New Finite Difference Methods Based on IIM for Inextensible Interfaces in Incompressible Flows

(1)

New Finite Difference Methods Based on IIM for Inextensible

Interfaces in Incompressible Flows

Zhilin Li1 and Ming-Chih Lai2

Zhilin Li: zhilin@unity.ncsu.edu; Ming-Chih Lai: mclai@math.nctu.edu.tw

1_{Center for Research in Scientific Computation & Department of Mathematics, North Carolina}

State University, Raleigh, NC 27695-8205, USA.

2_{Center of Mathematical Modeling and Scientific Computing & Department of Applied}

Mathematics, National Chiao Tung University, 1001 Ta Hsueh Road, Hsinchu 300, Taiwan.

Abstract

In this paper, new finite difference methods based on the augmented immersed interface method (IIM) are proposed for simulating an inextensible moving interface in an incompressible two-dimensional flow. The mathematical models arise from studying the deformation of red blood cells in mathematical biology. The governing equations are incompressible Stokes or Navier-Stokes equations with an unknown surface tension, which should be determined in such a way that the surface divergence of the velocity is zero along the interface. Thus, the area enclosed by the interface and the total length of the interface should be conserved during the evolution process. Because of the nonlinear and coupling nature of the problem, direct discretization by applying the immersed boundary or immersed interface method yields complex nonlinear systems to be solved. In our new methods, we treat the unknown surface tension as an augmented variable so that the augmented IIM can be applied. Since finding the unknown surface tension is essentially an inverse problem that is sensitive to perturbations, our regularization strategy is to introduce a controlled tangential force along the interface, which leads to a least squares problem. For Stokes equations, the forward solver at one time level involves solving three Poisson equations with an interface. For Navier-Stokes equations, we propose a modified projection method that can enforce the pressure jump condition corresponding directly to the unknown surface tension. Several numerical experiments show good agreement with other results in the literature and reveal some interesting phenomena.

Keywords

Inextensible interface; incompressible flow; Stokes equations; Navier-Stokes equations; immersed interface method; inverse problem; regularization; augmented immersed interface method

1. Introduction

In this paper, we develop some new finite difference methods based on the augmented immersed interface method (cf. for example [2, 7, 10, 12]) for simulating an inextensible moving interface in an incompressible shear flow. The problem involves finding an unknown surface tension σ (s, t) such that the surface divergence of the velocity is zero along the interface. Since the fluid is incompressible and the interface is inextensible, both the area enclosed by the interface and the total length of the interface should be conserved. The mathematical model has been used to describe the deformation of erythrocytes, also called red blood cells in the field of bio-rheology (cf. [4, 13–19] and the references therein for the bio-mathematical applications and other related information).

NIH Public Access

Author Manuscript

East Asian J Applied Math. Author manuscript; available in PMC 2013 June 19.

Published in final edited form as:

East Asian J Applied Math. 2011 January 1; 1(2): 155–171. doi:10.4208/eajam.030510.250910a.

NIH-PA Author Manuscript

(2)

The fluid equations can be formulated by either the Stokes equations (the inertial term is neglected)

(1.1) or the Navier-Stokes equations

(1.2)

with the fluid incompressibility constraint

(1.3)

Here, we assume that the interface motion is under a shear flow u = γ̇ y e1 along the boundary of a finite domain Ω as illustrated in Fig. 1, where γ̇ is the shear rate (a fixed number). The force term F(x, t) has the form

(1.4)

where X(s, t) is a parametric representation of the moving interface Γ, and {n, τ} are the unit normal and tangential directions of the moving interface, respectively. The bending force fb acting along the normal direction can be described as

(1.5)

where cb is a bending coefficient, and κ is the interface curvature. The detailed derivation of the above bending force can be found in [18]. However, unlike the previous literature in which a filter is used [3, 18], we add a tangential regularization force g(s, t) to the original problem such that we actually solve a modified perturbed problem. This technique is quite common for inverse problems.

The surface tension σ(s, t) on the interface is part of the unknown solutions and should be determined so that the inextensibility of the interface

(1.6)

