4. Implementation details
5.2. Convergence test for the Stokes flow with an inextensible interface enclosing a solid particle. In this subsection, we perform a convergence study
Ω
pe(x) dx.
So our initial guess p0ijin the GMRES iteration can be chosen as p0ij=#
Ωpe(x) dx/|Ω|.
In those tests, the tolerance of the residual is chosen as 10−8.
Table 5.1 shows the maximum errors between the exact and numerical solutions for different grid resolutions. One can see that the velocity field has clear second-order accuracy, while the pressure has first-order accuracy. Table 5.2 shows the efficiency of the present Stokes solver. One can see that the number of iterations increases slightly even we double the grid size, and the CPU time for a 512× 512 grid is just a few seconds.
5.2. Convergence test for the Stokes flow with an inextensible interface enclosing a solid particle. In this subsection, we perform a convergence study for the present numerical algorithm to the Stokes flow with an inextensible interface enclosing a solid particle. Here, we put an inextensible interface Γ and particle P with initial configuration X(s) = (0.25 cos(s), 0.5 sin(s)) and Y(s) = (0.1 cos(s), 0.1 sin(s)) under a shear flow (u, v) = (χy, 0) in a fluid domain Ω. The dimensionless shear rate χ is chosen to be χ = 1. We also choose the different grid size as m = n = 32, 64, 128, 256, 512 in which the corresponding mesh width is h = 2/m. We also set the Lagrangian mesh widths to be Δs≈ h/2 and Δα ≈ h/2, and the time step size to Δt = h/4.
Since the analytical solution is not available in this test, we choose the result obtained from the finest mesh m = n = 512 as our reference solution and compute the
Downloaded 04/28/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
B704 MING-CHIH LAI, WEI-FAN HU, AND WEN-WEI LIN Table 5.3
The mesh refinement results for the perimeter of the interfaceLh, the interface configuration Xh, and the velocityuhandvh.
m = 32 m = 64 rate m = 128 rate m = 256 rate
|Lh− L0|/L0 3.511e-2 2.089e-2 0.75 1.220e-2 0.78 6.749e-3 0.85
Xh− Xref∞ 3.413e-3 1.079e-3 1.66 4.438e-4 1.28 1.622e-4 1.45
uh− uref∞ 6.552e-2 2.592e-2 1.34 1.053e-2 1.30 4.047e-3 1.38
vh− vref∞ 1.133e-1 3.660e-2 1.63 1.554e-2 1.24 4.968e-3 1.64
maximum error between the reference and the numerical solutions. All the numerical solutions are computed up to time T = 0.0625. Since the interface is inextensible, the perimeter of the interface should remain a constant theoretically as time evolves.
Let L0 and Lh be the perimeters of the interface at the initial time and final time T = 0.0625, respectively. The relative error of the perimeter is defined as |Lh− L0|/L0. Table 5.3 shows the relative errors of the perimeter, the maximum errors of the interface configuration, and the maximum errors for the fluid velocity field.
Note that the fluid variables are defined at the staggered grid, so when we refine the mesh, the numerical solutions will not coincide with the same grid locations.
In these runs, we simply use a linear interpolation to compute the solutions at the desired locations. Due to the fact that the IB formulation has the singular forcing term in the equations, regularizing the singular term by smoothing the discrete delta function causes the method to be first-order accurate. As shown in [3], even though the calculated pressure has O(1) error near the interface, its gradient has the same discretization error as the singular delta force term near the interface. Thus, the overall accuracy of the velocity field in the IB method is first-order accurate in the L∞ norm. The numerical results shown in Table 5.3 are consistent with what we expect from theory.
5.3. Tank-treading to tumbling motion under shear flow. Unlike the pre-vious subsection where we focus on the numerical convergence test for our present scheme, here we consider the physical transient deformation of an inextensible inter-face with or without an enclosing rigid particle in a simple shear flow. As mentioned before, the motivation of this test is to simulate the compound vesicle dynamics, which have a lot of applications in bio-fluid problems [26].
In the present numerical experiments, we choose the residual tolerance for GM-RES as 10−4, which is larger than the pure Stokes flow. This is because the elastic tension σ tends to oscillate along the interface, which makes the GMRES method solving (4.11) difficult to converge if the residual tolerance is too small. It will be more appealing if one can find a good preconditioner for solving linear system (4.11) so that GMRES method can converge faster. However, this issue has not yet been resolved.
We first consider the case of an inextensible interface without enclosing a solid particle in a simple shear flow. It is wellknown that the equilibrium dynamics of an inextensible interface or vesicle under a simple shear flow will undergo a tanking-treading motion if the viscosity contrast is under a certain threshold [8]. Here, by tank-treading motion we mean that the configuration of the interface remains stationary while there is a tangential motion along the interface. As studied in previous literature [28, 9, 10, 25, 7], the motion of this steady interface can be characterized by both the inclination angle θ between the long axis of interface and the flow direction, and the
Downloaded 04/28/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
AN INEXTENSIBLE INTERFACE WITH A SOLID PARTICLE B705
Fig. 5.1. The inclination angles θ/π (left) and the tank-treading frequency f = 2π/
Γ dl u · τ (right) versus reduced areasV with different dimensionless shear rates for the tank-treading motion of an inextensible interface in a shear flow.
tank-treading frequency f = 2π/#
Γ dl
uτ of the revolution, where uτ is the tangential velocity component. The inclination angle has been found to be strongly dependent on the reduced area V = 4πA0
L20 , where A0 is the enclosed area of the interface and L0 is the total length of interface. (Notice that, by the above definition, a circle has the reduced area V = 1, while an ellipse with larger aspect ratio has a smaller reduced area.) However, the inclination angle is independent of the dimensionless shear rate χ. This behavior is verified in the left panel of Figure 5.1, which shows the steady inclination angle (θ/π) versus the reduced area (V ) with different shear rates χ = 1, 5, 10, and it is in a good agreement with previous studies [28, 9, 10, 25, 7]. As the reduced area increases, the inclination angle increases as well. The right panel of Figure 5.1 shows the tank-treading frequency f versus the reduced area V for different shear rates. One can see that as the dimensionless shear rate increases, the tangential motion becomes stronger; thus, the frequency becomes larger. Moreover, by fixing the shear rate, if the interface has larger reduced area, then it has larger frequency as well. Again, our numerical results are in good agreement with previous studies in the literature.
Now we consider the case of an inextensible interface enclosing a solid particle (a compound interface) under a shear flow. As shown in [26], the compound interfacial dynamics will have the transition from tank-treading (TT) to tumbling (TB) if the inclusion effect is strong enough. That is, if the filling fraction φ = a2π/A0 of the particle (a is the inclusion particle radius, A0is the interface enclosing area) is above some critical threshold, then the interface will start to tumble rather than being sta-tionary. This can be explained as follows. By including a solid particle, the energy dissipation enhances, so the compound interface behaves like an inclusion-free inter-face that encapsulates a more viscous fluid. The larger the inclusion is, the higher the viscosity will be. Figures 5.2 and 5.3 show the time evolutionary plots for interfacial configurations with different filling fractions, φ = 0.08 and φ = 0.42, respectively. One can indeed see that in the smaller filling fraction case in Figure 5.2, the compound interface has TT motion, while in the larger filling fraction case in Figure 5.3, it has the TB motion.
Meanwhile, in the TT regime, as the filling fraction increases, both the inclination angle and the TT frequency will decrease. Besides, as in the inclusion-free case, the
Downloaded 04/28/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
B706 MING-CHIH LAI, WEI-FAN HU, AND WEN-WEI LIN
Fig. 5.2. The motion of an compound interface in a shear flow with initial configuration X(s) = (0.25 cos(s), 0.5 sin(s)) and Y(s) = (0.1 cos(s), 0.1 sin(s)). φ = 0.08.
Fig. 5.3. The motion of an compound interface in a shear flow with initial configuration X(s) = (0.25 cos(s), 0.5 sin(s)) and Y(s) = (0.23 cos(s), 0.23 sin(s)). φ = 0.42.
compound interface with larger reduced area has larger inclination angle and TT frequency when the filling fraction is small. Both behaviors can be confirmed in Figure 5.4 and are consistent qualitatively with those in [26].
We also investigate the critical value of filling fraction versus the reduced area for the TT motion to TB transition, as in Figure 5.5. Above the critical value, the interface motion will transit from TT to TB. One can easily see that as the reduced area increases, the critical filling fraction increases too. This is exactly the similar behavior as the inclusion-free case in which the critical viscosity contrast increases with the reduced area.
Downloaded 04/28/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
AN INEXTENSIBLE INTERFACE WITH A SOLID PARTICLE B707
0 0.05 0.1 0.15 0.2
filling fraction φ
V=0.84 V=0.95
0 0.05 0.1 0.15 0.2
0.25 0.3 0.35 0.4 0.45 0.5
filling fraction φ
f
V=0.84 V=0.95
Fig. 5.4. The inclination angles θ/π (left) and the tank-treading frequency f = 2π/
Γ dl u · τ versus filling fractionφ the in tank-treading regime.
0.6 0.65 0.7 0.75 0.8 0.85 0.9
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55
reduced area
filling fraction φ
TB
TT
Fig. 5.5. The critical filling fraction for the TT to TB transition versus the reduced area.
6. Conclusions. Biological fluid mechanics often require one to study the dy-namics of a lipid membrane encapsulating some cellular contents. For instance, a compound vesicle usually consists of a lipid bilayer membrane enclosing a fluid with a suspended particle. Since the bilayer membrane resists area dilation and its thick-ness is on the order of a nanometer, one can regard the membrane as an inextensible (incompressible) interface. In this paper, we consider the problem of an inextensi-ble interface enclosing a solid particle in Stokes flow. We first write the governing equations in the immersed boundary (IB) framework, in which the interface and solid boundary are treated as force generators in the fluid. In additional to solving for the fluid variables, the present problem involves finding an extra unknown elastic ten-sion such that the surface divergence of the velocity is zero along the interface, and an extra unknown particle surface force such that the velocity satisfies the no-slip boundary condition along the particle surface. We show that the spreading operator of the tension and the surface divergence operator of the velocity are skew-adjoint mathematically. While the interface moves with local fluid velocity, the enclosed par-ticle hereby undergoes a rigid body motion, and the system is closed by the force-free and torque-free conditions along the particle surface.
Downloaded 04/28/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
B708 MING-CHIH LAI, WEI-FAN HU, AND WEN-WEI LIN
The governing equations are then discretized by standard centered difference schemes on a staggered grid, and the interactions between the interface and parti-cle and the fluid are discretized using discrete delta function as in the IB method.
The discrete spreading operator of the tension and the discrete surface divergence op-erator of the velocity used in the present scheme preserve the skew-adjoint property numerically. The resultant linear system of equations is therefore symmetric and can be solved by fractional steps so that only fast Poisson solvers are involved. The present method can be extended to Navier–Stokes flow with moderate Reynolds number by treating the nonlinear advection terms explicitly for the time integration. It is also important to mention that, unlike our previous work using the penalty approach [7], we are able to estimate the local error of inextensibility for two successive time steps in the present scheme. In addition, since there are no penalty parameters introduced as in [7], the time step size can be significantly increased.
We study the tank-treading (TT) and tumbling (TB) dynamics for an inextensible interface enclosing a solid particle with different filling fractions under shear flow. We have found that by the increase of the filling fraction, the interface tends to transit from TT to TB. In the TT regime (when the filling fraction is small), the inclination angle and TT frequency will increase as the reduced area increases. Those angle and frequency will decrease as the filling fraction increases. We also have found that the critical filling fraction (from TT to TB) will increase as the reduced area increases, which is qualitatively consistent with the results in [26].
Appendix A. In this appendix, we first give a simple derivation for the rigid particle motion equation (2.5) and the torque-free condition (2.7) in twodimensions.
These can be derived directly from the following equations in threedimensions, V = Vc+ω × (Y − Yc),
∂P
F× (Y − Yc) dP = 0,
by setting the third component of the vector-valued functions to be zero. For instance, we set V = (V1, V2, 0), so the vorticity vectorω = (0, 0, ω).
The symmetric property of the (3.6)–(3.8) can be seen as follows:
V1k = V1c− ω(Y2k− Y2c), V2k = V2c+ ω(Y1k− Y1c),
We can scale out the coefficient Δα so that the above three equations can be written in the matrix form
Downloaded 04/28/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
AN INEXTENSIBLE INTERFACE WITH A SOLID PARTICLE B709 Appendix B. In this appendix, we derive that the matrix obtained from the discrete spreading operator Sh of σh and the matrix obtained from discrete surface divergence operator ∇sh of U are transpose with each other. To see this, we first rewrite the operator Sh(σh) as
Sh(σh) =
M−1
k=0
Ds(στ )k δh(x− Xk) Δs
=
M−1
k=0
σk+1/2τk+1/2− σk−1/2τk−1/2
Δs δh(x− Xk)Δs
=
M k=1
σk−1/2τk−1/2δh(x− Xk−1)−
M k=1
σk−1/2τk−1/2δh(x− Xk)
− σ−1/2τ−1/2δh(x− X0) + σM−1/2τM−1/2δh(x− XM)
=
M k=1
(δh(x− Xk−1)− δh(x− Xk))τk−1/2σk−1/2.
Note that the last two terms are cancelled out due to the periodicity of the interface.
Now we can write the discrete operator∇sh as
∇sh· Uk =Uk− Uk−1
Δs · τk−1/2
|DsX|k−1/2
=− h2
Δs |DsX|k−1/2
x
(δh(x− Xk−1)− δh(x− Xk))τk−1/2· ui,j. Since the discrete surface divergence operator of the velocity is zero, we can scale out the coefficient −Δs |DshX|2
k−1/2 so that the resultant matrices obtained from Sh and
∇sh· are transpose to each other.
Acknowledgments. We thank Xiaofan Li for helping with the formulation of the motion of the solid particle.
REFERENCES
[1] J. Adams, P. Swarztrauber, and R. Sweet, Fishpack, a package of Fortran subpro-grams for the solution of separable elliptic partial differential equations, 1980; see http://www2.cisl.ucar.edu/resources/legacy/fishpack.
[2] J. H. Bramble, J. E. Pasciak, and A. T. Vassilev, Analysis of the inexact Uzawa algorithm for saddle point problems, SIAM J. Numer. Anal., 34 (1997), pp. 1072–1092.
[3] K.-Y. Chen, K.-A. Feng, Y. Kim, and M.-C. Lai, A note on pressure accuracy in immersed boundary method for Stokes flow, J. Comput. Phys., 230 (2011), pp. 4377–4383.
[4] F. H. Harlow and J. E. Welsh, Numerical calculation of time-dependent viscous incompress-ible flow of fluid with a free surface, Phys. Fluids, 8 (1965), pp. 2181–2189.
[5] D. Goldstein, R. Handler, and L. Sirovich, Modeling a no-slip flow boundary with an external force field, J. Comput. Phys., 105 (1993), pp. 354–366.
[6] K. H. de Haas, C. Blom, D. van den Ende, M. H. G. Duits, and J. Mellema, Deformation of giant lipid bilayer vesicles in shear flow, Phys. Rev. E, 56 (1997), pp. 7132–7137.
[7] Y. Kim and M.-C. Lai, Simulating the dynamics of inextensible vesicles by the penalty im-mersed boundary method, J. Comput. Phys., 229 (2010), pp. 4840–4853.
[8] S. R. Keller and R. Skalak, Motion of a tank-treading ellipsoidal particle in a shear flow, J. Fluid Mech., 120 (1982), pp. 27–47.
[9] V. Kantsler and V. Steinberg, Orientation and dynamics of a vesicle in tank-treading motion in shear flow, Phys. Rev. Lett., 95 (2005), 258101.
Downloaded 04/28/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
B710 MING-CHIH LAI, WEI-FAN HU, AND WEN-WEI LIN
[10] M. Kraus, W. Wintz, U. Seifert, and R. Lipowsky, Fluid vesicles in shear flow, Phys. Rev.
Lett., 77 (1996), pp. 3685–3688.
[11] Z. Li and M.-C. Lai, New finite difference methods based on IIM for inextensible interfaces in incompressible flows, East Asian J. Appl. Math., 1 (2011), pp. 155–171.
[12] M.-C. Lai and C. S. Peskin, An immersed boundary method with formal second order accuracy and reduced numerical viscosity, J. Comput. Phys., 160 (2000), pp. 705–719.
[13] M.-C. Lai, Y.-H. Tseng, and H. Huang, An immersed boundary method for interfacial flow with insoluble surfactant, J. Comput. Phys., 227 (2008), pp. 7279–7293.
[14] H. Noguchi and G. Gompper, Shape transitions of fluid vesicles and red blood cells in capillary flows, Proc. Natl. Acad. Sci. USA, 102 (2005), pp. 14159–14164.
[15] J. B. Perot, An analysis of the fractional method, J. Comput. Phys., 108 (1993), pp. 51–58.
[16] C. S. Peskin, The immersed boundary method, Acta Numer., 11 (2002), pp. 479–517.
[17] C. Pozrikidis, The axisymmetric deformation of a red blood cell in uniaxial straining Stokes flow, J. Fluid Mech., 216 (1990), pp. 231–254.
[18] J. Peters, V. Reichelt, and A. Reusken, Fast iterative solvers for discrete Stokes equations, SIAM J. Sci. Comput., 27 (2005), pp. 646–666.
[19] E. M. Saiki and S. Biringen, Numerical simulation of a cylinder in uniform flow: Application of a virtual boundary method, J. Comput. Phys., 123 (1996), pp. 450–465.
[20] M. J. Stevens, Coarse-grained simulations of lipid bilayers, J. Chem. Phys., 121 (2004), pp. 11942–11948.
[21] G. W. Schmid-Schonbein, Y. Y. Shih, and S. Chien, Morphometry of human leukocytes, Blood, 56 (1980), pp. 866–875.
[22] J. S. Sohn, Y.-H. Tseng, S. Li, A. Voigt, and J. S. Lowengrub, Dynamics of multicompo-nent vesicles in a viscous fluid, J. Comput. Phys., 229 (2010), pp. 119–144.
[23] J. C. Strikwerda, An iterative method for solving finite difference approximations to the Stokes equations, SIAM J. Numer. Anal., 21 (1984), pp. 447–458.
[24] K. Taira and T. Colonius, The immersed boundary method: A projection approach, J. Com-put. Phys., 225 (2007), pp. 2118–2137.
[25] S. K. Veerapaneni, D. Gueyffier, D. Zorin, and G. Biros, A boundary integral method for simulating the dynamics of inextensible vesicles suspended in a viscous fluid in 2D, J.
Comput. Phys., 228 (2009), pp. 2334–2353.
[26] S. K. Veerapaneni, Y.-N. Young, P. M. Vlahovska, and J. Blawzdziewicz, Dynamics of a compound vesicle in shear flow, Phys. Rev. Lett., 106 (2011), 158103.
[27] X. Yang, X. Zhang, Z. Li, and G.-W. He, A smoothing technique for discrete delta func-tions with application to immersed boundary method in moving boundary simulafunc-tions, J.
Comput. Phys., 228 (2009), pp. 7821–7836.
[28] H. Zhou and C. Pozrikidis, Deformation of liquid capsules with incompressible interfaces in simple shear flow, J. Fluid Mech., 283 (1995), pp. 175–200.
Downloaded 04/28/14 to 140.113.38.11. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php