Injection of Deformable Capsules in a Reservoir: A Systematic Analysis

Abstract

Motivated by red blood cell dynamics and injectable capsules for drug delivery, in this paper, a computational study of capsule ejection from a narrow channel into a reservoir is undertaken for a combination of varying deformable capsule sizes and channel dimensions. A mass-spring membrane model is coupled to an Immersed Boundary–Lattice Boltzmann model solver. The aim of the present work is the description of the capsules’ motion, deformation and the response of the fluid due to the complex particles’ dynamics. The interactions between the capsules affect the local velocity field and are responsible for the dynamics observed. Capsule membrane deformability is also seen to affect inter-capsule interaction. We observe that the train of three particles locally homogenises the velocity field and the leading capsule travels faster than the other two trailing capsules. Variations in the size of reservoir do not seem to be relevant, while the ratio of capsule diameter to channel diameter as well as the ratio of capsule diameter to inter-capsule spacing play a major role. This flow set-up has not been covered in the literature, and consequently we focus on describing capsule motion, membrane deformation and fluid dynamics, as a preliminary investigation in this field.

1. Introduction

Haemodynamcis in large arteries is commonly described by the incompressible Newtonian Navier–Stokes equations, hence modelling whole blood to have a constant density and viscosity. While this is acceptable for larger arteries at high flow rates, it is more appropriate to adopt a thixotropic non-Newtonian shear-thinning rheological model for viscosity when a larger variation of shear rates is apparent as evident in smaller vessels or slower flows. This would then take into account the presence of the erythrocytes and other constituents of whole blood [1,2] in a continuum model. However, when smaller vessels of the cardiovascular system are considered the dimension of the conduits and the circulating cells are of similar scales, it is therefore necessary to discretise and model whole blood as a multi-component medium. At these smaller scales, the properties of the cells (such as material constitutive laws for the membrane), the inter-cellular flow interactions (such as flow wakes) and other biochemical and biological phenomena (such as tethering or remodelling), must be considered and together describe a complex physical interplay.

Experimental works have been at the forefront and driving much of the research and our understanding of haemodynamic micro-circulation for many years [3], and it is only with the advance of commodity computational resources that new numerical methods have been developed and with them computational simulations of micro-circulation have been possible. Numerical simulations provide a fine spatial and temporal resolution of physical variables, which enable for a quantitative analysis and allow for different mathematical models and hypotheses to be tested. Some fundamental studies of cells or capsules have been undertaken using computational simulations, investigating the importance of cell shape and deformability, concentration and apparent viscosity, transport and migration, providing important insight into micro-circulation dynamics [4,5,6,7,8,9,10,11].

Research in the field of blood micro-circulation has seen a range of applications and interests. For example, specially designed micro-channel geometries have been used to separate or sort suspended cells. Such designs include simple, sudden expansions which promote cell focusing [12,13,14], or alternatively repeated sections of hyperbolic micro-channels [15], wavy channels [16], guiding grooves [17], multi-stage micro-fluidic devices involving bends and siphoning [18,19], though a range of different micro-fluidic device configurations exist [20,21]. These largely make use of inertial forces of the cells [22,23], as well as cell deformability [7,13,24,25]. Interestingly, however, micro-channels may also be designed in a very similar fashion to enhance mixing of the flow [26], and a review of low cost fabrication devices is presented in [14].

While inertial effects have been predominantly used to sort cells in micro-fluidic devices, it is known that cell deformability and shape play important roles in their transport dynamics [7,8,19,27,28,29,30]. The volumetric concentration of suspended particles in flow is also known to affect the apparent viscosity [5,9,11] of the medium, and the resulting inter-cellular flow interactions have been observed to affect transport of the cells through different micro-channel configurations [31,32,33,34,35]. The motion of suspended particles in micro-channels are also known to induce a pattern of wall shear stress variation along the wall [31,36,37], which is not only important in mechanotransduction and signalling pathways, but also in cell adhesion mechanics [38]. The effect of particle suspensions of different sizes has also been investigated, with relevance to leukocyte radial margination [39] and micro- and nanoparticles on drug delivery [40,41].

In the present work, we investigate ejection of capsules from a narrow channel to a reservoir, comparing different size ratios of channel and capsule diameters. Specifically, the aim of the present work is to detail the dynamics of circular capsules when navigating across a geometric discontinuity. The presence of capsules dragged by the flow locally increases the apparent viscosity of the fluid in the region immediately close to and inside the membrane itself. This causes the homogenisation of the velocity field disturbed by the presence of the particles, and observe that a train of three particles will tend to act as a single larger body (because of the inter-particle interaction). Finally, the effect of the local increased viscosity also causes the leading capsule to move faster than the other two trailing capsules. We perform numerical simulations of micro-fluidic particulate flow of deformable capsules in discontinuous geometries, with relevance to capsule injection in applications such as drug delivery. This flow set-up has not been covered in the literature, and consequently we focus on describing capsule motion, membrane deformation and fluid dynamics, as a preliminary investigation in this field.

