A well-conditioned augmented system for solving Navier-Stokes equations in irregular domains

(1)

A well-conditioned augmented system for solving Navier–Stokes

equations in irregular domains

Kazufumi Ito

a

, Ming-Chih Lai

b

, Zhilin Li

a,* a

Center for Research in Scientiﬁc Computation and Department of Mathematics, North Carolina State University, Raleigh, NC 27695-8205, United States b

Department of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsinchu 30050, Taiwan

a r t i c l e

i n f o

Article history:

Received 31 October 2007

Received in revised form 20 October 2008 Accepted 17 December 2008

Available online 31 December 2008

AMS classiﬁcation: 65M06 65M12 76T05 Keywords: Navier–Stokes equations Embedding technique Immersed interface method Irregular domain Augmented system Projection method Level set representation

a b s t r a c t

An augmented method based on a Cartesian grid is proposed for the incompressible Navier–Stokes equations in irregular domains. The irregular domain is embedded into a rectangular one so that a fast Poisson solver can be utilized in the projection method. Unlike several methods suggested in the literature that set the force strengths as unknowns, which often results in an ill-conditioned linear system, we set the jump in the normal derivative of the velocity as the augmented variable. The new approach improves the condition number of the system for the augmented variable signiﬁcantly. Using the immersed interface method, we are able to achieve second order accuracy for the velocity. Numerical results and comparisons to benchmark tests are given to validate the new method. A lid-driven cavity ﬂow with multiple obstacles and different geometries are also presented.

1. Introduction

In this paper, we consider the following incompressible Navier–Stokes equations in an irregular domain RnX:

q

@u @tþ ðu

r

Þu þ

r

p ¼

l

D

u þ G; x 2 R n

X

; ð1:1Þ

r

u ¼ 0; x 2 R n

X

; ð1:2Þ uj@R¼ u@R; uj@X¼ u@X; BC; ð1:3Þ uðx; 0Þ ¼ u0; IC; ð1:4Þ

where Gðx; y; tÞ is an external forcing term, R is a rectangular domain, andXis a set of inclusions (obstacles), seeFig. 1for an illustration. We assume that

q

¼ 1 in this paper.

Since the inclusionXis irregular inside a regular computational domain R, one interesting numerical issue is how to im-pose the prescribed velocity condition at the immersed boundary @X, particularly, when the boundary is moving. While

* Corresponding author.

E-mail address:zhilin@math.ncsu.edu(Z. Li).

Contents lists available atScienceDirect

Journal of Computational Physics

j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / j c p

(2)

there are a number of numerical methods for solving Navier–Stokes equations in irregular domains, we are interested in Cartesian grid based ﬁnite difference methods. One of the main motivations using Cartesian grid methods is to avoid mesh generation for moving boundary problems. Cartesian grid methods also enable us to use some commonly used fast solvers such as FFT for solving Poisson/Helmholtz equations.

Among Cartesian grid based ﬁnite difference methods, one common approach is to use a streamfunction–vorticity formu-lation, see for example,[3,20,16]. In these methods, the domain is embedded into a rectangular one so that the immersed interface method can be used to get second order accuracy for the vorticity. But this formulation is often restricted to two-dimensional problems only.

A second type of finite difference methods is based on the projection method for small to medium Reynolds number flow, see for example,[9,21]. This type of method has been modified to compute large eddy simulations as well, see for example,

[19]for a brief review and the references therein. In these methods, the boundary is treated as an immersed interface that exerts forces to the surrounding fluid. The force density is chosen so that the boundary condition is satisfied. The problem in finding such force strengths is an inverse problem and it is ill-conditioned. Different strategies have been developed to reg-ularize the problem. In[22], the method was coupled with an optimal control; in[10], singular value decomposition is used. In most of these methods, the boundary is represented by a set of control points (or a cubic spline). The condition number of the system of equations is often very large and depends on the number of control points. The ill-conditioned system may affects the accuracy of the computed solution in an adversely way.

Other related work can be found in[5,17]. The method in[5]used a direct discretization for a similar but different problem. Advantages of the method proposed in[5]are its simplicity and symmetry of the system of equations on the entire domain. The PCG method was applied to the entire system with OðN2

Þ unknowns at every time step. One of disadvantages is that a fast Pois-son solver can not be used. The stair-case approximation of the boundary proposed in[5]is less accurate than our approach that takes into account of the curvature of the boundary. In[18], the authors proposed a Cartesian method to simulate the incom-pressible flows around stationary or moving immersed boundaries. The finite difference discretization near the interface or boundary is modified so that an algebraic multigrid solver can be applied. Again, a fast Poisson solver based on the FFT can not be applied there. In an earlier work of Ye et al.[23], the authors used a finite volume method for a similar problem.

In order to overcome the difﬁculty, many attempts have been made to utilize a Helmholtz/Poisson solver on irregular do-mains based on a second order augmented immersed interface method in[6,12,15]. In an augmented immersed interface method, an augmented variable (often is one dimension lower than that of the solution) is introduced so that the IIM can be easily implemented. The solution then is a functional of the augmented variable and should satisfy the boundary or inter-face condition, which is the augmented equation. In this paper we develop such a method for the incompressible Navier– Stokes equations in irregular domains.