is satisfied. Further, note that we can get

(1.7)

which means that the tension force affects both tangential and normal directions.

The difficulties in solving the above interfacial problem include (1) both the area and total length of the interface should be conserved simultaneously, (2) the problem is known to be very stiff, which may require small time steps for numerical methods, and (3) the problem may have no solution, as for example when the initial shape is a circle. In previous literature, most related work is based on boundary integral methods – cf. for example [13, 16, 18, 19]. However, boundary integral methods generally assume infinite domains, and cannot be generalized to full Navier-Stokes equations since there is no corresponding Green function.

NIH-PA Author Manuscript

(3)

Until recently, Kim and Lai [3] have applied a penalty immersed boundary method to simulate the dynamics of inextensible vesicles. By introducing two different kind of Lagrangian markers, the authors are able to decouple the fluid and vesicle dynamics so that the computation can be performed more efficiently. One potential difficulty of this approach is that the time step depends on the penalty number and must be chosen smaller as the penalty number becomes larger. There is no such dependency in the methods proposed in this paper, since no penalty numbers are involved.

Using the idea of the augmented immersed interface method, it is natural to set the surface tension σ(s, t) as an augmented variable. Once we know the surface tension, we know all the jump conditions related to the velocity and the pressure. The flow can be solved easily, as with our previously developed solver in [11]. The augmented equation is the surface divergence free condition (1.6) along the interface. We obtain usual discretized Stokes or Navier-Stokes equations with corrections near or on the interface, plus a much smaller system of equations for the augmented variable. This technique can be applied to full Navier-Stokes equations in both 2D and 3D. The difference is the size of the discrete system.

The rest of the paper is organized as follows. In the next section, we describe our method (three Poisson equations approach) for the model of Stokes equations. A regularization technique is introduced to control the magnitude of the artificial tangential force g, and the augmented approach is explained. In Section 3, we propose our modified projection method so that the pressure jump condition proportional to the unknown surface tension can be implemented. In Section 4, we show some numerical simulations and compare our results with others obtained in the literature. Some conclusions are made in the last section.

2. The Numerical Method for the Stokes Equations Model

We assume that the domain Ω is a rectangle [a, b] × [c, d]. The spatial increment is chosen as hx = (b−a)/M, hy = (d − c)/N, where M and N are the number of grid points in the x and y directions, respectively. We use a standard uniform Cartesian grid and the cubic spline package [8] to represent the moving interface. In this representation, the interface Γ is represented as a periodic cubic spline (X(s), Y(s)) in terms of the arc-length parameter s. We denote the number of control points for the cubic spline by Nb.

In the following, we introduce the three Poisson equations approach for solving the Stokes equations. By applying the divergence operator to the momentum equation (1.1) and reformulating the singular force by the jump conditions [5], we obtain the pressure equations

(2.1)

Notice that we use the zero Neumann boundary condition ∂p/∂n|_∂_Ω = 0 for simplicity, even though the Stokes equations do not impose the pressure boundary condition explicitly. Once the pressure is computed, we can apply the momentum equation (1.1) to obtain the velocity by solving the equations

(2.2)

NIH-PA Author Manuscript

(4)

See for example [5], for a discussion of the jump conditions; and we note that the pressure equation (2.1) and the velocity equation (2.2) are coupled via the unknown surface tension σ(s, t), which should be determined by the inextensibility constraint (1.6). Since the system of equations consisting of (2.1)–(2.2) and (1.6) is complete, we can discretize those equations directly by the immersed interface method as follows.

Before proceeding, let be the column vector whose components are Uij and Pij, the approximate velocity and the pressure on a Cartesian grid at one particular time step. Let Q be the vector of the discrete unknown surface tension σ at some control points on the interface. The IIM discretization gives the matrix-vector equation

(2.3)

for some vector F1 and sparse matrices A11, A12, and A13. It requires solving three Poisson equations with different jump conditions to get . Note that the cost of solving A11 = F1 is about the cost of calling a fast Poisson solver three times [1, 6].

Along the interface, the inextensible condition (1.6) is discretized by a least squares interpolation scheme

