Numerical study of viscosity and inertial effects on tank-treading and tumbling
motions of vesicles under shear flow
Yongsam Kim*
Department of Mathematics, Chung-Ang University, Dongjakgu Heukseokdong, Seoul 156-756, Korea Ming-Chih Lai†
Center of Mathematical Modeling and Scientific Computing & Department of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsinchu 300, Taiwan
(Received 22 December 2011; revised manuscript received 27 September 2012; published 26 December 2012) An inextensible vesicle under shear flow experiences a tank-treading motion on its membrane if the viscosity contrast between the interior and exterior fluids is small. Above a critical threshold of viscosity contrast, the vesicle undergoes a tumbling bifurcation. In this paper, we extend our previous work [Kim and Lai,J. Comput. Phys. 229, 4840 (2010)] to the case of different viscosity and investigate the transition between the tank-treading and tumbling motions in detail. The present numerical results are in a good agreement with other numerical and theoretical studies qualitatively. In addition, we study the inertial effect on this transition and find that the inertial effect might inhibit the tumbling motion in favor of the tank-treading motion, which is observed recently in the literature. The critical viscosity contrast for the transition to the tumbling motion usually increases as the reduced area increases in the Stokes regime. However, we surprisingly observe that the critical viscosity contrast decreases as the reduced area increases to some point in the flow of slightly higher Reynolds number. Our numerical result also shows that the inertial effect has stronger inhibition to tumbling motion when the reduced area is small. DOI:10.1103/PhysRevE.86.066321 PACS number(s): 47.63.−b, 83.80.Lz
I. INTRODUCTION
We introduce an immersed boundary (IB) method [1–4] to simulate the dynamics of inextensible vesicles in an incompressible viscous fluid with different viscosity contrasts. The dynamics of an individual vesicle in a flow has been widely investigated by mathematical studies [5–8] and numerical simulations [9–12] as well as experiments [13–16] (and the ref-erences therein). This interest in vesicle dynamics is motivated by the strong resemblance to many biological systems [9] such as red blood cells [17,18] and drug-carrying capsules [19].
It has been known for a while that a vesicle in a shear flow has two types of motion: tank-treading and tumbling. This transition is of fundamental importance since it can alter rheological properties of a vesicle solution by reducing dissipation [6,12]. The selection between these two types of vesicle motion depends on the reduced volume (a flattened vesicle would tumble more easily than a swollen one) and the viscosity contrast (the more viscous the internal fluid in comparison to the external one is, the easier is the tumbling) [5–8,12,15,16]. Several attempts to understand the transition between tank-treading and tumbling motions have been made in the literature. The most prominent and pioneering result is the work of Keller and Skallak [5] who have derived an evo-lution equation of the vesicle inclination angle in a shear flow. The dynamics of moving vesicles is determined by their membrane elasticity (bending resistance), inextensibility, and hydrodynamical force. In Ref. [4], we have applied an IB method to simulate the dynamics of inextensible vesicles. We have used a penalty idea [20] to enforce the inextensibility con-straint of the vesicle, verified the numerical method by a series
*http://cau.ac.kr/˜kimy; kimy@cau.ac.kr
†mclai@math.nctu.edu.tw
of convergence study, and shown that the vesicle in a simple shear flow undergoes a tank-treading tangential motion at its equilibrium state. The simulation results in Ref. [4] were com-parable to other theoretical and numerical works in literature. Since the method used in Ref. [4] has assumed that the fluid viscosities inside and outside the vesicle are the same, i.e., the viscosity contrast is always 1, we have observed only the tank-treading motion of the vesicle. In this paper, we extend the method used in Ref. [4] to handle the case in which the fluid has a viscosity contrast other than 1. Since the viscosity is discontinuous and has two different constant values across the boundary of the vesicle, we adopt the approach of Tryggvason
et al. [21] and introduce an indicator function to smooth out
the discontinuity of the viscosity. In this way we can regard the fluid inside and outside of the vesicle as the same one with a nonuniform viscosity. The details shall be presented in Sec.III. In order to show that the present IB method is a robust and efficient numerical tool to handle inextensible vesicles with a nonunity viscosity contrast, we first perform a convergence study, which confirms the consistent first-order convergence for the velocity field in L2 norm. Then we simulate several two-dimensional problems with a single vesicle under shear flow. From those runs, we show that the simulation results obtained by the IB method are comparable to those in the previous literature [5–8,12,15,16]. We also investigate the tank-treading/tumbling transitions in a flow with Reynolds number of order 1, and show that the inertia effect may inhibit the tumbling motion in favor of the tank-treading motion. This inhibition is stronger especially when the reduced area is small.
II. MODEL EQUATIONS
Consider a two-dimensional viscous incompressible fluid in a domain which contains an immersed, inextensible, mass-less vesicle. Our mathematical formulation and computational
In the pIB method [4], we use a dual representation of the immersed vesicle boundary. One of its representatives, denoted by X(s,t), interacts with the fluid and moves at the local fluid velocity, just as the traditional immersed boundary formulation. The other representative, denoted by Y(s,t), is elastic and linked to X(s,t) by a stiff spring system. The boundary points of Y(s,t) are not coupled directly to the fluid, and they move as if in a vacuum. Then the equations of motion of the pIB formulation are as follows:
ρ ∂u ∂t + u · ∇u = −∇p + ∇ · [μ(x,t)(∇u + ∇uT )]+ f, (1) ∇ · u = 0, (2) f(x,t)= F(s,t)δ(x− X(s,t))ds, (3) ∂X ∂t (s,t)= U(s,t) = u(x,t)δ(x− X(s,t))dx, (4) F= K(Y(s,t) − X(s,t)) + R(V(s,t) − U(s,t)), (5) K(Y(s,t)− X(s,t)) + R(V(s,t) − U(s,t)) = ∂ ∂s(σ (s,t)τ(s,t)) − cb ∂4Y(s,t) ∂s4 , (6) ∂V ∂s · ∂Y ∂s = 0. (7)
Equations(1)and(2)are the familiar Navier-Stokes equa-tions for an incompressible viscous fluid with a nonuniform viscosity. The fluid density ρ is a constant, but the viscosity μ(x,t) is a piecewise constant function which has two different values inside and outside the vesicle. The unknown variables in the fluid equations are the fluid velocity, u(x,t), and the fluid pressure, p(x,t), where x= (x,y) are fixed Cartesian coordinates, and t is the time.
Equation (3) and (4) both involve the two-dimensional Dirac delta function δ(x)= δ(x)δ(y), which expresses the local character of the interaction between the fluid and the immersed boundary. Equation(3)simply expresses the relation between the two corresponding force densities f(x,t)dx and F(s,t)ds, and Eq.(4)is the equation of motion of the vesicle boundary which indicates that the boundary moves at the local fluid velocity.
Equations(5)–(7)are the vesicle boundary equations which are written in Lagrangian form. Equation(5)defines the force density F which is transmitted by the boundary X(s,t) to the fluid. It includes only the force density generated by the stiff springs that link the two representations of the immersed boundary. Here, V= ∂Y/∂t is the velocity of the boundary Y(s,t), and K and R are the penalty parameters. The larger the penalty parameters are, the greater is the penalty energy paid to separate the positions and velocities of those two boundary representations.
Equation(6)is the equation for the position and velocity of the boundary Y(s,t) which describes that the penalty force F is instantaneously balanced by the elastic force of the boundary
Here the function τ(s,t) is the unit tangent vector to the boundary Y(s,t); i.e.,τ(s,t) =∂Y∂s/|∂Y∂s|. Whereas the bending coefficient cb is a constant, the elastic tension σ (s,t) is an unknown function which plays the role of Lagrange multiplier to enforce the inextensibility constraint of the vesicle boundary in Eq.(7).
Consider the case in which the penalty parameters K and R go to infinity. Then the boundaries X(s,t) and Y(s,t) and the velocities U(s,t) and V(s,t) will tend to coincide, and F in Eq.(5)approaches the elastic force density Feof the form
Fe(s,t)= ∂
∂s(σ (s,t)τ(s,t)) − cb
∂4X(s,t)
∂s4 , (8)
which is the standard formula for a vesicle under tension (first term) and bending resistance (second term). If we use the force density Fe(s,t) instead of Eqs.(5)and(6), the system of equations becomes the traditional IB formulation. Since the surface tension σ (s,t) is an unknown function, the traditional IB method should be solved in an implicit way, which requires an iterative method to solve the whole system. By setting the dynamics for the boundary Y apart from the whole system (especially the fluid equations), the pIB method can find the unknown function σ (s,t) more easily. The dynamics of the boundary Y is related to the whole system only through the penalty force F.
In practice, we choose K and R sufficiently large to keep the positions and velocities of the two representations of the boundary close to each other. Notice that the present pIB idea shares the same spirit with the “virtual boundary method” [22,23] for simulating a flow past a rigid boundary. Both methods introduce the penalty parameters K and R to penalize the difference in the desired and numerical velocities and positions of the immersed boundaries. The difference is that the virtual boundary method is developed to impose the no-slip boundary condition along a rigid boundary, whereas the present pIB method is to impose the inextensibility condition on the immersed boundary.
III. NUMERICAL IMPLEMENTATION
For the numerical implementation of the system(1)–(7), we use a first-order IB method, generalized to take into account the boundary Y that is linked to the boundary X by stiff springs. In the time level n t, the configurations of the two boundary representatives X and Y and the velocities u and V are all given. We update these values in the next time level (n+ 1) t by the following procedure.
Step 1. Use Eqs.(5)and(3)to calculate the force density
F defined on Lagrangian grid points and spread it into the Eulerian grid points to obtain the force density f in the fluid equations.
Step 2. Find the indicator function I (x) to define the viscosity function. Since the indicator function has the value 1 (I = 1) inside the vesicle boundary X and the value zero (I = 0) outside, it can be calculated by the following procedure [21]. Let inrepresent the interior of the vesicle and n be the unit outward normal to the vesicle; then the indicator function
−4 −2 0 2 4 −4 −3 −2 −1 0 1 2 3 4 x/R 0 x/R 0 y/R 0 y/R 0 θ −4 −2 0 2 4 −3 −2 −1 0 1 2 3
FIG. 1. The configuration of a single vesicle suspended in the shear flow u0= (u0,v0)= γ (y,0) at t = 0 (left) and t = 72 τ (right) where τ
is the intrinsic time scale. The effective radius of the vesicle R0is used as the length scale. Here the reduced area of the vesicle is V = 0.7, the
Reynolds number is Re= 0.008, and the viscosity contrast is λ = 3.5. The vesicle boundary together with streamlines shows a tank-treading tangential motion with a steady inclination angle θ between the long axis of the vesicle and the flow direction.
can be represented by
I(x)=
in
δ(x− ˜x) d ˜x. (9) By taking the gradient first to I (x) and divergence operator next to the resultant equation, we have
∇I(x) = − X δ(x− X)n ds, (10) I(x)= −∇ · X δ(x− X)n ds. (11) Thus, the indicator function can be obtained by solving the Poisson equation(11)with a singular source using a standard finite difference method. Once the indicator function is found, the viscosity can be defined as
μ(x)= μout+ (μin− μout)I (x), (12) where μin and μout are the viscosities inside and outside the vesicle, respectively. In practice, the viscosity is a smoothing function rather than a piecewise constant function across the vesicle boundary.
Step 3. Given the computed force density f and the viscosity μ, we solve the discretized version of the fluid equations(1) and(2)to update the velocity and pressure fields.
Step 4. Update the velocity U using Eq.(4)and the velocity
V using Eqs.(6)and(7). The combination of Eqs.(6)and(7)
induces a modified Helmholtz equation for the unknown function σ (s,t), which can be solved by a standard finite-difference method; see [4] for the detail.
Step 5. As the last step, we update the configuration of the two representations X and Y of the boundary using the forward Euler method with the velocities U and V. This completes the evolution over one time step.
IV. NUMERICAL RESULTS
In this section, we conduct a series of numerical tests for a single vesicle suspended in a simple shear flow. This problem has been explored by several researchers through
numerical simulations [9–12] in the Stokes flow regime. We first present the initial setup for our vesicle model and display the physical and computational parameters used in our numerical experiments.
Consider the computational domain of a rectangular box which is filled with an incompressible viscous fluid and contains an elastic vesicle in the form of a closed curve; see the left panel of Fig.1which shows the initial configuration of a vesicle. Following the analysis and scaling in Refs. [9,11,12], we use the length scale R0= L/2π where L is the perimeter of the vesicle boundary. The computational domain in our numerical study is chosen as (−D/2,D/2) × (−D,D) where
D= 8R0.
In order to apply the shear flow, it is natural to impose the background velocity field u0= (u0,v0)= γ (y,0) to the whole domain, where y ranges from−D to D and γ is the shear rate. In this paper, however, since we use the periodic boundary condition for the computational domain, we impose the following shear flow:
u0(y)= ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ γ y, if |y| H, γ HDD−H−y, if H < y D, −γ HD+y D−H, if −D y < −H, (13)
and v0= 0. The set of points {(x,y) | |x| D/2, y = ±H } represents two walls parallel to the x axis in the domain. Even though the imposed shear flow is not smooth, it generates almost the same flow as the natural shear flow u0= (u0,v0)=
γ(y,0) in|y| H; see the left panel of Fig.1.
The distance 2H between the two parallel walls is chosen as 12R0. We have checked that when the distance 2H varied from 11R0up to 14R0, the desirable shear flow is generated well and the simulation results are all similar. We define the time scale as τ = μoutR30/cbwhere μoutis the fluid viscosity outside the vesicle and cb is the bending coefficient of the vesicle. Then the time step is t= 1.5 × 10−5 in the intrinsic time units, and the spatial mesh width is h= x = y = R0/16. The number of node points representing the vesicle boundary
0 0.5 1 −0.4 −0.2 0 0.2 0.4 x/R 0 y/R 0 N=64 N=128 N=256 N=512 0.5 1 1.5 2 2.5 3 1.95 2 2.05 2.1 2.15 t/τ error ratio
FIG. 2. The configuration of the right half vesicle at t= 3 τ for the four N’s (left) and convergence ratios (defined in the text) of the computed velocity field u(x,t) (right). In the right panel, the dashed line is the convergence ratio obtained using the grids N= 64,128,256, while the solid line is the ratio obtained with N= 128,256,512. The vesicle configurations are very close for the four N’s, and the convergence ratios for the velocity are near 2 (first-order accuracy).
is 200, which makes the boundary mesh width s to be approximately half the fluid mesh width h. Using τ μoutR0as the scale of mass, the two dimensionless penalty parameters are chosen as K= 4 × 105and R= 4 × 103.
We vary the initial configuration of the vesicle, the shear rate γ , and the fluid viscosity inside the vesicle μin to see the dependence of the vesicle dynamics on these parameters. These variations, which control the interaction between the fluid flow and the vesicle motion, introduce the following dimensionless parameters:
(i) The viscosity contrast, λ= μin/μout. (ii) The reduced area,
V = A
π R02 =
4 A π
L2 , (14)
where A is the area of the vesicle. The reduced area V , which indicates how much a vesicle is deflated, is 1 for a circular vesicle (maximally swollen one) and is less than 1 for a deflated one.
(iii) The Reynolds number, Re= ργ R
2 0
μout
, (15)
which measures the inertial force versus viscous force. Throughout this paper, we vary the Reynolds number by simply adjusting the shear rate γ while keeping other physical parameters fixed.
A. Convergence study
In order to validate that the present pIB method correctly solves the inextensible vesicle dynamics in an incompressible fluid with different viscosity contrasts, we first perform a convergence study. We consider a single vesicle in a shear flow with the reduced area V = 0.7, the viscosity contrast λ = 6, and the Reynolds number Re= 0.8. Now we choose the mesh size (N,2N ) of the domain where N = 64, 128, 256, and 512 so that the corresponding mesh width is h= 8R0/N. We also choose s and t proportional to h, so that each factor of 2 in refinement of the fluid mesh width is accompanied by the same factor of refinement for the boundary mesh and the time step
duration. The penalty parameters are adjusted with the mesh size as K = 105N/64 and R= 103N/64; i.e., both parameters increase as the mesh width and time step are refined.
The left panel of Fig. 2 shows the configuration of the right-half vesicle at t = 3 τ for the four N’s. We can see that the vesicle configurations are close for the four cases of N . Especially, the difference of the vesicle configurations between the cases of N = 256 and 512 is smaller than that between the cases of other pairs with coarser resolutions, which might imply the convergence behavior of the solutions.
To get a more quantitative measure of convergence, we compare the velocity fields computed on the four different mesh widths. The right panel of Fig.2shows the convergence ratios of the computed fluid velocity u(x,t). Since we do not have the exact solution for the problem, the estimation of the convergence ratio requires three numerical solutions for three consecutive grid sizes N . We first define the discrete L2 norm of a scalar valued function ψ defined on the Cartesian grid as||ψ||2= (
i,j|ψi,j|2h2)1/2. Let (uN,vN) be the velocity field for an N× 2N Cartesian grid; then the left panel of Fig.2shows the convergence ratios (||uN− u2N||22+ ||vN− v2N||22)1/2/(||u2N− u4N||22+ ||v2N− v4N||22)1/2versus time for each of the cases N = 64 (dashed line) and N = 128 (solid line).
One can see from the figure that the convergence ratios for the fluid velocity are around 2, which indicates that the present method is first-order accurate. This is the typical behavior in accuracy for the IB method as applied to problems with thin elastic boundaries in which both the pressure and velocity are nonsmooth along the immersed boundary [24]. For a second-order IB method in the case of an immersed elastic structure of finite thickness, the interested readers can refer to [20,25,26].
B. Single vesicle under a shear flow in Stokes regime It is well known that when the viscosity contrast λ is small, the vesicle in a simple shear flow undergoes a tank-treading tangential motion at its equilibrium configuration [5,9,12–15]. This can be observed in the simulation results of Fig. 3 in which the vesicle boundary together with streamlines shows a
−4 −2 0 2 4 −2 0 2 x/R 0 x/R0 x/R0 y/R 0 y/R 0 y/R 0 −4 −2 0 2 4 −2 0 2 −4 −2 0 2 4 −2 0 2
FIG. 3. Vesicle boundary and streamlines at t= 72 τ. The reduced area of the vesicle is V = 0.7, and the Reynolds number is Re = 0.008. The viscosity contrast is λ= 1 (left), 2 (middle), and 3 (right). The vesicle boundaries show a tank-treading tangential motion with a steady inclination angle θ which depends inversely on the viscosity contrast λ.
tank-treading tangential motion with a steady inclination angle θat time t= 72 τ. Whereas the reduced area of the vesicle is V = 0.7 and the Reynolds number is Re = 0.008 in all the three cases, the viscosity contrasts are different: λ= 1 (left), 2 (middle), and 3 (right). Notice that the inclination angle θ between the long axis of the vesicle and the flow direction depends inversely on the viscosity contrast λ.
The dependence of the inclination angle θ on the viscosity contrast and the reduced volume has long been studied. Under the assumption of a fixed ellipsoidal vesicle boundary, Keller and Skalak (KS theory) have derived an evolutional equation for the inclination angle θ [5] as
dθ
dt = A + B cos(2θ), (16)
where A and B depend on the shear rate γ , the viscosity contrast λ, and the shape of the ellipsoid; see [5] for the exact formula for A and B. Even though the determination of these coefficients A and B is quite complicated in a particular situation, Eq. (16)simply states that there are two types of motions [5–8]; namely, (1) when A/B < 1, a steady state can exist, in which case the inclination angle can be calculated by θ = cos−1(−A/B)/2 (tank-treading motion); (2) when
A/B >1, the motion is no longer stationary, and the vesicle
tumbles and rotates in a shear flow (tumbling motion). We can compare the tank-treading motions of three vesicles in Fig.3and a tumbling motion of a vesicle in Fig.4. While the other parameters are kept the same, the former has a small viscosity contrast up to λ= 3, and the latter has a larger viscosity contrast, λ= 10. The left panels of Fig.4depict the configurations of a single vesicle at some chosen times which clearly shows the tumbling motion of the vesicle. Note that the time scale in Fig.4is the reciprocal of the shear rate, i.e., 1/γ , which is different from the previous figures. This time scale will be used from here throughout the paper in order to remove the dependence of the tumbling frequency on the shear rate γ and to see the inertia effect more clearly; see the next section.
When a vesicle with a fixed ellipsoidal boundary tumbles, the general solution of(16)can be found as
θ(t)= arctan A+ B √ A2− B2tan[ A2− B2(t− t 0)] . (17)
The right panel presents the time evolution of the inclination angle θ of the tumbling vesicle shown in the left panels (solid line) which agrees well with a solution of the KS theory given
−2 0 2 −2 0 2 γt=2.304 −2 0 2 −2 0 2 γ t=4.608 −2 0 2 −2 0 2 γt=6.912 −2 0 2 −2 0 2 γt=9.216 −2 0 2 −2 0 2 γt=11.52 −2 0 2 −2 0 2 γt=13.824 −2 0 2 −2 0 2 γt=16.128 −2 0 2 −2 0 2 γt=18.432 0 10 20 30 40 −π/2 0 π/2 γ t angle θ
FIG. 4. The configuration of a single vesicle suspended in a shear flow at some chosen times (left panels) and the time evolution of the inclination angle θ (right panel). The figures are scaled using the time scale 1/γ and the length scale R0. The reduced area of the vesicle is V = 0.7, the Reynolds number is Re = 0.008, and the viscosity contrast is λ = 10. We can observe a tumbling motion of the vesicle (left
panels). The time evolution of the inclination angle θ of the numerical result (solid line) agrees well with a solution of the KS theory given as(17)(dashed line).
0 1 2 3 4 5 6 7 8 0 5 10 15 20 25 30 35 θ degrees viscosity contrast λ V=0.9 V=0.8 V=0.7 V=0.6 V=0.5 Keller Skalak 2D
FIG. 5. The inclination angle θ as functions of the viscosity contrast λ for various values of the reduced area V . The simulation results (solid lines) are compared with the KS theory in 2D (dashed lines).
as(17)(dashed line) in which A and B are determined using the ellipse which has the same area and perimeter as the vesicle shown in the left panels of Fig. 4. The small deviation of our result from the KS theory in the right panel is due to the fact that the KS theory has the assumption of a fixed elliptical vesicle boundary, which is apparently not satisfied in our vesicle model.
According to the KS theory, the coefficient B is a function of the viscosity contrast λ and the ellipsoidal shape which determines the reduced area V (volume in 3D). In the tank-treading regime, when the viscosity contrast λ increases or the reduced area V decreases, the coefficient B decreases, so the steady inclination angle θ decreases. This fact can be verified from Fig.5in which we plot the steady inclination angle in terms of the viscosity contrast λ for various reduced areas. In the figure, we compare our simulation results (solid lines) and the KS theory in 2D (dashed lines). Although both results capture the generic dependence of the inclination angle θ on viscosity contrast λ and reduced area V , there are
some quantitative discrepancies between the theory and the numerical results. Despite this, Fig.5is quite comparable to Fig.5in Ref. [12] which uses a different numerical method but observes the same type of discrepancies.
For each fixed reduced area V , as the viscosity contrast λ increases in the tank-treading regime, the inclination angle θ decreases as shown in Fig.5. If λ increases further to reach some critical value λc, there occurs a transition from tank-treading to tumbling. Figure6shows the prediction of λcusing the numerical simulations (solid line) in comparison with the KS theory (dashed line). The KS theory(16)predicts that the critical viscosity contrast is achieved at A= −B, and also that, close to this transition point, the angle θ is related to the viscosity contrast λ with a square root law. We use the numerical results shown in Fig. 5to interpolate around θ = 0 with a square root function, and to extrapolate it at θ = 0 to obtain the critical viscosity contrast λc. Though both the predictions produce qualitatively the same trend of the critical viscosity contrast versus the reduced area, it is a bit
0.5 0.6 0.7 0.8 0.9 1 2 3 4 5 6 7 8 reduced area V
critical viscosity contrast
λ c
Keller Skalak 2D Numerical result
−2 0 2 −2 −1 0 1 2 γt=9.60 −2 0 2 −2 −1 0 1 2 γt=11.52 −2 0 2 −2 −1 0 1 2 γt=13.44 −2 0 2 −2 −1 0 1 2 γt=15.36 −2 0 2 −2 −1 0 1 2 γt=17.28 −2 0 2 −2 −1 0 1 2 γt=19.20 −2 0 2 −2 −1 0 1 2 γt=21.12 −2 0 2 −2 −1 0 1 2 γt=23.04
FIG. 7. The configuration of a single vesicle in a shear flow at some chosen times with the time scale 1/γ and the length scale R0. The
physical parameters are all the same as those used in Fig.4except for the Reynolds number. The Reynolds number here is Re= 0.8 which is 100 times higher than that of Fig.4. One can see a clear tumbling motion with an ample vesicle deformation.
underestimated by the theory; see Fig.6 which is also quite comparable to Fig.6in Ref. [12].
C. Single vesicle under a shear flow with O(1) Reynolds number Most of the theoretical and numerical studies on the transi-tion between tank-treading and tumbling motransi-tions of a vesicle have been so far restricted to the Stokes regime by neglecting the inertia effect [5–12]. Also available experimental data on the vesicle dynamics correspond to a very low Reynolds number limit [13–16]. The vesicle dynamics can be a model of a red blood cell interacting with blood flow which may have the Reynolds number of order 1 [27,28]. In this case, one cannot ignore the inertia effect and needs to consider the full Navier-Stokes equations. Notice that the present pIB method is not restricted to Stokes flow and can handle flows with a high Reynolds number which can be realized simply by increasing the shear rate; see Eq.(15).
Figure7depicts a vesicle motion in a flow with 100 times higher Reynolds number than the previous simulations. The reduced area of the vesicle is V = 0.7 and the viscosity contrast is λ= 10. These two values and other parameters are the same as those of the vesicle shown in Fig.4 except for the Reynolds number. Though both vesicles go through the same qualitative behavior of tumbling motion, the vesicle with Re= 0.8 (Fig.7) shows an ample vesicle deformation as compared to the one with Re= 0.008 (Fig.4).
Figure8shows the inclination angle θ in terms of time for various Reynolds numbers. Since the vesicle may not maintain an elliptical boundary in a high Reynolds number flow, we define the inclination angle as the angle between the x axis and the longest line segment among the line segments from the center of the vesicle to the vesicle boundary. In order to remove the dependence of the frequency on the shear rate γ
and to see the inertia effect more clearly, we use as the time scale the reciprocal of the shear rate 1/γ . The reduced area of the vesicle is V = 0.7 and the viscosity contrast is λ = 15.
The upper panel of Fig.8shows that when the Reynolds numbers are low, i.e., Re= 0.008, 0.4, and 1.2, the vesicle undergoes a tumbling motion, with the tumbling frequency getting lower as the Reynolds number gets higher. The large difference of the tumbling frequencies without collapsing indicates that there exists an apparent inertia effect which is to prevent a vesicle from tumbling. When the Reynolds number becomes even higher, i.e., Re= 1.6 and 2 (lower panel), the inertia effect gets higher and the vesicle stops tumbling and returns to the tank-treading motion with an almost steady inclination angle θ . Also the steady inclination angle θ in these high Reynolds number flows is generally larger than those of tank-treading vesicles in the Stokes regime; compare the inclination angles in Figs.8and5.
Table I shows the tumbling frequency when the vesicle with V = 0.7 is in the tumbling regime (TB) and the steady inclination angle θ when it is in the tank-treading regime (TT) for various Reynolds numbers and viscosity contrasts. In the tumbling regime, as the Reynolds number gets lower, the tumbling frequency becomes higher; i.e., the vesicle tumbles more frequently. In fact, a higher Reynolds number implies a higher shear rate and thus a higher moment generated, which forces the vesicle to rotate more frequently in the real time. However, when we use the dimensionless time γ× time to remove the dependence on the shear rate, as the Reynolds number gets higher, the inertia effect gets stronger and the tumbling frequency becomes lower as shown in TableI.
Also in the tumbling regime, as the viscosity contrast gets larger, the tumbling frequency becomes higher. When the internal fluid gets more viscous, the vesicle behaves more as if quasirigid, which makes the torque for rotational motion
0 10 20 30 40 50 −π/2 0 γt angle θ Re=1.2 0 20 40 60 80 100 120 0 π/8 π/4 3π/8 π/2 γ t angle θ Re=1.6 Re=2
FIG. 8. The inclination angle θ in terms of dimensionless time for various Reynolds numbers. The reduced area of the vesicle is V= 0.7 and the viscosity contrast is λ= 15. When the Reynolds numbers are 0.008, 0.4, and 1.2 (upper panel), the vesicle undergoes a tumbling motion, with the frequency getting lower as the Reynolds number gets higher. When the Reynolds numbers become even higher, i.e., Re= 1.6 and 2 (lower panel), the tumbling motion is inhibited, and the vesicle returns to tank-treading motion with an almost steady inclination angle θ . more effective. For all the viscosity contrasts considered
here, however, we surprisingly observe that a further increase of Reynolds number inhibits tumbling motion in favor of tank-treading motion; see each row of TableI. This inhibition of tumbling motion in a flow with high Reynolds number agrees well with the recent results in Ref. [28].
Although a tank-treading motion is favorable in a high Reynolds number flow, as the viscosity contrast λ increases to reach some critical value λc, there occurs a transition from tank-treading to tumbling even in a high Reynolds number flow. Figure 9 compares the tumbling motions of vesicles with the reduced area V = 0.5 in two different situations: The vesicle in the upper panels has Re= 0.008 and λ = 2.75, and the vesicle in the lower panels has Re= 1.6 and λ = 55. The times are selected when the vesicle has similar inclination angles for the two cases. Both the vesicles initially have the
same configuration which includes a concave part around their center; see the upper panels. Whereas the concave part is kept in the low Reynolds number flow (upper), it is almost straightened in the high Reynolds number flow (lower). Even though vesicles with a high Reynolds number and a large viscosity contrast go through a large deformation as shown in the lower panels, their length and interior area are preserved well, and the velocity field maintains the desirable characteristic of the shear flow. This implies that the the simulations with a high Reynolds number and a large viscosity contrast are stable.
Figure 10 shows the numerical prediction of critical viscosity contrast λc versus the reduced area V for various Reynolds numbers. Given a reduced area, as the Reynolds number increases, the critical viscosity contrast λc for the tank-treading/tumbling transition increases as well. Whereas
TABLE I. The tumbling frequency in the tumbling regime (TB) and the steady inclination angle θ in the tank-treading regime (TT) for various Reynolds numbers Re and viscosity contrasts λ.
Re= 0.008 0.4 0.8 1.6 2.0
λ= 5 5.36× 10−2(TB) 10.5 (TT) 21.4 (TT) 23.9 (TT) 25.6 (TT)
10 8.79× 10−2(TB) 6.79× 10−2(TB) 12.8 (TT) 17.9 (TT) 20.6 (TT)
15 9.65× 10−2(TB) 8.10× 10−2(TB) 4.40× 10−2(TB) 11.5 (TT) 17.1 (TT)
−2 0 2 −2 −1 0 1 2 γt=6.336 Re=0.008; λ =2.75 −2 0 2 −2 −1 0 1 2 γt=29.376 −2 0 2 −2 −1 0 1 2 γt=46.08 −2 0 2 −2 −1 0 1 2 γt=47.232 −2 0 2 −2 −1 0 1 2 γt=74.88 Re=1.6; λ =55 −2 0 2 −2 −1 0 1 2 γt=176.64 −2 0 2 −2 −1 0 1 2 γt=180.48 −2 0 2 −2 −1 0 1 2 γt=182.40
FIG. 9. The comparison of the tumbling motion of a vesicle with the reduced area V = 0.5 at the times when the vesicles have similar inclination angles. The vesicle in the upper panels has Re= 0.008 and λ = 2.75, and the vesicle in the lower panels has Re = 1.6 and λ = 55. Whereas the vesicle in a low Reynolds number flow keeps its initial concave part, the vesicle in a high Reynolds number flow loses it and experiences a large deformation.
the increment is so small that the critical viscosity contrasts are almost indistinguishable when the Reynolds number is below 0.08, it becomes pronounced with Re > 0.08, in which case, the increment of critical viscosity contrast is steeper for vesicles with small reduced areas than for vesicles with large reduced areas.
It is also interesting to note that for Reynolds numbers Re= 0.8, 1.2, and 1.6, the critical viscosity contrast decreases as the reduced area increases up to about V = 0.9. This behavior is completely different from that in the case of Reynolds numbers
Re= 0.008,0.08, and 0.4, where the critical viscosity contrast is an increasing function of the reduced area. The latter behav-ior was generally observed in previous studies on the transition of vesicle dynamics in a shear flow [5–12,15,16]. We attribute this significantly different behavior between the present simulation results and the previous results in the literature to the difference of Reynolds numbers and thus the inertial effects.
In summary, whereas a flattened vesicle with a small reduced area V tumbles more easily than a swollen one
0.5 0.6 0.7 0.8 0.9 1 5 10 15 20 25 30 35 40 45 50 55 reduced area V
critical viscosity contrast
λc Re=0.008 Re=0.08 Re=0.4 Re=0.8 Re=1.2 Re=1.6
FIG. 10. The critical viscosity contrast λc versus the reduced area for various Reynolds numbers. For a low Reynolds number flow
(Re 0.4), the critical viscosity contrast increases as the reduced area increases. However, for Reynolds number Re 0.8, the critical viscosity contrast decreases as the reduced area increases up to about V = 0.9. For each reduced area, the higher the Reynolds number gets, the higher the critical viscosity contrast is.
[O(1)] Reynolds number flow. This unexpected change of the tumbling tendency with respect to the reduced area is discovered in our simulations for the first time, to the best of our knowledge. Even though the inhibition of tumbling motion in a flow with high Reynolds number was recently observed in Ref. [28], the authors in Ref. [28] focused on and simulated a vesicle with one reduced area V = 0.82 only. Here, however, we not only verify the inhibition of tumbling motion by the inertia effect for vesicles with various reduced areas from V = 0.5 to 0.95, but we also show that the inhibition effect is much stronger for flattened vesicles than for swollen vesicles in a high [O(1)] Reynolds number flow.
V. CONCLUSIONS
We have extended the penalty immersed boundary (pIB) method to simulate the transitional dynamics of inextensible vesicle under a shear flow with different viscosity contrast. We validate the method by performing a convergence study, which confirms the consistent first-order accuracy for the velocity field. We have also simulated several tests for two-dimensional
present pIB method are comparable to those in literature. In order to consider more realistic situations, we have investigated the transitional dynamics of a vesicle under shear flow with O(1) Reynolds number. The inertial effect inhibits the tumbling motion in favor of the tank-treading motion. This inhibition effect by the inertia is stronger especially when the reduced area is small.
Future work will include the extension of the present methodology to study the three-dimensional vesicle dynamics. The application of the present method to some biological systems such as red blood cells and drug-carrying capsules will also be the subject of future work.
ACKNOWLEDGMENTS
Y.K. is supported by a National Research Foundation of Korea grant funded by the Korean Government (Grant No. 2012R1A1A2043238). M.-C.L. is supported in part by the National Science Council of Taiwan under Research Grant No. NSC-98-2115-M-009-014-MY3 and by the NCTS in Taiwan.
[1] C. S. Peskin,Acta Numerica 11, 479 (2002).
[2] C. S. Peskin and D. M. McQueen,J. Comput. Phys. 81, 372 (1989).
[3] C. S. Peskin and D. M. McQueen, in Case Studies in
Mathemat-ical Modeling: Ecology, Physiology, and Cell Biology (Prentice
Hall, Englewood Cliffs, NJ, 1996), pp. 309–337.
[4] Y. Kim and M.-C. Lai,J. Comput. Phys. 229, 4840 (2010). [5] S. R. Keller and R. Skalak,J. Fluid Mech. 120, 27 (1982). [6] F. Rioual, T. Biben, and C. Misbah,Phys. Rev. E 69, 061914
(2004).
[7] C. Misbah,Phys. Rev. Lett. 96, 028104 (2006).
[8] M.-A. Mader, H. Ez-Zahraouy, C. Misbah, and T. Podgorski, Eur. Phys. J. E 22, 275 (2007).
[9] M. Kraus, W. Wintz, U. Seifert, and R. Lipowsky,Phys. Rev. Lett. 77, 3685 (1996).
[10] H. Zhou and C. Pozrikidis,J. Fluid Mech. 283, 175 (1995). [11] S. K. Veerapaneni, D. Gueyffier, D. Zorin, and G. Biros,
J. Comput. Phys. 228, 2334 (2009).
[12] J. Beaucourt, F. Rioual, T. Seon, T. Biben, and C. Misbah,Phys. Rev. E 69, 011906 (2004).
[13] K. H. de Haas, C. Blom, D. van den Ende, M. H. G. Duits, and J. Mellema,Phys. Rev. E 56, 7132 (1997).
[14] V. Kantsler and V. Steinberg,Phys. Rev. Lett. 95, 258101 (2005).
[15] M.-A. Mader, V. Vitkova, M. Abkarian, A. Viallat, and T. Podgorski,Eur. Phys. J. E 19, 389 (2006).
[16] V. Kantsler and V. Steinberg,Phys. Rev. Lett. 96, 036001 (2006). [17] H. Noguchi and G. Gompper,Proc. Natl. Acad. Sci. USA 102,
14159 (2005).
[18] C. Pozrikidis,J. Fluid Mech. 216, 231 (1990). [19] M. J. Stevens,J. Chem. Phys. 121, 11942 (2004). [20] Y. Kim and C. S. Peskin,Phys. Fluids 19, 053103 (2007). [21] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi,
W. Tauber, J. Han, S. Nas, and Y.-J. Jan,J. Comput. Phys. 169, 708 (2001).
[22] D. Goldstein, R. Handler, and L. Sirovich,J. Comput. Phys. 105, 354 (1993).
[23] E. M. Saiki and S. Biringen, J. Comput. Phys. 123, 450 (1996).
[24] K. Chen, K. Feng, Y. Kim, and M.-C. Lai,J. Comput. Phys. 230, 4377 (2011).
[25] M.-C. Lai and C. S. Peskin,J. Comput. Phys. 160, 705 (2000). [26] B. E. Griffith and C. S. Peskin,J. Comput. Phys. 208, 75 (2005). [27] Y.-C. Fung, Biomechanics: Motion, Flow, Stress, and Growth
(Springer-Verlag, New York, Berlin, Heidelberg, 1990). [28] A. Laadhari, P. Saramito, and C. Misbah, Phys. Fluids 24,