By investigating the source of the stiffness, we have found that the stiffness is from the force in the normal direction that corresponds to the jump in the pressure. Based on the fact that the pressure in Navier–Stokes equations is not a free variable, we propose a new method that sets the jump in the normal derivative of the velocity as augmented variables. This new method improves the condition of the Schur complement system for the augmented variables signiﬁcantly compared to the method that sets the force density as the augmented variable.

In an augmented method, it is relatively easy to choose and test different augmented variables with a few modiﬁcations, see[11,12]. In[13], the jump in the velocity was chosen as the augmented variable for the two-phase incompressible Stokes equations. We choose the jump in the normal derivative of the velocity for the Navier–Stokes equations because it is easier to deal with the nonlinear term since the velocity is continuous.

We use a zero level set of a Lipschitz continuous function to represent the boundary which is more flexible for multi-con-nected domains. At each time step, the linear system of equations for the augmented variable can be solved by the GMRES method, or by LU decomposition if we form the coefficient matrix at the first time step. The latter approach is especially effi-cient for stationary boundary problems and long time runs.

n

R

n

Ω

(3)

The rest of the paper is organized as follows. In Section2, we present an augmented method for solving(1.1)–(1.4)within the framework of the projection method. The novelty of our approach and the implementation issues are described in detail. In Sec-tion3, numerical results for non-trivial examples are presented and analyzed. Conclusions will be given in Section4. 2. The numerical algorithm

Our numerical method is based on the projection method for solving the incompressible Navier–Stokes equations. There are several versions of the projection method for solving the incompressible Navier–Stokes equations, see for example,

[1,7,8]and many others. The projection method that we used in our method is the one described in[2]which is based on the pressure increment formulation of[1,7].

As mentioned before, we assume that the domain R is a rectangle ½a; b ½c; d with holesX. The spatial spacing is chosen as hx¼ ðb aÞ=M, hy¼ ðd cÞ=N, where M and N are the number of grid points in the x and y-directions, respectively. We use

a standard uniform Cartesian grid for simplicity. The time integration from tk_{to t}kþ1_{can be written as:}

u_uk

D

t þ ðu

r

uÞ kþ1 2_¼

r

p k1 2þl 2ð

D

uþ

D

u k_{Þ þ G}kþ1 2 if x 2 R n

_X

l 2ð

D

u _þ

_D

_uk_Þ _{if x 2}

_X