(2.4)

where A21, A21, and A23 are some sparse matrices, and F2 is a vector (cf. [2] for more details).

When G = 0 the system of equations (2.3)–(2.4) is complete, since the number of equations and the number of unknowns are the same with A13 = 0 and A23 = 0. The system has a unique solution for both and Q. However, the condition number of the Schur complement for Q is quite large, reflecting that the problem is nearly ill-posed. The error would

accumulate and affects the accuracy of the computed area eventually, due to the fact that the inextensible condition is enforced at the interface and the incompressible condition is enforced at the grid points simultaneously. Our regularization strategy is to introduce a tangential force g, or G in the discrete case, which will not alter the interface motion much – and at the same time we can control the magnitude of G to be as small as possible.

Assume that the dimension of the vector Q is Nb, so the dimension of G is also Nb, and we have an additional Nb degrees of the freedom. Since the area enclosed by the interface Γ must be conserved, we enforce the condition

(2.5)

in the discrete case. Again, the velocity on the interface is interpolated from a least squares interpolation scheme using the velocity at the grid points, hence the discrete analogue of above equation (2.5) can be written as

(2.6) where A31,A32,A33 are row vectors and F3 is a number.

The Schur complement for [ Q G ]T is

NIH-PA Author Manuscript

(5)

(2.7)

where

(2.8)

Note that [ Q G ]T has dimension of 2Nb, so we still have Nb − 1 degrees of the freedom. Meanwhile, we do want to control the magnitude of the regularized tangential force term G, so we put

(2.9)

where E is a regularization parameter, which should be chosen large enough. However, by adding the above equation (2.9), the system of equations (2.3), (2.4), and (2.9) is now over-determined, for the total number of equations is one more than the total number of

unknowns. Thus we seek the least squares solution of (2.7) and (2.9), which is equivalent to solving the following normal equations (Tickhonov regularization):

(2.10)

We can choose the regularization parameter E to control the magnitude of G. In fact, the larger E is, the smaller is ‖G‖.

In the linear system of equations (2.3) and (2.7), we solve (2.7) first since it is O(Nb) dimensional lower than that of . Once we have Q and G, we can get by applying the one-step Stokes solver once. There are other components of the augmented method, such as how to obtain the matrix-vector multiplication for the Schur complement – cf. [2, 9].

3. The Numerical Method for the Navier-Stokes Equation Model

Numerical simulations of an inextensible interface in an incompressible flow have been reported in [15, 16, 19], using the boundary integral method for Stokes equations in which the Green function is available. Another novelty of this paper is the development of the finite difference method for the Navier-Stokes equation model that may be more realistic, since the inertial effect is taken into account. In addition to the difficulties that we mentioned for the Stokes model, another difficulty is how to enforce the pressure jump condition in the projection method. This turns out to be a crucial step in our numerical method. Here, we propose a modified projection method to enforce the pressure jump condition. Similar to the three Poisson equations approach for Stokes equations, we apply the divergence operator to the momentum equation (1.2) and reformulate the singular force term by the jump conditions to obtain

(3.1)

(3.2)

Thus, our modified projection method from time tk_{to t}k+1_{can be written as}

NIH-PA Author Manuscript

(6)

(3.3)

(3.4)

(3.5)

(3.6)

(3.7) In addition, the inextensible condition (1.6) expressed as

and the area conserving condition (2.5) expressed as

should be satisfied as well. The boundary condition for u* is the shear flow condition on the rectangular domain. Once we have the velocity uk+1, we can move the control points of the interface to new positions and form a new interface by our cubic spline representation. The remaining discretization details are as in Section 2, except that A₁₁ now corresponds to the Navier-Stokes instead of the Stokes solver. It is worth mentioning that in our simulations with the Navier-Stokes model we often obtain a better computed divergence-free velocity, compared with that obtained from the Stokes equations, since the projection step enforces the divergence free (or incompressible) condition at grid points.

4. Numerical Results

