An augmented approach for Stokes equations with a discontinuous viscosity and singular forces

(1)

An augmented approach for Stokes equations with a

discontinuous viscosity and singular forces

Zhilin Li

a,*

, Kazufumi Ito

a

, Ming-Chih Lai

b

a_{Center for Research in Scientiﬁc Computation & Department of Mathematics, North Carolina State University, Raleigh, NC 27695-8205, USA} b_{Department of Applied Mathematics, National Chiao Tung University, Hsinchu 30050, Taiwan}

Received 5 September 2005; received in revised form 4 February 2006; accepted 13 March 2006 Available online 27 June 2006

Abstract

For Stokes equations with a discontinuous viscosity across an arbitrary interface or/and singular forces along the interface, it is known that the pressure is discontinuous and the velocity is non-smooth. It has been shown that these discontinuities are coupled together, which makes it difficult to obtain accurate numerical solutions. In this paper, a new numerical method that decouples the jump conditions of the fluid variables through two augmented variables has been developed. The GMRES iterative method is used to solve the Schur complement system for the augmented variables that are only defined on the interface. The augmented approach also rescales the Stokes equations in such a way that a fast Poisson solver can be used in each iteration. Numerical tests using examples that have analytic solutions show that the new method has average second order accuracy for the velocity in the infinity norm. An example of a moving interface problem is also presented.

1. Introduction

The incompressible Stokes or Navier–Stokes equations with a discontinuous viscosity and singular forces arise from many important applications in fluid and bio-fluid mechan-ics. One particular example is Peskin’s immersed boundary (IB) model that was introduced to simulate the blood flow in a human’s heart[28,29,37]. The idea of IB formulation is to treat the complicated immersed boundary (such as a heart valve) as a force generator in the fluid domain, or mathematically, a Dirac delta function force distribution along the immersed boundary. In this paper, we consider the following two-dimensional stationary Stokes equations: rp ¼ r l ru þ ðruÞ Tþ g þ Z C fðsÞd2ðx XðsÞÞ ds; x2 X; ð1:1Þ r u ¼ 0; x2 X; ð1:2Þ

where u = (u, v) is the velocity vector, p is the pressure, C is an arbitrary interface parameterized by the arc-length s, and g(x) is an external force. The force density f is only de-ﬁned on the interface C1that separates two ﬂuid regions X and X+. The viscosity l is assumed to be a piecewise constant l¼ l þ_; _{if x}_{2 X}þ_; l_; _{if x}_{2 X}_: ð1:3Þ We refer the readers toFig. 1(a) for an illustration of the problem. The existence of the solution to the system

(1.1)–(1.2)can be found in[35].

There are at least two diﬃculties in solving (1.1)–(1.2)

numerically using finite difference methods. The first one is to deal with the singular force term. A simple way is to use Peskin’s discrete delta function approach to distribute the singular force to nearby grid points. Such a discretiza-tion is typically first order accurate and will smooth out the

*

Corresponding author.

E-mail addresses:zhilin@math.ncsu.edu(Z. Li),kito@math.ncsu.edu

(K. Ito),mclai@math.nctu.edu.tw(M.-C. Lai).

1 _{The singular source term}R

CfðsÞd2ðx XðsÞÞ ds can also be written as

((f Æ n)n + (f Æ s)s)d(C), or in the form of ((f Æ n)n + (f Æ s)s)j$ujd(u), where uis a level set function whose zero level set represents the interface C.

(2)

solution. The second diﬃculty is how to deal with the dis-continuous viscosity. A simple smoothing method may introduce large errors, see[21]for a one-dimensional exam-ple there.

In the case of a continuous viscosity with a Dirac delta function source distribution along an interface, various methods have been developed in the literature. We refer the readers to [4,5,9,18,26,37] for various methods and the references therein. The diﬃculty with a discontinuous viscosity is that the jump conditions for the pressure and velocity are coupled together, see(2.7)–(2.10), which makes it diﬃcult to discretize the system accurately.

In this paper, we develop a new second-order sharp interface method that uses the exact jump conditions for the two two-phase Stokes equations(1.1)–(1.3)with a dis-continuous viscosity. The idea is to introduce two aug-mented variables that are defined only along the interface so that the jump conditions can be decoupled and the immersed interface method (IIM) [17,19] can be applied. The GMRES iterative method is applied to the Schur com-plement system for the discrete augmented variables. Fur-thermore, our approach rescales the original problem and enables us to use a fast Poisson solver in the iterative pro-cess. Each GMRES iteration requires solving the rescaled Stokes equations, which can be done by calling a fast Pois-son solver three times, and an interpolation scheme to eval-uate the residual of the Schur complement. While augmented approaches have been developed for elliptic interface problems or problems defined on irregular domain problems in [3,8,10,20,23,24,40], and for dealing with the pressure boundary condition in [22], the aug-mented approach proposed in this paper provides a way to get a second order immersed interface method using a finite difference discretization to solve the Stokes equations

(1.1)–(1.3)with a discontinuous viscosity.

The main focus of this paper is to test and implement our proposed augmented approach to the two-phase Stokes equations. Our implementation and results reported in this paper are based on the cubic spline representation of the interface [19]. In order to apply the method for prob-lems on multi-connected domains and moving interface problems, the level set method is preferred, see, for exam-ple,[24,27,32]. The level set method using the augmented

approach for the two-phase Stokes equations has yet to be developed.

A few other ﬁnite diﬀerence methods may be applicable for solving the Stokes equations(1.1)–(1.3) with a discon-tinuous viscosity. A simple approach is the smoothing method that treats the viscosity as a continuous function using

lðxÞ ¼ l_{þ ðl}þ_l_ÞH

ðuðxÞÞ; ð1:4Þ

where u(x) is a two-dimensional Lipschitz continuous function whose zero level set is the interface C, H is a

smoothed Heaviside function, see, for example, [33,34, 36]. In this way, a standard Stokes solver on a rectangular region can be used. The disadvantage is that the solution is smeared in a neighborhood of the interface and it is usually first order accurate. Among sharp interface methods, that is, methods for which the computed solution satisfies the jump conditions across the interface either exactly or approximately, the fast Stokes solver on irregular domains based on an integral equation coupled with a finite differ-ence discretization[2,6]may be applicable to the two-phase Stokes flow discussed here. In [25], an immersed interface method that requires three fast Poisson solvers within each iteration was developed assuming that the jump in the viscosity is small so that the jump conditions can be approximated with those for a continuous viscosity, that is, the coupling terms in (2.7)–(2.10) are ignored. Such approximation may not be valid anymore if the jump in the viscosity is large. The ghost fluid method developed in[14]for the two-phase Navier–Stokes equations is appli-cable to the two-phase Stokes equations with first order accuracy. It is also possible to generalize the finite volume methods developed for elliptic and Navier–Stokes equa-tions on irregular domains in[1,12]to the two-phase Stokes flow.

Motivated by the augmented approach (fast IIM) [20]

for elliptic interface problems with piecewise constant coef-ﬁcients, Bube and Wiegmann developed the explicit immersed interface method (EJIIM) for elliptic interface problems in[38,40]. The EJIIM introduces unknown jumps in the solution and its up to second order derivatives along the interface. The GMRES method is then applied for the unknown jumps. The EJIIM has been applied to nonlinear 1D interface problems in [39]and two-dimensional elastic equations in shape design in conjunction with the level set method [31]. While we believe that the EJIIM should be applicable to the two-phase Stokes ﬂow, it has not yet been done and the derivation and implementation are not trivial. Furthermore, the augmented system of equations using EJIIM will be larger than the augmented system in this paper because more unknowns will be introduced as augmented variables.

