1. Introduction
Capsules, membrane-enclosing fluid objects, play pivotal roles across many natural and industrial domains, including the biological, pharmaceutical, cosmetic and food industries (Barthes-Biesel Reference Barthes-Biesel2016; Luo & Bai Reference Luo and Bai2019; Sun et al. Reference Sun, Fu, Peng, Jiao, Liu, Du and Zhang2021). Capsules also serve as simplified models for examining complex biological cells, such as red blood cells (RBCs) and circulating tumour cells, through numerical analysis (Krüger Reference Krüger2012; Gekle Reference Gekle2016; Balogh et al. Reference Balogh, Gounley, Roychowdhury and Randles2021; Huet & Wachs Reference Huet and Wachs2023). Given their significant role and applications, understanding and manipulating capsules require advanced theoretical, experimental and numerical techniques. Recent advances have highlighted the utility of microfluidic devices in separating cells by size and deformability, enriching specific cell concentrations, and facilitating precise cell characterisation (Häner et al. Reference Häner, Vesperini, Salsac, Le Goff and Juel2021; Wang et al. Reference Wang, Sui, Salsac, Barthès-Biesel and Wang2016). Enhanced computing capabilities have further bolstered the design and optimisation of these devices, enabling sophisticated numerical analyses that mimic experimental conditions (Zhu et al. Reference Zhu, Rorai, Mitra and Brandt2014; Gubspun et al. Reference Gubspun, Gires, de Loubens and Barthes-Biesel2016). Notably, innovative microfluidic designs, such as a semi-circular pillar geometry developed and optimised numerically by Zhu et al. (Reference Zhu, Rorai, Mitra and Brandt2014), have demonstrated efficacy in sorting cells based on membrane stiffness, confirmed through both numerical studies and experimental validations (Häner et al. Reference Häner, Vesperini, Salsac, Le Goff and Juel2021).
Over the last years, considerable research efforts have been directed towards modelling and understanding the deformation and dynamics of capsules in non-inertial flows, a focus driven largely by its relevance to microcirculatory systems such as capillary blood vessels and conventional microfluidic applications. Based on the thin-shell theory, Barthès-Biesel and Rallison initially provided an analytical framework for the temporal evolution of elastic capsule deformation under infinite, slow-moving shear flows, assuming infinitesimal deformations (Barthes-Biesel & Rallison Reference Barthes-Biesel and Rallison1981). Subsequently, Pozrikidis expanded on this by applying the boundary integral method to accommodate finite deformations (Pozrikidis Reference Pozrikidis1995), a method further applied to analyse significant deformation in capsules with distinct internal and external fluid viscosities (Ramanujan & Pozrikidis Reference Ramanujan and Pozrikidis1998), and to incorporate the effects of bending stresses (Pozrikidis Reference Pozrikidis2001). In terms of the capsule dynamics, it was shown that the migration and diffusion of a capsule suspension do differ from those of hard sphere suspensions due to the deformable capsule–capsule interactions and the dissipating inner fluid (Bishop et al. Reference Bishop, Popel, Intaglietta and Johnson2002; Kulkarni & Morris Reference Kulkarni and Morris2008; Gross et al. Reference Gross, Krüger and Varnik2015; Rahmani et al. Reference Rahmani, Hammouti and Wachs2018). More recently, Balogh and Bagchi (Balogh & Bagchi Reference Balogh and Bagchi2017, Reference Balogh and Bagchi2018, Reference Balogh and Bagchi2019) employed the front-tracking method (FTM) to delve into the behaviour and structural changes of RBCs within the complex vessel networks reminiscent of human microcirculation. Apart from the migration of capsules in channels of different geometries, the mechanical and rheological properties of suspensions of purely elastic capsules have been thoroughly studied in the non-inertial regime. The rheological behaviour of suspensions, especially for the rigid particles, was addressed in the pioneering work of Batchelor (Reference Batchelor1970), establishing that stresses due to the presence of particles are formulated using a particle stress tensor, which can be expressed as a summation of stresslets in a domain. Recently, Takeishi et al. (Reference Takeishi, Rosti, Imai, Wada and Brandt2019) conducted a comprehensive investigation of the rheology of a suspension of RBCs in a wall-bounded shear flow, highlighting that the dynamic motion of RBCs has an important influence on the suspension rheology. Aouane et al. (Reference Aouane, Scagliarini and Harting2021) reveal that the strain-hardening capsules share rheological features with both soft and solid particles depending on the ratio of the area dilatation to shear elastic modulus and the effective volume fraction. Guglietta et al. (Reference Guglietta, Pelusi, Sega, Aouane and Harting2023) studied the suspension of viscoelastic capsules by considering the effects of membrane viscosity on the transient dynamics of the suspension rheology. They showed that the influence of the membrane viscosity strongly depends on the capsule volume fraction.
Inertial particle microfluidics has recently emerged as a cutting-edge technique for high-precision and efficient manipulation of microparticles and biological cells. This technology finds application in fields ranging from microbiology to biotechnology, offering a promising avenue to expedite labour-intensive, non-destructive diagnostic processes such as cell sorting and cell characterisation (Krüger et al. 2014; Son et al. Reference Son, Samuel, Gale, Carrell and Hotaling2017; Yin et al. Reference Yin, Wang, Li, Lin, Shao, Cao and Hou2019; Owen et al. Reference Owen2023). Given the substantially higher flow velocities in inertial microfluidics compared with their non-inertial counterparts, traditional analytical theories and methods relying on the Stokes conditions such as the boundary integral method fail to account for the advective components of the fluid momentum equation (Amini et al. Reference Amini, Lee and Di Carlo2014). Initial investigations into the inertial dynamics of capsules, such as the work of Doddi & Bagchi (Reference Doddi and Bagchi2008), who explored the interactions and spiralling motions of two capsules in a shear flow using the FTM, laid the groundwork for this field. Research into the dynamics of deformable capsules in inertial flows has revealed distinct behaviours depending on channel geometry. Subsequent studies in straight microchannels identified multiple equilibrium positions for deformable capsules, diverging from the channel centreline towards the cross-section diagonals, and primarily influenced by inertial forces: wall repulsion, shear-gradient lift and drag forces in secondary flow fields (Kilimnik et al. Reference Kilimnik, Mao and Alexeev2011; Raffiee et al. Reference Raffiee, Dabiri and Ardekani2017; Owen et al. Reference Owen2023). Building on this, Ebrahimi & Bagchi (Reference Ebrahimi and Bagchi2021) extended the analysis to curved channels, systematically varying parameters such as channel Reynolds number, capsule deformability and channel geometry to uncover the effects of curvature on the capsule dynamics. Parallel advancements have focused on other deformable particles, notably droplets, whose behaviours under finite inertia have been extensively explored. These studies highlight insights into the dynamics of emulsions, emphasising the effects of inclination angles on rheological behaviours (Li & Sarkar Reference Li and Sarkar2005b ; Olapade et al. Reference Olapade, Singh and Sarkar2009; Singh & Sarkar Reference Singh and Sarkar2009; Singh & Sarkar Reference Singh and Sarkar2011, Reference Singh and Sarkar2015). For elastic capsule suspensions specifically, Banaei et al. examined nucleated capsules under finite inertia (Banaei et al. Reference Banaei, Loiseau, Lashgari and Brandt2017, Reference Banaei, Shahmardi and Brandt2021), showing significant effects on deformation and suspension rheology, including increased relative viscosity. In our previous work (Huet et al. Reference Huet, Morente, Gai and Wachs2024), we further investigated the dynamics of elastic capsule trains navigating sharp corners under inertial flow, identifying mechanisms behind interaction and deformation in these complex geometries.
Despite the existing body of knowledge previously summarised, significant gaps still remain in our understanding of the deformation, diffusion and rheology of elastic capsule suspensions in an inertial shear flow. Current research has yet to comprehensively address the intricate interactions and behaviours that arise under varying conditions of inertia and flow dynamics. We aim to answer the following questions using high-fidelity particle-resolved numerical simulations: (i) How does flow inertia impact the deformation and dynamics of a single capsule in a shear flow? (ii) What are the connections between the behaviour of a single capsule and that of a dense suspension? (iii) What are the major differences between capsule suspensions in the non-inertial and inertial regimes?
This manuscript is organised as follows: § 2 provides the governing equations and membrane models for the deformable capsule. In § 3, we introduce the numerical set-ups for a single capsule and a capsule suspension in a shear flow. The novel findings for a single capsule in a shear flow are presented in § 4. In § 5, we present the numerical results on the deformation, lateral migration and rheology of the capsule suspension and discuss their relationship with the behaviour observed in single capsule studies. Finally, § 6 provides the conclusions and perspectives for future investigation.
2. Numerical method and dimensionless numbers
2.1. Governing equation and membrane models
The capsule membrane
$\tilde {\Gamma }$
is assumed infinitely thin and is surrounded by an incompressible Newtonian fluid with constant viscosity and density. In this study, the viscosity ratio between the inner and outer fluids of the capsule may differ from unity. The fluid dynamics is governed by the mass and momentum conservation equations
\begin{align} &\tilde {\nabla } \cdot \boldsymbol {\tilde {u}} = 0, \end{align}
\begin{align} \frac {\partial \boldsymbol {\tilde {u}}}{\partial \tilde {t}} + \boldsymbol {\tilde {u}} \cdot \tilde {\nabla } \boldsymbol {\tilde {u}} &= \frac {1}{\tilde {\rho }} \tilde {\nabla } \tilde {p} + \tilde {\nu } \tilde {\nabla }^2 \boldsymbol {\tilde {u}} + \frac {1}{\tilde {\rho }} \boldsymbol {\tilde {f_b}}, \end{align}
where
$\boldsymbol {\tilde {u}}$
denotes the velocity field,
$\tilde {p}$
the pressure field and
$\tilde {\rho }$
the density. The kinematic viscosity is given by
$\tilde {\nu } = \tilde {\mu }/\tilde {\rho }$
, where
$\tilde {\mu }$
is the dynamic viscosity. The term
$\boldsymbol {\tilde {f_b}}$
represents a body force that accounts for the influence of the capsule membrane on its surrounding fluid. Dimensional quantities are indicated by the
$\sim$
symbol. The membrane can exhibit both elasticity and bending stresses, and its localised action on the fluid is described by the following expression for
$\boldsymbol {\tilde {f_b}}$
:
\begin{equation} \boldsymbol {\tilde {f_b}} = \left ( \boldsymbol {\tilde {F}_{{elastic}}} + \boldsymbol {\tilde {F}_{{bending}}} \right ) \tilde {\delta }(\boldsymbol {\tilde {x}} - \boldsymbol {\tilde {x}_{0}}), \end{equation}
where
$\tilde {{\boldsymbol {x}}} = (\tilde {x}, \tilde {y}, \tilde {z})$
is the position vector and
$\tilde {\delta }(\boldsymbol {\tilde {x}} - \boldsymbol {\tilde {x}_{0}})$
is a Dirac delta distribution defined at position
$\boldsymbol {\tilde {x}_{0}}$
, which is non-zero only on the surface of the membrane since
$\boldsymbol {\tilde {x}_{0}}\in \tilde {\Gamma }$
. Please note that only the elastic force is considered in the current research, while bending stresses are neglected.
The shear and area-dilatation membrane stresses are described using the thin-shell theory, which is briefly summarised here. For more detailed information, readers are referred to the analytical study by Barthes-Biesel & Rallison (Reference Barthes-Biesel and Rallison1981). In this work, we adopt the Skalak law to describe the elastic membrane, with the expression of the surface strain energy function
$\tilde {W_s}$
expressed as
\begin{equation} \tilde {W_s} = \frac {\tilde {E_s}}{4} \left ( I_1^2 + 2I_1 - 2I_2 + CI_2^2\right ). \end{equation}
The two invariants are defined as
$I_1 = \lambda _1^2 + \lambda _2^2 - 2$
and
$I_2 = \lambda _1^2\lambda _2^2 - 1$
, where
$\lambda _{1,2}$
are the principal stretches in the two tangential directions, and
$\tilde {E_s}$
is the shear modulus. The constant
$C = 1$
represents the area-dilatation modulus, which prevents significant area changes (Pozrikidis Reference Pozrikidis1995). Subsequently, the membrane elastic stresses
$\tilde {\sigma }_{i,j}$
are calculated using the surface strain energy function
$\tilde {W_s}$
\begin{equation} \tilde {\sigma }_{i,j} = \frac {1}{\lambda _j} \frac {\partial \tilde {W_s}}{\partial \lambda _i}, \qquad i, j \in \{1,2\}, \quad i \neq j. \end{equation}
Once the elastic stress is determined, the elastic force exerted by the membrane on the fluid is calculated as follows:
\begin{equation} \boldsymbol {\tilde {F}_{{elastic}}} = \tilde {\nabla } \cdot \boldsymbol {\tilde {\sigma }}. \end{equation}
The adaptive FTM is employed to solve the aforementioned equations. A concise overview of the numerical approach is provided in Appendix A. Equations (2.1) and (2.2) are spatially discretised using the finite volume method on an adaptive octree grid, or on a homogeneous Cartesian grid for capsule suspensions; and are implemented within the open-source software Basilisk (Popinet Reference Popinet2015). A comprehensive description of the numerical implementation is available in Huet & Wachs (Reference Huet and Wachs2023).
2.2. Physical quantities and dimensionless numbers
Our problem is governed by several dimensionless numbers that describe and control the flow configuration. These include the Reynolds number
$\mathcal {R}e$
, which measures the ratio of inertial to viscous forces; the capillary number
$\mathcal {C}a$
, representing the ratio of viscous to elastic stresses; the capsule volume fraction
$\phi$
, indicating the ratio of volume occupied by the capsules; and the fluid viscosity ratio
$\lambda$
of the internal and external fluids of the capsules. Please note that we adopt
$\tilde {r}_0$
and
$\tilde {\dot {\gamma }}\tilde {r}_0$
as length scale and velocity scale, respectively. For brevity, we summarise the key input parameters of the flow configuration in table 1, which presents the mathematical expression and range of values of each parameter, together with associated symbols and the corresponding descriptions.
Flow control parameters for numerical simulations: mathematical expressions, range of values and detailed symbol descriptions.
As for the output parameters, we investigate various physical quantities of capsules in an inertial shear flow, including the Taylor deformation factor (
$D$
), particle stress (
$\Sigma$
), relative viscosity (
$\mu _r$
) and capsule diffusivity (
$D_\psi$
). Below, we briefly outline the methods used to compute and normalise these physical quantities.
The information on the deformation of the capsules is retrieved from the moment of inertia tensor. We consider the ellipsoid with the same inertia tensor
$\tilde {\mathcal {I}}$
as the capsule. Once we obtain the three eigenvalues (
$\tilde {\zeta }_1 \lt \tilde {\zeta }_2 \lt \tilde {\zeta }_3$
) of
$\tilde {\mathcal {I}}$
, we can use them to compute the length of the three semi-axes :
$\tilde {r}_1, \tilde {r}_2$
and
$\tilde {r}_3$
(Ramanujan & Pozrikidis Reference Ramanujan and Pozrikidis1998; Guglietta et al. Reference Guglietta, Pelusi, Sega, Aouane and Harting2023). Here,
$\tilde {r}_1$
and
$\tilde {r}_3$
are the longest and shortest lengths of the semi-axis in the shear plane (
$x{-}y$
plane), and
$\tilde {r}_2$
denotes the radius along the vorticity direction along the axis normal to the shear plane (
$z$
-axis). Subsequently, the Taylor deformation factor is defined as
\begin{equation} D = \frac {\tilde {r}_1-\tilde {r}_3}{\tilde {r}_1+\tilde {r}_3} . \end{equation}
Sequential illustration of a spherical capsule’s deformation in the shear plane. Transitioning from blue to red depicts the deformation dynamics of a spherical capsule. Here,
$r_1$
,
$r_3$
and the inclination angle
$\theta$
are annotated for clarity.
Figure 1 presents a brief sketch illustrating the two main semi-axes
$r_1$
,
$r_3$
and the inclination angle
$\theta$
of the capsule in shear flow on a
$x{-}y$
cut plane crossing the capsule centroid. The transition of the line colours from blue to red depicts the deformation process of the capsule, transitioning from an initial spherical shape to an oval shape as it reaches the steady state. The inclination angle is defined as the angle between the longest semi-axis and the streamwise direction along the
$x$
-axis and is computed using the eigenvectors, as depicted in figure 1
\begin{equation} \theta = arctan\left(\frac {\tilde {v}_y}{\tilde {v}_x}\right). \end{equation}
The bulk excess stress of the capsule suspension can be estimated using Batchelor’s formulation (Batchelor Reference Batchelor1970), considering two contributions to the total excess stress:
$\Sigma ^p$
and
$\Sigma ^f$
. The particle stress tensor, which depends on the interaction between capsules
$\Sigma ^p$
, can be computed by incorporating the viscosity ratio and the Lagrangian nodal forces acting on the capsule surface, as described in the methodologies by Krüger (Reference Krüger2012) and Gross et al. (Reference Gross, Krüger and Varnik2014)
\begin{equation} \tilde {\Sigma }_{\alpha \beta }^p = -\frac {1}{\tilde {V}} \sum _{i = 1}^{N}\sum _{j}^{N_{ele}}\frac {1}{2}\!\left( \!\tilde {F}_\alpha ^{i,j}\tilde {r}_\beta ^{i,j} + \tilde {F}_\beta ^{i,j}\tilde {r}_\alpha ^{i,j}\right) + \frac {1}{\tilde {V}} \sum _{i = 1}^{N}\oint _{\tilde {\mathcal {A}}_i}\!\! \left(\tilde {\lambda }_v - 1\!\right)\tilde {\mu }_0\left(\tilde {u}_\alpha n_\beta + \tilde {u}_\beta n_\alpha \right) {\rm d}\mathcal {A}, \end{equation}
where
$\alpha$
and
$\beta$
are indices referring to the Cartesian directions (
$x$
,
$y$
or
$z$
),
$N$
denotes the number of capsules considered in the control volume
$\tilde {V}$
and
$N_{ele}$
is the number of triangular elements on the capsule;
$\tilde {\mathcal {A}}_i$
is the area of the capsule
$i$
and
$\tilde {\mathcal {A}}_0$
is the surface area of the initial spherical capsule;
$\tilde {F}^{i,j}$
is Lagrangian force on the node with position vector
$\tilde {\boldsymbol {r}}^{i,j}$
. In the inertial regime, an additional Reynolds-stress-like term, denoted as
$\Sigma ^{f}$
, needs to be considered to account for the contribution of flow fluctuations to the bulk stress. For a dilute system, the following formulation can be applied (Raja et al. Reference Raja, Subramanian and Koch2010; Srivastava et al. Reference Srivastava, Malipeddi and Sarkar2016; Mwasame et al. Reference Mwasame, Wagner and Beris2017):
\begin{equation} \tilde {\Sigma }^f_{\alpha \beta } = -\frac {1}{\tilde {V}} \int _{\tilde {V}} \tilde {\rho } \tilde {u}'_{\alpha } \tilde {u}'_{\beta } , {\rm d}V, \end{equation}
where
$\tilde {u}' = \tilde {u} - \overline {\tilde {u}}$
is the flow fluctuation relative to the reference (or undisturbed) velocity of the imposed shear flow
$\overline {\tilde {u}} = \tilde {\dot {\gamma }}\tilde {y}$
.
The rheology of the suspension is assessed using the relative viscosity and the normal stress differences. For example, the relative viscosity of the suspension is defined as
\begin{equation} \mu _r = 1 + \Sigma _{xy}^{tot} = 1 + \frac {\tilde {\Sigma }_{xy}^p + \tilde {\Sigma }_{xy}^f}{\tilde {\mu }_0{\tilde {\dot {\gamma }}}}, \end{equation}
where
${\Sigma }_{xy}^{tot}$
denotes the dimensionless total excess particle shear stress. The dimensionless first and second normal stress differences are defined as
\begin{equation} N_1^* = \frac {\tilde {N}_1^*}{\tilde {\mu }_0\tilde {\dot {\gamma }}} = \frac {\tilde {\Sigma }_{xx}^* - \tilde {\Sigma }_{yy}^*}{\tilde {\mu }_0\tilde {\dot {\gamma }}}, \quad N_2^* = \frac {\tilde {N}_2^*}{\tilde {\mu }_0\tilde {\dot {\gamma }}} = \frac {\tilde {\Sigma }_{yy}^* - \tilde {\Sigma }_{zz}^*}{\tilde {\mu }_0\tilde {\dot {\gamma }}}, \end{equation}
where * indicates the contribution to the particle stress that arises from viscosity ratio and particle interaction (
$p$
) or from the flow fluctuation (
$f$
). Similar approaches are applied to non-dimensionalise other quantities using components of the stress tensor.
In addition to suspension rheology, we also investigate particle lateral migration by analysing the statistical motion of the capsules, which can be described using the mean-square displacement (MSD)
\begin{equation} \big\langle \unicode{x1D6E5} _{\psi }(t)^2 \big\rangle := \left\langle (\tilde {\boldsymbol {X}}_{i,\psi }\left(\tilde {t} + \unicode{x1D6E5} \tilde {t}) - \tilde {\boldsymbol {X}}_{i,\psi }(\tilde {t})\right)^2 \right\rangle _{i,t}, \end{equation}
where
$\boldsymbol {X}$
denotes the position of the capsule centroid and
$\psi$
is the diffusion direction of the capsule centroid, i.e. the velocity gradient direction (
$y$
-axis) or the vorticity direction (
$z$
-axis) in the shear flow. The ensemble average is taken over time and over all freely moving capsules in the bulk flow, ensuring that both
$\tilde {t}$
and
$\tilde {t} + \Delta \tilde {t}$
are within the steady-state interval. Then, the diffusivity along
$\psi$
direction can be defined as the slope of the temporal evolution of the MSD in the diffusion regime. Conventionally, the dimensionless diffusivity of the capsules is scaled as
\begin{equation} D_\psi = \frac {\tilde {D}_\psi }{\tilde {\dot {\gamma }} \tilde {r}_0^2}, \,\,\, \mbox {with} \,\,\, \tilde {D}_\psi = \frac {\left\langle \unicode{x1D6E5} _{\psi }(t)^2 \right\rangle }{\tilde {d}_0^2\tilde {\dot {\gamma }}\tilde {t}}. \end{equation}
It is important to note that
$\Delta \tilde {t}$
here represents a time window used for the computation of the MSD and
$D_\psi$
. This time window is distinct from, and can be significantly larger than, the numerical time step.
3. Numerical set-up
3.1. Numerical set-up of a single capsule
Illustration of a single capsule within an octree grid in a three-dimensional cubic computational domain
$\tilde {\mathcal {D}}$
. The red box shows a zoomed view of the capsule at
$\phi = 8\times 10^{-4}$
in simple shear flow. The blue box indicates the position of the capsule with the adaptive grid in the
$x{-}y$
cut plane at
$z=0$
.
As depicted in figure 2, we consider a simple shear flow along the
$x$
axis characterised by a reference velocity
$\tilde {U}_0$
within a cubic domain
$\tilde {\mathcal {D}}$
with an edge length
$\tilde {L}$
and centred at
$(0,0,0)$
. The boundaries of the cubic computational domain are referred to as left and right in the
$x$
direction, top and bottom in the
$y$
direction and front and back in the
$z$
direction, collectively denoted as
$\partial \tilde {\mathcal {D}} = \text {left} \cup \text {right} \cup \text {top} \cup \text {bottom} \cup \text {front} \cup \text {back}$
. The flow velocity
$\boldsymbol {\tilde {u}}$
satisfies homogeneous Dirichlet boundary conditions on the top and bottom boundaries, and periodic boundary conditions on the left/right (along
$x$
-axis) and front/back (along
$z$
-axis) boundaries. The complete set of boundary and initial conditions is as follows:
\begin{align} &\quad\tilde {{\boldsymbol {u}}}(\tilde {{\boldsymbol {x}}}, \tilde {t}) = \left ( \frac { \tilde {U}_0}{2}, 0, 0 \right ) \,\,\, \text {on the top boundary}, \end{align}
\begin{align} &\tilde {{\boldsymbol {u}}}(\tilde {{\boldsymbol {x}}}, \tilde {t}) = \left ( - \frac {\tilde {U}_0}{2}, 0, 0 \right ) \,\,\, \text {on the bottom boundary} ,\end{align}
\begin{align} &\qquad\qquad\tilde {{\boldsymbol {u}}}(\tilde {{\boldsymbol {x}}}, 0) = \left ( \frac {\tilde {U}_0}{\tilde {L}} \tilde {y}, 0, 0 \right ) \,\,\, \text {in} \,\,\, \tilde {\mathcal {D}}.\end{align}
We begin with the numerical simulation of a single elastic capsule at
$0.01 \leqslant {\mathcal {C}a} \leqslant 1.0$
. The capsule is strategically positioned at the centre of the cubic domain
$\tilde {\mathcal {D}}$
, as shown in figure 2. This set-up provides a sufficiently large domain for the accurate resolution of the flow dynamics and flow–capsule interactions. Initially at rest, the capsule at
$t=0$
is subject to a simple shear flow in both the quasi-non- to low inertial regime (
${\mathcal {R}e} = 0.1, 1$
), hereafter called the non-inertial regime for simplicity, and the inertial regime (
${\mathcal {R}e} = 10, 20$
). Several viscosity ratios are considered:
$0.1 \leqslant \lambda \leqslant 10$
across the capsule membrane.
To achieve better accuracy in modelling capsule deformation while maintaining reasonable computational cost, we apply an octree adaptive mesh refinement. This strategy, implemented in the open-source platform Basilisk (Popinet Reference Popinet2015), involves partitioning a parent cube cell into eight sub-cubes for local mesh refinement in specific regions of interest (Huet & Wachs Reference Huet and Wachs2023). At each time step, the Cartesian octree grid is dynamically adapted: it is refined in regions with strong gradient variations in any field of interest and coarsened in regions with weak gradient variations. The primary field of interest in this study is the flow velocity. We ensure that the region near the capsule surface always has the finest grid resolution, guaranteeing that the stencils of the Lagrangian points on the capsule membrane are correctly resolved, which is crucial for the effective functioning of the immersed boundary method (IBM) (Appendix A). The hierarchical grid is constructed so that the cell size between two successive levels differs by a factor of two. Consequently, the smallest cell size is
$\tilde {\unicode{x1D6E5} } = \tilde {L}/2^{n_E}$
, where
$n_E$
denotes the maximum refinement level of the Eulerian octree grid. In figure 2, the locally refined grid around the capsule is highlighted. The blue box shows a cut plane at
$z=0$
across the capsule centroid, revealing the instantaneous particle position in the grid and illustrating the size difference between the particle and the computational domain. We conducted
$128$
simulations of a single capsule in shear flow, varying
$\mathcal {R}e$
,
$\mathcal {C}a$
and
$\lambda$
as depicted in table 1. Each simulation was run until a steady state of deformation and dynamics was reached, for a minimum duration.
3.2. Flow configuration of a capsule suspension
Illustration of the capsule suspension in a simple shear flow with increasing volume fractions (a)
$\phi =0.1$
, (b)
$\phi =0.4$
in the three-dimensional cubic computational domain
$\tilde {\mathcal {D}}$
with the upper and lower wall along the
$y$
-axis depicted. The periodic boundary conditions are applied along the
$x$
and
$z$
directions.
Based on the knowledge gained from simulating a single capsule in shear flow, we extend our simulations to investigate capsule suspensions with volume fractions ranging from
$\phi = 0.1$
to
$\phi = 0.4$
. Again, both the non-inertial regime (
${\mathcal {R}e} = 0.1, 1$
) and inertial regime (
${\mathcal {R}e} = 10, 20$
) are considered. Figure 3 depicts the capsule suspension in a simple shear flow at
${\mathcal {R}e} = 1$
and
${\mathcal {C}a}=0.1$
. To ensure consistent resolution across all simulations, all capsules have the same diameter and resolution as in the single capsule case, with a number of points per diameter of the initial spherical capsule
$2\tilde {r_0}/\tilde {\unicode{x1D6E5} } = 30$
. The increase in volume fraction is achieved by adding more capsules to the computational domain. Please note that in the simulation of capsule suspensions, the high number of capsules in the computational domain can significantly increase the computational cost, especially when using octree adaptive mesh refinement. This is due to the extensive refinement required around each capsule at every time step. To manage this computational expense effectively, we employ a uniform Cartesian grid for the simulation of all capsule suspensions presented in this work. This approach balances the need for spatial resolution and computational efficiency, ensuring that the simulations remain feasible without compromising accuracy. Figures 3(a) and 3(b) show capsule suspensions at
$\phi = 0.1$
and
$\phi = 0.4$
, respectively. They clearly illustrate that, as the capsule volume fraction increases, the frequency of strong short-range capsule interactions (in the sense of capsules being in close proximity) rises markedly. To gain a comprehensive understanding of the dynamics, rheology and diffusion of the capsule suspension, we conducted more than 355 numerical simulations with varying
$\mathcal {R}e$
,
$\mathcal {C}a$
,
$\phi$
and
$\lambda$
. Each simulation utilised 64 processors and ran for approximately 21 days to achieve a steady state of average deformation and dynamics, extending up to
$t=200$
. For capsule diffusion, specific cases were designed with smaller time steps and longer simulation times to accurately capture the diffusion coefficient and identify the critical
$\mathcal {R}e$
for regime transitions.
We thoroughly validated our solver by examining the deformation factor, semi-axis lengths, particle shear stress and the first/second normal stress differences. All our results are consistent with data reported in the literature (Bagchi & Kalluri Reference Bagchi and Kalluri2010; Aouane et al. Reference Aouane, Scagliarini and Harting2021; Guglietta et al. Reference Guglietta, Pelusi, Sega, Aouane and Harting2023). Notably, a mesh refinement test in figure 18(g–h) indicates that our choice of smaller grid size enhances the accuracy of numerical results compared with existing literature. For brevity, the detailed numerical validation of our solver for a single capsule and capsule suspensions in the shear flow is provided in Appendix B. In line with methodologies used in previous studies (Aouane et al. Reference Aouane, Scagliarini and Harting2021; Guglietta et al. Reference Guglietta, Pelusi, Sega, Aouane and Harting2023), a repulsive force is also introduced when the Lagrangian nodes of adjacent capsules are in close proximity. For more details, please refer to Appendix C.
4. Single capsule results
In this section, we delve into the effects of flow inertia on the deformation, dynamics and particle stress of a single elastic capsule in the simple shear flow. We aim to provide a comprehensive understanding of how flow inertia influences capsule behaviour in different flow conditions and to shed light on the interplay between these parameters.
4.1. Inertial effects on capsule deformation
Deformation of capsules in the simple shear flow on the
$x{-}y$
plane across the capsule centroid in the cases (a)
${\mathcal {R}e}=1$
and (b)
${\mathcal {R}e}=10$
. The capsules are coloured with the maximum elastic stress
$\sigma _{2,max}$
on its surface when the flow is fully developed. Capsule outlines on the cut plane
$x{-}y$
across the capsule centroid highlighted with pink dot showing (c–d) effects of
$\mathcal {R}e$
on the capsules at
$\lambda =0.1$
and
$1$
; (e–f) effects of
$\lambda$
on the capsules at
${\mathcal {R}e}=0.1$
and
$10$
.
We first examine the deformation of a single capsule within the shear plane (
$x{-}y$
plane) across the capsule centroid, as depicted in figure 4. Understanding the capsule deformation is critical for determining its mechanical response in the shear flow. The maximum stress experienced by the capsule is of particular interest to experimentalists, as it can be used to predict a priori whether a capsule will undergo plastic deformation or even break up in the flow (Häner et al. Reference Häner, Vesperini, Salsac, Le Goff and Juel2021). Specifically, the largest eigenvalue of the membrane elastic stress tensor
${\sigma }_{2}$
provides crucial insight into the membrane integrity.
In figures 4(a) and 4(b), we elucidate the impact of flow inertia (
$\mathcal {R}e$
), membrane elasticity (
$\mathcal {C}a$
) and viscosity ratio (
$\lambda$
) on capsule deformation and the maximum stress distribution. To facilitate comparison, we present the outlines of the most deformable capsule (
${\mathcal {C}a}=1$
) projected on the
$x{-}y$
cut plane across the capsule centroid in figure 4(c–f). At
${\mathcal {C}a} \leqslant 0.1$
, the deformation due to high inertia is not markedly visible. However, a significant increase of the maximum stress,
$\sigma _{2,max}$
with
$\mathcal {R}e$
, is observed. For all capsules featured in figure 4,
$\sigma _{2,max}$
consistently manifests at the centre of the membrane on the
$x{-}y$
plane, at the capsule extreme along the
$z$
-axis. An increase in
$\mathcal {C}a$
results in a relaxation of
$\sigma _{2,max}$
as the capsule membrane becomes more deformable.
The capsule exhibits clearly increased elongation when flow inertia is increased in figure 4(c–d). Conversely, the increase of
$\lambda$
acts to preserve the spherical shape of the capsule as a result of the internal resistance to deformation, especially at
${\mathcal {R}e} \leqslant 1$
. Illustrated in figures 4(a) and 4(e), the capsule at
${\mathcal {C}a}=1$
and
$\lambda =0.1$
adopts a spindle shape, gradually transitioning to a more blunt and oval shape with increasing
$\lambda$
, eventually achieving an inclined oblate spheroid shape at
$\lambda =10$
. In the inertial regime, depicted at
${\mathcal {R}e}=10$
in figures 4(b) and 4(f), the influence of
$\lambda$
on capsule morphology becomes notably less pronounced, where the outlines of the capsules tend to overlap.
4.1.1. Capsule deformation and inclination angle
Effects of flow inertia on the deformation of a single capsule in a shear flow at
${\mathcal {R}e} = 0.1 \sim 20$
, and
$\lambda = 0.1 \sim 10$
. (a–c) Evolution of the Taylor deformation factor
$D$
as a function of
$\mathcal {C}a$
. (d–f) Evolution of the inclination angle
$\theta$
in the
$x{-}y$
shear plane as a function of
$\mathcal {C}a$
, with black dashed lines denoting
$\theta =\pi /4$
. (g–i) Capsule angular velocity of tank treading,
$\omega$
.
Figure 5 shows the impact of flow inertia and viscosity ratio on capsule deformation as a function of
$\mathcal {C}a$
. As depicted in figure 5(a), one immediate effect of inertia is the enhancement of the Taylor deformation factor
$D$
, which increases with
$\mathcal {R}e$
across the range of
$\mathcal {C}a$
values explored. For capsules at
${\mathcal {C}a} \geqslant 0.1$
, elevating the flow inertia from
${\mathcal {R}e}=0.1$
to
${\mathcal {R}e}=20$
leads to a consistent increase in the deformation factor, with an increment of
$\unicode{x1D6E5} D \approx 0.15$
. From another point of view, we fix
$\mathcal {R}e$
to reveal the effects of
$\lambda$
on the capsule deformation in the non-inertial flow at
${\mathcal {R}e}=1$
in figure 5(b) and the inertial flow at
${\mathcal {R}e}=10$
in figure 5(c). In the non-inertial regime, the resistance to deformation becomes more pronounced as
$\mathcal {C}a$
increases. For example, at
${\mathcal {R}e}=1$
and
${\mathcal {C}a}=1$
, the deformation factor drops from
$D=0.54$
at
$\lambda =0.1$
to
$D=0.17$
at
$\lambda =10$
in figure 5(b). However, the impact of
$\lambda$
is much less pronounced in the inertial regime in figure 5(c). A plausible explanation is that in the non-inertial regime (low
$\mathcal {R}e$
), the flow is dominated by viscous forces. At high
$\mathcal {R}e$
, the inertial forces exert such a strong influence on the capsule deformation that the internal viscosity becomes a secondary factor, leading to only slight changes in
$D$
. Additionally, the evolution of capsule surface area
$\mathcal {A}$
closely mirrors the behaviour of the deformation parameter
$D$
under the influence of both
$\mathcal {R}e$
and
$\lambda$
. At
${\mathcal {C}a} = 1$
and
${\mathcal {R}e} = 20$
,
$\mathcal {A}$
reaches a maximum value approximately 30 % larger than the initial surface area of the capsule.
In terms of the inclination angle,
$\theta$
exhibits a higher value for less deformable capsules (low
$\mathcal {C}a$
). We see from figure 5(d) that increasing the flow inertia from
${\mathcal {R}e}=0.1$
to
${\mathcal {R}e}=10$
increases the inclination angle for all
$\mathcal {C}a$
values at
$\lambda =1$
. Capsules with
${\mathcal {C}a} \geqslant 0.2$
exhibit inclination angles
$\theta \lt \pi /4$
(extensional axis of shear), indicating that more deformable capsules preferentially elongate along the streamwise direction. For low
$\mathcal {C}a$
, the increased flow inertia raises
$\theta$
above
$\pi /4$
, consistent with observations for a single droplet in shear flow (Li & Sarkar Reference Li and Sarkar2005a
; Srivastava et al. Reference Srivastava, Malipeddi and Sarkar2016). A maximal
$\theta \approx 0.9$
is observed at
${\mathcal {C}a}=0.01$
, which significantly affects the interfacial stresses. At
${\mathcal {R}e}=20$
,
$\theta$
is very close to the case at
${\mathcal {R}e}=10$
, indicating that the increase in
$\theta$
due to inertia saturates. Similar to the deformation factor
$D$
, the increase of
$\lambda$
tends to decrease the inclination angle
$\theta$
for all
$\mathcal {C}a$
simulated, as shown in figure 5(e–f). Higher inertia impedes the alignment of the capsule with the flow direction, increasing
$\theta$
. Conversely, increasing the viscosity ratio
$\lambda$
reduces
$\theta$
and enhances the alignment of the capsule with the shear flow. This behaviour highlights the interplay between flow inertia and viscosity ratio in determining capsule orientation, which in turn influences the characteristics of normal stress differences – a topic we will explore further in the following discussion.
In a fully developed shear flow, the capsule exhibits a tank-treading motion, with the angular velocity
$\omega$
decreasing as
$\mathcal {C}a$
increases. Here,
$\omega$
is defined as the time-averaged value across all surface Lagrangian points, after the flow has achieved a fully developed state. Figure 5(g) shows that the increase of
$\mathcal {R}e$
impedes the tank-treading velocity, leading to a smaller magnitude of
$\omega$
(note that the sign of
$\omega$
is determined by the shear flow direction and the chosen coordinate system). This deceleration is primarily due to the capsule being stretched and elongated by inertial forces, which impedes capsule rotation. Similarly, the viscosity ratio
$\lambda$
also contributes to slowing down the angular velocity of the tank treading, as shown in figure 5(h–i). A higher
$\lambda$
means greater internal viscosity, which resists the internal circulation of the fluid within the capsule. An exception is noted at
${\mathcal {R}e}=1$
and
$\lambda =10$
in figure 5(h), where the angular velocity is increased compared with the case at
$\lambda =5$
for
${\mathcal {C}a} \geqslant 0.2$
. This behaviour can be attributed to the very high internal viscosity (
$\lambda =10$
) providing sufficient resistance to prevent excessive stretching, allowing the capsule to maintain a higher rotational speed at
${\mathcal {R}e}=1$
. At
${\mathcal {R}e}=10$
, the acceleration effect observed at
$\lambda =10$
is no longer present in figure 5(i).
4.1.2. Surface principal stretches
Evolution of principal stretches (
$\lambda _1$
,
$\lambda _2$
) on the capsule surface. (a) Sketch of the rotating capsule with the two critical points highlighted,
$\mathcal {C}_1$
(
) denotes a point on the capsule outline in the shear plane
$x{-}y$
across the capsule centroid and
$\mathcal {C}_2$
(
) is an extreme point along
$z$
-axis. (b) A snapshot of the principal stretch distribution in the
$\lambda _1$
$\lambda _2$
space (
) with the two critical points highlighted in corresponding colours. (c–f) Temporal evolution of the principal stretches on the two critical points
$\mathcal {C}_1$
and
$\mathcal {C}_2$
, with red point darkness denoting the time evolution (from light to dark).
To gain deeper insights into the deformation of the membrane and capsule dynamics, we explore the evolution of the principal stretches of the elastic membrane, denoted as
$\lambda _1$
and
$\lambda _2$
, presented in figure 6. We focus on two specific regions on the capsule surface that undergo significant deformation. The first point,
$\mathcal {C}_1$
, is located on the capsule outline in the shear
$x{-}y$
plane initially located at
$(x,y,z)=(-1, 0, 0)$
. The second point,
$\mathcal {C}_2$
, is found at an extreme of the capsule initially at
$(x,y,z)=(0, 0, 1)$
, where the membrane elastic stress reaches its maximum values, as shown in figure 4. The specific locations of the two critical points are illustrated in figure 6(a) for a better clarity.
In figure 6(b), we present a snapshot at
$t=90$
that showcases the phase diagram of principal stretches,
$\lambda _1$
and
$\lambda _2$
, on the whole membrane of a single capsule at
$\lambda =10$
and
${\mathcal {C}a}=1.0$
in the flow at
${\mathcal {R}e}=1$
. The red dots () represent the principal stretches at the point
$\mathcal {C}_1$
, which follow a periodic trajectory loosely resembling an infinity symbol
$\infty$
, as depicted in figure 6(d). This trajectory highlights a significant variation in stretches at the point
$\mathcal {C}_1$
, exhibiting values of
$\Delta \lambda _1 \approx 0.25$
and
$\Delta \lambda _2 \approx 0.2$
. In contrast, the evolution of principal stretches at the point
$\mathcal {C}_2$
is marked by a stationary dot (), indicating constant principal stretches throughout the simulation, as illustrated in figures 6(b) and 6(d). Point
$\mathcal {C}_2$
consistently shows the maximum
$\lambda _2$
value (also the maximum elastic stress
$\sigma _{2,max}$
) during the simulation. The remaining surface points are represented by blue dots () in figure 6(b). These blue dots are positioned between the stationary yellow dot (
$\mathcal {C}_2$
) and the dynamic trajectory of the red dots (
$\mathcal {C}_1$
). The principal stretches at all points on the capsule surface are governed by the dynamic behaviour of the two critical points,
$\mathcal {C}_1$
and
$\mathcal {C}_2$
.
To elucidate the impact of flow inertia and viscosity ratio, we examine additional cases of the trajectories of the principal stretches at critical points
$\mathcal {C}_1$
and
$\mathcal {C}_2$
, as illustrated in figure 6(c–f). The varying brightness of the red dots, from light to dark, illustrates the temporal evolution of the principal stretches at
$\mathcal {C}_1$
. The staggered colour darkness indicates that point
$\mathcal {C}_1$
follows a periodic trajectory, revealing multiple cycles of evolution of the principal stretches. A comparison of figures 6(c) and 6(d) reveals that an increase in
$\mathcal {C}a$
significantly broadens the range of principal stretches of both
$\lambda _1$
and
$\lambda _2$
, indicating more pronounced deformation. Similarly, increasing
$\mathcal {R}e$
leads to a higher maximum
$\lambda _2$
, from
$1.09$
to
$1.2$
, as shown in figures 6(c) and 6(f). Furthermore, figures 6(e) and 6(f) show that an increase in
$\lambda$
results in a slight decrease in
$\lambda _2$
, while maintaining a similar range for
$\lambda _1$
.
In summary, these trajectory patterns are influenced by the interplay of the flow parameters
$\mathcal {C}a$
,
$\mathcal {R}e$
and
$\lambda$
, which collectively dictate the deformation and dynamic behaviour of the capsule. Our analysis yields three important novel insights: (i) the trajectory of
$\mathcal {C}_1$
displays periodic behaviour, characterised by extensive variability and complex patterns at different
$\mathcal {R}e$
,
$\mathcal {C}a$
and
$\lambda$
; (ii) the trajectory of
$\mathcal {C}_2$
remains static within the principal stretch (
$\lambda _1$
,
$\lambda _2$
) space, suggesting that the tank-treading motion of the capsule maintains stability in the examined cases; (iii) the principal stretch evolution at points
$\mathcal {C}_1$
and
$\mathcal {C}_2$
governs the capsule overall deformation, with all other surface points exhibiting principal stretches that lie strictly between those of
$\mathcal {C}_1$
and
$\mathcal {C}_2$
.
4.1.3. Correlation of the deformation factor
Based on the previous discussion, it is clear that the capsule deformation is significantly influenced by three primary factors: flow inertia, elasticity of the capsule and viscosity ratio. A quantitative understanding of their roles is essential for optimising applications in biomedical engineering, where precise control over capsule behaviour is critical. In pursuit of this goal, we have formulated an empirical correlation to predict the capsule deformation, by incorporating
$\mathcal {R}e$
,
$\mathcal {C}a$
and
$\lambda$
\begin{align} \hat {D} &= \underbrace {0.21({\mathcal {C}a})^{0.23}}_{f_1}\underbrace {\left (0.73\log ({\mathcal {C}a}) + 0.027{\mathcal {R}e}\lambda ^{-0.097} + 3.6\right )}_{f_2}\underbrace {\exp \left (-\frac {0.078\lambda }{\exp (0.32{\mathcal {R}e})}\right )}_{f_3} \notag\\&\quad - \underbrace {0.23{\mathcal {C}a}}_{f_4} ,\end{align}
where
$\hat {D}$
denotes the predicted Taylor deformation factor. The first factor,
$f_1$
, represents the base influence of
$\mathcal {C}a$
, while the mixed factor,
$f_2$
, describes the combined nonlinear effects of
$\mathcal {C}a$
,
$\mathcal {R}e$
and
$\lambda$
. Overall, the deformation
$\hat {D}$
increases with higher values of both
$\mathcal {R}e$
and
$\mathcal {C}a$
;
$\lambda$
acts as a damping factor on the deformation, with its effect being most significant at low
$\mathcal {R}e$
and
$\mathcal {C}a$
, where surface elasticity and viscous resistance dominate, limiting deformation. The exponential decay term
$f_3$
shows that the damping effect due to
$\lambda$
decreases as
$\mathcal {R}e$
increases, allowing more significant deformation at higher flow inertia, which captures the transition to the inertial regime. The final term
$f_4$
, which subtracts a linear function of
$\mathcal {C}a$
, provides a stabilising effect that adjusts the overall deformation. This empirical correlation allows us to quantify the complex interplay between fluid dynamics and material properties that govern the capsule deformation.
Correlation prediction of the Taylor deformation factor
$\hat {D}$
in comparison with our computed results. (a–d) Evolution of
$D$
in function of
$\mathcal {C}a$
at
${\mathcal {R}e}=0.1 \sim 20$
; (e) direction comparison between the prediction
$\hat {D}$
and numerical results
$D$
.
Figure 7(a–d) presents a meticulous comparison between our numerical results (obtained by high-fidelity interface-resolved simulations, illustrated with coloured lines) and the results predicted by the correlation (as formulated in (4.1) and indicated by red dots) for the evolution of
$D$
of a single capsule. The comparison reveals the remarkable accuracy of the correlation in predicting the deformation factor
$D$
across the entire range of parameter space (
${\mathcal {R}e}, {\mathcal {C}a}, \lambda$
) explored. The panels distinctly illustrate the increase of
$D$
as a function of
$\mathcal {R}e$
, while the variations of
$D$
with respect to
$\mathcal {C}a$
are also accurately reproduced. Additionally, the influence of
$\lambda$
on the decrease of
$D$
is comprehensively delineated in figure 7(a–d). In general, the correlation in (4.1) demonstrates high precision in predicting the deformation factor
$D$
as depicted in figure 7(e), achieving an average relative error of
$\overline {\varepsilon }_r=7.21\%$
and a determination coefficient of
$r^2=0.9918$
. This level of accuracy underscores the robustness of our empirical correlation in capturing the complex dynamics governing capsule deformation. Moreover, it serves as a valuable tool for experimentalists and engineers, providing a reliable means to estimate capsule deformation and facilitating better design of microfluidic devices.
4.2. Inertial effects on particle stress
Effects of inertia on the particle stress of a single capsule in a shear flow at
${\mathcal {R}e}=0.1\sim 20$
and
$\lambda = 0.1\sim 10$
. (a-c) Particle shear stress
$\Sigma _{xy}^p/\phi$
, (d-f) first normal stress difference
$N_1^p/\phi$
, (g-i) second normal stress difference
$N_2^p/\phi$
.
In addition to capsule deformation, the particle stress exerted on the capsule membrane is also of vital importance. Understanding the behaviour of a single capsule provides valuable insights into the dynamics and rheology of suspensions of capsules. Figure 8 illustrates the evolution of particle shear stress
$\Sigma _{xy}^p$
, first and second normal stress differences (
$N_1^p$
,
$N_2^p$
) as a function of
$\mathcal {C}a$
for a single capsule. The influence of flow fluctuations is not accounted for in this section due to the minimal volume fraction. We divide the particle stress by its volume fraction
$\phi$
in figure 8 for a better visualisation.
From figure 8(a), we observe that an increase in
$\mathcal {R}e$
leads to a clear rise in
$\Sigma _{xy}^p$
. At
${\mathcal {R}e} \leqslant 1$
, higher
$\mathcal {C}a$
values tend to reduce the shear stress, whereas for
${\mathcal {R}e} \geqslant 10$
,
$\mathcal {C}a$
has an opposite effect. This behaviour mainly results from changes in two key factors affecting
$\Sigma _{xy}^p$
: (
$i$
) the capsule elongation, which is proportional to
$r_1$
and (
$ii$
) its inclination angle
$\theta$
. A longer
$r_1$
and a higher
$\theta$
contribute to increase the shear stress
$\Sigma _{xy}^p$
. For instance, the increase in
$\Sigma _{xy}^p$
with
$\mathcal {R}e$
in figure 8(a) is mainly due to the more elongated capsule at higher inertia. As the capsule deformation intensifies with increasing
$\mathcal {C}a$
, there is a corresponding rise in
$r_1$
, while
$\theta$
undergoes a noticeable reduction. In the non-inertial regime, such as at
${\mathcal {R}e}=0.1$
, the increase in
$r_1$
is modest while the effect of decreasing
$\theta$
is more pronounced, leading to a decrease in shear stress with
$\mathcal {C}a$
as shown in figure 8(a). At higher flow inertia (
${\mathcal {R}e}=10$
),
$r_1$
increases more significantly (
$r_1=2.0$
at
${\mathcal {C}a}=1$
), which outweighs the effect of the decreasing
$\theta$
, causing shear stress to rise with
$\mathcal {C}a$
. The interplay of
$r_1$
and
$\theta$
becomes clearer when considering the effects of
$\lambda$
on the particle shear stress, as shown in figures 8(b) and 8(c). We see in figure 5 that a higher
$\lambda$
hinders capsule elongation, while simultaneously reducing the inclination angle
$\theta$
, which results in a clear decreasing trend of
$\Sigma _{xy}^p$
with
$\lambda$
for
${\mathcal {C}a} \geqslant 0.1$
. For capsules at
${\mathcal {C}a}=0.01$
and
${\mathcal {C}a}=0.05$
,
$\Sigma _{xy}^p$
increases with
$\lambda$
due to the stronger internal viscous resistance and corresponding limited deformation.
Figure 8(d–f) illustrates the evolution of the first normal stress difference
$N_1^p = \Sigma _{xx}^p - \Sigma _{yy}^p$
as a function of
$\mathcal {C}a$
at various
$\mathcal {R}e$
and
$\lambda$
. In figure 8(d), we observe that at
${\mathcal {C}a} \leqslant 0.1$
, the flow inertia tends to decrease
$N_1^p$
, while at
${\mathcal {C}a} \geqslant 0.2$
, higher
$\mathcal {R}e$
leads to an increase in
$N_1^p$
in the shear flow. These two opposing behaviours can also be explained by the interplay between capsule elongation and inclination. At low
$\mathcal {C}a$
, the high inclination angle
$\theta$
(as depicted in figure 5
e) promotes normal stress along the
$y$
direction
$\Sigma _{yy}^p$
. For the same elastic force along
$r_1$
, a larger
$\theta$
results in a greater component projected onto the
$y$
-direction. Conversely, at
${\mathcal {C}a} \geqslant 0.1$
, the more significant elongation of the capsule together with the smaller
$\theta$
, enhances normal stress along the streamwise direction (
$x$
-axis)
$\Sigma _{xx}^p$
. Different from rigid particle suspensions (Haddadi & Morris Reference Haddadi and Morris2014), the sign of
$N_1^p$
changes in the capsule suspension. This change of sign aligns closely with the variation of the capsule inclination angle
$\theta$
, similar to what is observed in studies of droplet emulsions (Srivastava et al. Reference Srivastava, Malipeddi and Sarkar2016). We see that the reduction in
$\theta$
from
$\pi /4$
(extension axis of shear flow) in figure 5(d–f) drives the transition of
$N_1^p$
from negative to positive values for elastic capsules. For emulsions of droplets, with a constant surface tension, the stresslet term in (2.9) is purely geometric, which means that the sign change arises directly from the inclination angle exceeding
$\pi /4$
(Singh et al. Reference Singh, Li and Sarkar2014). Please note that while similar mechanism applies, the sign change in
$N_1^p$
is not purely geometric for elastic capsules in figure 8(d–f), due to the distinct nature of their membrane. In figure 8(e–f), the viscosity ratio
$\lambda$
clearly mitigates
$N_1^p$
. Higher internal viscosity helps maintain a more balanced stress distribution between the
$x$
and
$y$
directions, thereby reducing the difference
$N_1^p$
.
Similarly, figures 8(g–i) depicts the evolution of the second normal stress difference (
$N_2^p = \Sigma _{yy}^p - \Sigma _{zz}^p$
). An increase in
$\mathcal {R}e$
forces higher deformation in the
$x$
$y$
plane and hence higher stress along the
$y$
-axis relative to the
$z$
-axis. Conversely, an increase in
$\mathcal {C}a$
leads to smaller inclination angle
$\theta$
, which reduces
$\Sigma _{yy}^p$
and consequently decreases
$N_2^p$
, even to negative values. At
${\mathcal {R}e}=1$
in figure 8(h), the high internal viscosity promotes a more balanced stress distribution in the shear plane and enhances
$\Sigma _{yy}^p$
with respect to
$\Sigma _{zz}^p$
. However,
$\lambda$
has a very limited impact on
$N_2^p$
at
${\mathcal {R}e}=10$
in figure 8(i).
5. Capsule suspension results
In this section, we explore the behaviour of capsule suspensions in the shear flow. We begin by examining how flow inertia affects capsule deformation and dynamics, including microstructure and lateral diffusion of the capsule suspension. Then, we analyse the impact of inertia on the rheological properties of the suspension, by establishing a connection with the previous results of a single capsule.
5.1. Inertial effects on capsule deformation and dynamics
In a suspension, the interaction between capsules is the primary factor influencing its dynamics and rheological behaviour. In a very dilute system, capsules rarely encounter each other and behave similarly to single capsules in the flow. Hence, in this section, we investigate capsule suspensions with volume fractions in the range
$0.1 \leqslant \phi \leqslant 0.4$
in both the non-inertial regime (
${\mathcal {R}e} = 0.1, 1$
) and the inertial regime (
${\mathcal {R}e} = 10, 20$
).
5.1.1. Bridge structures
Snapshots of the capsule suspension at
$\phi =0.4$
projected on the shear plane
$z=0$
with the flow field coloured by the pressure distribution (red indicates regions of high pressure, while blue represents areas of low pressure) and the perturbed flow streamlines; with viscosity ratios (a–d)
$\lambda =1$
, (e–h)
$\lambda =5$
.
We first examine the capsule interactions in dense suspensions with volume fraction of
$\phi =0.4$
at different
$\mathcal {R}e$
and
$\mathcal {C}a$
in figure 9. The outlines of the capsules are projected onto the shear plane at
$z=0$
in black contours. Additionally, figure 9 illustrates the pressure distribution in the same shear plane, with regions of low and high pressure indicated by blue and red colours, respectively. In shear flow, the streamwise velocity along the
$x$
-axis is significantly greater than its transverse and spanwise components. To better illustrate the flow dynamics, we subtract the imposed shear flow velocity from the total velocity and present the streamlines of the perturbed velocity field
$\tilde {u}' = \tilde {u} - \tilde {U}_0 \tilde {y}/\tilde {L}$
in figure 9.
Figures 9(a) and 9(b) show that unsurprisingly capsules become more elongated at higher
$\mathcal {C}a$
. In very dilute systems, the capsules often form pairs. As
$\phi$
increases, these pairs tend to cluster together, forming longer structures. In such regions, an intriguing bridge structure emerges, reminiscent of the rouleaux structures observed in experimental studies of Chinchilla et al. (Reference Chinchilla, Armstrong, Mehri, Savoia, Fenech and Franceschini2021) and Lee & Paeng (Reference Lee and Paeng2021). The formation of the bridge structures is more prevalent at higher values of
$\mathcal {R}e$
and
$\mathcal {C}a$
, as shown in figures 9(c) and 9(d). When the internal fluid viscosity is increased to
$\lambda =5$
, the bridge-like structures in the flow become more prominent, as demonstrated in figures 9(e) to 9(h). It is important to note that these structures do not adhere together as real rouleaux structures would, since no adhesion model is applied to the surface of the capsules. Nonetheless, these observations provide strong evidence that pure hydrodynamic interactions can significantly contribute to the formation of rouleaux structures.
Additionally, all panels in figure 9 demonstrate a marked increase in local fluid pressure during capsule interactions, driven by compressive forces. The primary axis of the bridge structures are typically perpendicular to the extension axis of shear flow, contributing to the anisotropy of the microstructure and enhancing interactions. The anisotropy induced by bridge structures leads to a redistribution of stresses among the different directions. This impacts the particle stress tensor, altering the rheological properties of the suspension, including the first and second normal stress differences. As an additional discussion, the inhomogeneous distribution of membrane elastic stresses on the capsule surface in the bridge structure is presented and analysed in Appendix D.
5.1.2. Average deformation and inclination angle
Effects of flow inertia on the deformation and dynamics of capsule suspensions at
$\lambda =1$
,
${\mathcal {R}e}=0.1\sim 20$
and
$\phi =0.1\sim 0.4$
. (a–c) Ensemble-averaged Taylor deformation factor
$\langle D \rangle$
, with black dashed lines denoting
$\theta =\pi /4$
; (d–f) inclination angle
$\langle \theta \rangle$
; (g–i) angular velocity of tank-treading
$\langle \omega \rangle$
.
We then study the effects of inertia on the deformation of capsules at viscosity ratio
$\lambda =1$
in a suspension with increasing volume fraction. Simulations of suspensions with a higher viscosity ratio
$\lambda =5$
are also performed and analysed. Detailed results and discussions for these cases are provided in Appendix E. In the following, the symbol
$\langle \cdot \rangle$
represents the ensemble average of physical quantities once the flow regime has fully developed, such as
$\langle x \rangle :=\frac {1}{N}\sum _{i=1}^Nx_i$
, with
$N$
the total number of capsules.
In figure 10, we depict the evolution of the ensemble-averaged Taylor deformation factor
$\langle D \rangle$
, inclination angle
$\langle \theta \rangle$
and angular velocity
$\langle \omega \rangle$
as a function of
$\mathcal {C}a$
. In general, capsules exhibit similar deformation evolution with
$\mathcal {C}a$
as a single capsule in shear flow does. At
$\phi =0.4$
,
$\langle D \rangle$
increases with
$\mathcal {C}a$
and
$\mathcal {R}e$
in figure 10(a). At
${\mathcal {R}e} \leqslant 1$
, the effects of flow inertia on variations of
$\langle D \rangle$
are minimal, such that data points nearly overlap for
${\mathcal {R}e}=0.1$
and
${\mathcal {R}e}=1$
at the scale presented. Further increase of
$\mathcal {R}e$
enhances the deformation factor
$\langle D \rangle$
in the inertial regime. In both figures 10(b) and 10(c), the influence of
$\phi$
on
$\langle D \rangle$
is more pronounced for cases at low
$\mathcal {C}a$
. This is attributed to the more intense short-range capsule interactions at low
$\mathcal {C}a$
.
The variation of the inclination angle
$\langle \theta \rangle$
at
$\phi =0.4$
is depicted in figure 10(d). It is clear that flow inertia generally leads to an increase in
$\langle \theta \rangle$
. Conversely, in figure 10(e–f), the decrease of
$\langle \theta \rangle$
with
$\phi$
can be attributed to the reduced free space between capsules along the
$y$
direction. Regarding the average angular velocity
$\langle \omega \rangle$
, an increase in
$\mathcal {R}e$
leads to a lower magnitude of
$\langle \omega \rangle$
at
$\phi =0.4$
, as seen in figure 10(g). In comparison with figure 5(g),
$\langle \omega \rangle$
exhibits similar values as those of the single capsule in the non-inertial regime at
${\mathcal {R}e} \leqslant 1$
. We observe that the impact of
$\phi$
on
$\langle \omega \rangle$
is limited at
${\mathcal {R}e}=1$
in figure 10(h). Interestingly at
${\mathcal {R}e} \geqslant 10$
, the angular velocity of the dense suspension is observed to exceed that of a single capsule. The increase of volume fraction enhance the particle shear stress, which in turn accelerates the tank-treading of the capsules as shown in figure 10(i).
5.1.3. Lateral transport and diffusion
Evolution of the capsule suspension migration in a simple shear flow. (a) Temporal evolution of 50 capsule trajectories along the
$y$
-axis in the case of
$\phi =0.4, \lambda = 1, {\mathcal {R}e}=10, {\mathcal {C}a} = 1.0$
; (b) probability density function (
$f$
, PDF) of the capsule distribution along the
$y$
-axis in the suspension of increasing volume fraction
$\phi$
at
$\lambda =1, {\mathcal {R}e}=10$
and
${\mathcal {C}a}=0.1$
. (c, d) Evolution of the PDF of the standard deviation of the capsule motion along the
$y$
-axis under the effects of
$\phi$
,
$\mathcal {R}e$
,
$\lambda$
and
$\mathcal {C}a$
.
We now examine the diffusion dynamics of the capsules. We focus on how the capsules migrate in the direction of the velocity gradient (along the
$y$
-axis) and in the spanwise direction (along the
$z$
-axis). To provide an intuitive illustration, we present the trajectories of 50 initially randomly distributed capsules in a suspension at
$\phi =0.4$
, as shown in figure 11(a). This figure shows the initial positions of the 50 capsules and their trajectories along the
$y$
-axis over
$t \in [0, 70]$
in the shear flow. It can be observed that the capsules exhibit varying patterns of migration along the
$y$
-axis: some oscillate within a limited region, while others migrate as far as half the computational domain over
$t \in [0, 70]$
. By computing the time average of the number of capsules along the
$y$
-axis after the flow is fully developed, we plot the probability density function (PDF)
$f_y$
of the capsule spatial distribution. Figure 11(b) reveals two depletion layers adjacent to the walls, with small peaks in capsule concentration at
$\tilde {y}/\tilde {L} = \pm 0.4$
. An increase in volume fraction
$\phi$
appears to have only minor effects on the distribution
$f_y$
.
Additionally, we investigate how
$\lambda$
,
$\mathcal {R}e$
and
$\phi$
affect the standard deviation
$\sigma _y$
of capsule trajectory fluctuations along the
$y$
-axis. The corresponding PDFs, denoted
$f_{\sigma _y}$
, are plotted in figure 11(c–d). We observe that an increase in
$\mathcal {R}e$
flattens
$f_{\sigma _y}$
, and raises
$\sigma _{y,peak}$
, the value of
$\sigma _y$
corresponding to the peak of
$f_{\sigma _y}$
, indicating that flow inertia enhances capsule migration. However, increasing
$\lambda$
has the opposite effect, reducing
$\sigma _{y,peak}$
, though its impact on the shape of
$f_{\sigma _y}$
remains minor, as shown in figure 11(c). Figure 11(d) illustrates the evolution of
$f_{\sigma _y}$
with increasing
$\phi$
from
$0.1$
to
$0.4$
. An increase in
$\phi$
leads to more short-range capsule interactions and restricts the free motion of the capsules along the
$y$
-axis. This results in a narrower distribution and a lower
$\sigma _{y,peak}$
.
Subsequently, we take a closer look at capsule migration by computing the MSD using (2.13). The MSD represents the average squared distance that a particle has moved along a direction
$\psi$
over a time interval
$\Delta t$
. In an unsheared Brownian solid sphere system, without a suspending fluid, the MSD is quadratic (
$\propto \Delta t^2$
) for small
$\Delta t$
, known as the ballistic regime, and linear (
$\propto \Delta t$
) for large
$\Delta t$
, known as the diffusive regime. In a simple shear flow, the Brownian contributions are limited, and particle diffusion becomes primarily shear induced and anisotropic (Eckstein et al. Reference Eckstein, Bailey and Shapiro1977; Leighton & Acrivos Reference Leighton and Acrivos1987; Breedveld Reference Breedveld2000). Since the shear rate is constant in a simple shear flow, particle diffusion is solely due to their interactions (Krüger Reference Krüger2012). With the presence of neighbouring particles, straight-line motion is impeded, resulting in non-affine displacement and diffusive motion. In the following, we investigate the diffusion properties of elastic capsules and reveal particular properties of the capsule diffusivity
$D_{\psi }$
with the increase of flow inertia and capsule volume fraction. To wall effects, we will focus mainly on the diffusion behaviour of the capsules along the
$z$
-axis in the suspension.
Diffusivity of the capsule suspension as a function of
$\mathcal {C}a$
and
$\mathcal {R}e$
. (a-b) Comparison between
$D_y$
and
$D_z$
; (c-d) effects of flow inertia
$\mathcal {R}e$
and volume fraction
$\phi$
.
We begin by presenting the overall behaviour of capsule diffusivity
$D_{\psi }$
, computed using (2.14). The evolution of
$D_{\psi }$
along the
$y$
- and
$z$
-axes as a function of
$\mathcal {C}a$
in a suspension at
$\phi =0.3$
is shown in figures 12(a) and 12(b), respectively. The high flow inertia leads to an apparent increase in diffusivity in both directions, while the effects of
$\mathcal {C}a$
remain minor. Depending on the flow regime, the volume fraction
$\phi$
has two opposite effects on the capsule diffusion. In the non-inertial regime, the increase of
$\phi$
promotes the capsule diffusion
$D_z$
in figure 12(c), which simply results from enhanced interactions between capsules. Conversely, at
${\mathcal {R}e} \geqslant 10$
,
$D_z$
decreases with
$\phi$
, indicating that a high
$\phi$
impedes capsule diffusion along the
$z$
-axis. One possible explanation is that, at higher inertia, the streamline movement of capsules tends to dominate, reducing the influence of lateral motion, especially in more densely packed suspensions where inter-particle interactions are stronger and hinder transverse movements. A similar tendency is also observed at
${\mathcal {C}a}=1.0$
in figure 12(d), the most deformable case in this study. There is a clear regime transition in the capsule suspension between the non-inertial regime
${\mathcal {R}e}=1$
and the inertial regime
${\mathcal {R}e}=10$
.
(a) Evolution of MSD of the capsule suspension for increasing
$\phi$
; (b) evolution of capsule diffusivity
$D_z$
within the Reynolds number range
$1 \leqslant {\mathcal {R}e} \leqslant 10$
, showing the transition from the non-inertial to the inertial regime.
Figure 13 presents the MSD along the
$z$
-axis of the capsule suspension for
$1 \leqslant {\mathcal {R}e} \leqslant 10$
and
$0.1 \leqslant \phi \leqslant 0.4$
. In the
$\log -\log$
plots, the MSD evolution shows a slope
$\propto \Delta t^2$
in the ballistic regime and
$\propto \Delta t$
in the diffusion regime, consistent with the behaviour observed in solid particle suspensions. This transition process exists across all the parameter ranges investigated. In particular, in figure 13(a), increasing the capsule volume fraction leads to a decrease of MSD at
${\mathcal {R}e}=10$
, which supports the observations of figure 12(c). The decrease of MSD is more pronounced in the diffusion regime than in the ballistic regime. This trend is opposite to the case at
${\mathcal {R}e} \leqslant 1$
, confirming the existence of a transition from the non-inertial regime to the inertial regime. To better capture this novel transition process, we perform additional simulations and plot the diffusivity
$D_z$
as a function of
$\phi$
for a capsule suspension at
${\mathcal {C}a}=0.1$
in figure 13(b). At other values of
$\mathcal {C}a$
, the evolution of
$D_z$
mirrors that at
${\mathcal {C}a}=0.1$
, and therefore, the results are not plotted to eliminate redundancy. Figure 13(b) clearly shows that
$D_z$
exhibits an increasing trend with
$\phi$
at
${\mathcal {R}e}=1$
, similar to the behaviour observed in solid particle suspensions. As
$\phi$
increases, the probability of strong short-range capsule interactions rises, enhancing mixing and their lateral migration along the
$z$
-axis. Generally, increasing flow inertia leads to higher capsule diffusivity due to greater lateral motion following capsule interactions at all
$\phi$
levels. The enhancement of lateral motion by
$\mathcal {R}e$
along the
$z$
-axis is more evident at low volume fractions (e.g.
$\phi =0.1$
). However, a further increase in
$\phi$
impedes the lateral motion of capsules, as the probability of encountering a neighbouring capsule with opposite lateral motion rises. Thus, the steric interactions between neighbouring capsules result in a decrease in
$D_z$
with
$\phi$
, which is observed for
${\mathcal {R}e} \geqslant 5$
in figure 13(b). At
${\mathcal {R}e}=2.5$
the capsule diffusivity exhibits similar values across the range
$\phi = 0.1$
to
$0.4$
. Based on the current dataset, the critical
${\mathcal {R}e}_c$
for the non-inertial regime to inertial regime transition is identified to be between
${\mathcal {R}e}=1$
and
${\mathcal {R}e}=2.5$
.
Additionally, knowing that in the absence of capsule interactions, a single capsule at steady state exhibits little lateral motion. It is presumed that there exists a critical volume fraction (
$\phi _c \leqslant 0.1$
) at which
$D_z$
reaches its maximum value for the suspension at
${\mathcal {R}e}=10$
, before it decreases as depicted in figure 13(b).
5.2. Inertial effects on suspension rheology
Building on the insights gained from the deformation and diffusion of capsule suspensions, we now investigate the effects of flow inertia on their rheology. We begin by examining the relative viscosity, followed by an analysis of the normal stresses.
5.2.1. Relative viscosity
Relative viscosity
$\mu _r$
of the capsule suspension as a function of
$\phi$
. Einstein’s correlation is shown as a solid black line, Batchelor’s as a black dashed line and the effects of
$\mathcal {C}a$
are illustrated from light blue (low
$\mathcal {C}a$
) to dark blue (high
$\mathcal {C}a$
).
Figure 14 shows the evolution of relative viscosity defined in (2.11) as a function of
$\phi$
at
$0.1 \leqslant {\mathcal {R}e} \leqslant 20$
. The simulation results for capsule suspension, at
$0.01 \leqslant {\mathcal {C}a} \leqslant 1.0$
, are depicted using colours from light to dark blue. In figure 14(a), it is clear that most data points at
${\mathcal {R}e}=0.1$
lie between the correlations of Einstein (Reference Einstein1911) and Batchelor & Green (Reference Batchelor and Green1972). Einstein’s equation is a simple model for predicting the relative viscosity of a dilute suspension of non-interacting, rigid spherical particles. It assumes that the particle concentration is very low, so interactions between particles are negligible. Batchelor’s equation includes a quadratic term to account for hydrodynamic interactions and Brownian motion in a more concentrated suspension. We observe that less deformable capsules (at low
$\mathcal {C}a$
) exhibit relative viscosities close to Batchelor’s correlation. As the deformability of the capsule increases (at a higher
$\mathcal {C}a$
) the value of
$\mu _r$
tends to decrease for all volume fractions investigated (
$0.1 \leqslant \phi \leqslant 0.4$
). One possible reason is that more deformable capsules experience reduced hydrodynamic interactions because they can adapt to the flow and align more smoothly with the fluid streamlines, thus reducing
$\mu _r$
. As
$\mathcal {R}e$
increases, we have higher
$\mu _r$
especially at small
$\mathcal {C}a$
. In figure 14(b), the relative viscosity
$\mu _r$
of the suspension at
${\mathcal {C}a} \leqslant 0.2$
largely exceeds the values predicted by Batchelor’s correlation. In the inertial regime, at
${\mathcal {R}e} \geqslant 10$
, the values of
$\mu _r$
further increase, and both correlations fail to capture the numerical data of
$\phi \geqslant 0.2$
in figure 14(c–d).
Considering flow inertia is necessary for more accurate predictions and better control over the rheological properties of capsule suspensions within our parameter space. In the following, we establish several new correlations to better predict the relative viscosity
$\mu _r$
of the capsule suspension in the inertial regime. We first extend Batchelor’s correlation and propose the following empirical correlation based on the entire numerical dataset within the parameter space
$0.1 \leqslant \phi \leqslant 0.4$
,
$0.1\leqslant {\mathcal {R}e} \leqslant 20$
and
$0.01 \leqslant {\mathcal {C}a} \leqslant 1.0$
using a multivariate regression
\begin{equation} \mu _r = 1.34 - 6.68\phi + \left(-0.23 + 11.13\phi ^{1.19}\right)\mathcal {C}a^{-0.066}\left(1 + 0.081\mathcal {R}e^{0.62}\right). \end{equation}
Instead of a simple quadratic term, our multivariate correlation explicitly incorporates the two additional variables
$\mathcal {C}a$
and
$\mathcal {R}e$
to account for capsule deformability and flow inertia, respectively. The correlation shows good agreement with our computed results, accurately capturing the increasing trend of
$\mu _r$
with
$\phi$
and the decreasing trend with
$\mathcal {C}a$
. Although there exist a few discrepancies at low
$\mathcal {C}a$
and low
$\mathcal {R}e$
, the model demonstrates an average relative error
$\langle \epsilon \rangle =4.9\,\%$
and a determination coefficient of
$r^2 = 0.9916$
.
A second way to predict the relative viscosity of the capsule suspension is to incorporate the deformability of the capsules into the volume fraction
$\phi$
and apply an Eilers-type model (Aouane et al. Reference Aouane, Scagliarini and Harting2021)
\begin{equation} \mu _r = \left (1 + \frac {\alpha \phi _e}{1 - \beta \phi _e} \right )^2, \end{equation}
where
$\alpha$
and
$\beta$
are coefficients of the model, and
$\phi _e$
denotes the effective volume fraction of the capsule suspension. In the non-inertial regime, Rosti et al. (Reference Rosti, Brandt and Mitra2018) proposed the following effective volume fraction, considering an effective sphere:
\begin{equation} \phi _{e, Rosti} = \frac {4}{3}\pi \langle r_3 \rangle ^3\frac {N}{V}. \end{equation}
However, this definition is not well suited for elongated capsules in inertial flow, especially for
${\mathcal {R}e} \geqslant 10$
. Therefore, we propose a novel alternative definition of the effective volume fraction, considering an effective spheroid
\begin{equation} \phi _e = \frac {4}{3}\pi \langle r_1 \rangle \langle r_3 \rangle ^2\frac {N}{V}. \end{equation}
This definition emphasises the role of particle elongation under flow inertia and proves to be more effective for the parameter space investigated, particularly for capsules in inertial flow. The following analysis is based on the novel spheroidal effective volume fraction defined by (5.4).
(a–c) Relative viscosity
$\mu _r$
of the capsule suspension as a function of
$\phi _e$
. The panels illustrate cases with increasing flow inertia
$\mathcal {R}e$
; results at different
$\phi$
are represented using distinct markers. (d) Prediction of
$\mu _r$
using Eilers-type model.
In figure 15, we present the relative viscosity
$\mu _r$
of the capsule suspension as a function of the effective volume fraction
$\phi _e$
for
$0.1 \leqslant {\mathcal {R}e} \leqslant 20$
. To enhance clarity, we show the cases separately at different
$\mathcal {R}e$
in the panels. For a single capsule in the shear flow, we have
$\mu _r = 1$
as
$\phi _e \to 0$
. In the non-inertial regime, at
${\mathcal {R}e} \leqslant 1$
, the values of
$\mu _r$
collapses to a master curve as a function of
$\phi _e$
, as shown in figure 15(a–b). In the inertial regime, the influence of
$\mathcal {C}a$
on
$\mu _r$
becomes more pronounced, as illustrated in figure 15(c). Although the numerical results still approximate a master curve, the data points are grouped by
$\phi$
and show more outliers at small
$\mathcal {C}a$
. Consistent with the lateral migration of capsules discussed above, the behaviour of the capsule suspension in the inertial regime differs from that in the non-inertial regime. To accurately model the relative viscosity of the capsule suspension, it is essential to consider the effects of inertia and the transition between these two regimes.
Given that the values of
$\mu _r$
as a function of
$\phi _e$
fall close to a master curve, we can establish new correlations for
$\mu _r$
using the Eilers-type models in (5.2). For each
$\mathcal {R}e$
investigated, we determine the coefficients
$\alpha$
and
$\beta$
using the least-square method, as listed in table 2. The prediction results are presented and compared with our computed results in figure 15(d). The first observation is that all the numerical data for
${\mathcal {R}e}=0.1$
and
${\mathcal {R}e}=1$
collapse onto the same master curve (solid black line in figure 15
d). Consequently, we use the same coefficient set (
$\alpha$
,
$\beta$
) for modelling
$\mu _r$
for these cases in the non-inertial regime, as shown in table 2. The obtained correlation exhibits high accuracy, with an average relative error
$\langle \epsilon \rangle =2.3\,\%$
for
${\mathcal {R}e}=0.1$
and
${\mathcal {R}e}=1$
. In the inertial regime, the effects of
$\mathcal {R}e$
on the viscosity evolution become more pronounced, necessitating separate modelling for each
$\mathcal {R}e$
. The best fit yields a relative error of
$\langle \epsilon \rangle =4.1\,\%$
for
${\mathcal {R}e}=10$
and
$\langle \epsilon \rangle =7.8\,\%$
for
${\mathcal {R}e}=20$
, using the parameters listed in table 2.
Correlation coefficients
$\alpha$
,
$\beta$
, relative error
$\langle \epsilon \rangle$
and determination coefficient for Eilers-type models in (5.2).
Prediction of models based on the effective volume fraction
$\phi _e$
using the inertia-tuned model.
Even though Eilers-type models show a promising ability to fit relative viscosity data for suspensions in the non-inertial regime, the dependence on
$\mathcal {R}e$
of the model coefficients also makes it less convenient to use and extrapolate at
${\mathcal {R}e} \geqslant 10$
. To address this, we extend the Eilers-type model by introducing an inertia factor of the form
$A+B{\mathcal {R}e}^C$
, resulting in the inertia-tuned model
\begin{equation} \mu _r = 1 + \left (\frac {2.8\phi _e}{1 + 2.59\phi _e}\right )^2 \left (7.92 + 1.65\mathcal {R}e^{0.63}\right ). \end{equation}
Now, with the inertia-tuned model in (5.5), the entire data set is well captured with a single set of coefficients, as shown in figure 16(a). The inclusion of the inertia factor allows for accurate prediction of the dependence of
$\mu _r$
on
$\mathcal {R}e$
. As shown in figure 16(b), the inertia-tuned model in (5.5) achieves a satisfactory determination coefficient of
$r^2 = 0.9811$
.
Overall, the inertia-tuned model is more physically interpretable and has fewer empirical coefficients than the multivariate correlation in (5.1). However, it is important to note that (5.5) requires the a priori computation of the effective volume fraction
$\phi _e$
defined in (5.4). In contrast, the multivariate correlation in (5.1) depends directly on the physical properties of the system:
$\mathcal {R}e$
,
$\mathcal {C}a$
and
$\phi$
, and does not require the computation of
$\phi _e$
. In general, the correlations and models derived in this section can be of specific interest in different applications, depending on the particular requirements.
5.2.2. Normal stress of the suspension
Effects of inertia on the particle stress of capsule suspensions at
$\lambda =1$
,
${\mathcal {R}e}=0.1\sim 20$
and
$\phi =0.1\sim 0.4$
. Panels (a–c) depict the first normal stress difference
$N_1^p$
, panels (d–f) illustrate the second normal stress difference
$N_2^p$
and panels (g–i) present the bulk stress due to flow fluctuation
$N_1^f$
and
$N_2^f$
.
In addition to relative viscosity, we present the evolution of first and second normal stress differences of capsule suspension at
$\lambda =1$
in figure 17. We begin by analysing the particle stress contributions arising from capsule interactions,
$N_1^p$
and
$N_2^p$
. Subsequently, we take a closer examination of the stress resulting from flow fluctuations,
$N_1^f$
and
$N_2^f$
.
For the dense suspension at
$\phi =0.4$
, the first normal stress difference,
$N_1^p$
, predominantly shows positive values except in cases where
${\mathcal {R}e} \geqslant 10$
and
${\mathcal {C}a}=0.01$
, as depicted in figure 17(a). Both
$\mathcal {R}e$
and
$\mathcal {C}a$
contribute to the increase of
$N_1^p$
, with the effect of
$\mathcal {C}a$
being more pronounced in the inertial regime. This trend closely resembles the behaviour of a single capsule in shear flow plotted in figure 8(d). The high similarity suggests that the evolution of
$N_1^p$
with
$\mathcal {C}a$
in the suspension results from the collective effects of flow inertia on individual capsules. Figures 17(b) and 17(c) illustrate that a higher
$\phi$
significantly increases
$N_1^p$
in both the non-inertial and inertial regimes. Compared with figure 10(d–f), the observed sign change of
$N_1^p$
aligns closely with the variation in the capsule inclination angle
$\langle \theta \rangle$
, similar to the single capsule case discussed before. The reduction in
$\langle \theta \rangle$
from
$\pi /4$
(extension axis of shear flow) drives the transition of
$N_1^p$
from negative to positive values, particularly pronounced at
${\mathcal {R}e} = 10$
. This phenomenon mirrors a similar mechanism identified in both dilute and dense droplet emulsions, as reported by Li & Sarkar (Reference Li and Sarkar2005
a) and Srivastava et al. (Reference Srivastava, Malipeddi and Sarkar2016). Another possible contribution is the inhomogeneous microstructure of the suspension. Bridge structures often indicate regions of strong interaction and deformation of capsules, which can contribute to the considerable enhancement of particle stress,
$\Sigma _{xx}^p$
, along the streamwise direction in inertial flow, thus increasing
$N_1^p$
.
Regarding
$N_2^p$
, we observe a decreasing trend as a function of
$\mathcal {C}a$
at
$\phi =0.4$
in figure 17(d). At
${\mathcal {R}e} \leqslant 1$
,
$N_2^p$
exhibits negative values at high
$\phi$
, different from the positive values observed in the single capsule case in figure 8(g). The enhanced flow inertia at
${\mathcal {R}e} \geqslant 10$
results in positive
$N_2^p$
at small
$\mathcal {C}a$
. In figure 17(e), we note that a high volume fraction
$\phi$
leads to a significant decrease in
$N_2^p$
at
${\mathcal {R}e}=1$
. The negative values of
$N_2^p$
indicate that the increase in normal stress along the
$z$
-axis is more significant than along the
$y$
-axis as
$\phi$
increases. From the point of view of the suspension microstructure, the capsules forming the bridge structures provide ample free space for other capsules in the shear plane, potentially reducing the overall normal stress
$\Sigma _{yy}^p$
in the gradient direction. In contrast, along the
$z$
-axis, the normal stress
$\Sigma _{zz}^p$
increases with
$\phi$
due to the rise in binary interactions. At
${\mathcal {R}e}=10$
in figure 17(f), the sign change of
$N_2^p$
occurs at
$0.2\leqslant {\mathcal {C}a} \leqslant 0.4$
, distinct from that of
$N_1^p$
, as it does not directly correlate with the inclination angle
$\langle \theta \rangle$
.
The evolution of bulk stress due to flow fluctuations is presented in figures 17(g)–17(i). The shear stress
$\Sigma _{xy}^f$
is significantly smaller compared with
$\Sigma _{xy}^p$
and is therefore not plotted here. Flow fluctuations predominantly impact the first normal stress difference,
$N_1^f$
, whose magnitude increases with
$\phi$
, as shown in figure 17(g). Unlike
$N_1^p$
,
$N_1^f$
exhibits negative values for
$\phi$
in the range of
$0.1$
to
$0.3$
. Intuitively, increasing
$\mathcal {R}e$
amplifies the magnitude of
$N_1^f$
. At
$\phi =0.3$
, figure 17(h) illustrates the evolution of
$N_1^f$
as a function of
$\mathcal {C}a$
. Higher capsule deformability results in a sharp rise in the magnitude of
$N_1^f$
, which can surpass
$N_1^p$
, as shown in figures 17(b) and 17(c). Conversely, increasing
$\mathcal {C}a$
reduces the magnitude of the second normal stress difference,
$N_2^f$
, as shown in figure 17(i). This behaviour arises from the reorientation of capsules along the
$x$
-axis in more deformable cases, which mitigates stress along the
$y$
axis.
6. Conclusions
We conducted a comprehensive investigation into the deformation, dynamics and rheology of a single and suspensions of immersed elastic capsules in an inertial shear flow. We performed high-fidelity interface-resolved simulations using octree and multigrid methods to accurately capture the behaviour of elastic capsule suspensions and create an extensive dataset for numerical analysis and modelling.
For a single capsule in the shear flow, we elucidate the promoting effects of flow inertia on capsule deformation and the internal resistance to deformation brought by the viscosity ratio, highlighting the interplay of these two competing factors. We reveal the mechanism behind the stretching of the capsule surface during tank-treading motion. Specifically, we demonstrate that certain critical surface points exhibit constant principal stretches, while others undergo periodic changes. Our study shows that the first and second normal stress differences exhibit a change of sign with increasing flow inertia, indicating a complex interaction between the capsules and the surrounding fluid.
In examining capsule suspensions, the effects of increasing volume fraction on the capsule average deformation, diffusion and rheology are thoroughly discussed. We find the formation of bridge structures due to the pure hydrodynamic interactions among elastic capsules in the shear plane. These interactions enhance the inhomogeneity of the microstructure within the capsule suspension, altering the capsule surface stress distribution and providing new insights into the formation of rouleaux structures in inertial flows. In terms of capsule lateral migration, we identify a critical range of
$\mathcal {R}e$
marking the transition of capsule diffusion from the non-inertial regime to the inertial regime, providing a deeper understanding of the underlying mechanisms at play. Furthermore, we reveal close connections between the behaviour of individual capsules and dense suspensions and demonstrate that inertia significantly affects the deformation and rheology of capsule suspensions. We propose multiple new empirical correlations to predict the capsule deformation factor of a single capsule and the relative viscosity of the capsule suspension, which can serve as valuable tools for future research and practical applications.
Our findings enhance the understanding of elastic capsule behaviour in inertial shear flows, paving the way for advances in inertial microfluidic applications. The proposed empirical correlations offer predictive tools to guide the selection of flow parameters, channel geometries and capsule properties, helping to avoid costly and time-intensive trial-and-error processes. These insights enable the customisation of microchannel geometries to meet specific needs, such as isolating tumour cells or bacteria from whole blood, and optimising device designs to reduce clogging. By understanding capsule deformation and orientation under various shear stresses, device parameters can be precisely tuned to improve separation efficiency and sorting accuracy. The novel insights obtained in this study lay a foundation for future advancements, fostering the development of next-generation inertial microfluidic systems for specific, high-precision applications in healthcare and biotechnology. Future research could extend this work to explore polydisperse capsule suspensions and the behaviour of capsules in channels with complex geometries.
Acknowledgement
G.G. expresses his gratitude to the Pacific Institute of Mathematical Sciences for their support via his PIMS-CNRS postdoctoral fellowship. He also thanks Dr Z. Ge at University of British Columbia, Professor N. Takeishi at Kyushu University and Professor T. Krüger at University of Edinburgh for the insightful conversations.
Funding
The authors greatly appreciate the financial support of the Natural Sciences and Engineering Research Council of Canada (NSERC) via Anthony Wachs New Frontiers in Research Fund grant NFRFE-2018–01922. This research was enabled by support provided by Digital Research Alliance of Canada (https://alliancecan.ca/en) through Anthony Wachs’s 2023–2024 Computing Resources for Research Groups allocation qpf-764-ac. This research was also supported in part through computational resources and services provided by Advanced Research Computing at the University of British Columbia (https://arc.ubc.ca).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Immersed boundary method
In this study, the membrane is discretised using an unstructured triangulation. Equation (2.5) is resolved via a linear finite element method. Communication between the membrane triangulation and the grid utilises the IBM (Peskin Reference Peskin1977, Reference Peskin2002), where the regularisation of the Dirac distribution in (2.3) employs a smoothed 2-point formulation (Yang et al. Reference Yang, Zhang, Li and He2009)
\begin{equation} \tilde {\delta }_h(\boldsymbol {x} - \boldsymbol {x_0}) = \begin{cases} \begin{aligned} & \frac {1}{\tilde {\unicode{x1D6E5} }^3} \prod _{i = 1}^3 \left ( \frac {3}{4} - |x_i - x_{0,i}|^2 \right ) \, \text {if} \quad |x_i - x_{0,i}| \leqslant 0.5 \\ & \frac {1}{\tilde {\unicode{x1D6E5} }^3} \prod _{i = 1}^3 \left ( \frac {8}{9} - \frac {3}{2}|x_i - x_{0,i}| + \frac {1}{2}|x_i - x_{0,i}|^2 \right ) \, \text {if} \quad 0.5 \lt |x_i - x_{0,i}| \leqslant 1.5 \\ & 0 \quad \text {otherwise} \end{aligned} \end{cases}, \end{equation}
where
$\boldsymbol {x_0} = [x_{0,1}, \; x_{0,2}, \; x_{0,3}]$
represents the location of a Lagrangian node on the membrane surface,
$x = \tilde {x}/\tilde {\unicode{x1D6E5} }$
and
$\tilde {\unicode{x1D6E5} }$
is the local mesh size of the Eulerian grid. The prefactor
$1/\tilde {\unicode{x1D6E5} }^3$
ensures that the discrete integral over the entire space,
$\int _\Omega \tilde {\delta }_h(\boldsymbol {x_0} - \boldsymbol {x}) \, {\rm d}\boldsymbol {x}$
, remains equal to one.
The velocity
$\boldsymbol {\tilde {V}_0}$
of a Lagrangian node located at
$\boldsymbol {\tilde {x}_0}$
is interpolated from the Eulerian velocity field
$\boldsymbol {\tilde {u}}$
using
\begin{equation} \boldsymbol {\tilde {V}_0} = \sum _{i\in {supp}(\tilde {\delta }_h(\boldsymbol {x_0} - \boldsymbol {x_i}))} \boldsymbol {\tilde {u}_i}\tilde {\delta }_h(\boldsymbol {x_0} - \boldsymbol {x_i})\tilde {\unicode{x1D6E5} }^3, \end{equation}
where
$\text {supp}(\tilde {\delta }_h)$
denotes the support of the regularised Dirac delta function. Similarly, the membrane force
$\boldsymbol {\tilde {F}_0}$
at a Lagrangian node
$\boldsymbol {\tilde {x}_0}$
is distributed to the Eulerian force density function
$\boldsymbol {\tilde {f}}$
using (2.3). Once the Lagrangian velocities of all the capsule nodes are interpolated from the Eulerian velocity field using (A2), the position of each node is updated using an explicit forward Euler scheme (Krüger Reference Krüger2012; Aouane et al. Reference Aouane, Scagliarini and Harting2021).
Additionally, when the fluid viscosity inside and outside the capsule membrane differ, we adopt the approach developed by Tryggvason and co-workers (Unverdi & Tryggvason Reference Unverdi and Tryggvason1992; Tryggvason et al. Reference Tryggvason, Bunner, Esmaeeli, Juric, Al-Rawahi, Tauber, Han, Nas and Jan2001). A discrete indicator
$\tilde {\chi }$
is computed from a discrete grid-gradient field
$\boldsymbol {\tilde {G}}$
\begin{equation} \boldsymbol {\tilde {G}(x)} = \sum _{i\in \mathcal {T}}\tilde {\mathcal {A}}_{i}\tilde {\delta }_h(\boldsymbol {x_0} - \boldsymbol {x_i})\boldsymbol {n_i}, \end{equation}
where
$\mathcal {T}$
denotes the triangular elements on the capsule surface,
$\mathcal {\tilde {A}}_{i}$
is the surface area of triangle
$i$
,
$\boldsymbol {x_i}$
represents the coordinates of its centroid and
$\boldsymbol {n_i}$
is its inward normal vector. By definition,
$\boldsymbol {\tilde {G}}$
is non-zero only on the supports of
$\tilde {\delta }_h(\boldsymbol {x_0} - \boldsymbol {x_i})$
. The discrete indicator function
$\tilde {\chi }$
can be computed by solving a Poisson problem
\begin{equation} \nabla ^2 \tilde {\chi } = \nabla \cdot \boldsymbol {\tilde {G}}. \end{equation}
Since
$\tilde {\chi }$
is a regularised step function, it remains constant away from the membranes.
Our current implementation of the FTM requires all cells belonging to the supports of
$\tilde {\delta }_h(\boldsymbol {x_0} - \boldsymbol {x_i})$
to be the same size
$\tilde {\unicode{x1D6E5}}$
. In other words, all Eulerian cells surrounding a membrane must be of uniform size corresponding to the smallest grid size
$\tilde {\unicode{x1D6E5} }$
, through adaptive mesh refinement. Note that the volume of the capsule is forced to be conserved in all cases investigated.
Appendix B. Numerical validation of the solver
B.1. Validation of a single capsule in shear flow
Validation of the deformation and particle stress of a single capsule in a shear flow at
$\dot {\gamma } = 1$
,
$\lambda =1$
and
${\mathcal {R}e}=0.01$
as a function of
$\mathcal {C}a$
; (a) Taylor deformation factor, (b) inclination angle, (c-e) three semi-axis lengths of the moment of inertia equivalent ellipsoid
$r_1$
,
$r_2$
,
$r_3$
, (f) particle shear stress
$\Sigma _{xy}^p/\phi$
and (g-h) particle first/second normal stress difference
$N_1^p/\phi$
,
$N_2^p/\phi$
.
Figure 18 illustrates the validation of our solver in comparison with the literature for the simulation of a single spherical elastic capsule subjected to a shear flow at
${\mathcal {R}e} = 0.01$
. The viscosity ratio of the capsule studied is unity. We investigate the evolution of the Taylor deformation factor,
$D$
, as defined in (2.7) and illustrated in figure 18(a). Our simulations align closely with the data reported by (Lac et al. Reference Lac, Barthes-Biesel, Pelekasis and Tsamopoulos2004; Aouane et al. Reference Aouane, Scagliarini and Harting2021; Guglietta et al. Reference Guglietta, Pelusi, Sega, Aouane and Harting2023). To enhance understanding, we also compute the Taylor deformation factor using an alternate, yet popular, definition
\begin{equation} D_{m} = \frac {r_{max} - r_{min}}{r_{max} + r_{min}}, \end{equation}
where
$r_{max}$
and
$r_{min}$
represent the maximum and minimum radial distances from the capsule’s membrane to its centroid within the shear plane. Remarkably, as shown in figure 18(a), applying these two definitions at
${\mathcal {R}e} = 0.01$
results in overlapping curves, demonstrating their consistency. Similar observations extend to the inclination angle of the capsule within the shear plane, as demonstrated in figure 18(b). Despite minor discrepancies in
$\theta _m$
determined by the directions of
$\boldsymbol {r}_{max}$
and
$\boldsymbol {r}_{min}$
in comparison with the definition by (2.8), the fundamental physical principles continue to hold. In the remainder of the current study, we apply definitions of
$D$
and
$\theta$
using (2.7) and (2.8). In figure 18(c–e), we plot the evolution of the lengths of three semi-axis (
$r_1$
,
$r_2$
,
$r_3$
) of the inertia-equivalent ellipsoid in a function of
$\mathcal {C}a$
. Our numerical results show excellent agreement with existing data reported in the literature (Li & Sarkar Reference Li and Sarkar2008; Aouane et al. Reference Aouane, Scagliarini and Harting2021; Guglietta et al. Reference Guglietta, Pelusi, Sega, Aouane and Harting2023).
In terms of validation of forces acting on the capsule membrane, we analyse the particle shear stress
$\Sigma _{xy}^p$
, and the first/second normal stress differences (
$N_1^p$
,
$N_2^p$
), illustrated in figure 18(f–h). Our results indicate that the particle shear stress aligns closely with existing studies (Bagchi & Kalluri Reference Bagchi and Kalluri2010; Aouane et al. Reference Aouane, Scagliarini and Harting2021) as shown in figure 18(f). A mesh convergence analysis on the first and second normal stress differences
$N_1^p$
and
$N_2^p$
as functions of
$\mathcal {C}a$
is presented in figure 18(g–h). We denote
$n_E$
and
$n_L$
as the refinement levels of the Eulerian (fluid phase) and Lagrangian (capsule membrane) meshes, respectively. We note that simulations with a coarse mesh (
$n_E=7$
,
$n_L=3$
) reveal good agreement with the literature. Yet, increasing the mesh resolution modifies the outcomes; for instance, at
$n_E=8$
, we observe discrepancies for
$N_1^p$
and
$N_2^p$
in comparison with the reported findings by (Bagchi & Kalluri Reference Bagchi and Kalluri2010; Aouane et al. Reference Aouane, Scagliarini and Harting2021), underscoring their sensitivity to the mesh size. Further elevation of the Eulerian mesh resolution to
$n_E=9$
demonstrates incremental variance in
$N_1^p$
and
$N_2^p$
. Given this, we opt for
$n_E=8$
and
$n_L=4$
for subsequent simulations of a single capsule in shear flow, a compromise ensuring both high precision and computational efficiency.
B.2. Validation of capsule suspensions in a shear flow
Validation of the deformation and rheology of capsule suspensions in a shear flow at
$\dot {\gamma } = 1$
,
$\lambda =1$
,
${\mathcal {R}e}=0.01$
and
$\phi =0.1$
as a function of
$\mathcal {C}a$
; (a-c) three semi-axis lengths of the moment of inertia equivalent ellipsoid
$\langle r_1 \rangle$
,
$\langle r_2 \rangle$
,
$\langle r_3 \rangle$
; (d) relative viscosity of the suspension
$\mu _r$
; (e-f) particle first/second normal stress difference
$N_1^p$
and
$N_2^p$
.
Now that the solver is validated for single capsule scenarios, we proceed with simulations to examine the interactions between capsules in a suspension. Figure 19 illustrates the validation of our FTM for a suspension in a shear flow at
$\phi =0.1$
and
${\mathcal {R}e}=0.01$
in comparison with data reported in the literature (Aouane et al. Reference Aouane, Scagliarini and Harting2021; Guglietta et al. Reference Guglietta, Pelusi, Sega, Aouane and Harting2023). Our numerical results show a very good agreement with the data reported in Aouane et al. (Reference Aouane, Scagliarini and Harting2021) for the three semi-axis lengths of the moment of inertia equivalent ellipsoid,
$\langle r_1 \rangle$
,
$\langle r_2 \rangle$
and
$\langle r_3 \rangle$
, as seen in figure 19(a–c). The symbol
$\langle \cdot \rangle$
represents the ensemble average of physical quantities after the flow regime has fully developed. The relative viscosity
$\langle \mu _r \rangle$
, an important aspect of the suspension rheology, is also well captured, as shown in figure 19(d). Similar to the single capsule case shown in figure 18(g–h), our simulation results for
$N_1^p$
and
$N_2^p$
show slight differences with those reported in the literature. However, given the higher spatial resolution of our simulations, we are confident that our results are accurate and thoroughly validated.
Appendix C. Repulsive model and its impacts
To prevent capsules from overlapping in a relatively dense suspension, we introduce a repulsive force to simulate the normal component of the lubrication force. This force is activated when the Lagrangian nodes of two neighbouring capsules come into close proximity, specifically when the numerical stencils of the Lagrangian points begin to overlap
\begin{equation} \displaystyle { \boldsymbol {\tilde {f}}_{r,ij} = \left \{ \begin{array}{ll} \mathcal {K} \left [ \left (\frac {2\tilde {\unicode{x1D6E5} } }{\tilde {d}_{ij}} \right )^2 -1 \right ]\hat {\boldsymbol {d}}_{ij} & \text {if} \,\, \tilde {d}_{ij} \lt 2\tilde {\unicode{x1D6E5} }, \\ 0 & \text {if} \,\, \tilde {d}_{ij} \geqslant 2\tilde {\unicode{x1D6E5} }. \\ \end{array} \right . } \end{equation}
The repulsive force acts on nodes
$i$
and
$j$
, which belong to two neighbouring capsules. Let
$\tilde {d}_{ij}$
be the distance between these two nodes and
$\hat {\boldsymbol {d}}_{ij}$
be the unit vector connecting them. The repulsive coefficient
$\mathcal {K} = ||\boldsymbol {\tilde {u}}||/8$
is chosen to be proportional to the local fluid velocity. By implementing this repulsive interaction, we ensure that the capsules maintain a minimum separation distance, preventing any undesired physical overlap.
Effects of the repulsive force on the deformation and relative viscosity of the capsule suspension. Panels (a) and (b) illustrate the effects of the repulsive force on a dilute suspension at
$\phi =0.1$
in the non-inertial regime
${\mathcal {R}e}=0.01$
, while panels (c) and (d) show the results for a dense suspension with volume fraction
$\phi =0.4$
in the inertial regime at
${\mathcal {R}e}=10$
.
In figure 20, we show that the macroscopic behaviour of the suspension remains unaffected by the additional repulsive nodal force. We present the temporal evolution of the Taylor deformation factor
$\langle D \rangle$
and the relative viscosity
$\mu _r$
in a dilute suspension of capsules at
$\phi =0.1$
and
${\mathcal {R}e}=0.01$
. As illustrated in figure 20(a–b), the repulsive force model has no impact on the temporal evolution of both
$\langle D \rangle$
and
$\mu _r$
. In this dilute system, the likelihood of interactions between the capsules is low. The ample free space allows the lubrication effects to be well resolved when two capsules encounter each other, effectively preventing them from approaching too closely. As a result, the repulsive model is rarely activated.
In a dense suspension with
$\phi =0.4$
in an inertial flow at
${\mathcal {R}e}=10$
, interactions among the capsules become significant. As shown in figure 20(c), the repulsive force has minor effects on the temporal evolution of the Taylor deformation factor
$\langle D \rangle$
, while maintaining the same time-averaged value. Similarly, figure 20(d) demonstrates that the repulsive force has minor effects only on the temporal evolution of the relative viscosity
$\mu _r$
, with a relative error of less than
$6\,\%$
compared with the case without the repulsive force.
To conclude, the repulsive force effectively prevents membrane overlapping and accounts for lubrication effects, thereby enhancing numerical stability. Importantly, this is achieved without altering the dynamics or rheology of the capsule suspension in both the non-inertial and inertial regimes.
Appendix D. Capsule elastic stress in the bridge structure
Distribution of the maximum elastic stress
$\sigma _{2,{max}}$
in a capsule suspension at
$\phi = 0.3$
,
${\mathcal {C}a} = 0.05$
and
${\mathcal {R}e} = 0.1$
. (a) Location of the cut plane in the cubic computation domain, (b) snapshot of capsule contours coloured by
$\sigma _{2,{max}}$
.
We plot the capsule elastic stress distribution on the contours of capsules on a cut plane in a suspension with
$\phi = 0.3$
in figure 21. The cut plane, perpendicular to the
$z$
-axis, is located at
$z = 0$
as shown in figure 21(a). The capsule contours in this shear plane are coloured based on
$\sigma _{2,{max}}$
, with red indicating high stresses and blue indicating low values. From figure 21(b), we observe the formation of bridge structures in the suspension. The capsule membrane in the lubrication layer exhibits low
$\sigma _{2,{max}}$
values. The maximum elastic stress,
$\sigma _{2,{max}}$
, is observed at the two extremities of the capsule in the shear plane. The bridge structures are temporary because there is no adhesive force holding the capsules together. Without adhesion, the structures are easily broken apart by the shear flow. Despite their temporary nature, the conditions in the shear flow repeatedly bring capsules into close proximity, leading to the frequent formation of bridge structures, due to the hydrodynamic forces in the suspension.
Appendix E. Effects of viscosity ratio on capsule suspension
Effects of viscosity ratio on the capsule suspensions in a shear flow at
${\mathcal {R}e}=10$
and
$\phi =0.4$
. (a) Ensemble-averaged Taylor deformation factor
$\langle D \rangle$
, (b) inclination angle
$\langle \theta \rangle$
, (c) tank-treading velocity
$\langle \omega \rangle$
, (d) relative viscosity
$\mu _r$
, (e) first normal stress difference
$N_1^p$
and (e) second normal stress difference
$\Sigma _{xy}^p$
.
We briefly present and discuss the effects of the viscosity ratio
$\lambda$
on capsule suspensions at
$\phi = 0.4$
and
${\mathcal {R}e} = 10$
in figure 22. Generally, the effects of
$\lambda$
on a dense suspension are minor compared with the effects of
$\phi$
discussed in § 5. As shown in figure 22(a), increasing the viscosity ratio from
$\lambda = 1$
to
$\lambda = 5$
barely changes the average deformation
$\langle D \rangle$
of the suspension. However, the inclination angle
$\langle \theta \rangle$
and the magnitude of the tank-treading velocity
$\langle \omega \rangle$
are reduced due to the increased resistance in the internal fluid of the capsules. Regarding particle stresses, we observe a slight increase in the relative viscosity
$\mu _r$
in figure 22(d), which is opposite to the effect of
$\lambda$
in the single capsule case, as depicted in figure 8(c). This increase can be attributed to the enhanced short-range interactions among the suspended capsules due to the higher internal viscosity at a higher
$\lambda$
. In terms of the normal stress differences,
$\lambda$
shows only a minor impact on the evolution of
$N_1^p$
in figure 22(e), although a slight decrease in
$N_2^p$
is observed at a higher viscosity ratio
$\lambda = 5$
in figure 22(f). In the non-inertial regime, we observe similar effects of
$\lambda$
, so the numerical results are not presented for brevity.
Appendix F. Fluctuations of the shear stress
Shear stress fluctuation of the capsule suspension. (a-c) Temporal evolution of the shear stress under effects of
$\mathcal {C}a$
,
$\phi$
and
$\mathcal {R}e$
, respectively; (d-e) probability distribution function
$f$
of particle shear stress fluctuations at
$\phi =0.1$
and
$\phi =0.4$
; a reference normal distribution is presented in black solid line. Here,
$\sigma _{\Sigma _{xy}^p}$
denotes the standard deviation of the particle shear stress
$\Sigma _{xy}^p$
.
Due to the inhomogeneous microscopic structure of the capsule suspension, the particle shear stress consistently fluctuates around the average value, once the flow is fully developed
\begin{equation} \Sigma _{xy}^p = \overline {\Sigma }_{xy}^p + \delta \Sigma _{xy}^p ,\end{equation}
where
$\overline {\Sigma }_{xy}^p$
and
$\delta \Sigma _{xy}^p$
denote the time-averaged shear stress and its fluctuations, respectively. Figure 23(a–c) presents the temporal evolution of the shear stress fluctuations at different
$\mathcal {C}a$
,
$\phi$
and
$\mathcal {R}e$
. We observe that higher deformability of the capsules leads to increased shear stress fluctuations
$\delta \Sigma _{xy}^p$
at
${\mathcal {R}e} = 10$
and
$\phi =0.4$
in figure 23(a). In the inertial regime,
$\phi$
has minor effects on
$\delta \Sigma _{xy}^p$
as shown in figure 23(b). However, an increase in
$\mathcal {R}e$
results in a significant rise in shear stress fluctuations in figure 23(c). Generally,
$\delta \Sigma _{xy}^p$
represents less than 15 % of the total shear stress
$\Sigma _{xy}^p$
even at
${\mathcal {R}e} = 20$
.
Regarding the distribution of stress fluctuations, figure 23(d) shows that at low Reynolds number (
${\mathcal {R}e} = 1$
), low volume fraction (
$\phi = 0.1$
) and low capillary number (
${\mathcal {C}a} = 0.01$
),
$\delta \Sigma _{xy}^p$
follows a normal distribution with an average value of zero. Increasing
$\mathcal {C}a$
causes deviations from this normal distribution, leading to a lower peak value as seen in figure 23(d). In a dense suspension (
$\phi = 0.4$
), there is a clear deviation from the normal distribution for all investigated
$\mathcal {C}a$
values, as shown in figure 23(e). Higher
$\mathcal {R}e$
also causes deviations from the normal distribution, although these results are not shown here for brevity.
Open access
