Deformation of a Capsule in a Power-Law Shear Flow

An immersed boundary-lattice Boltzmann method is developed for fluid-structure interactions involving non-Newtonian fluids (e.g., power-law fluid). In this method, the flexible structure (e.g., capsule) dynamics and the fluid dynamics are coupled by using the immersed boundary method. The incompressible viscous power-law fluid motion is obtained by solving the lattice Boltzmann equation. The non-Newtonian rheology is achieved by using a shear rate-dependant relaxation time in the lattice Boltzmann method. The non-Newtonian flow solver is then validated by considering a power-law flow in a straight channel which is one of the benchmark problems to validate an in-house solver. The numerical results present a good agreement with the analytical solutions for various values of power-law index. Finally, we apply this method to study the deformation of a capsule in a power-law shear flow by varying the Reynolds number from 0.025 to 0.1, dimensionless shear rate from 0.004 to 0.1, and power-law index from 0.2 to 1.8. It is found that the deformation of the capsule increases with the power-law index for different Reynolds numbers and nondimensional shear rates. In addition, the Reynolds number does not have significant effect on the capsule deformation in the flow regime considered. Moreover, the power-law index effect is stronger for larger dimensionless shear rate compared to smaller values.


Introduction
Flow induced deformation of a capsule consisting of a membrane enclosing an internal medium such as a gel or a liquid is an important problem in fundamental research as well as bioengineering applications. For example, a capsule in shear flow is a fundamental process that is related to erythrocytes (or red blood cells), leukocytes (or white blood cells), and platelets in blood flow [1][2][3][4][5][6]. Deformation is essential for red blood cells to perform their physiological functions in the circulation of capillary blood vessels and thus affects the rheology of the blood [6][7][8]. The deformations of white blood cells and red blood cells can, respectively, affect the immune response and the oxygen load release [9,10]. The synthetic microcapsules with polymerized interfaces are designed for drug delivery, cosmetic production, and other technical usages [11,12]. Therefore, great effort has been made to study this problem (e.g., [1,4,6,8,10,[12][13][14]).
Both experimental and numerical methods have been conducted to observe capsule behaviors and the relevant underneath fluid-structure interaction physics. Schmid-Schönbein and Wells [15] and Goldsmith [16] observed that red blood cells tumble like rigid particles at low shear rates while they deform to a steady configuration and direction after which the membrane rotates around the internal liquid (tank-treading movement) at high shear rates. Later, Goldsmith and Marlow [17] and Keller and Skalak [18] found that the viscosity ratio between the liquids inside and outside the cell may also affect the type of behaviors. A higher viscosity inside would cause unsteady tumbling-rotating motion, while a smaller viscosity inside would lead to the tank-treading movement with a stationary shape. These phenomena were captured by Xu et al. [14]. More recently, Dupire et al. [19] reported rolling motion in addition to other behaviors. A hysteresis cycle and two transient dynamics driven by the shear rate (i.e., an intermittent regime during the "tanktreading-to-flipping" transition and a Frisbee-like "spinning" regime during the "rolling-to-tank-treading" transition) were highlighted.
There are several numerical methods that have been used to study capsule dynamics. Examples are the boundary element method (e.g., [20]), arbitrary Lagrangian-Euler method (e.g., [21][22][23]), immersed finite element method (e.g., [24]), and immersed boundary method (IBM) (e.g.,  [12][13][14][25][26][27][28][29][30][31][32][33][34]). Specifically, Zhou and Pozrikidis [20] studied the transient and large deformation of capsules with position-dependent membrane tension. Choi and Kim [21] simulated the motion of red blood cells freely suspended in shear flow to investigate the nature of pairwise interception of red blood cells using a fluid-particle interaction method based on the arbitrary Lagrangian-Eulerian method. Villone et al. [22,23] studied the effect of the non-Newtonian fluid on flexible particle deformation and migration in shear and channel flows by using the arbitrary Lagrangian-Eulerian method. The Navier-Stokes equations and cell-cell interaction were coupled in the framework of the immersed finite element method and mesh-free method by Y. Liu and W. K. Liu [24] to model complex blood flows with deformable red blood cells within micro and capillary vessels in three dimensions. The transient deformation of a liquid-filled elastic capsule in simple shear flow was studied by Sui et al. [1,4,35,36]. The fluid inertia on the dynamics of deformable particles has been studied by Krüger et al. [32] and Kaoui and Harting [34]. More recently, optical force based separation of particles/capsules was simulated by Chang et al. [37][38][39]. Still, as far as known to us, the existing numerical simulations seldom consider the non-Newtonian rheology effects on the capsule behaviors, while blood and most fluids involved in biomedical engineering are non-Newtonian fluids [6,8,40,41].
Following the work by Sui et al. [1] and Xu et al. [14], we develop an immersed boundary-lattice Boltzmann method (IB-LBM) to study the non-Newtonian effects on the deformation of a capsule in a shear flow. As a typical rheology, the power-law fluid is used. In this method, the capsule dynamics and the fluid dynamics are coupled by using the IBM, and the incompressible viscous power-law fluid motion is acquired by solving the lattice Boltzmann equation (LBE).
The rest of this paper is organized as follows. Section 2 briefly introduces the governing equations of the fluid and solid structures and describes the numerical approach. Section 3 presents the numerical results. Final conclusions are given in Section 4.

Physical Model and Mathematical Formulation.
In this work, a two-dimensional liquid capsule enclosed by an elastic membrane and immersed in an incompressible non-Newtonian fluid is considered, as illustrated in Figure 1 where is the arch length coordinate, n denotes the surface normal that points into the outer fluid, t denotes the tangent unit vector that points to the increasing arc length, and 0 is the velocity applied at both top and bottom walls to form a simple shear flow. The incompressible non-Newtonian fluid dynamics is achieved by using LBM [42,43]. Great effort has been made in using LBM to solve the complex flows (see several reviews [42][43][44] for the effort). Many publications have presented the details of LBM; thus we just provide a brief description in this paper and discuss the extension for non-Newtonian fluids. The details of LBM and its applications are referred to the references provided. Using the IB-LBM, the lattice Boltzmann equation (LBE) that governs the viscous flow dynamics and incorporates the traction jump across the interface due to the elastic membrane is written as [1,14,42,43,45,46] (x + e Δ , + Δ ) − (x, ) where (x, ) is the distribution function for particles with velocity e at position x and time , Δ is the size of the time step, eq (x, t) is the equilibrium distribution function, LB represents the dimensionless relaxation time, is the term representing the body force effect on the distribution function, are the weights, u = ( , V) is the velocity of the fluid, is the speed of sound defined by = Δ / √ 3Δ with Δ being grid spacing, f is the body force acting on the fluid, ΔF( , ) is the Lagrangian force density on the fluid by the elastic boundary, X is the position vector on the membrane, and (x − X(s, )) is Dirac's delta function.
In the present work, a two-dimensional nine-speed (D2Q9) model is used, as shown in Figure 2. In this model, the nine possible particle velocities are given by The values of e ensure that, within one time step, a particle moves to one of the eight neighboring nodes as shown in Figure 2 or stays at its current location. The weights, , are given by 0 = 4/9 and = 1/9 for where ] = / with being the dynamic viscosity of the ambient fluid and being the fluid density. When the particle density distributions are known, the fluid density, velocity, and pressure are then computed from Theoretically the LBM introduced above simulates the compressible viscous flow instead of incompressible viscous one, because the spatial density variation is not zero in LBM simulations. In the applications, the Mach number (Ma = 0 / ) should be low (e.g., Ma ≤ 0.3) so that the incompressible viscous flow can be correctly simulated. The deduction process from LBE to the incompressible viscous flow governing equations can be found in [47]. The dynamics viscosity is a constant for a Newtonian fluid, while it is dependent on the local shear rate for a non-Newtonian fluid. Without loss of generality, the power-law fluid is taken as a representation of non-Newtonian fluids in the present paper. The rheological equation of state for powerlaw fluids is defined by [48] =̇− 1 , where is the power-law consistency index, is the powerlaw fluid behavior index,̇is the shear rate, anḋis the minimum shear rate that is applied to avoid the numerical singularity caused by the zero shear rate. The power-law fluids of < 1, > 1 and = 1 are, respectively, the shear-thinning, shear-thickening, and Newtonian fluids. In (9), the Einstein summation convention is applied. In LBM implementation, can be either calculated macroscopically by using the finite difference method or locally in mesoscopic scale by using eq and f [49]. To achieve the non-Newtonian rheology, a shear rate-dependant relaxation time is used which can be obtained by applying the effective viscosity determined by (8) in (6).
Because of the deformation, the membrane develops a transverse shear tension and a bending moment . In addition, due to the stretching motion, a tension, , is induced. Consider the force balance of membranes; we acquire Please note that ΔF is the Lagrangian force on the fluid exerted by the elastic body boundary and is opposite to the fluid force on the boundary. To evaluate and for the thin membrane, we use Hooke's law which is a relatively simpler constitutive law for modeling small deformation of capsules. Hooke's law states that the tension and the bending moment are linearly related to the stretch and the curvature, respectively. It can be written in the form where is the bending coefficient, is the stretching coefficient, 0 is the initial arch length, is the curvature of membrane, and 0 is the curvature in the resting configuration. If is large so that the stretching deformation is small, Hooke's law works well. Following the work by Sui et al. [1], the capsule membrane is assumed to be infinitely thin so that the bending effect is neglected; that is, = 0. Actually, the effect of is similar to that of when is small compared to [1,35]. If is large, the capsule may undergo flipping motion [35].
The velocity of a point on the capsule is interpolated from the flow field, and the position of the capsule is updated explicitly; that is, where U( , ) is the velocity of the capsule.
In this work, we choose the flow shear rate (e.g., 0 /ℎ), density, and the radius of the capsule to nondimension the governing equations and obtain two dimensionless parameters: the Reynolds number Re and dimensionless shear rate , which are defined by where is the radius of the undeformed capsule. measures the ratio of shear force to the elastic force. For applications where inertia force is important, we can also use Re ⋅ to nondimension the elastic property, which measures the ratio of fluid inertial forces to stretching elastic forces. Please note that the two-dimensional model is used in this work, while red blood cell deformation is a three-dimensional problem; however, the results obtained in this research should show some features common with the three-dimensional simulations, as demonstrated in [1].