Our method is a new finite difference method that re-scales the problem, has second order accuracy in the max-imum norm, and takes advantage of a fast Poisson solver. These are a few merits about the new method that is competitive against other finite difference methods such Fig. 1. (a) A diagram for the incompressible Stokes equations defined on a

domain X with an interface C across which the viscosity l is discontin-uous. (b) Force density decomposition in which f1and f2are the force

density in the x- and y-directions, while ^f1and ^f2are the force density in

(3)

as the one that uses smoothing technique and a discrete delta function.

While augmented methods have some similarities to boundary integral methods to ﬁnd a source strength, the augmented methods have a few special features: (1) no Green function is needed, and therefore no need to evalu-ate singular integrals; (2) no need to set up the system of equations for the augmented variable explicitly; (3) appli-cable to general PDEs with or without source terms; (4) the process does not depend on boundary conditions. On the other hand, we have less information about the condi-tion number of the Schur complement system and how to apply pre-conditioning techniques. A boundary integral method requires less computation compared with the aug-mented method if only the velocity on the interface is needed. The augmented method proposed here requires less computation compared with a boundary integral method if we need the velocity in the entire domain as in many applications.

The paper is organized as follows. In the next section, we discuss the jump conditions of the Stokes equations with a discontinuous viscosity along an interface. We will see how the jump conditions for the pressure and velocity are coupled together, and explain how to decouple the jump conditions by introducing two augmented variables and two augmented equations. In Section 3, we present the new algorithm in detail. We explain how the GMRES iterative method can be applied to the discrete augmented variables without explicitly forming the coeﬃcient matrix. In Section4, we present some numerical experiments using examples that have analytic solutions to check the conver-gence and the performance of the new method. An example of moving interface is also presented there.

2. Jump conditions for the Stokes equations with a discontinuous viscosity and singular forces

Referring to Fig. 1(a) and Eqs.(1.1)–(1.3), we assume the interface C is a smooth curve. We use n and s to repre-sent the unit normal and tangential directions at a point (X, Y) on the interface. We use ^f1ðsÞ ¼ fðsÞ n and

^

f2ðsÞ ¼ fðsÞ s to represent the force density in the normal

and tangential directions. If C is smooth, then in a neigh-borhood of C, the distance function d(x, C) is also a smooth function. The normal and tangential directions of C can be extended to the neighborhood, for example, n = $d(x, C)/ j$d(x, C)j. Therefore the quantities such as u Æ n, ou

on, ou os,

etc., are well deﬁned in the neighborhood of C.

The interface conditions for the two-phase Stokes equa-tions with a discontinuous viscosity and a Dirac delta func-tion singular force are derived in[11]. The idea is from the immersed interface method to express all the quantities in the local Cartesian coordinate at a point (X, Y) on the interface C as,

n¼ ðx X Þ cos h þ ðy Y Þ sin h; g¼ ðx X Þ sin h þ ðy Y Þ cos h;

ð2:5Þ

where h is the angle between the x-axis and the normal direction at the point (X, Y), seeFig. 1(b) for an illustra-tion. In the local coordinate system, the interface can parameterized by n = v(g), g = g. Note that v(0) = 0, v0_{(0) = 0, v}0_{(0) = j, the curvature of the interface at the}

point (X, Y). The unit normal and tangential directions are n¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 1þ ðv0_ðgÞÞ2 q ; v 0_ðgÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1þ ðv0_ðgÞÞ2 q 0 B @ 1 C A; s¼ v 0_ðgÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1þ ðv0_ðgÞÞ2 q ; _{ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi}1 1þ ðv0_ðgÞÞ2 q 0 B @ 1 C A: ð2:6Þ

We also express the velocity components in the tangential and normal directions as ^u¼ u n and ^v ¼ u s.

The jump conditions at a ﬁxed point (X, Y) (or (0, 0) in the local coordinate system) are summarized in the follow-ing theorem.

Theorem 1. Assume C(s)2 C2, ^f₁ðsÞ 2 C1, and ^f₂ðsÞ 2 C1. Let p and u be the solution to the Stokes equations(1.1)– (1.2). We have the following jump conditions across the interface C at the fixed point (X, Y) in terms of the local Cartesian coordinate system

½p ¼ 2 loû on þ ^f1; ð2:7Þ op on ¼ ½g n þo^f2 og þ 2 l o2û og2 4j lo^v og ; ð2:8Þ lo^v on þ loû og þ ^f2¼ 0; ð2:9Þ ½lr u ¼ 0; ð2:10Þ ½u ¼ 0; ð2:11Þ

where the jump [ Æ ] of a quantity, for example, [p] at a point X is defined as

½pj_X2Cdef¼ lim

x!X;x2XþpðxÞ _x!X;x2Xlim pðxÞ: ð2:12Þ Note that, although the jump condition op_on is not needed to close the system mathematically, the condition is used in our numerical algorithm to obtain second-order discretization for the pressure. The jump condition is not arbitrary but determined from the governing equations and the force strength. This is similar to the discussion of the boundary condition of the pressure when a no-slip boundary condition is prescribed for the velocity.

We will omit the subscript X2 C in the rest of the paper for simplicity. The ﬁrst jump condition(2.7)is the result of balancing force in the normal direction while the third one

(2.9)is the result of balancing force in the tangential direc-tion. The jump condition (2.10) is obtained from the incompressibility condition $ Æ u = 0. The jump condition

(4)

for the normal derivative of the pressure (2.8) can be obtained by applying the divergence operator to the momentum equation(1.1),

rp ¼ r lðru þ ðruÞTÞ þ g

excluding the interface, see[11]for the details.

From the incompressibility condition, one can easily prove that o^u on h i ¼ ou on n ¼ 0. If l is continuous, then lo^u on h i ¼ 0 and 2 lo2_^_u og2 h i 4j lo^v og h i

¼ 0, and the jump condi-tions for the pressure and the velocity are decoupled in

(2.7)–(2.10). In this case, we recover the jump conditions derived and used in [16,18]. A second order accurate immersed interface method has been developed in[18,19]

if the viscosity is continuous.

Note that using the invariance of the ﬁrst order deriva-tives, the jump conditions(2.7) and (2.9)can also be writ-ten as ½p ¼ 2 lou on n þ ^f1; ð2:13Þ lou on s þ lou os n þ ^f2¼ 0: ð2:14Þ

These jump conditions can be found in the literature. To get a second order accurate algorithm based on the immersed interface method for the Stokes equations with a discontinuous viscosity, our strategy is to introduce two intermediate, or augmented variables, along the interface so that the jump conditions can be decoupled. We also need two augmented equations to close the system of the equations.

There are diﬀerent ways to introduce augmented vari-ables so that the jump conditions can be decoupled. For example, one could introduce lou

on

as an intermediate unknown so that the jump conditions in the pressure and the velocity are de-coupled. Different augmented variables and equations will lead to different algorithms. Note that for a discontinuous viscosity, there are two different scales corresponding to the two-phase flow. Generally, if the vis-cosity ratio max{l, l+}/min{l, l+} is large, then the velocity field with smaller viscosity has larger gradient than that with larger viscosity. Furthermore, with a constant vis-cosity and a periodic boundary condition, we can use a fast Poisson solver to solve the stationary Stokes equations, see

[18,19]. Based on these two considerations, we introduce the jumps [lu](s) and [lv](s) along the interface as two aug-mented variables. The advantages and details can be seen in the rest of the paper.

Using this local coordinate system, we have o on¼ o onand o os¼ o

og. We can rewrite the last two jump conditions

(2.9)–(2.10) in terms of the augmented variables [lu](s) and [lv](s) as follows.

Lemma 1. Let p, u, and v be the solution to the Stokes equations(1.1)–(1.2). Define

