- Email: [email protected]
Algorithmic and theoretical results on computation of incompressible viscous flows by finite element methods
Computers& FluidsVol.13,No. 3,pp. 361-373,1985
0045-7930/85 $3.00+ .00 © 1985PergamonPressLtd.
Printedin the U.S.A.
ALGORITHMIC
AND
COMPUTATION FLOWS
BY
OF
THEORETICAL
RESULTS
INCOMPRESSIBLE
FINITE
ELEMENT
ON
VISCOUS
METHODS
M. D. GUNZBURGER and R. A. NICOLAIDES Department of Mathematics, Carnegie-Mellon University, Pittsburgh, PA 15213, U.S.A. and C. H. LIU Low Speed AerodynamicsDivision, NASA-LangleyResearch Center, Hampton, VA 23665, U.S.A. (Received 14 December 1984; in revised form 12 February 1985)
Abstract--Primitive variable as well as streamfunction-vorticity and pure streamfunction formulations are discussed. For the primitive variable case alternative choices of the viscous stress term are shown to produce natural boundary conditions which are well suited for matching to various far field conditions. For the other cases recent analytical results, including error estimates are described, and an optical algorithm for pressure recovery as well as treatment for multiply connected domains are given. I. INTRODUCTION Many different weak formulations are used as a basis for finite element solutions of the incompressible Navier-Stokes equations, ranging from fourth order streamfunction techniques through second order streamfunction-vorticity methods, velocity-vorticity-pressure systems and, most important perhaps, the standard primitive variable approach. For a survey, see [1]. Of these methods, the primitive variable technique has been most thoroughly analyzed [2, 3] although the streamfunction-vorticity method (finite element style) has also received attention [3-8]. Most of the analysis deals with interior flows driven by body forces. In applications, external flows are also important, and in particular the treatment of far field boundary conditions is often a major difficulty. In such cases, the natural boundary conditions generated by finite element schemes often have a valuable role to play (see e.g. [9]). Even for a given system of variables the natural boundary conditions can vary according to how the viscous and convection terms are defined. In this paper we shall introduce some new results for the streamfunction and streamfunction-vorticity formulations, including error estimates, optimal recovery of the pressure field, the question of multiply connected regions, and also possible choices for natural boundary conditions. For the primitive variable case three types of natural boundary conditions will be considered arising from different formulations of the viscous terms in the equations. Additionally, the effect on the natural boundary conditions of alternate treatments of the convection terms is considered, along with the manner in which a parabolized Navier-Stokes system interacts with the natural boundary conditions in general. 2. PRIMITIVE VARIABLE FORMULATIONS 2.1. For stationary incompressible flows the governing equations are u. Vu + XTp = div a,
(2.1)
where u is the velocity field, p the pressure and a is the viscous stress tensor, defined by o- = P(VI! + VuT). 361
(2.2)
362
M.D. GUNZBURGERet al.
The constant density is incorporated into both p and a and of course we also have div u = 0.
(2.3)
In view of (2.3), whenever v is constant, div a = div (v(Vu + VuT))
(2.4)
= .(Au)
(2.5)
= - v curl curl u.
(2.6)
Although in the strong sense these are equivalent, in the weak sense (2.4)-(2.6) each yield distinct natural boundary conditions which we shall discuss. Let ~ denote a bounded domain of R 2, not necessarily simply connected, with boundary I'. Standard Galerkin technique yields a weak problem of the form
a(u, v) + b(p, v) + c(u, u, v) = (f, v)
b(q, u)
(2.7) (2.8)
= O,
where (u, p) denotes the sought velocity-pressure variables and (v, q) are test functions. Below it will be necessary to consider several choices for spaces containing these functions. In all cases
(f, v) = fn f" v, and
b(q, v)
=
-fq
div v.
(2.9)
However, (2.4)-(2.6) give, respectively,
ff.i(Vu+VAT). Vv a(u, v) =
v
(2.10)
~Tu. Vv
(2.1 1)
v f curl u curl v. dn
(2.12)
(We note that henceforth, since we are dealing with domains in R 2, curl u will denote the single non-trivial component of the curl of u, i.e. curl u is a scalar.) We note at the outset that although (2.12) has some useful properties associated with it, it suffers from the serious drawback of requiring the use of divergence free elements. There is also some flexibility in the choice of c(u, u, v). The first possibility is the obvious one, c(u, u, v) = f~ (u- ~Tu). v.
(2.13)
The second is the skew-symmetrized form [2] c(u, u, v) = ~
(u. Vu). v - (u. Vv). u
.
(2.14)
We note that (2.14) can be used, and indeed is equivalent to (2.13), only when div
Computation of incompressible viscous flows
363
u = 0 and either u. n = 0 or u is specified on t2. Otherwise a boundary integral must be added to (2.7). See (2.18) below. Associated with these weak forms are certain essential and natural boundary conditions. These can be deduced from the boundary integrals omitted during the formation of (2.9)(2.12) and (2.14). Specifically we have respectively r v. [pn - v(Vu + VuT) • n]
r v.
[pn -- vVu. n]
f r [pv. n - v curl u(v. r)],
(2.15)
(2.16)
(2.17)
where n and ~- respectively denote the unit outer normal and unit tangent vectors to F. If (2.14) is used, to each of these must be added
1f [ (u. n)(u. v).
2 Jr
(2.18)
For every case (2.15)-(2.18) specification of u . n and/or u . 7 , i.e. u, implies that v . n l r = 0 and/or v. fir = 0. For example, if Ulr is given, then Vlr = 0 and all the boundary integrals (2.15)-(2.18) vanish. In this case, there is no obvious particular advantage of one weak form over any other. On the other hand, for other boundary conditions some of the natural conditions can be used to advantage. Unfortunately (2.1 l) which is the most widely used formulation, has the disadvantage that the associated natural condition, which from (2.16) is t a n t a m o u n t to specifying - p n + vVu. n appear to have no physical significance, and in particular is not the stress vector at the boundary. To have specification of the stress vector - n p + v(Vu + VuT) • n be a natural condition, formulation (2.10) with (2.13) should be used if the stress vanishes. Otherwise (2.15) appears as data. Use of (2.10), (2.13) and (2.15) also allows, for example, the use of the normal velocity as an essential condition paired with the tangential stress as natural. Physically, being able to specify stress components is often necessary at interfaces between different fluids and in some free boundary calculations. Another possibility is to use (2.12), (2.13) and (2.17). Here, the natural conditions involve p and ~" -- curl u. For example, both m a y be specified as natural conditions or we m a y pair p with the tangential velocity or ~" with the normal velocity. Such combinations are useful for blending fully viscous flows with both inviscid and boundary layer type flows. This is discussed in Section 2.2. Following the remarks above we see that (2.14) can best be used when u. nit = 0 or u specified on all boundaries so that V[r = 0. Otherwise it is necessary to use (2.13) since no physical meaning seems to attach to (2.18). 2.2. In this section we shall discuss in more detail the use of (2.12), (2.13) and (2.17) to simulate far field conditions at a computational boundary. Specifically three cases that are practically important are inflow, outflow and matching to given inviscid flows. Consider first the inflow conditions. In some cases, e.g. channel flows, both components of u m a y be accurately known and so they can be specified as essential boundary conditions. However, in other cases, for example with a boundary layer inflow, the
364
M.D. GUNZBURGERel al.
boundary layer equations correctly determine u, whereas even in the simplest case of a flat plate the Blasius solution for v does not match the inviscid flow as y --~ ~ [10]. Thus, some method for imposing the inflow condition which does not involve knowledge of v would be preferable. It is natural to impose u as an essential condition, and, in view both of the inviscid flow at great distances from the boundary and our previous discussion of natural boundary conditions, to take as the other specified variable to be the vorticity ~. Then both u and ~" can be made arbitrarily close to their inviscid flow values by making the computational boundary sufficiently distant. To impose the vorticity condition, since ~"is not zero but is computed from the boundary layer equations, e.g. the Blasius solution, it is necessary to incorporate the boundary integral (2.17) into the weak form of the equations after substituting v. n = 0 and the known function curl u = ~'. This is illustrated in a specific case in subsection 2.3. Turning next to computational boundary segments which are nominally in inviscid flow regions, typically the known inviscid solution is imposed. One method is simply to specify u. An alternative is to specify p and ~" or some combination of these two possibilities. Usually the p, ~"choice is the preferred one because these variables approach their asymptotic limit faster than does u [ 11 ]. Once again, (2.17) with the known values p and ~" substituted must be restored to the weak form (2.12). Again, see Section 2.3 below. Next we consider the notoriously difficult question of finding satisfactory outflow conditions. Ideally these should have no observable upstream influence. Were the situation one of simple convection-diffusion, satisfactory conditions can be fairly easily given in most cases [12, 13]. However the fact of incompressibility with its attendant global coupling of the pressure field renders the Navier-Stokes case rather more difficult. Often, for example for wakes, boundary layer or channel flows, the outflow has the property that the diffusion in the streamwise direction is negligible compared with both the convection and the cross-stream diffusion. Mathematically, if the x-axis is aligned with the flow as we shall assume in the sequel, then viscous terms involving x-differentiations may be neglected. This observation is systematically exploited in the parabolized NavierStokes equations (PNS) as derived, e.g., by [14, 15]. Instead of imposing outflow conditions for the full Navier-Stokes equations, we shall assume that in the vicinity of the outflow boundary, the flowfield is described by the PNS equations. Thus the boundary conditions for PNS equations are what determine the outflow conditions. For example, for the PNS equations (2.12), (2.17) are replaced by
u fe (curl u)vl,y
v = (Vl, v2)
(2.19)
and
frP
V. n.
(2.20)
Thus, only p at the outflow must be specified. No vorticity or velocity is required. The need to impose the pressure arises from its previously mentioned global character. 2.3.
Here we give a specific formulation of a problem illustrating the above considerations. We consider boundary layer flows over obstacles of size comparable with the boundary layer thickness. A simple example of such a flow region is shown in Fig. 1. Here AOBC represents a solid boundary at which u = 0. AE is an inflow boundary at which a boundary layer profile is imposed. We assume ED is sufficiently far removed from the solid boundary so that it may be regarded as being in the inviscid regime. CD represents the outflow boundary. We shall assume that the inflow on AE is given by the Blasius profile, with a leading edge at x = Xo < xA. Following the discussion of the previous section, u and ~"will be prescribed (i.e. [10])
Computation of incompressibleviscousflows
365
D
E
B
C
t,y 0
A
J-
x
Fig. 1. Computationaldomain for flowover a step.
on AE
u~. = U ~ f ' ( n )
and fi~ = U o ~ f " ( n ) [ U o o / v ( x A -
Xo)]1/2
on AE,
where ~ = y [ U o o / v ( x A -- x0)] ~/2 and f is the Blasius function. Since on ED the flow is effectively inviscid, we may use the potential flow solution to get a pressure boundary condition and set ~"= 0. In general, the potential solution for p must be found numerically (e.g. by a panel method), although in the step case being considered a closed form solution can be found from the Schwartz-Christoffel transformation. In fact, on ED, we have for a step of height d, i.e. YB = d, 1 [ ( ~ + l ) ( f f + 1)11/2 Ppot = P~o + ~ ~'2UL -- 1)(w-- ~ J ~'~o, = O,
where w is obtained from z = x + i y e by
Z =
d
--
[(W 2
_ 1)l/2 - - cosh -I w ] + i d
and where we recall that we have absorbed the constant density into the pressure variable. We note that lim ~'in~
~'pot =
0
for
YE < ~ "
Y~YE
However, the discrepancy can be made arbitrarily small by making Y e sufficiently large. In practice, f ~ decays sufficiently rapidly for this not to be a practical issue. For example i f R e = Uood/v = 1000, then with YE = 4yB = 4d we have ~'i,(YE) "" l O - S U o ~ / d . Consider next the outflow DC. Through the use of PNS only p needs to be specified. Unfortunately in most cases p is either unknown or difficult to find. Two possible choices are (a) to use experimental data or (b) to assume that the flow at DC is of boundary layer type so that the pressure is constant on DC and is given by pou,(y) =
p~,(x~, yz).
To conclude, we shall now give the appropriate weak form incorporating this particular set of boundary conditions, and the PNS equations in the vicinity of the outflow. Let
366
M.D. GUNZBURGERet al. H~([2) : {u E
n'(f01UlAo~c -- 0,
UIAE
-----
Uin}
H6(f~) = {u ~ H'(f~)luLAosc = 0, UIAE = 0}, where H1(t2) is the usual Sobolev space. T h e n seek u ~ Hel(~2) and p @ L2(f/), f~ v curl U(l)l,y -- O/V2,x) q- fi2 (U'VU)" V --
~t p d i v v
-- ~AE V~'in(V'T)-- fED Pp°tV'n--
fnq h div u h =
fcDP°mV'n
for all
0
qh E
for all v E H6(~2),
L2(f~).
(2.21) (2.22)
In (2.21) a = 1 in the main flow field, and a = 0 near the outflow boundary with, if desired, a smooth transition. This simply activates the PNS approximation at appropriate field points. N o rigorous analysis of (2.21)-(2.22) is available. However in practice element pairs suitable for wall bounded flows in the standard formulation appear to be satisfactory in the above described case. See, e.g. [1-3, 16-21], for examples of such element pairs and for a discussion of error estimates.
3. S T R E A M F U N C T I O N - V O R T I C I T Y
FORMULATIONS
3.1. The stationary Navier-Stokes equations, written in terms of the streamfunction ~b and the vorticity ~', are given by Aft = --~'in ~2
(3.1)
and
where in (3. l) and (3.2) v is the kinematic viscosity, f the body force, and 12 is a bounded region in R 2. If ~2 is multiply connected, we denote by Fo its exterior boundary and by Fi, i = 1. . . . . m, the remaining parts of the boundary. (See Fig. 2.) m
Suppose that at these boundaries the velocity u is specified, i.e. u = g on F = U I?i.
i=0
I"
0
Fig. 2. A multiply connected domain ft.
367
Computation of incompressible viscous flows In order for a streamfunction to exist, we must have that [3] for
fr g.n=0
i=0 .....
(3.3)
m,
i
i.e. there is no net mass flow through any o f the b o u n d a r y pieces. It is well k n o w n that the streamfunction is unique up to an additive constant, and we fix its value by specifying it to be zero at an arbitrary point x0 on £o. Now, let q denote a function such that
cgq Or
g . n on £
q(xo) O.
and
(3.4)
T h e n the b o u n d a r y conditions for the system (3. l) and (3.2) are given by
~=q+aionri,
=qon£o,
i = 1. . . . .
m
(3.5)
and &b -r = -g'r
on £,
On
(3.6)
where ai, i = l . . . . . m are constants to be determined as part o f the solution. These constants are fixed by the requirement that the pressure be single valued, i.e. the change in the pressure p a r o u n d any closed curve surrounding any o f l~i, i = l . . . . . m, is zero. Since the m o m e n t u m equation stipulates that 1 -
Vp = - u . Vu + vAu + f,
(3.7)
P where a is the constant density, we have that, for i = l . . . . .
m,
o = ~ f r V( P~ 'Aru=-fur' V u + f ) ' r ' ~ Since A u - z = 0f/0n whenever div u = 0 and u . Vu. r = f u - n + O[(u. u)/2]/Oz, we then have that
fr
(V Of/an - f g . n + f- r) = 0
for
i = 1. . . . .
m.
(3.8)
i
T h u s (3. l), (3.2), (3.5), (3.6) and (3.8) are the governing equations and side conditions which are to determine the functions ~k and f and the constants ai, i -- 1. . . . . m. As in the primitive variable case, other b o u n d a r y conditions are o f interest. For example, if all or part o f F0 is a c o m p u t a t i o n a l b o u n d a r y in the inviscid potential flow regime one m a y want to specify f = Of/On = 0 there, or if part o f Fo represents a b o u n d a r y layer inflow, one m a y want to specify ~b and ~. Here we will only consider the b o u n d a r y conditions (3.5), (3.6) and (3.8). Other b o u n d a r y conditions m a y also be i m p l e m e n t e d in a straightforward m a n n e r (see [5]). The weak formulation o f the streamfunction vorticity equation requires the use o f some n o n s t a n d a r d function spaces. For details concerning these spaces see [3-5]. Here we will proceed formally. The most direct weak formulation is to seek ~b E Se, f E Z such that
fn fco-fv~b. Vo~=frO~g.r -v
Vf . V~P+
f
0tP
Ox o~b Ox
=-
for all ~o E Z f-curltP
forall~P~S0.
(3.9) (3.10)
368
M.D. GUNZBURGER el aL
Here Z = Hl(fl), Se is a subset of HI(f/) whose elements satisfy (3.5), and So is a subspace of Hl(~2) whose elements satisfy (3.5) with q = 0. In deriving (3.9) we have replaced the boundary integrals
o~ f r co ~nn
(3.11)
with the boundary integral on the right hand side of (3.9). Likewise, in deriving (3.10), the boundary integral 0~"
0ff
r~
(3.1 2)
!
appears. However, note that ~P = 0 on F0 and ~ = constant on I?i so that, with (3.8), (3.12) vanishes. (For details see [3-5].) We note that no boundary conditions on ~"need be contrived on solid walls. Also, the conditions (3.8) are natural to the weak formulation, as is (3.6) while (3.5) is an essential condition. One may view (3.9), (3.10) as a mixed method for the fourth order streamfunction equation. Other mixed methods for the linear case are discussed in, e.g. [3, 22, 23], and for the nonlinear case in [3]. The particular method given by (3.9), (3.10) is discussed, in the linear case, in [23-26] and for the nonlinear case in [3-5]. With regard to finite element discretizations of (3.9), (3.10), one is particularly interested in using merely continuous finite element spaces for both ~b and ~. A sample error estimate is given in the following result [24]. THEOREM 3. l Suppose piecewise continuous linear finite element spaces are used to approximate both ~k and L Then there exists a constant C independent of the grid size h such that I[Vff - XTffhl[0< Ch
K' II0 ~
and
Ch ~n
Here ~kh and ~-h denote the discrete solution. These estimates display the general trend, i.e. the gradient of streamfunction, and therefore the velocity field, is optimally approximated, while there is a loss of accuracy in the vorticity approximation. We remark that if one uses continuously differentiable finite element spaces for ~k and ~, then both ~k and ~"are optimally approximated [24]. 3.2.
The main difficulty which multiply connected domains present is that the streamfunction may be arbitrarily specified on only one part of the boundary. Thus, one sets and ~bh = 0 at one point on Fo, but these are determined only up to unknown constants ai and a~, respectively, on the remaining parts Fi, i = 1. . . . . m, of the boundary. In the discretization of (3.9), (3.10), the constants a h are determined as part of the solution process. However, that method requires the use of N+M,
E
%(x),
i = 1. . . . .
m
(3.13)
j=N+Mi-j+I
as basis functions associated with the a~, i = 1. . . . . m. Here tpj, j = 1. . . . . N represents basis functions associated with interior unknowns, and ~Pj, j = N + Mi-i + 1. . . . . N + Mi, represents basis function associated with unknowns on Fi, i = 1. . . . . m, with M0 = 0. In general, these basis functions are only "semi-local" in the sense that they couple all the points on a boundary F~. In particular, in the discrete set of equations, a h is coupled with all the unknowns associated with elements whose closure intersects with Fi. The result of this is that the bandwidth of the linear systems encountered will be greatly increased, resulting in both a computer storage and time penalty over simply connected domain problems. Some of this penalty may be mitigated by employing a numbering
Computation of incompressible viscous flows
369
scheme resulting in banded-bordered matrices wherein the onerous couplings are found near the bottom and the right of the matrices encountered. This, of course, requires the development of special programs to solve the linear systems. An alternative to using basis functions such as (3.13) is to guess the value of the constants ai, i = 1. . . . . m, appearing in (3.5). Then one may solve for ~ and ~" from (3.1), (3.2), (3.5) and (3.6). However, in general, (3.8) will not be satisfied. At this point one m a y change the guesses for the a~, repeating the process until convergence is achieved, i.e. until (3.8) is satisfied. In more detail, we proceed as follows. Given guesses a~, i = 1, . . . . m, k E Z +, we define ~kh E Sk and ~-h C Z h [where
ZhCZ and
s~c& and where now the constants ai = a/*, i.e. are assumed known] to be the solution of the discrete forms of (3.9), (3.10) where now (3.10) holds for all ¢ph E S0h = S h fq H~(~2). The discrete equations and unknowns may be ordered in such a way that the bandwidth of the discrete system is no larger than that obtainable in an analogous discretization of a problem posed in the simply connected domain bounded by Po- Having obtained ~-h from a~, i = 1. . . . . m, we m a y substitute into (3.8). In practice, we would evaluate the integrals appearing in (3.8) by a numerical quadrature formula. Let us denote such a n approximation to each equation in (3.8) by Ii(ak), i = 1. . . . . m, which in general will not be small. O f course, each Ii(. ) is a nonlinear function of a k = (a~ . . . . . akin) through the discrete vorticity ~-h. We now need an algorithm to generate new g u e s s e s a k+l, from which ffkh+l and ~'h+1 m a y be computed and for which Ii(a k+l) is closer to zero than was Ii(ag). The most practical way to accomplish this is to use a secant type method. For example, if m = 1, i.e. f / i s doubly connected, we m a y then define a k+' = a~ -
a~ - a~ ' ii(alk). Ii(a~) -- Ii(a~ -l)
For m > be used. as (3.14) practical
(3.14)
1, some generalized secant method in R m such as the Wolfe secant method can (See [4, 5] for details about such methods.) Note that the use of formulas such requires two starting guesses a ° and a ~. We also note that in general it is not to use Newton's method to update a k since evaluating the Jacobian of I(a) = ( I 1 , - . - , Ira) requires the solution of linear problems of the same size as those for qjh and ~-h for each pair O~lOai, Od//Oai. Thus, at the price of solving a sequence of simpler problems, one may avoid the bandwidth problems engendered by the use of basis functions such as (3.13). Whether or not the alternate method is more efficient depends, of course, on how m a n y of the simpler problems must be solved as well as on how m u c h cheaper it is to solve the simpler problems. Both of these factors, for a particular geometry, will depend on the n u m b e r of degrees of freedom, i.e. the gridsize, and the Reynolds number. Therefore it is clear that the performance of iterations such as (3.13) play a central role in the overall efficiency of the alternate method of treating multiply connected domains. Unfortunately, the sequence a k generated by secant type methods are not guaranteed to converge for arbitrary values of the initial guesses, especially at high values of the Reynolds n u m b e r (for details see [4, 5]). Similar approaches for treating multiply connected domain may be used for the pure streamfunction formulation. 4. STREAMFUNCTION FORMULATION 4.1.
The differential equation is the well-known
M.D. GUNZBURGERet al.
370
and natural boundary conditions are obtained after two formal integrations with a test function ×; thus,
o
o×
og,
,, fo ,x,ax +,, fr ~, (a~,)×-,, fr (a~,) T + fo ,X~,(xx~,y-xy~,~)- f r x ~ ,X~,= f o b , O 0 when ~nn and ~rr represent outer normal and tangential directional derivatives. For an interior problem X and X, will be zero on F, the boundary of the domain 12, and so each boundary integral vanishes. We are assuming in this section that u, the velocity vector, in specified on all boundaries. For multiply connected domains, ~b is specified on only one bounding surface and is determined up to an arbitrary constant on the others, which of course must be found as discussed for the streamfunction-vorticity formulation above. The natural conditions are thus
= -~
=0
and ('0__~+ (u. n)~" = 0.
On Clearly then, these give on F
at" i-=~=o. These are useful for matching to an irrotational external flow when computing flows for which the flow domain includes both rotational and irrotational zones. For an example of such a case see Section 2. Of course it is also possible to make combinations of natural and essential conditions, e.g. we may specify the normal velocity as an essential condition and the vorticity as natural. In any case the weak problem is the following: find ~b @ H~(f~) such that
u f, A~Ax+Ax(xxg/y-Xyd/x) = fafx+ N(g, X)
VXE
H2e,o(12),
(4.1)
where N(g, ×) represents whatever boundary integrals are kept, with data represented by g, and He2(~2) is a subset of H2(~) satisfying the imposed essential boundary conditions and elements of He2,o(ft) satisfy homogeneous essential conditions. If g = 0, i.e., the data for the natural boundary conditions is zero, then N ( . , . ) is identically 0. To approximate ~b we proceed in the standard way to find ~bh ~ S* C H~2(f~), where now fl is a bounded domain of R 2. Equation (4. l) is replaced by the following find ~bh E S h such that v x h ~ Soh where Soh is the analogous subspace of H2,o(12). The basic error estimate for this procedure is given in the following theorem. THEOREM 4.1 A constant C exists, independent of Seh, such that for sufficiently large u the following estimate is valid 11¢ - ¢h112 ~ C inf I1~ - ~hll=.
¢*cs,*
For the proof of this result and details of the constants see [27].
Computation of incompressible viscous flows
371
Another weak formulation for the streamfunction formulation is given in [28]. The present formulation appears more direct however. Theorem 2.1 applies, of course, to conforming finite elements. Here, the element choices are somewhat limited by the requirement of C 1 inter-element continuity. The lowest degree conforming finite element on triangles is the Clough-Tocher cubic macroelement. A possible alternative is the Argyris 21 degree of freedom quintic elements, free of course, of the need for special macro constructions. However it is desirable to use the lowest order of polynomials possible for several reasons. First, apart from exceptional circumstances, the exact solution is not in a sufficiently high order Sobolev space for the full approximation capabilities of the quintics to be realized. Second, the costs of the numerical integration for the assembly process are determined by the need to integrate 6th degree polynomials exactly for the quintic case but only 2nd degree polynomials for the cubic macroelement case. For these reasons we recommend the Clough-Tocher element. In this case we obtain from Theorem 2. l the following error estimate: I1¢ - ¢h112 ~ C h 2 - S l l ¢ l h - s
s = 0, 1.
(4.2)
In general, however, we must expect to have to take s > 0 as shown by results of [29]. However for a square domain with f E L2([2) we can take s = 0, although even in this relatively favorable case we cannot improve on (4.2) with the quintic element. Equation (4.2) used in conjunction with the standard duality method gives estimates in lower norms, for example, we have the estimate (4.3) for a square domain. This of course gives an O(h 3) rate for [[u - uh[[o in the most favorable case, where u and u h are the velocity fields defined from ~b and ~bh. In practice, numerical integrations are used to evaluate all contributions to stiffness matrices and right hand side terms. One valuable observation is that one needs to evaluate only the viscous, i.e. biharmonic, terms exactly, in order to preserve the estimates (4.2), (4.3), provided that the same quadrature rule is used for all the other terms. For example, for the Clough-Tocher element, it is sufficient to use a rule exact for quadratic polynomials, even though the convection term involves 5 th degree polynomials. This has been verified by the authors in practical computations and anyway is predicted at least for linear problems by finite element theory [30]. 4.2.
A pressure field, of optimal accuracy for the given velocity approximation may be computed (without specifying any explicit boundary conditions) from the primitive variable m o m e n t u m equation, by discretization of the equation
-re p div v = g(v)
V v E V,
(4.4)
where g(. ) denotes the linear functional
g(v)=f(v)-ufVu. Vv-f(u.Vu)v
u = curl ~b
(4.5)
which is bounded on H~([2) if f E L2(Q) and ~b E Ho2(fl). Since u h = curl ~bh is known, rather than the exact velocity u solving (4.4) differs from the usual case in that the data functional (4.5) must be discretized. Then the discrete problem is to find ph E S h C L02(~2) such that _ f a ph div vh = gh(Vh)
V vh @ Vh C H~(f~)
(4.6)
M.D. GUNZBURGER el al.
372
where if g = g ( f ~k), gh = g ( f ~bh). V h must be chosen so that the div-stability condition holds [17]. For the Clough-Tocher element, it is shown in [27] that choosing v h = {vh ~ H~(~2)lvhle, is quadratic; vh ~ C((~)} where ei are the microelements of the Clough-Tocher triangulation, and choosing S h = {qh G_ L2(f~)lqhlE, is linear}
where Ei are the macroelements of the triangulation gives a ,oh satisfying lip
-
PhlIL I )
QIIp-/~IIL2<~>+
c211
-
hl12
where C = C(f, v), and 13h is an arbitrary element of S h, so that for the elements being used, lip - Phil ~ C([Ip[12 + 1l~&ll4)h z
assuming the necessary regularity. This is the optimal result in the sense that using lower order pressure elements would decrease the power of h, whereas using higher order elements would not increase it. It is also optimal in the approximation theoretic sense. Algebraically, we usually do not want to impose explicitly the condition
foe--0 Instead, no normalizing condition on ph is enforced a priori. Then we get a matrix problem from (4.6), BXp = G,
(4.7)
where B x is a matrix which may be thought of as a discrete gradient operator. Although B T is rectangular in general (4.7) will have a unique solution, modulo an additive constant, as follows from the div-stability of (4.6). To compute p it is simplest to solve the positive semidefinite equation = 86
where the nullity of B B x is 1. Similar considerations can be applied to pressure recovery from other Navier-Stokes formulations not explicitly involving pressure, for example streamfunction-vorticity cases. Acknowledgments--This work was supported in part by the National Aeronautics and Space Administration under contract number NAS-I-17070 while the authors were in residence at the Institute for Computer Applications in Science and Engineering. REFERENCES I. F. Thomasset, Implementation of Finite Element Methods for Navier-Stokes Equations. Springer ( 1981 ). 2. P. Jamet and P. Raviart, Numerical solution of the stationary Navier-Stokes equations by finite element methods, in Proceedings of lnternational Symposiums Organized by IRIA-LABORIA I, 192-223. Springer (1973). 3. V. Girault and P. Raviart, Finite Element Approximations of the Navier-Stokes Equation. Springer (1979). 4. M. Gunzburger and J. Peterson, On the finite element approximation of the streamfunction-vorticity equations, in Advances in Computer Methods for Partial Differential Equations, V 47-56 IMACS (1984). 5. M. Gunzburger and J. Peterson, Treatment of boundary conditions and multiply connected domains in the finite element approximation of the streamfunction-vorticity equations, to appear. 6. A. Campion-Renson and M. Crochet, On the streamfunction-vorticity finite element solutions of NavierStokes equations. Int. J. Num. Meth. Eng. 12, 1809-1818 (1978). 7. K. Barrett, A variational principle for the strcamfunction-vorticity formulation of the Navier-Stokes equations incorporating no-slip conditions. J. Comp. Phys. 26, 153-161 ( ! 978). 8. G. Dhatt, B. Forno, and C. Bourque, A ~b-to finite element formulation for the Navier-Stokes equations. Int. J. Num. Meth. Eng. 17, 199-212 (1981).
Computation of incompressible viscous flows
373
9. Philip M. Gresho, S. T. Chan, C. D. Upson and R. L. Lee, A modified finite element method for solving the time-dependent, incompressible Navier-Stokes equations. Int. J. Num. Meth. F1. (1983). 10. H. Sehlichting, Boundary Layer Theory, 7th Ed. McGraw-Hill (1984). I I. J. Sen-in, Mathematical principles of classical fluid mechanics, in Handbuch der Physik, VII/l, 125-263. Springer (1959). 12. P. Roache, Computational Fluid Dynamics. Hermosa (1972). 13. G. Fix and M. Gunzburger, Downstream boundary conditions for viscous flow problems, Comput. Math. Appl. 3, 53-63 (1977). 14. S. G. Rubin, A review of marching procedures for parabolized Navier-Stokes equations, Proc. Syrup. on Num. Phys. Aspects of Aero. Flows. Springer-Verlag, New York ( 1981). 15. R. Chilukuri and R. H. Pletcher, Numerical solutions to the partially parabolized Navier-Stokes equations for developing flow in a channel. Numer. Heat Transfer 3, 169-188 (1980). 16. M. Crouzeix and P. Raviart, Conforming and non-conforming finite element methods for the stationary Stokes equation. RAIRO R-3, 33-76 (1973). 17. J. M. Boland and R. Nicolaides, Stability of finite elements under divergence constraints. SIAM J. Numer. Anal. 20, 722-731 (1983). 18. J. M. Boland and R. A. Nicolaides, Stable and semi-stable low order finite elements for viscous flows. SIAM ,L Num. An. 22, 474 (1985). 19. M. Gunzburger and J. Peterson, On conforming finite element methods for the inhomogeneous stationary Navier-Stokes equations. Numer. Math. 42, 173-194 (1983). 20. C. Taylor and P. Hood, A numerical solution of the Navier-Stokes equations using the finite element method. Compt. Fluids, l, 73-100 (1973). 2 I. M. Gunzburger and R. Nicolaides, Some aspects of finite element approximations of incompressible viscous flows, to appear in Recent Advances in Numerical Methods in Fluids, III. 22. F. Brezzi and P. Raviart, Mixed finite element methods for 4th order elliptic problems, in Topics in Numerical Analysis III. Academic (1977). 23. R. Falk and J. Osborn, "Error estimates for mixed methods," RAIRO, 249-261 (1980). 24. G. Fix, M. Gunzburger, R. Nicolaides, and J. Peterson, Mixed finite element approximations for the biharmonic equation, in Proceedings of the 5th International Symposium on Finite Elements and Flow Problems, 281-285. Univ. of Texas (1984). 25. R. Scholz, Approximation von Sattlepunkten mit finite elements. Bonnet Math. Schriften, 58-67 (1976). 26. R. Scholz, A mixed method for 4th order problems using linear finite elements. RAIRO R-3, 85-90 0978). 27. E. M. Cayco and R. A. Nicolaides, Recovery of pressure field from finite element streamfunction NavierStokes calculations, submitted to Num. Math. 28. M. D. OIsen and S. Y. Tuann, New finite element results for the square cavity. Comps. and Fluids 7, 123135 (1978). 29. H. Blum and R. Rannacher, On boundary value problems of the hiharmonic operator on domains with angular comers. Math. Meth. Appl. Sci. 2, 556-581 (1980). 30. P. Ciarlet, The Finite Element Method for Elliptic Problems. North-Holland (1978).
0045-7930/85 $3.00+ .00 © 1985PergamonPressLtd.
Printedin the U.S.A.
ALGORITHMIC
AND
COMPUTATION FLOWS
BY
OF
THEORETICAL
RESULTS
INCOMPRESSIBLE
FINITE
ELEMENT
ON
VISCOUS
METHODS
M. D. GUNZBURGER and R. A. NICOLAIDES Department of Mathematics, Carnegie-Mellon University, Pittsburgh, PA 15213, U.S.A. and C. H. LIU Low Speed AerodynamicsDivision, NASA-LangleyResearch Center, Hampton, VA 23665, U.S.A. (Received 14 December 1984; in revised form 12 February 1985)
Abstract--Primitive variable as well as streamfunction-vorticity and pure streamfunction formulations are discussed. For the primitive variable case alternative choices of the viscous stress term are shown to produce natural boundary conditions which are well suited for matching to various far field conditions. For the other cases recent analytical results, including error estimates are described, and an optical algorithm for pressure recovery as well as treatment for multiply connected domains are given. I. INTRODUCTION Many different weak formulations are used as a basis for finite element solutions of the incompressible Navier-Stokes equations, ranging from fourth order streamfunction techniques through second order streamfunction-vorticity methods, velocity-vorticity-pressure systems and, most important perhaps, the standard primitive variable approach. For a survey, see [1]. Of these methods, the primitive variable technique has been most thoroughly analyzed [2, 3] although the streamfunction-vorticity method (finite element style) has also received attention [3-8]. Most of the analysis deals with interior flows driven by body forces. In applications, external flows are also important, and in particular the treatment of far field boundary conditions is often a major difficulty. In such cases, the natural boundary conditions generated by finite element schemes often have a valuable role to play (see e.g. [9]). Even for a given system of variables the natural boundary conditions can vary according to how the viscous and convection terms are defined. In this paper we shall introduce some new results for the streamfunction and streamfunction-vorticity formulations, including error estimates, optimal recovery of the pressure field, the question of multiply connected regions, and also possible choices for natural boundary conditions. For the primitive variable case three types of natural boundary conditions will be considered arising from different formulations of the viscous terms in the equations. Additionally, the effect on the natural boundary conditions of alternate treatments of the convection terms is considered, along with the manner in which a parabolized Navier-Stokes system interacts with the natural boundary conditions in general. 2. PRIMITIVE VARIABLE FORMULATIONS 2.1. For stationary incompressible flows the governing equations are u. Vu + XTp = div a,
(2.1)
where u is the velocity field, p the pressure and a is the viscous stress tensor, defined by o- = P(VI! + VuT). 361
(2.2)
362
M.D. GUNZBURGERet al.
The constant density is incorporated into both p and a and of course we also have div u = 0.
(2.3)
In view of (2.3), whenever v is constant, div a = div (v(Vu + VuT))
(2.4)
= .(Au)
(2.5)
= - v curl curl u.
(2.6)
Although in the strong sense these are equivalent, in the weak sense (2.4)-(2.6) each yield distinct natural boundary conditions which we shall discuss. Let ~ denote a bounded domain of R 2, not necessarily simply connected, with boundary I'. Standard Galerkin technique yields a weak problem of the form
a(u, v) + b(p, v) + c(u, u, v) = (f, v)
b(q, u)
(2.7) (2.8)
= O,
where (u, p) denotes the sought velocity-pressure variables and (v, q) are test functions. Below it will be necessary to consider several choices for spaces containing these functions. In all cases
(f, v) = fn f" v, and
b(q, v)
=
-fq
div v.
(2.9)
However, (2.4)-(2.6) give, respectively,
ff.i(Vu+VAT). Vv a(u, v) =
v
(2.10)
~Tu. Vv
(2.1 1)
v f curl u curl v. dn
(2.12)
(We note that henceforth, since we are dealing with domains in R 2, curl u will denote the single non-trivial component of the curl of u, i.e. curl u is a scalar.) We note at the outset that although (2.12) has some useful properties associated with it, it suffers from the serious drawback of requiring the use of divergence free elements. There is also some flexibility in the choice of c(u, u, v). The first possibility is the obvious one, c(u, u, v) = f~ (u- ~Tu). v.
(2.13)
The second is the skew-symmetrized form [2] c(u, u, v) = ~
(u. Vu). v - (u. Vv). u
.
(2.14)
We note that (2.14) can be used, and indeed is equivalent to (2.13), only when div
Computation of incompressible viscous flows
363
u = 0 and either u. n = 0 or u is specified on t2. Otherwise a boundary integral must be added to (2.7). See (2.18) below. Associated with these weak forms are certain essential and natural boundary conditions. These can be deduced from the boundary integrals omitted during the formation of (2.9)(2.12) and (2.14). Specifically we have respectively r v. [pn - v(Vu + VuT) • n]
r v.
[pn -- vVu. n]
f r [pv. n - v curl u(v. r)],
(2.15)
(2.16)
(2.17)
where n and ~- respectively denote the unit outer normal and unit tangent vectors to F. If (2.14) is used, to each of these must be added
1f [ (u. n)(u. v).
2 Jr
(2.18)
For every case (2.15)-(2.18) specification of u . n and/or u . 7 , i.e. u, implies that v . n l r = 0 and/or v. fir = 0. For example, if Ulr is given, then Vlr = 0 and all the boundary integrals (2.15)-(2.18) vanish. In this case, there is no obvious particular advantage of one weak form over any other. On the other hand, for other boundary conditions some of the natural conditions can be used to advantage. Unfortunately (2.1 l) which is the most widely used formulation, has the disadvantage that the associated natural condition, which from (2.16) is t a n t a m o u n t to specifying - p n + vVu. n appear to have no physical significance, and in particular is not the stress vector at the boundary. To have specification of the stress vector - n p + v(Vu + VuT) • n be a natural condition, formulation (2.10) with (2.13) should be used if the stress vanishes. Otherwise (2.15) appears as data. Use of (2.10), (2.13) and (2.15) also allows, for example, the use of the normal velocity as an essential condition paired with the tangential stress as natural. Physically, being able to specify stress components is often necessary at interfaces between different fluids and in some free boundary calculations. Another possibility is to use (2.12), (2.13) and (2.17). Here, the natural conditions involve p and ~" -- curl u. For example, both m a y be specified as natural conditions or we m a y pair p with the tangential velocity or ~" with the normal velocity. Such combinations are useful for blending fully viscous flows with both inviscid and boundary layer type flows. This is discussed in Section 2.2. Following the remarks above we see that (2.14) can best be used when u. nit = 0 or u specified on all boundaries so that V[r = 0. Otherwise it is necessary to use (2.13) since no physical meaning seems to attach to (2.18). 2.2. In this section we shall discuss in more detail the use of (2.12), (2.13) and (2.17) to simulate far field conditions at a computational boundary. Specifically three cases that are practically important are inflow, outflow and matching to given inviscid flows. Consider first the inflow conditions. In some cases, e.g. channel flows, both components of u m a y be accurately known and so they can be specified as essential boundary conditions. However, in other cases, for example with a boundary layer inflow, the
364
M.D. GUNZBURGERel al.
boundary layer equations correctly determine u, whereas even in the simplest case of a flat plate the Blasius solution for v does not match the inviscid flow as y --~ ~ [10]. Thus, some method for imposing the inflow condition which does not involve knowledge of v would be preferable. It is natural to impose u as an essential condition, and, in view both of the inviscid flow at great distances from the boundary and our previous discussion of natural boundary conditions, to take as the other specified variable to be the vorticity ~. Then both u and ~" can be made arbitrarily close to their inviscid flow values by making the computational boundary sufficiently distant. To impose the vorticity condition, since ~"is not zero but is computed from the boundary layer equations, e.g. the Blasius solution, it is necessary to incorporate the boundary integral (2.17) into the weak form of the equations after substituting v. n = 0 and the known function curl u = ~'. This is illustrated in a specific case in subsection 2.3. Turning next to computational boundary segments which are nominally in inviscid flow regions, typically the known inviscid solution is imposed. One method is simply to specify u. An alternative is to specify p and ~" or some combination of these two possibilities. Usually the p, ~"choice is the preferred one because these variables approach their asymptotic limit faster than does u [ 11 ]. Once again, (2.17) with the known values p and ~" substituted must be restored to the weak form (2.12). Again, see Section 2.3 below. Next we consider the notoriously difficult question of finding satisfactory outflow conditions. Ideally these should have no observable upstream influence. Were the situation one of simple convection-diffusion, satisfactory conditions can be fairly easily given in most cases [12, 13]. However the fact of incompressibility with its attendant global coupling of the pressure field renders the Navier-Stokes case rather more difficult. Often, for example for wakes, boundary layer or channel flows, the outflow has the property that the diffusion in the streamwise direction is negligible compared with both the convection and the cross-stream diffusion. Mathematically, if the x-axis is aligned with the flow as we shall assume in the sequel, then viscous terms involving x-differentiations may be neglected. This observation is systematically exploited in the parabolized NavierStokes equations (PNS) as derived, e.g., by [14, 15]. Instead of imposing outflow conditions for the full Navier-Stokes equations, we shall assume that in the vicinity of the outflow boundary, the flowfield is described by the PNS equations. Thus the boundary conditions for PNS equations are what determine the outflow conditions. For example, for the PNS equations (2.12), (2.17) are replaced by
u fe (curl u)vl,y
v = (Vl, v2)
(2.19)
and
frP
V. n.
(2.20)
Thus, only p at the outflow must be specified. No vorticity or velocity is required. The need to impose the pressure arises from its previously mentioned global character. 2.3.
Here we give a specific formulation of a problem illustrating the above considerations. We consider boundary layer flows over obstacles of size comparable with the boundary layer thickness. A simple example of such a flow region is shown in Fig. 1. Here AOBC represents a solid boundary at which u = 0. AE is an inflow boundary at which a boundary layer profile is imposed. We assume ED is sufficiently far removed from the solid boundary so that it may be regarded as being in the inviscid regime. CD represents the outflow boundary. We shall assume that the inflow on AE is given by the Blasius profile, with a leading edge at x = Xo < xA. Following the discussion of the previous section, u and ~"will be prescribed (i.e. [10])
Computation of incompressibleviscousflows
365
D
E
B
C
t,y 0
A
J-
x
Fig. 1. Computationaldomain for flowover a step.
on AE
u~. = U ~ f ' ( n )
and fi~ = U o ~ f " ( n ) [ U o o / v ( x A -
Xo)]1/2
on AE,
where ~ = y [ U o o / v ( x A -- x0)] ~/2 and f is the Blasius function. Since on ED the flow is effectively inviscid, we may use the potential flow solution to get a pressure boundary condition and set ~"= 0. In general, the potential solution for p must be found numerically (e.g. by a panel method), although in the step case being considered a closed form solution can be found from the Schwartz-Christoffel transformation. In fact, on ED, we have for a step of height d, i.e. YB = d, 1 [ ( ~ + l ) ( f f + 1)11/2 Ppot = P~o + ~ ~'2UL -- 1)(w-- ~ J ~'~o, = O,
where w is obtained from z = x + i y e by
Z =
d
--
[(W 2
_ 1)l/2 - - cosh -I w ] + i d
and where we recall that we have absorbed the constant density into the pressure variable. We note that lim ~'in~
~'pot =
0
for
YE < ~ "
Y~YE
However, the discrepancy can be made arbitrarily small by making Y e sufficiently large. In practice, f ~ decays sufficiently rapidly for this not to be a practical issue. For example i f R e = Uood/v = 1000, then with YE = 4yB = 4d we have ~'i,(YE) "" l O - S U o ~ / d . Consider next the outflow DC. Through the use of PNS only p needs to be specified. Unfortunately in most cases p is either unknown or difficult to find. Two possible choices are (a) to use experimental data or (b) to assume that the flow at DC is of boundary layer type so that the pressure is constant on DC and is given by pou,(y) =
p~,(x~, yz).
To conclude, we shall now give the appropriate weak form incorporating this particular set of boundary conditions, and the PNS equations in the vicinity of the outflow. Let
366
M.D. GUNZBURGERet al. H~([2) : {u E
n'(f01UlAo~c -- 0,
UIAE
-----
Uin}
H6(f~) = {u ~ H'(f~)luLAosc = 0, UIAE = 0}, where H1(t2) is the usual Sobolev space. T h e n seek u ~ Hel(~2) and p @ L2(f/), f~ v curl U(l)l,y -- O/V2,x) q- fi2 (U'VU)" V --
~t p d i v v
-- ~AE V~'in(V'T)-- fED Pp°tV'n--
fnq h div u h =
fcDP°mV'n
for all
0
qh E
for all v E H6(~2),
L2(f~).
(2.21) (2.22)
In (2.21) a = 1 in the main flow field, and a = 0 near the outflow boundary with, if desired, a smooth transition. This simply activates the PNS approximation at appropriate field points. N o rigorous analysis of (2.21)-(2.22) is available. However in practice element pairs suitable for wall bounded flows in the standard formulation appear to be satisfactory in the above described case. See, e.g. [1-3, 16-21], for examples of such element pairs and for a discussion of error estimates.
3. S T R E A M F U N C T I O N - V O R T I C I T Y
FORMULATIONS
3.1. The stationary Navier-Stokes equations, written in terms of the streamfunction ~b and the vorticity ~', are given by Aft = --~'in ~2
(3.1)
and
where in (3. l) and (3.2) v is the kinematic viscosity, f the body force, and 12 is a bounded region in R 2. If ~2 is multiply connected, we denote by Fo its exterior boundary and by Fi, i = 1. . . . . m, the remaining parts of the boundary. (See Fig. 2.) m
Suppose that at these boundaries the velocity u is specified, i.e. u = g on F = U I?i.
i=0
I"
0
Fig. 2. A multiply connected domain ft.
367
Computation of incompressible viscous flows In order for a streamfunction to exist, we must have that [3] for
fr g.n=0
i=0 .....
(3.3)
m,
i
i.e. there is no net mass flow through any o f the b o u n d a r y pieces. It is well k n o w n that the streamfunction is unique up to an additive constant, and we fix its value by specifying it to be zero at an arbitrary point x0 on £o. Now, let q denote a function such that
cgq Or
g . n on £
q(xo) O.
and
(3.4)
T h e n the b o u n d a r y conditions for the system (3. l) and (3.2) are given by
~=q+aionri,
=qon£o,
i = 1. . . . .
m
(3.5)
and &b -r = -g'r
on £,
On
(3.6)
where ai, i = l . . . . . m are constants to be determined as part o f the solution. These constants are fixed by the requirement that the pressure be single valued, i.e. the change in the pressure p a r o u n d any closed curve surrounding any o f l~i, i = l . . . . . m, is zero. Since the m o m e n t u m equation stipulates that 1 -
Vp = - u . Vu + vAu + f,
(3.7)
P where a is the constant density, we have that, for i = l . . . . .
m,
o = ~ f r V( P~ 'Aru=-fur' V u + f ) ' r ' ~ Since A u - z = 0f/0n whenever div u = 0 and u . Vu. r = f u - n + O[(u. u)/2]/Oz, we then have that
fr
(V Of/an - f g . n + f- r) = 0
for
i = 1. . . . .
m.
(3.8)
i
T h u s (3. l), (3.2), (3.5), (3.6) and (3.8) are the governing equations and side conditions which are to determine the functions ~k and f and the constants ai, i -- 1. . . . . m. As in the primitive variable case, other b o u n d a r y conditions are o f interest. For example, if all or part o f F0 is a c o m p u t a t i o n a l b o u n d a r y in the inviscid potential flow regime one m a y want to specify f = Of/On = 0 there, or if part o f Fo represents a b o u n d a r y layer inflow, one m a y want to specify ~b and ~. Here we will only consider the b o u n d a r y conditions (3.5), (3.6) and (3.8). Other b o u n d a r y conditions m a y also be i m p l e m e n t e d in a straightforward m a n n e r (see [5]). The weak formulation o f the streamfunction vorticity equation requires the use o f some n o n s t a n d a r d function spaces. For details concerning these spaces see [3-5]. Here we will proceed formally. The most direct weak formulation is to seek ~b E Se, f E Z such that
fn fco-fv~b. Vo~=frO~g.r -v
Vf . V~P+
f
0tP
Ox o~b Ox
=-
for all ~o E Z f-curltP
forall~P~S0.
(3.9) (3.10)
368
M.D. GUNZBURGER el aL
Here Z = Hl(fl), Se is a subset of HI(f/) whose elements satisfy (3.5), and So is a subspace of Hl(~2) whose elements satisfy (3.5) with q = 0. In deriving (3.9) we have replaced the boundary integrals
o~ f r co ~nn
(3.11)
with the boundary integral on the right hand side of (3.9). Likewise, in deriving (3.10), the boundary integral 0~"
0ff
r~
(3.1 2)
!
appears. However, note that ~P = 0 on F0 and ~ = constant on I?i so that, with (3.8), (3.12) vanishes. (For details see [3-5].) We note that no boundary conditions on ~"need be contrived on solid walls. Also, the conditions (3.8) are natural to the weak formulation, as is (3.6) while (3.5) is an essential condition. One may view (3.9), (3.10) as a mixed method for the fourth order streamfunction equation. Other mixed methods for the linear case are discussed in, e.g. [3, 22, 23], and for the nonlinear case in [3]. The particular method given by (3.9), (3.10) is discussed, in the linear case, in [23-26] and for the nonlinear case in [3-5]. With regard to finite element discretizations of (3.9), (3.10), one is particularly interested in using merely continuous finite element spaces for both ~b and ~. A sample error estimate is given in the following result [24]. THEOREM 3. l Suppose piecewise continuous linear finite element spaces are used to approximate both ~k and L Then there exists a constant C independent of the grid size h such that I[Vff - XTffhl[0< Ch
K' II0 ~
and
Ch ~n
Here ~kh and ~-h denote the discrete solution. These estimates display the general trend, i.e. the gradient of streamfunction, and therefore the velocity field, is optimally approximated, while there is a loss of accuracy in the vorticity approximation. We remark that if one uses continuously differentiable finite element spaces for ~k and ~, then both ~k and ~"are optimally approximated [24]. 3.2.
The main difficulty which multiply connected domains present is that the streamfunction may be arbitrarily specified on only one part of the boundary. Thus, one sets and ~bh = 0 at one point on Fo, but these are determined only up to unknown constants ai and a~, respectively, on the remaining parts Fi, i = 1. . . . . m, of the boundary. In the discretization of (3.9), (3.10), the constants a h are determined as part of the solution process. However, that method requires the use of N+M,
E
%(x),
i = 1. . . . .
m
(3.13)
j=N+Mi-j+I
as basis functions associated with the a~, i = 1. . . . . m. Here tpj, j = 1. . . . . N represents basis functions associated with interior unknowns, and ~Pj, j = N + Mi-i + 1. . . . . N + Mi, represents basis function associated with unknowns on Fi, i = 1. . . . . m, with M0 = 0. In general, these basis functions are only "semi-local" in the sense that they couple all the points on a boundary F~. In particular, in the discrete set of equations, a h is coupled with all the unknowns associated with elements whose closure intersects with Fi. The result of this is that the bandwidth of the linear systems encountered will be greatly increased, resulting in both a computer storage and time penalty over simply connected domain problems. Some of this penalty may be mitigated by employing a numbering
Computation of incompressible viscous flows
369
scheme resulting in banded-bordered matrices wherein the onerous couplings are found near the bottom and the right of the matrices encountered. This, of course, requires the development of special programs to solve the linear systems. An alternative to using basis functions such as (3.13) is to guess the value of the constants ai, i = 1. . . . . m, appearing in (3.5). Then one may solve for ~ and ~" from (3.1), (3.2), (3.5) and (3.6). However, in general, (3.8) will not be satisfied. At this point one m a y change the guesses for the a~, repeating the process until convergence is achieved, i.e. until (3.8) is satisfied. In more detail, we proceed as follows. Given guesses a~, i = 1, . . . . m, k E Z +, we define ~kh E Sk and ~-h C Z h [where
ZhCZ and
s~c& and where now the constants ai = a/*, i.e. are assumed known] to be the solution of the discrete forms of (3.9), (3.10) where now (3.10) holds for all ¢ph E S0h = S h fq H~(~2). The discrete equations and unknowns may be ordered in such a way that the bandwidth of the discrete system is no larger than that obtainable in an analogous discretization of a problem posed in the simply connected domain bounded by Po- Having obtained ~-h from a~, i = 1. . . . . m, we m a y substitute into (3.8). In practice, we would evaluate the integrals appearing in (3.8) by a numerical quadrature formula. Let us denote such a n approximation to each equation in (3.8) by Ii(ak), i = 1. . . . . m, which in general will not be small. O f course, each Ii(. ) is a nonlinear function of a k = (a~ . . . . . akin) through the discrete vorticity ~-h. We now need an algorithm to generate new g u e s s e s a k+l, from which ffkh+l and ~'h+1 m a y be computed and for which Ii(a k+l) is closer to zero than was Ii(ag). The most practical way to accomplish this is to use a secant type method. For example, if m = 1, i.e. f / i s doubly connected, we m a y then define a k+' = a~ -
a~ - a~ ' ii(alk). Ii(a~) -- Ii(a~ -l)
For m > be used. as (3.14) practical
(3.14)
1, some generalized secant method in R m such as the Wolfe secant method can (See [4, 5] for details about such methods.) Note that the use of formulas such requires two starting guesses a ° and a ~. We also note that in general it is not to use Newton's method to update a k since evaluating the Jacobian of I(a) = ( I 1 , - . - , Ira) requires the solution of linear problems of the same size as those for qjh and ~-h for each pair O~lOai, Od//Oai. Thus, at the price of solving a sequence of simpler problems, one may avoid the bandwidth problems engendered by the use of basis functions such as (3.13). Whether or not the alternate method is more efficient depends, of course, on how m a n y of the simpler problems must be solved as well as on how m u c h cheaper it is to solve the simpler problems. Both of these factors, for a particular geometry, will depend on the n u m b e r of degrees of freedom, i.e. the gridsize, and the Reynolds number. Therefore it is clear that the performance of iterations such as (3.13) play a central role in the overall efficiency of the alternate method of treating multiply connected domains. Unfortunately, the sequence a k generated by secant type methods are not guaranteed to converge for arbitrary values of the initial guesses, especially at high values of the Reynolds n u m b e r (for details see [4, 5]). Similar approaches for treating multiply connected domain may be used for the pure streamfunction formulation. 4. STREAMFUNCTION FORMULATION 4.1.
The differential equation is the well-known
M.D. GUNZBURGERet al.
370
and natural boundary conditions are obtained after two formal integrations with a test function ×; thus,
o
o×
og,
,, fo ,x,ax +,, fr ~, (a~,)×-,, fr (a~,) T + fo ,X~,(xx~,y-xy~,~)- f r x ~ ,X~,= f o b , O 0 when ~nn and ~rr represent outer normal and tangential directional derivatives. For an interior problem X and X, will be zero on F, the boundary of the domain 12, and so each boundary integral vanishes. We are assuming in this section that u, the velocity vector, in specified on all boundaries. For multiply connected domains, ~b is specified on only one bounding surface and is determined up to an arbitrary constant on the others, which of course must be found as discussed for the streamfunction-vorticity formulation above. The natural conditions are thus
= -~
=0
and ('0__~+ (u. n)~" = 0.
On Clearly then, these give on F
at" i-=~=o. These are useful for matching to an irrotational external flow when computing flows for which the flow domain includes both rotational and irrotational zones. For an example of such a case see Section 2. Of course it is also possible to make combinations of natural and essential conditions, e.g. we may specify the normal velocity as an essential condition and the vorticity as natural. In any case the weak problem is the following: find ~b @ H~(f~) such that
u f, A~Ax+Ax(xxg/y-Xyd/x) = fafx+ N(g, X)
VXE
H2e,o(12),
(4.1)
where N(g, ×) represents whatever boundary integrals are kept, with data represented by g, and He2(~2) is a subset of H2(~) satisfying the imposed essential boundary conditions and elements of He2,o(ft) satisfy homogeneous essential conditions. If g = 0, i.e., the data for the natural boundary conditions is zero, then N ( . , . ) is identically 0. To approximate ~b we proceed in the standard way to find ~bh ~ S* C H~2(f~), where now fl is a bounded domain of R 2. Equation (4. l) is replaced by the following find ~bh E S h such that v x h ~ Soh where Soh is the analogous subspace of H2,o(12). The basic error estimate for this procedure is given in the following theorem. THEOREM 4.1 A constant C exists, independent of Seh, such that for sufficiently large u the following estimate is valid 11¢ - ¢h112 ~ C inf I1~ - ~hll=.
¢*cs,*
For the proof of this result and details of the constants see [27].
Computation of incompressible viscous flows
371
Another weak formulation for the streamfunction formulation is given in [28]. The present formulation appears more direct however. Theorem 2.1 applies, of course, to conforming finite elements. Here, the element choices are somewhat limited by the requirement of C 1 inter-element continuity. The lowest degree conforming finite element on triangles is the Clough-Tocher cubic macroelement. A possible alternative is the Argyris 21 degree of freedom quintic elements, free of course, of the need for special macro constructions. However it is desirable to use the lowest order of polynomials possible for several reasons. First, apart from exceptional circumstances, the exact solution is not in a sufficiently high order Sobolev space for the full approximation capabilities of the quintics to be realized. Second, the costs of the numerical integration for the assembly process are determined by the need to integrate 6th degree polynomials exactly for the quintic case but only 2nd degree polynomials for the cubic macroelement case. For these reasons we recommend the Clough-Tocher element. In this case we obtain from Theorem 2. l the following error estimate: I1¢ - ¢h112 ~ C h 2 - S l l ¢ l h - s
s = 0, 1.
(4.2)
In general, however, we must expect to have to take s > 0 as shown by results of [29]. However for a square domain with f E L2([2) we can take s = 0, although even in this relatively favorable case we cannot improve on (4.2) with the quintic element. Equation (4.2) used in conjunction with the standard duality method gives estimates in lower norms, for example, we have the estimate (4.3) for a square domain. This of course gives an O(h 3) rate for [[u - uh[[o in the most favorable case, where u and u h are the velocity fields defined from ~b and ~bh. In practice, numerical integrations are used to evaluate all contributions to stiffness matrices and right hand side terms. One valuable observation is that one needs to evaluate only the viscous, i.e. biharmonic, terms exactly, in order to preserve the estimates (4.2), (4.3), provided that the same quadrature rule is used for all the other terms. For example, for the Clough-Tocher element, it is sufficient to use a rule exact for quadratic polynomials, even though the convection term involves 5 th degree polynomials. This has been verified by the authors in practical computations and anyway is predicted at least for linear problems by finite element theory [30]. 4.2.
A pressure field, of optimal accuracy for the given velocity approximation may be computed (without specifying any explicit boundary conditions) from the primitive variable m o m e n t u m equation, by discretization of the equation
-re p div v = g(v)
V v E V,
(4.4)
where g(. ) denotes the linear functional
g(v)=f(v)-ufVu. Vv-f(u.Vu)v
u = curl ~b
(4.5)
which is bounded on H~([2) if f E L2(Q) and ~b E Ho2(fl). Since u h = curl ~bh is known, rather than the exact velocity u solving (4.4) differs from the usual case in that the data functional (4.5) must be discretized. Then the discrete problem is to find ph E S h C L02(~2) such that _ f a ph div vh = gh(Vh)
V vh @ Vh C H~(f~)
(4.6)
M.D. GUNZBURGER el al.
372
where if g = g ( f ~k), gh = g ( f ~bh). V h must be chosen so that the div-stability condition holds [17]. For the Clough-Tocher element, it is shown in [27] that choosing v h = {vh ~ H~(~2)lvhle, is quadratic; vh ~ C((~)} where ei are the microelements of the Clough-Tocher triangulation, and choosing S h = {qh G_ L2(f~)lqhlE, is linear}
where Ei are the macroelements of the triangulation gives a ,oh satisfying lip
-
PhlIL I )
QIIp-/~IIL2<~>+
c211
-
hl12
where C = C(f, v), and 13h is an arbitrary element of S h, so that for the elements being used, lip - Phil ~ C([Ip[12 + 1l~&ll4)h z
assuming the necessary regularity. This is the optimal result in the sense that using lower order pressure elements would decrease the power of h, whereas using higher order elements would not increase it. It is also optimal in the approximation theoretic sense. Algebraically, we usually do not want to impose explicitly the condition
foe--0 Instead, no normalizing condition on ph is enforced a priori. Then we get a matrix problem from (4.6), BXp = G,
(4.7)
where B x is a matrix which may be thought of as a discrete gradient operator. Although B T is rectangular in general (4.7) will have a unique solution, modulo an additive constant, as follows from the div-stability of (4.6). To compute p it is simplest to solve the positive semidefinite equation = 86
where the nullity of B B x is 1. Similar considerations can be applied to pressure recovery from other Navier-Stokes formulations not explicitly involving pressure, for example streamfunction-vorticity cases. Acknowledgments--This work was supported in part by the National Aeronautics and Space Administration under contract number NAS-I-17070 while the authors were in residence at the Institute for Computer Applications in Science and Engineering. REFERENCES I. F. Thomasset, Implementation of Finite Element Methods for Navier-Stokes Equations. Springer ( 1981 ). 2. P. Jamet and P. Raviart, Numerical solution of the stationary Navier-Stokes equations by finite element methods, in Proceedings of lnternational Symposiums Organized by IRIA-LABORIA I, 192-223. Springer (1973). 3. V. Girault and P. Raviart, Finite Element Approximations of the Navier-Stokes Equation. Springer (1979). 4. M. Gunzburger and J. Peterson, On the finite element approximation of the streamfunction-vorticity equations, in Advances in Computer Methods for Partial Differential Equations, V 47-56 IMACS (1984). 5. M. Gunzburger and J. Peterson, Treatment of boundary conditions and multiply connected domains in the finite element approximation of the streamfunction-vorticity equations, to appear. 6. A. Campion-Renson and M. Crochet, On the streamfunction-vorticity finite element solutions of NavierStokes equations. Int. J. Num. Meth. Eng. 12, 1809-1818 (1978). 7. K. Barrett, A variational principle for the strcamfunction-vorticity formulation of the Navier-Stokes equations incorporating no-slip conditions. J. Comp. Phys. 26, 153-161 ( ! 978). 8. G. Dhatt, B. Forno, and C. Bourque, A ~b-to finite element formulation for the Navier-Stokes equations. Int. J. Num. Meth. Eng. 17, 199-212 (1981).
Computation of incompressible viscous flows
373
9. Philip M. Gresho, S. T. Chan, C. D. Upson and R. L. Lee, A modified finite element method for solving the time-dependent, incompressible Navier-Stokes equations. Int. J. Num. Meth. F1. (1983). 10. H. Sehlichting, Boundary Layer Theory, 7th Ed. McGraw-Hill (1984). I I. J. Sen-in, Mathematical principles of classical fluid mechanics, in Handbuch der Physik, VII/l, 125-263. Springer (1959). 12. P. Roache, Computational Fluid Dynamics. Hermosa (1972). 13. G. Fix and M. Gunzburger, Downstream boundary conditions for viscous flow problems, Comput. Math. Appl. 3, 53-63 (1977). 14. S. G. Rubin, A review of marching procedures for parabolized Navier-Stokes equations, Proc. Syrup. on Num. Phys. Aspects of Aero. Flows. Springer-Verlag, New York ( 1981). 15. R. Chilukuri and R. H. Pletcher, Numerical solutions to the partially parabolized Navier-Stokes equations for developing flow in a channel. Numer. Heat Transfer 3, 169-188 (1980). 16. M. Crouzeix and P. Raviart, Conforming and non-conforming finite element methods for the stationary Stokes equation. RAIRO R-3, 33-76 (1973). 17. J. M. Boland and R. Nicolaides, Stability of finite elements under divergence constraints. SIAM J. Numer. Anal. 20, 722-731 (1983). 18. J. M. Boland and R. A. Nicolaides, Stable and semi-stable low order finite elements for viscous flows. SIAM ,L Num. An. 22, 474 (1985). 19. M. Gunzburger and J. Peterson, On conforming finite element methods for the inhomogeneous stationary Navier-Stokes equations. Numer. Math. 42, 173-194 (1983). 20. C. Taylor and P. Hood, A numerical solution of the Navier-Stokes equations using the finite element method. Compt. Fluids, l, 73-100 (1973). 2 I. M. Gunzburger and R. Nicolaides, Some aspects of finite element approximations of incompressible viscous flows, to appear in Recent Advances in Numerical Methods in Fluids, III. 22. F. Brezzi and P. Raviart, Mixed finite element methods for 4th order elliptic problems, in Topics in Numerical Analysis III. Academic (1977). 23. R. Falk and J. Osborn, "Error estimates for mixed methods," RAIRO, 249-261 (1980). 24. G. Fix, M. Gunzburger, R. Nicolaides, and J. Peterson, Mixed finite element approximations for the biharmonic equation, in Proceedings of the 5th International Symposium on Finite Elements and Flow Problems, 281-285. Univ. of Texas (1984). 25. R. Scholz, Approximation von Sattlepunkten mit finite elements. Bonnet Math. Schriften, 58-67 (1976). 26. R. Scholz, A mixed method for 4th order problems using linear finite elements. RAIRO R-3, 85-90 0978). 27. E. M. Cayco and R. A. Nicolaides, Recovery of pressure field from finite element streamfunction NavierStokes calculations, submitted to Num. Math. 28. M. D. OIsen and S. Y. Tuann, New finite element results for the square cavity. Comps. and Fluids 7, 123135 (1978). 29. H. Blum and R. Rannacher, On boundary value problems of the hiharmonic operator on domains with angular comers. Math. Meth. Appl. Sci. 2, 556-581 (1980). 30. P. Ciarlet, The Finite Element Method for Elliptic Problems. North-Holland (1978).