2. Computational Method

Computational methods to model and solve for multi-component micro-circulation have developed immensely in the last decades. Numerical methods which discretise the domain as lumped volumes (or masses) of fluid, typically denoted as particle methods have been popular, including dissipative particle dynamics (DPD), smoothed particle hydrodynamics (SPH), moving particle semi-implicit method (MPS), and multiparticle collision dynamics (MCP) [42,43,44,45,46,47]. These methods are based on expressing the governing equations in a moving reference frame, which is well suited to flows with deformable bodies and moving boundaries. Here, we adopt a mixed approach, in which the fluid is solved on a fixed grid, while the capsule membranes are described in a moving reference frame. The solution to the membrane forces is then interpolated to the fixed grid. In doing so, we adopt an immersed boundary method, and employ the Lattice Boltzmann method as the fluid solver.

2.1. Lattice Boltzmann Method

The evolution of the fluid is defined in terms of a set of N discrete distribution functions, $[f_i],(i=0,\cdots ,N-1)$, which obey the dimensionless Boltzmann equation

$$f_i(\overrightarrow{x}+{\overrightarrow{e}}_i\Delta t,t+\Delta t)-f_i(\overrightarrow{x},t)=-\frac{\Delta t}{\tau }[f_i(\overrightarrow{x},t)-f_i^{eq}(\overrightarrow{x},t)],$$

(1)

in which $\overrightarrow{x}$ and t are the spatial and time coordinates, respectively; $[{\overrightarrow{e}}_i],(i=0,\dots ,N-1)$ is the set of N discrete velocities; $\Delta t$ is the time step; and $\tau$ is the relaxation time given by the unique non-null eigenvalue of the collision term in the BGK-approximation [48]. The kinematic viscosity of the flow is strictly related to $\tau$ as $\nu =c_s^2(\tau -\frac{1}{2})\Delta t$ being $c_s=\frac{1}{\sqrt{3}}\frac{\Delta x}{\Delta t}$ the reticular speed of sound. The moments of the distribution functions define the fluid density $\rho ={\sum }_i{f}_i$, velocity $\overrightarrow{u}={\sum }_i{f}_i{\overrightarrow{e}}_i/\rho$, and pressure $p=c_s^2\rho =c_s^2{\sum }_i{f}_i$. The local equilibrium density functions $[f_i^{eq}](i=0,\dots ,N-1)$ are expressed by the Maxwell–Boltzmann distribution:

$$f_i^{eq}(\overrightarrow{x},t)={\omega }_i\rho \left[1+\frac{1}{c_s^2}({\overrightarrow{e}}_i·\overrightarrow{u})+\frac{1}{2c_s^4}{({\overrightarrow{e}}_i·\overrightarrow{u})}^2-\frac{1}{2c_s^2}(\overrightarrow{u}·\overrightarrow{u})\right].$$

(2)

On the two-dimensional square lattice with $N=9$ speeds $\left(D2Q9\right)$ [49], the set of discrete velocities is given by

