- Email: [email protected]
Stabilized finite element methods: II. The incompressible Navier-Stokes equations
Computer Methods in Applied Mechanics and Engineering North-Holland
99 (1992) 209-233
CMA 258
Stabilized finite element methods: II. The incompressible Navier-Stokes equations* Leopold0 Laboradrio
P. Franca
National de ComputaGiio CientljLica (LNCCICNPq), 22290 Rio de Janeiro, RJ, Brazil
Rua Lauro Miller 455,
SCrgio L. Frey Departamento
de Computaqa”o, lnstituto de Matemhtica, Universidade Federal Fluminense PraGa do Valonguinho slno., 24210 Niterdi, RJ, Brazil
CUFF),
Received 21 August 1991
Stabilized methods are proposed and analyzed for a linearized form of the incompressible Navier-Stokes equations. The methods are extended and tested for the nonlinear model. The methods combine the good features of stabilized methods already proposed for the Stokes flow and for advective-diffusive flows. These methods also generalize previous works restricted to low-order interpolations, thus allowing any combination of velocity and continuous pressure interpolations. A careful design of the stability parameters is suggested which considerably simplifies these generalizations.
1. Introduction Simulations of the incompressible Navier-Stokes equations with the classical Galerkin method may suffer from spurious oscillations arising from two main sources. The first is due to the advective-diffusive character of the equations, and as in the scalar advective-diffusive equation, oscillations may contaminate the velocity field for high advection, or, in this context, for high Reynolds number. The second source has to do with the mixed formulation character of the equations, which limits the choice of combinations of the finite element interpolations used to approximate the velocity and pressure fields. Employment of combinations that do not satisfy a compatibility condition (or inf-sup, or BabuSka-Brezzi [l, 21 condition) may yield undesirable pathologies in the approximation of pressure and velocity. To deal with these potential shortcomings, finite element researchers have proposed to simulate the incompressible Navier-Stokes problem with a pair of interpolations that fulfills the inf-sup condition, and modify the advective operator to include some ‘upwinding’ effect. While this seems a reasonable procedure, if applied without further thought, instabilities may still prevail, or excessive numerical diffusivity may affect the final accuracy of the results. Correspondence
to: Leopold0 P. Franca, LNCC, Rua Lauro Miiller 4.55, 22290 Rio de Janeiro, * For Part I, see [16].
0045-7825/92/$05.00
0
1992 Elsevier Science Publishers
B.V. All rights reserved
RJ, Brazil.
210
L. P. Franca, S. L. Frey, Stabilized jinite element methods: II
For some of the Galerkin methods and related methodologies employed in this equation, the interested reader is referred to [3-S] and references therein. In this paper we continue our search for improved stabilized methods, now applied to the incompressible Navier-Stokes equations. In Section 2, we first address a linearized model for these equations, propose stabilized methods, make error convergence studies and show an interesting relationship with the Galerkin method using the mini-element. The methods studied deal with the inf-sup condition requirements in the same manner as in the methods introduced in [9-121 for the Stokes equations. Namely, by modifying the variational equations in a careful way, we circumvent the need to satisfy the BabuSka-Brezzi condition. Also, the methods stabilize the advective operator in the streamline direction, as in the methods for the scalar advective-diffusive equation (SUPG [13, 141, Galerkin-least-squares [15] or as in [16], referred to as Part I henceforth). While one version of these methods has already been studied recently [17, 181, these authors limited their presentation to low-order interpolation combinations, such as Pl/Pl or Ql/Ql for the velocity-pressure pair. To generalize these methods for high order interpolations, the design of the stability parameter is readdressed as in Part I, using a recent computation of inverse estimate constants [19]. Although our analysis is performed for a linearized model, in Section 3 we extend the methods to the incompressible Navier-Stokes equations and all the features anticipated in the analysis presented in Section 2 are confirmed numerically. A predictor multi-corrector scheme is employed to integrate the methods to the steady-state solutions and concomitantly to deal with the nonlinearity of the equations, closely following all the suggestions proposed in [18, 201. The preliminaries are the same as in Part I; that is, the linearized model and the incompressible Navier-Stokes equations are defined on a bounded domain fl C RN, N = 2, 3, with a polygonal or polyhedral boundary r. A partition %‘hof fi into elements consisting of triangles (tetrahedrons in R”) or convex quadrilaterals (hexahedrons) is performed in the usual way (i.e., no overlapping is allowed between any two elements of the partition; the union of all element domains K reproduces fin, etc.) and a combination of triangles and quadrilaterals for the two-dimensional case can be accommodated. Quasiuniformity is nut assumed. In what follows C, C,, C,, C,, . . . denote various positive constants independent of viscosity and of any element diameter h,, KE GCh. As usual, L’(R) is the space of square-integrable functions in R ; Lo, the space of L2-functions with zero mean value in 0; (- , * ) denotes the L2-inner product in 0 and 1)* /lo, the L*(0)-norm. We also employ to denote the L2-inner product and norm in the element domain K, (* 7*)K andII *ll0.K respectively. Also, %‘“(a) is the space of continuous functions in a, HA(a) is the Sobolev space of functions with square-integrable value and derivatives in fl with zero value on the boundary r, and the H’-norm is denoted by )I * 1)1. For convenience, we adopt the following notation: if K is a triangle or tetrahedron, if K is a quadrilateral or hexahedron where for each integer m 20,
P,,, and Q, have the usual meaning.
,
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
2. A linearized incompressible Navier-Stokes
211
model
2.1. Stabilized methods A steady-state,
linearized,
incompressible
Navier-Stokes
model is given by
(Vu)a - ~YV. E(U) + Vp =f
in L? ,
(14
v-u=0
in J2 ,
(lb)
Zf=O
on r,
(14
where a is a given velocity field, u is the unknown velocity, p is the pressure, v is the viscosity, E(U) is the symmetric part of the velocity gradient andf is the body force. For the discussion below, the homogeneous Dirichlet boundary condition (lc) suffices (for an interesting account on various boundary conditions for the incompressible Navier-Stokes, see [213). The finite clement spaces we wish to work with are given by Vh =
{v E H;(o)”
( VI, E R,(lqN,
K E Vh}
(2)
)
where the integers k and 1 are greater than or equal to 1. In the methods that follow, any combination of k and 1 is allowed, a possibility ruled out, in general, in the standard Galerkin method. The methods we wish to study can be written as: find u, E Vh and ph E Ph such that
with B,(u,
p;
v, q) = ((Vu)a,
v) +
(Zm(u),
E(U)) - (V* v, p) - (V. u, 4) + (V* u, av* v)
+ Kzg ((Vu)a + VP - 2s
44
~.(x, Re,(~))((W~
h
-vq
2 2247. E(v)))K
(5)
and
U-% d = (A u) + K;v(f, 4x, Re,(x))((Vv)a -Vq +-m7- +WK , h
where the stability parameters
S and r are defined as
(6)
212
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
(11)
mk = min{ 4, 2C,}
(12)
,
‘k K,q,, ‘:llV*E(U)IIi,~ d II+‘)II; ,
u E V, ,
(13)
and A > 0 is a positive parameter. For comparison purposes, let us also review the SUPG method for the model [13, 17, 18, 22, 231. Find u,, E V,, and ph E Ph such that B S”Pa&h> phi u, q) = FS”PG(U, q) 7 (‘3 q) E
‘h
’
‘/,
(14)
3
with B S”P&’
p; u, 4) = (Wb,
4 + (244
E(U)) - (V* UTP) - (V* u7 4)
+ (V. u, 6V. u) + KFz,, ((9~
+ VP - 2vV.
44,4x7
Re&))((Wa
-
WK (15)
and (.L 7(x, Re.(x))((Wa
-W)K
.
(16)
h
REMARK
1. In the Stokes limit, i.e., when a = 0 above, the ‘plus’ method reduces to the Galerkin-least-squares method of [ll], the ‘minus’ method reduces to the method introduced in [lo], and the SUPG method reduces to the Brezzi and Pitkaranta method [24] up to a correction on the right-hand-side (for generalization of the latter in the Stokes context, see also [9]). For these methods, there is no need to satisfy the standard BabuSka-Brezzi condition, and equal-order interpolation can be used for velocity and pressure, a result that is preserved for a # 0, as shown in the next section.
REMARK 2. If for a moment, we drop the pressure term from the above formulations (i.e. set p = q = 0 throughout), then for the advective-diffusive equations that remain, the ‘plus’
method reduces to the method analyzed in Part I, the ‘minus’ method reduces to the Galerkin-least-squares method of [15], and the SUPG method is just the original one proposed in [13], up to the new definition of the stability parameter introduced in Part I. As thoroughly discussed in Part I, the stability considerations for each method are similar in connection with the advective-diffusive model. However, with respect to the Stokes operator
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
213
part, the method inspired in the Douglas and Wang method (the ‘minus’ formulation above) is more insensitive to the choice of the stability parameter. Different values of mk (eq. (12)) still produce desirable numerical solutions for this formulation, as anticipated in the stability analysis, presented in the next subsection.
REMARK
3. The ‘plus’ and ‘minus’ formulations
are direct extensions of the stabilized methods that we have been investigating for some years now, but apparently not introduced before for the incompressible Navier-Stokes equations. They both reduce to the SUPG formulation for the linearly interpolated velocity approximation. In this case, with an equally linear pressure, Hansbo and Szepessy [17] have made an analysis and some computations of the SUPG method, with a design of the stability parameter T, which is always 0(/z,) in their case. By assuming that the Reynolds number is always sufficiently high, they do not distinguish the locally diffusive flow cases in their choice of T. Tezduyar et al. [18] have also studied a similar formulation, making a proper distinction between locally advective-dominated and locally diffusive-dominated flows. One of the designs of T in [lS] reduces to our formulae (8)-(11) by setting mk = f and by dropping inequality (13). While this is satisfactory for linear interpolations, it breaks down for high order interpolations. A related formulation is also proposed by Lube and Auge [25]. In the next section we show that the ‘plus’ formulation is unstable for larger values of 7 than those proposed above when we employ high-order interpolations for the velocity field. In particular, for equal-order serendipity interpolation, using a ‘wrong’ value of m = 4 instead of mkz2 = 0.0814814. . . as suggested by recent computation of inverse estimates [19], the numerical simulations present spurious oscillations for both velocity and pressure in the ‘plus’ formulation (see Section 3.3).
2.2. Erroranalysis In this section we study the convergence of the ‘plus’ and ‘minus’ methods given by (4)-(13). Let us first note that the stability parameter T is bounded by a constant in each element domain K. By definition 7(x, Re&))
= 2,a;:),
2
Re,(x) 2 1 ,
(17a)
P 7(x3 Re,&N
Therefore,
for Re.(x) 7(X, Re,(x))
e.4 = -jg-
,
OaRe,(x)<
1.
(17b)
3 1, =
hK 2l+)l,
1 Re,(x)
mMx)l,h, 4~
(18) and combining with the definition (17b), we conclude that the bound (18) is valid for all values of Re,(x).
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
214
We are now ready to establish the stability characteristics LEMMA 2.1. (Stability). Assuming (I) V*a(x)=O, (II) v(x) = v = constant > 0 , then for all (u, p) E V,, X Ph , B,(u,
PROOF.
p; u, -p)
the given data a(x) and v(x) to satisfy
2 $(24IE(4Il(:
uE
ll~1’2(P4~ + VP>lli>+ w2v
+
First note that by integration
((Vu>a,u) = 0,
of each method.
by parts, by assumption
4l:, ’
(I) and by (2),
v,,,
since u=O on rfor uEV,. Next, for each formulation we have (i) The ‘plus’ formulation: By (5), B+(u, p; u, -p> = 0 + 24E(U)II:,
+ VP)ll:, - Kg
+ ll~1’2((wJ
Combining
+ llg1’2v. UlK
l17”22~v’ 44lL
(19)
.
(18) to estimate the last term of (19) gives Kzq,, 1171’22vV*e(~)Ili,~ d mky”’
2
h;llV. &(~)ll,?;,~ (by (II) and (18))
KE%;I,
(by(13)) (20)
(by (12)) * Combining this estimate with (19) yields the desired result. (ii) The ‘minus’ formulation: Let p > 1. By (5), B_(u,
p; u, -p)
= 2V((E(U)((; + ((PV* + K&
IIT
ulic: + IlT”‘((v~)a
+ Vp)ll,:
“*2Vv E(U)Il(:.K- 2 K& (2Vv’ +)7
>2V((E(U)((; + ((PV.
T((Vu)a + %))K
p-‘>~~~“‘<~ + Vp>((:,
z& + (l-
+ (l - P) Kzq,, ~~T1’22vv’ E(U)Il:,.K 2
(1
+
;(I
-
k9&4+411(:
+
(1
-
fi-‘)h”‘((v+
+ )(61’2v. ull; .
For the specific choice of p = 2, the stability estimate follows.
+
vP)ll:, (21)
Cl
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
215
Next, let us recall that by standard approximation theory [26], for a regular family of elements, there exist interpolants u”hlK E R,(QN and p”hlK E R,(K) such that JJU- U”hlKJJm.K S Ch~l-mlU~k+,,K
vu E Hk+‘(qN,
I(p - P”hlKllm,K =sChym(pI,+,.K
vp E HI+‘(K),
OG m d k + 1 ) 09 m f I+ 1 .
(The expression regular family of elements has the usual meaning The following interpolation estimate can now be established.
(22)
(23)
[26]).
LEMMA 2.2. Assume that the solution to (1) satisfies u E Hk”(fl)N fl H:,(R)N and p E H’+‘(O) rl L:(R). Denoting by T-J,= u”,,- u and 7, = ch - p the interpolation errors, for each KEY?,,, we have VxEK, then (a) if Re .al,
11~-“*%ll~.K +2410?uN~.K + 11~“*P?,)4I~.K + ll7”*2~v*o-L)ll~.K + IP”*v’ %ll~,K+ ll~-“*~,llL + ll~1’2v77pll:),K s Cc;~~l”l~hZKk+ll~I:.l.K (6)
if 0 S Re,(x)
+ SE; l4,‘h~+‘lPl:+l.K)
;
(24)
< 1, V’xE K, then
Il~-“z~Pll~.K +24Ie?.)llL + Il~“*
(2%
Therefore,
ll~-“*11,ll~ +241+7U)II~ + II~“2w?u)41~ + lP”*\J*sll,: + Kzs [lb“*m7-+~)ll~.~ + HO% - ~W”2~,Il:).K + H(;- Re,W-‘lh,ll&J
+ ll~“*vrlpll~
d C[ c hZ,kl&+‘,K(H(ReK - l)h, KEY&
+
sup /al,, + H(1
-
Re,)2v)
XEK
hZ,‘l~lf,‘,dHWK - W, ;Ec ld,’ + HU - &)h:@+‘)]
KTe
, (26)
h
where H( * ) is the Heuviside function given by
(27)
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
216
PROOF.
The velocity estimates related to the momentum equation follow as in Part I, noting that the obvious extensions apply for a vector field. For pressure and velocity (related to the continuity equation) estimates, let us divide the proof into two cases. (a) Let Re,(x) 3 1, Vx E K. Then 116li2v* %ll;,K +
ll~-1’2.t7,11~,R + ll~1’2Vrlpll;.K = II(M~)l,~K)1’2V* ?IuIlL
+ lI(~l~(~)l,~,)-1’217,11~.K + ll~21~~~~~~~1’2~~,:.~ -h/c ;~p,bl,llV* s.lli.K+ ;‘J; d
(28)
C sup lalphZKk+lkl:+,,R + sup l~~,‘h;+‘l~lf+,,~). ( XEK
(b) Let 0 d Re,(x) Ils1’2v’
Therefore,
l~l,l(~-‘h,‘~~~pll~,K + h,cIhpll;,K)
XEK
< 1, VX E K. Then dli,K
+
G
2 2vllv’
d
C@‘h2,kbd:+~.K
(2v)-1b?p~~~.K
+
%IIi.K
(24)-(26)
+
+
Ib1’2vt7pll;l.K
(2V)-1(b?pll:,K
+
h~llv~pll~,K)
(29)
(2v)-1h2,‘+2bt+1.K).
follows from (28), (29).
q
We may now establish the following convergence
estimate.
THEOREM 2.1, Under the same hypotheses as in Lemmas 2.1 and 2.2, the solutions (u,, , p,,) converge to (u, p), solutions of (1) as follows: of the methods given by (4)-(13)
4 i 2Wuh
-
u>lli+ ll~1’2(w% - u>a + VP,
s c KTz (h~bl~+I,K(H(R% - ‘lhK ;EF
-
b/p
P>>lI3 + lW2V*(Uh- 4lli +
HU
-
ReKP)
h
+
h2K PI~+I.AH(R~K -
l)h,
sup [al,’ + H(l - Re,)h:(2u)-l))
,
XEK
where H( * ) is defined in (27). PROOF. Let e: = u,, - kh, ei = p,, - @,, and ey = ei + Q, proof is as follows:
ep = ei + vp. The convergence
L.P. Frwma, S. I.. Frey, Stabilized finite element methods: II
217
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
218
Before continuing,
note that
$ 3 H(Re,
d
+
$
H(Re,
Combining with the previous yS = 1 f H(Re, - l), gives
- 1)//S 1’2V. eiII:_, + $1
- l)ila -i’2n,lI:.~ + ; estimate,
~(2vl(iz(e~)l)~ f l171'2((Ve:)a
and choosing
-I- Ve~)ll~) +
H(1 - Re.)2++j:)~l;,,
H(1 -
Re,)(2d-‘h,//;,,]
‘y, = 3 + 2NH(l - Re,),
.
‘yz= 4 and
1116 “*V*ef:lli
-117-1’2rf&+K~~h[f(3+2h'H(1-Re~))2vllF(r),)ll~.~ + (7 + 2NH(l - Re,))ll7”‘(V~,>all~,~
e
+ t(l f H(Re,
+-2(7 + 2NH(l-
Re,))(fr”*2vV.
C 4(3 +2N)H(l
- ReK)(2v)-‘jj77,(j&
- l))llS”2V.
E( ~,)I/~,, + H(Re,
rl,lli,K
- l)ll~-“~~~ll~~~
+ 2(7 + 2NH(l - Re,))~/~“*V~,~~~,~]
c(ll~-1’211,11~ + WI&(rl”)lli + ll~“2mubII,:+ IP *‘2v*s,lG &f~~)lli,K+ ff@eK- l>llS-1’2~p /Ii\, + & trlT1’22~V.
d C
2 (~~lul~+~,K(H(ReK - l)h,
sup 101~-t H(l - Re,)2v) XEK
KE%h
+ GIPI?+,.AWR~K - l)h, The last inequality
;EJ la/,’ -t H(l
-
(30)
Re,)hi(2v)-‘)) .
follows from Lemma 2.2, which also gives that
aP4147?u)ll: + Il~““~P~,~~ + b,)Il3 + 4ll~“2vef?,II~ s c KzT
(~~~~l~+~,K(H~ReK
-
l)h,
sup lalp + H(1-
Re,)2v)
XEK
h
+
and the result follows inequality. q
h:‘l!#+l,K(H(ReK
by redefining
-
l)hK
constants
,sF;
la/,’
+
ff(1
-
ReK)h;&‘)-I))
,
(31)
in (30) and (31) and using the triangle
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
2.3. Relationship
219
with the Galerkin method using the mini-element
Pierre [27, 281 and Bank and Welfert [29] have shown, for the Stokes flow, the equivalence between the Galerkin method employing the mini-element of Arnold et al. [30] with the stabilized methods of Brezzi and Pitkaranta [24] and Hughes et al. [12], for a particular choice of the stability parameter emanating from the specific definition of the bubble function used. Recently, Brezzi et al. [31] have indicated that for advective-diffusive equations a similar equivalence holds. They have shown that the bubbles not only help in circumventing the compatibility of space requirements, but also stabilize the advective operator, as established for a linearized model of the compressible Navier-Stokes equations. In this section we use similar arguments to show the relation between the Galerkin method employing the minielement and the SUPG method, (14)-(16), w h en applied with equal-order linear elements for velocity and pressure to the two-dimensional linearized incompressible Navier-Stokes equations. Let us start by recalling the finite element spaces for the mini-element:
vi={uE H#qN
) UJK E (P,(K)@
qzqy,K E
%,}
(32)
and
where B(K) denotes standard bubble functions (such as the cubic bubble in [30]) defined on K with the following properties:
I cj K
do=
C,,hi
= C,, meas(
llWll”,K= c3,3
$1, E B(K), K E Vi’,, ,
(34)
AK E W), KE Q% 3
(35)
where CIK, C,, and C,, are independent of the element parameter h,. Then the Galerkin method reads: find u,, E Vi and p,, E Pi such that
B,(“,>
ph;
u; 4) = &@, 4))
@, 4) E
v;’
p:, ,
(36)
with B,(u,
p;
K&
4) = (J-Yu) *
u, 4) = ((VU)&
4
+ (244,
+.J))
- (V* u, p> - (V* u, q>
(37)
and (38)
Following the demonstration in [31], the first step is to ‘statically condensate’ the bubble from the momentum equation. To this end, let us decompose the trial function for the velocity as ‘h@)
=
%(d + K& +(-+,
3
h
where ur(x) E V; = {u E HA(0))N ( u(, E P,(ZC)~, K E %,}, ub is the bubble nodal value and
AK E B(K) c fwo
Next let us take u = 4ei, on K, for i = 1, . . . , N, and zero elsewhere, where ei is the unity Cartesian vector with components (ei)j = 6, for j = 1, . . . , N. For this choice of weighting
L. P. Franca, S. L. Frey, Stabi~ize~~n~te element methods: II
220
functions
and using (39) in (36) with 4 = 0, we obtain (tV$)%
cbe,)~ + (2VE(#Ub), e(&&
In obtaining (40), we have used integration and f are constant in K) to obtain ((W&M,
+ (VP, rbe,)~ = (.L
#ei)x
(40)
e
by parts (and we have assumed that the data a, Y
(41)
(~VE(Z~I),
E(+ei))K
(42)
= 0 -
Note that
(43)
=Az+,-e,. Therefore
the ‘bubble equation’ 2vAu, *e, = -((Vu,)a
(40) becomes + Vp -f,
(44)
4ei)K .
Or, using the piecewise constant data assumption, z+, = - ;v A-‘((Vu&
+ Vp -f)l,
I, 4 d0 .
(45)
Having computed the bubble nodal value, let us examine its effect on the Galerkin method with u = ut E VL. Thus by using the decomposition (39) in (37), we obtain ‘C&z,
Ph;
‘1,
d
=
((vu,)a~
%>
+
Kzq
(@buda,
+C
+
(2VE(Ulh
E(Ul))
h
-
(V- Ul,
Ph)
- P*
U1r
(46)
4) -
where we used the fact that (2ve(4uu,), e(u,)), = 0. Note that integrating combining the two terms above under the sum, we obtain
B&h,
ph; u,7
4) = &t”,
2
Ph; ul,
d -
K;qtub+,tvulb
-‘dK
7
by parts and
(47)
k
which is just the Galerkin operator for equal-order linear elements plus a perturbation term. This additional term can be further worked out with the help of (34), (35) and (45) to yield
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
PT = - &
221
(u,+, (Vv,)a -V&C h
Kzsh $,A-‘((Vu,b+ VP-f)i.>((V+
=
=c
‘l;:k
h; A-‘((Vu&z + Vp -f,
-oq),(/K 4 da)*
(Vu,)u -V&
(48)
.
KE’Gh
Therefore the perturbation terms that emanate by condensation are in fact adding a term as in the SUPG method [13] with a stability parameter matrix defined as ClkC2khi 'mini
=
A-l
2u
REMARK
4. Up to the matrix character, the stability parameter is form-identical to the formula (17b) for locally diffusive dominated behavior (Re, < 1). However, for locally advective dominated flows, this formula has been shown to be ineffective and to recover the O(h,) for the T parameter. In this case, either one modifies the condensation procedure in ad-hoc fashion, or a new definition of the bubble function must be introduced. Both alternatives are discussed in [31]. REMARK 5. A relation between stabilized methods and the Galerkin method with bubble functions has also been noticed by Arnold [32] and Pitkaranta [33] for the Reissner-Mindlin model of plate bending. In particular, Arnold [32] shows the equivalence between the Galerkin method introduced by Arnold and Falk [34] and the Galerkin-least-squares method of Franca and Stenberg [35].
3. The incompressible Navier-Stokes equations 3.1.
Stabilized methods
Let us now consider the incompressible
Navier-Stokes
equations,
$+(VU)~-~VV.E((~)+V~=~
innX(O,T),
v*u=o
in 0 X (0, T) ,
u=g
on rg X (0, T) ,
an=h
on r, x (0, T) ,
u = 240
inRatt=O,
given by
(50)
where in addition to the notation introduced in Section 2.1, rg and r, are complementary subsets of r to accommodate given Dirichlet (g) and Neumann (h) boundary conditions, and
L. F. Franca,
222
S. L. Frey, Stabilized finite element methods: II
u is the stress tensor given by o=-pz+/e!v&(u).
Using the finite element
spaces given by (2) and (3), and defining
we can now extend the methods presented follows. Find uI, E Vf and ph E P, such that
Bty(u,, ph;
0,
q) =
in Section 2.1, in a semi-discrete
fashion,
F:Yu>4) > tu*4) E 62x Ph
as
(53)
’
with BY{% p; u, 4) =
2
+ (Vz+, u) f (22+),
E(fJ)) - (V* u, p) - (V. @>4)
u + Vp - 287. E(U), 7(x, Re,(x))((Vu)u -v~t2~v.~(u)))K+(V.u,6V.v)
(54)
and J?=%, q) = (f, r.9 + (k u)~~ + c
(f, 7(x, Re&))((V@
-Vq I: 2~7. ~(0)~
KEqk (55)
where 6 and the stability parameter
T are defined as
(57)
L. P. Frunca, S. L. Frey, Stabilized finite element methods: II
REMARK 6. The SUPG 0. e(u) = 0 in (54)-(55).
method
is obtained
by setting
223
in the weighting
function
slot,
REMARK 7. The methods above are form-identical to the methods introduced in Part I and studied in detail for the linearized set of Navier-Stokes equations in Section 2. Here, we do not attempt to perform a convergence analysis for the nonlinear model. This is done for SUPG employing equal-order linear elements in [17]. We should emphasize that all the main numerical difficulties in the simulation of the nonlinear model are also present in the linearized model of Section 2, namely the simulation of high advective flows and the compatibility of the velocity and pressure spaces. There it is shown that these stabilized methods deal with these features appropriately and simultaneously. We understand that this is an important result per se, and it is applicable to various high order interpolations, with different stability features for the different methods studied, which is confirmed numerically for the nonlinear model in Section 3.3. In the next subsection we show how the nonlinearity is dealt with as part of the general predictor multi-corrector algorithm employed in [18, 201 to march in time to steady solutions. 3.2. Time integration algorithm Discretization of (53) is carried out by expanding the trial functions { uh, ph} and {u, 4) in terms of their finite element basis or shape functions. This leads to the following set of semi-discrete equations:
[M +
MT]a +
N(u) + N”,(u) + [K + K", + K”,]u +
[M”,]a + N;(u) + [G' + K4,]u + [GCI,]p = E, ,
[G
+
Gy]p
=
F
+ F, , (63)
where u and p are the vectors of degrees of freedom for uh and p,, , respectively and a is the vector of degrees of freedom for au,, / at. The matrices [Ml, [K] and [G] emanate from time-dependent, viscous and pressure terms, respectively. The nonlinear vector N(u) is derived from the advective term. The other matrices arise from the additional T- and a-terms. To define the solution procedure to (63), let us introduce the following notation: n ,_, the maximum number of time steps, i max, the number of corrector passes through the algorithm, (Y, a parameter governing the approximation of a. Following [20] and [18], we consider a one-step generalized trapezoidal method given by the following predictor multi-corrector algorithm. ALGORITHM 3.1 Step 1. Initialize ~“0,p”, and ai. Step 2. Set n:=O and i:=O. Step 3. Predictor phase: u(i) *n+1*-
uw .n+l
‘-
PCi) n+l'_--
.;) + At(1 - a)~;) ,
0, p(i) n
*
(64)
L. P. Franca, S. L. Frey, Stabilized
224
Step 4. Form and factorize the incremental
finite element methods:
II
matrices [M*], [G*], [G*‘], [G,] defined by
[M*] := [M] + [M;] + a At[N*] ,
(UpI,)] +[Z (&,I +[K] + [K;] + [K;]
[N*] := [g [G*]:= [G] CC*‘] :=
+
)
[G;],
[M;] + a At( j[G'] +
[z (u;;,)]+[K;]).
(65)
Step 5. Form the residual vectors R(‘)
.-
s(i)
._
n+1*-
n+1--
Fm+,- ([M + iqay;,+N(uy:,) + Ny(u;j,) + [K + KY + Kfi’]z(;, + [G + G;]&) ,
F n+l +
E
7n+,
([M;]a;;, + N’f(z&) + [G' + K’f]u;;, + [G;]p;;,)
-
Step 6. Solve the incremental
.
(66)
system
(67) Step 7. Corrector
phase:
(68)
Step 8. If i < i,,,,
then i := i + 1 and return to Step 4; otherwise,
u(i)
._
$1
:= &wJ
n+1.-
n+l
P@) ??+I*-
u(i,,,)
n+l
n+l
. _
i := 0 and
p(i,,,) n+l
3
(6%
)
.
Step 9. If n < nmax, then II := y1+ 1 and return to Step 3; otherwise stop. REMARK suggestions
8. The matrices given by
defined in (65) are simplified
following Tezduyar
et al.‘s [18]
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
[iv”] := [M] [G"]
:=
[G]
)
(70)
,
:= a At[G’]
[G*']
225
.
With these simplifications, the matrices in (65) do no longer depend on uzi 1. This implies that the matrices can be formed and factorized once and for all, before entering the loops of time step and iterations. In this case, one may consider that Step 4 is performed before the time loop starts. Another advantage is that the system given in Step 6 becomes symmetric with these simplifications. 3.3. Numerical results In this section we present some numerical experiments Navier-Stokes equations using the methods described in serendipity approximations were employed for all tests performed on SUN Sparcstations l+ of LNCC using the MODULEF/INRIA and VISU2/Laval. 3.3.1.
with the nonlinear incompressible Sections 3.1 and 3.2. Equal-order reported. All computations were graphics facilities of the programs
The leaky cavity problem
The problem consists of the flow in a unit square with velocity boundary conditions: u, = 1, u,=Oaty=l(O
Next, we perform the computations for the refined mesh employing the ‘minus’ formulation with Re = 400. The number of nodes are practically the same as used by Tezduyar et al. [18] with a Ql/Ql combination of interpolations, and their results are qualitatively similar to those in Fig. 3. In Fig. 3, the flow is locally diffusive dominated, since by definition (see (58)) the element Reynolds number Re,, in the refined mesh for the flow with maximum unit velocity is equal to Re
=
K
Considering
wlww 4(1/400)
= 0.52 < 1 .
a flow with Re = 5000, locally advective dominated
(71) flows are now present and the
L. P. Franca. S. L. Frey, Stabilized finite element methods: II
226
Velocity
Stream
Horizontal
function
velocity
Pressure
elevation
at x = 0.50
0.
-0.
0. 705 P
=I
0. 428
-0.
0. 150
-0.
0.0
contours
Pressure
at P = 0.50
0. 982
-0, 127
Pressure
-0. 0.2
0.5
0.7
;Y Fig. 1. ‘Minus’ formulation: serendipity elements).
steady-state
1.0
Y leaky cavity flow for Re = 400 employing a coarse mesh (8 x 8 Q2S/Q2S-
227
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
“Minus”
“Plus”
Velocity
Velocity
ir
..,**,*,.I..
. * t + .
b
’
r.o~~t\~.,
* .
.
3 .
.
.
,
,
. .
-
.
,
.
,
, c
,
.
.
c
. .
.
. .
. .
,
I
.
.
.
.
.
. .
.
.
+
.
1
,,,,
+,t
,
sLeCC,.
.
.
c
-
..I
.
.
-
.,t.*\\.-._&_
4
..LC’/...
.
*
l
.
.
-
///I
.I\\-,
v,;
.,t\tt,.*,
.
. .
-.
.
.
c
,
,
.
-
.
*
.
(
+
\.
.
.
-
.
h.
*.
.
*
_
*.
.
..*
_
Pressure contours
Pressure contours
0 -
Pressure elevation
.
.
( +. . . . . . . . . . . . . . * . * 4
.
/
-F-C,,
c
*
.
/
.
1
n
Pressure elevation
Fig. 2. Comparison between ‘minus’ and ‘plus’ formulation with a ‘wrong ’ value of m (m = 4): steady-state cavity tlow for Re = 400 employing 8 X 8 Q2S/Q2S_serendipity elements.
leaky
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
228
Pressure contours
Velocity I
- f--?-Y-?-7+.--:-~?+-----’
. L L
.C
I...
,:
,v..
.:
.
.
.
LICCC#. .
.
.**CCLCa+CC*.
.
.
.
.,
.
.
.
.
.
.,
,
*....***.....<
*
. .
.
. .
.
.
.
.
.
. .
. .
.
.
.
.
.
. .
. .
.
.
. .
. .
.
. .
.
. .
.
.
.
.
.
Pressure elevation
Stream function
Pressure at 5 = 0.50
Horizontal velocity at x = 0.50
Ul
1. 000
0
0. 697
a
a. 394
-0
091
-a
025
-0
037
211
0.0
0.2
0.5
0.8
Fig. 3. ‘Minus’ formulation: steady-state Q2S/Q2S_serendipity elements).
1.
leaky
cavity
flow for
0.0
0.2
Re = 400 employing
0.5
0.8
a refined
1.0
mesh
(16 x 16
L.P. Franca, S. L. Frey, Stabilized finite element methods: iI
229
Pressure contours
St ream function
Pressure elevation
0 ~ I
d
-
Pressure at 2 = 0.50
l3orizontal velocity at T = 0.50
Ul
1.000
0.004
0.717
0.001
P
0.433
0.150 -oat33
-0.002 -0.005
a.0
0.3
0.5
0.8
Y Fig. 4. ‘Minus’ formulation: elements.
steady-state
1.0
-0.008
r
0.0
0.3
0.5
0.6
rlo
Y leaky cavity flow for Re = 5000 employing16 X 16 Q2S/Q2S+erendipity
230
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
ihl
611
Fig. 5.
The backward step flow: (a) problem statement; (b)
mesh employed
(224 Q2S/Q2S_serendipity
elements).
results employing the minus formulation are shown in Fig. 4. (We take A = 1 in the 6 formula, (56).) The results obtained are similar to those found in [36]. 3.3.2. The backward step flow The problem statement and the mesh employed are shown in Fig. 5. We take Re = 60 with respect to the inflow section. For the Q2S element, this is basically a diffusive dominated flow, by computing Re, similarly as in (71). The results for the ‘minus’ and ‘plus’ formulations are shown in Fig. 6 and 7, respectively. In Fig. 8, the ‘minus’ formulation is tested with a free-traction outflow condition. Note that the results are similar to the Dirichlet outflow test, except for a region near the outflow. In all computations, we observe the recirculation at the step, which is absent for a Stokes flow.
Fig. 6. ‘Minus’ formulation for the steady-state elevation; (c) pressure contours.
backward
step flow with Re = 60: (a) velocity
vectors;
(b) pressure
L.p.
Franca, S. L. ‘Frey, Stabilized finite element methods: 11
Fig. 7. ‘Plus’ formulation of the steady-state elevation; (c) pressure contours.
231
backward step flow with Re = 60: (a) velocity vectors; (b) pressure
Fig. 8. ‘Minus’ formulation for the steady-state backward step flow with Re = 60 using free-traction condition: (a) velocity vectors; (b) pressure elevation; (c) pressure contours.
as outflow
232
L. P. Franca,
S. L. Frey, Stabilized
finite element methods:
II
Acknowledgment The authors wish to thank Tayfun Tezduyar for fruitful discussions and useful suggestions and Michel Fortin, Michel Lapointe and Roger Pierre from Lava1 University for the graphics post-processing code VISU2. During the course of this work Leopold0 P. Franca has been partially supported by CNPq Procs. 300419/90-2 and 402598/90-3.
References [l] I. BabuSka, The finite element method with Lagrangian multipliers, Numer. Math. 20 (1973) 179-192. [2] F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from Lagrange multipliers, RAIRO Ser. Rouge 8 (1974) 129-151. [3] 0. Pironneau, Finite Element Methods for Fluids (Wiley, New York, 1989). [4] R.L. Sani, P.M. Gresho, R.L. Lee and D.F. Griffiths, The cause and cure (!) of the spurious pressures generated by certain FEM solutions of the incompressible Navier-Stokes equations: Part l., Internat. J. Numer. Methods Fluids 1 (1981) 17-43. [5] R.L. Sani, P.M. Gresho, R.L. Lee, D.F. Griffiths and M. Engelman, The cause and cure (!) of the spurious pressures generated by certain FEM solutions of the incompressible Navier-Stokes equations: Part 2, Internat. J. Numer. Methods. Fluids 1 (1981) 171-204. [6] A. Soulaimani, M. Fortin, Y. Ouellet, G. Dhatt and F. Bertrand, Simple continuous pressure elements for two- and three-dimensional incompressible flows, Comput. Methods Appl. Mech. Engrg. 62 (1987) 47-69. [7] C. Taylor and P. Hood, Numerical solution of the Navier-Stokes equations using the finite element technique, Comput. & Fluids 1 (1973) l-28. [8] R. Verfiirth, Finite element approximation of incompressible Navier-Stokes equations with slip boundary condition, Numer. Math. 50 (1987) 697-721. [9] F. Brezzi and J. Douglas, Stabilized mixed methods for Stokes problem, Numer. Math. 53 (1988) 225-236. [lo] J. Douglas and J. Wang, An absolutely stabilized finite element method for the Stokes problem, Math. Comp. 52 (1989) 495-508. [ll] T.J.R. Hughes, L.P. Franca, A new finite element formulation for computational fluid dynamics: VII. The Stokes problem with various well-posed boundary conditions: Symmetric formulations that converge for all velocity/pressure spaces, Comput. Methods Appl. Mech. Engrg. 65 (1987) 85-96. [12] T.J.R. Hughes, L.P. Franca and M. Balestra, A new finite element formulation for computational fluid dynamics: V. Circumventing the BabuSka-Brezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations, Comput. Methods Appl. Mech. Engrg. 59 (1986) 85-99. [13] A.N. Brooks and T.J.R. Hughes, Streamline upwind/Petrov-Galerkin formulations for convective dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1982) 199-259. [14] C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method (Cambridge Univ. Press, Cambridge, 1987). [15] T.J.R. Hughes, L.P. Franca and G.M. Hulbert, A new finite element formulation for computational fluid dynamics: VIII. The Galerkin-least-squares method for advective-diffusive equations, Comput. Methods Appl. Mech. Engrg. 73 (1989) 173-189. [16] L.P. Franca, S.L. Frey and T.J.R. Hughes, Stabilized finite element methods: I. Application to the advective-diffusive model, Comput. Methods Appl. Mech. Engrg. 95 (1992) 253-276. [17] P. Hansbo and A. Szepessy, A velocity-pressure streamline diffusion finite element method for the incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 84 (1990) 175-192. [18] T.E. Tezduyar, R. Shih, S. Mittal and S.E. Ray, Incompressible flow using stabilized bilinear and linear equal-order-interpolation velocity-pressure elements, Report UMSI 90/165, University of Minnesota/Supercomputer Institute, Minneapolis, 1990.
L. P. Franca, S. L. Frey, Stabilized finite element methods: II [19] I. Harari and T.J.R.
233
Hughes, What are C and h?: Inequalities for the analysis and design of finite element methods, Comput. Methods Appl. Mech. Engrg. 97 (1992) 157-192. [20] T.E. Tezduyar and T.J.R. Hughes, Development of time-accurate finite element techniques for first-order hyperbolic systems with particular emphasis on the compressible Euler equations, Report prepared under NASA-Ames University Consortium Interchange No. NCAZ-OR745104, 1982. [21] P.M. Gresho, Incompressible fluid dynamics: Some fundamental formulation issues, Annu. Rev. Fluid Mech. 23 (1991) 413-4.53. [22] T.E. Tezduyar, M. Behr and J. Liou, A new strategy for finite element computations involving moving boundaries and interfaces - The deforming-spatial-domain/space-time procedure: 1. The concept and the preliminary numerical tests, Comput. Methods Appl. Mech. Engrg. 94 (1992) 339-351. [23] T.E. Tezduyar, M. Behr, S. Mittal and J. Liou, A new strategy for finite element computations involving moving boundaries and interfaces - The deforming-spatial-domain/space-time procedure: II. Computation of free-surface flows, two-liquid flows, and flows with drifting cylinders, Comput. Methods Appl. Mech. Engrg. 94 (1992) 353-371. [24] F. Brezzi and J. Pitkaranta, On the stabilization of finite element approximations of the Stokes problem, in: W. Hackbusch, ed., Efficient Solutions of Elliptic Systems, Notes on Numerical Fluid Mechanics, Vol. 10 (Viewig, Wiesbaden, 1984) 11-19. [25] G. Lube and A. Auge, Regularized mixed finite element approximations of incompressible flow problems. II. Navier-Stokes flow, Otto von Guericke University, Preprint 15/91, 1991. [26] P.G. Ciarlet, The Finite Element Method for Elliptic Problems (North-Holland, Amsterdam, 1978). [27] R. Pierre, Simple C” approximations for the computation of incompressible flows, Comput. Methods Appl. Mech. Engrg. 68 (1988) 205-227. [28] R. Pierre, Regularization procedures of mixed finite element approximations of the Stokes problem, Numer. Methods Partial Differential Equations 5 (1989) 241-258. [29] R.E. Bank and B.D. Welfert, A comparison between the mini-element and the Petrov-Galerkin formulations .for the generalized Stokes problem, Comput. Methods Appl. Mech. Engrg. 83 (1990) 61-68. [30] D.N. Arnold, F. Brezzi and M. Fortin, A stable finite element for the Stokes equations, Calcolo 23 (1984) 337-344. [31] F. Brezzi, M.-O. Bristeau, L.P. Franca, M. Mallet and G. RogC, A relationship between stabilized finite element methods and the Galerkin method with bubble functions, Comput. Methods Appl. Mech. Engrg. 96 (1992) 117-129. [32] D.N. Arnold, Innovative finite element methods for plates, Report No. AM 61, Department of Mathematics, PennState, 1990; to appear in Mat. Apl. Comput. (331 J. Pitkaranta, Analysis of some low-order finite element schemes for Mindlin-Reissner and Kirchhoff Plates, Numer. Math. 53 (1988) 237-254. [34] D.N. Arnold and R.S. Falk, A uniformly accurate finite element method for the Reissner-Mindlin plate, SIAM J. Numer. Anal. 26 (1989) 1276-1290. [35] L.P. Franca and R. Stenberg, A modification of a low-order Reissner-Mindlin plate bending element, in: J.R. Whiteman, ed., The Mathematics of Finite Elements and Applications VII, MAFELAP 1990 (Academic Press, London, 1991) 425-436. [36] P.M. Gresho and S.T. Chan, On theory of semi-implicit projection for viscous incompressible flow and its implementation via a finite element method that also introduces a nearly consistent mass matrix. Part 2: Implementation, Internat. J. Numer. Methods. Fluids 11 (1990) 621-659.
99 (1992) 209-233
CMA 258
Stabilized finite element methods: II. The incompressible Navier-Stokes equations* Leopold0 Laboradrio
P. Franca
National de ComputaGiio CientljLica (LNCCICNPq), 22290 Rio de Janeiro, RJ, Brazil
Rua Lauro Miller 455,
SCrgio L. Frey Departamento
de Computaqa”o, lnstituto de Matemhtica, Universidade Federal Fluminense PraGa do Valonguinho slno., 24210 Niterdi, RJ, Brazil
CUFF),
Received 21 August 1991
Stabilized methods are proposed and analyzed for a linearized form of the incompressible Navier-Stokes equations. The methods are extended and tested for the nonlinear model. The methods combine the good features of stabilized methods already proposed for the Stokes flow and for advective-diffusive flows. These methods also generalize previous works restricted to low-order interpolations, thus allowing any combination of velocity and continuous pressure interpolations. A careful design of the stability parameters is suggested which considerably simplifies these generalizations.
1. Introduction Simulations of the incompressible Navier-Stokes equations with the classical Galerkin method may suffer from spurious oscillations arising from two main sources. The first is due to the advective-diffusive character of the equations, and as in the scalar advective-diffusive equation, oscillations may contaminate the velocity field for high advection, or, in this context, for high Reynolds number. The second source has to do with the mixed formulation character of the equations, which limits the choice of combinations of the finite element interpolations used to approximate the velocity and pressure fields. Employment of combinations that do not satisfy a compatibility condition (or inf-sup, or BabuSka-Brezzi [l, 21 condition) may yield undesirable pathologies in the approximation of pressure and velocity. To deal with these potential shortcomings, finite element researchers have proposed to simulate the incompressible Navier-Stokes problem with a pair of interpolations that fulfills the inf-sup condition, and modify the advective operator to include some ‘upwinding’ effect. While this seems a reasonable procedure, if applied without further thought, instabilities may still prevail, or excessive numerical diffusivity may affect the final accuracy of the results. Correspondence
to: Leopold0 P. Franca, LNCC, Rua Lauro Miiller 4.55, 22290 Rio de Janeiro, * For Part I, see [16].
0045-7825/92/$05.00
0
1992 Elsevier Science Publishers
B.V. All rights reserved
RJ, Brazil.
210
L. P. Franca, S. L. Frey, Stabilized jinite element methods: II
For some of the Galerkin methods and related methodologies employed in this equation, the interested reader is referred to [3-S] and references therein. In this paper we continue our search for improved stabilized methods, now applied to the incompressible Navier-Stokes equations. In Section 2, we first address a linearized model for these equations, propose stabilized methods, make error convergence studies and show an interesting relationship with the Galerkin method using the mini-element. The methods studied deal with the inf-sup condition requirements in the same manner as in the methods introduced in [9-121 for the Stokes equations. Namely, by modifying the variational equations in a careful way, we circumvent the need to satisfy the BabuSka-Brezzi condition. Also, the methods stabilize the advective operator in the streamline direction, as in the methods for the scalar advective-diffusive equation (SUPG [13, 141, Galerkin-least-squares [15] or as in [16], referred to as Part I henceforth). While one version of these methods has already been studied recently [17, 181, these authors limited their presentation to low-order interpolation combinations, such as Pl/Pl or Ql/Ql for the velocity-pressure pair. To generalize these methods for high order interpolations, the design of the stability parameter is readdressed as in Part I, using a recent computation of inverse estimate constants [19]. Although our analysis is performed for a linearized model, in Section 3 we extend the methods to the incompressible Navier-Stokes equations and all the features anticipated in the analysis presented in Section 2 are confirmed numerically. A predictor multi-corrector scheme is employed to integrate the methods to the steady-state solutions and concomitantly to deal with the nonlinearity of the equations, closely following all the suggestions proposed in [18, 201. The preliminaries are the same as in Part I; that is, the linearized model and the incompressible Navier-Stokes equations are defined on a bounded domain fl C RN, N = 2, 3, with a polygonal or polyhedral boundary r. A partition %‘hof fi into elements consisting of triangles (tetrahedrons in R”) or convex quadrilaterals (hexahedrons) is performed in the usual way (i.e., no overlapping is allowed between any two elements of the partition; the union of all element domains K reproduces fin, etc.) and a combination of triangles and quadrilaterals for the two-dimensional case can be accommodated. Quasiuniformity is nut assumed. In what follows C, C,, C,, C,, . . . denote various positive constants independent of viscosity and of any element diameter h,, KE GCh. As usual, L’(R) is the space of square-integrable functions in R ; Lo, the space of L2-functions with zero mean value in 0; (- , * ) denotes the L2-inner product in 0 and 1)* /lo, the L*(0)-norm. We also employ to denote the L2-inner product and norm in the element domain K, (* 7*)K andII *ll0.K respectively. Also, %‘“(a) is the space of continuous functions in a, HA(a) is the Sobolev space of functions with square-integrable value and derivatives in fl with zero value on the boundary r, and the H’-norm is denoted by )I * 1)1. For convenience, we adopt the following notation: if K is a triangle or tetrahedron, if K is a quadrilateral or hexahedron where for each integer m 20,
P,,, and Q, have the usual meaning.
,
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
2. A linearized incompressible Navier-Stokes
211
model
2.1. Stabilized methods A steady-state,
linearized,
incompressible
Navier-Stokes
model is given by
(Vu)a - ~YV. E(U) + Vp =f
in L? ,
(14
v-u=0
in J2 ,
(lb)
Zf=O
on r,
(14
where a is a given velocity field, u is the unknown velocity, p is the pressure, v is the viscosity, E(U) is the symmetric part of the velocity gradient andf is the body force. For the discussion below, the homogeneous Dirichlet boundary condition (lc) suffices (for an interesting account on various boundary conditions for the incompressible Navier-Stokes, see [213). The finite clement spaces we wish to work with are given by Vh =
{v E H;(o)”
( VI, E R,(lqN,
K E Vh}
(2)
)
where the integers k and 1 are greater than or equal to 1. In the methods that follow, any combination of k and 1 is allowed, a possibility ruled out, in general, in the standard Galerkin method. The methods we wish to study can be written as: find u, E Vh and ph E Ph such that
with B,(u,
p;
v, q) = ((Vu)a,
v) +
(Zm(u),
E(U)) - (V* v, p) - (V. u, 4) + (V* u, av* v)
+ Kzg ((Vu)a + VP - 2s
44
~.(x, Re,(~))((W~
h
-vq
2 2247. E(v)))K
(5)
and
U-% d = (A u) + K;v(f, 4x, Re,(x))((Vv)a -Vq +-m7- +WK , h
where the stability parameters
S and r are defined as
(6)
212
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
(11)
mk = min{ 4, 2C,}
(12)
,
‘k K,q,, ‘:llV*E(U)IIi,~ d II+‘)II; ,
u E V, ,
(13)
and A > 0 is a positive parameter. For comparison purposes, let us also review the SUPG method for the model [13, 17, 18, 22, 231. Find u,, E V,, and ph E Ph such that B S”Pa&h> phi u, q) = FS”PG(U, q) 7 (‘3 q) E
‘h
’
‘/,
(14)
3
with B S”P&’
p; u, 4) = (Wb,
4 + (244
E(U)) - (V* UTP) - (V* u7 4)
+ (V. u, 6V. u) + KFz,, ((9~
+ VP - 2vV.
44,4x7
Re&))((Wa
-
WK (15)
and (.L 7(x, Re.(x))((Wa
-W)K
.
(16)
h
REMARK
1. In the Stokes limit, i.e., when a = 0 above, the ‘plus’ method reduces to the Galerkin-least-squares method of [ll], the ‘minus’ method reduces to the method introduced in [lo], and the SUPG method reduces to the Brezzi and Pitkaranta method [24] up to a correction on the right-hand-side (for generalization of the latter in the Stokes context, see also [9]). For these methods, there is no need to satisfy the standard BabuSka-Brezzi condition, and equal-order interpolation can be used for velocity and pressure, a result that is preserved for a # 0, as shown in the next section.
REMARK 2. If for a moment, we drop the pressure term from the above formulations (i.e. set p = q = 0 throughout), then for the advective-diffusive equations that remain, the ‘plus’
method reduces to the method analyzed in Part I, the ‘minus’ method reduces to the Galerkin-least-squares method of [15], and the SUPG method is just the original one proposed in [13], up to the new definition of the stability parameter introduced in Part I. As thoroughly discussed in Part I, the stability considerations for each method are similar in connection with the advective-diffusive model. However, with respect to the Stokes operator
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
213
part, the method inspired in the Douglas and Wang method (the ‘minus’ formulation above) is more insensitive to the choice of the stability parameter. Different values of mk (eq. (12)) still produce desirable numerical solutions for this formulation, as anticipated in the stability analysis, presented in the next subsection.
REMARK
3. The ‘plus’ and ‘minus’ formulations
are direct extensions of the stabilized methods that we have been investigating for some years now, but apparently not introduced before for the incompressible Navier-Stokes equations. They both reduce to the SUPG formulation for the linearly interpolated velocity approximation. In this case, with an equally linear pressure, Hansbo and Szepessy [17] have made an analysis and some computations of the SUPG method, with a design of the stability parameter T, which is always 0(/z,) in their case. By assuming that the Reynolds number is always sufficiently high, they do not distinguish the locally diffusive flow cases in their choice of T. Tezduyar et al. [18] have also studied a similar formulation, making a proper distinction between locally advective-dominated and locally diffusive-dominated flows. One of the designs of T in [lS] reduces to our formulae (8)-(11) by setting mk = f and by dropping inequality (13). While this is satisfactory for linear interpolations, it breaks down for high order interpolations. A related formulation is also proposed by Lube and Auge [25]. In the next section we show that the ‘plus’ formulation is unstable for larger values of 7 than those proposed above when we employ high-order interpolations for the velocity field. In particular, for equal-order serendipity interpolation, using a ‘wrong’ value of m = 4 instead of mkz2 = 0.0814814. . . as suggested by recent computation of inverse estimates [19], the numerical simulations present spurious oscillations for both velocity and pressure in the ‘plus’ formulation (see Section 3.3).
2.2. Erroranalysis In this section we study the convergence of the ‘plus’ and ‘minus’ methods given by (4)-(13). Let us first note that the stability parameter T is bounded by a constant in each element domain K. By definition 7(x, Re&))
= 2,a;:),
2
Re,(x) 2 1 ,
(17a)
P 7(x3 Re,&N
Therefore,
for Re.(x) 7(X, Re,(x))
e.4 = -jg-
,
OaRe,(x)<
1.
(17b)
3 1, =
hK 2l+)l,
1 Re,(x)
mMx)l,h, 4~
(18) and combining with the definition (17b), we conclude that the bound (18) is valid for all values of Re,(x).
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
214
We are now ready to establish the stability characteristics LEMMA 2.1. (Stability). Assuming (I) V*a(x)=O, (II) v(x) = v = constant > 0 , then for all (u, p) E V,, X Ph , B,(u,
PROOF.
p; u, -p)
the given data a(x) and v(x) to satisfy
2 $(24IE(4Il(:
uE
ll~1’2(P4~ + VP>lli>+ w2v
+
First note that by integration
((Vu>a,u) = 0,
of each method.
by parts, by assumption
4l:, ’
(I) and by (2),
v,,,
since u=O on rfor uEV,. Next, for each formulation we have (i) The ‘plus’ formulation: By (5), B+(u, p; u, -p> = 0 + 24E(U)II:,
+ VP)ll:, - Kg
+ ll~1’2((wJ
Combining
+ llg1’2v. UlK
l17”22~v’ 44lL
(19)
.
(18) to estimate the last term of (19) gives Kzq,, 1171’22vV*e(~)Ili,~ d mky”’
2
h;llV. &(~)ll,?;,~ (by (II) and (18))
KE%;I,
(by(13)) (20)
(by (12)) * Combining this estimate with (19) yields the desired result. (ii) The ‘minus’ formulation: Let p > 1. By (5), B_(u,
p; u, -p)
= 2V((E(U)((; + ((PV* + K&
IIT
ulic: + IlT”‘((v~)a
+ Vp)ll,:
“*2Vv E(U)Il(:.K- 2 K& (2Vv’ +)7
>2V((E(U)((; + ((PV.
T((Vu)a + %))K
p-‘>~~~“‘<
z& + (l-
+ (l - P) Kzq,, ~~T1’22vv’ E(U)Il:,.K 2
(1
+
;(I
-
k9&4+411(:
+
(1
-
fi-‘)h”‘((v+
+ )(61’2v. ull; .
For the specific choice of p = 2, the stability estimate follows.
+
vP)ll:, (21)
Cl
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
215
Next, let us recall that by standard approximation theory [26], for a regular family of elements, there exist interpolants u”hlK E R,(QN and p”hlK E R,(K) such that JJU- U”hlKJJm.K S Ch~l-mlU~k+,,K
vu E Hk+‘(qN,
I(p - P”hlKllm,K =sChym(pI,+,.K
vp E HI+‘(K),
OG m d k + 1 ) 09 m f I+ 1 .
(The expression regular family of elements has the usual meaning The following interpolation estimate can now be established.
(22)
(23)
[26]).
LEMMA 2.2. Assume that the solution to (1) satisfies u E Hk”(fl)N fl H:,(R)N and p E H’+‘(O) rl L:(R). Denoting by T-J,= u”,,- u and 7, = ch - p the interpolation errors, for each KEY?,,, we have VxEK, then (a) if Re .al,
11~-“*%ll~.K +2410?uN~.K + 11~“*P?,)4I~.K + ll7”*2~v*o-L)ll~.K + IP”*v’ %ll~,K+ ll~-“*~,llL + ll~1’2v77pll:),K s Cc;~~l”l~hZKk+ll~I:.l.K (6)
if 0 S Re,(x)
+ SE; l4,‘h~+‘lPl:+l.K)
;
(24)
< 1, V’xE K, then
Il~-“z~Pll~.K +24Ie?.)llL + Il~“*
(2%
Therefore,
ll~-“*11,ll~ +241+7U)II~ + II~“2w?u)41~ + lP”*\J*sll,: + Kzs [lb“*m7-+~)ll~.~ + HO% - ~W”2~,Il:).K + H(;- Re,W-‘lh,ll&J
+ ll~“*vrlpll~
d C[ c hZ,kl&+‘,K(H(ReK - l)h, KEY&
+
sup /al,, + H(1
-
Re,)2v)
XEK
hZ,‘l~lf,‘,dHWK - W, ;Ec ld,’ + HU - &)h:@+‘)]
KTe
, (26)
h
where H( * ) is the Heuviside function given by
(27)
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
216
PROOF.
The velocity estimates related to the momentum equation follow as in Part I, noting that the obvious extensions apply for a vector field. For pressure and velocity (related to the continuity equation) estimates, let us divide the proof into two cases. (a) Let Re,(x) 3 1, Vx E K. Then 116li2v* %ll;,K +
ll~-1’2.t7,11~,R + ll~1’2Vrlpll;.K = II(M~)l,~K)1’2V* ?IuIlL
+ lI(~l~(~)l,~,)-1’217,11~.K + ll~21~~~~~~~1’2~~,:.~ -h/c ;~p,bl,llV* s.lli.K+ ;‘J; d
(28)
C sup lalphZKk+lkl:+,,R + sup l~~,‘h;+‘l~lf+,,~). ( XEK
(b) Let 0 d Re,(x) Ils1’2v’
Therefore,
l~l,l(~-‘h,‘~~~pll~,K + h,cIhpll;,K)
XEK
< 1, VX E K. Then dli,K
+
G
2 2vllv’
d
C@‘h2,kbd:+~.K
(2v)-1b?p~~~.K
+
%IIi.K
(24)-(26)
+
+
Ib1’2vt7pll;l.K
(2V)-1(b?pll:,K
+
h~llv~pll~,K)
(29)
(2v)-1h2,‘+2bt+1.K).
follows from (28), (29).
q
We may now establish the following convergence
estimate.
THEOREM 2.1, Under the same hypotheses as in Lemmas 2.1 and 2.2, the solutions (u,, , p,,) converge to (u, p), solutions of (1) as follows: of the methods given by (4)-(13)
4 i 2Wuh
-
u>lli+ ll~1’2(w% - u>a + VP,
s c KTz (h~bl~+I,K(H(R% - ‘lhK ;EF
-
b/p
P>>lI3 + lW2V*(Uh- 4lli +
HU
-
ReKP)
h
+
h2K PI~+I.AH(R~K -
l)h,
sup [al,’ + H(l - Re,)h:(2u)-l))
,
XEK
where H( * ) is defined in (27). PROOF. Let e: = u,, - kh, ei = p,, - @,, and ey = ei + Q, proof is as follows:
ep = ei + vp. The convergence
L.P. Frwma, S. I.. Frey, Stabilized finite element methods: II
217
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
218
Before continuing,
note that
$ 3 H(Re,
d
+
$
H(Re,
Combining with the previous yS = 1 f H(Re, - l), gives
- 1)//S 1’2V. eiII:_, + $1
- l)ila -i’2n,lI:.~ + ; estimate,
~(2vl(iz(e~)l)~ f l171'2((Ve:)a
and choosing
-I- Ve~)ll~) +
H(1 - Re.)2++j:)~l;,,
H(1 -
Re,)(2d-‘h,//;,,]
‘y, = 3 + 2NH(l - Re,),
.
‘yz= 4 and
1116 “*V*ef:lli
-117-1’2rf&+K~~h[f(3+2h'H(1-Re~))2vllF(r),)ll~.~ + (7 + 2NH(l - Re,))ll7”‘(V~,>all~,~
e
+ t(l f H(Re,
+-2(7 + 2NH(l-
Re,))(fr”*2vV.
C 4(3 +2N)H(l
- ReK)(2v)-‘jj77,(j&
- l))llS”2V.
E( ~,)I/~,, + H(Re,
rl,lli,K
- l)ll~-“~~~ll~~~
+ 2(7 + 2NH(l - Re,))~/~“*V~,~~~,~]
c(ll~-1’211,11~ + WI&(rl”)lli + ll~“2mubII,:+ IP *‘2v*s,lG &f~~)lli,K+ ff@eK- l>llS-1’2~p /Ii\, + & trlT1’22~V.
d C
2 (~~lul~+~,K(H(ReK - l)h,
sup 101~-t H(l - Re,)2v) XEK
KE%h
+ GIPI?+,.AWR~K - l)h, The last inequality
;EJ la/,’ -t H(l
-
(30)
Re,)hi(2v)-‘)) .
follows from Lemma 2.2, which also gives that
aP4147?u)ll: + Il~““~P~,~~ + b,)Il3 + 4ll~“2vef?,II~ s c KzT
(~~~~l~+~,K(H~ReK
-
l)h,
sup lalp + H(1-
Re,)2v)
XEK
h
+
and the result follows inequality. q
h:‘l!#+l,K(H(ReK
by redefining
-
l)hK
constants
,sF;
la/,’
+
ff(1
-
ReK)h;&‘)-I))
,
(31)
in (30) and (31) and using the triangle
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
2.3. Relationship
219
with the Galerkin method using the mini-element
Pierre [27, 281 and Bank and Welfert [29] have shown, for the Stokes flow, the equivalence between the Galerkin method employing the mini-element of Arnold et al. [30] with the stabilized methods of Brezzi and Pitkaranta [24] and Hughes et al. [12], for a particular choice of the stability parameter emanating from the specific definition of the bubble function used. Recently, Brezzi et al. [31] have indicated that for advective-diffusive equations a similar equivalence holds. They have shown that the bubbles not only help in circumventing the compatibility of space requirements, but also stabilize the advective operator, as established for a linearized model of the compressible Navier-Stokes equations. In this section we use similar arguments to show the relation between the Galerkin method employing the minielement and the SUPG method, (14)-(16), w h en applied with equal-order linear elements for velocity and pressure to the two-dimensional linearized incompressible Navier-Stokes equations. Let us start by recalling the finite element spaces for the mini-element:
vi={uE H#qN
) UJK E (P,(K)@
qzqy,K E
%,}
(32)
and
where B(K) denotes standard bubble functions (such as the cubic bubble in [30]) defined on K with the following properties:
I cj K
do=
C,,hi
= C,, meas(
llWll”,K= c3,3
$1, E B(K), K E Vi’,, ,
(34)
AK E W), KE Q% 3
(35)
where CIK, C,, and C,, are independent of the element parameter h,. Then the Galerkin method reads: find u,, E Vi and p,, E Pi such that
B,(“,>
ph;
u; 4) = &@, 4))
@, 4) E
v;’
p:, ,
(36)
with B,(u,
p;
K&
4) = (J-Yu) *
u, 4) = ((VU)&
4
+ (244,
+.J))
- (V* u, p> - (V* u, q>
(37)
and (38)
Following the demonstration in [31], the first step is to ‘statically condensate’ the bubble from the momentum equation. To this end, let us decompose the trial function for the velocity as ‘h@)
=
%(d + K& +(-+,
3
h
where ur(x) E V; = {u E HA(0))N ( u(, E P,(ZC)~, K E %,}, ub is the bubble nodal value and
AK E B(K) c fwo
Next let us take u = 4ei, on K, for i = 1, . . . , N, and zero elsewhere, where ei is the unity Cartesian vector with components (ei)j = 6, for j = 1, . . . , N. For this choice of weighting
L. P. Franca, S. L. Frey, Stabi~ize~~n~te element methods: II
220
functions
and using (39) in (36) with 4 = 0, we obtain (tV$)%
cbe,)~ + (2VE(#Ub), e(&&
In obtaining (40), we have used integration and f are constant in K) to obtain ((W&M,
+ (VP, rbe,)~ = (.L
#ei)x
(40)
e
by parts (and we have assumed that the data a, Y
(41)
(~VE(Z~I),
E(+ei))K
(42)
= 0 -
Note that
(43)
=Az+,-e,. Therefore
the ‘bubble equation’ 2vAu, *e, = -((Vu,)a
(40) becomes + Vp -f,
(44)
4ei)K .
Or, using the piecewise constant data assumption, z+, = - ;v A-‘((Vu&
+ Vp -f)l,
I, 4 d0 .
(45)
Having computed the bubble nodal value, let us examine its effect on the Galerkin method with u = ut E VL. Thus by using the decomposition (39) in (37), we obtain ‘C&z,
Ph;
‘1,
d
=
((vu,)a~
%>
+
Kzq
(@buda,
+C
+
(2VE(Ulh
E(Ul))
h
-
(V- Ul,
Ph)
- P*
U1r
(46)
4) -
where we used the fact that (2ve(4uu,), e(u,)), = 0. Note that integrating combining the two terms above under the sum, we obtain
B&h,
ph; u,7
4) = &t”,
2
Ph; ul,
d -
K;qtub+,tvulb
-‘dK
7
by parts and
(47)
k
which is just the Galerkin operator for equal-order linear elements plus a perturbation term. This additional term can be further worked out with the help of (34), (35) and (45) to yield
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
PT = - &
221
(u,+, (Vv,)a -V&C h
Kzsh $,A-‘((Vu,b+ VP-f)i.>((V+
=
=c
‘l;:k
h; A-‘((Vu&z + Vp -f,
-oq),(/K 4 da)*
(Vu,)u -V&
(48)
.
KE’Gh
Therefore the perturbation terms that emanate by condensation are in fact adding a term as in the SUPG method [13] with a stability parameter matrix defined as ClkC2khi 'mini
=
A-l
2u
REMARK
4. Up to the matrix character, the stability parameter is form-identical to the formula (17b) for locally diffusive dominated behavior (Re, < 1). However, for locally advective dominated flows, this formula has been shown to be ineffective and to recover the O(h,) for the T parameter. In this case, either one modifies the condensation procedure in ad-hoc fashion, or a new definition of the bubble function must be introduced. Both alternatives are discussed in [31]. REMARK 5. A relation between stabilized methods and the Galerkin method with bubble functions has also been noticed by Arnold [32] and Pitkaranta [33] for the Reissner-Mindlin model of plate bending. In particular, Arnold [32] shows the equivalence between the Galerkin method introduced by Arnold and Falk [34] and the Galerkin-least-squares method of Franca and Stenberg [35].
3. The incompressible Navier-Stokes equations 3.1.
Stabilized methods
Let us now consider the incompressible
Navier-Stokes
equations,
$+(VU)~-~VV.E((~)+V~=~
innX(O,T),
v*u=o
in 0 X (0, T) ,
u=g
on rg X (0, T) ,
an=h
on r, x (0, T) ,
u = 240
inRatt=O,
given by
(50)
where in addition to the notation introduced in Section 2.1, rg and r, are complementary subsets of r to accommodate given Dirichlet (g) and Neumann (h) boundary conditions, and
L. F. Franca,
222
S. L. Frey, Stabilized finite element methods: II
u is the stress tensor given by o=-pz+/e!v&(u).
Using the finite element
spaces given by (2) and (3), and defining
we can now extend the methods presented follows. Find uI, E Vf and ph E P, such that
Bty(u,, ph;
0,
q) =
in Section 2.1, in a semi-discrete
fashion,
F:Yu>4) > tu*4) E 62x Ph
as
(53)
’
with BY{% p; u, 4) =
2
+ (Vz+, u) f (22+),
E(fJ)) - (V* u, p) - (V. @>4)
u + Vp - 287. E(U), 7(x, Re,(x))((Vu)u -v~t2~v.~(u)))K+(V.u,6V.v)
(54)
and J?=%, q) = (f, r.9 + (k u)~~ + c
(f, 7(x, Re&))((V@
-Vq I: 2~7. ~(0)~
KEqk (55)
where 6 and the stability parameter
T are defined as
(57)
L. P. Frunca, S. L. Frey, Stabilized finite element methods: II
REMARK 6. The SUPG 0. e(u) = 0 in (54)-(55).
method
is obtained
by setting
223
in the weighting
function
slot,
REMARK 7. The methods above are form-identical to the methods introduced in Part I and studied in detail for the linearized set of Navier-Stokes equations in Section 2. Here, we do not attempt to perform a convergence analysis for the nonlinear model. This is done for SUPG employing equal-order linear elements in [17]. We should emphasize that all the main numerical difficulties in the simulation of the nonlinear model are also present in the linearized model of Section 2, namely the simulation of high advective flows and the compatibility of the velocity and pressure spaces. There it is shown that these stabilized methods deal with these features appropriately and simultaneously. We understand that this is an important result per se, and it is applicable to various high order interpolations, with different stability features for the different methods studied, which is confirmed numerically for the nonlinear model in Section 3.3. In the next subsection we show how the nonlinearity is dealt with as part of the general predictor multi-corrector algorithm employed in [18, 201 to march in time to steady solutions. 3.2. Time integration algorithm Discretization of (53) is carried out by expanding the trial functions { uh, ph} and {u, 4) in terms of their finite element basis or shape functions. This leads to the following set of semi-discrete equations:
[M +
MT]a +
N(u) + N”,(u) + [K + K", + K”,]u +
[M”,]a + N;(u) + [G' + K4,]u + [GCI,]p = E, ,
[G
+
Gy]p
=
F
+ F, , (63)
where u and p are the vectors of degrees of freedom for uh and p,, , respectively and a is the vector of degrees of freedom for au,, / at. The matrices [Ml, [K] and [G] emanate from time-dependent, viscous and pressure terms, respectively. The nonlinear vector N(u) is derived from the advective term. The other matrices arise from the additional T- and a-terms. To define the solution procedure to (63), let us introduce the following notation: n ,_, the maximum number of time steps, i max, the number of corrector passes through the algorithm, (Y, a parameter governing the approximation of a. Following [20] and [18], we consider a one-step generalized trapezoidal method given by the following predictor multi-corrector algorithm. ALGORITHM 3.1 Step 1. Initialize ~“0,p”, and ai. Step 2. Set n:=O and i:=O. Step 3. Predictor phase: u(i) *n+1*-
uw .n+l
‘-
PCi) n+l'_--
.;) + At(1 - a)~;) ,
0, p(i) n
*
(64)
L. P. Franca, S. L. Frey, Stabilized
224
Step 4. Form and factorize the incremental
finite element methods:
II
matrices [M*], [G*], [G*‘], [G,] defined by
[M*] := [M] + [M;] + a At[N*] ,
(UpI,)] +[Z (&,I +[K] + [K;] + [K;]
[N*] := [g [G*]:= [G] CC*‘] :=
+
)
[G;],
[M;] + a At( j[G'] +
[z (u;;,)]+[K;]).
(65)
Step 5. Form the residual vectors R(‘)
.-
s(i)
._
n+1*-
n+1--
Fm+,- ([M + iqay;,+N(uy:,) + Ny(u;j,) + [K + KY + Kfi’]z(;, + [G + G;]&) ,
F n+l +
E
7n+,
([M;]a;;, + N’f(z&) + [G' + K’f]u;;, + [G;]p;;,)
-
Step 6. Solve the incremental
.
(66)
system
(67) Step 7. Corrector
phase:
(68)
Step 8. If i < i,,,,
then i := i + 1 and return to Step 4; otherwise,
u(i)
._
$1
:= &wJ
n+1.-
n+l
P@) ??+I*-
u(i,,,)
n+l
n+l
. _
i := 0 and
p(i,,,) n+l
3
(6%
)
.
Step 9. If n < nmax, then II := y1+ 1 and return to Step 3; otherwise stop. REMARK suggestions
8. The matrices given by
defined in (65) are simplified
following Tezduyar
et al.‘s [18]
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
[iv”] := [M] [G"]
:=
[G]
)
(70)
,
:= a At[G’]
[G*']
225
.
With these simplifications, the matrices in (65) do no longer depend on uzi 1. This implies that the matrices can be formed and factorized once and for all, before entering the loops of time step and iterations. In this case, one may consider that Step 4 is performed before the time loop starts. Another advantage is that the system given in Step 6 becomes symmetric with these simplifications. 3.3. Numerical results In this section we present some numerical experiments Navier-Stokes equations using the methods described in serendipity approximations were employed for all tests performed on SUN Sparcstations l+ of LNCC using the MODULEF/INRIA and VISU2/Laval. 3.3.1.
with the nonlinear incompressible Sections 3.1 and 3.2. Equal-order reported. All computations were graphics facilities of the programs
The leaky cavity problem
The problem consists of the flow in a unit square with velocity boundary conditions: u, = 1, u,=Oaty=l(O
Next, we perform the computations for the refined mesh employing the ‘minus’ formulation with Re = 400. The number of nodes are practically the same as used by Tezduyar et al. [18] with a Ql/Ql combination of interpolations, and their results are qualitatively similar to those in Fig. 3. In Fig. 3, the flow is locally diffusive dominated, since by definition (see (58)) the element Reynolds number Re,, in the refined mesh for the flow with maximum unit velocity is equal to Re
=
K
Considering
wlww 4(1/400)
= 0.52 < 1 .
a flow with Re = 5000, locally advective dominated
(71) flows are now present and the
L. P. Franca. S. L. Frey, Stabilized finite element methods: II
226
Velocity
Stream
Horizontal
function
velocity
Pressure
elevation
at x = 0.50
0.
-0.
0. 705 P
=I
0. 428
-0.
0. 150
-0.
0.0
contours
Pressure
at P = 0.50
0. 982
-0, 127
Pressure
-0. 0.2
0.5
0.7
;Y Fig. 1. ‘Minus’ formulation: serendipity elements).
steady-state
1.0
Y leaky cavity flow for Re = 400 employing a coarse mesh (8 x 8 Q2S/Q2S-
227
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
“Minus”
“Plus”
Velocity
Velocity
ir
..,**,*,.I..
. * t + .
b
’
r.o~~t\~.,
* .
.
3 .
.
.
,
,
. .
-
.
,
.
,
, c
,
.
.
c
. .
.
. .
. .
,
I
.
.
.
.
.
. .
.
.
+
.
1
,,,,
+,t
,
sLeCC,.
.
.
c
-
..I
.
.
-
.,t.*\\.-._&_
4
..LC’/...
.
*
l
.
.
-
///I
.I\\-,
v,;
.,t\tt,.*,
.
. .
-.
.
.
c
,
,
.
-
.
*
.
(
+
\.
.
.
-
.
h.
*.
.
*
_
*.
.
..*
_
Pressure contours
Pressure contours
0 -
Pressure elevation
.
.
( +. . . . . . . . . . . . . . * . * 4
.
/
-F-C,,
c
*
.
/
.
1
n
Pressure elevation
Fig. 2. Comparison between ‘minus’ and ‘plus’ formulation with a ‘wrong ’ value of m (m = 4): steady-state cavity tlow for Re = 400 employing 8 X 8 Q2S/Q2S_serendipity elements.
leaky
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
228
Pressure contours
Velocity I
- f--?-Y-?-7+.--:-~?+-----’
. L L
.C
I...
,:
,v..
.:
.
.
.
LICCC#. .
.
.**CCLCa+CC*.
.
.
.
.,
.
.
.
.
.
.,
,
*....***.....<
*
. .
.
. .
.
.
.
.
.
. .
. .
.
.
.
.
.
. .
. .
.
.
. .
. .
.
. .
.
. .
.
.
.
.
.
Pressure elevation
Stream function
Pressure at 5 = 0.50
Horizontal velocity at x = 0.50
Ul
1. 000
0
0. 697
a
a. 394
-0
091
-a
025
-0
037
211
0.0
0.2
0.5
0.8
Fig. 3. ‘Minus’ formulation: steady-state Q2S/Q2S_serendipity elements).
1.
leaky
cavity
flow for
0.0
0.2
Re = 400 employing
0.5
0.8
a refined
1.0
mesh
(16 x 16
L.P. Franca, S. L. Frey, Stabilized finite element methods: iI
229
Pressure contours
St ream function
Pressure elevation
0 ~ I
d
-
Pressure at 2 = 0.50
l3orizontal velocity at T = 0.50
Ul
1.000
0.004
0.717
0.001
P
0.433
0.150 -oat33
-0.002 -0.005
a.0
0.3
0.5
0.8
Y Fig. 4. ‘Minus’ formulation: elements.
steady-state
1.0
-0.008
r
0.0
0.3
0.5
0.6
rlo
Y leaky cavity flow for Re = 5000 employing16 X 16 Q2S/Q2S+erendipity
230
L. P. Franca, S. L. Frey, Stabilized finite element methods: II
ihl
611
Fig. 5.
The backward step flow: (a) problem statement; (b)
mesh employed
(224 Q2S/Q2S_serendipity
elements).
results employing the minus formulation are shown in Fig. 4. (We take A = 1 in the 6 formula, (56).) The results obtained are similar to those found in [36]. 3.3.2. The backward step flow The problem statement and the mesh employed are shown in Fig. 5. We take Re = 60 with respect to the inflow section. For the Q2S element, this is basically a diffusive dominated flow, by computing Re, similarly as in (71). The results for the ‘minus’ and ‘plus’ formulations are shown in Fig. 6 and 7, respectively. In Fig. 8, the ‘minus’ formulation is tested with a free-traction outflow condition. Note that the results are similar to the Dirichlet outflow test, except for a region near the outflow. In all computations, we observe the recirculation at the step, which is absent for a Stokes flow.
Fig. 6. ‘Minus’ formulation for the steady-state elevation; (c) pressure contours.
backward
step flow with Re = 60: (a) velocity
vectors;
(b) pressure
L.p.
Franca, S. L. ‘Frey, Stabilized finite element methods: 11
Fig. 7. ‘Plus’ formulation of the steady-state elevation; (c) pressure contours.
231
backward step flow with Re = 60: (a) velocity vectors; (b) pressure
Fig. 8. ‘Minus’ formulation for the steady-state backward step flow with Re = 60 using free-traction condition: (a) velocity vectors; (b) pressure elevation; (c) pressure contours.
as outflow
232
L. P. Franca,
S. L. Frey, Stabilized
finite element methods:
II
Acknowledgment The authors wish to thank Tayfun Tezduyar for fruitful discussions and useful suggestions and Michel Fortin, Michel Lapointe and Roger Pierre from Lava1 University for the graphics post-processing code VISU2. During the course of this work Leopold0 P. Franca has been partially supported by CNPq Procs. 300419/90-2 and 402598/90-3.
References [l] I. BabuSka, The finite element method with Lagrangian multipliers, Numer. Math. 20 (1973) 179-192. [2] F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from Lagrange multipliers, RAIRO Ser. Rouge 8 (1974) 129-151. [3] 0. Pironneau, Finite Element Methods for Fluids (Wiley, New York, 1989). [4] R.L. Sani, P.M. Gresho, R.L. Lee and D.F. Griffiths, The cause and cure (!) of the spurious pressures generated by certain FEM solutions of the incompressible Navier-Stokes equations: Part l., Internat. J. Numer. Methods Fluids 1 (1981) 17-43. [5] R.L. Sani, P.M. Gresho, R.L. Lee, D.F. Griffiths and M. Engelman, The cause and cure (!) of the spurious pressures generated by certain FEM solutions of the incompressible Navier-Stokes equations: Part 2, Internat. J. Numer. Methods. Fluids 1 (1981) 171-204. [6] A. Soulaimani, M. Fortin, Y. Ouellet, G. Dhatt and F. Bertrand, Simple continuous pressure elements for two- and three-dimensional incompressible flows, Comput. Methods Appl. Mech. Engrg. 62 (1987) 47-69. [7] C. Taylor and P. Hood, Numerical solution of the Navier-Stokes equations using the finite element technique, Comput. & Fluids 1 (1973) l-28. [8] R. Verfiirth, Finite element approximation of incompressible Navier-Stokes equations with slip boundary condition, Numer. Math. 50 (1987) 697-721. [9] F. Brezzi and J. Douglas, Stabilized mixed methods for Stokes problem, Numer. Math. 53 (1988) 225-236. [lo] J. Douglas and J. Wang, An absolutely stabilized finite element method for the Stokes problem, Math. Comp. 52 (1989) 495-508. [ll] T.J.R. Hughes, L.P. Franca, A new finite element formulation for computational fluid dynamics: VII. The Stokes problem with various well-posed boundary conditions: Symmetric formulations that converge for all velocity/pressure spaces, Comput. Methods Appl. Mech. Engrg. 65 (1987) 85-96. [12] T.J.R. Hughes, L.P. Franca and M. Balestra, A new finite element formulation for computational fluid dynamics: V. Circumventing the BabuSka-Brezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations, Comput. Methods Appl. Mech. Engrg. 59 (1986) 85-99. [13] A.N. Brooks and T.J.R. Hughes, Streamline upwind/Petrov-Galerkin formulations for convective dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1982) 199-259. [14] C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method (Cambridge Univ. Press, Cambridge, 1987). [15] T.J.R. Hughes, L.P. Franca and G.M. Hulbert, A new finite element formulation for computational fluid dynamics: VIII. The Galerkin-least-squares method for advective-diffusive equations, Comput. Methods Appl. Mech. Engrg. 73 (1989) 173-189. [16] L.P. Franca, S.L. Frey and T.J.R. Hughes, Stabilized finite element methods: I. Application to the advective-diffusive model, Comput. Methods Appl. Mech. Engrg. 95 (1992) 253-276. [17] P. Hansbo and A. Szepessy, A velocity-pressure streamline diffusion finite element method for the incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 84 (1990) 175-192. [18] T.E. Tezduyar, R. Shih, S. Mittal and S.E. Ray, Incompressible flow using stabilized bilinear and linear equal-order-interpolation velocity-pressure elements, Report UMSI 90/165, University of Minnesota/Supercomputer Institute, Minneapolis, 1990.
L. P. Franca, S. L. Frey, Stabilized finite element methods: II [19] I. Harari and T.J.R.
233
Hughes, What are C and h?: Inequalities for the analysis and design of finite element methods, Comput. Methods Appl. Mech. Engrg. 97 (1992) 157-192. [20] T.E. Tezduyar and T.J.R. Hughes, Development of time-accurate finite element techniques for first-order hyperbolic systems with particular emphasis on the compressible Euler equations, Report prepared under NASA-Ames University Consortium Interchange No. NCAZ-OR745104, 1982. [21] P.M. Gresho, Incompressible fluid dynamics: Some fundamental formulation issues, Annu. Rev. Fluid Mech. 23 (1991) 413-4.53. [22] T.E. Tezduyar, M. Behr and J. Liou, A new strategy for finite element computations involving moving boundaries and interfaces - The deforming-spatial-domain/space-time procedure: 1. The concept and the preliminary numerical tests, Comput. Methods Appl. Mech. Engrg. 94 (1992) 339-351. [23] T.E. Tezduyar, M. Behr, S. Mittal and J. Liou, A new strategy for finite element computations involving moving boundaries and interfaces - The deforming-spatial-domain/space-time procedure: II. Computation of free-surface flows, two-liquid flows, and flows with drifting cylinders, Comput. Methods Appl. Mech. Engrg. 94 (1992) 353-371. [24] F. Brezzi and J. Pitkaranta, On the stabilization of finite element approximations of the Stokes problem, in: W. Hackbusch, ed., Efficient Solutions of Elliptic Systems, Notes on Numerical Fluid Mechanics, Vol. 10 (Viewig, Wiesbaden, 1984) 11-19. [25] G. Lube and A. Auge, Regularized mixed finite element approximations of incompressible flow problems. II. Navier-Stokes flow, Otto von Guericke University, Preprint 15/91, 1991. [26] P.G. Ciarlet, The Finite Element Method for Elliptic Problems (North-Holland, Amsterdam, 1978). [27] R. Pierre, Simple C” approximations for the computation of incompressible flows, Comput. Methods Appl. Mech. Engrg. 68 (1988) 205-227. [28] R. Pierre, Regularization procedures of mixed finite element approximations of the Stokes problem, Numer. Methods Partial Differential Equations 5 (1989) 241-258. [29] R.E. Bank and B.D. Welfert, A comparison between the mini-element and the Petrov-Galerkin formulations .for the generalized Stokes problem, Comput. Methods Appl. Mech. Engrg. 83 (1990) 61-68. [30] D.N. Arnold, F. Brezzi and M. Fortin, A stable finite element for the Stokes equations, Calcolo 23 (1984) 337-344. [31] F. Brezzi, M.-O. Bristeau, L.P. Franca, M. Mallet and G. RogC, A relationship between stabilized finite element methods and the Galerkin method with bubble functions, Comput. Methods Appl. Mech. Engrg. 96 (1992) 117-129. [32] D.N. Arnold, Innovative finite element methods for plates, Report No. AM 61, Department of Mathematics, PennState, 1990; to appear in Mat. Apl. Comput. (331 J. Pitkaranta, Analysis of some low-order finite element schemes for Mindlin-Reissner and Kirchhoff Plates, Numer. Math. 53 (1988) 237-254. [34] D.N. Arnold and R.S. Falk, A uniformly accurate finite element method for the Reissner-Mindlin plate, SIAM J. Numer. Anal. 26 (1989) 1276-1290. [35] L.P. Franca and R. Stenberg, A modification of a low-order Reissner-Mindlin plate bending element, in: J.R. Whiteman, ed., The Mathematics of Finite Elements and Applications VII, MAFELAP 1990 (Academic Press, London, 1991) 425-436. [36] P.M. Gresho and S.T. Chan, On theory of semi-implicit projection for viscous incompressible flow and its implementation via a finite element method that also introduces a nearly consistent mass matrix. Part 2: Implementation, Internat. J. Numer. Methods. Fluids 11 (1990) 621-659.