The Immersed Interface Method for the Navier–Stokes Equations
with Singular Forces
Zhilin Li and Ming-Chih Lai
Center for Research in Scientific Computation & Department of Mathematics, North Carolina State University, Raleigh, North Carolina 27695-8205; and Department of Mathematics,
Chung Cheng University, Minghsiung, Chiayi 621, Taiwan E-mail: zhilin@math.nscu.edu; mclai@math.ccu.edu.tw
Received September 12, 2000; revised April 5, 2001
Peskin’s Immersed Boundary Method has been widely used for simulating many fluid mechanics and biology problems. One of the essential components of the method is the usage of certain discrete delta functions to deal with singular forces along one or several interfaces in the fluid domain. However, the Immersed Boundary Method is known to be first-order accurate and usually smears out the solutions. In this paper, we propose an immersed interface method for the incompressible Navier–Stokes equations with singular forces along one or several interfaces in the solution domain.
The new method is based on a second-order projection method with modifications only at grid points near or on the interface. From the derivation of the new method, we expect fully second-order accuracy for the velocity and nearly second-order accuracy for the pressure in the maximum norm including those grid points near or on the interface. This has been confirmed in our numerical experiments. Furthermore, the computed solutions are sharp across the interface. Nontrivial numerical results are provided and compared with the Immersed Boundary Method. Meanwhile, a new version of the Immersed Boundary Method using the level set representation of the interface is also proposed in this paper. ° 2001 Academic Pressc
Key Words: Navier–Stokes equations; interface; discontinuous and nonsmooth so- lution; immersed interface method; immersed boundary method; projection method;
level set method.
1. INTRODUCTION
In this paper, we consider the incompressible Navier–Stokes equations in a bounded domainÄ ⊂ R2that contains one or several interfaces0(t):
ρ µ∂u
∂t + (u · ∇)u
¶
+ ∇ p = µ1u + G + F, x ∈ Ä, (1.1)
822 0021-9991/01 $35.00
Copyright c° 2001 by Academic Press All rights of reproduction in any form reserved.
∇ · u = 0, (1.2)
u|∂Ä= ub, BC, (1.3)
u(x, 0) = u0, IC. (1.4)
Here we write F and G separately to distinguish different irregularities. The singular force F has support only on the immersed interface0(t) and has the form
F(x, t) = Z
0(t)
f(s, t)δ2(x − X(s, t)) ds, (1.5) where X(s, t) is the arc-length parameterization (s is the arc-length parameter) of the inter- face, and f(s, t) is the force strength. The above integral is over the entire interface, and δ2
is the two-dimensional Dirac function,δ2(x) = δ(x)δ(y) with x = (x, y). The term G may have a finite jump across the interface, but is bounded and piecewise continuous in the entire domain. Throughout this paper, we simply assume that the densityρ ≡ 1 and the viscosity µ is a constant. Figure 1 is an illustration of the geometry and the local coordinates for the problems discussed in this paper.
Although we assume that the density and the viscosity are two constants, the model (1.1)–
(1.5) has many applications, for example, cardiac blood flow [32–34], platelet aggregation and coagulation [9, 10], swimming microorganisms [8], biofilm formation [6], suspension flows [38], and others; see also [35] for a brief review.
In the 1970–1980s, Peskin [32, 33] developed the Immersed Boundary Method (IB method) which used the model (1.1)–(1.5) to simulate the blood flow through heart valves.
The valve leaflet (the interface) exerts some force into the blood (the fluid) and moves along with the fluid simultaneously. Peskin’s IB method has become a standard numerical method for interface problems that involve singular forces along one or several interfaces.
The method has been applied successfully to many fluid and biological problems.
However, it is also known that generally the IB method is only first-order accurate for nonsmooth but continuous quantities. Different variations of IB method have been proposed to improve the accuracy. For example, Lai and Peskin [18] proposed a new formally second- order IB method with reduced numerical viscosity and applied the new method to simulate
FIG. 1. A diagram of the geometry for the interface problems discussed in this paper. We use n andτ to denote the unit normal and tangential directions of the interface, respectively; andθ is the angle between the normal direction and the x-axis.
the flow around a circular cylinder. By formal second-order accuracy, the authors meant that the scheme would be second-order accurate if the delta function is replaced by a fixed smooth function, independent of the mesh. The numerical result has a good agreement with the experimental data; however, the fully second-order accuracy in the maximum norm, which is the goal of this paper, has not been achieved there. More recently, Cortez and Minion [5] proposed another high-order scheme using the smoothed blob projection to compute the force. Their method shows better accuracy in L1and L2norms and preserves the volume better than the original IB method. Nevertheless, both approaches described in [5] and [18] use the discrete delta function which usually smears out the solution across the interface.
In contrast to the discrete delta function approach, the immersed interface method (IIM) [20–23, 25], which incorporates the jump conditions into the finite difference scheme, was developed to improve the accuracy of the IB method and to obtain sharp interface solutions. The IIM method has been shown to be second-order accurate in the maximum norm for elliptic problems [26] and has been applied to Stokes flow [21] and other moving interface/free boundary problems [14, 27, 28]. In this paper, we extend the IIM method to the full Navier–Stokes equations (1.1)–(1.5) and compare the new method with different versions of the IB method.
The Ghost Fluid Method (GFM), which is also a sharp interface method that uses the jump conditions in the solution and the flux, has been developed by Fedkiw, Liu, Kang and their group for elliptic interface problems [29]. The GFM is simpler than IIM and the resulting linear system of equations for the self-adjoint elliptic equations is symmetric positive definite. The trade-off is that the solution is first-order accurate in the maximum norm. The method has been applied to multiphase incompressible flow [16] and two phase incompressible flames simulations [30].
Another sharp interface method in which the jump conditions in the coefficients was taken into account in constructing the algorithm was developed in [11] and [12] for the incompressible Navier–Stokes equations.
The method developed in this paper is for sharp interface models. Besides the differ- ences in the methodology, our method distinguishes from [16, 29, 30] in accuracy, and distinguishes from [11, 12] in dealing with different problems. There are some other mod- els and methods for interface problems, notably, the phase field model and the finite volume method, for example, [36]. Each model has advantages and disadvantages, and may have specific applications. It is easy to compute surface tension and curvature with sharp interface models. It is also an advantage to use the level set method with sharp interface models.
The rest of the paper is organized as follows. In Section 2, we introduce the jump condi- tions for (1.1)–(1.5). In Section 3, we present the IIM scheme for (1.1)–(1.5) which modifies a projection method by adding some correction terms. In Section 4, numerical results for nontrivial examples including one with a moving interface are presented and analyzed. A new version of the IB method with the level set formulation is also proposed there. Some conclusions will be given in Section 5.
2. JUMP CONDITIONS ACROSS THE INTERFACE
Because of the singular forces, the velocity of the solution to (1.1)–(1.5) is typically non- smooth, and the pressure may be discontinuous, across the interface. The jump conditions of the velocity and the pressure are quantitatively described in the following theorem.
THEOREM1. Let(X, Y ) be a point on the interface. Let 0 be a smooth closed interface within the domain. Let the unit outward normal direction be n= (cos θ, sin θ), where θ is the angle between the outward normal direction and the x-axis. Then we have the jump conditions
[u]= 0, [µun]= − ˆf2τ, (2.6)
[ p]= ˆf1, [pn]= ∂ ˆf2
∂s + [G] · n, (2.7)
whereτ = (−sin θ, cos θ) is the unit tangent direction, and ˆf1and ˆf2are the force strengths in the normal and tangential directions, respectively
fˆ1(s, t) = f1(s, t) cos θ + f2(s, t) sin θ,
(2.8) fˆ2(s, t) = − f1(s, t) sin θ + f2(s, t) cos θ.
The jump [·] is defined as the difference of the limiting value from outside of the interface and that from the inside, and s is the arc-length parameter of the interface.
The proof can be found in [19, 21] with minor modifications.
2.1. Additional Interface Relations
In this subsection, we derive some additional interface relations by differentiating the known jump conditions in Theorem 1 along the interface, and by making use of the mo- mentum equation (1.1). These interface relations are summarized in Theorem 2 and will be used in our new numerical scheme. For interface problems, most physical quantities along/across the interface are associated with the normal and tangential directions. Thus, we introduce a local coordinate system which lies along those directions, and derive more interface relations in the local coordinate system.
THEOREM2. Under the same notations and assumptions as in Theorem 1, we define the local coordinates at(X, Y ), a point on the interface, as
ξ = (x − X) cos θ + (y − Y ) sin θ,
(2.9) η = −(x − X) sin θ + (y − Y ) cos θ.
Then the interface0 can be represented by ξ = χ(η) in the neighborhood of (ξ, η) = (0, 0), which satisfiesχ(0) = 0, χ0(0) = 0, and χ00(0) = κ, the curvature of the interface at (0, 0).
The following interface relations are true at(X, Y ):
[ p]= ˆf1, [pξ]= ∂ ˆf2
∂η , [u]= 0, [µuξ]= − ˆf2τ, [uη]= 0,
(2.10) [µuηη]= κ ˆf2τ, [µuξη]= −∂ ˆf2
∂ητ − κ ˆf2n, [µuξξ]= −[µuηη]+ [pξ] n+ [pη]τ + [uξ] u· n − [G].
The proof of the theorem is quite technical, long, and tedious. However, it is straightfor- ward if we know the key ideas. Therefore, we will provide only the outline of the proof.
Sketch of the proof. The first four equalities are directly copied from the jump conditions in Theorem 1. Let the interface0 be ξ = χ(η) in the neighborhood of (ξ, η) = (0, 0).
Then we know thatχ(0) = 0, χ0(0) = 0, and χ00(0) = κ. Differentiating [u] = 0 along the interface and usingχ0(0) = 0, we get the fifth equality. By differentiating [u] = 0 along the interface twice and using χ0(0) = 0, and χ00(0) = κ, we get the sixth equality. By differentiating the fourth equality along the interface and using∂τ∂η = κn, we get the seventh equality.1Finally, by expressing the momentum equation (1.1) in the local coordinate (2.9), we get the last equality.
2.2. Projection of the Jump Relations on x- and y-Direction
In order to be numerically useful, those interface relations in (2.10) in the local coordi- nate (2.9) are transformed into the jump relations in the Cartesian coordinates. Sinceµ is continuous, we have
[ux]= [uξ] cosθ − [uη] sinθ, [uy]= [uξ] sinθ + [uη] cosθ,
(2.11) [ux x]= [uξξ] cos2θ − 2[uξη] cosθ sin θ + [uηη] sin2θ,
[uyy]= [uξξ] sin2θ + 2[uξη] cosθ sin θ + [uηη] cos2θ.
Using these identities and the interface relations in (2.10), we can easily write down [u], [ux], [uy], [ux x], [uyy], [ p], [ px], and [ py] at any point on the interface in terms of the force strength and its derivatives, and the geometric information of the interface.
3. THE NUMERICAL ALGORITHM
Our numerical method is based on the projection method for solving the incompressible Navier–Stokes equations with careful treatment at grid points that are near or on the interface.
There are several versions of the projection method for solving the incompressible Navier–
Stokes equations; see for example, [1, 15, 17] and many others. The projection method that we used in this paper is the one described in [2] which is based on the pressure increment formulation of [1, 15] with an additional accurate pressure correction (in [2]; the method is denoted by Pm II). Our key modification to the projection method (Pm II) is to add some correction terms at grid points near or on the interface. Those correction terms are determined from the interface relations that we stated in the previous sections.
For simplicity, we assume that the domainÄ is a rectangle [a, b] × [c, d], and the spatial spacing is h= (b − a)/M = (d − c)/N, where M and N are the number of grid points in the x and y directions, respectively. Here, a standard uniform Cartesian grid is used. Our method from time tnto tn+1can be written as
u∗− un
1t + (u · ∇hu)n+12 = ∇hpn−12+µ
2(1hu∗+ 1hun) + Gn+12+ Cn1,
(3.12) u∗|∂Ä= unb+1,
1In our notation, a circle has negative curvature.
where(u · ∇hu)n+12 is approximated using (u · ∇hu)n+12 =3
2(un· ∇h)un−1
2(un−1· ∇h)un−1+ Cn2. (3.13) The projection step is the following:
1hφn+1= ∇h· u∗
1t + Cn3, ∂φn+1
∂n
¯¯¯¯
∂Ä= 0, (3.14)
un+1 = u∗− 1t∇hφn+1+ Cn4. (3.15)
∇hpn+12 = ∇hpn−12+ ∇hφn+1+ Cn5. (3.16) In the expressions above,∇hand1hare the standard central difference operators regardless of the interface, thus,
1hui j =
µui+1, j− ui−1, j
2h ,ui, j+1− ui, j−1
2h
¶ , 1hui j = ui−1, j+ ui, j−1+ ui+1, j + ui, j+1− 4ui j
h2 .
It is clear that the modifications that we made to the projection method are those correction terms Cknat grid points near or on the interface.
3.1. Determining the Correction Terms
At a regular grid point (xi, yj), meaning that all the grid points in the standard centered five-point stencil are from the same side of the interface, we use the standard central difference scheme to discretize Eqs. (1.1)–(1.4) without any correction; that is Ckn ≡ 0.
At an irregular grid point (xi, yj), we need to determine the correction terms so that the finite difference scheme can be as accurate as possible. We assume that the interface cuts a grid line between two grid points no more than once. This is guaranteed if the interface is expressed in terms of a level set function.
First, we establish the following lemma, which is the basis for determining the correction terms.
LEMMA1. Let u(x) be a piecewise twice differentiable function. Assume that u(x) and its derivatives have finite jumps [u], [ux], and [ux x], at x∗= x + αh, −1 ≤ α ≤ 1, then the following relations hold,
u(x + h) − u(x − h)
2h =
u0(x) + C(x, α)
2h + O(h2), if 0 ≤ α ≤ 1, u0(x) − C(x, α)
2h + O(h2), if −1 ≤ α < 0,
(3.17)
u(x + h) − 2u(x) + u(x − h)
h2 = u00(x) +C(x, α)
h2 + O(h), (3.18) where
C(x, α) = [u] + [ux](1 − |α|)h + [ux x](1 − |α|)2h2
2 , (3.19)
and the jumps are defined as
[u]=
(limx→x∗+u(x) − limx→x∗−u(x), if 0 ≤ α ≤ 1,
limx→x∗−u(x) − limx→x∗+u(x), if −1 ≤ α < 0. (3.20) [ux]=
(limx→x∗+ux(x) − limx→x∗−ux(x), if 0 ≤ α ≤ 1,
limx→x∗−ux(x) − limx→x∗+ux(x), if −1 ≤ α < 0. (3.21)
[ux x]=
(limx→x∗+ux x(x) − limx→x∗−ux x(x), if 0 ≤ α ≤ 1,
limx→x∗−ux x(x) − limx→x∗+ux x(x), if −1 ≤ α < 0. (3.22) Proof. Without loss of generality, we assume that 0< α ≤ 1. Therefore, x − h and x are on the same side while x+ h is on the other side. We use the Taylor expansion twice for u(x + h) at x∗, then at x, as follows,
u(x + h) = u(x∗+ (1 − α)h) = u(x∗+) + ux(x∗+)(1 − α)h + ux x(x∗+)(1 − α)2h2
2 + O(h3)
= u(x∗−) + ux(x∗−)(1 − α)h + ux x(x∗−)(1 − α)2h2
2 + C(x, α) + O(h3)
= u(x) + ux(x)(x∗− x) + ux x(x)(x∗− x)2
2 + ux(x)(1 − α)h + ux x(x)(x∗− x)(1 − α)h + ux x(x)(1 − α)2h2
2 + C(x, α) + O(h3).
The Taylor expansion for u(x − h) is
u(x − h) = u(x) − hux(x) + ux x(x)h2
2 + O(h3).
After substituting these two expansions into the left-hand sides of (3.17) and (3.18) and with some manipulations, we get the desired equalities.
Remark 1. The lemma above is the basis of our numerical algorithm for determining the correction terms.
• The correction terms given in the lemma can be found in the literature. For example, (3.18) is a special case in the expressions (2.22) and (2.24) in [22]. It was also derived independently in [29] for the ghost fluid method. Equation (3.17) is a special case in the expressions (26) and (39) in [24]. Nevertheless, in the lemma we have “uniform” formulae for both first- and second-order derivatives; see below for the explanation.
• Lemma 1 tells us that the correction term is simply the product of the finite difference coefficient, corresponding to the grid point from different sides of the interface in reference to the master grid point, and C(x, α) which corresponds to the singular source term. For example, if−1 < α < 0, in Lemma 1, which means x − h and x are on different sides of the interface, then the correction terms for ux is the product of the coefficient−1/(2h) and C(x, α). Therefore, it is very easy to modify the finite difference scheme to maintain a certain order of accuracy.
• If there are two interfaces between (x − h, x + h), say, x1∗= x + α1h and x2∗= x + α2h, withα1< 0 and α2≥ 0, then we just need to add another correction term, for example,
u(x + h) − 2u(x) + u(x − h)
h2 = u00(x) + C(x, α1)
h2 +C(x, α2)
h2 + O(h), u(x + h) − u(x − h)
2h = u0(x) + C(x, α1)
2h −C(x, α2)
2h + O(h2).
3.2. Correction Terms to the Projection Method
In our numerical scheme (3.12)–(3.16), there are several correction terms to be deter- mined. For the sake of simplicity, we will only explain how to evaluate some of the correction terms of the x-component of Cn1. The other correction terms can be treated in the same way so we omit the details here. Assume that x∗, xi ≤ x∗≤ xi+1, is an intersection of the in- terface and the x-axis. The contribution to Cn1 in the x-direction from this interface point includes the following:
• the correction term for unDx,hunfrom the nonlinear term(un· ∇h)un:
−3 4h
µ£unx¤
(xi+1− x∗) +£
unx x¤(xi+1− x∗)2 2
¶
un(x∗, yj);
• the correction term for un−1Dx,hun−1from the nonlinear term(un−1· ∇h)un−1:
1 4h
µ£unx−1¤
(xi+1− x∗) +£
unx x−1¤(xi+1− x∗)2 2
¶
un−1(x∗, yj);
• the correction term for Dx,hpn−12 from∇hpn−12:
−1 2h
³£pn−12¤ +h
pn−
1
x 2
i(xi+1− x∗)´
;
• the correction term for µ(u∗i−1, j − 2u∗i, j+ u∗i+1, j)/(2h2) from µ1hu∗/2:
− µ 2h2
µ£unx+1¤
(xi+1− x∗) +£
unx x+1¤(xi+1− x∗)2 2
¶
;
• the correction term for µ(uni−1, j − 2uni, j+ uni+1, j)/(2h2) from µ1hun/2:
− µ 2h2
µ£unx¤
(xi+1− x∗) +£
unx x¤(xi+1− x∗)2 2
¶ .
We use the bilinear interpolation to approximate the values of u andv at (x∗, yj) from the values of u andv at the neighboring four grid points. Such approximation is at least first- order accurate since the velocity is continuous. Note that we have used the jump conditions at tn+1 to approximate the jump conditions at t∗ corresponding to the same treatment for the velocity boundary condition in the projection method.
Similarly, if the interface cuts through between xi−1 and xi, we need to add the corre- sponding correction terms as well. In this way, we take care of the correction terms in the x-direction for Cn1.
The spatial local truncation error of our scheme is O(h2) at regular grid points, but O(h) at irregular grid points. However, since the number of irregular grid points usually is one dimensional lower than the total number of grid points, we expect the global accuracy is still second order for the velocity, and nearly second order for the pressure. A proof of second-order convergence for an elliptic interface problem can be found in [26]; see also the numerical experiments in Section 4.
It is worth pointing out that the modifications that we made to the projection method only affect the right-hand sides of the original projection method. Therefore, our method does not change the stability nature of the original projection method, which is stable if the CFL condition is satisfied.
Once we have computed those correction terms, we need to solve two Helmholtz equations to get u∗ in (3.12) and one Poisson equation to getφn+1 in (3.14). The discrete systems for the Helmholtz and Poisson equation are symmetric and positive definite, and are solved using the fast solver HWSCRT from Fishpack.
3.3. Further Correction Near the Boundary and the Interface
Using the projection method with the correction terms added at irregular grid points, we have observed that the velocity is second-order accurate but the pressure is only first order.
The largest error in the pressure appears near the boundary or the interface. As analyzed in [2, 7, 37, 39], this is due to the inherent numerical boundary layer in the projection method.
In [2], see also [17], a correction scheme is proposed which is
pn+12 = pn−12+ φn+1−µ
2(∇h· u∗+ C6n), (3.23) whereφn+1, in our notation, is the solution of the Poisson equation of (3.14)–(3.15) from tn to tn+1, C6n again is a correction term only needed at irregular grid points near or on the interface. This correction does dramatically improve the accuracy for the pressure (see Section 4).
In our implementation, we use the centered difference scheme to compute the correction term∇h· u∗ in interior grid points. Then we use an extrapolation scheme to extend the quantity to the boundary of the rectangle, for example, at the boundary x = a,
∇h· u∗(0, j) = 3∇h· u∗(1, j) − 3∇h· u∗(2, j) + ∇h· u∗(3, j). (3.24)
Near the interface, even though the velocity is second-order accurate, the term∇h· u∗/1t may have only zeroth order accuracy since the error may not be smooth near the interface, and1t ∼ h. This will reduce the accuracy of pressure to first order. Our remedy to this problem is to replace∇h· u∗/1t at an irregular grid point by the same quantity at the nearest regular grid point. This again improves the accuracy of the pressure. The same technique was also used independently in [4] where a second-order accurate symmetric sharp interface method for the Poisson equation with Dirichlet boundary conditions was proposed.
We will call the numerical method described in this section as the projection-IIM method.
4. NUMERICAL EXAMPLES
In this section, we present some numerical experiments for solving the velocity and presure in (1.1)–(1.4) at some fixed time t . The velocity, and the presure up to a constant, are uniquely determined by the initial conditions u(x, 0), p(x, 0), and 0(s, 0); the initial interface; the boundary condition u(x, t)|∂Ä; and the force strength( ˆf1, ˆf2) in the normal and tangential directions along the interface.
All the calculations in this section were perfomed at North Carolina State University using Sun workstations. Most simulations are done within minutes to a couple of hours depending on the number of grid points. Unless specified differently, we use a level set function to represent the interface; see (4.40).
4.1. Example 1: A Noninterface Problem
We start our numerical tests by checking the accuracy of the projection method for a regular problem, and use the result as the basis for comparing different algorithms for the interface examples. The following exact solution is taken from [40]:
u(x, y, t) = −sin2(πx) sin(2πy) cos t,
(4.25) v(x, y, t) = sin2(πy) sin(2πx) cos t.
Define a functionψ by
ψ =
·
(x − x2) sin(πx)(y − y2) sin(πy) − 16 π6
¸
cos t. (4.26)
The pressure then is
p(x, y, t) = −∂ψ
∂t + ν1ψ. (4.27)
The force term G is determined from the exact solution. The domain is the unit square Ä = [0, 1] × [0, 1]. In this example, we have∂p∂n 6= 0 on the boundary. Thus, the numerical boundary layer might deteriorate the accuracy of the pressure.
Table I shows the grid refinement analysis at t= 10. Throughout this paper, the error computed in the maximum norm is up to the boundary, meaning that the error at the boundary is also taken into account. The order of accuracy is defined as
order=log[E(N)/E(2N)]
log 2 . (4.28)
The pressure has been adjusted by some constant that minimizes the error in the maximum norm. We can see the velocity is second-order accurate, but the pressure is first order using the original projection method. However, if we add the correction term in (3.23), then the pressure is very close to second-order. The last column in Table I is the errors of the pressure with the correction term in (3.23). While it is not fully second-order accurate, it is better than the original projection method without the correction. Away from the boundary, the pressure is second-order accurate.
TABLE I
The Grid Refinement Analysis of the Projection Method withµ = 1, ∆t = h, for Example 1 at t = 10
Projection method without (3.23) With (3.23)
N kE1(u)k∞ Order kE1(p)k∞ Order kE2(p)k∞ Order
32 2.7021 × 10−4 0.0608 3.2871 × 10−2
64 8.3429 × 10−5 1.6955 0.0269 1.1785 9.2432 × 10−3 1.8303
128 2.2471 × 10−5 1.8925 0.0109 1.2975 3.2825 × 10−3 1.4936
256 5.7841 × 10−6 1.9580 0.0045 1.2665 1.1341 × 10−4 1.5333
Note. The second columnkE1(u)k∞, and the fourth columnkE1(p)k∞are the errors of the velocity and pressure using the projection method without correction in (3.23), respectively. While the velocity is nearly second order accurate, the pressure is roughly first order. The sixth columnkE2(p)k∞is the error for the pressure using the projection method with correction in (3.23). While it is not entirely second order, the accuracy of the pressure has been improved. The velocity (not listed) remains second-order accurate with smaller errors.
In practice, it is difficult to get fully second-order accurate results for pressure for some problems using projection methods with a single Cartesian grid. It may be due to the difficulty in dealing with the numerical boundary layer and the four corners of the rectangular domain. We did notice that, however, the full second-order accuracy was achieved for pressure in [2] for the channel flow and a problem with homogeneous boundary conditions.
We think it may be possible to achieve full second-order accuracy for pressure if a suitable extrapolation scheme is used to evaluate the correction term at the boundary, particularly, at the four corner points. Usually, it is more accurate to use the MAC grid than to use a single Cartesian grid. So we believe that it is possible to achieve full second-order accuracy using the MAC grid. The trade-off is the increased complexity in the implementation.
4.2. Example 2: An Interface Problem with a Constant Jump in the Pressure In this example, we consider a fixed interface,
(x − 0.5)2
a2 +(y − 0.5)2
b2 = 1, (4.29)
in the domain [0, 1] × [0, 1]. The normal and tangential force strength are ˆf1= 10 and fˆ2= 0, respectively. We also set G1= 0 and G2= 0. The solution now is given by
u= 0, v = 0, p =
C, if(x − 0.5)a2 2 +(y − 0.5)b2 2 − 1 > 0,
−10 + C, if(x − 0.5)a2 2 +(y − 0.5)b2 2 − 1 ≤ 0,
(4.30)
where C is an arbitrary constant. Here, we compare our method with Peskin’s IB method using the discrete cosine delta function. Depending on how the interface is represented, there are two different versions of the IB method.
4.2.1. The Lagrangian particle approach. In this traditional approach, the interface is represented by a set of Lagrangian points Xk, k = 1, 2, . . . , Nb, where Nb is the total number of the particle points and is usually chosen in such a way that
max1sk= max kXk+1− Xkk ∼ h.
The forces defined on those particle points are then distributed to the nearby grid points (x, y) by the formula
Nb
X
k=1
f(sk)δw(x − Xk)δw(y − Yk)1sk,
whereδwis the one-dimensional discrete delta function
δw(x) = ( 1
4w(1 + cos(πx/2w)), if |x| < 2w,
0, if|x| ≥ 2w. (4.31)
Although there are other discrete delta functions in the literature, the computational results are not much different for two- and three-dimensional problems. The most common choices ofw is h, the spatial mesh size. As we can see, when w gets larger, the cost to spread the singular force to the grid points increases significantly. Our numerical tests for this and other examples show that the accuracy of the computed velocity remains pretty much the same for different choices ofw. However, the pressure, if it is discontinuous, smears out more widely asw increases. Therefore, the best choice of w is h. The sixth column in Table II is the grid refinement analysis for the velocity which is clearly first-order accurate.
4.2.2. The level set representation. If the interface is represented by the zero level set of a two-dimensional functionϕ(x, y) as in (4.40), it is not straightforward to use the IB method directly. However, if the tangential component of the singular force is zero, then we can write
Z
0(t)f(s)δ2(x − X(s)) ds = f · nδ(ϕ)∇ϕ; (4.32) see [3] and others. It is easy to apply the discrete delta function in the above form.
TABLE II
The Grid Refinement Analysis and Comparison of the Three Projection-IB Methods for Example 2 withµ = 0.1, a = 0.35, b = 0.25, and ∆t = h
Level Set IB widthw = h Level Set IB widthw =√
h Lagrangian IB widthw = h
N kE(u)k∞ Order kE(u)k∞ Order kE(u)k∞ Order
32 1.2434 × 10−1 4.9254 × 10−3 5.1204 × 10−2
64 6.3619 × 10−2 0.96673 2.9594 × 10−3 0.7349 2.5839 × 10−2 0.9867 128 4.0730 × 10−2 0.64339 1.1062 × 10−3 1.4197 1.2968 × 10−2 0.9929 256 4.3059 × 10−2 0.09459 4.1558 × 10−4 1.4124 6.5055 × 10−3 0.9952
Note. The initial data is taken from the solution in Example 1. The second and fourth columns are the error of the computed velocity using the level set function with widthw being h and√
h, respectively. The sixth column is the error of the computed velocity using the original IB method withw being h. Generally, the velocity is first-order accurate. The pressure which is not listed in the table has zeroth order accuracy because it is smeared out.
Unfortunately, the choice of the width of the discrete delta function is crucial for the IB method. Ifw = h, the IB method barely converges; see the second column of Table II. T.
Hou [13] may be the first one to note this and suggest that the best choice of the width is w =√
h. This has also been confirmed in our numerical tests; see the fourth column of Table II. However,w =√
h is much larger than h if h is small. The IB method smears out the pressure and shows zero-order convergence for the pressure. Since the width has to be taken asw =√
h, we cannot guarantee the right pressure even at those grid points that are away from the interface.
A remedy to the problem is to use the projections of irregular grid points of a particular side, sayϕ(x, y) ≤ 0, on the interface as the control points. Then the Peskin’s original IB method can be applied with those control points. In [14, 28], we described in detail how to find the projection of an irregular grid point on the interface. Using this new approach, we still can use a thin layer(w = h) of the interface in the IB method, so that the solution will only be smeared out in the thin layer while we can keep the advantages of the level set method.
4.2.3. The grid refinement analysis and comparison. Table II lists the comparison of three different versions of the IB method using the particle approach and the level set func- tion. We have also tested the example with different but modestµ, different major axis a and minor axis b, and different initial conditions. All the results are qualitatively similar.
In Table III, we show the grid refinement analysis of the projection-IIM method devel- oped in this paper. Without the modification of the pressure equation (3.23), the velocity in the second column is second-order accurate, while the pressure in the fourth column is first-order accurate with the largest errors occurring near the boundary or the interface.
With the correction term in (3.23) added, we see second-order accuracy for both the ve- locity and the pressure. Since the pressure is piecewise constant in this case, the accuracy for the pressure is much better than second order. However, this is not true for general problems.
In a further numerical experiment, we revisit Example 1 but add a singular force term ˆf1= 10 to the Navier–Stokes equations in (1.1). The exact solution for the velocity is the same
TABLE III
The Grid Refinement Analysis of Projection-IIM Method for Example 2. The Parameters areµ = 0.1, a = 0.35, b = 0.25, and ∆t = h
Projection IIM without (3.23) Projection IIM with (3.23)
N kE(u)k∞ Order kE(p)k∞ Order kE(p)k∞ Order
32 1.6412 × 10−4 1.5785 × 10−2 9.9058 × 10−3
64 3.9857 × 10−5 2.0418 7.9281 × 10−3 0.9935 1.1533 × 10−3 3.1026 128 9.8987 × 10−6 2.0095 3.9676 × 10−3 0.9987 1.5946 × 10−4 2.8544 256 1.2367 × 10−6 2.0008 1.9842 × 10−3 0.9997 2.4832 × 10−5 2.6830
Note. The first and third columns are the errors of the computed velocity and the pressure without the correction term in (3.23), respectively. We can clearly see that the velocity is second-order accurate while the pressure is first-order because of the boundary layer. The sixth column is the error of the computed pressure using the same method but with the correction term for the pressure in (3.23).
The accuracy actually is better than second-order because of the fact that the pressure is piecewise constant. The computed velocity (not shown in this table) with the correction is also more accurate meaning a smaller error constant.
FIG. 2. Computed pressure for Example 1 with an added force term ˆf1= 10 using a 64 by 64 grid. (a) The result computed using the projection-IIM method which catches the jump in the pressure. (b) The result computed using the IB method with the interface represented by a level set function(w =√
h). The width of the discrete delta function isw =√
h. (c) The result computed using the original IB method(w = h). The error for the pressure is O(1) for both (b) and (c).
as in Example 1, but the pressure is the sum of the pressure in Example 1 and the piecewise constant in (4.30). The analysis and observations are similar to the discussions above.
Figure 2a shows the pressure computed using the projection-IIM method. The computed pressure has the right jump up to second order. Figures 2b and 2c show the pressure computed using the IB method with the level set representation of the interface(w =√
h) and the original IB method with particle representation of the interface(w = h), respectively. In Fig. 2b, we see that the pressure is smeared out by the widthw =√
h while in Fig. 2c, we see larger errors but a thinner layer(w = h) near the interface.
For this test problem, our projection-IIM method is faster compared with the IB method.
The reason, we believe, is the triple loop of the IB method to spread the forces to nearby grid points.
4.3. Example 3: An Interface Problem with Nonconstant Jump in the Pressure
In this example, we want to check how well our method can handle a nonconstant jump in the pressure and a finite jump in G. The interface0 and the domain Ä are the same as in Example 2.
The source term G is the gradient of a function q(x, y), which is the solution of the following
1q = 0, (x, y) ∈ Ä, ∂q
∂n
¯¯¯¯
∂Ä= 0, (4.33)
[q]|0= x + y,
·∂q
∂n
¸¯¯¯¯
0= 0. (x, y) ∈ 0. (4.34) The function q(x, y) and G = ∇q are computed using the IIM method [25] with a 512 by 512 fine grid. Note that from [∂q∂n]= 0 and [q] = x + y, it is easy to derive
[G1]= [qx]= sin2θ − sin θ cos θ, [G2]= [qy]= cos2θ − sin θ cos θ, where n= (cos θ, sin θ) is the unit normal direction of 0 pointing outward.
The full Navier–Stokes equations of Example 3 are (1.1)–(1.5) with the following parameters
u|∂Ä= 0, µ = 1, G = ∇q(x, y), fˆ2= 0, ˆf1= x + y, (x, y) ∈ 0, with initial conditions for u and p.
For this example, we do not know the exact solutions. However, the steady state solution for this example is
tlim→∞u(x, t) = 0, lim
t→∞p(x, t) = q(x).
In Table IV, we show the grid refinement analysis for the steady state solution obtained using the projection-IIM method at t= 10. The initial velocity is taken from (4.25) in Example 1, while the pressure is
p0(x, y) = cos πx cos πy + q(x, y).
In this way, we guarantee the initial pressure has the right jump, [ p]= ˆf1. We see that the projection-IIM method gives a second-order accurate solution for both the velocity and the pressure. The results are qualitatively the same if different initial conditions are used, as
TABLE IV
The Grid Refinement Analysis of Projection-IIM Applied to Example 3 withµ = 1, ∆t = h, at t = 10
N kEI I M(u)k∞ Order kEI I M(p)k∞ Order
32 2.4930 × 10−5 5.1614 × 10−4
64 3.5197 × 10−6 2.8244 8.4365 × 10−5 2.6131
128 4.6660 × 10−7 2.9152 1.6308 × 10−5 2.3710
256 3.5073 × 10−8 2.9332 5.7023 × 10−6 1.5160
Note. The initial velocity is taken from (4.25) in Example 1. The initial pressure is taken as cosπx cos πy + q(x, y).
TABLE V
The Grid Refinement Analysis of Projection-IIM Applied to Example 4 withµ = 0.02, ∆t = h, at t = 10; Second-Order Convergence Is Achieved
N kEIIM(u)k∞ Order kEIIM(p)k∞ Order
32 2.4215 × 10−3 1.1513 × 10−2
64 5.3547 × 10−4 2.1771 3.2255 × 10−3 1.8357
128 1.4970 × 10−4 1.8388 9.1307 × 10−4 1.8357
256 3.6173 × 10−5 2.0491 1.9727 × 10−4 2.2106
long as the jump condition in p is initially satisfied. Note that cosπx cos πy satisfies the homogeneous Neumann boundary condition∇ p · n = 0 on the boundary. If we take the initial pressure from (4.27) in Example 1, then the computed pressure is only first-order accurate since [ p]6= ˆf1initially while the velocity is still second-order accurate.
4.4. Example 4: Circular Flow with a Line Force
Now we consider an example with nonsmooth velocity. The exact solution is
u(x, y, t) =
½h(t)¡y
r − 2y¢
if r >12,
0 r ≤ 12, (4.35)
v(x, y, t) =
½h(t)¡
−xr + 2x¢
if r> 12,
0 r ≤ 12, (4.36)
p(x, y, t) = 0, (4.37)
where r=p
x2+ y2. The interface is the circle r= 12and the solution domain is [−1, 1] × [−1, 1]. It is easy to verify that the velocity satisfies the incompressibility constraint, and it is continuous but has a finite jump in the normal direction
[un]= [ur]= −2h(t) sin θ, [vn]= [vr]= 2h(t) cos θ, across the interface. Thus, the normal and tangential force strength are
fˆ1 = 0, ˆf2= −2h(t).
Note that, ut· n is not zero on the boundary ifdhdt 6= 0, which may affect the performance of the projection method.2Outside of the interface r= 0.5, the nonzero and bounded source term G is derived directly from the exact solution using Maple. There is also a finite jump in G.
With h(t) = 1, our method converges to the steady state solution. The grid refinement analysis reveals similar results as in Example 3. In our time-dependent test, we take h(t) = 1− e−t. We present the numerical results at t = 10 in Table V and plot the x-component
2It seems that the projection method used here requires ut· n to be zero along the boundary. Numerical tests showed that if ut· n is not zero, then the velocity still converges but not with second-order accuracy; the pressure may not converge at all for the projection method if utis large on the boundary. However, some other projection methods may work even if the condition ut· n = 0 is violated.
FIG. 3. Plot of the x-component of the velocity−u(x, y) of Example 4 at t = 10. The results are computed with a 64× 64 grid and µ = 0.02. The time step size is 1t = h.
of the velocity in Fig. 3. One can easily see that the velocity is second-order accurate, and the pressure is nearly second-order accurate in Table V.
We also tested a slightly different example by setting
G= 0, ˆf1= 0, ˆf2= 10 u|∂Ä= 0.
The domain and the interface are the same. The initial velocity and the pressure are taken to be zero on the rectangular domain. Since the force is only along the tangential direction, the pressure is continuous, but the normal derivative of the velocity has a nonconstant jump across the interface. While the exact solution is difficult to find, the motion of the steady state is nothing but a simple rotation along the interface. In Fig. 4a, we plot the x-component
FIG. 4. The computed steady state velocity for the modified example in Example 4 with a circular interface r= 0.5 with a 64 × 64 grid. Other parameters are G = 0, ˆf1= 0, ˆf2= 10, and µ = 0.02. The time step size is 1t = h. (a) Plot of −u(x, y) at t = 10 which is not smooth. (b) The plot of the velocity field at t = 10, the flow approaches a steady rotation along the interface.
of the velocity−u(x, y) at t = 10 computed using a 64 by 64 grid. We can clearly see the nonconstant jump in the normal derivative. In Fig. 4b, we plot the velocity field at t= 10.
The flow approaches a steady rotation along the circular interface.
4.5. Example 5: A Moving Interface Problem
Now we consider a moving interface problem which has uniform density and viscosity.
The initial interface is given by
r(θ) = r0+ ² sin(kθ), 0 ≤ θ ≤ 2π. (4.38) The initial velocity and the pressure are all set to be zero. The only driven force of the fluid motion is the surface tension which is proportional to the curvatureκ, that is,
fˆ1(0, t) = ²κ, ˆf2(0, t) = 0. (4.39)
The velocity is smooth but the pressure has a nonconstant jump [ p]|0(t)= ²κ, [pn]|0(t)= 0.
In our test, we take r0= 0.5, k = 5, ² = 0.05, and µ = 0.1.
We use the level set method [31] to update the interface. We refer the readers to [14, 27, 28] for the details about how to combine the immersed interface method with the level set method. In the level set method, the interface is the zero level set of a two-dimensional function
ϕ(x, y, t)
<0, if (x, y) is in the inside of the interface,
=0, if (x, y) is on the interface,
>0, if (x, y) is in the outside of interface.
(4.40)
At each time step, we apply the projection-IIM method described in Section 3 to compute the velocity. Then the velocity is used to solve the Hamilton–Jacobian equation
ϕt+ Eu·∇ ϕ = 0· (4.41)
to obtain the new level set function at t= tn+1. Since we use the level set method to update the interface, the time step is chosen as
1t = min
½ h
2kuk∞, h
¾
(4.42)
to maintain the stability for the level set method and the accuracy for the projection-IIM method.
In Fig. 5, we plotted the computed interface at different times with a 160 by 160 grid.
Using the immersed interface method, we know the intersections of the interface and the grid lines. The set of intersections, which is on the interface at a particular time, is
An animation of the evolution of the interface in this example can be found at http://www4.ncsu.edu/
˜zhilin/research.html.