In this section, we present some simulations of the motion of an inextensible interface in a shear flow. All the computations were performed at the North Carolina State University using either notebook or desktop computers. We use the same notations as described at the beginning of Section 2. We use the cubic spline package [8] to represent the moving interface. In this representation, the interface is represented as a periodic cubic spline (X(s), Y(s)) in terms of the arc-length parameter s, and Nb denotes the number of control points for the cubic spline. In most simulations, the viscosity coefficient chosen is μ = 20.

NIH-PA Author Manuscript

(7)

4.1. Accuracy Check Against an Exact Solution

We first check our numerical algorithm and test the accuracy against an exact solution with a stationary interface r = 0.5 in the domain [−1, 1] × [−1, 1]. The Dirichlet boundary and initial conditions are from the following exact solution:

(4.1)

(4.2)

(4.3)

where and C is a constant. The force term F is derived directly from the exact solution.

In Table 1, we show a grid refinement analysis to check the order of accuracy for our method. We set the following

as the errors of the velocity and the pressure, respectively, at time T = 0.5. The number order is the approximated order of accuracy from the two consecutive errors for the solution. Second order accuracy is clearly seen for the velocity. The pressure seems to be second order accurate because ∂p/∂n = 0 is satisfied exactly for the solution.

4.2. The Computed Surface Tension and the Tank-treading Motion

Next, we show a grid refinement analysis for the computed surface tension of an

inextensible interface when it reaches the quasi-equilibrium state. In the quasi-equilibrium state, there is only a tangential motion that behaves like tank-treading. Fig. 2 shows the computed surface tension with different mesh sizes when the bending coefficient is zero, so that we can compare with the results obtained in [16, 19]. The initial interface is an ellipse with major and minor axes being 0.5 and 0.3, respectively. In Fig. 2, we have observed that there are some sawtooth oscillations in the computed surface tension. This oscillatory behavior has also been found in the previous study [19], where a point-wise collocation method is applied. In order to remove the oscillations, one might apply numerical smoothing so that the high frequency modes can be eliminated. Here, however, these sawtooth

oscillations generally have low amplitudes as we refine the mesh, and the surface tension remains with sufficient accuracy as long as the interface is smooth enough. The results also agree well with those obtained in [16].

In Fig. 3 (a), we plot an interface (x/0.3)2_{+ (y/0.5)}2_{= 1 at initial time (t = 0), and at time t =} 6.1529 when the interface has reached its quasi-equilibrium state. The initial length of the interface and the area enclosed by the interface can easily be calculated as (2.5523, 0.47124). The regularization parameter is chosen as E = 103_{and the shear rate is γ̇ = 0.5.} We also carry out the grid refinement analysis for the computed length and the area at t = 2.

NIH-PA Author Manuscript

(8)

The errors with M = N = 80, Nb = 80 are (−8.7566 × 10−4, 7.4570 × 10−4); and (−2.1310 × 10−4, 6.4505 × 10−4) when M = N = 160, Nb = 160 is used. The results show second order accuracy in both length and area for this case. In fact, the errors in the length and area are within the discrete error of the cubic spline interpolation in discretizing the interface. In Fig. 3 (b), we plot the velocity field near some portion of the interface at t = 6.1529 after the interface reaches its quasi-equilibrium state. We can see that the direction of the velocity is along the tangential direction of the interface, indicating that the interface undergoes a tank-treading motion.

Note that when the initial shape is a circle, the coefficient matrix for the discrete surface tension is singular. If there is no flow (i.e. the shear rate is zero), then there are infinite solutions of the unknown surface tension. If the shear rate is non-zero, then no such surface tension exists. Mathematically, we know that that a circle has largest area among all geometric objects with fixed circumference. This was confirmed in our numerical tests.

4.3. Effect of the Cell Circularity

It is stated in [16, 19] that the final shape and the inclination angle (angle with x axis) of an inextensible interface at the quasi-equilibrium depend on the circularity defined by

(4.4)

where L and A are the length and area of the interface, respectively. Note that c = 1 for a circle, and for smaller circularity the inclination angle at the quasi-equilibrium is larger. If the circularity is large enough, then the interface at the quasi-equilibrium state has a biconcave shape.