Numerical
Method. Similar to [1], the capsule is discretized by nodal points which are initially distributed with equal distances. The position of the th node at time level is denoted by X . To compute the stretching force at th node, a finite difference scheme is used; that is, where Δ is the Lagrangian grid spacing on the membrane and the tension and tangent vector, t = X/ , at the segment center, + 1/2, are both computed using a secondorder central difference scheme. The time integration of (14) is calculated according to In the IBM, a smooth approximation [50] of Dirac's delta function, ℎ , is used, In the present simulations, Δ = Δ = Δ (in lattice units) is used. Now, the computational algorithm can be summarized as follows. Given all values at time step , the values at time step + 1 can be undated by the following: (1) Calculate the Lagrangian force density ΔF +1 from X by using (11)-(12).
In the present work, the above-mentioned computational simulation algorithm is implemented in the Fortran 90 programming language.

Validation.
The IB-LBM in this work has been validated and verified in our previous studies (see, e.g., [14,46]) and has been used to study filament(s) flapping in viscous fluids [51][52][53], sperm swimming, and cell/particle flows [10,54]. In the present work, we focus on the validation of non-Newtonian flow by considering a power-law flow in a straight channel which is one of the benchmark problems to validate an in-house computational fluid dynamics solver. As in our previous work [41], we consider a two-dimensional steady laminar developing flow of power-law fluid with a uniform incoming velocity ∞ through a rectangular channel of height ℎ and length , as shown in Figure 3. The physically realistic initial and boundary conditions are given as   Figure 4 with the corresponding analytical solution for fully developed velocity profile [41,48] for power-law fluid flow in a channel which is given as From Figure 4, it is found that the present numerical results show a good agreement with the analytical solutions for various values of power-law index, giving us confidence in the reliability and accuracy of the present numerical solution procedure. It is noted from Figure 4 that the shear layer is thinned for < 1 and thickened for > 1 compared to the Newtonian fluid case ( = 1).