~

u¼ lu; ~v¼ lv; ~u¼ ð~u; ~vÞ: ð2:15Þ

Then we have the following jump relations for ~u and ~v o~u on ¼ ^f2þ o~u os n sin h o~u os s cos h; ð2:16Þ o~v on ¼ ^f2þ o~u os n cos h o~u os s sin h; ð2:17Þ o~u on n ¼ o~u os s : ð2:18Þ

Proof. Note that n = (cos h, sin h) and s = (sin h, cos h). Re-write the incompressibility condition [l$ Æ u] = 0 in the local coordinates, we have

o~u on cos h o~u os sin hþ o~v on sin hþ o~v os cos h¼ 0; ð2:19Þ which is o~u on cos hþ o~v on sin h¼ o~u os s : ð2:20Þ

Re-write the interface relation (2.14) in the local coordi-nates, we have o~u on sin hþ o~v on cos h¼ ^f2 o~u os n : ð2:21Þ

From the two equalities (2.20) and (2.21)above, we solve

o~u on and o~v on

to get(2.16) and (2.17). The last equality is ver-iﬁed by substituting o~u on with(2.16), and o~v on with(2.17)in the following: o~u on n ¼ o~u on cos hþ o~v on sin h¼ o~u os s :

3. The numerical algorithm

Our numerical method is based on the following theorem.

Theorem 2. Let p, u, and v be the solution to the Stokes equations (1.1)–(1.2). Let q₁ðsÞ ¼ ½~uðsÞ ¼ ½luðsÞ, q₂ðsÞ ¼ ½~vðsÞ ¼ ½lvðsÞ, and q(s) = (q1(s), q2(s)). Then ~u, ~v, p, q1(s),

q2(s) are the solution of the following augmented system of