In Fig. 4, we show the evolution of an inextensible interface with relatively large circularity, c = 1.7770. The initial interface is an ellipse (x/0.1)2_{+ (y/0.55)}2_{= 1. In Fig. 4 (a), the results} are obtained using the Stokes model; while in Fig. 4 (b), the results are obtained using the Navier-Stokes model. Compared with Fig. 3 (a), Fig. 5 (c) and (d), and Fig. 6, the quasi-equilibrium interface with a larger circularity is more likely to have a biconcave shape and a smaller inclination angle. Note that, due to the large curvature near the tip area, an extensive smoothing may cause the area loss.

4.4. The Stokes Equations Model Versus the Navier-Stokes Model and the Bending Effect

When the angle of the initial interface (an ellipse) with the x-axis is less than 90 degrees, the results for the interface quasi-equilibrium are almost the same for both the Stokes and the Navier-Stokes models, and smooth in either case. However, if the angle of the initial interface with the x-axis is larger than 90 degrees, we can see that some oscillations along the interface develop in the beginning and disappear after the angle becomes less than 90 degrees – cf. Fig. 5. We also see stronger oscillations for the Navier-Stokes equations with μ = 20, compared with that for the Stokes equations, due to the convection effect of the Navier-Stokes equations. In Fig. 5 (a) and (b), the initial interface is an ellipse x2_/0.12₊ y2_/0.552_{= 1 with major axis rotated by 45 degrees counter-clockwise. The circularity is c =} 1.7770. The computations are performed using a 160 × 160 grid with Nb = 160. Both the area and length are almost constants after the inclination angle becomes less than 90 degrees, and the errors in the length and the area are about 5%. We observe the biconcave shape when the interface reaches its quasi-equilibrium state. In Fig. 5 (c) and (d), the initial interface is an ellipse x2_/0.32_{+ y}2_/0.52_{= 1 with major axis rotated by 45 degrees} counter-clockwise, whose circularity is c = 1.0491 (smaller than the previous one). Again, we

NIH-PA Author Manuscript

(9)

observe some oscillations for the Navier-Stokes equations model compared with almost no oscillations for the Stokes model.

We also test the effect of the bending force. In Fig. 6, we show two different bending coefficients cb = 10−6 and 10−2 used in the Navier-Stokes model. We have observed that the stronger the bending coefficient the smoother the interface, during the transition process.

4.5. More General Initial Shapes

We show some simulations with more general initial configurations using the Navier-Stokes model. The results are comparable with those presented in [16, 19] using the Stokes

equations. In Fig. 7, we show the evolution process of an inextensible interface with the initial configuration

(4.5)

where α = 0.2451 and β = 0.1473. The computational domain is [−1, 1] × [−1, 1]. In the next test, we show the results with an initial configuration from the equation

(4.6)

where α = 0.207161, β = 2.002558, and γ = −1.122762, see [19]. The initial configuration is biconcave like a dumbbell. As in [19], we put the interface in two different orientations; one is in the flow direction (Fig. 8); and another is inclined parallel with the x-axis (Fig. 9). In the first case, the interface stretches gradually until it reaches its quasi-equilibrium state. Again, we see the biconcave shape due to the large circularity. In the second case, however, the interface rotates and stretches simultaneously. It rotates immediately so that the interface is aligned in the direction of the flow before gradually relaxing to the quasi-equilibrium state.

5. Conclusions

In this paper, we propose new finite difference methods for simulating the motion of an inextensible interface in an incompressible flow using either the Stokes or the Navier-Stokes equations. The unknown surface tension should be determined in such a way that the surface divergence of the velocity must be zero along the interface. For the Stokes equations model, we use a three Poisson approach. For the Navier-Stokes model, we propose a modified projection method, which can enforce the pressure jump condition that corresponds directly to the unknown surface tension. A regularization technique is useful for the interfaces with large curvature, or when the orientation of the interface is at a large angle with the flow direction.

