Numerical simulation of red blood cell suspensions behind a moving interface in a capillary
aa r X i v : . [ phy s i c s . f l u - dyn ] A p r Numerical simulation of red blood cell suspensionsbehind a moving interface in a capillary
Shihai Zhao and Tsorng-Whay PanDepartment of Mathematics, University of Houston, Houston, Texas 77204-3008, USA
Abstract : Computational modeling and simulation are presented on the motion of red bloodcells behind a moving interface in a capillary. The methodology is based on an immersedboundary method and the skeleton structure of the red blood cell (RBC) membrane ismodeled as a spring network. The computational domain is moving with either a designatedRBC or an interface in an infinitely long two-dimensional channel with an undisturbed flowfield in front of the domain. The tanking-treading and the inclination angle of a cell in asimple shear flow are briefly discussed for the validation purpose. We then present the resultsof the motion of red blood cells behind a moving interface in a capillary, which show thatthe RBCs with higher velocity than the interface speed form a concentrated slug behind theinterface.
Key words: red blood cells, moving domain, immersed boundary method.
The rheological property of the red blood cells (RBCs) is a key factor of the blood flow char-acteristics at the microchannel level, especially the particulate nature of the blood becomessignificant when studying blood drop through a glass capillary within miniature blood diag-nostic kit. The penetration of the blood suspension in a perfectly wettable capillary has beenanalyzed in [1, 2]. The failure of such penetration is attributed to three RBCs segregationmechanisms: (i) corner deflection at the entrance, (ii) the intermediate deformation-inducedradial migration and (iii) shear-induced diffusion within a packed slug at the meniscus. Thekey mechanism responsible for penetration failure is the deformation-induced radial migra-tion, which endows the blood cells with a higher velocity than the meniscus to form theconcentrated slug behind the meniscus (see Figure 1). The results in [1, 2] shed light onmaking the smallest microfluidic kit and loading microneedle that require the least amountof blood sample.Nowadays in silico mathematical modeling and numerical study of RBC rheology haveattracted growing interest (see, e.g., [3, 4]). The immersed boundary method developed byPeskin, e.g, [5, 7, 6], has been one of the popular methodologies for numerically studying theRBC rheology due its distinguish features in dealing with the problem of fluid flow inter-acting with a flexible fluid/structure interface. For example, in [8]-[17], immersed boundarymethods have been combined with different RBC membrane models to simulate the motion Figure 1: Schematics of the BRCs moving behind a meniscus.of RBCs and vesicles in fluid flow. In [15, 16, 17], we have successfully combined an im-mersed boundary method with a spring model developed in [18] to simulate the motion ofRBCs in shear flows and Poiseuille flows. In this paper we have generalized the aforemen-tioned methodology to simulate the RBCs aggregation behind a moving interface consideredin [1, 2] by having the computational domain moving with an interface in an infinitely longtwo-dimensional channel with an undisturbed flow field in front of the domain since thetypical periodic boundary condition in the direction of the channel wall is not well suitedanymore. To mimic the motion of the RBCs behind a meniscus in a capillary, we haveconsidered a flat interface moving with a given constant speed in this paper. The simulatingresults of the motion of red blood cells behind a moving interface show that the RBCs withhigher velocity than the interface speed form the concentrated slug behind the interface,which resembles the motion of the RBCs observed in [1, 2]. The structure of this paper is asfollows: We discuss the elastic spring model and numerical methods in Section 2. In Section3, the tanking-treading and the inclination angle of a cell in a simple shear flow are brieflydiscussed for the validation purpose. We then present the results of the motion of red bloodcells behind a moving interface in a capillary. The conclusions are summarized in Section 4.
Let Ω be a bounded rectangular domain filled with blood plasma which is incompressible,Newtonian, and contains RBCs with the viscosity of the cytoplasm same as that of the bloodplasma (see Figure 2). For some
T >
0, the governing equations for the fluid-cell system are ρ (cid:20) ∂ u ∂t + u · ∇ u (cid:21) = − ∇ p + µ ∆ u + f in Ω , t ∈ (0 , T ) , (2.1) ∇ · u = 0 in Ω , t ∈ (0 , T ) (2.2)where u and p are the fluid velocity and pressure, respectively, ρ is the fluid density, and µ is the fluid viscosity, which is assumed to be constant for the entire computational domain. Ω x x Γ Γ Γ Γ Figure 2: An example of the computational domain with a cell.In (2.1), f is a body force which accounts for the force acting on the fluid/cell interface.Equations (2.1) and (2.2) are completed by the following boundary and initial conditions: u = g on Γ d , (2.3) µ ∂ u ∂ n − n p = on Γ n , (2.4) u (0) = u (2.5)where the domain Ω is taken from an infinitely long channel with its boundary denoted byΓ = ∪ i =1 Γ i . In the simulations, we have considered two types of boundary conditions: (i)Γ n = ∅ and Γ d = Γ, (ii) Γ n = Γ and Γ d = Γ ∪ Γ ∪ Γ with g having the profile of eitherPoiseuille flow or simple shear flow on Γ d A two-dimensional elastic spring model used in [18] is considered in this paper to describe thedeformable behavior of the RBCs. Based on this model, the RBC membrane can be viewedas membrane particles connecting with the neighboring membrane particles by springs, asshown in Figure 3. Elastic energy stores in the spring due to the change of the length l ofthe spring with respected to its reference length l and the change in angle θ between twoneighboring springs. The total elastic energy of the RBC membrane, E = E l + E b , is the sumof the total elastic energy for stretch/compression and the total energy for bending which,in particular, are E l = k l N X i =1 ( l i − l l ) (2.6)and E b = k b N X i =1 tan ( θ i / . (2.7)In equations (2.6) and (2.7), N is the total number of the spring elements, and k l and k b are spring constants for changes in length and bending angle, respectively. θ l Figure 3: The elastic spring model of the RBC membrane
Remark . In the process of creating the initial shape of RBCs described in [18], the RBCis assumed to be a circle of radius R = 2 . µm initially. The circle is discretized into N = 76membrane particles so that 76 springs are formed by connecting the neighboring particles.The shape change is stimulated by reducing the total area of the circle through a penaltyfunction Γ s = k s s − s e s e ) (2.8)where s and s e are the time dependent area of the RBC and the equilibrium area of theRBC, respectively, and the total energy is modified as E + Γ s . Based on the principle ofvirtual work the force acting on the i th membrane particle now is F i = − ∂ ( E + Γ s ) ∂ r i (2.9)where r i is the position of the i th membrane particle. When the area is reduced, each RBCmembrane particle moves on the basis of the following equation of motion: m ¨ r i + γ ˙ r i = F i (2.10)Here, ˙() denotes the time derivative; m and γ represent the membrane particle mass and themembrane viscosity of the RBC. The position r i of the i th membrane particle is solved bydiscretizing (2.10) via a second order finite difference method. The total energy stored inthe membrane decreases as the time elapses. The final shape of the RBC is obtained as thetotal elastic energy is minimized (please see [19]). The area of the final shape has less than0 . s e and the length of the perimeter of thefinal shape has less than 0 . s ∗ = s e /πR . Remark . When simulating the case involving a moving interface, we have applied arepulsive force to prevent the overlapping between cell and wall. The repulsive force isobtained from the following Morse potential (e.g., see [20]) φ ( d ) = k r (1 − e − ( d − d ) ) where the parameter d is the shortest distance between the membrane particle and the walland d is the range of the repulsive force (when the distance d is greater than d , there is norepulsive force). The parameter k r is a constant for the strength of the potential. The immersed boundary method developed by Peskin, e.g, [5, 7, 6], is employed in this studybecause of its distinguish features in dealing with the problem of fluid flow interacting with aflexible fluid/structure interface. Over the years, it has demonstrated its capability in studyof computational fluid dynamics including blood flow. Based on the method, the boundaryof the deformable structure is discretized spatially into a set of boundary nodes. The forcelocated at the immersed boundary node r i = ( r i, , r i, ) affects the nearby fluid mesh nodes x = ( x , x ) through a 2D discrete δ -function D h ( x − r i ): f ( x ) = X F i D h ( x − r i ) f or | x − r i | ≤ h, (2.11)where h is the uniform finite element mesh size and D h ( x − r i ) = δ h ( x − r i, ) δ h ( x − r i, ) (2.12)with the 1D discrete δ -functions being δ h ( z ) = h (cid:16) − | z | /h + p | z | /h − | z | /h ) (cid:17) , | z | ≤ h, h (cid:16) − | z | /h − p − | z | /h − | z | /h ) (cid:17) , h ≤ | z | ≤ h, , otherwise. (2.13)The velocity of the immersed boundary node r i is also affected by the surroundingfluid and therefore is enforced by summing the velocities at the nearby fluid mesh nodes x weighted by the same discrete δ -function: U ( r i ) = X h u ( x ) D h ( x − r i ) f or | x − x i | ≤ h. (2.14)After each time step, the position of the immersed boundary node is updated by r n +1 i = r ni + ∆ t U ( r ni ) . (2.15) Concerning the finite element based space approximation of { u , p } in problem (2.1)-(2.5),we use the P - iso - P and P finite element approximation (e.g., see [22] (Chapter 5)). For arectangular computational domain Ω ⊂ R , let T h be a finite element triangulation of Ω forvelocity and T h a twice coarser triangulation for pressure where h is a space discretizationstep. We introduce the finite dimensional spaces: W h = { v h | v h ∈ C (Ω) , v h | T ∈ P × P , ∀ T ∈ T h } W h = { v h | v h ∈ W h , v h = on Γ d } ,L h = { q h | q h ∈ C (Ω) , q h | T ∈ P , ∀ T ∈ T h } ,L h, = { q h | q h ∈ L h , Z Ω q h d x = 0 } where P is the space of polynomials in two variables of degree ≤
1. We apply the
Lie’sscheme [21, 22] with the above finite elements to equations (2.1)-(2.5) with the backwardEuler method in time for some subproblems and obtain the following sequence of fractionalstep subproblems (some of the subscripts h have been dropped): u = u is given; for n ≥ u n being known, we compute the approximate solution viathe following fractional steps:1. Update the position of the membrane by (2.14) and (2.15) and then compute the force f n based on the fluid/cell interface by (2.9) and (2.11).2. Solve Z Ω ∂ u ( t ) ∂t · v d x + Z Ω ( u n · ∇ ) u ( t ) · v d x = 0 , on ( t n , t n +1 ) , ∀ v ∈ W − h , u ( t n ) = u n , u ( t ) ∈ W h , u ( t ) = g ,h on Γ − × ( t n , t n +1 ) , (2.16)and set u n +2 / = u ( t n +1 ).3. Finally solve ρ Z Ω u n +1 − u n +2 / △ t · v d x + µ Z Ω ∇ u n +1 : ∇ v d x − Z Ω p n +1 ( ∇ · v ) d x = Z Ω f n · v d x , ∀ v ∈ W h , Z Ω q ∇ · u n +1 d x = 0 , ∀ q ∈ L h , u n +1 ∈ W h , u ( t ) = g ,h on Γ d ; p n +1 ∈ L h ( p n +1 ∈ L h, if Γ n = ∅ ) . (2.17)In eq. (2.16), we have Γ − = { x | x ∈ Γ, g ,h ( x ) · n ( x ) < } and W − h = { v h | v h ∈ W h , v h = on Γ − } . The quasi-Stokes problem (2.17) is solved by a preconditioned conjugate gradientmethod (see, e.g., [22]). The subproblem (2.16) is an advection type subproblem. It is solvedby a wave-like equation method, which is described in detail in [23] and [24]. Remark . In simulations, the computational domain Ω moves to the right with eitherthe mass center of a RBC or the interface (see, e.g., [25, 26] and references therein foradjusting the computational domain according to the position of the particle). Due to theuse of structured and uniform mesh in our simulations, it is relatively easy to have thecomputational domain moving with a designated cell. Generally when the mass center of aRBC moves to the right in an infinitely long channel, we add one vertical grid line to theright end of the computational domain if the cell mass center crosses one vertical grid lineafter we predict its new position and at the same time we drop one vertical grid line at theleft end of the computational domain. In the mean time at these new grid points addedat the right end, we assign the values of velocity field according to either Poiseuille flow orsimple shear flow depending on the test cases. When following an interface moving to theright with a constant speed, we have applied the same strategy.
Swelling ratio I n c li na t i on ang l e case Icase IIcase III([15])([13]) Swelling ratio T an k − t r ead i ng v e l o c i t y case Icase IIcase III([15])([13]) Figure 4: (Color online). Steady inclination angle versus the cell swelling ratio (left) andmembrane tank-treading velocity (scaled by γR /
2) versus the cell swelling ratio (right) incomparison to Shi. et al.[15] and Kaoui. et al.[13] in different cases. Case I: 112 µ m × µ mdomain with Dirichlet boundary conditions, Case II: 80 µ m × µ m domain with Dirichletboundary conditions, Case III: 80 µ m × µ m domain with Neumann inflow condition andDirichlet outflow condition.We have first validated the computational methodology with two types of boundaryconditions discussed in Section 2 by comparing the inclination angle and the tank-treadingfrequency of a single RBC in shear flow. Here are the parameters used in the simulations: Thevalues of parameters for modeling cells are same with [15, 16, 17] as follows: The bendingconstant is k b = 5 × − N · m, the spring constant is k l = 5 × − N · m, the penaltycoefficient is k s = 10 − N · m, the repulsive force coefficient is k r = 10 − N · m, and the rangeof the repulsive force is d = 2 h where h is the mesh size for the flow velocity field. The cells
10 20 30 40 50 60 70 800510 x ( µ m) y ( µ m ) t=0.01ms
20 30 40 50 60 70 80 900510 x ( µ m) y ( µ m ) t=0.5ms70 80 90 100 110 120 130 1400510 x ( µ m) y ( µ m ) t=2.5ms130 140 150 160 170 180 190 2000510 x ( µ m) y ( µ m ) t=4.5ms170 180 190 200 210 220 230 2400510 x ( µ m) y ( µ m ) t=6.25ms Figure 5: (Color online). The positions of 12 cells in a capillary behind a moving interfaceat t = 0.01, 0.5, 2.5, 4.5 and 6.25 ms and the velocity field with 12 cells at t =6.25 ms (fromtop to bottom).
20 40 60 80 100 120 140 16001020 x ( µ m) y ( µ m ) t=0.01ms50 100 150 20001020 x ( µ m) y ( µ m ) t=1.8ms150 200 250 30001020 x ( µ m) y ( µ m ) t=5.37ms280 300 320 340 360 380 400 42001020 x ( µ m) y ( µ m ) t=10ms Figure 6: (Color online). The position of 68 cells in a capillary behind a moving interfaceat t = 0.01, 1.8, 5.37 and 10 ms and the velocity field with 68 cells at t = 5.37 and 10 ms(from top to bottom).0are suspended in blood plasma which has a density ρ = 1 . / cm and a dynamical viscosity µ = 0 . / (cm · s). The viscosity ratio which describes the viscosity contrast of the innerand outer fluid of the RBC membrane is fixed at 1.0. The dimensions of the computationaldomain are 112 µ m × µ m and 80 µ m × µ m. Then the degree of confinement (2 R /H ) is 0.8where H is the height of the channel. The grid resolution for the computational domain is80 grid points per 10 µ m. The time step ∆ t is 1 × − ms. The initial position of the masscenter of the cells are (56, 3.5) and (40, 3.5) for the longer domain and the shorter domain,respectively. To have a shear flow, a Couette flow driven by two walls at the top and bottomwhich have the same speed U/ U is given by U = γ ∗ H with a given shear rate γ . The shearrate used in the simulation is γ = 275/s. The steady inclination angles of the tank-treadingfor four values of s ∗ =0.6, 0.7, 0.8 and 0.9 are presented in Figure 4, which show the verygood agreement with the lattice-Boltzmann simulation results in [13] and those previouslyobtained with periodic boundary conditions in [15]. The membrane tank-treading velocity(scaled by γR /
2) is also in good agreement with the results in [13, 15]. The results showthat there is no significant difference when having the Dirichlet boundary conditions on Γwith the length L = 112 and 80 µ m or the conditions (2.3) and (2.4) on the boundary of theshorter domain. For the cases involving a moving interface in a capillary, we have considered the one movingto the right with constant speed U to mimic the motion of the RBCs behind a meniscus ina capillary. Then the associated boundary condition in (2.3) on Γ d is g = on Γ ∪ Γ and g = ( U, t on Γ and the boundary condition (2.4) is satisfied on Γ . We have keptall the related parameters the same except the following. We have first considered the caseof 12 cells of swelling ratio s ∗ =0.481 in a capillary of the height 10 µ m. The computationaldomain Ω is 80 µ m × µ m. The interface speed is U = 8 / s ∗ =0.481 in a capillaryof the height 20 µ m, we have considered the computational domain Ω = 160 µ m × µ m. Theinterface speed is U = 8 / In summary, we have developed computational modeling and methodologies for simulatingthe motion of many RBCs in a capillary behind a moving interface in this paper. Themethodology is based on an immersed boundary method and the skeleton structure of thered blood cell (RBC) membrane is modeled as a spring network. The computational domainis moving with either a designated RBC or an interface in an infinitely long two-dimensionalchannel with an undisturbed flow field in front of the domain. The tanking-treading andthe inclination angle of a cell in a simple shear flow are briefly discussed for the validationpurpose. The results of the motion of red blood cells behind a moving interface in a capillaryshow that the RBCs with higher velocity than the interface speed form a concentrated slugbehind the interface, which is consistent with the results in [1, 2]. The lateral migration isalso a key factor for the formation of a slug behind the moving interface.
Acknowledgments
The authors acknowledge the support of NSF (grant DMS-0914788). We acknowledge thehelpful comments of James Feng, Ming-Chih Lai and Sheldon X. Wang.
References [1]
H.-C. Chang, R. Zhou , Capillary penetration failure of blood suspensions , J. Col-loid Interface Sci. 287 (2005), pp. 647–656.[2]
R. Zhou, J. Gordon, A.F. Palmer, H.-C. Chang , Role of erythrocyte deforma-bility during capillary wetting , Biotechnology and Bioengineering 93 (2006), pp. 201-211.[3]
V. Cristini, G.S. Kassab , Computer odeling of red blood cell rheology in the mi-crocirculation: a brief overview. , Ann. Biomed. Eng. 33 (2005), pp. 1724–1727.[4]
C. Pozrikidis , Modeling and simulation of capsules and biological cells . Chapman& Hall/CRC: Boca Raton, 2003.[5]
C.S. Peskin , Numerical analysis of blood flow in the heart , J. Comput. Phys., 25(1977), pp. 220–252.[6]
C.S. Peskin , The immersed boundary method , Acta Numer., 11 (2002), pp. 479–517.[7]
C.S. Peskin, D.M. McQueen , Modeling prosthetic heart valves for numerical anal-ysis of blood flow in the heart , J. Comput. Phys., 37 (1980), pp. 113-32.2 [8]
C. Eggleton, A. Popel , Large deformation of red blood cell ghosts in a simpleshear flow , Phys. Fluids, 10 (1998), pp. 1834–1845.[9]
P. Bagchi, P. Johnson, A. Popel , Computational Fluid Dynamic Simulation ofAggregation of Deformable Cells in a Shear Flow , J. Biomech. Eng., 127 (2005), pp.1070–1080.[10]
P. Bagchi , Mesoscale simulation of blood flow in small vessels , Biophys. J., 92(2007), pp. 1858–1877.[11]
J. Zhang, J. Johnson, A.S. Popel , Effects of erythrocyte deformability and ag-gregation on the cell free layer and apparent viscosity of microscopic blood flows ,Microvasc. Res., 77 (2009), pp. 265–272.[12]
L.M. Crowl, A.L. Fogelson , Computational model of whole blood exhibiting lat-eral platelet motion induced by red blood cells , Int. J. Numer. Meth. Biomed. Engng.,26 (2010), pp. 471–487.[13]
B. Kaoui, J. Harting, C. Misbah , Two-dimensional vesicle dynamics under shearflow: Effect of confinement , Phys. Rev. E, 83 (2011), 066319.[14]
Y. Kim, M.-C. Lai , Numerical study of viscosity and inertial effects on tank-treadingand tumbling motions of vesicles under shear flow , Phys. Rev. E, 86 (2012), 066321.[15]
L. Shi, T.-W. Pan, R. Glowinski , Deformation of a single blood cell in boundedPoiseuille flows , Phys. Rev. E, 85 (2012), 016307.[16]
L. Shi, T.-W. Pan, R. Glowinski , Lateral migration and equilibrium shape andposition of a single red blood cell in bounded Poiseuille flows , Phys. Rev. E, 86 (2012),056308.[17]
L. Shi, T.-W. Pan, R. Glowinski , Numerical simulation of lateral migration ofred blood cells in Poiseuille flows , Int. J. Numer. Methods Fluids, 68 (2012), pp.1393–1408.[18]
K. Tsubota, S. Wada, T. Yamaguchi , Simulation study on effects of hematocriton blood flow properties using particle method , J. Biomech. Sci. Eng., 1 (2006), pp.159–170.[19]
T. Wang, T.-W. Pan, Z. Xing, R. Glowinski , Numerical simulation of rheologyof red blood cell rouleaux in microchannels , Phys. Rev. E, 79 (2009), 041916.[20]
A. Alexeev, R. Verberg, A.C. Balazs , Modeling the interactions between de-formable capsules rolling on a compliant surface , Soft Matter 2 (2006), pp. 499–509.[21]
A.J. Chorin, T.J.R. Hughes, M.F. McCracken, J.E. Marsden , Product for-mulas and numerical algorithms , Comm. Pure Appl. Math., 31 (1978), pp. 205–256.3[22]
R. Glowinski , Finite element methods for incompressible viscous flow , in Handbookof Numerical Analysis, Vol. IX, Ciarlet PG and Lions JL (Eds.). North-Holland,Amsterdam (2003), pp. 7-1176.[23]
E.J. Dean, R. Glowinski , A wave equation approach to the numerical solutionof the Navier-Stokes equations for incompressible viscous flow , C.R. Acad. Sc. Paris,S´erie 1, 325 (1997), pp. 783–791.[24]
E.J. Dean, R. Glowinski, T.-W. Pan , A wave equation approach to the numericalsimulation of incompressible viscous fluid flow modeled by the NavierStokes equations ,in Mathematical and Numerical Aspects of Wave Propagation, De Santo JA (Ed.).SIAM: Philadelphia(1998), pp. 65–74.[25]
H.H. Hu, D.D. Joseph, M.J. Crochet , Direct simulation of fluid particle mo-tions , Theoret. Comput. Fluid Dynamics 3 (1992), pp. 285–306.[26]