Injection of Deformable Capsules in a Reservoir, a Systematic Analysis

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 particle dynamics. The interactions between the capsules affect the local velocity field significantly and are responsible for the dynamics observed. Capsule membrane deformability is also seen to affect inter capsule interaction, and we observe that the train of three particles locally homogenizes the velocity field and the leading capsule travels faster than the other two trailing capsules. On the contrary, variations in size of the reservoir do not seem to be relevant, while the ratio of capsule diameter with respect to channel diameter plays a major role as well as the ratio of capsule diameter to inter capsule spacing. 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.


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 are apparent as evident is 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 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 [24,7,13,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, 27-30, 8, 19]. The volumetric concentration of suspended particles in flow is also known to affect the apparent viscosity [9,5,11] of the medium, and the resulting inter-cellular flow interactions have been observed to affect transport of the cells through different microchannel 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 [36,37,31], 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 nano-particles 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 2 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 near 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.

Computational Method
Computational methods to model and solve for multi-component micro-circulation has 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), multiparticle collision dynamics (MCP) [42? ? -45]. 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 a immersed boundary method, and employ the lattice Boltzmann method as the fluid solver.

Lattice Boltzmann Method
The evolution of the fluid is defined in terms of a set of N discrete distribution functions,[ f i ],(i = 0,...,N − 1), which obey the dimensionless Boltzmann equation in which ⃗ x and t are the spatial and time coordinates, respectively; [⃗ e i ],(i = 0,...,N − 1) is the set of N discrete velocities; ∆t is the time step; and τ is the relaxation time given by the unique non-null eigenvalue of the collision term in the BGK-approximation [46]. The kinematic viscosity of the flow is strictly related to τ as ν = c 2 s (τ − 1 2 )∆t being c s = 1 velocity space of the equilibrium distribution based on the Hermite polynomial expansion of this distribution [48].

Immersed Boundary Treatment
Deforming body models are commonly based on continuum approaches using strain energy functions to compute the membrane response [49][50][51]. However, a particle-based model governed by molecular dynamics has emerged due to its mathematical simplicity while providing consistent predictions [52,10,53,54]. In this work, a particlebased 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 [ ⃗ F i ](i = 0,..., 8), accounting for the immersed boundary, is included as an additional contribution on the right-hand side of Eq. (1): ⃗ F i is expanded in term of the reticular Mach number, ⃗ e i c s , resulting in: where ⃗ f ib is a body force term. Due to the presence of the forcing term, the mass density and the momentum density Within this parametrisation, the forced Navier-Stokes equations is recovered with a second order accuracy [55-4 59]. The external boundaries of the computational domain are treated with the known velocity bounce back conditions by Zou and He [60]. The IBM procedure, extensively proposed and validated by Coclite and colleagues [61,28,29,27,62], is here adopted and the moving-least squares reconstruction by Vanella et al. [63] is employed to exchange all LBM distribution functions between the Eulerian lattice and the Lagrangian chain. Finally, the body force term in Eq.(5), ⃗ f ib , is evaluated through the formulation by Favier et al [64].
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 where ⃗ r i, j = ⃗ r i − ⃗ 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 being k b the bending constant, k v the current local curvature in the v-th vertex, 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 (θ − θ 0 ), with θ 0 the angle in the stress free configuration. Given this, the forces on the nodes v le f t , v, and v right are obtained as where l right and l le f t are the length of the two adjacent left and right edges, respectively, and ⃗ 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 = 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 term of the reference pressure p re f and directed along the normal inward unity vector of the l-th element 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 = ∫ x l dy l − y l dx l . Within this formulation k a = 1 returns a perfectly incompressible membrane. Note that, ⃗ F a l is evenly distributed to the two vertices connecting the l-th element (v le f t and v right ) as ⃗ 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 ∆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 ⃗ n − 1 and is given by: Hydrodynamics Stresses Pressure and viscous stresses exerted by the l-th linear element are: whereτ l and p l are the viscous stress tensor and the pressure evaluated in the centroid of the element, respectively; ⃗ n l is the outward normal unit vector while l l is its length. The pressure and velocity derivatives in Eq.s (11) and (12) are computed using a probe in the normal positive direction of each element, being the probe length 1.2∆x [63,65]

FLUID-STRUCTURE INTERACTION
Particles dynamics is determined by dynamics IB technique described in [62], using the solution of the Newton equation for each Lagrangian vertex, accounting for both internal, Eq.s (6), (8), (9), and (10), and external stresses, Eq.s (11) and (12). Then, no-slip boundary conditions are imposed using a weak coupling approach [61]. The total force ⃗ F tot v (t) 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, The Newton equation of motion is integrated by using the Verlet algorithm. Specifically, a first tentative velocity is considered into the integration process,⃗ x v,0 (t), obtained interpolating the fluid velocity from the surrounding lattice nodes then, the velocity at the time level t + ∆t is computed as It should be noted that the present formulation is unconditionally stable for small deformation of the capsule membrane and for small velocity variations applied, as previously demonstrated by the authors [61,27,62].

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 x−axis or co-axial direction, and the vertical axis is denoted by y−axis or radial direction. The analysis is based on simulations either a single or three inline 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, 7 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∆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 [66]. The Reynolds number is fixed equal to 0.1 and is given by Re = u max d ν ; 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 ν the kinematic viscosity of the fluid (ν = 1.2 × 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 = ρνu max k s , where ρ is the density of the fluid in which capsules are immersed (ρ = 1000 kg/m 3 ) in and k s is the elastic constant used for the worm-like chain composing the membranes (see Eq (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 [25,7]. Note that, the typical stiffness for a red blood cell is k s = 1.5 × 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 τ = 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.

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.  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 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 [36,37,31].
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 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 which 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. 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 ≈ 8. We observe that the perimeter length change δ 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 are larger in the narrower channel due to higher shear rates, more so with 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 proportion 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 is 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.

13
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 intercapsule interaction, due to the relatively large capsule separation distance and the small capsule size which 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 capsule 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.

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 important mechanism which 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 were investigated. The 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 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 threedimensional 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 is 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 threedimensional modelling is certainly of interest and will be the object of future investigations.