Numerical Results
We first consider the power-law index effect on the deformation of a cylindrical capsule in a shear flow. The Reynolds number is 0.05, which is in the range of normal physiological conditions. The dimensionless shear rate is 0.04. The computational domain ranges from 0 to 20 in bothaxis and -axis. The capsule is at the center of the domain, and its membrane is equally discretized into 80 Lagrangian nodes. The grid resolution is Δ = Δ = Δ = 0.1 .
The characteristic velocity is set as 0 = 0.05 so that the dimensionless relaxation time is 0.5 < LB < 3.0. Such setup is consistent with that used in [1]. To study the power-law index effect, is set in the range of 0.2 < < 1.8, covering the shearthinning, Newtonian, and shear-thickening fluids. In order to quantify the deformation of the capsule, the Taylor shape parameter is introduced [1], where and B are, respectively, the length and width of a cross-section of the cylindrical capsule. Figure 5 shows the deformation of the flexible capsule in a shear flow for Reynolds number of Re = 0.05, dimensionless shear rate of = 0.04, and power-law index of = 0.2 to 1.8. There are several interesting observations from Figure 5. First, the capsule deforms to a steady shape and then the membrane rotates around the liquid inside (tank-treading motion), which is further indicated by the streamlines in Figure 6. Second, the deformation increases with the powerlaw index. When the fluid is shear-thinning (i.e., < 1.0), the deformation is smaller compared to the Newtonian fluid case ( = 1.0), while the deformation is larger compared to the Newtonian fluid case for the shear-thickening fluid; that is, > 1.0. This can be explained by the power-law rheology. When < 1.0, the effective viscosity near the capsule is smaller compared to the Newtonian fluid, while the effective viscosity near the capsule is higher than that of Newtonian fluid for > 1.0. Based on the definition of in (16), the local is larger for > 1.0 and smaller for < 1.0 compared to that of Newtonian fluid. As presented by Sui et al. [1], a larger corresponds to a larger , that is, larger deformation of the capsule. Third, the Taylor shape parameter , which is used to quantify the deformation, increases with the powerlaw index. Finally, it is noted that y is approximately linear function of , as shown in Figure 5(b).
In order to study the Reynolds number effect on the deformation of the capsule, we simulate two additional Reynolds numbers, Re = 0.1 and 0.025. Figure 7 shows the deformation of the flexible capsule in a shear flow for Reynolds number of Re = 0.1 and 0.025, dimensionless shear rate of = 0.04, and power-law index of 0.2 ≤ ≤ 1.8. It is found that the deformation ( ) for Re = 0.1 is larger compared to the cases of Re = 0.05 and 0.025. However, the difference is quite small, implying that, in the low Reynolds number regime, for example, Re ≤ 0.1 in this work, the deformation of the capsule is not significantly affected by the Reynolds numbers used, as the inertial force is ignorable here, and the shear forces and capsule elastic forces are dominant. Therefore, the dimensionless shear rate ( ) should significantly affect the deformation of the capsule, which will be further verified by the simulations shown below by varying .
Finally, we study the shear rate effect on the deformation of the capsule by using  observations are obtained. First, the capsule deformation is sensitive to the dimensionless shear rate. This can be explained by the definition of in (16): measures the ratio of shear (viscous) forces to the stretching elastic forces, which is the dominant physical process here. A change of this ratio would cause significant difference in the capsule deformation. Second, the power-law index effect is stronger for larger , as indicated by the slopes of the functions shown in Figure 8(c). This can be explained by the fact that the physical for Re = 0.05 is shown in (c) for comparison. process changes from a shear force dominant to an elasticforce dominant process when varies from 0.1 to 0.004. For low , for example, 0.004, the elastic forces are dominant, and thus the shear force change caused by the change of the non-Newtonian rheology is smaller compared to that for larger , for example, 0.1. Finally, the deformed capsule is obviously biased from elliptical cylinder for large and ; for example, = 0.1 and ≥ 1.2. This is caused by the shear-induced torque on the deformed capsule and the decrease of effective bending resistance caused by the shear-induced elongation.
To further discuss the non-Newtonian effect, = − , which measures the local shear stress, is introduced [40]. compared to that for = 0.6. This is a further explanation of larger deformation for larger .