partial differential equations: Dp¼ r g; ½p ¼ ^f1 2oq_os s; op_on ¼ ½g n þo^f2 os þ wðqÞ; ( ð3:22Þ D~u¼ px g1; ½~u ¼ q1; o~onu ¼ ^f2þoq_os n sin hoq_os s cos h; ( ð3:23Þ D~v¼ py g2; ½~v ¼ q2; o~von ¼ ^f2þoq_os n cos h oq os s sin h; ( ð3:24Þ ~ u l ¼ 0; ~v l ¼ 0; ð3:25Þ

(5)

where wðqÞ ¼ 2o 2_q_^ 1 og2 2j o ^q₂ og ; ð3:26Þ ^

q₁¼ q n, and ^q₁¼ q s in the local Cartesian coordinate system.

The proof of the theorem is straightforward from the Stokes equations(1.1)–(1.2), the jump conditions in Theo-rem 1and Lemma 1. The periodic boundary condition is used here so that we are not introducing additional bound-ary condition for the pressure.

The existence and uniqueness of the solution to the system above is the same as the original incompressible Stokes equa-tions(1.1)–(1.2). This is because if (u, v) and p are the solu-tion to the original Stokes equasolu-tions, then they are also the solution to the system inTheorem 2according to the deﬁni-tion of (u, v),Theorem 1, andLemma 1. On the other hand, if (u, v) and p are the solution to the system(3.22)–(3.25)above plus the periodic boundary condition, then they satisfy all the equations in(1.1)–(1.2)and the incompressibility condi-tion. So they are also solution to the original problem.

Notice that if we know q, then the jump conditions for the pressure are all known and we can solve the pressure independently. After the pressure is solved, we can solve the velocity from(3.23) and (3.24). The three equations with the given jump conditions can be solved using the immersed interface method[17–19]in which a fast solver is invoked with modiﬁed right hand sides at grid points near or on the interface. This observation is the basis of our new method. The compatibility condition for the two augmented variables are the two equations in(3.25)(which means that the velocity is continuous across the interface). It is also important to mention that the incompressibility condition is used to obtain the pressure Poisson equation of(3.22).

Once the augmented variables ([lu] and [lv]) and the augmented equations (the two equations in(3.25)) are cho-sen, the success of the numerical algorithm depends on how eﬃciently we can solve the augmented variables. Note that the augmented approaches have been developed for elliptic interface problems with a piecewise constant coeﬃcient

[8,20], and the fast algorithms for Poisson and biharmonic equations on irregular domains[3,10,24,23].

We assume that the domain X is a rectangle: [a, b]· [c, d], and use a uniform Cartesian grid

xi¼ a þ ihx; i¼ 0; 1; . . . ; m; hx¼ b a M ; y_j¼ c þ jhy; j¼ 0; 1; . . . ; n; hy ¼ d c N :

For simplicity of discussion, we also assume that M = N. We ﬁrst choose a set of interface points, sometimes also called the control points, {Xk} = {(Xk, Yk)}, k = 1,2, . . . ,

Nb, on the interface.2The auxiliary variable q(s) = (q1(s),

q2(s)) is deﬁned, and the augmented equations (3.25) are

discretized, at {Xk}. We use upper case letters such as Pij,

Uij, Vij, Qk, for the discrete approximations at grid points

and at the interface points, respectively. We use the bold face upper case letters without subscripts to represent the vectors formed by those discrete components.

Given an initial guess of Q at the interface points, we can approximate its ﬁrst and second order tangential deriv-atives oq_os and o_os2q2. Thus all the jump conditions in (3.22)–

(3.24)are known and we can solve (3.22)–(3.24)using the immersed interface method [17,18]. Note that the jump conditions for the pressure and the velocity are decoupled if we know Q. Since the solution depends on Q, the solu-tion can be written as P(Q), eUðQÞ and eVðQÞ.

If the computed eUðQÞ and eVðQÞ satisfy the two equa-tions in (3.25), then P(Q), eUðQÞ=l and eVðQÞ=l are an approximate solution to the original system (1.1)–(1.2). Otherwise, we use an accurate linear interpolation scheme to evaluate the residual of the two equations in (3.25)

which will be explained in Section3.2.

3.1. The discrete system of equations in the matrix–vector form

Given a discrete approximation of (q1, q2) at {Xk}, we

can solve the ﬁrst three equations (3.22)–(3.24) using the immersed interface method [18] to get an approximate solution: the pressure P(Q), the scaled velocity eUðQÞ and

e

VðQÞ. Generally the computed velocity ð eUðQÞ; eVðQÞÞ do not satisfy the two augmented equations in(3.25), that is, ðU; VÞ ¼ ð eU=l; eV=lÞ may not be continuous across the interface.

Let us assemble the discrete solution {Pij}, {Uij}, and

{Vij} together as a big vector fU whose dimension is

3MN. We denote also the vector of the discrete values of (q1, q2) at the interface points {Xk} by Q whose dimension

is 2Nb. Then the discrete solution of(3.22)–(3.24) given Q

can be written as

AfU þ BQ ¼ F1 ð3:27Þ

for some vector F1and sparse matrices A and B. It requires

solving three Poisson equations with diﬀerent force terms and jump conditions to get fU.

Once we know the solution fU given Q, we can use ð eU; eVÞ and the jump conditions oeU

on h i and oeV on h i which also depend on Q, to get ½UðQÞ ¼ ½ eUðQÞ=l and ½VðQÞ ¼ ½ eVðQÞ=l at those interface points {Xk}, 1 6 k 6 Nb. If

bothk[U(Q)]k and k[V(Q)]k are smaller than a given toler-ance, then the method has already converged and Q, eU=l,

e

V=l are the approximate solution. The interpolation scheme to get ½ eUðQÞ=l and ½ eVðQÞ=l, which will be explained in detail in the next sub-section, depends fU, Q linearly. Therefore, we can write

½UðQÞj_C ¼ ð½ eU=l; ½ eU=lÞT¼ SfU þ EQ F2; ð3:28Þ

2 _{In a front tracking method, {X}

k} is a set of points that represents the

interface, see[19], for examples in which a cubic spline is used. In a level set method, the interface points can be chosen as the orthogonal projections of irregular grid points on the interface, see[8,20], for example.

(6)

where S and E are two sparse matrices, and F2is a vector.

The matrices depend on the interpolation scheme but do not need to be actually constructed in our algorithm. We want to choose such a vector Q that the continuity condi-tion for the velocity is satisﬁed along the interface C. If we put the two matrix–vector equations(3.27) and (3.28) to-gether we get A B S E _f U Q " # ¼ F1 F2 : ð3:29Þ

Note that Q is deﬁned only on a set of interface points {Xk}

on the interface while fU is deﬁned at all grid points. The Schur complement for Q is

ðE SA1BÞQ ¼ F2 SA1F1¼ F: ð3:30Þ

After we solve(3.30)for Q, we can get fU easily. Because the dimension of Q is much smaller than that of fU, we ex-pect to get a reasonably fast algorithm for the two-phase Stokes equations.

In implementation, we use the GMRES [30] iterative method to solve(3.30). The GMRES method only requires the matrix–vector multiplication, and thus the Schur com-plement, maybe a dense matrix, of Q does not have to be explicitly constructed. We explain below how to evaluate the right hand side F of the Schur complement, and how to evaluate the matrix–vector multiplication needed in the GMRES iteration.

3.1.1. Evaluation of the right hand side of the Schur complement

First we set Q = 0 and solve the de-coupled system

(3.22)–(3.24), or (3.27) in the discrete form, to get fUð0Þ which is A1F1from(3.27). From the interpolation scheme

(3.28), we also have

½Uð0Þj_C¼ SfUð0Þ þ E0 F2¼ SfUð0Þ F2: ð3:31Þ

Note that the residual of the Schur complement for Q = 0 is

Rð0Þ ¼ ðE SA1BÞ0 F ¼ F ¼ ðF2 SA1F1Þ

¼ F2þ SfUð0Þ ¼ ½Uð0ÞjC ð3:32Þ

which gives the right hand side of the Schur complement system with an opposite sign.

3.1.2. Evaluation of the matrix–vector multiplication The matrix–vector multiplication of the Schur comple-ment system given Q is obtained from the following two steps:

Step 1: Solve the coupled system(3.22)–(3.24), or(3.27)in the discrete form, to get fUðQÞ.

Step 2: Interpolate fUðQÞ using (3.28) to get [U(Q)]jC.

Then the matrix–vector multiplication is ðE SA1_{BÞQ ¼ ½UðQÞj}

C ½Uð0ÞjC: ð3:33Þ

This is because

ðE SA1BÞQ ¼ EQ SA1BQ

¼ EQ SðA1F1 fUðQÞÞ ðfrom ð3:27ÞÞ;

¼ EQ þ SfUðQÞ F2þ F2 SA1F1

¼ ½UðQÞjCþ F ðfrom ð3:28ÞÞ;

¼ ½UðQÞjC ½Uð0ÞjC ðfrom ð3:32ÞÞ:

Now we can see that a matrix–vector multiplication is equivalent to solving the coupled system (3.22)–(3.24), or

(3.27)in the discrete form, to get fU, and using an interpo-lation scheme(3.28)to get [U(Q)]jCat the interface points.

Since we know the right hand side of the linear system of equations and the matrix–vector multiplication of the coef-ﬁcient matrix, it is straightforward to use the GMRES or other iterative methods.

3.2. The least squares interpolation scheme to compute the residual

The interpolation scheme (3.28)to evaluate ½ eU=l and ½ eV=l is crucial to the efficiency of the method. This is because the interpolation coefficients will affect the entries of the Schur complement, so its condition numbers. To reduce the number of iterations, it is important to couple the solutions on both sides of the interface using the jump conditions. Although the least squares interpolation scheme is not a new idea anymore, the details vary with problems. We explain the least squares interpolation scheme for computing(3.28)to see why we have the second matrix–vector equation in expression (3.29).

Given an approximation to the augmented variable Q, we can solve the pressure, and then the velocity ðU; VÞ ¼ ð eU=l; eV=lÞ from(3.22)–(3.24). Since

e U l " # ¼Ue þ lþ e U l ;

we need to evaluatef eUþ_{g and f e}_U_{g at all interface points}

to get the vector ½ eU=l. To explain the idea, however, we just need to explain the interpolation scheme for eUðXÞ at a point X on the interface. The interpolation scheme can be written as e UðXÞ ¼X ks1 k¼0 c_kUeiþik;jþjk þ C; ð3:34Þ where ksis the number of grid points involved in the

inter-polation scheme,ðxi; y

jÞ is the closest grid point to X, and C is a correction term. We should point it out that a one-sided interpolation scheme works poorly in the sense that the convergence speed is slow for the GMRES iteration. Below we discuss how to determine the coeﬃcients ckand

the correction term C using the information from both sides of the interface. Note that ck and C depend on X,

but for simplicity of the notation, we have omitted the dependency.

(7)

We use an un-determined coeﬃcients method to deter-mine the coeﬃcients ckby minimizing the truncation error

of (3.34)when eUi_þi

k;jþjk is substituted by the exact solu-tion ~uðxi_þi

k; yj_þj

kÞ. Using the local coordinates system centered at the point X, see (2.5), and denoting the local coordinates ofðxi_þi

k; yj_þj

kÞ as (nk, gk), we have the follow-ing from the Taylor expansion:

~ uðxiþk; yj_þkÞ ¼ ~uðn_k;g_kÞ ¼ ~uþ nk~un þ gk~u g þ 1 2n 2 k~u nn þ nkgk~ungþ 1 2g 2 k~u ggþ Oðh 3_Þ; _ð3:35Þ

where the + or sign is chosen depending on whether (nk, gk) lies on the + or side of C, ~u, ~un; . . . ; ~uggare

eval-uated at the local coordinates (0, 0), or X = (X, Y) in the original coordinates system (Fig. 2). Note that we should have used something like ~uðX ; Y Þ ¼ ~uð0; 0Þ to distinguish the two coordinate systems. However, we omit the bars and use the same notation ~uðX ; Y Þ ¼ ~uð0; 0Þ for simplicity. We carry out this expansion at all the grid points used in the interpolation scheme and plug(3.35)into(3.34). After collecting and arranging terms, we can write

~

uðXÞ a1~uþ a2~uþþ a3~un þ a4~uþn þ a5u~g þ a6~uþg

þ a7~unnþ a8~uþnnþ a9u~ggþ a10~uþggþ a11~ung

þ a12~ungþ C; ð3:36Þ

where the ai’s are given by

a1¼ X k2K ck; a2¼ X k2Kþ ck; a3¼ X k2K nkck; a4¼ X k2Kþ nkck; a5¼ X k2K gkck; a6¼ X k2Kþ gkck; a7¼ 1 2 X k2K n2_k; c_ka8¼ 1 2 X k2Kþ n2_kc_k a9¼ 1 2 X k2K g2_kck; a10¼ 1 2 X k2Kþ g2_kck; a11¼ X k2K nkgkck; a12¼ X k2Kþ nkgkck: ð3:37Þ Note that ~uþ_{¼ ~}_u_{þ q} 1 and o~onu is known from (2.16). From [18,19], we also have the following interface relations:3 ~ uþ_g ¼ ~u_g þdq1 dg; ~ uþ_nn¼ ~u_nnþ j o~u on d 2 q₁ dg2 þ ½px g1; ~ uþ_gg¼ ~u_gg j o~u on þd 2 q₁ dg2 ; ~ uþ_ng¼ ~u_ngþ jdq1 dgþ d dg o~u on ; ð3:38Þ

where j is the curvature of the interface at X. Therefore we can express all the quantities from + side in(3.36)in terms of those from side and the known quantities. Thus

(3.34), when eUiþik;jþjk is substituted by the exact solution ~ uðxiþik; yj_þj kÞ, can be written as ~ uðXÞ X k ck~uðxi_þi k; yjþjkÞ þ C ¼ a1~uþ a2~uþþ a3~un þ a4~uþn þ a5~ug þ a6~uþg þ a7~unn þ a8~uþnnþ a9~uggþ a10u~þggþ a11~ungþ a12~uþngþ C

¼ ða1þ a2Þ~uþ ða3þ a4Þ~un þ ða5þ a6Þ~ug

þ ða7þ a8Þ~unnþ ða9þ a10Þ~uggþ ða11þ a12Þ~ung

þ a2½~u þ a4½~un þ a6½~ug þ a8½~unn

þ a10½~ung þ a12½~ung þ C:

To minimize the local truncation error, we should set the following linear system of equations for the coeﬃcients ck

by matching the terms of ~u; ~u_n; . . . ; ~u_ng: a1þ a2¼ 1; a3þ a4¼ 0;

a5þ a6¼ 0; a7þ a8¼ 0;

a9þ a10¼ 0; a11þ a12¼ 0:

ð3:39Þ

The system of equations for {ck} is independent of jumps

which means we can calculate {ck} outside of the GMRES

iteration. Once we have the coeﬃcients, the correction term is C¼ ða2½~u þ a4½~un þ a6½~ug þ a8½~unn þ a10½~ung þ a12½~ungÞ ¼ a2q1 a4 o~u on a6 dq₁ dg a8 o~u on jd 2 ½~u dg2 þ ½px g1 a10 d2q₁ dg2 o~u on j a12 d½~u dg jþ d dg o~u on : ð3:40Þ We choose a neighborhood of {X} that contains more than six diﬀerent grid points so that we have an

+ ( 1 1) X Ω Ω η η − ξ ξ ,

Fig. 2. A diagram of the local coordinates used in the interpolation(3.35).

3 _{Note that q}

1now is a quantity deﬁned only on the interface, and we

can only take its derivatives along the interface which is called surface derivative in the literature.

(8)

under-determined system. In our numerical tests, we choose ks= 12, that is, we selected 12 closest grid points

to X = (X, Y) as the interpolation stencil. We use the singu-lar value decomposition (SVD) to ﬁnd the least squares solution, which also has the least l 2 norm among all the solutions. In this way, the magnitude of the coeﬃcients ckis controlled and balanced. The least squares

interpola-tion plays an important role in the stability of the algorithm.

The only trade-oﬀ of the least squares interpolation is that we have to solve an under-determined system of equa-tions. However, the size of the linear system is small and the coeﬃcients can be pre-determined before the GMRES iteration. The extra time needed in dealing with the inter-face is usually less that 5% of the total CPU time and the percentage decreases as the mesh size (h) decreases. Remark. By setting a1+ a2= 0 and a3+ a4= 1 while

keeping other equations unchanged in(3.39), we can easily get the normal derivative of the solution u. This is the method that we used in Section4for the convergence rate analysis of the computed normal derivative of the velocity. There are also a few issues about how to evaluate pxand

pyat all grid points before we can solve the velocity. Taking

px as an example, if (xi1, yj), (xi, yj), (xi+1, yj) are on the

same side of the interface, we can approximate px by the

standard central difference, ðpxÞij 1 2hðPiþ1;j Pi1;jÞ: Otherwise, we use ðpxÞij 1 hðPi;j Pi1;jÞ

if ðxi1; yjÞ and ðxi; yjÞ are on the same side; 1

hðPiþ1;j Pi;jÞ

if ðxi; yjÞ and ðxiþ1; yjÞ are on the same side;

Pij Plj ½p ½pxðxl XijÞ ðxi xlÞ otherwise: 8 > > > > > > > > > > > > < > > > > > > > > > > > > :

In the last case the interface cuts through both between xi1

and xi, and xiand xi+1. LetðXij; yjÞ be one of the two

inter-sections between the interface and the grid line y = yj. The

subscript l = i 1 or l = i + 1 is chosen such that jxl Xijj

is the smaller one. The sign in the expression above de-pends on which side of the interface the point (i, j) is on, and [px] (or [py] in the y-direction) are calculated at

ðX ij; yjÞ from ½px ¼ op on cos h op os sin h; ½py ¼ op on sin hþ op os cos h:

The analysis of such approximations to px (or py) is

dis-cussed in[19].

4. Numerical experiments and analysis

In this section, we present several numerical tests to check the order of accuracy and efficiency of the aug-mented method proposed in this paper. In the first three examples, we choose the interface as a unit circle and we know the analytic solution under different situations. It is important to mention that the reason for choosing the cir-cular interface is just for the sake of constructing the exact solutions that satisfy Eqs.(1.1)–(1.3) and the jump condi-tions(2.7)–(2.11). The circular interface is not a restriction of our method. In fact, we will also show the results for more complicated interfaces in this section.

In all the numerical tests, the interface points are chosen in such a way that the interface mesh size Ds has the same order of magnitude as the Cartesian mesh size h. In other words, the ratio Ds/h is a constant. When we reduce the grid mesh h by half during a grid reﬁnement analysis pro-cess, we also reduce the interface mesh size Ds by half.

How many interface points that are needed to represent the interface depends on the complexity of the interface. The detailed analysis about the eﬀects of the number of interface points can be found in [20]. The surface deriva-tives needed in the algorithms are computed using the cubic spline generated by the interface points {Xk}, see[19]. The

main cost of the algorithm is from the fast Poisson solver that requires O(N2log N) operations, where N is the num-ber of grid lines in the x- and y-directions. Therefore, each GMRES iteration requires about O(3N2log N) operations. All the simulations are done with double precision.

Most of the computations are done on workstations or a Notebook PC’s within a few seconds to a few minutes for stationary Stokes problems. The tolerance of GMRES iter-ation is set as = 106. The GMRES is set to re-started after Nb iterations, where Nb is the number of interface

points where Q is deﬁned. In almost all the cases, the GMRES converges in much fewer steps than Nb. The

ini-tial guess for the GMRES iteration is chosen as zero vector unless 0 is the exact solution, in which case a random vec-tor is used. Throughout this section, the computational domain is X = [2, 2] · [2, 2] unless it is stated otherwise. Example 4.1. We start with a simple example where the velocity is smooth and the pressure is discontinuous across the interface. The exact velocity and the pressure are given by u¼ yðx2_{þ y}2_1Þ; _{ðx; yÞ 2 X;} _ð4:41Þ v¼ xðx2_{þ y}2_1Þ; _{ðx; yÞ 2 X;} _ð4:42Þ p¼ 1; if x 2_{þ y}2₆_1; 0; if x2_{þ y}2_>_1: ð4:43Þ The interface is the unit circle. The viscosity is

l¼ 1; if x 2_{þ y}2₆_1; 1 2; if x 2_{þ y}2_>_1: ( ð4:44Þ

(9)

The external forcing term g is given by g₁¼ 8y; if x 2_{þ y}2₆_1; 4y; if x2_{þ y}2_>_1; ð4:45Þ g₂¼ 8x; if x 2_{þ y}2₆_1; 4x; if x2_{þ y}2_>_1; ð4:46Þ which has a ﬁnite jump across the interface. The normal and tangential force density are

^ f1¼ ½p 2 l ou on n ¼ 1; ð4:47Þ ^ f2¼ l ou on s lou os n ¼ 1; ð4:48Þ

calculated from (2.7) and (2.9), respectively.

InTable 1, we show the result of the convergence rate analysis. We use the Dirichlet boundary condition for the velocity u, and the Neumann boundary condition op_on¼ 0 for the pressure when we solve the three Poisson equations

(3.22)–(3.24). Similarly, the Dirichlet boundary condition for the velocity and Neumann boundary condition for the pressure are used in other examples when the exact solution is not periodic in this section. In[7,13], the correct boundary conditions for the pressure has been derived. The Neumann boundary conditions for the pressure is coupled with the velocity. Since the augmented approach proposed in this paper is an iterative method, we can use the updated pressure and velocity to approximate the boundary condi-tion as proposed in[7,13]. The fast Poisson solver that we use if from the Fishpack which allows Dirichlet, Neumann, and periodic boundary conditions along each side of a rect-angular domain. The errors inTable 1are measured in the maximum norm at all grid points, for example,

EuðN Þ ¼

1

2 06i;j6Nmax jUij uðxi; yjÞj þ max06i;j6NjVij vðxi; yjÞj

; ð4:49Þ where u(xi, yj) is the exact solution at (xi, yj) while Uijis the

approximate solution and so on. In all the tables in this sec-tion, N is the number of grid lines used in both x- and y-directions. The ratio

p-order¼logðEpðN Þ=Epð2N ÞÞ

log 2 ; ð4:50Þ

is an indication of the order of accuracy for the pressure. On average, we observe second order convergence for all the quantities. The last column is the number of iterations

(No.) of the GMRES method. We can see that the number of iterations remains roughly the same as we halve the mesh size h.

The errors usually do not decrease monotonically as we reﬁne the grid, see[20]. To be more precise, we use linear regression analysis to ﬁnd the approximate order of accu-racy. InFig. 3, we show the error plot on a log–log scale for the pressure and the velocity versus the grid spacing h = hx= hy, which shows that, from the slopes, the average

convergence rate of the pressure and the velocity are 1.9187 and 1.9416, respectively. The mesh size varies from N = 100 to N = 320 according to N = 100 + 5k, k = 0, 1, . . . , 44. The order of accuracy for the normal velocity ou/on is 2.1928 from the linear regression analysis. Example 4.2. In the ﬁrst example, while the pressure is discontinuous, the velocity is smooth and vanishes at the interface. Thus the exact solution to the augmented variable q = [lu] is zero which makes the problem easier to compute. In this example, we keep all the quantities u, v, p in(4.41)–(4.48) inside the unit circle unchanged, but set both pressure and the velocity outside the circle to be zero. Therefore, the periodic boundary conditions are satisﬁed.

Table 1

Numerical results and convergence analysis forExample 4.1

N Ep p-order Eu u-order Eou/on ouon-order No.

32 8.2573· 103 _6.5931_{· 10}3 _2.4605_{· 10}2 ₈

64 3.0540· 103 _1.4353 _1.7372_{· 10}3 _1.9242 _6.0768_{· 10}3 _2.0176 ₉

128 9.4747· 104 _1.6883 _3.9504_{· 10}4 _2.1367 _1.6713_{· 10}3 _1.8623 ₈

256 2.6866· 104 _1.8183 _8.2274_{· 10}5 _2.2635 _4.1920_{· 10}4 _1.9953 ₇

512 7.4314· 105 _1.8541 _2.5053_{· 10}5 _1.7155 _1.0445_{· 10}4 _2.0048 ₇

Fig. 3. Linear regression analysis of the convergence order in log–log scale for the pressure and the velocity forExample 4.1. The order of accuracy for the pressure and the velocity are 1.9187 and 1.9416, respectively.

(10)

The velocity is non-smooth and the jump in the normal velocity is not a constant along the interface. Now the normal force density component is still ^f1¼ 1, but the tangential force density component is

^ f2¼ l ou on s lou os n ¼ 2: ð4:51Þ

The external force is also adjusted to g = 0 outside of the circle. Thus, there is a ﬁnite jump in g across the interface as well. InTable 2, we show the result of the convergence analysis. Once again we observe second order accuracy on average. Only a few iterations are needed in the GMRES iteration and the number changes little as N in-creases. The results of the linear regression analysis for the pressure and the velocity are given inFig. 4which con-ﬁrms average second order accuracy for both the pressure and the velocity.

Example 4.3. In previous examples, the force density com-ponents are constants, and ½lo2^u

os2 ¼ 0. In this example, we construct the exact solutions in such a way that all the jumps and their derivatives along the interface are non-constant functions. The exact velocity and the pressure are given by u¼ y 4; if x 2_{þ y}2₆_1; y 4ðx 2_{þ y}2_Þ; _{if x}2_{þ y}2_>_1; 8 > < > : ð4:52Þ v¼ x 4ð1 x 2_Þ; _{if x}2_{þ y}2₆_1; xy 2 4 ; if x 2_{þ y}2_>_1; 8 > < > : ð4:53Þ p¼ 3 4x 3_þ3 8x y; if x2_{þ y}2₆_1; 0; if x2_{þ y}2_>_1: 8 < : ð4:54Þ

The external forcing term g is

g₁¼ 9 4x 2_þ3 8 y; if x2_{þ y}2₆_1; 2lþ_y _{if x}2_{þ y}2_>_1; 8 < : ð4:55Þ g₂¼ 3 4x 3_þ3 8x 3l 2 x; if x 2_{þ y}2₆_1; lþ 2 x; if x 2_{þ y}2_>_1; 8 > < > : ð4:56Þ

which is discontinuous across the interface. The force density components in the normal and tangential directions are ^ f1¼ 3 4 cos 3_h3 8 cos h sin h3 2½l cos 3_h_{sin h;} _ð4:57Þ ^ f2¼ 1 2l þ_þ3 4½l cos 2_{hð1 2 cos}2_hÞ; _ð4:58Þ

respectively. All the jump conditions(2.7)–(2.10)are satis-ﬁed. We use the exact Dirichlet boundary condition for the velocity and the homogeneous Neumann boundary condi-tion for the pressure.

In Table 3, we show the convergence rate analysis for different jump in l. We scale the problem such that max{l, l+} = 1, and test our results for l/l+= 10, 103, and 103. While the accuracy does depend on l/l+, the average convergence rates are about the same (second order accurate). Note that, there are two very different scales for the problems inTable 3(b) and (c). The number of iterations seems to be dependent on the ratio l/l+but not on the mesh size N. Note that the number of iterations of the GMRES method in Table 3(c) is many more than that inTable 3(b). One intuitive explanation is the follow-ing. When l+ l, the solution inside the circle is much smoother than that outside. In other words, it is more dif-ficult to resolve the solution outside than that inside. Since there are more grid points outside of the interface, there are Table 2

32 8.4430· 103 _3.4549_{· 10}3 _2.8308_{· 10}2 ₉ 64 2.8405· 103 _1.5716 _8.8800_{· 10}4 _1.9600 _6.1798_{· 10}3 _2.1956 ₁₀ 128 8.0952· 104 _1.8110 _2.2666_{· 10}4 _1.9700 _1.8260_{· 10}3 _1.7589 ₁₁ 256 2.5417· 104 _1.6713 _4.7693_{· 10}5 _2.2487 _5.3612_{· 10}4 _1.7681 ₁₂ 512 1.4086· 105 _2.1296 _1.4086_{· 10}5 _1.7595 _1.2538_{· 10}4 _2.0962 ₁₃ –3.2 –3 –2.8 –2.6 –2.4 –2.2 –2 –1.8 –1.6 –1.4 –4.4 –4.2 –4 –3.8 –3.6 –3.4 –3.2 –3 –2.8 –2.6 log 10h log 10 E p order=1.9168 u_order=2.1519

Fig. 4. Linear regression analysis forExample 4.2. The average conver-gence order for the pressure and the velocity are porder= 1.9168,

(11)

more iterations needed in this case. If we increase the radius of the interface from r = 1 to r = 1.5, then the num-ber of the GMRES iterations is down to 22 from 36 in

Table 3(c) for N = 32. This is consistent with our explana-tion above. Another explanaexplana-tion is that: In the absence of information from the outer region, the inner solution would be arbitrary up to an additive constant. That con-stant is determined by consistency with the outer solution and takes some extra eﬀort by the solver to pin down. This problem does not arise for l+ lbecause in this case the outer solution is (relatively) small so that it is not diﬃcult to approach the correct constant. In any case, the issue is less sensitive because of the scaling.

Example 4.4. In the examples above, the interface is chosen as a circle so that we can construct exact solutions to compare with the computed solution. In this example, we test our method for a complex geometry. The interface is given by

q¼ 0:35 þ 0:1 sinð6hÞ; 0 6 h 6 2p; ð4:59Þ

in polar coordinates in the rectangular domain [1, 1] · [1, 1], seeFig. 5. The curvature of the interface varies both in the magnitude and the sign, and the source term is g = 0. The force density is ^f1¼ 0:1j and ^f2¼ 0:1, respectively,

where j is the curvature. Periodic boundary conditions are used for all the variables. Since the analytic solution is not available, we compare the computed solution with the solu-tion that is obtained from the ﬁnest grid, 512 by 512. In

Table 4, we show the convergence rate analysis for the pressure and the velocity. We also listed the number of iterations, the CPU time in seconds, and the ratio of two consecutive errors. For a second order method, the ratio is between 4 and 5. The justiﬁcation of such an analysis is given in [19]. The interface points are initially chosen

from (4.59) with equally spaced Dh. Then they are re-distributed with equally spaced arc-length by the cubic spline package developed in[19].

Table 3

(a) l= 1, l+= 0.1 32 6.8928· 102 _4.6299_{· 10}2 _2.0170_{· 10}2 ₁₇ 64 5.6851· 103 _3.5998 _3.4079_{· 10}3 _3.7640 _3.8114_{· 10}3 _2.4038 ₁₇ 128 2.2966· 103 _1.3077 _1.2068_{· 10}3 _1.4977 _8.7465_{· 10}4 _2.1235 ₁₇ 256 5.4715· 104 _2.0695 _2.6908_{· 10}4 _2.1651 _2.2514_{· 10}4 _1.9579 ₁₄ 512 1.5365· 104 _1.8323 _6.4921_{· 10}5 _2.0513 _5.6411_{· 10}5 _1.9968 ₁₄ (b) l_{= 0.001, l}+_{= 1} 32 1.3803· 102 8.1811· 101 13.323 12 64 4.1261· 103 1.7421 2.2177· 101 2.0017 2.8062 2.2371 13 128 1.0414· 103 1.9863 6.2257· 102 1.8328 7.7665· 101 1.8533 11 256 3.5892· 104 1.5368 1.4046· 102 2.1481 1.8162· 101 2.0969 11 512 7.0865· 105 _2.3405 _2.8175_{· 10}3 _2.3177 _4.7867_{· 10}2 _1.9293 ₉ (c) l= 1, l+= 0.001 32 0.6950 42.0260 1.4587 36 64 1.4356· 102 _2.0017 _9.4294_{· 10}1 _5.4780 _2.6627_{· 10}1 _2.4537 ₃₀ 128 6.5307· 103 _1.1363 _3.1469_{· 10}1 _1.5832 _5.9376_{· 10}2 _2.1649 ₃₁ 256 1.1757· 103 _2.4737 _4.6464_{· 10}2 _2.7598 _1.4749_{· 10}2 _2.0093 ₃₆ 512 3.0160· 104 _1.8260 _1.1697_{· 10}3 _1.9900 _3.4355_{· 10}3 _2.1020 ₃₅

The CPU time for (c) are 1.5527, 2.8925, 6.1523, 17.275, and 44.248 s, respectively, for the 5 diﬀerent grids.

Fig. 5. The interface forExample 4.4. Table 4

Numerical results and convergence analysis forExample 4.4with l= 1 and l+= 0.01

N Ep p-ratio Eu u-ratio CPU (s) No.

32 16.149 6.6068· 101 _3.5957 ₄₈

64 2.8549 5.6566 9.4029· 102 _7.0263 _8.9805 ₄₁

128 6.9682· 101 _4.0970 _2.5461_{· 10}2 _3.6931 _21.986 ₃₉

(12)

Compared with previous examples, the number of GMRES iterations is larger. This is not surprising since the curvature of the interface is large. Nevertheless, the expected second order accuracy on average is still observed.

Example 4.5. As a ﬁnal example, we show the motion and deformation of a liquid drop in simple shear ﬂow in two space dimensions. We refer the readers to [15] and the references therein for detailed description of the problem and physical explanations.

Fig. 6. Contour plot of stream-function of the shear ﬂow with various viscosity ratio k and Ca. In all of cases Ca = 1.3898, c = 0.4605 except for (d) in which Ca = 0.4695, c = 1.4605. (a) k = 0.04 and t = 0. (b) k = 0.04 and t = 0.75. (c) k = 0.08 and t = 0.75. (d) k = 1 and t = 0.75. (e) k = 6.4 and t = 0. (f) k = 6.4 and t = 0.75. The drop is the solid-dot line in the plots.

(13)

The set-up of the problem is as follows. A circular viscous drop with viscosity l2is placed in a shear ﬂow with

viscosity l1. The velocity of the shear ﬂow without drop is

u = (Gy, 0), where G is the shear rate. The computational domain is X = [2, 2] · [2, 2]. The Dirichlet boundary condition uoX= (Gy, 0) is used along oX; and Neumann

boundary conditionop/on = 0 is used for the pressure along oX. The capillary number is deﬁned as Ca = Gr0l1/c, where

r0is the radius of the circular drop, c is the surface tension.

The source strength is ^f1¼ f n ¼ cj, where j is the

curva-ture. The ratio of the viscosity is k = l2/l1. A front tracking

method is used to evolve the interface. The crucial para-meters are the capillary number Ca and the viscosity ratio k. The long-term behavior of the interface has different con-figurations depending on the values of Ca and k. Typical long-term configurations include internal circulations and an asymptotic equilibrium state as explained below.

InFig. 6we present a sequence of contour plots of the stream function of the velocity for the two-phase Stokes ﬂow in comparison with the results in[15]. In our simula-tion, the initial drop is a circle with r0= 0.75 and the shear

rate is G = 0.9143. The parameters are chosen to mimic the set-up in[15]. We try to re-produce the results there except now it is in two space dimensions. Our results agree with those presented in [15] qualitatively. A quantitative com-parison is not available. The stream-function is obtained by solving the Poisson equation Dw = vx uy with the

boundary condition determined from the shear flow. When k is small, we observe two regions of recirculating fluid within the drop. Generally, for smaller k, we observe more deformation but larger angle (less tilting) between the longer axis of the drop and the x-axis. As k gets larger, we see less deformation but smaller angle (more tilting) between the longer axis of the drop and the x-axis. All these agree with the results and analysis discussed in[15]. Note that for various parameters, internal circulations do not always occur. For instance, for circular drop at t = 0 with k = 1, the solution of the pressure is piecewise con-stant and the solution of the velocity is simply the shear flow itself: u = (Gy, 0). While we do observe internal circu-lation for Ca = 0.4695 and c = 1.4605 as presented in

Fig. 6(d), this is not true for Ca = 1.3898 and c = 0.4605. Also at t = 0.75, for all the simulations inFig. 6we have observed that the drop has reached its steady state. 5. Conclusion

In this paper, a new second order accurate finite differ-ence method has been developed for incompressible sta-tionary Stokes equations with a discontinuous viscosity in which the jump conditions for the pressure and the velocity are coupled together. The idea is to introduce two augmented variables that are only defined along the interface so that the jump conditions can be decoupled. The GMRES iterative method then is used to solve the Schur complement system for the augmented variables. The main cost in one step of the GMRES iteration is three

calls to a fast Poisson solver. Numerical examples demon-strate the eﬃciency and accuracy of the method. The idea should be applicable to other interface problems with cou-pled jump conditions. One of the remaining open questions is how to use a preconditioning technique for the Schur complement system.

Acknowledgments

The ﬁrst author is partially supported by NSF grants DMS-0201094 and DMS-0412654. The ﬁrst and second author are supported by an ARO under grant number 43751-MA. The third author is partially supported by the National Science Council of Taiwan under research grant NSC-92-2115-M-009-012. We would like to thank Drs. X.-B. Lin and S. Lubkin of North Carolina State Univer-sity for useful discussions. We would like also to thank the referees for their useful suggestions.

References

[1] Almgren A, Bell J, Colella P, Marthaler T. A Cartesian grid projection method for the incompressible Euler equations in complex geometries. SIAM J Sci Comput 1997;18:1289–309.

[2] Biros G, Ying L, Denis Z. An embedded boundary integral solver for the Stokes equations. J Comput Phys 2004;196:591–626.

[3] Chen G, Li Z, Lin P. A fast ﬁnite diﬀerence method for biharmonic equations on irregular domains. CRSC-TR04-09, North Carolina State University; 2004.

[4] Cortez R. The method of regularized Stokeslets. SIAM J Sci Comput 2001;23:1204–25.

[5] Fogelson AL, Peskin CS. Numerical solution of the three dimensional Stokes equations in the presence of suspended particles. In: Proceed-ings of the SIAM conference on multi-phase ﬂow. SIAM, June 1986. [6] Greengard L, Kropinski MC, Mayo A. Integral equation methods for Stokes ﬂow and isotropic elasticity in the plane. J Comput Phys 1996; 125:403–14.

[7] Guermond JL, Shen J. A new class of truly consistent splitting schemes for incompressible ﬂows. J Comput Phys 2003;192:262–76. [8] Hou T, Li Z, Osher S, Zhao H. A hybrid method for moving interface

problems with application to the Hele–Shaw ﬂow. J Comput Phys 1997;134:236–52.

[9] Hou TY, Stredie VG, Wu TY. A 3D numerical method for studying vortex formation behind a moving plate. Commun Comput Phys 2006;1:207–28.

[10] Hunter J, Li Z, Zhao H. Autophobic spreading of drops. J Comput Phys 2002;183:335–66.

[11] Ito K, Li Z. Interface relations for Stokes equations with discontin-uous viscosity and singular sources. Appl Math Lett 2006;19(3): 229–34.

[12] Johansen H, Colella P. A Cartesian grid embedded boundary method for Poisson equations on irregular domains. J Comput Phys 1998;147: 60–85.

[13] Johnston H, Liu J. Accurate, stable and eﬃcient Navier–Stokes solvers based on explicit treatment of the pressure term. J Comput Phys 2004;188:221–59.

[14] Kang M, Fedkiw R, Liu X. A boundary condition capturing method for multiphase incompressible ﬂow. J Sci Comput 2000;15:323–60. [15] Kenedy MR, Pozrikidis C, Skalak R. Motion and deformation of

liquid drops and the rheology of dilute emulsions in shear ﬂow. Comput Fluids 1994;23:251–78.

[16] Lai M-C, Li Z. A remark on jump conditions for the three-dimensional Navier–Stokes equations involving an immersed moving membrane. Appl Math Lett 2001;14:149–54.

(14)

[17] LeVeque RJ, Li Z. The immersed interface method for elliptic equations with discontinuous coeﬃcients and singular sources. SIAM J Numer Anal 1994;31:1019–44.

[18] LeVeque RJ, Li Z. Immersed interface method for Stokes ﬂow with elastic boundaries or surface tension. SIAM J Sci Comput 1997;18: 709–35.

[19] Li Z. The Immersed interface method — a numerical approach for partial diﬀerential equations with interfaces. PhD thesis, University of Washington; 1994.

[20] Li Z. A fast iterative algorithm for elliptic interface problems. SIAM J Numer Anal 1998;35:230–54.

[21] Li Z. An overview of the immersed interface method and its applications. Taiwan J Math 2003;7:1–49.

[22] Li Z, Wan X, Ito K, Lubkin S. An augmented pressure boundary condition for a Stokes ﬂow with a non-slip boundary condition. Commun Comput Phys 2006;1(5):874–85.

[23] Li Z, Wang C. A fast ﬁnite diﬀerence method for solving Navier– Stokes equations on irregular domains. J Commu Math Sci 2003;1: 180–96.

[24] Li Z, Zhao H, Gao H. A numerical study of electro-migration voiding by evolving level set functions on a ﬁxed Cartesian grid. J Comput Phys 1999;152:281–304.

[25] Lubkin SR, Li Z. Force and deformation on branching rudiments: cleaving between hypotheses. Biomech Model Mechanobiol 2002;1: 5–16.

[26] Mayo A, Peskin CS. An implicit numerical method for ﬂuid dynamics problems with immersed elastic boundaries. Contemp Math 1991;141: 261–77.

[27] Osher S, Fedkiw R. Level set methods and dynamic implicit surfaces. New York: Springer; 2002.

[28] Peskin CS. Numerical analysis of blood ﬂow in the heart. J Comput Phys 1977;25:220–52.

[29] Peskin CS. Lectures on mathematical aspects of physiology. Lect Appl Math 1981;19:69–107.

[30] Saad Y. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J Sci Stat Comput 1986; 7:856–69.

[31] Sethian J, Wiegmann A. Structural boundary design via level set and immersed interface methods. J Comput Phys 2000;163:489–528. [32] Sethian JA. Level set methods and fast marching methods. 2nd

ed. Cambridge University Press; 1999.

[33] Sussman M, Almgren A, Bell JB, Colella P, H Howell L, Welcome ML. An adaptive level set approach for incompressible two-phase ﬂows. J Comput Phys 1999;148:81–124.

[34] Sussman M, Smereka P, Osher S. A level set approach for computing solutions to incompressible two-phase ﬂow. J Comput Phys 1994;114: 146–59.

[35] Teman R. Navier–Stokes equations. New York: North-Holland; 1977.

[36] Tornberg A-K, Engquist B. Segment projection method for interface tracking. Commu Pure App Math 2003;56:47–79.

[37] Tu C, Peskin CS. Stability and instability in the computation of ﬂows with moving immersed boundaries: a comparison of three methods. SIAM J Sci Stat Comput 1992;13:1361–76.

[38] Wiegmann A. The explicit jump immersed interface method and interface problems for diﬀerential equations. PhD thesis, University of Washington; 1998.

[39] Wiegmann A, Bube K. The immersed interface method for nonlinear diﬀerential equations with discontinuous coeﬃcients and singular sources. SIAM J Numer Anal 1998;35:177–200.

[40] Wiegmann A, Bube K. The explicit jump immersed interface method: ﬁnite diﬀerence methods for PDEs with piecewise smooth solutions. SIAM J Numer Anal 2000;37:177–200.