When the circularity and the Reynolds number are modest and the initial angle between the interface and the flow direction is less than 90 degrees, our numerical methods can preserve the area enclosed by the interface and the total length of the interface very well, even if no regularization is used, (i.e. g = 0). We have analyzed the effects of the circularity, flow directions, and different initial configurations on the motion of the inextensible interface.

Acknowledgments

We thank Drs. Peng Song (UCI), Pingwen Zhang (Peking University), and Kazufumi Ito (NCSU) for valuable discussions.

NIH-PA Author Manuscript

(10)

The authors also acknowledge the following sources of financial support. The first author is partially supported by the US ARO grants 56349MA-MA, and 550694-MA, the AFSOR grant FA9550-09-1-0520, the US NSF grant DMS-0911434, the NIH grant 096195-01, and the C-NSF grant 11071123. The second author is partially supported by the Taiwan National Science Council under grant NSC-97-2628-M-009-007-MY3 and MoE-ATU project.

References

1. Adams, J.; Swarztrauber, P.; Sweet, R. Fishpack: Efficient Fortran subprograms for the solution of separable elliptic partial differential equations. http://www.netlib.org/fishpack/.

2. Ito K, Li Z, Lai M-C. An augmented method for the Navier-Stokes equations on irregular domains. J. Comput. Phys. 2009; 228:2616–2628.

3. Kim Y, Lai M-C. Simulating the dynamics of inextensible vesicles by the penalty immersed boundary method. J. Comput. Phys. 2010; 229:4840–4853.

4. Kraus M, Wintz W, Seifert U, Lipowsky R. Fluid vesicles in shear flow. Physical Review Letters. 2005; 95

5. Lai M-C, Li Z. A remark on jump conditions for the three-dimensional Navier-Stokes equations involving an immersed moving membrane. Applied Mathematics Letters. 2001; 14, vol 2:149–154. 6. LeVeque RJ, Li Z. Immersed interface method for Stokes flow with elastic boundaries or surface

tension. SIAM J. Sci. Comput. 1997; 18:709–735.

7. Li Z, Lai M-C, He G, Zhao H. An augmented method for free boundary problems with moving contact lines. Computers and Fluids. 2010; 39:1033–1040.

8. Li, Z. IIMPACK: A collection of fortran codes for interface problems. Anonymous ftp at

ftp.ncsu.edu under the directory: /pub/math/zhilin/Package and http://www4.ncsu.edu/~zhilin/IIM, last updated: 2008.

9. Li, Z.; Ito, K. The Immersed Interface Method – Numerical Solutions of PDEs Involving Interfaces and Irregular Domains. SIAM Frontier Series in Applied mathematics, FR33. 2006.

10. Li Z, Ito K, Lai M-C. An augmented approach for Stokes equations with a discontinuous viscosity and singular forces. Computers and Fluids. 2007; 36:622–635.

11. Li Z, Lai M-C. The immersed interface method for the Navier-Stokes equations with singular forces. J. Comput. Phys. 2001; 171:822–842.

12. Li Z, Wan X, Ito K, Lubkin S. An augmented pressure boundary condition for a Stokes flow with a non-slip boundary condition. Commun. Comput. Phys. 2006; 1:874–885.

13. Lowengrub JS, Ratz A, Voigt A. Phase-field modeling of the dynamics of multicomponent vesicles: Spinodal decomposition, coarsening, budding, and fissio. Physical Review E. 2009; 79 14. Markvoort AJ, Van Santen RA, Hilber PAJ. Vesicle shapes from molecular dynamics simulation.

J. Physical Chemistry B. 2006; 110

15. Sohn JS, Tseng Y-H, Li S, Voigt A, Lowengrub JS. Dynamics of multicomponent vesicles in a viscous fluid. J. Comput. Phys. 2010; 229:119–144. [PubMed: 20808718]

16. Song, P. PhD thesis. Peking University; 2008. Numerical simulations of clindrical membranes. 17. Veerapaneni SK, Gueyffier D, Zorin D, Biro G. A numerical method for simulating the dynamics

of 3d axismmetric vesicles suspended in viscous flow. J. Comput. Phys. 2009; 228:7233–7249. 18. Veerapaneni SK, Gueyffier D, Zorin D, Biros G. A boundary integral method for simulating the

