1. Introduction
Microfluidic devices rely inherently on confined flow dynamics for a wide range of applications, including disease diagnostics (Gharib et al. Reference Gharib, Bütün, Muganlı, Kozalak, Namlı, Sarraf, Ahmadi, Toyran, van Wijnen and Koşar2022), single-cell encapsulation and analysis (Ling et al. Reference Ling, Geng, Chen, Du and Xu2020), droplet-based microreactors (Liu, Xiang & Ni Reference Liu, Xiang and Ni2020), and the fabrication of tailored particulate matter (Li et al. Reference Li2018). Often, shallow rectangular channels are utilized in which suspended objects flow at roughly the same depth, and hence the same optical focus, thus allowing microscopy-based characterization and monitoring. In such confined flows, it is crucial to understand the subtle hydrodynamic interactions that control the dynamics and positioning of particles, relative to one another and relative to the channel walls, in order to optimize the performance of such devices.
In confined flows, at volume fractions above the dilute limit, interparticle hydrodynamic interactions can lead to the emergence of collective behaviour and particle ordering. Early studies focusing on rigid spherical particles within quasi-one-dimensional (quasi-1-D) channels showed that reductions in drag forces occur due to shielding interactions, which vanish at distances beyond the channel diameter (Wang & Skalak Reference Wang and Skalak1969; Leichtberg, Pfeffer & Weinbaum Reference Leichtberg, Pfeffer and Weinbaum1976; Cui, Diamant & Lin Reference Cui, Diamant and Lin2002). In quasi-two-dimensional (quasi-2-D) flows (i.e. planar Poiseuille flow), rigid particles induce a dipolar disturbance to the flow field that is unscreened in the lateral dimensions (Janssen et al. Reference Janssen, Baron, Anderson, Blawzdziewicz, Loewenberg and Wajnryb2012). This dipolar field, which is superimposed on the bare flow field, arises due to the confinement (Liron & Mochon Reference Liron and Mochon1976; Cui et al. Reference Cui, Diamant, Lin and Rice2004), is maximal at the centreplane, and vanishes at the walls. Collective phenomena have been shown to arise due to these dipolar fields in prearranged particle lattices, including particle-displacement waves and dislocation propagation (Baron, Bławzdziewicz & Wajnryb Reference Baron, Bławzdziewicz and Wajnryb2008; Blawzdziewicz et al. Reference Blawzdziewicz, Goodman, Khurana, Wajnryb and Young2010), as well as reduced friction coefficients (Bhattacharya Reference Bhattacharya2008; Kohale & Khare Reference Kohale and Khare2010) when particles are propagated by a prescribed velocity rather than by the surrounding fluid. However, in planar Poiseuille flow, arrays of rigid particles are unstable to lateral perturbations (Janssen et al. Reference Janssen, Baron, Anderson, Blawzdziewicz, Loewenberg and Wajnryb2012), therefore order-to-disorder transitions are observed (Blawzdziewicz et al. Reference Blawzdziewicz, Goodman, Khurana, Wajnryb and Young2010), while the opposite, i.e. disorder-to-order transitions, are not observed – although in confined shear flows, spontaneous ordering of rigid particles has been demonstrated (Cheng et al. Reference Cheng, Xu, Rice, Dinner and Cohen2012).
Deformable particles – including red blood cells (RBCs), droplets, vesicles and capsules – exhibit strikingly different flow behaviour compared to rigid particles, and in particular a much higher propensity for collective ordering in confined flow. The 1-D ordering of RBCs in cylindrical channels has been studied using computer simulations, revealing a flow-induced ordering with a stable and finite cell–cell separation (McWhirter, Noguchi & Gompper Reference McWhirter, Noguchi and Gompper2009), as well as the role of cell deformation on the apparent viscosity (Pozrikidis Reference Pozrikidis2005). The finite separations observed by McWhirter et al. (Reference McWhirter, Noguchi and Gompper2009) (and later by Janssen et al. (Reference Janssen, Baron, Anderson, Blawzdziewicz, Loewenberg and Wajnryb2012) for droplets in parallel-wall channels) are not observed in rigid suspensions in confinement, whereby particle–particle separations vanish to zero (Blawzdziewicz et al. Reference Blawzdziewicz, Goodman, Khurana, Wajnryb and Young2010; Janssen et al. Reference Janssen, Baron, Anderson, Blawzdziewicz, Loewenberg and Wajnryb2012) or very small distances due to short-range screened electrostatic repulsions (Cui et al. Reference Cui, Diamant and Lin2002). (It is noted that rigid particles flowing through rectangular channels at relatively lower confinement ratios do form stable cross-streamline pairs and trains with an equilibrium finite separation, e.g. Schaaf & Stark Reference Schaaf and Stark2020.) Experimental observations of ordered arrays of RBCs (Abkarian et al. Reference Abkarian, Faivre, Horton, Smistrup, Best-Popescu and Stone2008; Tomaiuolo et al. Reference Tomaiuolo, Lanotte, Ghigliotti, Misbah and Guido2012; Iss et al. Reference Iss, Midou, Moreau, Held, Charrier, Mendez, Viallat and Helfer2019) and droplets (Beatus, Tlusty & Bar-Ziv Reference Beatus, Tlusty and Bar-Ziv2006; Beatus, Bar-Ziv & Tlusty Reference Beatus, Bar-Ziv and Tlusty2007, Reference Beatus, Bar-Ziv and Tlusty2012; Hashimoto et al. Reference Hashimoto, Garstecki, Stone and Whitesides2008; Raven & Marmottant Reference Raven and Marmottant2009) have also been reported.
Janssen et al. (Reference Janssen, Baron, Anderson, Blawzdziewicz, Loewenberg and Wajnryb2012) proposed that particle deformation in confined planar Poiseuille flows results in an additional (superimposed) quadrupole field that gives rise to long-range relative motions (both lateral and longitudinal) and stable finite separations. The quadrupole interactions also stabilize 1-D arrays from lateral perturbations at high confinement (Janssen et al. Reference Janssen, Baron, Anderson, Blawzdziewicz, Loewenberg and Wajnryb2012), but are decreasingly effective in low confinement (Bryngelson & Freund Reference Bryngelson and Freund2016). The ordering of initially random monolayer distributions of deformable capsules and RBCs into 1-D arrays has been demonstrated in confined planar Poiseuille flows (Iss et al. Reference Iss, Midou, Moreau, Held, Charrier, Mendez, Viallat and Helfer2019). Interestingly, 1-D and 2-D arrays also develop in confined shear flows with deformable capsules and droplets (Shen et al. Reference Shen, Fischer, Farutin, Vlahovska, Harting and Misbah2018; Singha et al. Reference Singha, Malipeddi, Zurita-Gotor, Sarkar, Shen, Loewenberg, Migler and Blawzdziewicz2019; Ishida et al. Reference Ishida, Matsumoto, Matsunaga and Imai2022), driven also by quadrupolar flow disturbances arising from particle deformation (albeit, the particle shapes are different from those within planar Poiseuille flow). Three-dimensional (3-D) ordering of elastic capsules has also been observed in confined shear flows at higher particle volume fractions (Liao et al. Reference Liao, Wu, Chien, Huang and Chen2017). In rectangular channels at lower confinement, increasing deformability enhances the stability of particle pairs (Patel & Stark Reference Patel and Stark2021), as well as the probability of pair formation when an initial cross-streamline separation exists (Owen & Krüger Reference Owen and Krüger2022). In confined planar flow, the stable separation between a pair of deformable particles varies with the degree of confinement, and for a particular range of confinements (of rather low values), there appear to be two coexisting stable separations: a shorter one and a longer one (Aouane et al. Reference Aouane, Farutin, Thiébaud, Benyoussef, Wagner and Misbah2017).
Previous studies have yet to produce a systematic and quantitative analysis of the ordering of deformable particle suspensions in planar Poiseuille flow. Here, such an effort is made through a parametric study in which both the confinement ratio and the particle deformability are varied systematically within ranges wide enough to fully characterize the ordering. Large-scale computer simulations are employed that track the evolution of hundreds of spherical elastic capsules from an initially random distribution to an ordered state consisting of 1-D train assemblies. A numerical order parameter is proposed to quantify the degree of ordering within the suspension, and used to identify trends within the parameter space. Uniform suspensions (i.e. zero dispersity in capsule size and deformability) are considered first. An optimal confinement ratio is identified corresponding to a maximum degree of ordering. Next, the impact of dispersity in both the key parameters (confinement ratio and deformability) is studied, using both Gaussian and bimodal distributions. A sharp decrease in ordering is found at rather low values of dispersity in capsule size (i.e. confinement ratio), which is attributable to variations in capsule velocity across the distribution range. In comparison, a relatively smaller decrease in ordering is observed for increasing dispersity in capsule deformability, although this depends on the mean deformability for the suspension. The results of this study are applicable directly to microfluidic applications whereby particle ordering is desired within homogeneous and heterogeneous particle mixtures.
2. Methods
The simulations utilize a lattice Boltzmann method (LBM) coupled with an immersed boundary method (IBM) technique to model the deformable capsule suspensions, an approach used commonly for such non-stationary fluid–structure interaction problems (Krüger et al. Reference Krüger, Kusumaatmaja, Kuzmin, Shardt, Silva and Viggen2017). Full details and validation of the model can be found in a recent paper by the author (Millett Reference Millett2023) or other similar studies (e.g. Krüger, Kaoui & Harting Reference Krüger, Kaoui and Harting2014; Schaaf & Stark Reference Schaaf and Stark2017), but the important points are discussed below.
The fluid domain is discretized with a 3-D lattice consisting of 
$N_x \times N_y \times N_z$ lattice sites with uniform spacing 
$\Delta x$. The D3Q19 stencil is used, hence each site stores a set of 19 particle distribution functions 
$f_i({\boldsymbol {x}},t)$ that are updated with discrete time steps (
$\Delta t$) according to the equation
\begin{equation} f_i({\boldsymbol{x}} + {\boldsymbol{c}}_i\,\Delta t, t + \Delta t) - f_i({\boldsymbol{x}},t)= \varOmega_i ({\boldsymbol{x}},t) + F_i\,\Delta t,\end{equation}
where 
${\boldsymbol {c}}_i$ are the discrete velocity vectors, 
$F_i$ is the Guo forcing term (Guo, Zheng & Shi Reference Guo, Zheng and Shi2002), and 
$\varOmega _i$ is the Bhatnagar–Gross–Krook collision operator defined as
\begin{equation} \varOmega_i ({\boldsymbol{x}},t) ={-} \frac{\Delta t}{\tau} \left[ f_i({\boldsymbol{x}},t) - f_i^{eq}({\boldsymbol{x}},t) \right], \end{equation}
with 
$\tau$ a relaxation time, and 
$f_i^{eq}$ the equilibrium particle distribution functions for a given fluid velocity 
${\boldsymbol {u}}$ and density 
$\rho$. The kinematic fluid viscosity is related to the relaxation time 
$\tau$ by
\begin{equation} \nu = c_s^2 \left(\tau - \frac{\Delta t}{2} \right), \end{equation}
where 
$c_s$ is the lattice speed of sound. In all simulations herein, the time step size and lattice spacing are taken as unity: 
$\tau = 1\,\Delta t$, 
$c_s = (1/\sqrt {3})\,\Delta x/\Delta t$, and 
$\nu = (1/6)\,\Delta x^2/\Delta t$.
As shown in figure 1(a), the suspension is confined between two parallel walls that are normal to the 
$z$-direction. At these walls, halfway bounce-back rules are applied to enforce zero-slip conditions (Ladd Reference Ladd1994). In the 
$x$- and 
$y$-directions, periodic boundary conditions are implemented, hence the walls can be assumed to extend infinitely in these two directions. In LBM units, the dimensions of the domain are 
$N_x = N_y = 400\,\Delta x$, and 
$N_z$ is varied systematically over the range 16–28 times 
$\Delta x$. The channel depth is denoted by 
$2h$ (see figure 1b). The confinement ratio 
$a/h$ (where 
$a$ is the resting capsule radius) is therefore varied but kept in a range that preserves a monolayer of capsules.
Schematic of the simulation domain. (a) Stationary walls exist normal to the 
$z$-direction, and periodic boundaries are enforced in the 
$x$- and 
$y$-directions. Flow is in the 
$x$-direction. (b) The channel depth, denoted as 2
$h$, is varied systematically in this study. Each capsule membrane is discretized with 1280 faces and 642 nodes.
Flow within the channel is driven by a uniform body force density 
$f_x^b$ applied in the 
$x$-direction at every lattice site. Absent particles, this results in a parabolic flow profile
\begin{equation} u(z) = u_{max} \left[1 - \left( \frac{z-h}{h} \right)^2 \right],\end{equation}
assuming that the origin 
$z = 0$ exists at the bottom wall. The maximum velocity at the channel centreplane can be calculated as
\begin{equation} u_{max} = \frac{Re_0\,\nu}{h},\end{equation}
where 
$Re_0$ is the bare Reynolds number, which is assigned a value 
$Re_0 = 1$ in all simulations herein. Hence the largest value of 
$u_{max}$ in this study is 0.0208 
$\Delta x/\Delta t$, thereby ensuring negligible fluid compressibility (Krüger et al. Reference Krüger, Kaoui and Harting2014). The body force density for given 
$Re_0$ and 
$h$ is calculated by 
$f_x^b = (2u_{max}^2 \rho )/(Re_0\,h)$.
The capsules are modelled with the IBM, and the membrane of each capsule is discretized with a triangular mesh consisting of 1280 faces and 642 nodes (see figure 1b). At rest, the capsules assume a spherical shape. In most of the results below, the capsule radius is set to a value 
$a = 6\,\Delta x$ (the exception is when dispersities in capsule radii within a suspension are considered; see § 3.3, where values for 
$a$ are specified). The capsule membranes are considered to be hyperelastic shells and therefore deform in flow conditions. Considering a single capsule, the deformation energy is
\begin{equation} E = E_s + E_b + E_v, \end{equation}
with contributions from in-plane shear (
$E_s$), out-of-plane bending (
$E_b$), and volumetric changes (
$E_v$). Full details of these terms can be found in Millett (Reference Millett2023). Briefly, the in-plane shear energy is described by the Skalak model (Skalak et al. Reference Skalak, Tozeren, Zarda and Chien1973), which consists of a shear modulus 
$\kappa _s$ and an area-dilation modulus 
$\kappa _a$ (here, 
$\kappa _a/\kappa _s = 2$). A force is calculated for each node (positioned at 
$X$) on a capsule membrane, 
${\boldsymbol {F}}^m(X) = {\boldsymbol {F}}^m_{el}(X) + {\boldsymbol {F}}^m_{rep}(X)$, where 
${\boldsymbol {F}}^m_{el}(X) = -\partial E/\partial X$ is the elastic deformation energy, and 
${\boldsymbol {F}}^m_{rep}(X) = -\partial E_r/\partial X$ is a short-range soft repulsion force between nodes of neighbouring capsules and between nodes and the walls. Here, 
$E_r = A(r-D)^2/2$ is a harmonic energy with cutoff length 
$D = \Delta x$ and 
$A = 0.1$ in reduced units. Standard to the IBM, this node force (for every node) is extrapolated to nearby sites on the LBM lattice and included in the Guo forcing term in (2.1), thereby modifying the local flow field. In return, the fluid velocity field is interpolated to each IBM node in order to propagate its position for that time step. Here, a trilinear spread function (Krüger, Varnik & Raabe Reference Krüger, Varnik and Raabe2011) is used for both the extrapolation and interpolation couplings.
The deformability of a capsule is conveniently described by the capillary number
\begin{equation} Ca = \frac{\rho \nu u_{max} a}{h \kappa_s},\end{equation}
hence 
$Ca$ is inversely proportional to the in-plane shear modulus 
$\kappa _s$, which is calculated for a given 
$Ca$ using (2.7). To explore a span of capsule deformabilities, the results below use a range of capillary numbers from 
$Ca = 0.008$ (stiffest) to 
$Ca = 0.3$ (softest).
The initial state of each simulation is a random distribution of capsules within the channel (ensuring no overlap among neighbouring capsules or with the channel walls) and the fluid velocity profile specified by (2.4). The fluid density is 
$\rho = 1$ throughout the domain. In all simulations, the volume fraction 
$\phi = N(4{\rm \pi} a^3)/(3V)$ is held fixed at 
$\phi = 0.1$ (where 
$N$ is the number of capsules, and 
$V$ is the volume of the fluid domain). This corresponds to a semi-dilute concentration. Hence the number of capsules in a simulation depends on the channel depth 
$2h$. For the shallowest channels (
$2h = 16\,\Delta x$, or 
$a/h = 0.75$ assuming 
$a = 6\, \Delta x$), the number of capsules is 
$N = 283$, while for the deepest channels (
$2h = 28\, \Delta x$, or 
$a/h = 0.43$), the number of capsules is 
$N = 495$. The wide channel dimensions in the 
$x$- and 
$y$-directions (400 
$\Delta x$ or 66.67
$a$) minimize any potential effects of periodicity. In some of the results below, a bimodal distribution of capsule radii is simulated, and in these cases, 
$\phi$ remains fixed at 0.1 for each number fraction of small versus large capsules. Each simulation is executed for 
$2 \times 10^6$ time steps.
As demonstrated by Iss et al. (Reference Iss, Midou, Moreau, Held, Charrier, Mendez, Viallat and Helfer2019), 1-D trains are the preferred order arrangement for deformable particles in planar Poiseuille flow. Here, in order to explore fully the degree of ordering and its dependency on the channel and suspension conditions, it is necessary to identify which capsules belong to trains and which do not. For a particular capsule 
$i$, a search of neighbouring capsules is conducted to determine if any are located within a certain range in front of and behind capsule 
$i$, as shown in figure 2. A cutoff radius 
$r_{cut} = 3.5a$ is defined, as well as an angular range 
$\theta = \pm 15^{\circ }$ centred along the 
$x$-direction. The relative positions are calculated using the centre of mass of each capsule. Capsule 
$i$ is tagged as belonging to a train if it has two neighbouring capsules (one in front and one behind) located within these ranges. In order to also account for capsules at the front and rear of a train, an alternative criterion is defined: if capsule 
$i$ has one neighbour in this range, and that neighbour is in a train, then capsule 
$i$ is also tagged as belonging to a train. In the figures below, capsules belonging to trains are shaded grey, and capsules not belonging to trains are shaded blue. An order parameter is defined that characterizes the overall degree of ordering in the suspension by calculating the fraction of capsules belonging to a train:
\begin{equation} \varPhi_t = \frac{N_{train}}{N}, \end{equation}
where 
$N_{train}$ is the number of capsules belonging to a train. Here, 
$\varPhi _t$ is calculated at periodic intervals throughout every simulation, and after an equilibration time, its mean and standard deviation values are reported below.
Schematic illustrating the criteria for determining whether a capsule belongs to a train. (a) For each capsule, a search is conducted to find if other capsules are located within an axial range of 
$r_{cut}$ and a radial range of 2
$\theta$, both in front of and behind the capsule. (b) If two neighbouring capsules are within this range (one in front and one behind), then that capsule is flagged to be in a train. In addition, if a capsule has one neighbour in this range, and that neighbour is in a train, then the capsule is also flagged to be in a train. Herein, grey shading denotes that a capsule is in a train, and blue shading denotes that it is not.
It is noted that the chosen values of 
$r_{cut}$ and 
$\theta$ are somewhat arbitrary. The value of 
$r_{cut}$ should be larger than the equilibrium spacing of capsules in a train (shown below to be 2.3
$a$–2.8
$a$; see figure 8), but not so large that next-nearest capsules are counted. There is no quantitative rationale for using 
$\theta = \pm 15^{\circ }$ rather than another value between 
$0^{\circ }$ and 
$90^{\circ }$, other than providing directional bias to the train criteria along the flow direction (which would be eliminated if, for example, 
$\theta = \pm 90^{\circ }$).
3. Results
3.1. Dynamics of a single capsule
First, the dynamics of a single capsule in planar Poiseuille flow is investigated. In figure 3(a), the steady-state non-dimensionalized capsule velocity 
$\tilde {u}_c= u_c/u_{max}$ is plotted versus 
$a/h$ for varying values of deformability (where 
$u_c$ is the capsule velocity in the flow direction). First, it can be seen that 
$\tilde {u}_c$ decreases with increasing confinement. This is to be expected due to the increased viscous interactions with the walls and increased disturbance of flow with increasing confinement (also, it is evident that as 
$a/h$ approaches zero, 
$\tilde {u}_c$ should approach 1). Second, 
$\tilde {u}_c$ has very little dependency on 
$Ca$ for 
$a/h \leq 0.6$, and an increasing dependency for 
$a/h > 0.6$. A very similar result was found by Griggs, Zinchenko & Davis (Reference Griggs, Zinchenko and Davis2007) in their simulations of single drops between parallel plates. In particular, with increasing confinement, the decrease in velocity is greatest for small 
$Ca$, as stiffer capsules deform less and therefore their membrane surfaces remain closer to the walls. The solid red line in figure 3(a) is a second-order polynomial fit to the 
$Ca = 0.3$ data that intersects the 
$y$-axis at 
$\tilde {u}_c = 1$ (
$\tilde {u}_c = -0.3398(a/h)^2 - 0.0216(a/h) + 1$). This fit will be used for comparison in results in later sections.
(a) Velocity of a single capsule (
$\tilde {u}_c = u_c/u_{max}$) versus 
$a/h$ for varying values of 
$Ca$, revealing that 
$\tilde {u}_c$ has an inverse dependency on 
$a/h$, negligible dependency on 
$Ca$ for 
$a/h \leq 0.6$, and moderate dependency on 
$Ca$ for 
$a/h > 0.6$. The solid red line is a second-order polynomial fit to the 
$Ca = 0.3$ data. (b) Slice in the 
$x$–
$z$ plane bisecting the capsule used for the velocity plots in (c–f). (c) Induced flow field and (d) comoving flow field for 
$Ca = 0.3$ and 
$a/h = 0.6$. (e) Induced flow field and (f) comoving flow field for 
$Ca = 0.01$ and 
$a/h = 0.6$. The colour scale bars represent velocity magnitudes normalized by 
$u_{max}$.
Figures 3(c–f) provide views of two different flow fields on the 
$x$–
$z$ plane that bisects the capsule, for 
$a/h = 0.6$. For 
$Ca = 0.3$, the induced flow field is shown in figure 3(c) defined as the difference between the actual fluid velocity and the particle-free velocity profile, (2.4). Also shown in figure 3(d) is the comoving velocity field defined as the velocity relative to the frame moving with the capsule. Comparing the images for the two 
$Ca$ values, one does not see significant differences in the velocity fields, apart from slight differences due to capsule deformation (again, for this confinement, there is little dependency of velocity with 
$Ca$). Examining the induced flow fields (figures 3c,e), it can be seen that the presence of a capsule reduces the magnitude of velocity along the channel centreplane relative to the particle-free flow profile. Also, the presence of a capsule increases the magnitude of velocity at the top and bottom regions of the capsule relative to the particle-free flow profile (see the red arrows). Examining the comoving velocity fields (figures 3d,f) reveals that there is essentially no relative internal flow within the capsule, a result also shown by Aouane et al. (Reference Aouane, Farutin, Thiébaud, Benyoussef, Wagner and Misbah2017), who showed that on the other hand, drops in similar conditions exhibit internal recirculation zones.
3.2. Uniform suspensions
Next, the ordering behaviour of uniform suspensions (i.e. zero polydispersity in capsule size and deformability) is investigated. Here, both the confinement and deformability are systematically varied. The set of confinement ratios considered is 
$a/h = [0.43, 0.46, 0.50, 0.55, 0.60, 0.67, 0.75]$ (corresponding to 
$N_z = [28, 26, 24, 22, 20, 18,$ 
$16]\,\Delta x$, respectively, with a constant capsule radius 
$a = 6\,\Delta x$). Time is non- dimensionalized by dividing by the advection time, or the time required for an infinitesimal particle located at the centreplane to travel a distance 
$a$:
\begin{equation} \tilde{t} = t \left( \frac{u_{max}}{a} \right). \end{equation}
Each simulation is run for 
$2\times 10^6$ steps. Because 
$u_{max}$ varies inversely with 
$h$ (see (2.5)), the final non-dimensionalized time at the end of the simulation is 
$\tilde {t} = [3968, 4274, 4630, 5051, 5556, 6173, 6944]$ corresponding to 
$N_z = [28, 26, 24, 22, 20,$ 
$18, 16] \,\Delta x$, respectively.
Figure 4 displays top-down views of the capsule distributions at progressive instances of time for confinement ratio 
$a/h = 0.75$ and deformability 
$Ca = 0.3$ (i.e. the highest confinement and the softest capsules). In this simulation, the number of capsules is 
$N = 283$. The initial disordered state at 
$\tilde {t} = 0$ is shown in figure 4(a). Early on in the simulation (figures 4b,c), 1-D train assemblies begin to develop as evidenced by the grey-shaded capsules. At this preliminary stage, the trains are not too long and for the most part consist of less than 10 capsules. At the end of the simulation (figure 4(e), corresponding to 
$\tilde {t} = 6944$), the number of trains increases, as does the train length. At the end of the simulation, the capsules have travelled the length of the simulation box in the 
$x$-direction approximately 100 times. At this point, many of the trains extend the entire length of the channel in the 
$x$-direction. It can be seen that disorder remains in some regions of the channel where clusters of non-aligned blue-shaded capsules exist. Furthermore, very localized regions of defects can also be observed in the system, including dislocations associated with locations where two trains merge into one train resulting in a V-shaped capsule arrangement.
Top-down views illustrating the time evolution of a capsule suspension with 
$Ca = 0.3$ and 
$a/h = 0.75$. Snapshots correspond to progressive instances of time: (a) 
$\tilde {t} = 0$, (b) 
$\tilde {t} = 347$, (c) 
$\tilde {t} = 694$, (d) 
$\tilde {t} = 1736$, and (e) 
$\tilde {t} = 6944$. (f) Order parameter 
$\varPhi _t$ versus time 
$\tilde {t}$ with confinement 
$a/h = 0.75$ for each value of 
$Ca$ used in the study. In subsequent discussions, 
$\varPhi _t$ is taken as an average over the second half of the simulation, as indicated by the horizontal lines.
The order parameter 
$\varPhi _t$ is plotted versus time in figure 4(f) for each of the 
$Ca$ values. For each deformability value, 
$\varPhi _t$ exhibits a rapid early increase and eventually levels off to a steady-state value. For stiffer capsules, including 
$Ca = 0.01$ and 0.05, the steady-state value is rather small. For all results below, the first half of the simulation (
$10^6 \,\Delta t$) is treated as an equilibration period, and throughout the second half of the simulation, a running average of 
$\varPhi _t$ as well as its standard deviation is calculated. The horizontal lines in figure 4(f) represent the overall average throughout the second half of the simulation. Images of the final capsule configurations corresponding to the simulation data shown in figure 4(f) are shown in supplementary material available at https://doi.org/10.1017/jfm.2023.1052 (see figure S1).
Figure 5 shows representative snapshots of capsule configurations for 
$Ca = 0.05$ for varying confinement ratios: 
$a/h = 0.75$, 0.6 and 0.43. These images were taken at the end of the simulations. These capsules are relatively stiff and therefore are expected to exhibit less ordering compared to the results shown in figures 4(a–e). Figures 5(a–c) present top-down views of the system, and figures 5(d–f) present side-on views. For 
$a/h = 0.75$ (figures 5a,d), only a few short trains are observed, and certainly the degree of ordering is less than that seen in figure 4 for 
$Ca = 0.3$. However, for the slightly lower confinement ratio 
$a/h = 0.6$ (figures 5b,e), a greater degree of ordering appears to have developed. This confinement ratio corresponds to an intermediate channel depth (
$N_z = 20\,\Delta x$) within the range investigated in this study, and the number of capsules is increased to 
$N = 354$ in order to maintain the same value of 
$\phi$. Further decreasing the confinement ratio to 
$a/h = 0.43$ (figures 5c,f), however, does not yield an additional increase in ordering; rather, there is a clear reduction in train development. This confinement ratio corresponds to the deepest channel considered in this study (
$N_z = 28\,\Delta x$). As can be seen in figure 5(f), some capsules become displaced from the channel centreplane towards the top and bottom walls, thus precluding the formation of a single monolayer of capsules. These capsules possessed the second-lowest deformability value in our study, hence they did not exhibit significant flow-induced deformation, yet capsules displaced away from the centreplane towards the walls do exhibit an ellipsoidal shape (see figure 5f) associated with the higher levels of shear stress near the walls. At this lower level of confinement, capsules were observed to pass one another occasionally via mutual displacements in the 
$z$-direction.
Simulation snapshots of capsule configurations for 
$Ca = 0.05$: (a–c) top-down views, and (d–f) side-on views. All images were taken at the end of the simulation run. The degree of confinement is varied: (a,d) 
$a/h = 0.75$, corresponding to the shallowest channel; (b,e) 
$a/h = 0.6$, corresponding to an intermediate channel depth; and (c,f) 
$a/h = 0.43$, corresponding to the largest channel depth. Note that all simulations held 
$\phi = 0.1$, hence the appearance of more capsules in (c) compared to (a). Capsules in trains are shaded grey, and capsules not in trains are shaded blue.
The images in figure 5 suggest that an optimal value of confinement exists to maximize the ordering process. Similar images are shown in figure 6, although for suspensions with higher deformability, 
$Ca = 0.3$. Comparing figures 5 and 6, it is clear that the softer capsules exhibit a higher degree of ordering for each confinement ratio. Again, the confinement ratio 
$a/h = 0.6$ (figures 6b,e) appears to promote the highest degree of ordering. The side-on views in figures 6(d–f) reveal that the capsules assume a folded shape, and that the amount of deformation increases with increasing confinement. The capsules assemble into monolayers even for the deepest channel (
$a/h = 0.43$), as opposed to the case for 
$Ca = 0.05$ shown in figure 5(f). The top-down view of the configuration for 
$a/h = 0.6$ (figure 6b) shows that even though the degree of ordering is quite high, there still remain a few defects associated with dislocations in the train assemblages. This is to be expected, as self-assembling systems generally require special strategies (i.e. directed self-assembly strategies) to eliminate defects on length scales above an order of magnitude larger than the individual units (Li & Müller Reference Li and Müller2015).
Same as figure 5 except 
$Ca = 0.3$, i.e. softer capsules. A higher degree of alignment is observed compared to 
$Ca = 0.05$ for each channel depth.
The order parameter 
$\varPhi _t$ is plotted as a function of 
$a/h$ and 
$Ca$ in figure 7. The data points and error bars correspond to the averaged values and standard deviations, respectively, calculated throughout the second half of the simulations. Figure 7(a) shows 
$\varPhi _t$ plotted versus the confinement ratio, revealing that indeed 
$a/h = 0.6$ represents an optimal confinement ratio for particle ordering for each value of deformability. Increasing or decreasing the confinement ratio from this value results in a reduction in ordering. Figure 7(b) shows 
$\varPhi _t$ plotted versus the deformability, displayed separately for each confinement ratio. Note that the 
$x$-axes in figure 7(b) have a logarithmic scale. Figure 7(b) illustrates that the order–disorder transition occurs rather precipitously as the particle deformability is varied within the approximate range 
$Ca = 0.05\unicode{x2013}0.1$. Overall, for each confinement ratio, the ordering increases with increasing capsule deformability. In addition, it can be seen that the error bars are largest when the level of ordering is intermediate (i.e. when 
$\varPhi _t = 0.4\unicode{x2013}0.7$). In these conditions, it was observed that capsule trains form and dissolve transiently. Conversely, the error bars are relatively small when 
$\varPhi _t$ is close to 1 (whereby the capsule trains were quite stable with only a few isolated defects) or when 
$\varPhi _t$ is close to zero (whereby the conditions discourage train formation and the system remains disordered).
(a) Order parameter 
$\varPhi _t$ versus confinement ratio 
$a/h$ for each 
$Ca$ used in this study. (b) Plots of 
$\varPhi _t$ versus 
$Ca$ for each confinement. The tendency for ordering depends strongly on confinement and deformability. Maximal values of 
$\varPhi _t$ are observed for the softest capsules (
$Ca = 0.2, 0.3$) at confinement 
$a/h = 0.6$. Deviating above or below this confinement ratio results in a decrease in 
$\varPhi _t$. The error bars represent the standard deviation of 
$\varPhi _t$.
In order to ensure that periodicity in the flow direction is not impacting the values of 
$\varPhi _t$, three additional simulations were run with a channel length twice as large (see figures S2 and S3 in the supplementary material). The 
$\varPhi _t$ values and the general capsule configurations were similar for the two channel lengths, hence the effects of periodicity are assumed to be negligible. Furthermore, a test simulation was run with a Weeks–Chandler–Anderson short-range repulsion potential rather than the soft harmonic potential, and the 
$\varPhi _t$ values and the general capsule configurations were similar for the two repulsion models (see figure S4 in the supplementary material).
To formulate an explanation for the observed peak in 
$\varPhi _t$ at confinement 
$a/h = 0.6$, an analysis is first carried out on the attraction between a single pair of capsules aligned in the flow direction. Here, a set of simulations were performed consisting of two capsules placed initially on the centreplane of the channel, aligned in the flow direction. The lateral dimensions of the channel are 120 
$\Delta x$ (
$=20a$) in both the 
$x$- and 
$y$-directions (using periodic boundary conditions), with varying channel depth. Figure 8(a) displays the relative velocity in the flow direction between the two capsules, 
$\tilde {u}_{12} = (u_2-u_1)/u_{max}$, plotted versus the separation distance in the flow direction 
$\varDelta _{\parallel }$. In figure 8(a), the confinement ratio is fixed at 
$a/h = 0.6$ while the deformability is varied as indicated. The curves show that pairs of softer capsules experience a greater attraction in the flow direction for a given confinement. Quantitatively, these curves are quite similar to previous results for droplets in planar channel flow (Janssen et al. Reference Janssen, Baron, Anderson, Blawzdziewicz, Loewenberg and Wajnryb2012). In addition, close inspection of figure 8(a) shows that the equilibrium spacing is slightly larger for the softest capsules, a result also seen by Janssen et al. (Reference Janssen, Baron, Anderson, Blawzdziewicz, Loewenberg and Wajnryb2012), arising from the long-range hydrodynamic interactions.
Relative velocity between a pair of capsules aligned in the flow direction, 
$\tilde {u}_{12} = (u_2-u_1)/u_{max}$. (a) Relative velocity versus separation distance for confinement 
$a/h = 0.6$ and varying deformability. (b) Relative velocity versus separation distance for 
$Ca = 0.3$ and varying confinement, 
$a/h = 0.43\unicode{x2013}0.75$. The inset plots the maximum relative velocity (i.e. the depth of the well) versus 
$a/h$. (c) Relative velocity in the 
$y$-direction (
$\tilde {v}_{12}$) between a pair of capsules displaced initially in the 
$y$-direction.
Figure 8(b) displays a similar plot; however, the deformability is fixed at 
$Ca = 0.3$ and the confinement is increased systematically, revealing that the attraction in the flow direction increases with increasing confinement. Finally, in figure 8(c), the capsule pair is displaced initially in the 
$y$-direction (perpendicular to flow) as well as in the 
$x$-direction (parallel to flow), and the relative 
$y$-direction velocity 
$\tilde {v}_{12} = (v_2-v_1)/u_{max}$ is plotted versus lateral displacement 
$\varDelta _{\bot }$. Again, the lateral attraction between a pair of capsules increases with increasing confinement, albeit the magnitude of the lateral attraction appears to be less than that of the in-flow direction attraction shown in figure 8(b).
Overall, the results of figure 8 suggest that the tendency for ordering should be greatest in channels with the highest confinement, rather than the peak observed at 
$a/h = 0.6$ in figure 7. However, the results for uniform suspensions were carried out with constant volume fraction 
$\phi = 0.1$. In all the simulations, the capsules assembled into monolayers (with the exception of the hardest capsules in the deepest channels). At fixed 
$\phi$, once a monolayer of capsules has formed, the planar density of capsules (number of capsules per area) is inversely proportional to the confinement ratio: 
$PD (a^{-2}) = {3\phi }/({2{\rm \pi} (a/h)})$, as derived in Appendix A. Thus at higher confinements, the average spacing between capsules in the monolayer plane is higher, which conceivably would reduce the ordering tendency and therefore reduce 
$\varPhi _t$. To test this possibility, additional simulations similar to those shown in figures 5 and 6 were run for confinement ratios 
$a/h > 0.6$. However, here the number of capsules (and hence the planar density) was held equal to the values used in the 
$a/h = 0.6$ case. The data shown in figure 9 reveal that the drop in 
$\varPhi _t$ at confinement ratios 
$a/h > 0.6$ mostly disappears when the planar density of capsules is held fixed, rather than the volume fraction.
Order parameter 
$\varPhi _t$ versus confinement ratio 
$a/h$ for (a) 
$Ca = 0.3$, (b) 
$Ca = 0.2$, and (c) 
$Ca = 0.1$. Here, two cases are considered. First, the volume fraction is held constant (
$\phi = 0.1$) as displayed in figure 7. Second, the monolayer (areal) density is held constant for confinement ratios 
$a/h \geq 0.6$ (this is accomplished by keeping the number of capsules in the simulation the same as in the case 
$a/h=0.6$). These results indicate that the observed decrease in 
$\varPhi _t$ for 
$a/h > 0.6$ when volume fraction is held constant is largely due to the decrease in monolayer density with increasing confinement.
3.3. Dispersity in capsule size
The above results demonstrate that ordering is optimized with high deformability (
$Ca = 0.3$) and a particular value for the confinement ratio (
$a/h = 0.6$). Attention is now turned to understand how dispersity of properties within the suspension alters the ordering behaviour. Starting with these optimal ordering conditions (i.e. 
$Ca = 0.3$ and 
$a/h = 0.6$), increasing levels of dispersity will be introduced to ascertain how decreasing uniformity impacts ordering. First, dispersity in capsule size will be considered. Two types of distributions in capsule radii were simulated: (i) a Gaussian distribution of particle radii centred at 
$a/h = 0.6$ with increasing values of standard deviation in radius 
$\sigma _{\tilde {a}}$; and (ii) a bimodal distribution of particle radii with two sub-populations of capsules of radii 
$a_1$ and 
$a_2$ (here, the average of 
$a_1/h$ and 
$a_2/h$ is equal to 0.6). For both distribution types, the deformability is held fixed at 
$Ca = 0.3$, and the channel depth is held fixed at 
$2h = 20\,\Delta x$. Note that in this set of simulations, the shear modulus 
$\kappa _s$ for each capsule is calculated using 
$Ca = 0.3$ and the mean radius 
$\bar {a}$, rather than its individual radius. Hence every capsule has the same shear modulus despite the variability in radius; see (2.7).
First, the ordering behaviour for suspensions with a Gaussian distribution in capsule radii will be discussed. Note that the standard deviation values herein are associated with the confinement ratio 
$\tilde {a} = a/h$ (not radius 
$a$). Figure 10 shows simulation snapshots (after equilibration) of capsule configurations for increasing dispersity in radii. Three values of standard deviation in 
$a/h$ are shown: 
$\sigma _{\tilde {a}} = 0.01$ (figures 10a,d,g), 
$\sigma _{\tilde {a}} = 0.03$ (figures 10b,e,h) and 
$\sigma _{\tilde {a}} = 0.06$ (figures 10c,f,i). Increasing dispersity in capsule size leads to a decrease in ordering, despite the high deformability for these capsules. For the highest value of 
$\sigma _{\tilde {a}} = 0.06$ (figures 10c,f,i), the capsules appear to aggregate into disordered clusters. In some cases, these clusters partially consist of some small train assemblies that are generally much shorter than those observed in the more uniform suspensions.
Capsule configurations in systems with Gaussian distributions in capsule radii. Here, 
$Ca = 0.3$ and 
$\bar {a}/h = 0.6$, with increasing standard deviations of radius: (a,d,g) 
$\sigma _{\tilde {a}} = 0.01$, (b,e,h) 
$\sigma _{\tilde {a}} = 0.03$, and (c,f,i) 
$\sigma _{\tilde {a}} = 0.06$. (a–c) Histograms for capsule radii (here, 
$p$ is the probability density such that the sum area of the bars equals 1). (d–f) Top-down views of the suspension at the end of simulation. (g–i) Close-up views of a sub-portion of the domain (indicated by the red boxes in d–f), in which capsules are coloured by the principal membrane tension 
$T_1/\kappa _s$.
Figures 10(g–i) show close-up views of a sub-portion of the domain indicated by the red squares in figures 10(d–f). These close-up views colour the capsules according to the principal tension in the capsule membranes, 
$T_1$, dividing by 
$\kappa _s$ for non-dimensionalization. Examining figures 10(a,d,g), when the dispersity in capsule size is low, the membrane tension is quite uniform from capsule to capsule. For higher dispersity (figures 10c,f,i), it can been seen that the membrane tension is much more varied amongst the capsules. In particular, larger capsules exhibit higher values of 
$T_1$ due to the fact that they have a higher (local) confinement ratio and therefore a higher deformation. Smaller capsules experience a smaller confinement ratio, therefore they exhibit less deformation and lower membrane tension.
To quantify how the ordering varies with dispersity in capsule radius, figure 11(a) plots 
$\varPhi _t$ versus 
$\sigma _{\tilde {a}}$ for three values of deformability. This reveals another precipitous order–disorder transition with increasing 
$\sigma _{\tilde {a}}$, particularly in the range 
$\sigma _{\tilde {a}} = 0.00\unicode{x2013}0.02$. For the slightly less deformable suspensions with 
$Ca = 0.1$, the order–disorder transition occurs with even less dispersity in capsule size, in the approximate range 
$\sigma _{\tilde {a}} = 0.00\unicode{x2013}0.01$. Interestingly, after the significant initial drop in 
$\varPhi _t$ with increasing 
$\sigma _{\tilde {a}}$, the values of 
$\varPhi _t$ level off to a more constant lower value (approximately 
$\varPhi _t = 0.2\unicode{x2013}0.4$, depending on 
$Ca$) that does not change much with further increasing 
$\sigma _{\tilde {a}}$.
(a) Order parameter 
$\varPhi _t$ versus standard deviation of capsule radius 
$\sigma _{\tilde {a}}$ for 
$Ca = 0.1$, 0.2 and 0.3. The degree of ordering drops precipitously with 
$\sigma _{\tilde {a}}$. (b) For 
$\sigma _{\tilde {a}} = 0.01$ and 
$Ca = 0.3$, the probability density 
$p$, capsule velocity 
$\tilde {u}_c$, and 
$\varPhi _t$ are plotted as histograms versus 
$a/h$. The error bars for 
$\tilde {u}_c$ represent standard deviation for velocity within each bin. (c) Same as (b) except with 
$\sigma _{\tilde {a}} = 0.06$. In (b,c), the red lines in the middle row indicate the velocity of a single capsule (see figure 3a), and the red lines in the bottom row indicate the value of 
$\varPhi _t$ for the suspension as a whole, taken from (a).
A statistical analysis was done to understand the variation of capsule dynamics and ordering throughout the suspensions. Figures 11(b) and 11(c) provide three histogram plots versus 
$a/h$: the probability density displaying the dispersity in capsule radius (top), the average capsule velocity 
$\tilde {u}_c$ (middle), and the order parameter 
$\varPhi _t$ (bottom). Figure 11(b) displays data for a suspension with low dispersity, 
$\sigma _{\tilde {a}} = 0.01$, where a high degree of ordering was observed (
$\varPhi _t = 0.89$ for the suspension as a whole). For this system, the capsule velocity as well as 
$\varPhi _t$ was very uniform within the narrow range of 
$a/h$. On the other hand, figure 11(c) provides data for a suspension with higher dispersity, 
$\sigma _{\tilde {a}} = 0.06$, and a lower degree of ordering (
$\varPhi _t = 0.37$ for the suspension as a whole). In this system, a much greater variation in capsule velocity and ordering was observed in the range of 
$a/h$. Overall, in this suspension, larger capsules travelled with lower velocity and were more likely to be part of a train compared with the smaller capsules. Smaller capsules experience less confinement, hence they travel with higher velocity. It seems apparent that a wide distribution of capsule velocities works to hinder the ordering process and limit the development of train structures. Smaller capsules are often seen trailing larger capsules, but such structures exist only transiently, as the smaller capsules eventually pass the larger capsules. Pair formation for rigid particles has been studied recently, with some conflicting results: the 2-D simulations of Chen, Lin & Hu (Reference Chen, Lin and Hu2021) suggest pair stability only when a larger particle is leading, whereas the 3-D simulations of Thota, Owen & Krüger (Reference Thota, Owen and Krüger2023) demonstrated pair stability only when the smaller particle is leading. Note that, compared to the present work, those studies considered lower confinement ratios (
$a/h < 0.4$) in which the lateral position was an influencing factor. Here, the capsules shown in figure 10 are essentially confined to the centreplane.
Another observation from figures 11(b) and 11(c) is that the capsule velocities in the suspensions were below the curve associated with 
$\tilde {u}_c$ for a single capsule (the red lines in the middle row of figures 11b,c) throughout the ranges of 
$a/h$. This will be addressed in § 3.4.
Next, the ordering behaviour of suspensions with bimodal distributions in capsule radii will be discussed. Here, the capsule population is split between two sub-populations, one with radius 
$a_1$, and the other with radius 
$a_2$ (there is no dispersity about either of these radii, hence the probability density histograms consist of only two bars). The average of 
$a_1/h$ and 
$a_2/h$ is 0.6 for all simulations. Also, 
$a_1 < a_2$ throughout. The fraction of capsules with radius 
$a_2$ is represented by 
$X_2$, which has been varied systematically in the results below. For these bimodal suspensions, the overall volume fraction is fixed at 
$\phi = 0.1$ regardless of the value of 
$X_2$, hence suspensions with a greater fraction of capsules with radius 
$a_1$ will have a larger number of capsules than suspensions with a greater fraction of capsules with radius 
$a_2$.
Figure 12 displays simulation snapshots for suspensions with 
$a_1/h = 0.55$, 
$a_2/h = 0.65$ and 
$Ca = 0.3$. Here, the fraction of capsules with radius 
$a_2$ is varied: 
$X_2 = 0.0$ in figures 12(a,d,g), 
$X_2 = 0.3$ in figures 12(b,e,h), and 
$X_2 = 0.9$ in figures 12(c,f,i). For these systems, the extreme values of 
$X_2$ (i.e. 0 and 1) correspond to uniform suspensions and therefore are expected to undergo a high degree of ordering, given the high capsule deformability. This can be seen when examining figures 12(a,d,g). On the other hand, figures 12(c,f,i) show images for 
$X_2 = 0.9$ that also exhibit significant ordering, but the presence of a small quantity of smaller capsules with radius 
$a_1$ does disturb the ordering process. It is less intuitive, however, how more equal quantities of the two radii will impact ordering. Figures 12(b,e,h) provide images of a suspension with 
$X_2 = 0.3$, showing clearly a significant reduction in ordering. The close-up images in figure 12 reveal several important features. First, as seen above, the smaller capsules experience less membrane tension due to their reduced confinement. Second, for the bimodal suspensions, it was common to see a larger capsule leading a train followed by one or more smaller capsules. Third, the trailing smaller capsules were often seen to move around the larger capsules, hence these train structures were not long-lived.
Capsule configurations for bimodal distributions in capsule radius with two populations: 
$a_1/h = 0.55$ and 
$a_2/h = 0.65$. The deformability is 
$Ca = 0.3$. The fraction of capsules with radius 
$a_2$ is varied: (a,d,g) 
$X_2 = 0.0$, (b,e,h) 
$X_2 = 0.3$, and (c,f,i) 
$X_2 = 0.9$. The volume fraction is fixed at 
$\phi = 0.1$, therefore the simulation in (a,d,g) has more capsules than in (d–i). A visual inspection indicates that the degree of ordering is lowest in (b,e,h), where the mix of smaller and larger capsules is more even; (g–i) indicate that larger capsules (with radius 
$a_2$) experience higher levels of membrane tension.
The order parameter 
$\varPhi _t$ for bimodal radii suspensions is plotted versus 
$X_2$ in figure 13(a) for three different sets of radii 
$a_1$ and 
$a_2$, as indicated in the legend. For each set of 
$a_1$ and 
$a_2$, the data show that the ordering is highest when 
$X_2 = 0$ or 
$1$, corresponding to the uniform suspensions. The value of 
$\varPhi _t$ decreases for intermediate values of 
$X_2$ associated with mixtures of smaller and larger capsules. A greater reduction in ordering is seen with increasing discrepancy between 
$a_1$ and 
$a_2$, as would be expected. The smallest difference in size, 
$a_1/h = 0.59$ and 
$a_2/h = 0.61$, corresponds to only approximately a 3
$\,\%$ size difference, and there was not much of a drop in 
$\varPhi _t$ across the range of 
$X_2$. The largest difference in size, 
$a_1/h = 0.55$ and 
$a_2/h = 0.65$, corresponds to approximately an 18 % size difference, which results in a significant decrease in 
$\varPhi _t$, particularly when 
$X_2 = 0.3$. An interesting observation is that the curves in figure 13(a) are not symmetric about the 
$X_2 = 0.5$ isopleth. It seems that adding a small amount of larger capsules to a suspension of smaller capsules will reduce the tendency for ordering more than by adding a small amount of smaller capsules to a suspension of larger capsules.
(a) Order parameter 
$\varPhi _t$ versus 
$X_2$ in systems with bimodal distributions in capsule radius. Three different sets of radii were considered, as indicated in the legend. Relative to uniform distributions (
$X_2 = 0$ or 1), bimodal distributions can exhibit considerably decreased ordering. (b) For 
$a_1/h = 0.55$, 
$a_2/h = 0.65$ and 
$X_2 = 0.3$, the probability density 
$p$, capsule velocity 
$\tilde {u}_c$, and 
$\varPhi _t$ are plotted as histograms versus 
$a/h$. (c) Same as (b) except with 
$X_2 = 0.8$. Capsules with larger radii 
$a_2$ move at lower velocity, but are more likely to belong to a train.
Figures 13(b) and 13(c) show histogram plots similar to those in figure 11 for two different values of 
$X_2$, for 
$a_1/h = 0.55$ and 
$a_2/h = 0.65$. For both compositions, the larger capsules travel at a lower velocity and are more likely to belong to a train assembly. Similar to the Gaussian distributions, the capsule velocities in these suspensions are below the curve associated with a single capsule.
3.4. Capsule dynamics in train assemblies
To address the above results concerning the discrepancy between the velocity of an isolated capsule and the velocity of capsules in suspensions that assemble into trains, an additional investigation is presented. Here, initially capsules are placed in train structures with uniform in-line spacing (
$\varDelta _{\parallel }$) and lateral spacing (
$\varDelta _{\bot }$), as indicated in figure 14(a). Here, the deformability is fixed at 
$Ca = 0.3$ and the confinement ratio is fixed at 
$a/h = 0.6$. Figure 14(a) shows the steady-state non-dimensionalized capsule velocity 
$\tilde {u}_c$ versus 
$\varDelta _{\parallel }/a$ for varying values of 
$\varDelta _{\bot }/a$. Note that the 
$x$-axis in figure 14(a) is in logarithmic scale. The data show that the velocity of capsules within trains decreases with decreasing 
$\varDelta _{\parallel }$ and decreasing 
$\varDelta _{\bot }$. Janssen et al. (Reference Janssen, Baron, Anderson, Blawzdziewicz, Loewenberg and Wajnryb2012) observed a similar result for a pair of droplets, namely a decrease in droplet velocity with decreasing droplet separation, which they attributed to a decrease in the hydrodynamic drag exerted by the surrounding flow on the pair. Given that both droplets and capsules are deformable particles that experience similar particle–particle interactions in Poiseuille flow (Aouane et al. Reference Aouane, Farutin, Thiébaud, Benyoussef, Wagner and Misbah2017), this explanation is likely also valid for capsules in infinitely long trains. Figure 14(a) also reveals that decreasing distance between adjacent trains 
$\varDelta _{\bot }$ also reduces the train velocity. Figures 14(c) and 14(d) show the induced flow field and the comoving flow field for a capsule train, respectively. Due to the close proximity of the capsules, the velocity fields do not transition back to the particle-free flow profile at points increasingly removed from the surface of any particular capsule, as was seen for the case of an isolated capsule (see figure 3, particularly along the centreplane). Hence it seems plausible that the induced drag force on a capsule due to the surrounding flow is reduced within the train assemblies, and that it decreases further with decreasing capsule separation.
(a) Capsule velocities when capsules are positioned in trains for 
$Ca = 0.3$ and 
$a/h = 0.6$. The in-line spacing of capsules is denoted as 
$\varDelta _{\parallel }$, and the lateral spacing of capsules in adjacent trains is denoted as 
$\varDelta _{\bot }$. Note that the 
$x$-axis is in logarithmic scale. The plot indicates decreasing 
$\tilde {u}_c$ with decreasing 
$\varDelta _{\parallel }$ and 
$\varDelta _{\bot }$. (b) Slice in the 
$x$–
$z$ plane bisecting a capsule train used for the velocity plots in (c,d). (c) Induced flow field and (d) comoving flow field for 
$\varDelta _{\parallel }/a = 2.22$ and 
$\varDelta _{\bot }/a = 20$.
It was also observed that the principal tension within the capsule membranes also decreased for capsules within trains, relative to isolated capsules, as shown in figure 15. The maximum principal tension on a capsule surface 
$T_1^{max}$ was found to decrease with decreasing capsule spacing. In these planar channels, capsules assume the characteristic folded shape, and the location of 
$T_1^{max}$ on the capsule surface was observed to be midway between the centreplane and the top (or bottom) of the capsule. The reduction in membrane tension for capsules within trains likely is a result of the decreased velocity shown in figure 13, i.e. capsules travelling at lower velocity experience less deformation. However, this also raises the question as to whether a reduction in the mechanical deformation energy is an additional driving factor for capsule self-assembly (i.e. minimizing system energy), one that acts in addition to hydrodynamic forces. This question will be left for a future study.
(a) Maximum membrane principal tension for capsules in trains versus in-line spacing for varying lateral spacing values. Here, 
$Ca = 0.3$. Capsules in trains exhibit decreasing 
$T_1^{max}$ values with decreasing 
$\varDelta _{\parallel }$. (b) Rear, side and front views of a single capsule (
$\varDelta _{\parallel }/a = \varDelta _{\bot }/a = 20$) coloured by 
$T_1/\kappa _s$. (c) Similar views of capsules in a train (
$\varDelta _{\parallel }/a = 2.22$ and 
$\varDelta _{\bot }/a = 20$). The colouring reveals reduced values of principal membrane tension in (c) relative to (b).
3.5. Dispersity in capsule deformability
Attention is now turned to suspensions with dispersity in the deformability, using two types of distributions: (i) a Gaussian distribution in the 
$Ca$ values throughout the suspension; and (ii) a bimodal distribution whereby two sub-populations of capsules exist, each with a unique value of 
$Ca$. For all simulations in this section, the capsule radii are uniform and the confinement ratio is fixed at 
$a/h = 0.6$.
For suspensions with Gaussian distributions, two parameters are adjusted independently: the mean deformability 
$\overline {Ca}$, and the standard deviation in deformability 
$\sigma _{Ca}$. Figure 16 shows simulation snapshots for three different suspensions, each with the same standard deviation 
$\sigma _{Ca} = 0.1$ but with different mean deformabilities: 
$\overline {Ca} = 0.1$ (figures 16a,d,g), 
$\overline {Ca} = 0.2$ (figures 16b,e,h) and 
$\overline {Ca} = 0.3$ (figures 16c,f,i). Figures 16(a–c) show histograms of the probability distributions for 
$Ca$. A lower bound on 
$Ca$ is enforced (
$Ca_{min} = 0.008$), otherwise the distribution curves would in some cases extend into negative 
$Ca$ values, which is not physical. Inspection of figure 16 reveals that increasing order occurs with increasing 
$\overline {Ca}$, which may be expected for systems with overall higher mean deformability. In figure 16(a), for the suspension with 
$\overline {Ca} = 0.1$, the portion of capsules with 
$Ca$ values below 0.05 is not insignificant. In the results for uniform suspensions (see figure 6), the degree of ordering for 
$Ca < 0.05$ was rather low. Hence the presence of these harder capsules on the left tail of the distribution curve inhibits ordering for the suspension as a whole. On the other hand, in figure 16(c), for the suspension with 
$\overline {Ca} = 0.3$, the lower tail of the 
$Ca$ distribution curve essentially vanishes for 
$Ca < 0.1$, hence there are very few rigid capsules in this suspension, therefore resulting in a higher degree of ordering.
Capsule configurations in systems with Gaussian distributions in 
$Ca$. The confinement is 
$a/h = 0.6$, and the standard deviation is 
$\sigma _{Ca} = 0.1$, with varying values of mean 
$Ca$: (a,d,g) 
$\overline {Ca} = 0.1$, (b,e,h) 
$\overline {Ca} = 0.2$, and (c,f,i) 
$\overline {Ca} = 0.3$. (a–c) Histograms of the probability density for 
$Ca$. (d–f) Top-down views of the suspension at the end of simulation. (g–i) Close-up views of a sub-portion of the domain (indicated by the red boxes in d–f), in which capsules are coloured by the principal membrane tension 
$T_1/\kappa _s$. Harder capsules deform less and therefore have lower values of 
$T_1$, thus they are shaded a darker blue. In these simulations, a floor 
$Ca = 0.008$ was enforced.
Figure 17(a) plots 
$\varPhi _t$ versus 
$\sigma _{Ca}$ for the three values of mean deformability studied. The data show that for a mean deformability 
$\overline {Ca} = 0.3$, essentially no reduction in ordering is observed with increasing dispersity up to the maximum value used, 
$\sigma _{Ca} = 0.1$. For a smaller mean deformability, 
$\overline {Ca} = 0.1$, there is an observed decrease in ordering with increasing 
$\sigma _{Ca}$, particularly when the dispersity is increased beyond 
$\sigma _{Ca} > 0.04$. The statistical analysis for these suspension flows, shown in figures 17(b) and 17(c), reveals that despite the variability in deformability, the capsule velocities are very uniform throughout the suspensions. This is consistent with the earlier result shown in figure 3(a) that capsule velocity is not particularly dependent on deformability for confinement ratio 
$a/h = 0.6$. Finally, it can be seen that softer capsules are slightly more likely to belong to trains than harder capsules.
(a) Order parameter 
$\varPhi _t$ versus standard deviation of 
$Ca$ for three values of 
$\overline {Ca} = 0.1$, 0.2 and 0.3. The confinement is 
$a/h = 0.6$. (b) For 
$\overline {Ca} = 0.1$ and 
$\sigma _{Ca} = 0.04$, the probability density 
$p$, capsule velocity 
$\tilde {u}_c$, and 
$\varPhi _t$ are plotted as histograms versus 
$Ca$. (c) Same as (b) except with 
$\sigma _{Ca} = 0.1$. Capsule velocities are fairly uniform across the deformability range. Softer capsules have a higher probability of belonging to a train.
Finally, a bimodal distribution in deformability is investigated, whereby two sub-populations are assigned one of two values of deformability, 
$Ca_1$ and 
$Ca_2$. The fraction of capsules with 
$Ca_2$ is denoted 
$X_2$. In this study, 
$Ca_2$ is held fixed at value 0.3, and 
$Ca_1$ and 
$X_2$ are varied systematically. Again, the confinement ratio is fixed at 
$a/h = 0.6$. Figure 18 shows simulation snapshots for suspensions with 
$Ca_1 = 0.01$ (relatively rigid) and 
$Ca_2 = 0.3$ (soft) for three different values of 
$X_2$. Increasing the fraction of softer capsules increases the ordering, as expected. Also, figures 18(g–i) show that softer capsules experience greater levels of membrane tension due to higher deformation.
Capsule configurations for bimodal distributions in 
$Ca$ with two populations: 
$Ca_1 = 0.01$ and 
$Ca_2 = 0.3$. The confinement is 
$a/h = 0.6$. The fraction of capsules with 
$Ca_2$ is varied: (a,d,g) 
$X_2 = 0.1$, (b,e,h) 
$X_2 = 0.4$, and (c,f,i) 
$X_2 = 0.8$. A visual inspection indicates that the degree of ordering increases with increasing 
$X_2$, as would be expected. In (g–i), the softer capsules (with deformability 
$Ca_2$) experience higher levels of membrane tension.
For bimodal distributions in deformability, the plot of 
$\varPhi _t$ versus 
$X_2$ is shown in figure 19(a) for three different values of 
$Ca_1$. At the extremes of this plot, 
$X_2 = 0$ and 
$X_2 = 1$, the suspensions are uniform and the data values agree very well with those shown in figure 7. In particular, for 
$X_2 = 1$, the suspension consists of all soft capsules, and the value of 
$\varPhi _t$ is very close to 1. For 
$X_2 = 0$, the values of 
$\varPhi _t$ depend on the value of 
$Ca_1$, with decreasing 
$\varPhi _t$ as 
$Ca_1$ decreases from 0.1 to 0.05 to 0.01.
(a) Order parameter 
$\varPhi _t$ versus 
$X_2$ in systems with bimodal distributions in capsule deformability. Three different sets of 
$Ca$ were considered, as indicated in the legend (
$Ca_2 = 0.3$ in each set). The confinement is fixed at 
$a/h = 0.6$. (b) For 
$Ca_1 = 0.01$, 
$Ca_2 = 0.3$ and 
$X_2 = 0.3$, the probability density 
$p$, capsule velocity 
$\tilde {u}_c$, and 
$\varPhi _t$ are plotted as histograms versus 
$Ca$. (c) Same as (b) except with 
$X_2 = 0.8$. Softer capsules with deformability 
$Ca_2$ are more likely to belong to a train.
For mixtures of hard and soft capsules, the transition from 
$X_2 = 0$ to 
$X_2 = 1$ is monotonic for the case 
$Ca_1 = 0.01$, but not for the cases 
$Ca_1 = 0.05$ and 
$Ca_1 = 0.1$. For 
$Ca_1 = 0.05$ and 
$Ca_1 = 0.1$, it was found that adding 20 % to 30 % of the softer capsules into the suspension actually reduced the ordering compared to a uniform suspension (
$X_2 = 0$) of capsules with less deformability. The curves for 
$Ca_1 = 0.05$ and 
$0.1$ in figure 19(a) loosely resemble the curves shown in figure 13 for bimodal distributions of radii. Here, the data suggest that adding a small quantity of soft capsules into a suspension of medium-soft capsules reduces the ordering more than vice versa (or adding a small quantity of medium-soft capsules to a suspension of soft capsules). The histograms shown in figures 19(b) and 19(c) reveal that similar to the case of a Gaussian distribution in deformability (see figures 17b,c), the velocities of hard and soft capsules are essentially the same within the suspensions, and that softer capsules are more likely to belong to trains than harder capsules.
Overall, the above simulation results suggest that dispersities in capsule radius as well as capsule deformability can impact the hydrodynamic ordering process in planar Poiseuille flow. When comparing the reduction in ordering due to Gaussian dispersity in capsule size versus Gaussian dispersity in capsule deformability, based on the plots of figures 11(a) and 17(a), it appears that dispersity in capsule size leads to a greater reduction in ordering. However, confinement and deformability are both necessary factors leading to the 1-D assembly phenomenon.
4. Conclusions
The organization of deformable particles in confined flow conditions is an important phenomenon in the microcirculation of blood, and a useful mechanism for the manipulation of particle-laden suspensions within microfluidic devices. Here, three-dimensional LBM/IBM simulations were performed to carry out a systematic investigation into the dependency of ordering on two key parameters, capsule deformability and confinement ratio, for a suspension of elastic capsules in planar Poiseuille flow. Furthermore, dispersity in these two parameters was also investigated to understand how ordering may diminish in non-uniform suspensions.
The primary conclusions from this work are listed below.
(i) For uniform suspensions in planar flow, the degree of ordering was found to increase with increasing deformability. The degree of ordering also increases with increasing confinement ratio, but only to a certain point (
$a/h = 0.6$), beyond which the ordering decreases, therefore resulting in an optimal confinement ratio. At confinement ratios $a/h > 0.6$, at fixed volume fraction, the decrease in $\varPhi _t$ is attributed to a decrease in the in-plane density of capsules.(ii) Capsules that are assembled in trains are transported at a lower velocity and experience lower membrane tension compared with isolated capsules (i.e. capsules not in trains).
(iii) Dispersity in the capsule size leads to a sharp decrease in ordering, at standard deviations
$\sigma _{\tilde {a}} < 0.02$ (this is the standard deviation of the confinement ratio). In suspensions with dispersity in capsule size, larger capsules travel at lower velocities, experience higher membrane tension, and are more likely to belong to trains.(iv) Dispersity in deformability can lead to a decrease in ordering, depending on the mean deformability. For higher values of
$\overline {Ca}$ ($= 0.3$, for example), very little change in ordering was observed for increasing dispersity in deformability. For a lower value of $\overline {Ca}$ ($= 0.1$), a decrease in ordering was observed for increasing dispersity, due to an increased percentage of stiffer capsules within the suspension.(v) Overall, the results suggest that ordering is more sensitive to dispersity in capsule size relative to dispersity in deformability.
(vi) A bimodal distribution of capsules consisting of larger and smaller capsules also decreases the ordering behaviour, particularly with increasing discrepancy between the two radii. It was found that the
$\varPhi _t$ curve was asymmetric, so that within a binary mixture, a small quantity of larger capsules disrupts the ordering more than a small quantity of smaller capsules.
Finally, this study was limited to a single Reynolds number and particle volume fraction. Further work is needed to understand the role of inertia and concentration on the ordering behaviour. In addition, the capsule deformability is characterized only in terms of the in-plane shear modulus (
$\kappa _s$). Varying the in-plane area dilation modulus (
$\kappa _a$; see Aouane, Scagliarini & Harting Reference Aouane, Scagliarini and Harting2021) or the out-of-plane bending modulus, or using a different membrane mechanics model altogether, may lead to different ordering behaviour.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2023.1052.
Funding
The author gratefully acknowledges financial support from the 21st Century Endowed Professorship provided by the University of Arkansas. This research is also supported by the Arkansas High Performance Computing Center, which is funded through multiple National Science Foundation grants and the Arkansas Economic Development Commission.
Declaration of interest
The author reports no conflict of interest.
Appendix A
Within the planar channels, once a monolayer of capsules has developed, the planar density (
$PD$) can be calculated as
\begin{equation} PD = \frac{N}{L_x L_y}, \end{equation}
where 
$N$ is the number of capsules, and 
$L_x$ and 
$L_y$ are the dimensions of the box in the 
$x$- and 
$y$-directions, respectively. The number of capsules is a function of the volume fraction 
$\phi$:
\begin{equation} N = \frac{\phi L_xL_yL_z}{4/3{\rm \pi} a^3}. \end{equation}
Hence, using the equivalency 
$L_z = 2h$, the planar density in units of 
$a^{-2}$ is
\begin{equation} PD = \frac{3\phi}{2{\rm \pi}(a/h)}. \end{equation}
One can also approximate the average spacing between capsules, assuming that capsules are arranged in a square lattice with uniform spacing 
$\varDelta$ (in units of 
$a$):
\begin{equation} \varDelta = \sqrt{\frac{1}{PD}} = \sqrt{\frac{2{\rm \pi}(a/h)}{3\phi}}. \end{equation}
