doi: 10.4208/eajam xxx 201x
New Finite Difference Methods Based on IIM for Inextensible Interfaces in Incompressible Flows
Zhilin Li∗,1and Ming-Chih Lai2
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.
Received 3 May 2010; Accepted (in revised version) 25 September 2010 Available online xx XXXXX 2010
Abstract. In this paper, new finite difference methods based on the augmented im- mersed 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 govern- ing 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 nonlin- ear 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 perturba- tions, 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. Sev- eral numerical experiments show good agreement with others results in the literature and reveal some interesting phenomena.
AMS subject classifications: 65M06, 65M12, 76T05
Key words: Inextensible interface, incompressible flow, Stokes equations, Navier-Stokes equa- tions, immersed interface method, inverse problem, regularization, augmented immersed interface method.
∗Corresponding author.
http://www.global-sci.org/eajam 1 "201x Global-Science Pressc
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 un- known 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).
Ω Γ
The fluid equations can be formulated by either the Stokes equations (the inertial term is neglected)
∇p = µ ∆u + F(x, t), x∈ Ω, (1.1)
or the Navier-Stokes equations
∂ u
∂ t + u· ∇u + ∇p = µ ∆u + F(x, t), x∈ Ω, (1.2) with the fluid incompressibility constraint
∇ · u = 0. (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 termF(x, t) has the form
F(x, t) =
!
Γ(t)
"
∂
∂ s
#
σ(s, t) τ(s, t)
$
+ fbn + g(s, t)τ(s, t)
%
δ(x− X(s, t)) ds, (1.4)
whereX(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
fb=−cb
&
κss+κ3 2
'
, (1.5)
where cbis 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
(∇s· u)Γ= ∂ u
∂ τ· τ
*
*
*
*Γ
= 0 (1.6)
is satisfied. Further, note that we can get
∂
∂ sσ(s, t) τ(s, t) = ∂
∂ sσ(s, t)τ + σκn, (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 of 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. 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 regular- ization 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 pro- jection 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 com- pare 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,
∆p = 0, x∈ Ω \ Γ, ∂ p
∂ n
*
*
*
*∂ Ω
= 0,
[p]*
*Γ= σκ + fb,
+∂ p
∂ n ,
Γ
= ∂ g
∂ s. (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
µ∆u =∇p, [u]|Γ= 0, +
µ∂ u
∂ n ,
Γ
=− -∂ σ
∂ s + g .
τ. (2.2)
See for example [5], for a discussion of the jump conditions; and we note that the pres- sure 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 areUi j and Pi j, 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
A11! + A12Q + A13G = F1, (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! = F1is 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
A21! + A22Q + A23G = F2, (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 Nbdegrees of the freedom. Since the area enclosed by the interface Γ must be conserved, we enforce the condition
!
Γ
(u· n ) ds = 0, or
Nb
/
k=1
(Uk· nk) ∆sk= 0, (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
A31! + A32Q + A33G = F3, (2.6)
where A31, A32, A33are row vectors and F3 is a number.
The Schur complement for [ Q G ]T is
B
Q G
=
F2 F3
−
A21 A31
A−111F1=
F¯2 F¯3
, (2.7)
where
B =
4 A22 A23 A32 A33
5
−
A21 A21
A−1116
A12 A13 7
. (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
E G = 0 (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 number of total equations is one more than the number of total 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):
8BTB + E2I9
Q G
= BT
F¯2 F¯3
. (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 realis- tic, 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
∆p =− ∇ · (u · ∇u) , ∂ p
∂ n
*
*
*
*∂ Ω
= 0, (3.1)
[p]*
*Γ= σκ + fb,
+∂ p
∂ n ,
Γ
=∂ g
∂ s. (3.2)
Thus, our modified projection method from time tkto tk+1can be written as
∆˜pk+1=− ∇ · (u · ∇u)k, ∂ ˜pk+1
∂ n
*
*
*
*
*∂ Ω
= 0,
6˜pk+17
Γk = σk+1κ + fbk 4 ∂ ˜pk+1
∂ n 5
Γk
= ∂ gk+1
∂ s ,
(3.3)
u∗− uk
∆t + (u· ∇ u)k+∇ ˜pk+1= µ∆ u∗, [u∗]Γk = 0,
4 µ∂ u∗
∂ n 5
Γk
=−& ∂ σk+1
∂ s + gk+1(s) '
τ,
(3.4)
∆ φk+1= ∇ · u∗
∆t ,
∂ φk+1
∂ n
*
*
*
*
*∂ Ω
= 0, 6
φk+17
Γk= 0, 4 ∂ φk+1
∂ n 5
Γk
= 0,
(3.5)
uk+1= u∗− ∆t ∇ φk+1, (3.6)
pk+1= ˜pk+1+ φk+1. (3.7)
In addition, the inextensible condition (1.6) expressed as 8∇s· uk+19
Γk= ∂ uk+1
∂ τ · τ
*
*
*
*
*Γk
= 0, and the area conserving condition (2.5) expressed as
!
Γk
8uk+1· n9
ds = 0,
should be satisfied as well. The boundary condition foru∗ is the shear flow condition on the rectangular domain. Once we have the velocityuk+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 A11 now corre- sponds 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
T = 0.5 ' Eu'∞ or d er
M = N , Nb ' Eu'∞ or d eru ' Ep'∞ or d erp 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
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 Nbdenotes the number of control points for the cubic spline. In most simulations, the viscosity is chosen asµ = 20.
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:
u(x , y, t) =
sin(t)
#y
r − 2y$
, if r≥ 1/2, sin(t)
- r2− 1
4 .
y, otherwise,
(4.1)
v(x , y, t) =
sin(t)
#
−x r + 2x
$
, if r≥ 1/2,
− sin(t) -
r2−1 4
.
x , otherwise,
(4.2)
p(x , y, t) = C sin(t), (4.3)
where r =>
x2+ y2 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
' Eu'∞= max
i j {|Ui jk− u(xi, yj, T )|} + max
i j {|Vi jk− v(xi, yj, T )|}, ' Ep'∞= max
i j {|Pi jk− p(xi, yj, T )|},
(a) 80× 80
0 10 20 30 40 50 60 70 80
−2.7
−2.6
−2.5
−2.4
−2.3
−2.2
−2.1
−2
−1.9
−1.8
M=N=N1=80, E =1e+3
(b) 160× 160
0 20 40 60 80 100 120 140 160
−1.85
−1.8
−1.75
−1.7
−1.65
−1.6
−1.55
−1.5
−1.45
−1.4
M=N=N1=160, E =1e+3
(c) 320× 80
0 10 20 30 40 50 60 70 80
−2
−1.9
−1.8
−1.7
−1.6
−1.5
−1.4
−1.3
M=N=320,Nl=80, E =1e+3
(d) 320× 320
0 50 100 150 200 250 300 350
−1.95
−1.9
−1.85
−1.8
−1.75
−1.7
−1.65
−1.6
−1.55
−1.5
−1.45
M=N=N1=320, E =1e+3
as the errors of the velocity and the pressure, respectively, at time T = 0.5. The number or d er is the approximated order of accuracy from the two consecutive errors for the so- lution. 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 in- extensible 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
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8 1
t=0
6.1529
err=[1.1067e−03, 5.0365e−03]
M=128, Nb=120, E=103
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4
(a) (b)
t = 6.1529 1.1067×10−3 5.0365×10−3
method is applied. In order to remove the oscillations, one might apply numerical smooth- ing so that the high frequency modes can be eliminated. Here, however, these sawtooth oscillations are generally 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. 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
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8 1
t=0
t=6.0217
err=[1.0750e−02, 3.9102e−03]
M=160, Nb=120, E=103
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8 1
t=0
t=6.1341
err = [1.2254e−4, 1.7835e−3]
M=128, Nb=120, E=103
(a) (b)
c = 1.7770
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
c = L 2*
Aπ, (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 circu- larity, 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 de- grees, 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 oscilla- tions 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 convec- tion 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+ y2/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 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
? X (θ ) = α cos(θ ) + β cos(3θ ) cos(θ ),
Y (θ ) = α sin(θ ) + β cos(3θ ) sin(θ ), 0≤ θ < 2π, (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 Y (x ) =±1
2
>
1− x2#
α + β x2+ γx4
$
, (4.6)
where α = 0.207161, β = 2.002558, and γ =−1.122762, see [19]. The initial config- uration 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
−1 0 1
−1
−0.5 0 0.5 1
t = 0
−1 0 1
−1
−0.5 0 0.5 1
t=12.3048
−1 0 1
−1
−0.5 0 0.5 1
t=1.3708
−1 0 1
−1
−0.5 0 0.5 1
t=2.7416
−1 0 1
−1
−0.5 0 0.5 1
t=4.4550
−1 0 1
−1
−0.5 0 0.5 1
t=6.1685
−1 0 1
−1
−0.5 0 0.5 1
t=10.3816
−1 0 1
−1
−0.5 0 0.5 1
t=6.9119
−1 0 1
−1
−0.5 0 0.5 1
t = 0
−1 0 1
−1
−0.5 0 0.5 1
t=1.3494
−1 0 1
−1
−0.5 0 0.5 1
t=2.6987
−1 0 1
−1
−0.5 0 0.5 1
t=4.4513
(a) (b)
−1 −0.5 0 0.5 1
−1
−0.5 0 0.5 1
−1 −0.5 0 0.5 1
−1
−0.5 0 0.5 1
t = 0
−1 −0.5 0 0.5 1
−1
−0.5 0 0.5 1
t=1.0281
−1 −0.5 0 0.5 1
−1
−0.5 0 0.5 1
t=2.0562
t=6.1671 t=12.334
−1 −0.5 0 0.5 1
−1
−0.5 0 0.5
1 t=5.9958
t=12.162
−1 −0.5 0 0.5 1
−1
−0.5 0 0.5 1
t=1.0281
−1 −0.5 0 0.5 1
−1
−0.5 0 0.5 1
t = 0
−1 −0.5 0 0.5 1
−1
−0.5 0 0.5 1
t=2.0562
(c) (d)
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
−1 0 1
−1
−0.5 0 0.5 1
t=1.0364
−1 0 1
−1
−0.5 0 0.5 1
t=1.3664
−1 0 1
−1
−0.5 0 0.5 1
t=6.1467
−1 0 1
−1
−0.5 0 0.5 1
−1 0 1
−1
−0.5 0 0.5 1
−1 0 1
−1
−0.5 0 0.5 1
cb= 10−6
cb= 10−2
t = 0
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8 1
err=[5.6375e−02, 4.5774e−02]
t=0
t=18.494 M=128, Nb=120, E=103
−1 0 1
−1
−0.5 0 0.5 1
−1 0 1
−1
−0.5 0 0.5 1
−1 0 1
−1
−0.5 0 0.5 1
−1 0 1
−1
−0.5 0 0.5 1
−1 0 1
−1
−0.5 0 0.5 1
−1 0 1
−1
−0.5 0 0.5 1
−1 0 1
−1
−0.5 0 0.5 1
−1 0 1
−1
−0.5 0 0.5 1
−1 0 1
−1
−0.5 0 0.5 1
(a) (b)
t = 18.494 5.6375× 10−2
4.5774× 10−2
equations model, we use a three Poisson approach. For the Navier-Stokes model, we pro- pose 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.
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2
−2
−1.5
−1
−0.5 0 0.5 1 1.5 2
error=[−1.2656e−04, e6.4850−03], t=2.4986 M=128, Nb=120, E=103, µ=20
−2 −1 0 1 2
−2
−1 0 1 2
t=0
−2 −1 0 1 2
−2
−1 0 1 2
t=0.5556
−2 −1 0 1 2
−2
−1 0 1 2
t=1.3889
−2 −1 0 1 2
−2
−1 0 1 2
t=2.4986
(a) (b)
t = 0 t = 2.4986
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2
−2
−1.5
−1
−0.5 0 0.5 1 1.5 2
err=[−1.2760e−04,6.4350e−03], t=2.4914 M=128, N
b=120, E=103, µ=20
−2 −1 0 1 2
−2
−1 0 1 2
t=0
−2 −1 0 1 2
−2
−1 0 1 2
t=2.7778e−02
−2 −1 0 1 2
−2
−1 0 1 2
t=5.5556e−01
−2 −1 0 1 2
−2
−1 0 1 2
t=2.4986
(a) (b)
x
t = 0 t = 2.4914
When the circularity and the Reynolds number are modest and the initial angle be- tween 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 inex- tensible interface.
Acknowledgments
We thanks Drs. Peng Song (UCI), Pingwen Zhang (Peking University), and Kazufumi Ito (NCSU) for valuable discussions.
The authors also acknowledge the following sources of financial supports. The first au- thor is partially supported by the US ARO grants 56349MA-MA, and 550694-MA, the AF- SOR 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 National Sci- ence Council of Taiwan under grant NSC-97-2628-M-009-007-MY3 and MoE-ATU project.
References
[1] J. Adams, P. Swarztrauber, and R. Sweet. Fishpack: Efficient Fortran subprograms for the solution of separable elliptic partial differential equations. http://www.netlib.org/fishpack/.
[2] K. Ito, Z. Li, and M-C. Lai. An augmented method for the Navier-Stokes equations on irregular domains. J. Comput. Phys., 228:2616–2628, 2009.
[3] Y. Kim and M.-C. Lai, Simulating the dynamics of inextensible vesicles by the penalty im- mersed boundary method, J. Comput. Phys., 229:4840-4853, 2010.
[4] M. Kraus, W. Wintz, U. Seifert, and R. Lipowsky. Fluid vesicles in shear flow. Physical Review Letters, 95, 2005.
[5] M.-C. Lai and Z. Li, A remark on jump conditions for the three-dimensional Navier-Stokes equations involving an immersed moving membrane, Applied Mathematics Letters, 14, vol 2 , 149-154, (2001).
[6] R. J. LeVeque and Z. Li. Immersed interface method for Stokes flow with elastic boundaries or surface tension. SIAM J. Sci. Comput., 18:709–735, 1997.
[7] Z. Li, , M-C. Lai, G. He, and H. Zhao. An augmented method for free boundary problems with moving contact lines. Computers and Fluids, 39:1033–1040, 2010.
[8] Z. Li. 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] Z. Li and K. Ito. The Immersed Interface Method – Numerical Solutions of PDEs Involving Interfaces and Irregular Domains. SIAM Frontier Series in Applied mathematics, FR33, 2006.
[10] Z. Li, K. Ito, and M-C. Lai. An augmented approach for Stokes equations with a discontinuous viscosity and singular forces. Computers and Fluids, 36:622–635, 2007.
[11] Z. Li and M-C. Lai. The immersed interface method for the Navier-Stokes equations with singular forces. J. Comput. Phys., 171:822–842, 2001.
[12] Z. Li, X. Wan, K. Ito, and S. Lubkin. An augmented pressure boundary condition for a Stokes flow with a non-slip boundary condition. Commun. Comput. Phys., 1:874–885, 2006.
[13] J. S. Lowengrub, A. Ratz, and A. Voigt. Phase-field modeling of the dynamics of multicom- ponent vesicles: Spinodal decomposition, coarsening, budding, and fissio. Physical Review E, 79, 2009.
[14] A. J. Markvoort, R.A. Van Santen, and P.A.J. Hilber. Vesicle shapes from molecular dynamics simulation. J. Physical Chemistry B, 110, 2006.
[15] J. S. Sohn, Y.-H. Tseng, S. Li, A. Voigt, and J. S. Lowengrub. Dynamics of multicomponent vesicles in a viscous fluid. in press, 2010.
[16] P. Song. Numerical simulations of clindrical membranes. PhD thesis, Peking University, 2008.
[17] S. K. Veerapaneni, D. Gueyffier, D. Zorin, and G. Biro. A numerical method for simulating the dynamics of 3d axismmetric vesicles suspended in viscous flow. J. Comput. Phys., 228:7233- 7249, 2009.
[18] 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:2334-2353, 2009.
[19] H. Zhou and C. Pozrikidis. Deformation of capsules with incompressible interfaces in simple shear flow. J. Fluid Mech., 283:175–200, 1995.