dynamics of inextensible vesicles suspended in a viscous fluid in 2d. J. Comput. Phys. 2009; 228:2334–2353.

19. Zhou H, Pozrikidis C. Deformation of capsules with incompressible interfaces in simple shear flow. J. Fluid Mech. 1995; 283:175–200.

NIH-PA Author Manuscript

(11)

Figure 1.

A diagram of an inextensible interface in a shear flow.

NIH-PA Author Manuscript

(12)

Figure 2.

The computed surface tension with different mesh sizes and the number of control points. We have observed that there are some sawtooth oscillations for the surface tension when the interface is in its quasi-equilibrium state. However, by refining the mesh, the magnitude of the oscillations can be reduced significantly. The results agree with the those in [19].

NIH-PA Author Manuscript

(13)

Figure 3.

(a) An initial interface (red-curve) and its quasi-equilibrium state (blue-curve) at t = 6.1529. The errors of the total length and the area enclosed by the interface are 1.1067 × 10−3 and 5.0365 × 10−3, respectively. (b) Plot of the velocity field near a portion of the interface after it reaches the quasi-equilibrium state. The velocity is along the tangential direction.

NIH-PA Author Manuscript

(14)

Figure 4.

Evolution of an inextensible interface with relatively large circularity, c = 1.7770. (a) The results obtained using the Stocks equations. (b) The results obtained using the Navier-Stokes equations. Compared with Fig. 3, Fig. 5 (c), (d), and Fig. 6, the interface at the quasi-equilibrium state has a biconcave shape with a smaller inclination angle.

NIH-PA Author Manuscript

(15)

Figure 5.

The snap shots of an intextensible interface in a shear flow. (a) and (c): simulation results using the Stokes equations; (b) and (d): simulation results using the Navier-Stokes equations.

NIH-PA Author Manuscript

(16)

Figure 6.

Effect of the bending force – viz. the stronger the bending coefficient, the smoother the interface. The snapshots in the top row are the simulation results with cb = 10−6 in which we see some oscillations along the interface. The snapshots in the bottom row are the simulation results with larger bending coefficient cb = 10−2 in which the interface remains smooth during the evolution process. The red-curves in the very right column are the plots of the initial interface (t = 0).

NIH-PA Author Manuscript

(17)

Figure 7.

A rose shaped initial interface deforms under the shear flow. (a) The initial and final configurations (t = 18.494) of the interface. The errors of the length and the area are 5.6375 × 10−2 and 4.5774 × 10−2 respectively. (b) The snapshots at different times.

NIH-PA Author Manuscript

(18)

Figure 8.

The evolution of a biconcave interface in the orientation of the flow direction. (a) The initial state is at t = 0 and the quasi-equilibrium state is at t = 2.4986. (b) Some snapshots at different times. The motion is comparable with the one obtained in [19].

NIH-PA Author Manuscript

(19)

Figure 9.

The evolution of a biconcave initial interface in the horizontal (x) direction. (a) The initial state is at t = 0 and the quasi-equilibrium state is at t = 2.4914. (b) Some snapshots at different times. The interface starts with a rigid body rotation so that it can be aligned with the flow direction, an then deforms gradually to the quasi-equilibrium state.

NIH-PA Author Manuscript

(20)

NIH-PA Author Manuscript

Table 1

A grid refinement analysis against the exact solution at time T = 0.5. ‖ Eu ‖∞ is the maximal error in the velocity, and "order" is the approximated order of accuracy computed from the two consecutive errors.

M = N, Nb ‖ Eu‖∞ orderu ‖ Ep‖∞ orderp 32, 32 1.2850 10−3 _{1.0643 10}−2 64, 64 3.1483 10−3 _2.0291 _{2.5466 10}−3 _2.0633 128, 128 6.7513 10−4 _2.2213 _{9.2588 10}−4 _1.4596 256, 256 1.6127 10−5 _2.0657 _{2.7097 10}−4 _1.7727 512, 512 4.1006 10−6 _1.9756 _{7.0427 10}−5 _1.9440