Annual Reviews

(2.7)

1 72

GLOWINSKI & PIRONNEAU

n qJ +

n+ -qJ

( 1 - 28)k

n+8 +

+A ( 1 qJ

2

qJ

= ,

(2.8)

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

(2.9) Theoretical analysis and numerical experiments show that if e is properly chosen then the above scheme has much better properties than the Peace­ man-Rachford scheme concerning asymptotic behavior and the time inte­ gration of stiff systems of differential equations. In the context of the incompressible Navier-Stokes equations, the value e = 1 - 1/j2 seems to be near optimal. Another operator splitting scheme that has been widely used (par­ ticularly by the Russian applied mathematicians) is the one defined as follows:

Then for n

�

(2 . 10) 0,

n n qJ being known we compute cp +

"'n qJ 1+ -qJn n +A(" 1 qJ +l)=O, k

(2.11) (2.1 2)

Scheme (2. 1 0)-(2. 1 2) is first-order accurate and quite stable; its main drawback is that it is not well suited to capture steady state solutions, unless k is sufficiently small. A popular (and more accurate) variant of the above scheme is the one obtained by symmetrization of (2.1 0)-(2.1 2); it is defined by (2.1 0) and

-n n qJ + -qJ +A (-n+ 1 qJ k/2

"n -n qJ +1 - qJ ++A(-n+ 2 qJ k

n qJ + -qJ +A (mn+ 1 ) 1 y k/2

(2.13)

,

=

,

(2. 1 4)

0.

(2.15)

Its application to the numerical simulation of incompressible viscous flow is discussed in Beale & Majda (1981).

173

FINITE ELEMENTS FOR NAVIER-STOKES

2.2

Application to the Navier-Stokes Equations

We discuss now the application of the time discretization schemes described in the above sections to the solution of the time-dependent Navier-Stokes equations [( 1 . 1 ) and ( 1 .2)], with the initial-value condition ( 1 . 3); we suppose that the boundary conditions are of the mixed type given by ( 1 . 7) in Section 1 . If r I fO we need to have S r go' n dr o. We shall consider application of the {)-scheme only, since it produces the best results with regard to accuracy and convergence to steady states. We obtain therefore the following scheme =

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

=

Then for

n

2':

0 and starting from

u"

(2. 1 6) we solve (2. 17 a)

v

. U,,+8

= 0 in n,

(2. 1 7b) (2. 1 7c)

U,,+I-O_ U" +O -f3vAu"+ O U + 1 -8. V)u" + 1-8 fn+ 1 -8 I

_____

(l-2())k

un+

I-a

=

g�+ 1-8 on r 0,

+(

-

=

"

f3v aun+an1 -0

= g�-t- I-a

GUn+O " an

+nipn+o -IXV-- on r

un+ un+ g�-tv

I = 0 in n,

.

I

=

f3vAun+ 1 -0 _ (un+ 1-8. V)Un+ 1 -8 au"+I oun+ I -O 1-f3v on on f3 f3

in n,

+

I on r 0,

IXV--

_nipn+ 1= gn+ I

(2. 18b)

(2. 1 9a) (2. 1 9b)

on r I· (2. 1 9c)

A sensible choice for IX and is to take IX = ( l - 2())/(l-()) , = ()/(1 - ()); with such a choice many computer subprograms are common to both the linear and nonlinear subproblems, thereby saving a substantial amount of

174

GLOWINSKI & PIRONNEAU

core memory. Concerning (), numerical experiments show that () 11/)2seems to produce the best results,even in those situations where the Reynolds number is large. We observe that using the above 8-scheme we have been able to decouple the nonlinearity and the incompressibility in the Navier-Stokes equations [(1.1) and (1.2)]. We note that un+Oand un+ Iare obtained from the solution of linear problems very close to the steady Stokes problem. To conclude this section, we would like to mention that there is prac­ tically no loss in accuracy and stability by replacing (un+ 1-0. V)u"+ 1-0by (un+0 V)un+ 1-0in (2.1Sa). Remark 2.1: For () = 1/4, the convergence of scheme (2.16)-(2.19) is proved in Fernandez-Cara & Beltran (1989).

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

=

•

3.

NUMERICAL TREATMENT OF THE

ADVECTION 3.1

Introduction

The foJlowing advection equation (3.1) is a model problem for the Navier-Stokes equations because it is a sig­ nificant as a separate problem in fluid mechanics in convection dominated prob­ lems such as pollutant transport by fluids. Even though the advection equation is linear in

FINITE ELEMENTS FOR NAVIER-STOKES

1 75

technique,the least square Galerkin method,and the characteristic-Galer­ kin method. 3.2

A Third-Order Upwinding Scheme

Stability analysis of finite difference methods shows that upwinding schemes for the discretization of u . V¢ are much superior to centered schemes. The same idea can be used with the finite element method: u'V¢(x) is approximated by [¢(x)-¢(y)]/A where y x-Au, A> O. A third-order method of this type has been studied in Tabata & Fujima ( 1 99 1 )and can be described as follows.

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

=

1 . Choose 4 points x-2, X-I, Xl, x2 on the line [XO,xo + u(XO)]; x-l, X-I upwind of XO and x I, Xl downwind of XO with x-I and xI being the nearestupwind and downwind points to xO,respectively. 2. Denote by � 2, , � 2their coordinates on the line on a scale such that has coordinate �i i, i -2, -1,..., 2;call h = Ix-I-xol and

Xi

_

=

'i 3.

• • .

