This paper presents a new boundary integral equation (BIE) method for simulating particulate and multiphase flows through periodic channels of arbitrary smooth shape in two dimensions. The authors consider a particular system--multiple vesicles suspended in a periodic channel of arbitrary shape--to describe the numerical method and test its performance. Rather than relying on the periodic Green's function as classical BIE methods do, the method combines the free-space Green's function with a small auxiliary basis and imposes periodicity as an extra linear condition. As a result, we can exploit existing free-space solver libraries, quadratures, and fast algorithms and handle a large number of vesicles in a geometrically complex channel. Spectral accuracy in space is achieved using the periodic trapezoid rule and product quadratures, while a first-order semi-implicit scheme evolves particles by treating the vesicle-channel interactions explicitly. New constraint-correction formulas are introduced that preserve reduced areas of vesicles, independent of the number of time steps taken. By using two types of fast algorithms--(i) the fast multipole method for the computation of the vesicle-vesicle and the vesicle-channel hydrodynamic interaction, and (ii) a fast direct solver for the BIE on the fixed channel geometry--the computational cost is reduced to O(N) per time step, where N is the spatial discretization size. Moreover, the direct solver inverts the wall BIE operator at t = 0, stores its compressed representation, and applies it at every time step to evolve the vesicle positions, leading to dramatic cost savings compared to classical approaches. Numerical experiments illustrate that a simulation with N = 128;000 can be evolved in less than a minute per time step on a laptop. [ABSTRACT FROM AUTHOR]