Conclusion
A numerical approach combining the immersed boundary method and the lattice Boltzmann method has been developed for fluid-structure interactions involving non-Newtonian fluids. Without loss of generality, the power-law fluid is taken as a representation of non-Newtonian fluids to present the method. This method couples the flexible structure (e.g., capsule) dynamics and the fluid dynamics by using the immersed boundary method and calculates the incompressible viscous power-law fluid motion by solving the lattice Boltzmann equation. In order to achieve the non-Newtonian rheology, a shear rate-dependant relaxation time is employed. The non-Newtonian flow solver has been validated by conducting a power-law flow in a straight channel. The power-law index has been varied from 0.6 to 1.4. The present 8 Computational and Mathematical Methods in Medicine numerical results show a good agreement with the analytical solutions for various values of power-law index, giving us confidence in the reliability and accuracy of the present numerical solution procedure.
To study the non-Newtonian effects on the deformation of a capsule in a power-law shear flow, we have performed simulations by varying the Reynolds number from 0.025 to 0.1, dimensionless shear rate from 0.004 to 0.1, and power-law index from 0.2 to 1.8. It is found that the capsule deformation increases with the power-law index for different Reynolds numbers and nondimensional shear rates. In addition, the Reynolds number does not have significant effect on the capsule deformation in the flow regime considered. Finally, the power-law index effect is stronger for larger dimensionless shear rate compared to smaller values.