( ð2:5Þ u_j @R¼ u@Rðx@R;tkþ1Þ; ð2:6Þ ½u @X¼ 0; @u @n @X ¼ qkþ1_; _u_j @X¼ u@Xðx@X;tkþ1Þ; ð2:7Þ

D

/kþ1¼ru Dt ; x 2 R; @/kþ1 @n _@R¼ 0; ½/kþ1@X¼ 0; @/kþ1 @n h i @X¼ 0; 8 < : ð2:8Þ ukþ1_{¼ u}

_D

_t;

_r

_/kþ1_; _{x 2 R;} _ð2:9Þ

where ðu ruÞkþ1

2_{is approximated by} ðu

r

uÞkþ1 2_¼3 2ðu k

_r

_Þuk1 2ðu k1

_r

_Þuk1_: _ð2:10Þ

The solution u_{above is a functional of the augmented variable q}kþ1_{which should be determined by imposing the boundary}

condition u_j

@X¼ u@Xðx@X;tkþ1Þ: ð2:11Þ

The Eqs.(2.5)–(2.11)now become a complete system for ðu_;_qkþ1_;_/kþ1_;_ukþ1_{Þ. Once we have solved this system, then the}

pressure is determined from

pkþ1=2_{¼ p}k1=2_{þ /}kþ1_; _{x 2 R n}

_X

_: _ð2:12Þ

Table 1

A grid refinement analysis against the exact solution at final time T ¼ 5. kEuk1is the maximal error in the velocity, kEpk1is the maximal error in the pressure, cond1is the condition number of the coefficient matrix for the unknown jump in the normal derivative of the velocity, while cond2is the condition number of the coefficient matrix when we set the normal and tangential force strengths as the unknowns.

(a) Dirichlet BC on all sides,l¼ 0:5.

M N kEuk1 r1 kEpk1 r2 cond1 cond2 CPUðsÞ

64 32 4.8913 103 _{3.3494 10}2 _17.686 _{4.4674 10}3 _0:8656 128 64 1.1499 103 _2.0887 _{8.3891 10}3 _1.9973 _19.206 _{1.0545 10}4 5.2305 256 128 2.7371 104 2.0708 5.8427 103 0.5219 28.046 8.5386 105 42.072 512 256 6.6709 105 2.0367 2.5283 103 1.2085 46.338 1.0199 107 281.12 (b) Dirichlet BC on all sides,l¼ 0:005.

M N kEuk1 r1 kEpk1 r2 cond1 cond2 64 32 1.0918 101 _{3.1598 10}2 _58.888 _{6.4600 10}4 128 64 1.7605 102 _2.6326 _{6.5403 10}3 _2.2724 _19.356 _{1.1110 10}6 256 128 2.9241 103 _2.5899 _{1.2525 10}3 _2.3845 _13.447 _{3.1600 10}7 512 256 3.4913 104 3.0622 2.7237 104 2.2012 9.7787 1.006 1010

(c) Dirichlet BC for u on x ¼ a, for v on x ¼ a, y ¼ c, and y ¼ d. Neumann BC for u on x ¼ c, y ¼ c, and y ¼ d, forvon x ¼ b,l¼ 0:005. The condition numbers are almost the same as in table (b) above so we did not list them.

M N kEuk1 r1 kEpk1 r2 64 32 8.8951 102 _{1.3054 10}1 128 64 1.7636 102 _2.3344 _{3.4977 10}2 _1.9000 256 128 2.9041 103 _2.6024 _{1.4110 10}2 _1.3097 512 256 3.4473 104 3.0745 6.6851 103 1.0785

(4)

There are two novel ideas that are signiﬁcantly different from other methods in the literature. First, we set the jump in the normal derivative of the velocity as the augmented variable so that the system(2.5)–(2.7)is complete for ðu_;_qkþ1_{Þ. Here,}

qkþ1_{represents the augmented variable, and the augmented equation is u}_j

@X¼ u@Xðx@X;tkþ1Þ. Secondly, we apply the aug-mented approach for the coupled system of ðu_;_qkþ1_{Þ so that a fast Helmholtz solver can be applied when we use GMRES to}

solve the Schur-complement system for the augmented variable qkþ1_{, see Section}_2.3_{for more details.}

We also have tested the method by embedding the Navier–Stokes equations to the entire domain and selecting the nor-mal and tangential forces along the boundary as augmented variables as described in several papers in the literature. The resulting Schur-complement system is severely ill-conditioned, seeTable 1for an illustration. But our selection of the aug-mented variable ½@u

@n@X¼ qkþ1results in a very well-conditioned system. Since we are only interested in the solution in the region of R nX, the other terms such as ðu ruÞkþ1=2_and_r_pk1=2_{can be treated as a forcing term in the elliptic equations that}

do not need to be extended. These quantities at irregular grid points where the boundary @Xcuts through the standard cen-tral 5-point stencil, can be approximated by a one-sided interpolation scheme.

Note that Eq.(2.8)is also solved on the entire rectangular domain R so that the same fast Poisson solver can be used. Note also that, the above algorithm is equivalent to the following formulation without introducing any augmented variable in the domain R nX: u_uk

D

t þ ðu

r

uÞ kþ1 2_¼

_r

_pk12_þ

l

2ð

D

u _þ

_D

_uk Þ þ Gkþ12_; _{x 2 R n}

_X

u_j @R¼ u@Rðx@R;tkþ1Þ; uj@X¼ u@Xðx@X;tkþ1Þ

D

/kþ1¼

r

u

D

t ; x 2 R n

X

; @/kþ1 @n @R ¼ 0; ukþ1_{¼ u}

_D

_t

_r

_/kþ1_; _{x 2 R n}

_X

_;

r

ukþ1_{¼ 0;} _{x 2 R n}

_X

_:

Remark 1. The projection method that we currently used is a version based on the pressure increment formulation[1,2,7]in which the no-slip boundary conditions at the embedded boundary @Xand the physical boundary @R are imposed for the intermediate velocity u _{instead of u}kþ1_{. This kind of simplification is made to be consistent with the efficiency of the} projection method of Navier–Stokes solver and has also been used in[21]. Since we are interested in the solution outside the embedded boundary, (we disregard the computed solution inside the boundary), all we need to care is whether the solution satisfies the no-slip boundary condition in a suitable accuracy comparable with the fluid solver accuracy. As discussed in[4], the spurious slip velocity error of ukþ1_{at the boundary is of second order. One may think of such error could} cause the possible fluid penetration (leakage) into or out of the embedded boundary. In fact, we have also implemented the following alternative approach by setting

ukþ1_j

@X¼ u@Xðx@X;tkþ1Þ

as the constraint (the augmented equation) instead of using u_j

@X¼ u@Xðx@X;tkþ1Þ:

In such case, the entire system for solving ðukþ1_;_qkþ1_;_pkþ1_{Þ with the augmented variable ½}@u

@n@X¼ qkþ1are all coupled to-gether. This approach is more expensive since each iteration requires to solve one more Poisson equation for /kþ1. However, the results are almost the same. From the numerical experiments that we conduct in this paper, the effect of ﬂuid penetra-tion in present simulapenetra-tions caused by the enforcement of no-slip condipenetra-tion for u_{instead of u}kþ1_{is not that signiﬁcant. As a}

matter of fact, for the simulation of ﬂow past a circular cylinder in Section3, our present method obtains comparable results with others by comparing the drag, lift and Strouhal numbers.

2.1. Some implementation details

The method can be implemented with different representations of the boundary @X. If a front-tracking method is used, then the augmented variable qkþ1_{¼ ½@u}_{=@n is deﬁned at a set of control points (Lagrangian markers) as discussed in}

[21,10,22]. In our present scheme, the boundary @Xis represented by a zero level set of a Lipschitz continuous function (of-ten an approximated signed distance function). The augmented variable is deﬁned at those orthogonal projections of irreg-ular grid points1_{located outside the domain}_X_{. The level set representation is more ﬂexible to deal with multi-connected} domains, or moving objects with possible topological changes (merging and splitting).

1

(5)

2.2. Solving an elliptic interface problem with singular sources

As described before, our scheme consists of solving two generalized Helmholtz equations for the intermediate velocity u

in(2.5), and one Poisson equation for the pressure increment /kþ1in(2.8)with a given jump condition in the normal deriv-ative of the velocity. The details for solving Helmholtz/Poisson equations with jump conditions in the solution and its normal derivative (or the elliptic PDEs with singular sources) can be found in[11,14,12]. Here, we just give a brief sketch on how we solve such problems in an efﬁcient and accurate way. Without loss of generality, we consider the following generalized Helmholtz equation wxxþ wyy kw ¼ f ; ðx; yÞ 2 R; ½w@X¼ 0; @w @n @X ¼ q: ð2:13Þ

In our application, we have k ¼ 2=ð

l

DtÞ for the prediction step, and k ¼ 0 for /kþ1_in_(2.8)_{. Since k is a constant, a fast Poisson}

solver can be used.

The ﬁnite difference discretization using the immersed interface method can be simply written as

W

iþ1;jþ

W

i1;j 2

W

i;j

h2x

þ

W

i;jþ1þ

W

i;j1 2

W

i;j h2y

k

W

ij¼ f ðxi;yjÞ þ Cij; ð2:14Þ

where the correction term Cijis zero at regular grid points where the boundary @Xdoes not cut through the standard

cen-tered 5-point stencil. The correction term Cijat those irregular grid points can be determined in a dimension by dimension

fashion, see[14]for the formula of Cij.

2.3. The augmented method

The procedure of an augmented approach varies with different applications and can be found in[12,13]. Here, we give a brief description for our particular problem in this paper. From time tk_{to t}kþ1_{, given the jump ½@u}_{=@n ¼ q}kþ1_{, or Q}kþ1_{in the}

discrete form, the discrete solution U_{for u}_{satisﬁes a linear system}

AU

þ BQkþ1¼ F1; ð2:15Þ

where the matrix A corresponds to the prediction step (2.5)–(2.7), and the vector BQkþ1_{corresponds to the correction term}

Ci;j in (2.14) due to the jump condition ½@u

@n@X¼ qkþ1. The boundary condition (2.11) (or the augmented equation) u_j

@X¼ u kþ1

@X is discretized via a least squares interpolation, and can be written in the discrete form

CU_{þ DQ}kþ1_{¼ U}kþ1

@X : ð2:16Þ

If we put those two matrix–vector Eqs.(2.15) and (2.16)together, we get

A B C D _U Qkþ1 ¼ F1 Ukþ1 @X : ð2:17Þ

In practice, we do not need to form the matrices A, B, C, and D. Note that the dimension of Qkþ1_{is Oð2NÞ (assuming M N),}

which is much smaller than that of U_(Oð2N2

Þ).

Eliminating U_{from Eq.}_(2.17)_{, we get the Schur-complement system for Q}kþ1

; ðD CA1BÞQkþ1¼ Ukþ1@X CA 1 F1¼ def F2; or EQkþ1¼ F2: ð2:18Þ

Note that E is not symmetric in general. We can either use a direct method or an iterative method to solve the linear system

(2.18)for Qkþ1depending on different situations. If the boundary @Xis stationary and the time step is fixed, it is more effi-cient to use a direct method (say, Gaussian elimination with partial pivoting) since the coeffieffi-cient matrix E is a constant ma-trix. Thus, we can form the coefficient matrix E explicitly and apply the LU decomposition to it directly. This is also advantageous for long time runs since the costly LU decomposition needs to be done just one time in the beginning. For a moving boundary problem, the matrix E is not a constant so it is more efficient to use an iterative method. In this case, it is recommended to use GMRES iteration.

Three Helmholtz/Poisson solvers are needed to evaluate EQkþ1or A1F1. First, the prediction step for urequires a fast

Helmholtz solver on a rectangular domain for each velocity component. Secondly, the projection step requires a fast Poisson solver for /kþ1_{given u}_{. For stationary boundary @}_X_{and ﬁxed time step}_D_{t, all the matrices are constant matrices that}

de-pend on the fast Helmholtz/Poisson solver, and the interpolation scheme for the boundary condition.

Since we do not form the matrices A, B, C, and D explicitly, one question is how to use an iterative or direct method. This has been explained in detail in our recent work[12,13]. First, we set Qkþ1

¼ 0 and then solve the Navier–Stokes equations. The residual of the linear system(2.18)(or the difference between the exact and the computed boundary condition), is actu-ally the right hand side of the Schur complement with an opposite sign. Next, we explain how to ﬁnd the matrix–vector mul-tiplication given Q , a guess of Qkþ1_{. This again involves only two steps: (1) solving}_(2.15)_{by given Q , to get U}kþ1;

(6)

we use the superscript ðk þ 1; Þ to indicate that the solution is not the ﬁnal solution at time tkþ1 _{yet; (2) interpolating}

Ukþ1;

ðQ Þ at @Xvia the least squares interpolation. Once we know the matrix–vector multiplication, we can apply the GMRES or other iterative method easily.

Also, we can form the coefﬁcient matrix E of the Schur-complement by setting Qkþ1;¼ e‘, the ‘th unity vector, l ¼ 1; 2; . . ..

For a stationary interface and a ﬁxed time step, this is needed just initially.

The idea of the augmented method has been explained in[11,12]for elliptic interface problems with piecewise constant coefﬁcient, and Poisson equations on irregular domains. The Fortran code is available to the public either by request or by anonymous ftp.

3. Numerical examples

In this section, we present some numerical results for solving the velocity and pressure in(1.1)–(1.4)up to some ﬁxed time T. All the computations were performed at the North Carolina State University using notebook or desktop computers. Most simulations are done within minutes to a couple of hours depending on the mesh, the geometric properties of @X, and the ﬁnal time T.

Example 1. Validation of the method against an exact solution We ﬁrst consider an example with the exact solutions as

uðx; y; tÞ ¼ sinðtÞ y r 2y if r P 1=2 0; otherwise; ( ð3:19Þ

v

ðx; y; tÞ ¼ sinðtÞ x rþ 2x if r P 1=2 0; otherwise; ( ð3:20Þ pðx; y; tÞ ¼ cosð

p

xÞ cosð

p

yÞ if r > 1=2

0; otherwise;

ð3:21Þ where r ¼pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffix2_{þ y}2_{. The obstacle is a circular disk whose boundary is a circle r ¼ 1=2. The solution domain is}

R ¼ ½2; 2 ½1; 1. The source term G is derived directly from the exact solution.

InTable 1, we show the grid reﬁnement analysis for different viscosity and boundary conditions. Since we are interested in the computed solutions in the domain R nX, we set

kEuk1¼ maxx2 iþy

2

jP1=4fjUij uðxi;yj;TÞjg þ maxx 2 iþy

2

jP1=4fjVij

v

ðxi;yj;TÞjg to be the error in the velocity at time T, and

kEpk1¼ maxx2 iþy

2

jP1=4fjPij pðxi;yj;TÞjg;

to be the error in the pressure. The numbers r1and r2are the approximated order of accuracy from the two consecutive

er-rors for the velocity and pressure, respectively. The number cond1is the condition number of the coefﬁcient matrix for the

unknown jump in the normal derivative of the velocity proposed in this paper, while cond2is the condition number of the

coefﬁcient matrix that we set the normal and tangential force strengths as the unknowns. One can easily see that our pro-posed method has much better conditioned system than those that set the normal and tangential force strengths as the unknowns.

InTable 1(a) and (b), the velocity are prescribed along all four sides of the rectangular domain and @Xwith the ﬁnal time T ¼ 5. InTable 1(a), the viscosity is

l

¼ 0:5 which corresponds to a small Reynolds number while inTable 1(b), the viscosity is

l

¼ 0:005 which corresponds to a larger Reynolds number. InTable 1(c), we use the Dirichlet boundary condition from the exact solution for u and x ¼ a, for

v

on x ¼ a, y ¼ c, and y ¼ d. We use the Neumann boundary condition from the exact solu-tion for u on x ¼ c, y ¼ c, and y ¼ d, for

v

on x ¼ b. The viscosity is

l

¼ 0:005 and the ﬁnal time is T ¼ 5. The purpose of this simulation is to mimic the set-up for the next example, ﬂow past a circular cylinder.

In all cases, we clearly obtain second order accuracy for the velocity. The accuracy of the pressure is between ﬁrst and second order. Note that, the pressure correction scheme proposed in[2]does not seem to improve the accuracy of the pres-sure because the boundary @Xis arbitrary.

Example 2. Flow past a circular cylinder

We now present our numerical simulations for a classical benchmark problem of simulating the ﬂow around a circular cylinder and compare our results with those published in the literature. For small Reynolds number (Re < 47), the ﬂow structure remains symmetric with stationary recirculating vortices behind the cylinder. As the Reynolds number increases, the symmetry breaks down and the vortex starts to shed up and down alternately.

The problem set-up is similar to those in the literature[21]. The inﬁnite domain is truncated to a rectangular domain ½a; b ½c; d that contains a circular cylinder centered at ðx0;y0Þ. The far-ﬁeld boundary condition from the left is given by u ¼ U1;

v

¼ 0 for t > 0. We have tested a number of parameters and set-ups. Unless stated otherwise, the plots shown in this section for ﬂow past a circular cylinder are all calculated using the ﬁnest mesh 900 450 in the domain ½10; 10 ½5; 5

(7)

with the cylinder center ð5; 0Þ to minimize the effect from the approximated boundary conditions. We also ﬁx U1¼ 1 and the radius of the cylinder r0¼ 0:5, but vary the Reynolds number Re ¼ 2r0U1=

l

¼ 1=

l

. All of the results agree with those in the literature very well, particularly with the latest second order method using vorticity–streamfunction formulation[3]. We run our tests to a ﬁnal time T ¼ 60 or longer. In the GMRES approach (without forming the coefﬁcient matrix), the time step is chosen as ð

D

tÞk¼ min h 2; h 2 maxfjUkijj þ jV k ijkg ( ) ; h ¼ minfhx;hyg; ð3:22Þ where Uk ijand V k

ijare the computed velocity at a grid point ðxi;yjÞ at time level tk. For long time computations, it is faster to

form the coefficient matrix ifDt is fixed. We take the CFL condition as 0:5 for Re 6 100, and 1=6 for Re ¼ 150; 200. We show our result at the final time T ¼ 60 or longer. Most of computations were done within a few minutes to a couple of hours depending on the mesh size. For example, for Re ¼ 100 and the mesh size 900 450, the entire simulation took about 132 min.

InFig. 2, we show the plots of the vorticity contours. The plots show very good agreement with the results in[3]where the streamfunction–vorticity formulation is used. The contours in the wake of cylinder for Re ¼ 20 has slightly wider opening than that of Re ¼ 40. For Re ¼ 80, the vortex shedding behavior has been fully developed.

InFig. 3, we show a grid convergence analysis of the drag coefficient between 0 6 t 6 60 for the case Re ¼ 40. The dash-dotted line is the result obtained from the coarse grid 226 113; the dashed line is from a finer grid 450 225; and the solid line is the result from the finest grid 900 450. One can see that as the grid refines, the drag coefficient approaches to CD¼ 1:6622 which is also listed inTable 2.

InFig. 4, we show the plots of a few selected streamlines for Re ¼ 20 and Re ¼ 40. The ﬂow remains perfectly symmetric. The size and the length of the stationary recirculating vortices behind the cylinder are also in a very good agreement with known results in the literature.

InFig. 5, we show the vorticity contour plots for different Reynolds numbers Re ¼ 100; 150; 200. The characteristic vortex shedding has been fully developed for all the cases.

InTables 2 and 3, we list the drag and lift coefﬁcients, Strouhal numbers, and compare with some of recent results in the literature. The drag and lift coefﬁcients are computed according to the following formulas given in[3]:

Re=20,−2:5:0.2:2.5 −6 −4 −2 0 2 4 6 8 −2 −1 0 1 2 Re=40,−2:5:0.2:2.5 −6 −4 −2 0 2 4 6 8 −2 −1 0 1 2 Re=80,−2:5:0.2:2.5 −6 −4 −2 0 2 4 6 8 −2 −1 0 1 2

(8)

5 10 15 20 25 30 35 40 45 50 55 60 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 1.6622 1.6575 1.5611 time drag coef.

Fig. 3. A grid refinement analysis for the drag coefficient between 0 6 t 6 60 for the case Re ¼ 40. The solid line is the result from the finest grid; the dash– dotted line is the result from the coarsest grid.

Table 2

The comparison of drag coefﬁcients for different Reynolds numbers.

Re Present Su et al.[21] Calhoun[3] Xu et al.[22] Russell et al.[20] Le et al.[10]

20 2.2069 2.20 2.19 2.23 2.13 2.05 40 1.6622 1.63 1.62 1.48 1.60 1.56 80 1.4503 1.43 – – – – 100 1.4027 1.40 1.330 1.423 1.38 1.37 150 1.3494 1.39 – – – – 200 1.2722 – 1.172 1.42 1.29 1.34 −8 −6 −4 −2 0 2 4 6 −2 −1 0 1 2 Re=20 −8 −6 −4 −2 0 2 4 6 −2 −1 0 1 2 Re=40

(9)

CD¼

l

U21 Z2p 0

x

ðhÞ þ r@

x

ðhÞ @n sinðhÞdh; ð3:23Þ CL¼

l

U2 1 Z 2p 0

x

ðhÞ þ r@

x

ðhÞ @n cosðhÞdh; ð3:24Þ

where r is the radius of the circular obstacle, and h is the angle between the outward-directed normal to @Xand the x-axis. Here, the vorticity

x

and its normal derivative@xðhÞ

@n can be computed using the velocity as

x

¼

v

x uyand

@

x

ðhÞ

@n ¼

@

x

@x cosðhÞ þ @

x

@y sinðhÞ ¼ ð

vxx

uxyÞ cosðhÞ þ ð

vxy

uyyÞ sinðhÞ ¼ ð

v

xxþ

v

yyÞ cosðhÞ ðuxxþ uyyÞ sinðhÞ: Re=100,−2:5:0.2:2.5 −6 −4 −2 0 2 4 6 8 −2 −1 0 1 2 Re=150,−2:5:0.2:2.5 −6 −4 −2 0 2 4 6 8 −2 −1 0 1 2 Re=200,−2:5:0.2:2.5 −6 −4 −2 0 2 4 6 8 −2 −1 0 1 2

Fig. 5. Vorticity plots for Re ¼ 100; 150; 200 at time T ¼ 60.

Table 3

The comparison of maximum lift coefﬁcients (CD) and Strouhal number (St) for Reynolds numbers Re ¼ 100 and Re ¼ 200.

Re ¼ 100 Re ¼ 200 CL St CL St Present 0.3068 0.162 0.605 0.204 Su et al.[21] 0.34 0.168 – – Calhoun[3] 0.298 - 0.668 – Xu et al.[22] 0.34 0.171 0.66 0.202 Russell et al.[20] 0.300 0.169 0.67 0.195 Le et al.[10] 0.323 0.160 0.43 0.187

(10)

Note that, since the vorticity and its normal derivative is needed only on the interface, we use a second order one-sided least squares interpolation, for instance,

uxðX;YÞ ¼ Xns ðxi;yjÞ2RnX

c

ijUij; uxxðX;YÞ ¼ Xns ðxi;yjÞ2RnX

c

ijUij;

to approximate those partial derivatives at a number of equally spaced points ðX_;_Y

Þ on @X. The number of grid points (in the neighborhood of ðX_;_Y_{Þ) n}

sis taken to be 9 12 but only requires a second order accurate scheme (that is, matching up

to second partial derivatives). The linear system of equations for the interpolation coefﬁcients

c

ijand

c

ijis under-determined

but is solvable. We use the singular value decomposition (SVD) to solve for the coefﬁcients. Using this approach, the com-puted partial derivatives have almost the same order of accuracy as the interpolation function itself.

We can see good agreements of our results in the drag coefficients. The lift coefficient also agrees with those in the literature for Re ¼ 100. It is a little under-estimated for Re ¼ 200, more so as in the results presented in[10,20]. One of possible explanations may be that the projection method works well for small Reynolds numbers, but not for large ones, or/ and that the domain is not large enough. We have observed that the position of the cylinder has some effect on CLwhen Re ¼ 200. InFig. 6, we plot the time evolution of the drag coefficient for the case Re ¼ 200.

InFig. 7(a)–(c), we show the vorticity plots for the ﬂow past different dumbbell-shaped objects. In (a) and (b), we have the same geometry but different

l

,

l

¼ 0:05 for (a) and

l

¼ 0:01 for (b). In (c), we have the same

l

as in (b) but with a different geometry. It is clear that no vortex shedding is observed in this ﬂow. InFig. 7(d), we show the vorticity plot of the ﬂow around two circular cylinders located at different positions in the domain with

l

¼ 0:01. The radius of the front cylinder is r ¼ 0:5, while the second one is r ¼ 0:6. The ﬂow behind the second cylinder behaves similarly as the case where only one cylinder exists if those two cylinders are far apart enough.

We have also tested a moving boundary example for the ﬂow past an in-line oscillating cylinder as described in[21]. The prescribed velocity of the cylinder is u ¼ 0,

v

¼ 0:14 sinð2

p

fctÞ, where the cylinder oscillating frequency fcis about two times of the natural vortex shedding frequency fq. The average drag coefﬁcient is CD¼ 1:693; and the lift coefﬁcient is CL¼ 1:010, which are very close to the results given in[21].

Example 3. A lid-driven cavity ﬂow with multiple obstacles

We show our numerical experiments for a lid-driven cavity ﬂow with multiple obstacles and different geometries. Similar numerical results can be found in the literature[21]. Since the focus of this paper is on the new numerical method, we thus present numerical results for a few problems without too much elaborations about their physical meanings.

In Fig. 8, we show several plots of numerical simulations. The ﬂow is driven by the horizontal velocity along the boundary y ¼ d where u ¼ 1 and

v

¼ 0. No-slip boundary conditions are prescribed along other sides of the boundary. The mesh size is 256 256. InFig. 8(a), we show the steady velocity field of the flow with five circular cylinders. We obtain qualitatively similar flow behavior as shown in[21]for

l

¼ 0:01.Fig. 8(b) and (c) shows the steady velocity field of the flow with five rose-shaped obstacles at different locations.Fig. 8(d) shows the vorticity plot of the flow corresponding to the one inFig. 8(c).

150 155 160 165 170 175 180 185 190 195 200 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 time drag coef.

(11)

μ=0.05 −10 −8 −6 −4 −2 0 2 4 6 8 −4 −3 −2 −1 0 1 2 3 4 −10 −8 −6 −4 −2 0 2 4 6 8 −4 −3 −2 −1 0 1 2 3 4 −10 −8 −6 −4 −2 0 2 4 6 8 −4 −3 −2 −1 0 1 2 3 4 −10 −8 −6 −4 −2 0 2 4 6 8 −4 −3 −2 −1 0 1 2 3 4

a

b

c

d

Fig. 7. The vorticity plots of ﬂow past different objects at ﬁnal time T ¼ 60. In (a) and (b), we have the same geometry but differentl,l¼ 0:05 for (a) and l¼ 0:01 for (b). In (c), we have the samelas in (b) but with different geometry. In (d), we have two separated cylinders withl¼ 0:01.

(12)

4. Conclusions and acknowledgment

In this paper, we propose a new augmented method for solving Navier–Stokes equations in irregular domains. This new approach significantly improves the condition number of the linear system of equations for the augmented variable com-pared with some other methods. The efficiency of the new method has been illustrated by running several interesting exam-ples. The method is very efficient for fixed obstacles and long time simulations, since we can form the coefficient matrix in the beginning and perform the matrix decomposition outside of the time loop. For a moving boundary, or a short time sim-ulation with a fixed boundary, the GMRES iteration would be more efficient.

The author would like to acknowledge the following supports. The ﬁrst and third authors are partially supported by US-ARO Grant 49308-MA, and US-AFSOR Grant FA9550-06-1-0241. The third author is also partially supported by US-NSF Grant DMS0412654, and USA NSF-NIH Grant #0201094. The second author is partially supported by National Science Council of Taiwan under Grant NSC-95-2115-M-009-010-MY2 and MoE-ATU project.

References

[1] J.B. Bell, P. Colella, H.M. Glaz, A second-order projection method for the incompressible Navier–Stokes equations, J. Comput. Phys. 85 (1989) 257–283. [2] D.L. Brown, R. Cortez, M.L. Minion, Accurate projection methods for the incompressible Navier–Stokes equations, J. Comput. Phys. 168 (2001) 464. [3] D. Calhoun, A Cartesian grid method for solving the streamfunction–vorticity equation in irregular regions, J. Comput. Phys. 176 (2002) 231–275. [4] E. Weinan, J. Liu, Projection method: I. Convergence and numerical boundary-layers, SIAM J. Numer. Anal. 32 (1995) 1017–1057.

[5] E. Guendelman, A. Selle, F. Losasso, R. Fedkiw, Coupling water and smoke to thin deformable and rigid shells, SIGGRAPH, ACM TOG 24 (2005) 973–981. [6] J. Hunter, Z. Li, H. Zhao, Autophobic spreading of drops, J. Comput. Phys. 183 (2002) 335–366.

[7] J. Van Kan, A second-order accurate pressure-correction scheme for viscous incompressible ﬂow, SIAM J. Sci. Comput. 7 (1986) 870–891. [8] J. Kim, P. Moin, Application of a fractional-step method to incompressible Navier–Stokes equations, J. Comput. Phys. 59 (1985) 308–323.

−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2

a

c

_d

b

Fig. 8. A lid-driven cavity flow with multiple obstacles. The viscosity isl¼ 0:01. (a) Velocity field for the flow with five circular cylinders, see[21]. (b) and (c) velocity field for the flow with five rose-shaped obstacles at different positions. (d) The corresponding vorticity plot of (c).

(13)

[9] M.-C. Lai, C.S. Peskin, An immersed boundary method with formal second-order accuracy and reduced numerical viscosity, J. Comput. Phys. 160 (2000) 705–719.

[10] D.V. Le, B.C. Khoo, J. Peraire, An immersed interface method for viscous incompressible ﬂows involving rigid and ﬂexible boundaries, J. Comput. Phys. 220 (2006) 109–138.

[11] Z. Li, A fast iterative algorithm for elliptic interface problems, SIAM J. Numer. Anal. 35 (1998) 230–254.

[12] Z. Li, K. Ito, The immersed interface method – numerical solutions of PDEs involving interfaces and irregular domains, SIAM Front. Ser. Appl. Math. (2006). FR33.

[13] Z. Li, K. Ito, M-C. Lai, An augmented approach for Stokes equations with a discontinuous viscosity and singular forces, Comput. Fluids 36 (2007) 622– 635.

[14] Z. Li, M-C. Lai, The immersed interface method for the Navier–Stokes equations with singular forces, J. Comput. Phys. 171 (2001) 822–842. [15] Z. Li, H. Zhao, H. Gao, A numerical study of electro-migration voiding by evolving level set functions on a ﬁxed cartesian grid, J. Comput. Phys. 152

(1999) 281–304.

[16] M.N. Linnick, H.F. Fasel, A high-order immersed interface method for simulating unsteady incompressible ﬂows on irregular domains, J. Comput. Phys. 204 (2005) 157–192.

[17] F. Losasso, G. Irving, E. Guendelman, R. Fedkiw, Melting and burning solids into liquids and gases, IEEE TVCG 12 (2006) 343–352.

[18] S. Marella, S. Krishnan, H. Liu, H.S. Udaykumar, Sharp interface Cartesian grid method I: an easily implemented technique for 3D moving boundary computations, J. Comput. Phys. 210 (2005) 1–31.

[19] R. Mittal, G. Iaccarino, Immersed boundary methods, Ann. Rev. Fluid Mech. 37 (2005) 239–261.

[20] D. Russell, Z.J. Wang, A Cartesian grid method for modeling multiple moving irregular objects in 2D incompressible viscous ﬂow, J. Comput. Phys. 191 (2003) 177–205.

[21] S.-W. Su, M.-C. Lai, C.-A. Lin, An immersed boundary technique for simulating complex flows with rigid boundary, Comput. Fluids 36 (2007) 313–324. [22] S. Xu, Z.J. Wang, An immersed interface method for simulating the interaction of a fluid with moving boundaries, J. Comput. Phys. 216 (2006) 454–493. [23] T. Ye, R. Mittal, H.S. Udaykumar, W. Shyy, An accurate Cartesian grid method for viscous incompressible flows with complex immersed boundaries, J.