TI (�i-()+IXJ [ H{ i,O) = ((i-() D =I=

j-

i

=

, i ofO,

(3.3)

Define

(3.4 ) It can be shown that

(3.5) Thus the approximation is fourth-order accurate if IX = 0 and third-order accurate otherwise. If Xiiare symmetrical with respect to xO,the approxi­ mation is centered; if IX > 0 then the approximation is upwinded. When the Xiareuniformly distributed around XOand (X = 6, the finite difference scheme in Kawamura et al ( 1 985)is recovered. To approximate the advection equation (3. 1 ) one replaces u'V¢(x) by [u' V¢(x)h and uses a high order explicit Runge-Kutta scheme (RK) for the time derivatives. For example with RK2equation (3.1) is approximated by

1 76

GLOWINSKI & PIRONNEAU

(3.6) where 4>"+ 1/2is given by 1

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

k/2

[4>n+ 1/2(X) - 4>"(x)] +[u"· V4>"(x)h

=

j".

On a triangulation o f n good choices for the auxiliary points Xi are the intersections nearest to XU of the line [x,x+u(x)] with the edges of the triangles. The parameter (1can be chosen anywhere from 0to 6. 3.3

The Characteristic-Galerkin Method

The characteristic-Galerkin methods (Benque et al 1 982, Douglas & Russell 1 982, Pironneau 1 982) are derived from the analytic solution of the advection equation (Chorin 1 9 73)

04> . at +U' V4> = fill Q = n 4>(x,O) = 4>°(x),
=

x

(0, T),

(3 . 7)

VXEn,

(3.8)

g(x,t), VXE r - (t) = {x: u(x, t)· n(x) < 0, x E an}, Vt E (0, T), (3.9)

where u, J, cjJ0, g, n, r = anare given; n is the outer unit normal vector to r. In simple words h t e of ocjJ/ot. Thus (3 . 7) could be discretized in time by

1

k {4>"+ I(x) - 4>"[x- un(x)kJ} = j",

Vxdl

(3. 1 0)

To justify (3. 1 0)let {X(r)}O
; X(r) = u[X(r), r] if X(r)

E n(r)

and

0otherwise

(3. 1 1 )

with the fn i al

X(t)

=

x.

(3.12)

If u is the velocity of a fu l id then located at position xat time t. In (3.II)and (3. 1 2) xand tare parameters so X is written X(x, t; r); it is also the characteristic of the hyperbolic equation (3.7) . Noting that

FINITE ELEMENTS FOR NAVIER-STOKES

d 8¢ (x t)+u' V¢(x, t), -d ¢[X(X , t; r), r]1<=1 = ---;(jt , ,

1 77 (3. 1 3)

Equation (3.7)is rewritten as

d 8¢ at +u· V¢ = d, ¢[X(x, t; ,), ,]1.=1 Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

3.3. 1

APPROXIMATION IN TIME

(n+l )k]

=

=

f(x, t).

(3 . 1 4)

By using the fact that X[x, (n+ 1)k;

x, we can say that (3. 1 5)

where xn(x)is an approximation of X[x, (n+ l)k; nk]. Let us denote by X� an O(e)approximation of Xn(x) t[ he between the subscript of X and the precision order is due to the xn is an approximation of Xon a time intervalof size k; an O(ka) scheme gives a precision of O(ka+ I)]. For instance, choose pand {kq}�-' with � kq k; compute

1

=

xq-un(xq)kq, q = 0, . . . ,p-l. (3. 1 6) We will choose {kq} so that Xi(X)EO, i.e. Xi(O) c O. Then we derive a Xi(x) = xPwith XO

=

x; xq+

=

scheme for (3.7),namely:

(3. 1 7) where 4>oXdenotes the function x -+ 4> [X(x)] . The following result is easy to show

Proposition 3.1:

Ifu is regular and ijX'j(Q) c Qfor all n, scheme (3.17) is L2-stable. It is 0 (k) accurate when the distance between Qand X'; (Q) is 0 (k2) . Proof in the caseu . n

=

O.

We multiply (3. 1 7) by 4>n+ 1 and integrate over O. Denote by 14> I the L2 norm of ¢; then

(3. 1 8) Now by thechange of variable

y =

X,\(x)and starting from

178

GLOWINSKI & P IRONNEAU

(3 .19)

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

it is found that:

Therefore
(3.2 1 ) To obtain the error estimate w e proceed exactly i n the same way but on the difference en =
=

� [
�� (tn+ ) '

=

O(k).

So by (3 .21)

(3.23)

lenl:s; c[leol+IO(k)l]'

3. 3.2

APPROXIMATION IN SPACE

(3.22)

Let us discretize (3.15) with the fn i ite

element method;we obtain then

1


1


l

r+'whdx,

VWhEHOh,


(3.24)

where HOhis a space of continuous polynomial functions on a triangulation :Y" = {TIL of n, zero valued on the boundary r and where ghis a poly­ nomial approximation of ginside Q. In :Y", the TIare triangles such that

U, T, =Q.

Notice that (3.24)is a linear system of the type A n+ ,

=

(3.25)

Agh+b,

where

(3.26) and where {W } is a basis of H�t' and
=

i

L <1>7+ 'W . We have then the

FINITE ELEMENTS FOR NAVIER-STOKES

1 79

Proposition 3.2

Scheme (3.24) is L2(Q) stable.

The proof is done as in Proposition 3.1. We also have

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

Proposition 3.3 If HOh is the space of continuous piecewise affine functions on a triangulation of Q then the L2 (Q) error between CP'h the solution of (3.24) and the solution (pn of (3.17) is O(h2Ik+h). Therefore the scheme is O(h2Ik+h+k). Proof Let us subtract (3. 1 7) from (3.24) to obtain an L2 error. Denote +' - IIh'f' Aln +' , - Al enh+' 'f'hn

(3.27)

where IIh
1

e7,+lwhdx-

l

ehoX�Whdx =

1

HOh

of 4>"+ 1.

We

obtain

(
1

From (3.28), with Wh = e�+ we find

1

le;:+ 12

S

[ ( 1 +cluI1.ook)leJ:I +lcpn+ l_IT"cpn+ 11 + (1 +clu 11.ook)I
so that (3.29) 3.3.3 IMPL EMENTATION PROBLEMS Two points remain to be clarified, namely the computation ofX�(x), and the computation of r

=

1

CPhOX'hWhdx.

(3. 30)

Computation of (3.30):

The following Gauss quadrature formula is used (3.31) for instance with piecewise linear elements we can take {�k} = all the midpoints of the edges and wk = (h/3 in two dimensions ((h/4 in three dimensions) where (h is the area (volume) of the elements that contain ¢,k.

1 80

GLOWINSKI & PIRONNEAU

Computation of X" (x):

To compute ¢�[xn(�k)l one must address the following problem:

Given ;:k find I such that X"(�k)E T[ a triangle of the triangulation .9h of on.

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

This is not a simple problem. A possible solution is as follows: Compute all intersections of g\ X"(�k)} with the edges of the triangu­ lation, starting from �k in the direction u"(�k) and proceeding from triangle to neighboring triangle until xn(�k) is reached (see Figure I). The xi' in (3. 1 6) are these intersections; thus the numerical scheme for X"«(k) is an explicit Euler method applied within an element and when another element is reached the value of 0" on the new element is used. To compute the intersections, the barycentric coordinates {A;} of �k may be used. First on each element compute /1, such that

U(qb)

=

L

i�l, .. ,d+l

/1iqi,

L

i�l, .. ,d+l

/1i

=

0,

(3.32)

where {qi} is in the set of the vertices of the element, qh is its barycenter, and d is the dimension of the space. Then find m and p such that (3.33) To find m we assume that m is the index of the A;" that is zero, then compute p

Figure 1

Computation of XI: by following the streamline in the triangulation, On each

triangle for a given entry point (here 1=3),

(j (on

this figure j

=

2)

the exit point

(I must be computed

FINITE ELEMENTS FOR NAVIER-STOKES

p

_ =

Am . 11m

181 (3. 34)

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

If A;;;:: 0, Vi, m is the correct index. Thus most of the work goes into computing the l1i for all the elements. Because of round off errors, it may be difficult at times to decide which is the next triangle that the charac­ teristic will cross, for instance if xn(�k) is on a vertex; so careful pro­ gramming is needed. 3.4

A Least-Square Galerkin Method

Consider again the advection equation (3. 1 7). Choose a small parameter r and a smooth function w(x,t), multiply (3. 1 7) by w+r (ow/ot+u'Vw), and integrate over Q n x (0, T) to obtain: =

i [��

J[ (��

+u'V¢-I

w+r

)J

+U'VW

dXdt=O.

(3.3 5)

3 .4. 1 DISCRETIZATION Choose a time step k and construct a quad­ rangulation Th of n; denote tn = nk. In this way one obtains a proper quadrangulation of the space time slab (3.36) Finally define the space of piecewise linear functions in each component of the vector (x, t) W;:+' = {Wh E CO(Qn+ I): Wh piecewise affine in Xi' affine in ton Q;:+ '},

(3.37) (3.38)

¢

Notice that this construction gives an approximation for that is con­ tinuous in x but has jump discontinuities in time at the interfaces (tn, tn+ I) , therefore ¢h has two values at tn, which we have denoted by ¢J:(x, (n) and ¢'/, f ' (x, (n). The space-time least-square Galerkin method defines the approximate solution ¢i'+ , as a function in W'/,+' equal to gh on r and such that for all w" in Waf I we have:

1.+1 [O¢;t+

,

+u'V¢'/,+

'-IJ [ ( � Wh+r

a h +

+u· VWh

)J

r [¢;:+ '- ¢h] W" Jan

=

O.

(3.39)

In the last integral the domain of integration on is to be understood as

1 82 GLOWINSKI & PIRONNEAU + Qn ln {t tn}. The scalar T is a positive parameter of O(h+k). Notice the last term, which is a sort of penalty term to make sure that cp7,+1 � cp'/,. Then we have =

Proposition 3.4 Problem (3.39) has a unique solution.

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

Proof: Equation (3.39) is a linear system with respect to the values of cp7,+ 1 at n the vertices of the quadrangulation of Q + equations as unknowns so we only need to check that the kernel of the system is zero, i.e. that f = 9 = cpo = 0 implics cp/,+ 1 = 0 for all n. With these zero data and Wh = cp /, + I, (3.39) reduces to

Let us i ntegrate by parts the first two terms

Now the last i ntegral is positive because (3.4 1 ) yields

� r (CPh+ 2

Jon+

I

2

cp'/,+1

r (cp/,+1 )2+ r ( Jon [ cp;: Jon

(3 .42)

Define (3.43) then, using Schwartz' i nequality, (3.42) can also be written as a;;+ 1 + b;;+1 ::; 2anbn+ I, which implies i n turn that a;;+1 a;; + (bn+ 1 an)2 ::; O. So an+ that an= O. Concerning convergence the following result can be found in Johnson (1989): -

-

183

FINITE ELEMENTS FOR NAVIER-STOKES

Proposition 3.5 When L2(Q).

1"

is O(h+k) and 9 = gh the method is O[(k+ h)3/2} accurate in

NUMERICAL IMPLEM ENTATION At every time step a new non­ symmetric linear system needs to be solved; its size is twice the number of vertices for first order elements. It can either be solved by Gauss' elim­ ination or an iterative method like the Generalized Minimum Residual algorithm (GMRES) described in Saad & Schultz (1986). It may appear unnecessarily expensive to use a space-time formulation and indeed the method also works very well without it, but we have chosen to present it this way for two reasons. First the mathematical results are available in this formulation and second this formulation generaliz es without modi­ fication to moving domains and free boundaries. The elements are slanted instead of straight (see Tezduyar et al 1990).

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

3.4.2

4.

NUMERICAL TREATMENT OF THE STOKES

STEPS 4.1

Generalities

The operator splitting methods discussed in Section 2 (and also other methods) imply the solution at each time step of one or several systems of partial differential equations of the following type:

c(U-v�U+Vp

=

f inO,

(4.1)

V'U = 0 inn.

(4.2)

Using the notation of previous sections we shall take (4. 3) as boundary conditions. Here ex and v are two positive constants (ex 11k) and f, go, and g I are given functions. It can b e shown that if f, go, and g I are sufficiently smooth, then problem (4.1)-(4.3) has a unique solution if r I of ff (if r I = ff, i . e. r = r 0, u is still unique, but p is defined to within an arbitrary constant). Due to the incompressibility condition V' u = 0, problem (4.1)-(4.3) is a nontrivial one. However, suppose thatp is known in L 2(n), then we obtain u via the solution of an elliptic system whose variational formulation is given by '"

184

GLOWINSKI & PIRONNEAU

UEVg; vE o I;f

v we have a

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

wher e

Vo

=

{VIV E [H I (QW,

v=

0

lu'vdx+ v 1 VU'Vvdx= lr,vdx on ro},

Vg

=

{VIVE[H1(Q)]d, V = go on r o}.

(4.5)

Problem (4.4) can be easily solved by finite element or finite d ifference methods. This observation is the basis of powerful iterative methods for solving the Stokes problem (4. 1)-(4.3). One of these methods will be discussed in the following paragraph. Indeed these methods are sophis­ ticated variants of the following very simple algorithm:

pO E L 2(Q) is given;

(4.6)

then/or m 0, assuming that pm is known, we compute um andpm+ 1 via U"'EVg; I;fvEVo we havea lum'Vdx+vlVU""Vvdx= If'Vdx r pmV'vdx+ r g l 'vd r, In Jr, �

+

(4.7) (4.8)

Problem (4.7) is clearly equivalent to

OU'" +npm

va;:; = gl

on r I ' (4. 9)

Concerning now the convergence of algorithm (4.6)-(4.8), it follows from, e.g. Appendix 3 of Glowinski (1984) that if 0< then lim

p<

v 2d,

{um,pm} = {u,p} {u,p}

m--+ + <Xl

(4. 1 0)

in [H\Q)]d

x

L2(Q),

(4.1 1 )

where is a solution of (4. 1)--(4.3). Algorithm (4.6)-(4.8) may be slow in practice, particularly for flow at large Reynolds numbers where a 1 1k is taken very large (to follow the �

FINITE ELEMENTS FOR NAVIER-STOKES

18 5

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

fast dynamics of such flow) and where v is small. In the following Section 4.2 we shall describe a preconditioned conjugate gradient variant of algo­ rithm (4.6)-(4.8) whose convergence properties are quite good and uniform with respect to the values of v/rx. In practice, algorithm (4.6)-(4.8), and the one to be discussed in Section 4.2, are applied to finite element (or finite d ifference, or spectral) approxi­ mations of the Stokes problem (4. 1)-(4.3); we shall address in Section 4.3 the finite element approximation of problem (4. 1 )-(4.3). A Preconditioned Conjugate Gradient Algorithm for Solving the Stokes Problem (4.1)-(4.3)

For the mathematical justification of the preconditioned conjugate gradi­ ent variant of algorithm (4.6)-(4.8) described here, see, e.g., Cahouet & Chabard (1988), Bristeau et a1 ( 1 987), Glowinski & LeTallec (1989), and Glowinski (1990). The key idea (which is due to J. Cahouet) is to linearly combine two residuals, one active for IX/V» 1 (large Reynolds numbers), while the other is well suited to IX/V « 1 . Suppose that Sr, dr > 0 ; the preconditioned conjugate gradient variant of algorithm (4.6)-(4.8) is then described as follows:

pO E L 2(n) is given;

(4. 1 2)

solve o /Xu -v�u0 =f-VPoinn;

auo u0 =goonro,va;;=g,+np 0 onr" (4. 1 3)

and set (4.14)

Next, solve 0 = rO in n,' - �m 'f'

acpo 0 = Oon r , 0 on r 0, m 'f' , an =

(4.15)

and set (4.16) (4.17)

Then for m ::::: 0, assuming that pm, um, yin, gm, " " + m p " um+ r'n+ gm+', wm+ , as follows:

IV"'

are known, compute

186

GLOWINSKI & PIRONNEAU

Solve (4.18) and set pm

=

V

-OI11.

(4.19)

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

Compute (4.20) and then (4.21) (4.22) (4.23) Next, solve a

-m

:n = 0 on r

0,


(4.24)

and compute l11+ I g

= gl11 Pm(vfrn + a
(4.25)

(4.26) and then (4.27) Do m

=

m+1 andgo back to (4.18).

Algorithm (4.12)-(4.27) has proven to be quite efficient for solving Navier-Stokes problems at quite large Reynolds numbers. To conclude this section, we would like to make the following remarks:

Remark 4.1: In the case where ru = r, we should replace (4.13), (4.15), (4.18), (4.24) by

ocuO-vAuO

=

f-

FINITE ELEMENTS FOR NAVIER-STOKES

Vpo in Q;

UO

=

go

on r,

187 (4.13)' (4.15)' (4.18),

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

- Air = fm in Q;

a -m

In =

° on r,

rn ipm dx = 0, I

(4.24),

respectively. Remark 4.2: Each iteration of algorithm (4.12)-(4.27) requires the solution of one elliptic system for the operator !XI -vA. As already mentioned, for flow at large Reynolds number where oc 11k i s large and v is small the �

discrete analogues of operator oc/-vA are fairly well-conditioned matrices, making the iterative solution of these elliptic systems quite inexpensive. We also have to solve the Poisson problems (4.15) and (4.24) [or (4.15)' and (4.24)']; we shall discuss this aspect of the numerical implementation in Section 4.4. 4.3

Finite Element Approximation of the Stokes Problem

(4.1 )-(4.3)

4.3.1 GENERALITIES We have discussed in Sections 2 and 3 the time discretization of the Navier-Stokes equations (1.1) and (1.2) coupled to convenient initial and boundary conditions. The finite element treatment of the advection has also been treated in Section 3. Therefore, we still have to address the finite element approximation of the Stokes problem [(4.1) and (4.2)]. We assume that the boundary conditions are still those given by (4.3). The literature concerning the finite element approximation of the Navicr-Stokes equations is quite large (indeed, every issue of the International Journal of Numerical Methods in Fluids contains at least one article on these topics) ; concentrating on books, we shall mention Glowinski (1984), Thomasset (1981), Peyret & Taylor (1983), Temam (1977), Cuvelier et al (1986), Girault & Raviart (1986), Pironneau (1989), and Gunzburger (1989). In this article we concentrate on triangular elements since they are easier to implement and better suited to complicated geometries. It is a well accepted fact that the main difficulty related to the space approximation of the Navier-Stokes equations, in the pressure-velocity formulation, is the treatment of the incompressibility condition V' u = O. In Section 5.2 of Glowinski (1990) it was shown, using Fourier Analysis, that the mechanism producing spurious oscillations for the discrete velocity

188

GLOWINSKI & PIRONNEAU

and pressure is the fact that some discrete variant of the operator A, defined by

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

(4.28) too strongly damps those discrete pressure modes whose wave number is at least half the maximal wave n umber of the discrete velocity modes. An obvious cure to this trouble is to define the discrete pressure on a (finite difference or finite element) grid twice coarser than the one used to dis­ cretize the velocity. Indeed a finite element method for thc Stokes problem introduced in B ercovier & Pironneau (1 979) was (implicitly) based on this idea. In the above reference one proves the convergence of a finite element approximation of the Stokes problem (4.1)-(4.3) (with 0: = 0 and r 0 = r), where one uses a continuous piecewise linear pressure on a triangulation :Yh/2' twice finer than :T,,; this finite element approximation will be described in the following Section 4.3.2.

4.3.2 FUNDAMENTAL DISCRETE SPACES We suppose that n is a bounded polygonal domain ofR z. With!Y;, as a standard finite element triangulation ofn [see, e.g. Ciarlet (1978) and Raviart & Thomas (1983) for this no tion] , and h as the maximal length of the edges of :T", we introduce the following discrete spaces (with Pk = space of the polynomials in two variables of degree �k): Ph

=

{qhlqhECO(n), qhITEPI, VTE:T,,},

Vh = {vhlvhECO(n)

x

CO(n), V"ITEPz

(4.29)

x

Pz, VTE!Y;,}.

(4.30)

If the boundary conditions imply u = go on r 0, we shall need the space V h defined by O (4. 3 1 a) VOh = {vhlvhE Vh, Vh = 0 on q, ifTo = r, and by VOh = {vhlvhE Vh, Vh =

0

on

ro}, if

r dr > o. Jr,

(4 . 3 1 b)

If we are in the situation associated with ( 4.31 b) it is of fundamental importance to have the points at the interface of r 0 and r 1(= nr 0) as vertices of :T". Three useful variants of V" (and Vo,,) are obtained as follows: either Vh = {vhlvhECO(n)

x

CO(n), vhlTEPI

hlvhEco (n)

x

cO(n), vhlTEPh

V"

=

x

PI, VTE5!./z}, or

(4.32)

Ph, VTE5!.}, or

(4.33)

x

1 89

FINITE ELEMENTS FOR NAVIER-STOKES

Figure 2

To obtain the triangulation for

the velocity each triangle

of the tri­

angulation used for the pressure is divided

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

into 4 subtriangles by the mid-sides.

(4.34) In (4.32), !Y;,/2 is the triangulation of 0 obtained from f/;, by join ing the midpoints of the edges of TE!y;', as shown in Figure 2. In (4.34), �, is the triangulation of 0 obtained from!Y;, by joining the vertices to the centroid in TE�, as shown in Figure 3. Finally, in (4.33), PfT is the subspace of P3 defined as follows

Ph = {qlq = q] +A({JT, with q] EP],AER,

({JTE P3, ({JT = 0 on aT, ({JT(GT)

=

I},

(4. 35)

where, in (4.35), GT is the centroid of T. A function like ({JT is usually called a bubble function. The spaces Vh defined by (4.32), (4. 33), and (4.34) have been introduced in Bercovier & Pironneau ( 1 979), Arnold et al ( 1 984), and Pironneau ( 1 989), respectively. 4.3. 3 APPROXIMATION OF conditions are defined by

THE BOUNDARY CONDITIONS

Figure 3

If the boundary

To obtain the

triangulation for

the velocity each triangle

of the tri­ angulation used for the pressure is divided into

3 subtriangles by the center.

190

GLOWINSKI & PIRONNEAU

u =

g on r, with

r 0

g' dr

=

0,

(4.36)

it is of fundamental importance to approximate g by gh such that

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

f,

g"'Odr = o.

(4.37)

Let us discuss the construction of such a gil [we follow here Appendix 3 of Glowinski (1984)]. For simplicity, we suppose that g is continuous over r. We first define thc space 1'V" as (4.38)

i.e., yVh is the space of the traces on r of those functions Vh belonging to V". Actually, if Vh is defined by (4.38), then yVh is also the space of those functions defined over r (taking their values in R2) , continuous over r, and piecewise quadratic over the edges of � contained in r. Our problem is to construct an approximation gh of g such that

ghEYV,,,

1

9h'Odr

=

o.

(4.39)

If n"g is the unique element of yV,,, obtained by piecewise quadratic interpolation of g over r, i.e., obtained from the values taken by g at those nodes of � belonging to r, we usually have Sr nhg'0dr # O. To overcome this difficulty we may proceed as follows: (i) We define an approximation 0h of 0 as the solution of the following linear variational problem in yVh:

°hEYVh;

f,

Oh'llhdr

=

L

O'llhdr, VllhEYVh.

(4.40)

Problem (4.40) is equivalent to a linear system whose matrix is sparse, symmetric positive definite, extremely well conditioned, and easy to compute [also, problem (4.40) needs to be solved only once]. (ii) Define gh by

gh = 1l:hg-

(L

1l:"g' odr

IL 0

)0".

'Ohdr

It is easy to check that (4.40) and (4.41) imply (4.39). Now, if the boundary conditions are defined by either

(4.41)

191

FINITE ELEMENTS FOR NAVIER-STOKES

or

u= g

(4.4 2)

o

(4.43)

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

a well chosen variational formulation will a utomatically satisfy the boundary condition on rI; we shall see in the next paragraph that there is no difficulty associated with the Dirichlet condition on ro. 4.3.4 APPROXIMATION OF THE STOKES PROBLEM The problem that we consider first is defined by (4. 1)-(4.3) with ro using the spaces and defined by (4.29), (4.30) [or (4.32), (4.33), or (4.34)] and (4.3 1 a), respectively, we approximate (4. 1 )-(4.3) by:

Ph, Vh VOh

Find {Uh' Ph} E rx

Vh

X

Ph such that

LUh·v"dx+v LVu"oVv"dx- Lp"voVhdx= Lfhov"dx, 'v'VhEVO", (4.44) (4.45)

In (4.44) and (4.45),

fh gh and

(4.46)

f

are convenient approximations of and

g, respectively (see Section 4.3.3 for the construction of gil).

dr

If the boundary conditions are given by (4.3) with Sr. > 0, we approximate the problem (4. 1 )--{4.3) by the following variant of (4.44)­ (4.46):

Find {u",Ph} E Vh Ph such that LUhoVhdx+v LVUhoVv"dx- Lp"VOVhdX= Lf"OVhdX X

IY.

LVou"q"dx 'iq"EP", Uh = gOh on r = 0,

(4.48)

o.

(4.49)

There is no particular difficulty solving the above discrete Stokes prob-

192

GLOWINSKI & PIRONNEAU

lem by discrete variations of the conjugate gradient algorithm discussed in Section 4.2. 4.4

Remarks Concerning the Computer Solution of the

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

Discrete Problems

To summarize, the practical implementation of the Stokes solvers described in Section 4.2 and 4.3 requires the solution of two types of elliptic problems. If, for example, the boundary conditions for the velocity are of the Dirichlet type (i.e. if r 0 = r), we shall have to solve a problem of the type au-v�u = fin il, u =

g

on r,

(4. 50)

and a problem of the type ocp . -�cp =J;, mil, =gpon r.

(4. 5 1 )

an

Paradoxically, solving (4. 50) is not very expensive for flows at high Reynolds numbers; for such flows, the viscosity v is small, and their fast dynamics requires small �t, i.e., large values of ()(. The spectral analysis done in Dean et al ( 1 989) and Glowinski ( 1 990) shows that under such circumstances standard iterative methods (like overrelaxation and con­ jugate gradient) will have a very fast convergence when applied to the solution of the linear systems approximating problem (4. 50). On the other hand, solvin g (4. 5 1 ) may be quite time consuming, particularly for 3dimensional problem s, despite the fact that (4.51) has to be solved in the discrete pressure space whose dimension is much smaller than the dimension of the discrete velocity space. Fortunately iterative methods such as multigrid and preconditioned conjugate gradient algorithms have proven to be quite efficient for solving the discrete variants of the problem (4.51). The same conclusion holds for the other type of boundary con­ ditions considered in this article. For more details and complements, see the above references. 5.

NUMERICAL RESULTS

In the above sections we have discussed three upwinding techniques and three classes of finite element approximations. As examples of the utility of finite element methods, we shall now present some numerical results for complicated geometries. 5. 1

Flow Around a Car

The automobile industry is interested in cars with low drag because they require less gas. Also, for example, the design of the rear window is

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

FINITE ELEMENTS FOR NAVIER-STOKES

193

critical against raindrops coming from the road. Even though the Reynolds number for this problem is high, since the flow is detached and complex behind the vehicle, a direct simulation of Reynolds number around 1000 already gives a lot of information on the main eddies and on the drag. The result shown in Figure 4 is a pressure contour on the plane of symmetry of the vehicle. The computation was performed with the finite element method in the spaces defined by (4.29) and (4.33) and the Charac­ teristic-Galerkin method. The mesh was generated by a Delaunay-Voronoi mesh generator. It has approximately 50,000 vertices. The computing time for 1 00 time steps is a few hours on a CRAY 2. 5.2

Flow Over a Periodic Array of Cylinders

Some computer industries try to improve the cooling of circuit boards by intruding over the board an array of parallel rods perpendicular to the cooling-fan flow. For this problem the Reynolds number is moderate and a direct simulation is possible. After Re � 200 the flow becomes three dimensional. The program has been optimized for the Alliant FX80-8 and runs in 2 1 7 minutes for 40,000 vertices; the memory required is about 20 M bytes. The results shown in Figure 5 have been obtained with the finite element approximation associated with the spaces P I and VI> defined (in Section 4.3) by (4.29) and (4.32), respectively; the resulting nonlinear problems have been solved by one step of Newton's method [see, e.g. Dennis & Schnabel (1983) for a fairly complete description of Newton's method] mixed with a step of the characteristic-Galerkin method (see Buffat 1 99 1 ).

Figure 4

Flow around a car at Re

5000 and a viscosity of 2.5

x

10-4• The pressure is

visualized at the plane of symmetry of the car. (Courtesy of F. Hecht.) =

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

'-0 ..,.

C)

�

J-x �x

�

a

fil>

:g

�

� :»

e

F",",,5 Flow 'mood , "'ri odk 'm, of "liod "".ndan", Th, ,,, " R, � 6/10. OOOfig,,"lioo W",. 001'0", ,,1/ ;, "" """ "" ,;mol,red .ru/ panel. (Courtesy ERCOFTAC P''';od'' " odl of M. Buffat.) moellog;o "" llo", 'co ,ppliod ,"0"" M,,,h " ;01" .,d 199(). Sooood" 0'''" , 'ow;, ,ho wo ;0 'h, ',w ", right

FINITE ELEMENTS FOR NAVIER-STOKES

195

THE COMPRESSIBLE N AVIER-STOKES

6.

EQUATIONS

The Navier-Stokes equations describing the motion of Newtonian com­ pressible perfect fluids are given by: op

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

at + V . (pu) = 0,

a �� a pE --

(6. 1 ) 1

---at + V ' (pu ® u) + Vp - flAu - "3 flV(V ' u) (

at

)

=

f,

(6.2)

+ V ' (upE) + V ' (pu) = f · u

with the constitutive equation u2

E = 2" + e,

(6.4)

where fl is the viscosity, K the thermal diffusivity, f the external body forces, the i dentity tensor, and y the adiabatic constant ( = 1 .4 for air). A standard set of boundary conditions are:

I

u, e given on the boundary r; p given on r-

=

{X E r : u(x) · n(x) < O}. (6 . 5)

Relations (6.1 )-(6.4) define the velocity field u, the density p, the energy E, and temperature e. The pressure p is recovered from the law of perfect gas

p

=

(y - I)pe.

(6.6)

It follows from Matsumara & Nishida ( 1 982) that problem (6. 1 )-(6.4) has a unique solution on a time interval (0, t*) if the initial and boundary data are sufficiently small. In condensed form this system of equations is written as iJW

at + V ' F(W) - V ' G (W, VW) = 0,

(6.7)

where W [p, pu 1 0 PU2, PU 3, Ef and F (resp. G) is a nonlinear function of W (resp. W and VW). Without the G-term, (6. 7) are the Euler equations, which are clearly dominated by a dvection; they are also a system of =

1 96

GLOWINSKI & PIRONNEAU

conservation laws for which a number of sophisticated schemes have been proposed. Basically upwinding should be done in the direction of each characteristic on the diagonalized system. Without the F-term, (6.7) is a nonlinear partially parabolic system for which diffusion dominates. So we retain the idea of operator splitting and write (6.7) as

oW +A I (W)+AiW)

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

at

=

0,

(6.8)

for which a numerical scheme could be as follows: one time step of an explicit scheme for the Euler equation and one time step of an implicit scheme for the Navier-Stokes equations without the advection terms:

�[wn+ 1/2_wnl+ AI (wn)

=

0,

� [wn + I wn + 1/2] + A 2(Wn+ I ) _

7.

(6.9)

=

O.

(6.10)

THE ADVECTION STEP

There are several ways to carry out the advection step (Hughes et al 1 987, Johnson 1 989, Jameson & Baker 1 986); here we present a first-order method due to Dervieux ( 1 985). For simplicity we shall present it for two­ dimensional flows. The advection step (6.9) is in fact: (7. 1 ) T o discretize i t i n space we will use a Petrov-Galerkin method; equation (7. 1 ) is multiplied by a function that is not in the same function space as Following Dervieux ( 1 985) we multiply (7.1) by the characteristic function of a cell volume O'j and integrate; we obtain then

Wn•

I wn+ 1/2 Iwn-k IV ' F(Wn) Iwn +k L,F(Wn)'n, =

=

(7.2)

where the last formula was obtained by an integration by parts. Upwinding is introduced into (7.2) by a careful choice of O'i'

1. First construct a triangulation of Q and denote the vertices {ql 2. Then construct one cell (Ji per vertex qi delimited by the medians (resp. median planes in three dimensions) of the triangles passing through all the vertices qj neighboring qi. Then mass lumping is applied to the first two integrals in (7.2), implying:

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

FINITE ELEMENTS FOR NAVIER-STOKES

197

where I O"il denotes the area of O"i; F(W") has been replaced by an upwinded approximation, [F(W")]h, which is constructed by an adaptation of van Leer ( 1 973). The cell boundary i'JO"i is made of straight segments {oO"J H 30"i: j such that qi neighbors qi}. The art of upwinding is to choose a function cI>(W, W) such that F(W) = cI>(W, W),

(7.4)

and to set the last integral in (7.3) to

r [F(W")]h ' 0 = L: cI>(Wlui, W I,, ) ' oloai Jaa; j=/= i

n

(7.5)

oaJ

To construct cI> we notice that F being homogeneous in W, there exists a 4 x 4 matrix f4(W) such that F(W) ' n = f4(W, n)W. This matrix can be diagonalized, i.e. f4 T - 1 AT . Let us call A± diag( ± max{ ± Ah O}), f4± = T - 1 A ± T and 1 f4 1 = f4+ - gj - . Then Osher ( 1 982) suggested taking =

=

(7.6) the guiding idea being to get cI>(Wi, Wi)1 � F(Wi)1 if A/ > 0 and F(Wi) otherwise, where I denotes the components in the eigen direction cor­ responding to AI' Also when the method is applied to the advection equation, the simple upwinding [F(W)]hlaui = F(W upstream lau) is recovered. 8.

THE DIFFUSION STEP

Following (6. 1 0) we have to solve

1 98

GLOWINSKI & PIRONNEAU

�[En+ I pn+ 1 _ En + 1/2pn+ 1/2] _ V . {Kven + 1 + I {VU"+ 1 (VUn+ I )T +

This is a problem of the following type (with

=

l /k):

aup - J1Llu - 3" J1V(V ' u) = ctEp -V· {KV(E- �2)+J1[Vu+ VuT- �IV' uJU} p u Uh V� 1

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

a

(8.4)

f,

=

g,

(8.5)

where is known. Since (8.4) is linear in it is (classically) approximated by: Find E [a space of continuous piecewise linear functions (for ex­ ample), satisfying the boundary conditions] such that

ct In PUh 'Vh+ J1 In VUh' VVh+ � J1 In (V, Uh)(V 'Vh) In 'Vh Vh V� E =

f

,

(8.6)

for all E that are zero on the boundary. Once U is known, (8.5) is also linear in and can be solved in a similar way, namely:

ct lPEhwh+K 1 VEh'VWh = K 1 V�2 'VWh Wh Vh

for all E that are zero on the boundary. Excellent results were obtained at high Mach number by several authors with this method. At low Mach number the above method is not so robust because it treats the incompressibility condition in the convection step instead of the dissipation step. New ways are being explored (as in Bristeau et al 1 990b) where the full nonlinear system is discretized in time by an implicit method and the nonlinear set of equations is solved by a quasi-Newton algorithm like GMRES (Saad & Schultz 1 986). While a mixed approximation was necessary for incompressible flows, it is not clear whether it is necessary for compressible fluids but several indicators point in that direction:

FINITE ELEMENTS FOR NAVIER-STOKES

1 99

1 . Even in a compressible flow at high Mach number there are regions of low velocity where the fluid is essentially incompressible.

2. Checkerboard oscillations have been identified even at high Mach num­

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

bers (Bristeau et al 1 990b).

These oscillations can be damped by increasing the upwinding parameter or adding more artificial viscosity in the scheme but can also be eliminated by using a mixed approximation such as quadratic or linear plus bubble functions for the velocity (as in Section 4.3.2) and linear for the density and temperature. Theoretical results (Pironneau & Rappaz 1989, Fortin & Soulaimani 1 988) show that such mixed approximations are necessary on a simplified variant of (6. 1 )- (6.3), which can be obtained by assuming that the temperature is constant and that the advection velocity is small (high Reynolds number); this simplified system is given by V · (pu) = 0,

u lr

=

g,

(8.8)

p lr- = h .

(8 .9)

This problem could be named the compressible Stokes problem; it is a good test problem for the diffusion step of the full Navier-Stokes equations. 9.

ADAPTIVE MESH REFINEMENT

Obviously the performances of the computational methods are strongly influenced by the quality of the mesh. To a certain extent it is possible to monitor local mesh size when generated automatically by advancing fronts or Delaunay-Voronoi diagrams [see for example Peraire et al (1 988), Jameson & Baker ( 1 986), Georges et al ( 1 990)] . But most of the time it is only after the computation is finished that the mesh is found to be inappropriate. It is much better to adapt the mesh during the computation. The criteria for mesh refinement should be a local error estimate. For example, for the incompressible Navier-Stokes equations it is known that the difference between the computed solution Uh and the exact solution U in the case of the mini-element (see Section 4.3.2) is related to the mesh size h by (Bernardi & Raugel 1 987) (9. 1 ) Thus one way would be to divide the triangles when I u . V uI is larger than a threshold s (Erikson & Johnson 1 988).

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

200

GLOWINSKI & PIRONNEAU

Many other criteria have been proposed. In Babuska-Rheinboldt ( 1 978) and Verfuth ( 1 99 1 ) the mesh is refined when the norm of the strong residual error II Uh · VUh +VPh-VdUh II is large. In Babuska & Rheinboldt ( 1 978) and Verfuth ( 1 99 1 ) the local error is estimated by solving a local Stokes problem with the residual error on the right hand side and/or with a higher order method. In all cases the elements are subdivided. Thus the successive meshes are imbedded in each other and multigrid methods can be used to quickly solve the linear systems (Bank et al 1 987). But it is also possible to com­ pletely remesh the domain and redistribute the vertices according to the local errors as in Hecht & Vallet ( 1 99 1 ); then a complex interpolation procedure is needed. Multigrid methods can be used [at least theoretically (Scott 1 990)] but the main advantage lies in the possibility of adapting the mesh to nonisotropic singularities such as shocks and boundary layers. Time dependent singularities can even be followed by derefining/refining the mesh at every time step (Lohner 1 987). 1 0.

NUMERICAL RESULTS

The compressible Navier-Stokes equations are solved by the method described in Bristeau et al ( 1 990b), which is an adaptation to compressible flows of the methods described here for incompressible fluids. In dis­ cretizing (8. 1)-(8.3) the duality between p and V · u, is kept in mind and a different approximation is used for p, () and for u. Here p, () are chosen continuous piecewise linear on triangles and u is continuous piecewise linear on the triangulation obtained by dividing each triangle into 4 sub­ triangles by joining the mid-sides. No upwinding is used. The mesh of p, e is shown on Figure 6. The mesh is adapted to the result along the computational process with a Delaunay-Voronoi nonisotropic mesh gen­ erator for which the density of vertices can be prescribed in two directions (Hecht et al 1 99 1 ). The criterion for refinement is the second derivative matrix VVu; at each point the eigenvalues and eigenvectors are computed and the density of vertices is takcn proportional to the eigenvalue in the direction of each eigenvector. The mesh shown has 29, 1 52 triangles and 1 4,740 vertices. The method is applied to the computation of a separating flow over an airfoil in order to reproduce its stall at high angles of incidence (Figure 7). The profile is a NACA001 2 airfoil at 1 0 degrees of incidence; the Mach number is 0.2 and the Reynolds number is 1 0 5 . Although the flow is laminar, a separation bubble appears and the triangulation is correctly refined by the mesh generator at the separation bubble.

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

FINITE ELEMENTS FOR NAVIER-STOKES

Figure 6

20 1

Nonisotropic mesh refinement . (Courtesy of M. O. Bristeau and M. G. Vallet.)

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

202

GLOWINSKI & PIRONNEAU

.'

Figure 7 !1. =

..

-:'�--�­ .'

Computation of a laminar compressible separating flow at Re 1 0. (Courtesy of M. O. Bristeau and M. G. Vallet.)

=

I O j, Mach

=

0.2,

Literature Cited

Arnold, D., Brezzi, F., Fortin, M . 1 984. A stable finite element for the Stokes equa­ tions. Calcolo 2 1 : 337-44 Babuska, I., Rheinbolt, W. C. 1 978. A post­ eriori estimates for the finite element method. Int. 1. Numer. Methods Eng. 12: 1 597� 1 6 1 5 Baker, T . J. 1989. Automatic mesh gen­ eration for complex three dimensional regions using a constrained Delaunay tri­ angulation. Eng. Comput. 5: 1 61�75 Bank, R., Dupont, T., Yserentant, H. 1 987. The hierarchical basis multigrid meth­ od. Konrad-Zuse-Zentrum preprint SC v87-2

Beale, 1. T., Majda, A. 1 98 1 . Rates of convergence for viscous splitting of the Navier-Stokes equations. Math. Comput. 37: 243-60 Benque, J. P., Ibler, 8., Keramsi, A., Labadie, G. 1982. A new finite element method for Navier-Stokes equations coupled with a temperature equation. In Proc. 4th Int. Symp. on Finite elements in flow problems, ed. T. Kawai, pp. 295-30 1 .

Amsterdam: North-Holland Bercovier, M . , Pironneau, O. 1 979. Error estimates for finite element method solu­ tion of the Stokes problem in primitive variables. Numer. Math. 33: 2 1 1-24 Bernardi, c., Raugel, G. 1985. A conforming finite element method for the time-depen­ dent Navier-Stokes equations. SIAM J. Numer. Anal. 22: 455-73

Bristeau, M. 0., Bernardi, C., Pironneau, 0., Vallet, M. G. 1 990a. Numerical analysis for compressible isothermal flows. Proc. Symp. on Applied Math. Venice Bristean, M. 0., Glowinski, R., Mantel, 8., Periaux, J., Perrier, P. 1 985. Numerical methods for incompressible and com­ pressible Navier-Stokes problems. In Finite Elements in Fluids, ed. R. Gallagher, G. Carey, J. T. Oden, O. C. Zienkiewiez, Vol. 6, pp. 1-40. Chichester: Wiley Bristeau, M. 0., Glowinski, R., Dutto, L., Periaux, J., Roge, G. 1990b. Compressible flow calculations using compatible finite element approximations. Int. J. Numer. Methods Fluids 1 1 (6): 7 19-49 Bristeau, M . 0., Glowinski, R., Periaux, J. 1987. Numerical methods for the Navier­ Stokes equations. Compul. Phys. Rep. 6: 73- 1 87 Brown, P., Saad, Y. 1 990. Hybrid Krylov methods for nonlinear systems of equa­ tions. SIAM J. Sci. Stat. Comput. I I : 45081 Buffat, M . N . 1 99 1 . Numerical simulation of complex three dimensional flows. Int. J. Numer. Methods Fluids 1 2: 683-704 Cahouet, J., Chabard, J. P. 1 988. Some fast 3D solvers for the generalized Stokes problem. Int. J. Numer. Methods Fluids 8: 269-95 Chorin, A. 1. 1 973. Numerical study of slightly viscous flows. J. Fluid Mech. 57: 785-96

FINITE ELEMENTS FOR NAVIER-STOKES

1978. The Finite Element Method/or Elliptic Problems. Amsterdam;

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

Ciarlet, P. G.

North-Holland Cuvelier, c., Segal, A., Van Steenhaven, A . 1986. Finite Element Method and Navier­ Stokes Equations. Dordrecht; Reidel Dean, E. J., Glowinski, R., Li, C. H. 1 989. Supercomputer solution of partial differ­ ential equation problems in com­ putational fluid dynamics and control. Computer Physics Communications 53; 401-39 Dennis, J . E., Schnabel, R. B . 1 983 . Numeri­ cal Methods for Unconstrained Opti­ mization and Nonlinear Equations. Engle­ wood Cliffs, NJ; Prentice-Hall . Dervieux, A. 1985. Steady Euler Simulations Us ing Unstructured Meshes. Von Karman Lecture notes series 1 8 84-04 Dougla s, J. , Russell, T. F. 1 982. Numeri cal methods for advection dominated diffusion problems based on combining the method of characteristics with finite element methods or finite difference method. SIAM J. Numer. Anal. 19(5); 87 1-85 Erikson, K., Johnson, C. 1988. An adaptive finite element method for linear elliptic problems. Math. Comput. 50; 361-83 Fernandez-Cara, E . , Bel tran , M . M . 1 989. The convergence of two numerical schemes for the Navier-Stokes equations. Numerishe Mathematik 55; 3 3-60 Fortin, M . , Soulaimani, A. 1 988. Finite element approximation of compressible viscous flows. In Proc. Comput. Methods in Flow Analysis, ed. H. Niki, M. Kawa­ hara. Okayama; Univ. Sci. Press Gabay, D. 1 98 3. Application of the methods of multipl iers to variational inequalities. In Augmented Lagrangian Methods, ed. M. Fortin, R. Glowinski, pp. 299-33 1 . Amsterdam; North-Holland Georges, P . L., Hecht, F., Sallel, E. 1 990. Automatic 3D mesh generation with pre­ scribed mesh boundaries. IEEE Trans. Magn. 26(2); 7 7 1 -74 Girault, V., Raviart, P. A. 1 986. Finite

203

V . Balakrishnan, A. A. Doronitsyn, J. L. Lions, pp. 57-95. New York: Opti­ mization Software Glowinski, R. 1 989. Supercomputing and the finite element approximation of the Navier-Stokes equations for incom­ pressible viscous fluids. In Recent

A dvances in Computational Fluid Dyna­

mics, ed. C. C. Chao, S. A. Orszag, W. Shyy, Lecture Notes in Engineering, Vol. 43. Berlin; Springer-Verlag Glowinski, R . 1 990. Finite element methods for the numerical simulation of incom­ pressible viscous flow. Introduction to the control of the Navier-Stokes equations. Math. Dept. Research Report UH/MD-94. Univ. Houston, Houston, Texas Glowinski , R . , LeTallec, P. 1989. Augmented

Lagrangian and Operator Splitting Methods in Nonlinear Mechanics. Phila­

delphia: SIAM lecture notes Godlewsky, E. 1980. Methodes a pas mul­ tiples et de direction alternees pour la dis­ cretization d'equations d'evolution. These de 3e Cycle. Univ. P. et M . Curie, Paris Gunzburger, M. 1989. Finite Element Methods/or Viscous Incompressible Flows. Boston: Acad. Press Hecht, F., Vallet, M . G. 1 99 1 . Non isotro pic automatic mesh generation. In p repara­ tion Hughes, T. J. R., Franca, L. P., Mallet, M . 1 9 8 7 . A new finite element formulation for CFD. VI Convergence analysis of the generalized SUPG formulation for linear time dependent multi-dimensional advec­ tive diffusive systems. Comput. Methods Appl, Mech. Eng. 63: 97-1 1 2 Jameson, A . , Baker, J. 1 986. Improvement to aircraft Euler method. AIAA 25th Aero­

space Meet., Reno

York: Springer-Verlag Glowinski, R. 1 985. Vi scous flow simulation by finite element methods and related numerical techniques. In Progress and

Johnson, C. 1 989. Numerical Methods for Partiul Differential Equations. Cambridge; Univ. Press Johnson, C., Szepessy, A . , Hansbo, P. 1 990. On the convergence of shock-capturing stream diffusion FEM for hyperbolic con­ servation laws. Math. Comput. 54(1 89): 1 07-30 Kawamura, T . , Takami, H., Kuwahara, K . 1 985. New high order upwind scheme for incompressible Navier-Stokes equations. In Numerical Methods in Fluid Dynamics, Lecture Notes in Physics, Vol. 2 1 8, pp. 29 1-95. Berlin: Springer Kreiss, H. 0., Lorenz, J. 1 989. Initial-Bound­

banel, pp. 1 73-21 O. Boston: Birkhauser Glowinski, R. 1986. Splitting methods for the numerical solution of the incom­ pressible Navier-Stokes equations. In Vistas in Applied Mathematics, ed. A.

Ladysenskaya, O. A. 1 969. The Mathe­ matical Theory 0/ Viscous Incompressible Flows. New York: Gordon & Breach LeVeq ue, R . , Oli ger , J. 1 983 . Numerical Methods based on additive splitting for

Element Methods/or Navier-Stokes Equa­ tions. Berlin: Spri nger-Verlag Glowinski, R. 1 984. Numerical Methods /or Nonlinear Variational Problems. New

Supercomputing in Computational Fluid Dynamics, ed. E. M . Murman, S. S. Abar­

ary Value Problems and the Navier-Stokes Equations. Boston: Acad. Press

204

GLOWINSKI & PIRONNEAU

hyperbolic partial differential equations. Math. Comput. 40: 469�97 Lions, J. L. 1969. Quelques Methodes de

Annu. Rev. Fluid Mech. 1992.24:167-204. Downloaded from arjournals.annualreviews.org by Universidade Federal de Minas Gerais on 12/04/09. For personal use only.

Resolution des Problemes aux Limites Non Lineaires. Paris: Dunod Lions, P. L., Mercier, B. 1 979 . Splitting algo­

rithms for the sum of two nonlinear oper­ ators. SIAM 1. Numer. Anal. 16: 964�79 Lohner, R. 1 987. 3D grid generation by the advancing front method. In Laminar and Turbulent Flows, ed. C. Taylor, W. G. Habashi. Swansea, UK: Pinneridge Marchuk, G. I. 1 975. Methods of Numerical Mathematics. New York: Springer-Verlag Matsumura, A., Nishida, T. 1980. The initial value problem for the equations of motion of viscous and heat conductive gases. 1. Math. Kyoto Vniv. 20: 67�104 Osher, S., Chakravarthy, S. 1982. Upwind difference schemes for hyperbolic systems of conservation laws. Math. Comput. 38: 1 32�54 Peraire, J., Pciro, J. , Formaggia, L., Morgan, K., Zienkiewicz, O. C. 1 988. Finite element Euler computations in three dimensions. Int. 1. Numer. Methods Eng. 26: 2 1 3 5�59 Peaceman, D., Rachford, H. 1 955. The numerical solution of parabolic and ellip­ tic differential equations. 1. SIAM 3: 28� 41 Peyret, R . , Taylor, T. D . 1 983. Com­ putational Methods for Fluid Now. New York: Springer-Verlag Pironneau, O. 1 982. On the transport­ diffusion algorithm and its applications to the Navier-Stokcs equations. Numer. Math. 38: 309�32 Pironneau, O. 1989. Finite Element Methods (or Fluids. Chichester: Wiley Pironneau, 0., Rappaz, 1. 1989. Numerical analysis for compressible isentropic vis­ cous flows. IMPA CT I : 1 09�37

Raviart, P. A., Thomas, J. M. 1 983. Intro­ duction a rAnalyse Numerique des Equa­ tions aux Derivees Partielles. Paris: Mas­

son Saad, Y., Schultz, M. H. 1 986. GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM 1. Sci. Stat. Comput. 7: 856-69 Scott, R., Zhang, S. 1990. Higher dimen­ sional nonnested multigrid methods. Research Report VH/MD-84. Univ. Hou­ ston, Houston, Texas Strang, G. 1 968. On the construction and comparison of difference schemes. SIAM J. Numer. A nal. 5: 506- 1 7 Tabata, M . , Fuj ima, S. 1 99 1 . An upwind finite element scheme for high Reynolds number flows. Int. 1. Numer. Methods Fluids 1 2(4): 305�22 Tartar, L. 1 978. Topics in Nonlinear Analysis. Univ. Paris-Sud, Orsay Temam, R. 1 977. Navier-Stokes Equations. Amsterdam: North-Holland Tezduyar, T., Liou, J., Behr , M. 1 990. A new strategy for FEM computations involving moving boundaries and interfaces l. Thomasset, F. 1 98 1 . Implementation of Finite VMSI report 90/169

Element Methods for Navier-Stokes Equa­ tions. New York: Springer-Verlag van Leer, B. 1 973. Toward the ultimate con­ servative difference scheme. In LeclUre Notes in Physics, Vol. 1 8, pp. 1 63�68. New

York: Springer

Verfurth, R. 1 99 1 . A posteriori error esti­

mators and adaptive mesh-refinement techniques for the Navier-Stokes equa­ tions. In Incompressible CFD: Trends and A dvances, ed. M . D. Gunzburger, R. A. Nicolaides. Cambridge: Univ. Press Yanenko, N. N. 1 97 1 . The Method of Frac­ tional Steps. Berlin: Springer-Verlag