$${\overrightarrow{e}}_i=\left\{\begin{array}{cc}(0,0),\hfill & ifi=0,\hfill \\ \left(cos\left(\frac{(i-1)\pi }{2}\right),sin\left(\frac{(i-1)\pi }{2}\right)\right),\hfill & ifi=1-4,\hfill \\ \sqrt{2}\left(cos\left(\frac{(2i-9)\pi }{4}\right),sin\left(\frac{(2i-9)\pi }{4}\right)\right),\hfill & ifi=5-8,\hfill \end{array}\right.$$

(3)

with the weight, ${\omega }_i=1/9$ for $i=1-4$, ${\omega }_i=1/36$ for $i=5-8$, and ${\omega }_0=4/9$. Here, we adopt a discretisation in the velocity space of the equilibrium distribution based on the Hermite polynomial expansion of this distribution [50].

2.2. Immersed Boundary Treatment

Deforming body models are commonly based on continuum approaches using strain energy functions to compute the membrane response [51,52,53]. However, a particle-based model governed by molecular dynamics has emerged due to its mathematical simplicity while providing consistent predictions [10,54,55,56]. In this work, a particle-based model is employed by coupling the Immersed-Boundary (IB) technique with BGK-lattice Boltzmann solver. The immersed body is a worm-like chain of $nv$ vertices linked with $nl$ linear elements, whose centroids are usually called Lagrangian markers. A forcing term $[{\overrightarrow{\mathcal{F}}}_i](i=0,\dots ,8)$, accounting for the immersed boundary, is included as an additional contribution on the right-hand side of Equation (1):

$$f_i(\overrightarrow{x}+{\overrightarrow{e}}_i\Delta t,t+\Delta t)-f_i(\overrightarrow{x},t)=-\frac{\Delta t}{\tau }[f_i(\overrightarrow{x},t)-f_i^{eq}(\overrightarrow{x},t)]+\Delta t{\overrightarrow{\mathcal{F}}}_i.$$

(4)

${\overrightarrow{\mathcal{F}}}_i$ is expanded in terms of the reticular Mach number, $\frac{{\overrightarrow{e}}_i}{c_s}$, resulting in:

$${\overrightarrow{\mathcal{F}}}_i=\left(1-\frac{1}{2\tau }{\omega }_i\left[\frac{{\overrightarrow{e}}_i-\overrightarrow{u}}{c_s^2}+\frac{{\overrightarrow{e}}_i·\overrightarrow{u}}{c_s^4}{\overrightarrow{e}}_i\right]\right)·{\overrightarrow{f}}_{ib},$$

(5)

where ${\overrightarrow{f}}_{ib}$ is a body force term. Due to the presence of the forcing term, the mass density and the momentum density are derived as $\rho ={\sum }_i{f}_i$ and $\rho \overrightarrow{u}={\sum }_i{f}_i{\overrightarrow{e}}_i+\frac{\Delta t}{2}{\overrightarrow{\mathcal{F}}}_i$.

Within this parametrisation, the forced Navier–Stokes equations is recovered with a second order accuracy [57,58,59,60,61]. The external boundaries of the computational domain are treated with the known velocity bounce back conditions by Zou and He [62]. The IBM procedure, extensively proposed and validated by Coclite and colleagues [27,28,29,63,64], is adopted here and the moving-least squares reconstruction by Vanella et al. [65] is employed to exchange all LBM distribution functions between the Eulerian lattice and the Lagrangian chain. Finally, the body force term in Equation (5), ${\overrightarrow{f}}_{ib}$, is evaluated through the formulation by Favier et al. [66].

Elastic Membrane deformation. Elastic membranes are modelled by means of an elastic strain, bending resistance, and total enclosed area conservation potentials. Specifically, the nodal forces corresponding to the elastic energy for nodes 1 and 2 connected by edge l reads as

$$\left\{\begin{array}{cc}{\overrightarrow{F}}_1^s\hfill & =-k_s(l-l_0)\frac{{\overrightarrow{r}}_{1,2}}{l},\hfill \\ {\overrightarrow{F}}_2^s\hfill & =-k_s(l-l_0)\frac{{\overrightarrow{r}}_{2,1}}{l},\hfill \end{array}\right.$$

(6)

where ${\overrightarrow{r}}_{i,j}={\overrightarrow{r}}_i-{\overrightarrow{r}}_j$ with $r_i$ position vector of the node i.

The bending resistance related to the v-th vertex connecting two adjacent element is

$$V_v^b=\frac{1}{2}{k}_b{(k_v-k_{v,0})}^2,$$

(7)

$k_b$ being the bending constant, $k_v$ the current local curvature in the v-th vertex, and $k_{v,0}$ the local curvature in the v-th vertex for the stress-free configuration. The curvature is evaluated by measuring the variation of the angle between two adjacent elements ($\theta -{\theta }_0$), with ${\theta }_0$ the angle in the stress free configuration. Given this, the forces on the nodes $v_{left}$, v, and $v_{right}$ are obtained as

$$\left\{\begin{array}{cc}{\overrightarrow{F}}_{v_{left}}^b\hfill & =k_b(\theta -{\theta }_0)\frac{l_{left}}{l_{left}+l_{right}}{\overrightarrow{n}}_v,\hfill \\ {\overrightarrow{F}}_v^b\hfill & =-k_b(\theta -{\theta }_0){\overrightarrow{n}}_v,\hfill \\ {\overrightarrow{F}}_{v_{right}}^b\hfill & =k_b(\theta -{\theta }_0)\frac{l_{right}}{l_{left}+l_{right}}{\overrightarrow{n}}_v,\hfill \end{array}\right.$$

(8)

where $l_{right}$ and $l_{left}$ are the length of the two adjacent left and right edges, respectively, and ${\overrightarrow{n}}_v$ is the outward unity vector centred in v. Note that, in this context, the relation between the strain response constant $k_s$ and $k_b$ is $E_b=\frac{k_b}{k_s{r}^2}$, where r is the particle radius.

In order to limit the membrane stretching, an effective pressure force term is considered. Thus, the penalty force is expressed in terms of the reference pressure $p_{ref}$ and directed along the normal inward unity vector of the l-th element $n_l^{-}$, as

$${\overrightarrow{F}}_l^a=-k_a(1-\frac{A}{A_0})p_{ref}{\overrightarrow{n}}_l^{-}{l}_l,$$

(9)

with $l_l$ the length of the selected element, $k_a$ the incompressibility coefficient, A the current enclosed area, $A_0$ the enclosed area in the stress-free configuration. The enclosed area is computed using the Green’s theorem along the curve, $A=\int x_{l}d{y}_l-y_{l}d{x}_l$. Within this formulation, $k_a=1$ returns a perfectly incompressible membrane. Note that ${\overrightarrow{F}}_l^a$ is evenly distributed to the two vertices connecting the l-th element ($v_{left}$ and $v_{right}$) as ${\overrightarrow{F}}_l^a=0.5{\overrightarrow{F}}_{v_{left}}^a+0.5{\overrightarrow{F}}_{v_{right}}^a$.

Particle–Particle Interaction Two-body interactions are modelled with a repulsive potential centred in each vertex. The purely repulsive force is such that the minimum allowed distance between two vertices coming from two different particles is $\Delta x$. The impulse acting on vertex 1, at a distance $d_{1,2}$ from the vertex 2 of an adjacent particle, is directed in the inward normal direction identified by ${\overrightarrow{n}}_1^{-}$ and is given by:

$${\overrightarrow{F}}_1^{pp}=\frac{{10}^{-4}}{8\sqrt{2}}\sqrt{\frac{\Delta x}{d_{1,2}^5}}{\overrightarrow{n}}_1^{-}.$$

(10)

Hydrodynamics Stresses Pressure and viscous stresses exerted by the l-th linear element are:

$${\overrightarrow{F}}_l^p\left(t\right)=(-p_l{\overrightarrow{n}}_l)l_l,$$

(11)

$${\overrightarrow{F}}_l^{\tau }\left(t\right)=({\overline{\tau }}_l·{\overrightarrow{n}}_l)l_l,$$

(12)

where ${\overline{\tau }}_l$ and $p_l$ are the viscous stress tensor and the pressure evaluated in the centroid of the element, respectively; ${\overrightarrow{n}}_l$ is the outward normal unit vector while $l_l$ is its length. The pressure and velocity derivatives in Equations (11) and (12) are computed using a probe in the normal positive direction of each element, the probe length being $1.2\Delta x$ [65,67].

2.3. Fluid–Structure Interaction

Particles dynamics are determined by dynamics IB technique described in [64], using the solution of the Newton equation for each Lagrangian vertex, accounting for both internal, Equations (6) and (8)–(10), and external stresses, Equations (11) and (12). Then, no-slip boundary conditions are imposed using a weak coupling approach [63]. The total force ${\overrightarrow{F}}_v^{tot}\left(t\right)$ acting on the v-th element of the immersed body is evaluated in time and the position of the vertices is updated at each Newtonian dynamics time step considering the membrane mass uniformly distributed over the $nv$ vertices,

$$m_v{\dot{\overrightarrow{u}}}_v={\overrightarrow{F}}_v^{tot}\left(t\right)={\overrightarrow{F}}_v^s\left(t\right)+{\overrightarrow{F}}_v^b\left(t\right)+{\overrightarrow{F}}_v^a\left(t\right)+{\overrightarrow{F}}_v^{pp}\left(t\right)+{\overrightarrow{F}}_v^p\left(t\right)+{\overrightarrow{F}}_v^{\tau }\left(t\right).$$

(13)

The Newton equation of motion is integrated by using the Verlet algorithm. Specifically, a first tentative velocity is considered into the integration process, ${\dot{\overrightarrow{x}}}_{v,0}\left(t\right)$, obtained interpolating the fluid velocity from the surrounding lattice nodes

$${\overrightarrow{x}}_v(t+\Delta t)={\overrightarrow{x}}_v\left(t\right)+{\dot{\overrightarrow{x}}}_{v,0}\left(t\right)\Delta t+\frac{1}{2}\frac{{\overrightarrow{F}}_v^{tot}\left(t\right)}{m_v}\Delta t^2+O(\Delta t^3).$$

(14)

Then, the velocity at the time level $t+\Delta t$ is computed as

$${\overrightarrow{u}}_v(t+\Delta t)=\frac{\frac{3}{2}{\overrightarrow{x}}_v(t+\Delta t)-2{\overrightarrow{x}}_v\left(t\right)+\frac{1}{2}{\overrightarrow{x}}_v(t-\Delta t)}{\Delta t}+O(\Delta t^2).$$

(15)

It should be noted that the present formulation is unconditionally stable for a small deformation of the capsule membrane and for small velocity variations applied, as previously demonstrated by the authors [27,63,64].

2.4. Set-Up and Boundary Conditions

The simulations are performed for a two-dimensional domain as shown in Figure 1 and the fluid is considered to be water. The flow direction is left-to-right, the horizontal axis is denoted by the x-axis or co-axial direction, and the vertical axis is denoted by the y-axis or radial direction. The analysis is based on simulations either a single or three in-line spherical capsules, flowing from a channel of small diameter into that of a larger diameter, as shown in Figure 1. In subsequent discussion and presentation of results, we refer to the upstream direction as that closer to the inflow, and the downstream direction that closer to the outflow. One should, however, recognise that the capsule will be travelling faster than the bulk flow in the channel sections, since it is located furthest from the stationary walls. Consequently, in a moving reference frame following a capsule, its wake and disturbance it induces on the flow will in effect be in the upstream direction.

The computational domain is a rectangular channel by the height l and length $3l$ with no-slip boundary conditions along y. The grid resolution is such that l is discretised with $500\Delta x$. The section $x=0$ is parametrised as a velocity inlet section with a parabolic inlet profile, while the outlet section is located at $x=3l$ (see Figure 1) and is set as a convective condition [68]. The Reynolds number is fixed equal to 0.1 and is given by $Re=\frac{u_{max}d}{\nu }$; where $u_{max}$ is the maximum velocity for the plane Hagen–Poiseuille profile (parabolic profile) established in $x=0$, d is the narrow channel diameter, and $\nu$ the kinematic viscosity of the fluid ($\nu =1.2\times {10}^{-6}$ m${}^2$/s). Note that $Re$ represents the ratio of inertial to viscous forces and can also be interpreted as the ratio of viscous to convective time scales which act on the fluid. Here, the Reynolds number equals 0.1; consequently, viscous effects dominate, with viscous forces greater than inertial forces and the viscous time scale is smaller (hence acts faster, stabilising the flow) than the convective time scale. The capsules are initially at the rest, with no pre-stress applied, and of circular section with diameter r and stiffness modulated by the capillary number $Ca={10}^{-3}=\frac{\rho \nu u_{max}}{k_s}$, where $\rho$ is the density of the fluid in which capsules are immersed ($\rho =1000$ kg/m${}^3$) in and $k_s$ is the elastic constant used for the worm-like chain composing the membranes (see Equation (6)). $Ca$ represents the ratio of the viscous force to the elastic force consequently representing, within this definition, capsules slightly more rigid than in other studies [7,25]. Note that the typical stiffness for a red blood cell is $k_s=1.5\times {10}^{-4}$ Nm that would lead to Ca${}_{rbc}={10}^{-2}$ within this definition. To ensure the scheme stability, all the computations are performed with $\tau =1.0$ in the Lattice Boltzmann Method.

The present parametrisation deals with the deformation of a circular membrane in a rectangular two-dimensional channel injected into a reservoir. Given the narrow dimensions and axial symmetry of the problem and set-up, the results obtained are directly transferable to the analogous three-dimensional set-up. For a general case, the deformation and dynamics of initially spherical capsules in a three-dimensional circular capillary flowing into a reservoir will yield different quantitative results. However, the governing physics is unaltered and, while the mechanics in two or three dimensions is different, results from a two-dimensional investigation are transferable to a three-dimensional set-up and one can expect similar qualitative results and trends.

3. Results and Discussion

Numerical simulations for flow set-up outlined in Figure 1 were run for the following three cases: without capsules; with one capsule; with three capsules (numbered left-to-right). The available set-up combinations have resulted in a set of simulations, aimed at sampling the solution space in order to capture the physics of flow of capsules as they are injected into a reservoir. Videos of the simulations run are available online as Supplementary Files, organised by the capsule to channel ratio $r/d$.

The solution for flow without any capsules are shown in Figure 2 and Figure A1 (in Appendix A), for purpose of comparison. As expected, we see the flow profile develop from the parabolic inflow profile to a flattened paraboloid profile on approaching the reservoir. The flow accelerates in the radial (vertical) direction and decelerates in the co-axial (horizontal) direction as it approaches the end of the smaller channel before discharging into the reservoir. The radial acceleration is effected by a pressure gradient which drives the flow to turn at the geometric discontinuity. The flow separates at the geometric discontinuity, reattaching onto the horizontal walls of the larger channel.

We now turn our attention to the ejection of a single capsule. While the results and discussion are transferable to the other set-up combinations, we present results for the geometry set-up $r/d=0.50$ and $l/d=5.0$ in Figure 3 as these favour presentation and discussion. In Figure 3, we present the (normalised) co-axial flow component, as well as the difference in flow velocity without capsules and with one capsule (coloured by co-axial velocity component). In the following presentation of results, this velocity difference is termed the relative velocity or relative flow field.

Let us first focus on results for time snapshot $t{u}_{max}/d=5$ in Figure 3, when the capsule is located in the narrower channel. We observe that the presence of the capsule tends to create a more uniform velocity profile, by increasing the co-axial velocity near to the walls while decreasing it at the centre of the channel. The vector plot showing the relative velocity also shows that the presence of the flowing capsule results in flow around the capsule, from upstream to downstream. This is because the capsule in effect creates a resistance to the flow, causing it to flow more around it in comparison to the parabolic velocity profile. The relative flow field consequently appears as a vortex pair (or ring in three-dimensions) close to the wall, travelling with the capsule and aligned approximately at its maximum diameter. This vortex ring, identified in the relative flow field, induces a co-axial velocity which drives the capsule through the channel, and this can be seen by the higher velocity at the maximum diameter of the capsule. The wall shear stress is the tangential traction force that the flow exerts on the wall due to viscous effects, and we observe that, at the capsule’s maximum diameter, the co-axial velocity near the channel wall is the same as the parabolic flow, and hence the wall shear stress is effectively unaltered. However, slightly upstream and downstream of the capsule’s maximum diameter, we observe that the co-axial velocity near the wall is greater than for the parabolic profile, indicating a higher wall shear stress. This results in two locations (or rings in three-dimensions) of higher wall shear stress compared to the parabolic flow solution, travelling with the capsule and approximately aligned with the capsule’s leading and trailing edges. This is in contract to the results of the wall shear stress `footprint’ observed as red blood cells, which show a single higher band of wall shear stress as the deformable cells are in proximity to the walls [31,36,37].

We now turn our attention to results for time snapshot $t{u}_{max}/d=10$ in Figure 3, when the capsule has been ejected into the reservoir. We observe that, for this instance, the mean velocity of the capsule is the same as a solution when no capsule is present, while, in comparison, the flow is moving faster in front of the capsule and slower behind. We also observe that the relative flow field presents a vortex ring at the trailing side of the capsule, travelling with the capsule. This vortex ring is counter-rotating to the vortex ring observed in the narrower channel section, and is set up by the geometric discontinuity. The direction of rotation of this vortex ring induces a velocity that promotes the flow to turn around the geometric discontinuity and results in a smaller flow separation. The induced velocity of this vortex ring also acts to decelerate the capsule co-axial motion as it ejects into the reservoir.

Finally, we note that, for time snapshot $t{u}_{max}/d=20$, the flow of a single capsule in a large channel or reservoir has no marked influence on the flow field as compared to the above discussed time snapshots $t{u}_{max}/d=$ 5 and 10.

Similar flow fields and relative flow fields were observed in the other set-up combinations. We turn our attention therefore to the motion of the capsule and its change in perimeter length, presented in Figure 4, and investigate the effects of the geometric variations based on the ratio $l/d$ and $r/d$. Overall, the trends appear linear as one varies the geometric ratios $l/d$ and $r/d$; however, there are some small deviations that are worth highlighting and discussing.

For comparison purposes in this figure, the case for $l/d=1.0$ (straight channel without discontinuity) is also plotted. We observe that the capsule velocity and perimeter length change are constant (after rapid adjustment from the zero velocity initial conditions of the simulations); however, there are variations based on the $r/d$ ratio. At larger values of $r/d$, the capsule forms a larger blockage in the channel, and we observe that a more blunt velocity profile is caused by the $r/d=0.75$ capsule compared to the $r/d=0.25$ capsule. The greater perimeter length change is also observed for the $r/d=0.75$ capsule, since there are larger shear rates in proximity with the channel wall aligned to the capsule’s largest diameter as it flows. This was observed in Figure 3 from the relative flow field.

We now focus on the channel geometry set-up with $l/d=$ 2.5 and 5.0 within Figure 4. We observe that the traces of normalised co-axial position $x_{cm}/d$ for capsules with $r/d=$ 0.25 and 0.50 lie closer than for capsule $r/d=0.75$. The lag is seen to appear at the moment when the capsule is injected into the reservoir, with $t{u}_{max}/d\approx 7.5$. We observe that the traces of normalised co-axial velocity $u_{cm,x}/u_{max}$ are indistinguishable while the capsule lies within the narrower channel, and parallel based on their $r/d$ value. The capsule velocity then transits to the similar values once the capsule enters the reservoir occurs rapidly, now based on the value of $l/d$, while the ratio $r/d$ has little effect. For $l/d=5.0$ and $r/d=0.75$, we observe that there is a secondary peak in the velocity at $t{u}_{max}/d\approx 8$. We observe that the perimeter length change $\delta p$ is relatively constant as the capsule travels along the narrower channel but decreases and then increases sharply as the capsule is ejected into the reservoir. The decrease in perimeter length is due to the capsule leading edge slowing down as it reaches the reservoir, resulting in a decrease in membrane stress. The subsequent increase in perimeter length occurs as the capsule completes its transition into the reservoir, during which the anterior portion of the capsule is already in the reservoir and has a low velocity; however, the posterior portion of the capsule still has a larger velocity, causing the capsule to flatten (stretching radially). Once in the reservoir, the capsule membrane relaxes and tends to assume an undeformed shape. The change in perimeter length is larger in the narrower channel due to higher shear rates, more so with an increasing $r/d$ ratio, which was also observed from the relative flow field shown in Figure 3. For $l/d=5.0$ and $r/d=0.75$, we observe that the sudden decrease and subsequent increase in perimeter lengths were in proportions higher than other cases.

Focusing on the case with $l/d=5.0$ and $r/d=0.75$, we summarise that we observed a different behaviour as the capsule was ejected into the reservoir, compared to the other simulations. This led to a lag in its co-axial position, a second peak in the co-axial velocity and more pronounced change in the perimeter length. The reasons for these phenomena are principally due to the size of the capsule, which owing to the membrane have the effect of locally constraining the flow to be more uniform (hence a homogeneous velocity field). Capsule deformability and the elastic forces are also important, without which we would not obtain the second velocity peak, for example.

We now turn our attention to the results for three capsules flowing in the channel and ejecting into the reservoirs. The flow field and relative velocity for the geometry set-up $r/d=0.50$ and $l/d=5.0$ with three capsules are presented in Figure 5. In comparison to the flow of the single capsule (see Figure 3), we see that, within the narrow channel, the disturbance in the flow field extends to influence both upstream and downstream capsules, due to their relative proximity. Specifically, the capsules effect a blunter velocity profile with respect to the parabolic profile (obtained without capsules). This phenomenon is often described by a shear-thinning non-Newtonian model for the viscosity. At time $t{u}_{max}/d=10$, when the capsules are entirely in the reservoir, the leading capsules (capsules 2 and 3) are travelling faster than is the case without capsules, while the left-most capsule (capsule 1) is travelling slower. This is again due to the capsule acting to locally constrain the flow to be more uniform, and in this case extending the apparent jet formed by the flow ejecting from the narrower channel into the reservoir. At this time instance, we also note that a vortex ring is formed at the trailing section of capsule 1 only, further highlighting that the flow disturbances between the three capsules result in a local homogenisation of the velocity field, and the capsules therefore tend to act as a single larger object.

In Figure 6, the properties of the single capsule (also shown in Figure 4) and of the multiple particles are plotted as function of normalised time. We notice that the effect of the relative reservoir diameter, $l/d$, does not play a noticeable influence beyond a given size. Additionally, we observe a greater effect of the multiple capsules when they are larger, and therefore focus our presentation of results for the set-up $l/d=5.0$ and $r/d=0.75$. The discussion is amenable to the other set-up cases, and differences are presented. In general, for the smaller capsules with $r/d=$ 0.50 and 0.25 the same phenomena are present as for the larger capsules, but to a lesser extent. Indeed, for the smallest capsules $r/d=0.25$, the effect of the multiple cells is almost imperceivable. The reason for this is a reduced inter-capsule interaction, due to the relatively large capsule separation distance and the small capsule size that will not significantly affect the flow field.

In Figure 6, observing the traces for set-up $l/d=5.0$ and $r/d=0.75$, we first note that the capsules exhibit oscillating velocities and change in perimeter, more marked for capsule 3 and least for capsule 1. These oscillations are related to the ejection of capsules into the reservoir and we see, for example, that, as capsule 1 enters the reservoir, it induces an oscillation in both capsules 2 and 3. This chain effect is due to the capsules effectively constraining the fluid, and locally homogenising the velocity field. The apparent viscosity is locally higher, and the capsules, due to their close proximity, effectively act as a single larger body. Another phenomenon of particular interest is that capsule 3 tends to move overall faster, and, from the $x/d$ trace, we see that it is farther from capsules 1 and 2 towards later times. In order to explain this, we turn to Figure 5. We observe that the relative velocity is greatest ahead of capsule 3 when the capsules are in the reservoir (at time $t{u}_{max}/d=10$) and the flow is in the narrower channel a short distance ahead of capsule 3 is undisturbed and parabolic. Since capsule 3 is at the leading edge of the capsules, it is not affected by the wake and will be able to travel faster. This phenomenon of capsule spacing rearrangement was also observed in [16], though in very different geometries. Finally, we note that the perimeter stretching is greatest overall for capsule 1, due to its location in the end of the capsule train, inducing a more marked wake and consequently higher shear rate (velocity gradients) in the fluid and higher stresses in its membrane. In fact, we observe that the change in perimeter length is comparable to the set-up of the single capsule, since it has a similar wake flow field.

4. Conclusions

In this work, we investigate the dynamics of capsule ejection from a narrow channel into a reservoir, across a geometric discontinuity. We observe that inter-capsule interaction (due to the wakes of their motion) and the constraining of the fluid within the membranes are an important mechanism that affects the local apparent viscosity since the stress field must be continuous in the domain.

In order to span a meaningful parameter space, a combination of different configurations was investigated. The capsules were varied to have different sizes, namely $r/d=$ 0.25, 0.50 and 0.75, where r is the capsule diameter and d is the narrow channel diameter. Three configurations of channel geometries were investigated, namely $l/d=$ 1.0, 2.5 and 5.0, where l is the diameter of the reservoir. Additionally, three different configurations: no capsule, a single capsule, and three in-line capsules, were simulated and investigated. The Capillary number and Reynolds numbers were chosen to be $Ca={10}^{-3}$ and $Re={10}^{-1}$.

The simulations were investigated by observing the relative flow field that is the flow field resulting from capsule flow as compared to the no capsule solutions. This has proved to be an effective means of identifying where the flow field has altered, and, consequently, to identify the fluid mechanics phenomena causing the changes observed. Additionally, the trajectories, velocities and perimeters of the capsules were tracked during the simulations.

Overall, we have seen that the reservoir diameter has a negligible effect beyond a threshold, and, in the resent investigation, similar results were obtained for $l/d=2.5$ and $l/d=5.0$. The effect of capsule size was seen to be have a greater effect, with $r/d=0.75$ unsurprisingly resulting in the greatest deviation from a flow field with no capsule, however capsules with size $r/d=0.25$ were also seen to affect the flow field.

Capsule membranes constrain the flow internally, and since the stress field must be continuous across the capsule membrane, the effect is to locally homogenise (i.e., create greater uniformity) the velocity field. This can be seen as a local increase in apparent viscosity. When multiple capsules were investigated, the inter-capsule interaction caused the capsules to effectively act as a single larger body. This resulted in an increased apparent viscosity spanning the region of the capsules. This effect was clearly observed as the capsules flow in the narrow channel, for which the apparent viscosity resembled that of a shear-thinning non-Newtonian rheological model. An effect of the local increased viscosity is also the cause that the leading capsule tends to move faster than the other two trailing capsules.

The effect of the multiple capsules is to reduce the perimeter change, due to their wakes and inter-capsule interaction which reduces the shear rate (i.e., velocity gradients) of the fluid integrated over the capsule surface. This then leads to a decrease in overall strain for the capsule membrane. The capsule at the trailing edge however is not shielded and its wake promotes a vortex ring in the relative velocity field, and its perimeter change is the same as that of a single capsule flow.

Lastly, we highlight that, while the two-dimensional results reported here are representative of the analogous three-dimensional problem, due to the symmetry and regimes (based on capillary and Reynolds numbers) of the set-up, this is generally not the case. Indeed, in complex systems, such as a general set-up where deformable capsules are injected into a reservoir, not only are the mechanics of the jet collapse different between two and three dimensions, importantly also the specific stresses involved in the fluid–structure interactions will differ. This noted, two-dimensional simulations can still provide fruitful information on the regulating biophysical mechanisms without the inconvenience of the computationally intense three-dimensional simulations. The extension of the current physical problem to three-dimensional modelling is certainly of interest and will be the object of future investigations.