- Email: [email protected]
Accuracy of least-squares methods for the Navier-Stokes equations
Computers Fluids Vol. 22, No. 4/5, pp. 549-563, 1993 Printed in Great Britain. All rights reserved
ACCURACY
0045-7930/93 $6.00+ 0.00 Copyright © 1993 Pergamon Press Ltd
OF LEAST-SQUARES METHODS NAVIER-STOKES EQUATIONS
FOR
THE
PAVEL B. BOCHEV a n d MAX D . GUNZBURGER Department of Mathematics and Interdisciplinary Center for Applied Mathematics, Virginia Polytechnic Institute and State University, Blacksburg, VA 24061-0531, U.S.A.
(Receit,ed It February 1993) Abstract--Recently there has been substantial interest in least-squares finite element methods for velocity--vorticity pressure formulations of the incompressible Navier-Stokes equations. The main cause for this interest is the fact that algorithms for the resulting discrete equations can be devised which require the solution of only symmetric, positive definite systems of algebraic equations. On the other hand, it is well-documented that methods using the vorticity as a primary variable often yield very poor approximations. Thus, here we study the accuracy of these methods through a series of computational experiments, and also comment on theoretical error estimates. It is found, despite the failure of standard methods for deriving error estimates, that computational evidence suggests that these methods are, at the least, nearly optimally accurate. Thus, in addition to the desirable matrix properties yielded by least-squares methods, one also obtains accurate approximations.
1. I N T R O D U C T I O N
The approximate solution of the Navier-Stokes equations of incompressible flow has received tremendous attention from engineers and mathematicians [e.g. 1-3]. Among the more recent developments has been the use of least-squares ideas; see, for example, Ref. [4] for a recent survey of one such approach. Also, truly least-squares methods have been developed and applied [e.g. 5-14]. Here, a finite element method based on a least-squares variational principle is examined for the approximate solution of the stationary, incompressible Navier-Stokes equations. These equations are cast into a first-order system of partial differential equations involving the velocity, vorticity and pressure as dependent variables. In 3-D one has 7 unknown scalar fields. However, the application of a least-squares principle along with, for example, a Newton linearization, results in symmetric, positive definite linear systems, at least in a neighborhood of the solution. The influence of the Reynolds number on the positive definiteness of these linear systems is felt only through the size of the neighborhood. Thus, if properly implemented continuation (with respect to the Reynolds number) methods are used, one can expect to only encounter symmetric, positive definite linear systems in the solution procedure. A further advantage of this method is that a single piecewise polynomial finite element space may be used ,for all test and trial functions, i.e. one may use equal-order interpolation with respect to a single grid for all dependent variables and test functions. A final advantage resulting from the use of a least-squares principle is that, unlike some other methods involving the vorticity, no art!)Cicial numerical boundary conditions for the vorticity need be devised. The method discussed here is similar to those in Refs [5-12]; however, there are certain differences as well, especially in the order in which the least-squares, the discretization and the linearizations steps are taken. Furthermore, the analyses found in some of these papers are incorrect, leaving open the question of the accuracy of approximations. In Section 2, we define the least-squares finite element method. In Section 3, we discuss some practical and theoretical issues connected with the method described in Section 2. Then, in Section 4, we give the results of a computational study of the accuracy of the algorithm. Finally, in Section 5 we give some concluding remarks. 549
550
PAVEr. B, BOCHEV and MAX D. GUNZBURGER
2. T H E L E A S T - S Q U A R E S F I N I T E E L E M E N T M E T H O D
2.1. The velocity-vorticity-pressure equatzons Let the bounded set f~cN3 denote the flow domain and let F denote its boundary. The dimensionless equations governing the steady incompressible flow of a viscous fluid may be written in the form divu=0 curlu-to=0 vcurlto+to xu+gradp=f
in~,
(l~
inf~
(2)
in ~'2
(3)
in fL
(4)
and divto=0
where u, p and to denote the velocity, pressure and vorticity fields, respectively, v is the inverse of the Reynolds number and f is a given body force. The first-order system of partial differential equations (1)-(4) constitutes the velocity-vorticity-pressure form of the equations of steady, incompressible flow. Note that the terminology "pressure" is a little misleading since the variable p is actually the total head, i.e. p =/~ + ½1ul2, where/~ denotes the pressure. In view of equation (2), it seems that equation (4) is redundant. However. it is crucial to the stability of the least squares algorithm to explicitly require equation (4) [e.g. 5, 10]. The system (1)-(4) should be supplemented with boundary conditions. Here we examine two such boundary conditions. The first is to impose the velocity on the boundary, i.e. u=U~
on F,
(5)
where U1 denotes a given function defined along F. Note that equation (5) implies that the normal component of the vorticity is also known, i.e., to.n=n.curlU~
on F,
(6)
where n denotes the unit outer normal to ~q. To see this, one merely needs to observe that n • curl u involves only tangential derivatives of the tangential components of u and these may be deduced from equation (5). To these one must add a condition to fix the pressure; we choose to fix the pressure at a single point x0 in ~, i.e. p(x0) = Po,
(7)
where P0 is a given number. The second combination of boundary conditions we consider is the pressure and the normal component of the velocity, i.e. p = P
on F
(8)
and u-n=U;
on F,
(9)
where P and U2 denote given functions defined along F. The boundary conditions (8) and (9) are not all that useful as boundary conditions for the Navier-Stokes equations; we consider them here for two reasons. First, a complete, rigorous analysis of least-squares finite element approximations of equations (1)-(4), (8) and (9) can be given using standard techniques; this is not the case for equations (5)-(7). Furthermore, equations (1)-(4), (8) and (9) can be shown to be related to second-order elliptic partial differential equations; we will discuss these issues in more detail below.
2.2. The least-squares principle One can use equations (1)-(4) to define the least-squares functional
~-(u, to~ P) = f ( I c u r l u - to 12+ Idiv ul2 + Idiv to 12+ Iv curl to + to x u + grad p -- t]2) d~. J.
(lO)
551
Accuracy of least-squares methods for the Navier-Stokes equations
The least-squares principle then requires the minimization of this functional over appropriate function spaces. Then standard techniques from the calculus of variations may be used to deduce that the minimizers (u, p, o) of Y necessarily satisfy
s
[(curlu-o).curlv+divudivv+(vcurlo+gradp+w n
xu-f)*(o
xv)]dQ=O,
(11)
[divodivr-(curlu-o).i+(vcurlo+gradp+wxu-f).(vcurlr+rxu)]dn=O sn (12) and
s
(vcurlo+gradp+o R
xu-f)*gradqdfi=O.
(13)
In equations (1 1)-(13) the test functions (v, q, c) are required to belong to suitable function spaces; we do not go into detail here since we are primarily interested in finite element discretizations of these equations. Clearly any solution (u, 0,~) of, say equations (l)--(7), satisfies equations (1 lH13). 2.3. The 2-D case For planar flows we have that u = (u,, ul, O)Tand u, , u2 and p are functions of x, and x2 only. Then, we have that w = (0,0,w)’ where w = &,/ax, - au, /ax,. In this case, system (l)-(4) simplifies to $+$=O I au2 ----
84
ax,
ax,
o=o
am ap
vZ+--u2w=f, 2 ax,
am
ap
-vz+-+u,w=f, I
in Q,
(14)
in R,
(15)
in R
(16)
in a,
(17)
2
ax2
where f = (f, , f2, O)T.The boundary conditions (5)--(7) reduce to just equations (5) and (7) and the boundary conditions (8) and (9) remain unchanged. The functional (10) and the necessary conditions (1 lH13) also simplify in the obvious manner for 2-D problems. 2.4. Finite element methods Starting with the weak formulation (11 )(13), a conforming finite element method can be defined in a completely standard manner. We choose a finite space Sh parametrized by h. For example, for a given positive integer r, Sh could consist of continuous (over Q) piecewise polynomials of degree
on r}, on r}
552
PAVELB. BOCHEVand MAX D. GUNZBURGER
F o r the b o u n d a r y conditions (5)-(7), the discrete problem is defined as follows: seek u ~'~ S ~, to~ ~ S ~' and p~ e S ~ such that u ~ = U h and to~ . n = W h on F, ph(x°) = P0,
f [(curl u h -- toh) • curl vh + div u h div vh + (v curl to~' + grad p~ +tohXUh--f)'(to~XVh)]dfl=0
Vv~eV~,
(18)
V~heZ~
(19)
Vqh~Qh
(20)
n[div toh div ~ - (curl u h - toh). (h + (V curl toh + grad p~, +to~×u~--f)-(vcurl(~+~hxuh)]dfl=0 and
f
(vcurlto~+gradp
h+to~xu ~-f).gradqndfl=0
are satisfied. Here UIh and W h are approximations of the data U 1 and curl U~.n, respectively.
F o r example, we could define the former pair to be boundary interpolants of the latter pair. N o t e that all o f the discrete variables i.e. q~ and the c o m p o n e n t s o f u h and to~, are a p p r o x i m a t e d by the s a m e - d e g r e e continuous piecewise p o l y n o m i a l s defined with respect to a single grid.
The discrete problem for the boundary conditions (8) and (9) is also given by equations (l 8)-(20) except that now we have that u h - n = U h and ph = ph on F and the test functions spaces in equations (18)-(20) must be suitably redefined to account for the boundary conditions vh. u = 0 and q~ = 0 on F. Again, U~ and Ph are approximations to U2 and P, respectively. 3. P R A C T I C A L
AND THEORETICAL
ISSUES
3. I. N e w t o n ' s m e t h o d
The discrete equations (18)-(20) are a nonlinear system of algebraic equations that must be solved in an iterative manner. There are many methods that one might use for such a purpose; here we only consider Newton's method. However, we do remark that if a quasi-Newton method is used, it should be chosen so that it preserves symmetry and positive definiteness of the approximate Hessian matrices. Newton's m e t h o d for the solution of equations (18)-(20) is defined as follows. Given initial guesses u (°), to (°) and p(O) for u h, toh and ph, respectively, the sequence o f Newton iterates {u(k), to(k), p(k)}k>O is generated recursively by solving, for k = 1, 2 . . . . . the system fn[(curl u(k) _ to(k)), curl vh + div u (k) div yh + (v curl to (h) + grad p(k) + to~k) × U(k ]) + to(k - I) × U~k)). (to(k- ~) × Vh) + (V curl totk ~t + gradp(k
~ + took ~) × U~k-~1__ f ) . (to(kl+ Vh)]dr2
=((vcurltoIk--l)+gradptk
I)-~-2totk
I)×utkl)).(to(k-l)×vh)d~,"~
Vvh~V~,
(21)
¥~heZ~0
(22)
a[(di v totk) div ~h _ (curl u ok)-- to(kt). ~h + (v curl e) ok)+ grad p(kl + to~k~ × utk ~t + to(k ~) × U(k)). (V curl ~h) + (v curl e) ~k) + grad p~kl + to~k) × U(k- 1) + totk t) × U(k)). (~h × u~k 1)) + (V curl ~otk ~ + g r a d p fk i~+ to(k ~t × u(k- ~) _ f ) . (~h X U~k))]d92 f [ v ( f + took ~)× u~k , ) . curl~h + ( v c u r l t o ~k ~ ) + g r a d p ~ k - L ) + 2 t o k I ) × u ( k - I ) ) . ( ~ h × u ( k - ~ ) ) ] d ~
Accuracy of least-squares methods for the Navier-Stokes equations
553
and o(v curl ta ck)+ gradp (k) + ta ck) x u(k I) ..{_O.}(k- 1) X U(k)) " g r a d q* d D
= jn(f + took- t) x utk- l)). grad
qh d ~
Vq h E Qh.
(23)
The system of linear algebraic equations (21)-(23) that is used to determine the kth Newton iterate from the (k - 1)th looks rather formidable. However, it also has some very good features. First, it is easy to see that this system is symmetric. Moreover, in a neighborhood of a minimizer, the Hessian matrix for functional (10) is necessarily positive definite; but this Hessian matrix is exactly the coefficient matrix of equations (21)-(23). Thus, in a neighborhood of a solution of equations (18)-(20), the system (21)-(23) is symmetric and positive definite. This feature is independent of the value of v, i.e. of the value of the Reynolds number. These observations, along with the guaranteed local and quadratic convergence of Newton's method, are potentially very advantageous. 3.2. Continuation methods
At this point one may well ask what goes wrong with the method as the Reynolds number increases? Surely one difficulty is that the attraction ball for Newton's method, or for any other iterative method for solving the nonlinear equations, decreases in size with increasing Reynolds number. This is known to be true for other discretizations for the Navier-Stokes equations. A related observation is that the positive definiteness of the Hessian matrix is guaranteed only in a neighborhood of a minimizer; again the size of this neighborhood surely decreases with increasing Reynolds number. As a result, for an arbitrary initial guess we may have that Newton's method does not converge and/or that the coefficient matrix in system (21)-(23) is not positive definite. The former is, of course, unacceptable, while the latter occurrence would preclude the use of simple iterative methods for solving the linear systems (21)-(23) that define the Newton iterates. In order to determine an initial guess that is within the attraction ball of Newton's method, and also is such that the coefficient matrix in system (21)-(23) with k -- 1 is positive definite, one can use continuation or homotopy methods, among others. Here we describe a simple continuation method [e.g. 15-17]. Let us symbolically express system (18)-(20) in the form F(u h, tab, ph; Re) = O. Here Re = 1/v denotes the target Reynolds number, i.e. the value of the Reynolds number at which we want a solution. Now, suppose we have a sequence of increasing Reynolds numbers {Re," }~= t with ReM = Re. We denote the solution of system (18)-(20) for vm = 1/Re,. by (Urn,tam, p,"). For any value of m, we obtain (Urn,tom,p,") by Newton's method, i.e. we solve the sequence of linear systems
F,(u~-i,,to~-l),p~ l);Rein). ((u~)'tom.(k~,,~k)~,~,," ,_ ,'mt"~k-l), tO~- ,),ptk m l))) =-F(u~-~),to~
I),p~-';Re,")
fork=l,2
.....
(24)
Here, F' denotes the Jacobian of F with respect to (uh, to h, ph). We need to specify the initial guesses (u~),Wm''c°),.m"~°)~:to be used to start, for each m, the iteration with respect to k in iteration (24). Assume that Re~ is sufficiently small so that the iteration (24) converges if (u~°), to~°~, p~0)) is chosen to be the solution of a discrete linear Stokes problem. For example, we could choose Rej = 1 and take for (u~°), to]°),pl°) ) the solution of system (18)-(20) with all cross-product terms deleted and with v -- I. The latter problem involves a symmetric, positive definite linear algebraic system. The remaining initial guesses {(u~ ), to~), Fm ,~t°)~M )'Jm=2 are determined by "continuing along the tangent", i.e. by solving the linear algebraic system F'((Um-l),C"Om ,,Pm-l', Re,. - I) " ((u(m°),tf'l (0)-m , t"n(°) ," (~U -- , " - : ' 1, ~Om-l,P,.-I)) =-(Rez-Re,._t)FRe((um_~,co,._,,p,.
~;Re,,_~)).
(25)
554
PAVEL B. Boch~v and MAX D. GUNZBURGER
Here, (u,,_ ~, to,. ~,p ..... ~) denotes the converged iterates determined from iteration (24) at the (m - 1)th stage and FRo denotes the Frechet derivative of F with respect to the parameter Re. Note that the coefficient matrix of linear system (25) may be chosen to be the same as the coefficient matrix for the last iteration sequence of (24) at the (m - 1)th stage. The combined Newton-continuation method is now completely defined. If (Rein - Re,, ~) is sufficiently small, the use of system (25) should yield initial guesses that are within the attraction ball of Newton's method and such that the coefficient matrices in system (24) with k = 1 are positive definite. (They are always symmetric, of course.) In fact, since Newton's method is guaranteed to be locally convergent, i.e. its attraction ball is nontrivial, and since the neighborhood of a minimizer for which the Hessian matrix is positive definite is also nontrivial, by choosing ( R e , . - Rein ~) sufficiently small, one can guarantee that the combined Newton-continuation method should only deal with symmetric, positive definite matrices. The method can be made self-correcting. For example, suppose that the linear systems in equations (24) are solved by an iterative method, e.g. the conjugate gradient method, that works, or works well, only for symmetric positive definite matrices. Then, if we have chosen an increment (Re,. - Re,._ ~) that is too large, then either the Newton iteration or the conjugate gradient iteration will fail. In either case, one can restart the mth stage by choosing a smaller value of Rein in systems (24) and (25). We note that often the even simpler "continuation along a constant" method ,
,F"
I--~
"-l,tom-I
Pro-l)
(26)
for generating initial guesses has been found to work well in viscous flow calculations [18]. The disadvantage of equation (26) is that it breaks down in the vicinity of bifurcation or turning points, while system (25) can be amended so that it can handle such singular points [15-17]. We also note that there are modifications possible to Newton's method that circumvent the difficulty of not having positive definite Hessian matrices. For example, one simple such modification is to add to the diagonal entries of the coefficient matrix in system (21)-(23) a multiple of the magnitude of the residual of the previous Newton iterate. In the notation of equations (24), we would replace the Jacobian matrix F' on the left-hand side with ,,~k -,7., Re")l!, F'(u~-,7, --mt"(k-l;,r,, ntk "" Re,.) + ~ I[ F(y~ '), ta~- u ,~-m
(27)
where l[ ' L[denotes the Euclidean length. By choosing the constant V sufficiently large, we can make sure that the matrix represented by formula (27) is positive definite. However, as we approach the solution of the problem, the term multiplying V approaches zero, so that we recover the local convergence properties of Newton's method in the neighborhood of the solution. 3.3. Enhancing mass conservation
In many instances one is especially interested in conserving mass, i.e. satisfying the continuity equation (1) as well as possible. The method discussed here is not exactly mass conserving, i.e. div u h 4= 0 exactly. In fact, one can easily show that, at best, div u h ~ Ch r whenever the finite element space, restricted to any element, contains all polynomials of degre ~ 0 is a constant, then the constant C can be shown to be proportional to 1/x/~. Thus, by choosing a large value of ~, one can make div u h small. This is another potential advantage of the least-squares method. However, one must keep in mind that the larger the value of ~, the worse one satisfies the momentum equation relative to the continuity equation. 3.4. Theoretical observations concerning accuracy
Error estimates for the least-squares finite element approximations are derived through the following process. First, using the theory of Ref. [19] (see also Refs [1, 20], one can show that
Accuracy of least-squares methods for the Navier-Stokes equations
555
the error estimates for the nonlinear Navier-Stokes equations are essentially the same as those for the linear Stokes problem, at least away from singular points. Now for the Stokes problem divu=O
in [1,
(28)
in f~,
(29)
v curl ta + gradp = f
in f~
(30)
divta=0
in f~,
(31)
curlu-ta=0
standard techniques for estimating the error of least-squares finite element approximations require that this system be elliptic and that the boundary conditions satisfy the complementing condition of Ref. [21]; see, for example, Refs [6, 7, 22]. Moreover, in order to use finite element functions that are merely continuous across element boundaries, one must be able to bound the L2-norm of the derivatives of solutions of equations (28)-(31) in terms of the data of the problem. For the boundary conditions (5)-(7) this program cannot be carried out, i.e. system (28)-(31) with conditions (5)--(7) does not satisfy the complementing condition of Ref. [21] when one requires that the derivatives of u, ta and p be merely square integrable. Thus, standard error estimation techniques for least-squares finite element discretization of equations (1)-(7) cannot be used. On the other hand, one can easily show that system (28)-(31) with conditions (8) and (9) does satisfy all the requisite conditions of the theory of Ref. [21] and, therefore, one may conclude that in that case one may obtain optimal error estimates. For example, if continuous piecewise polynomials of degree ~
lu - uhh + IoJ - tobit + [p --phh <~Chr;
(32)
and, under mild restrictions on the domain, Ilu - u*llo + IIoJ - ¢ohllo+ liP -Philo ~ Ch r+ ~,
(33)
where h denotes a measure of the grid size and where II " II0 and l" I~ denote the L2(~)-norm and the L2(f~)-norm of the first derivatives, respectively. Nothing we have said necessarily means that least-squares finite element approximations of equations (1)-(7) do not satisfy conditions (32) and (33) as well. All we have concluded so far is that we cannot, using standard estimation techniques, derive these estimates. Indeed, this is the motivation for the computational study of Section 4.
3.5. The connection of conditions (8) and (9) with second-order elliptic problems Consider the problem Ap = g
in t~,
(34)
where A denotes the Laplacian operator and p=P
onF.
(35)
Now, let u=gradp
in fL
(36)
divu=g
in f~.
(37)
Then we have that
Thus, equations (36) and (37) are a first-order system equivalent to equation (34). However, it can be shown that this first-order system is not elliptic and that, in general, least-squares finite
556
PAVELB. BOCHEVand MAX D. GUNZBURGER
element methods for equations (36) and (37) are not stable [23]. However, if one considers system (37), curlta + g r a d p = u ,
c u r l u = 0 and divta = 0
in fL
(38)
then it can be shown that this system is elliptic. Clearly, i f p is a solution of system (37) and (38) with condition (35), then p is also a solution of system (34) with condition (35). Note that the principal part of system (37) and (38), i.e., the differentiated terms, is identical to those of system (1)-(4), and that the boundary condition (35) is the same as equation (8). We can then infer the optimal accuracy of least-squares finite element methods for appropriately defined first-order systems arising from second-order elliptic equations. 4. A C O M P U T A T I O N A L S T U D Y OF A C C U R A C Y The accuracy of least-squares finite element approximations for the Navier-Stokes equations is essentially the same as that for the linear Stokes equations. Moreover, the accuracy is independent of dimension. Thus, here we restrict attention to the Stokes equations in 2-D. We also take for our domain the unit square f~ = {0 ~
Ou
8v
~ + ~yy=
t~v 0u - -Ox 8y
-
•
=
g~ i n t'l,
(39)
g2 in fL
(40)
Oo~ 8p
--7 ~!) + ~x = f i in f~
(41)
and
ap - - 0--x~ +
~yy=f2 in
f~,
(42)
where g~, g2, fi and f2 are given functions. We consider the two sets of boundary conditions: u = U and v = V on F and p(0,0) = P0; L2 Error
for U;
(BCI)
L2 E r r o r
BCI/BC2
for V; BCI/BC2
o
0.01 0.005 0.00 0.0005
°'°°°1I
0.0001 2
3
L2 E r r o r
5
7
for W;
I0.
. 2
15.
BCI/BC2
0.01 0.001
0.01
.
.
.
3
l
. 5
L2 E r r o r 0.1
0.1
.
7
10.
15.
for P; BCI/BC2
. i
0.0001 0.001
0.00001 2
3
5
7
i0.
15.
2
3
Fig. 1. L 2 errors vs number of grid intervals in each direction for example (44)'. ..... , condition (BC2).
5
7
I0.
15.
, condition (BCI);
Accuracy of least-squares methods for the Navier-Stokes equations
HI Error for U; BCI/BC2
HI Error for V; BCI/BC2
0.5
0.5
0.2
0.2 0.i 0.05
0.1 0.05
557
0.02 0.01
0.02 0.01
2
3
5
7
I0.
15.
2
HI Error for W; BCI/BC2 2 1 0.5
3
5
7
i0.
15.
HI Error for P; BCI/BC2
°o15
O.
-.. "'..
0.2 0.i
.
"-,..
0.01
""-,,,.,.
0.005 0.001 2
3
5
7
i0.
15.
2
3
5
Fig. 2. H' errors vs number of grid intervals in each direction for example (44): . . . . . , condition (BC2).
7
i0.
15.
, condition (BCI);
and p = P a n d unl + vn2 = W on F,
(BC2)
where U, V, P a n d W are given functions defined o n F, P0 is a given n u m b e r a n d n,, n2 denote the c o m p o n e n t s of the u n i t outer n o r m a l . We will define the various data functions by choosing a n exact solution (u, v, o~,p) a n d then substituting into the above equations. We will measure the differences (L 2 error)
L2 Error for U; BCI/BC2 0.001 0.0005
0.001 0.000!
0.0001 0.00005
0.0003 0.00005
0.00001
Error
for
V;
2
3
5
7
BCI/BC2
0.00001
2
3
5
7
i0. 15.
L2 Error for W; BCI/BC2 0
L2
.
0
i0. 15.
L2 Error for P; BCI/BC2
1
o.oi ! 0.001
0.001
0.0001
0.0001
0.00001 2
,
,
,
3
5
7
"'.'"'""~,, i0. 15.
0.00001
,
2
g
,
.......
"> 16
Fig. 3. i 2 errors vs number of grid intervals in each direction for example (45) with s = 2.5: - condition (BC1); . . . . . , condition (BC2).
.".
15.
PAVEL B. BOCHEV and MAX D. GUNZBURGER
558
0.05 0.02 0.01 0.005 0.002 0.001:
HI Error for U; BCI/BC2
HI Error for V; BCI/BC2 0.02
00 . 0 1 0.005
i
0.002
0.001
2
3
5
HI Error
7
i0.
2
15.
5
7
i0.
0.1 0.05
1 .
0.021 o.ol I
' ........... .....
o.oo51
.... . .
00:0012[
2
3
5
7
i0,
15,
HI Error for P; BCI/BC2
for W; BCI/BC2
0.2 0.I 0.05 0.02 0.01 0.005 0.002 0.001
3
15.
.... .
2
3
5
7
.....)
I0.
15.
Fig. 4. H I errors vs number of grid intervals in each direction for example (45) with s = 2.5: - condition (BCI), . . . . . , condition (BC2).
and (H l error)
[¢
- -
~h[1 =
~X (¢
-
-
~h)
(~
..[_
_
~h)
dfl
,
w h e r e ~ c o u l d be a n y o f u, v, ~ , o r p. T h r o u g h o u t we use picccwise q u a d r a t i c finite e l e m e n t spaces; thus, for sufficiently s m o o t h s o l u t i o n s a n d for the d o m a i n f l = {0 ~< x ~< 1, 0 ~< y ~< 1}, we expect t h a t for the b o u n d a r y c o n d i t i o n s (BC2) we have t h a t
ll~
-
-
~hll0 = O(h3) a n d [~
-
~hh = O(h2),
L2 Error for V; BCI/BC2
L2 Error for U; BCI/BC2 0.00005
0
o.oooo~ 5. I0
0.0000 5. i0
-6 i. i0-7 5, i0
i. 105. i0 2
3
5
7
i0. 15,
L2 Error for W; BCI/BC2
.
0
2
0.00001 _61
"'....
0
0
3
0
Fig. 5.
3
L 2 errors
5
7
i0, 15.
5
5
7
10.
15,
2
0
1
........................
0.000I0-6 011.
i. i0 2
0
L2 Error for P; BCI/BC2
0 . 0.0001
0.0001
(43)
3
5
7
i0.
vs number of grid intervals ineachdirecfion ~rexampM ~5) w i t h s = 2 . 0 h - - , condition (BCI); . . . . . ,condition(BC2).
15.
559
Accuracy of least-squares methods for the Navier-Stokes equations
HI E r r o r
H1
for U; BCI/BC2
Error
for V;
BCI/BC2
5
i0.
0.001
0.001 0.0005
0.0005
0.0002
0.0002
0.0001
0.000
O. 00005
~,.. 2
3
5
H1 E r r o r
7
i0.
0.00005
15.
2
HI E r r o r
for W; BCI/BC2
0.005
0.005
0.001 0.0005
0.001 0.0005
0.0001 0.00005
"--
2
3
5
7
i0.
3
7
15.
for P; B C I / B C 2
0.0001 0.00005
"x
15.
2
3
5
7
i0.
15.
Fig. 6. H ~ errors vs number of grid intervals in each direction for example (45) with s = 2.01: - condition (BC1); . . . . . , condition (BC2).
w h e r e a g a i n ~ c o u l d be a n y o f u, v, o~ o r p. T h i s is c o n f i r m e d by the c o m p u t a t i o n s t h a t f o l l o w . O f special i n t e r e s t to us h e r e a r e the c o m p u t a t i o n a l l y d e t e r m i n e d rates o f c o n v e r g e n c e for the b o u n d a r y c o n d i t i o n s ( B C 1). T h e first e x a m p l e we p r e s e n t has the s m o o t h e x a c t s o l u t i o n u=
-- c o s n x s i n r~y + l - - y 3,
~o
L2 Error
av
-
v = s i n g x c o s Try + l -- x 3,
~u
~x
~ y , a n d p = sin y cos x + x y 2.
for U;
BCI/BC2
0.005
(44)
L2 Error
for V;
BCI/BC2
0.005
0.001 0.0005
00-001
.0005 k
0.0001 0.00005
0.0001 0.00005 2
3
L2 Error
5
7
for W;
i0.
15.
2
BCI/BC2
3
5
7
i0.
15.
L2 Error
for P; BCI/BC2
2
5
0.i 0.01
0.01
0.001
0.001
0.0001
\
2
3
5
7
i0.
15.
0.0001 3
7
i0.
Fig. 7. L 2 errors vs number of grid intervals in each direction for example (45) with s = 1.5: - condition (BCI); . . . . . , condition (BC2).
15.
560
PAVEL B. BOCHEV and MAX D. GUNZBURGER
HI E r r o r
for U;
BCI/BC2
0.05
HI E r r o r 0
0.02
.
0
5
for V;
BCI/BC2
~
0.02
\
0.01
0.01~
"~" '4
,
0.005
\/•
"4
0.0051 i
2
3
5
HI E r r o r
7
i0.
for W;
15.
BCI/BC2
0.2
--
0.1 0.05
3
5
7
10.
15.
HI E r r o r
for
P; B C I / B C 2
2
5
7
0.2 0.1 0.05
0.02
0.02 0.01 o
2
~ ~
0.01
.oo5!
0.005 2
3
5
7
i0.
15.
3
I0.
15.
Fig. 8. H ~ errors vs number of grid intervals in each direction for example (45) with s = 1.5: - - - , condition (BC1); ..... , condition (BC2). Figures 1 and 2 show log-log plots of the L 2 and H t errors vs the number of grid intervals in each direction for a uniform grid spacing. (Similar results have been obtained for nonuniform grid spacings.) F r o m these figures one can confirm that for the boundary conditions (BC2) one does indeed obtain the best approximation convergence rates of solution (44) (As is usually the case, computed L 2 rates are usually less reliable than their H t counterparts.) In fact, the asymptotic slopes of each of the dotted curves of Fig. 1 are approximately 3, while of Fig. 2 are approximately 2. Surprisingly, the slopes of the solid curves are asymptotically the same as those of the dotted curves; thus, for the boundary conditions (BC1) we seemingly obtain the best approximation rates of convergence of solution (44). Experiments with other smooth solutions yield similar results. L2 E r r o r
for U;
BCI/BC2
0.01 i 0. 005 ~
L2 E r r o r
for V;
BCI/BC2
2
5
i0.
0.02 0.01 0.005 0.002 0.001 0.0005
0.002 0.001 0.0005
0.0002
0.0002 2
3
L2 E r r o r
5
7
for W;
i0.
15.
BCI/BC2
0.i 0.05
3
7
15.
L2 E r r o r
for P; B C I / B C 2
2
5
0.1
0.01 0.005
0.01
0.001 0.00051
0.001 2
3
5
7
i0.
15.
3
7
i0.
15.
Fig. 9. L 2 errors vs number o f grid intervals in each direction for example (45) with s = 1,00h - - condition (BC1); . . . . . , condition (BC2).
561
A c c u r a c y o f l e a s t - s q u a r e s m e t h o d s for t h e N a v i e r - S t o k e s e q u a t i o n s
HI Error
for U;
BCI/BC2
HI E r r o r
0.2 0.15
0.2 0.15
0.i 0.07 0.05
0.i 0.07 0.05 x "'^~
0.03
0.02
0.02 3
HI
5
Error
7
i0.
for W;
BCI/BC2
^
0.03 2
for V;
15.
~,!,,
2
BCI/BC2
3
5
HI Error
0.5;
10.
7
V 15.
for P; BCI/BC2
.
0.2
0.
0.i
0.
0.05
0.0 2
3
5
7
i0.
i
15.
7
i0.
is.
Fig. 10. H t e r r o r s v s n u m b e r o f g r i d i n t e r v a l s in e a c h d i r e c t i o n f o r e x a m p l e (45) w i t h s = 1 . 0 0 h - c o n d i t i o n ( B C I ) ; . . . . . , c o n d i t i o n (BC2).
The next set o f examples uses the exact solution u = v = to = p = [(x - a) 2 + (y - b)2]s/2.
(45)
We use a = b = 0.1234. The exponent s m a y be chosen to adjust the smoothness o f the exact solution, and thus the rates o f convergence o f the best approximations out o f the finite element space used. Thus, except for some exceptional integer values, we have that I[~ - (hll0 =
O(h '+')
and
I~ - ~hh =
O(ht),
where t = min(s, 3).
(46)
and where again ~ denotes any o f u, v, co or p, and (h denotes the corresponding best approximation. Figures 3-10 are l o g - l o g plots o f L ~ and H j errors vs the n u m b e r o f grid intervals in each direction for a uniform grid spacing and for different values o f s. According to equations (46), the Table 1. Rates of convergence of the H ~ and L 2 errors in the best approximation (b.a.) and in the least-squares finite element solution (u, v, to, p )
Exact solution
Boundary condition type
(44)
BCI BC2
(45), s = 2.5
BC2 BC2
(45), s = 2.02
BC1 BC2
(45),s = 2.5
BCI BC2
(45),s = 2.001
BC2 BC2
CAF 22-4/~-K
Rates of convergence Error type
b.a.
u
v
~o
p
Ht L2 H~ L2 H~ L2 HI L2 H~ L2 Hi L2 H~ L2 H~ L: HI L2 H~ Lz
2 3 2 3 2 3 2 3 2 3 2 3 2.5 2.5 2.5 2.5 1.001 2.001 1.002 2.002
1.97 3.54 2.98 3.10 2.93 3.00 1.93 2.96 1.98 2.82 1.98 2.93 1.41 2.39 1.39 2.39 0.99 1.87 0.98 1.92
1.96 3.61 2.98 3.11 1.93 2.98 1.93 2.97 2.98 2.81 1.98 2.95 1.41 2.42 1.39 2.42 0.99 1.87 0.98 1.94
2.02 3.44 1.98 3.04 1.88 2.88 1.90 2.89 1.86 2.70 1.84 2.71 1,82 2.52 1.43 2.52 1.34 1.98 0.97 2.06
2.09 3.21 2.00 3.00 1.89 2.92 1.94 2.94 1.85 1.98 1.76 2.87 1.82 2.36 1.41 2.36 1.37 1.81 0.97 1.89
562
PAVELB. BOCHEVand MAX D. GUNZBURGER
best approximation for the cases s = 2.5 and 2.001 have L 2 and H ~ errors of O ( h 3) and O(h2), respectively. Figures 3-6 show that the least-squares finite element solutions for both boundary conditions (BC1) and (BC2) also seem to converge at approximately this rate for these values of s. Similar conclusions can be drawn from Figs 7 and 8, for which the best approximations and both finite element solutions have L 2 and H ~ errors close to O (h 25) and O (h ~5), respectively. Likewise, from Figs 9 and 10 we have that finite element solutions have L 2 and H ~errors close to O ( h 2) and O(h~), respectively; these rates are again those of the best approximation for s - 1.001~ The above conclusions drawn from Figs 1-10 can be quantized by computing the slope of a least-squares straight-line fit to the various curves in the figure. The results are given in Table 1~ The first column of Table 1 gives the exact solution used to obtain the remaining columns. The second column gives the error type, i.e. either H j or L ~, to which the entries of the remaining columns correspond. The third column gives the rate of convergence of the best approximations to the exact solutions out of the quadratic finite element space. Specifically, that column gives the value o f t for the H ~error and t + 1 for the L 2 error as determined from equations (46). The fourth column gives the exponent fl in the error formula ]iu - uhll = O(h'), where II • II denotes either the L 2 n o r m or the H ~ semi-norm, u is the exact x - c o m p o n e n t of the velocity and u h is its least-squares finite element approximation. The exponent fl is determined computationally from the data that was used to produce the plots for the errors in the approximation to u in Figs 1-10. Specifically, fl was determined as the slope of the best least-squares straight-line fit of the curves in those figures. The remaining three columns give similar data for v, o and p. Note that the behavior of the L ~ errors is much more erratic than that of the H ~ errors; this is usually the case. However, we see that the errors for condition (BCI) are not all that different from those for condition (BC2) and we may conclude that the errors in the former case are, at the least, nearly optimal. 5. C O N C L U D I N G
REMARKS
We have studied the accuracy of a least-squares finite element method for the Navier-Stokes equations based on a velocity-vorticity-pressure formulation. In particular, we have focused on velocity boundary conditions. In this case, standard techniques for error estimation fail so that a computational study of accuracy is called for. In spite of this failure, the computational experiments reported here indicate that the accuracy is, at the least, nearly optimal, i.e. that the rate of convergence as the grid size is reduced is seemingly the same as that for the best approximations out of the finite element spaces. When the least-squares finite element discretization algorithm is coupled with a Newton linearization with continuation, desirable discrete algebraic properties result. Thus, the overall method seems to provide a good combination of accuracy and efficiency. There remains substantial issues to study connected with the least-squares finite element method for incompressible flows. These include practical implementation issues such as the use of iterative linear solvers and theoretical issues such as the derivation of rigorous error estimates. We will address these issues in a forthcoming paper. Acknowledgements--This work was supported in part by the Air Force Office of Scientific Research under Grant No.
AFOSR-90-0179. M. D. Gunzburger was also supported by the Institute for Computational Methods in Propulsion at the NASA Lewis Research Center. REFERENCES I. V. Girault and P.-A. Raviart, Finite Element Methods for Navier-Stokes Equations. Springer-Verlag, Berlin (1986). 2. M. Gunzburger, Finite Element Methods for Viscous Incompressible Flows. Academic Press, Boston, MA (1989). 3. M. Gunzburger and R. Nicolaides (Eds), Incompressible Computational Fluid Dynamics: Trends and Advances. Cambridge Univ. Press, Cambridge (1993). 4. L. Franca, T. Hughes and R. Stenberg, Stabilized finite elements. In Incompressible Computational Fluid Dynamics: Trends and Advances (Edited by M. Gunzburger and R. Nicolaides). Cambridge Univ. Press, Cambridge (in press). 5. C. Chang, Least squares finite-element method for incompressible flow in 3-D (1993). 6. C. Chang, Piecewise linear approach to the Stokes equations in 3-D (in press).
Accuracy of least-squares methods for the Navier-Stokes equations
563
7. C. Chang and B.-N. Jiang, An error analysis of least-squares finite element methods of velocity-pressure-vorticity formulation for the Stokes problem. Comput. Meth. Appl. Mech. Engng 84, 247 (1990). 8. B.-N. Jiang, A least-squares finite element method for the incompressible Navie~Stokes equations. Int. J. Numer. Meth. Fluids 14, 843 (1992). 9. B.-N. Jiang and C. Chang, Least-squares finite elements for the Stokes problem. Comput. Meth. Appl. Mech. Engng 78, 297 (1990). 10. B.-N. Jiang, T. Lin and L. Povinelli, Large-scale computation of incompressible viscous flow by least-squares finite element method (in press). I I. B.-N. Jiang and L. Povinelli, Least-squares finite element method for fluid dynamics. Comput. Meth. Appl. Mech. Engng 81, 13 (1990). 12. B.-N. Jiang and V. Sonnad, Least-squares solution of incompressible Navie~Stokes equations with the p-version of finite elements (in press). 13. D. Lefebvre, J. Peraire and K. Morgan, Least squares finite element solution of compressible and incompressible flows (in press). 14. L. Tang and T. Tsang, A least-squares finite element for time-dependent incompressible flows with thermal convection (in press). 15. H. Keller, Numerical Methods in Bifurcation Problems. Springer-Verlag, Berlin (1987). 16. H. Keller, Global homotopies and Newton methods. In Recent Advances in Numerical Analysis (Edited by C. de Boor and G. Golub), pp. 73-79. Academic Press, New York (1978). 17. W. Rheinboldt, Solution fields of nonlinear equations and continuation methods. SIAM Jl Numer. Analysis 17, 221 (1980). 18. M. Gunzburger and J. Peterson, Predictor and steplength selection in continuation methods for the Navier-Stokes equations. Computers Math. Applic. 22, 72 (1991). 19. F. Brezzi, J. Rappaz and P.-A. Raviart, Finite-dimensional approximation of nonlinear problems, Part I: branches of nonsingular solutions. Numer. Math. 36, 1 (1986). 20. M. Crouzeix and J. Rappaz, On Numerical Approximation in Bifurcation Theory. Springer-Verlag, Berlin (1990). 21. S. Agmon, A. Douglis and L. Niremberg, Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions I, II. Commun Pure Appl. Math. 12, 623 (1959); 17, 35 (1964). 22. W. Wedland, Elliptic Systems in the Plane. Pittman, London (1979). 23. G. Fix, M. Gunzburger and R. Nicolaides, On finite element methods of the least squares type. Computers Math. Applic. 5, 87 (1979).
ACCURACY
0045-7930/93 $6.00+ 0.00 Copyright © 1993 Pergamon Press Ltd
OF LEAST-SQUARES METHODS NAVIER-STOKES EQUATIONS
FOR
THE
PAVEL B. BOCHEV a n d MAX D . GUNZBURGER Department of Mathematics and Interdisciplinary Center for Applied Mathematics, Virginia Polytechnic Institute and State University, Blacksburg, VA 24061-0531, U.S.A.
(Receit,ed It February 1993) Abstract--Recently there has been substantial interest in least-squares finite element methods for velocity--vorticity pressure formulations of the incompressible Navier-Stokes equations. The main cause for this interest is the fact that algorithms for the resulting discrete equations can be devised which require the solution of only symmetric, positive definite systems of algebraic equations. On the other hand, it is well-documented that methods using the vorticity as a primary variable often yield very poor approximations. Thus, here we study the accuracy of these methods through a series of computational experiments, and also comment on theoretical error estimates. It is found, despite the failure of standard methods for deriving error estimates, that computational evidence suggests that these methods are, at the least, nearly optimally accurate. Thus, in addition to the desirable matrix properties yielded by least-squares methods, one also obtains accurate approximations.
1. I N T R O D U C T I O N
The approximate solution of the Navier-Stokes equations of incompressible flow has received tremendous attention from engineers and mathematicians [e.g. 1-3]. Among the more recent developments has been the use of least-squares ideas; see, for example, Ref. [4] for a recent survey of one such approach. Also, truly least-squares methods have been developed and applied [e.g. 5-14]. Here, a finite element method based on a least-squares variational principle is examined for the approximate solution of the stationary, incompressible Navier-Stokes equations. These equations are cast into a first-order system of partial differential equations involving the velocity, vorticity and pressure as dependent variables. In 3-D one has 7 unknown scalar fields. However, the application of a least-squares principle along with, for example, a Newton linearization, results in symmetric, positive definite linear systems, at least in a neighborhood of the solution. The influence of the Reynolds number on the positive definiteness of these linear systems is felt only through the size of the neighborhood. Thus, if properly implemented continuation (with respect to the Reynolds number) methods are used, one can expect to only encounter symmetric, positive definite linear systems in the solution procedure. A further advantage of this method is that a single piecewise polynomial finite element space may be used ,for all test and trial functions, i.e. one may use equal-order interpolation with respect to a single grid for all dependent variables and test functions. A final advantage resulting from the use of a least-squares principle is that, unlike some other methods involving the vorticity, no art!)Cicial numerical boundary conditions for the vorticity need be devised. The method discussed here is similar to those in Refs [5-12]; however, there are certain differences as well, especially in the order in which the least-squares, the discretization and the linearizations steps are taken. Furthermore, the analyses found in some of these papers are incorrect, leaving open the question of the accuracy of approximations. In Section 2, we define the least-squares finite element method. In Section 3, we discuss some practical and theoretical issues connected with the method described in Section 2. Then, in Section 4, we give the results of a computational study of the accuracy of the algorithm. Finally, in Section 5 we give some concluding remarks. 549
550
PAVEr. B, BOCHEV and MAX D. GUNZBURGER
2. T H E L E A S T - S Q U A R E S F I N I T E E L E M E N T M E T H O D
2.1. The velocity-vorticity-pressure equatzons Let the bounded set f~cN3 denote the flow domain and let F denote its boundary. The dimensionless equations governing the steady incompressible flow of a viscous fluid may be written in the form divu=0 curlu-to=0 vcurlto+to xu+gradp=f
in~,
(l~
inf~
(2)
in ~'2
(3)
in fL
(4)
and divto=0
where u, p and to denote the velocity, pressure and vorticity fields, respectively, v is the inverse of the Reynolds number and f is a given body force. The first-order system of partial differential equations (1)-(4) constitutes the velocity-vorticity-pressure form of the equations of steady, incompressible flow. Note that the terminology "pressure" is a little misleading since the variable p is actually the total head, i.e. p =/~ + ½1ul2, where/~ denotes the pressure. In view of equation (2), it seems that equation (4) is redundant. However. it is crucial to the stability of the least squares algorithm to explicitly require equation (4) [e.g. 5, 10]. The system (1)-(4) should be supplemented with boundary conditions. Here we examine two such boundary conditions. The first is to impose the velocity on the boundary, i.e. u=U~
on F,
(5)
where U1 denotes a given function defined along F. Note that equation (5) implies that the normal component of the vorticity is also known, i.e., to.n=n.curlU~
on F,
(6)
where n denotes the unit outer normal to ~q. To see this, one merely needs to observe that n • curl u involves only tangential derivatives of the tangential components of u and these may be deduced from equation (5). To these one must add a condition to fix the pressure; we choose to fix the pressure at a single point x0 in ~, i.e. p(x0) = Po,
(7)
where P0 is a given number. The second combination of boundary conditions we consider is the pressure and the normal component of the velocity, i.e. p = P
on F
(8)
and u-n=U;
on F,
(9)
where P and U2 denote given functions defined along F. The boundary conditions (8) and (9) are not all that useful as boundary conditions for the Navier-Stokes equations; we consider them here for two reasons. First, a complete, rigorous analysis of least-squares finite element approximations of equations (1)-(4), (8) and (9) can be given using standard techniques; this is not the case for equations (5)-(7). Furthermore, equations (1)-(4), (8) and (9) can be shown to be related to second-order elliptic partial differential equations; we will discuss these issues in more detail below.
2.2. The least-squares principle One can use equations (1)-(4) to define the least-squares functional
~-(u, to~ P) = f ( I c u r l u - to 12+ Idiv ul2 + Idiv to 12+ Iv curl to + to x u + grad p -- t]2) d~. J.
(lO)
551
Accuracy of least-squares methods for the Navier-Stokes equations
The least-squares principle then requires the minimization of this functional over appropriate function spaces. Then standard techniques from the calculus of variations may be used to deduce that the minimizers (u, p, o) of Y necessarily satisfy
s
[(curlu-o).curlv+divudivv+(vcurlo+gradp+w n
xu-f)*(o
xv)]dQ=O,
(11)
[divodivr-(curlu-o).i+(vcurlo+gradp+wxu-f).(vcurlr+rxu)]dn=O sn (12) and
s
(vcurlo+gradp+o R
xu-f)*gradqdfi=O.
(13)
In equations (1 1)-(13) the test functions (v, q, c) are required to belong to suitable function spaces; we do not go into detail here since we are primarily interested in finite element discretizations of these equations. Clearly any solution (u, 0,~) of, say equations (l)--(7), satisfies equations (1 lH13). 2.3. The 2-D case For planar flows we have that u = (u,, ul, O)Tand u, , u2 and p are functions of x, and x2 only. Then, we have that w = (0,0,w)’ where w = &,/ax, - au, /ax,. In this case, system (l)-(4) simplifies to $+$=O I au2 ----
84
ax,
ax,
o=o
am ap
vZ+--u2w=f, 2 ax,
am
ap
-vz+-+u,w=f, I
in Q,
(14)
in R,
(15)
in R
(16)
in a,
(17)
2
ax2
where f = (f, , f2, O)T.The boundary conditions (5)--(7) reduce to just equations (5) and (7) and the boundary conditions (8) and (9) remain unchanged. The functional (10) and the necessary conditions (1 lH13) also simplify in the obvious manner for 2-D problems. 2.4. Finite element methods Starting with the weak formulation (11 )(13), a conforming finite element method can be defined in a completely standard manner. We choose a finite space Sh parametrized by h. For example, for a given positive integer r, Sh could consist of continuous (over Q) piecewise polynomials of degree
on r}, on r}
552
PAVELB. BOCHEVand MAX D. GUNZBURGER
F o r the b o u n d a r y conditions (5)-(7), the discrete problem is defined as follows: seek u ~'~ S ~, to~ ~ S ~' and p~ e S ~ such that u ~ = U h and to~ . n = W h on F, ph(x°) = P0,
f [(curl u h -- toh) • curl vh + div u h div vh + (v curl to~' + grad p~ +tohXUh--f)'(to~XVh)]dfl=0
Vv~eV~,
(18)
V~heZ~
(19)
Vqh~Qh
(20)
n[div toh div ~ - (curl u h - toh). (h + (V curl toh + grad p~, +to~×u~--f)-(vcurl(~+~hxuh)]dfl=0 and
f
(vcurlto~+gradp
h+to~xu ~-f).gradqndfl=0
are satisfied. Here UIh and W h are approximations of the data U 1 and curl U~.n, respectively.
F o r example, we could define the former pair to be boundary interpolants of the latter pair. N o t e that all o f the discrete variables i.e. q~ and the c o m p o n e n t s o f u h and to~, are a p p r o x i m a t e d by the s a m e - d e g r e e continuous piecewise p o l y n o m i a l s defined with respect to a single grid.
The discrete problem for the boundary conditions (8) and (9) is also given by equations (l 8)-(20) except that now we have that u h - n = U h and ph = ph on F and the test functions spaces in equations (18)-(20) must be suitably redefined to account for the boundary conditions vh. u = 0 and q~ = 0 on F. Again, U~ and Ph are approximations to U2 and P, respectively. 3. P R A C T I C A L
AND THEORETICAL
ISSUES
3. I. N e w t o n ' s m e t h o d
The discrete equations (18)-(20) are a nonlinear system of algebraic equations that must be solved in an iterative manner. There are many methods that one might use for such a purpose; here we only consider Newton's method. However, we do remark that if a quasi-Newton method is used, it should be chosen so that it preserves symmetry and positive definiteness of the approximate Hessian matrices. Newton's m e t h o d for the solution of equations (18)-(20) is defined as follows. Given initial guesses u (°), to (°) and p(O) for u h, toh and ph, respectively, the sequence o f Newton iterates {u(k), to(k), p(k)}k>O is generated recursively by solving, for k = 1, 2 . . . . . the system fn[(curl u(k) _ to(k)), curl vh + div u (k) div yh + (v curl to (h) + grad p(k) + to~k) × U(k ]) + to(k - I) × U~k)). (to(k- ~) × Vh) + (V curl totk ~t + gradp(k
~ + took ~) × U~k-~1__ f ) . (to(kl+ Vh)]dr2
=((vcurltoIk--l)+gradptk
I)-~-2totk
I)×utkl)).(to(k-l)×vh)d~,"~
Vvh~V~,
(21)
¥~heZ~0
(22)
a[(di v totk) div ~h _ (curl u ok)-- to(kt). ~h + (v curl e) ok)+ grad p(kl + to~k~ × utk ~t + to(k ~) × U(k)). (V curl ~h) + (v curl e) ~k) + grad p~kl + to~k) × U(k- 1) + totk t) × U(k)). (~h × u~k 1)) + (V curl ~otk ~ + g r a d p fk i~+ to(k ~t × u(k- ~) _ f ) . (~h X U~k))]d92 f [ v ( f + took ~)× u~k , ) . curl~h + ( v c u r l t o ~k ~ ) + g r a d p ~ k - L ) + 2 t o k I ) × u ( k - I ) ) . ( ~ h × u ( k - ~ ) ) ] d ~
Accuracy of least-squares methods for the Navier-Stokes equations
553
and o(v curl ta ck)+ gradp (k) + ta ck) x u(k I) ..{_O.}(k- 1) X U(k)) " g r a d q* d D
= jn(f + took- t) x utk- l)). grad
qh d ~
Vq h E Qh.
(23)
The system of linear algebraic equations (21)-(23) that is used to determine the kth Newton iterate from the (k - 1)th looks rather formidable. However, it also has some very good features. First, it is easy to see that this system is symmetric. Moreover, in a neighborhood of a minimizer, the Hessian matrix for functional (10) is necessarily positive definite; but this Hessian matrix is exactly the coefficient matrix of equations (21)-(23). Thus, in a neighborhood of a solution of equations (18)-(20), the system (21)-(23) is symmetric and positive definite. This feature is independent of the value of v, i.e. of the value of the Reynolds number. These observations, along with the guaranteed local and quadratic convergence of Newton's method, are potentially very advantageous. 3.2. Continuation methods
At this point one may well ask what goes wrong with the method as the Reynolds number increases? Surely one difficulty is that the attraction ball for Newton's method, or for any other iterative method for solving the nonlinear equations, decreases in size with increasing Reynolds number. This is known to be true for other discretizations for the Navier-Stokes equations. A related observation is that the positive definiteness of the Hessian matrix is guaranteed only in a neighborhood of a minimizer; again the size of this neighborhood surely decreases with increasing Reynolds number. As a result, for an arbitrary initial guess we may have that Newton's method does not converge and/or that the coefficient matrix in system (21)-(23) is not positive definite. The former is, of course, unacceptable, while the latter occurrence would preclude the use of simple iterative methods for solving the linear systems (21)-(23) that define the Newton iterates. In order to determine an initial guess that is within the attraction ball of Newton's method, and also is such that the coefficient matrix in system (21)-(23) with k -- 1 is positive definite, one can use continuation or homotopy methods, among others. Here we describe a simple continuation method [e.g. 15-17]. Let us symbolically express system (18)-(20) in the form F(u h, tab, ph; Re) = O. Here Re = 1/v denotes the target Reynolds number, i.e. the value of the Reynolds number at which we want a solution. Now, suppose we have a sequence of increasing Reynolds numbers {Re," }~= t with ReM = Re. We denote the solution of system (18)-(20) for vm = 1/Re,. by (Urn,tam, p,"). For any value of m, we obtain (Urn,tom,p,") by Newton's method, i.e. we solve the sequence of linear systems
F,(u~-i,,to~-l),p~ l);Rein). ((u~)'tom.(k~,,~k)~,~,," ,_ ,'mt"~k-l), tO~- ,),ptk m l))) =-F(u~-~),to~
I),p~-';Re,")
fork=l,2
.....
(24)
Here, F' denotes the Jacobian of F with respect to (uh, to h, ph). We need to specify the initial guesses (u~),Wm''c°),.m"~°)~:to be used to start, for each m, the iteration with respect to k in iteration (24). Assume that Re~ is sufficiently small so that the iteration (24) converges if (u~°), to~°~, p~0)) is chosen to be the solution of a discrete linear Stokes problem. For example, we could choose Rej = 1 and take for (u~°), to]°),pl°) ) the solution of system (18)-(20) with all cross-product terms deleted and with v -- I. The latter problem involves a symmetric, positive definite linear algebraic system. The remaining initial guesses {(u~ ), to~), Fm ,~t°)~M )'Jm=2 are determined by "continuing along the tangent", i.e. by solving the linear algebraic system F'((Um-l),C"Om ,,Pm-l', Re,. - I) " ((u(m°),tf'l (0)-m , t"n(°) ," (~U -- , " - : ' 1, ~Om-l,P,.-I)) =-(Rez-Re,._t)FRe((um_~,co,._,,p,.
~;Re,,_~)).
(25)
554
PAVEL B. Boch~v and MAX D. GUNZBURGER
Here, (u,,_ ~, to,. ~,p ..... ~) denotes the converged iterates determined from iteration (24) at the (m - 1)th stage and FRo denotes the Frechet derivative of F with respect to the parameter Re. Note that the coefficient matrix of linear system (25) may be chosen to be the same as the coefficient matrix for the last iteration sequence of (24) at the (m - 1)th stage. The combined Newton-continuation method is now completely defined. If (Rein - Re,, ~) is sufficiently small, the use of system (25) should yield initial guesses that are within the attraction ball of Newton's method and such that the coefficient matrices in system (24) with k = 1 are positive definite. (They are always symmetric, of course.) In fact, since Newton's method is guaranteed to be locally convergent, i.e. its attraction ball is nontrivial, and since the neighborhood of a minimizer for which the Hessian matrix is positive definite is also nontrivial, by choosing ( R e , . - Rein ~) sufficiently small, one can guarantee that the combined Newton-continuation method should only deal with symmetric, positive definite matrices. The method can be made self-correcting. For example, suppose that the linear systems in equations (24) are solved by an iterative method, e.g. the conjugate gradient method, that works, or works well, only for symmetric positive definite matrices. Then, if we have chosen an increment (Re,. - Re,._ ~) that is too large, then either the Newton iteration or the conjugate gradient iteration will fail. In either case, one can restart the mth stage by choosing a smaller value of Rein in systems (24) and (25). We note that often the even simpler "continuation along a constant" method ,
,F"
I--~
"-l,tom-I
Pro-l)
(26)
for generating initial guesses has been found to work well in viscous flow calculations [18]. The disadvantage of equation (26) is that it breaks down in the vicinity of bifurcation or turning points, while system (25) can be amended so that it can handle such singular points [15-17]. We also note that there are modifications possible to Newton's method that circumvent the difficulty of not having positive definite Hessian matrices. For example, one simple such modification is to add to the diagonal entries of the coefficient matrix in system (21)-(23) a multiple of the magnitude of the residual of the previous Newton iterate. In the notation of equations (24), we would replace the Jacobian matrix F' on the left-hand side with ,,~k -,7., Re")l!, F'(u~-,7, --mt"(k-l;,r,, ntk "" Re,.) + ~ I[ F(y~ '), ta~- u ,~-m
(27)
where l[ ' L[denotes the Euclidean length. By choosing the constant V sufficiently large, we can make sure that the matrix represented by formula (27) is positive definite. However, as we approach the solution of the problem, the term multiplying V approaches zero, so that we recover the local convergence properties of Newton's method in the neighborhood of the solution. 3.3. Enhancing mass conservation
In many instances one is especially interested in conserving mass, i.e. satisfying the continuity equation (1) as well as possible. The method discussed here is not exactly mass conserving, i.e. div u h 4= 0 exactly. In fact, one can easily show that, at best, div u h ~ Ch r whenever the finite element space, restricted to any element, contains all polynomials of degre ~
Error estimates for the least-squares finite element approximations are derived through the following process. First, using the theory of Ref. [19] (see also Refs [1, 20], one can show that
Accuracy of least-squares methods for the Navier-Stokes equations
555
the error estimates for the nonlinear Navier-Stokes equations are essentially the same as those for the linear Stokes problem, at least away from singular points. Now for the Stokes problem divu=O
in [1,
(28)
in f~,
(29)
v curl ta + gradp = f
in f~
(30)
divta=0
in f~,
(31)
curlu-ta=0
standard techniques for estimating the error of least-squares finite element approximations require that this system be elliptic and that the boundary conditions satisfy the complementing condition of Ref. [21]; see, for example, Refs [6, 7, 22]. Moreover, in order to use finite element functions that are merely continuous across element boundaries, one must be able to bound the L2-norm of the derivatives of solutions of equations (28)-(31) in terms of the data of the problem. For the boundary conditions (5)-(7) this program cannot be carried out, i.e. system (28)-(31) with conditions (5)--(7) does not satisfy the complementing condition of Ref. [21] when one requires that the derivatives of u, ta and p be merely square integrable. Thus, standard error estimation techniques for least-squares finite element discretization of equations (1)-(7) cannot be used. On the other hand, one can easily show that system (28)-(31) with conditions (8) and (9) does satisfy all the requisite conditions of the theory of Ref. [21] and, therefore, one may conclude that in that case one may obtain optimal error estimates. For example, if continuous piecewise polynomials of degree ~
lu - uhh + IoJ - tobit + [p --phh <~Chr;
(32)
and, under mild restrictions on the domain, Ilu - u*llo + IIoJ - ¢ohllo+ liP -Philo ~ Ch r+ ~,
(33)
where h denotes a measure of the grid size and where II " II0 and l" I~ denote the L2(~)-norm and the L2(f~)-norm of the first derivatives, respectively. Nothing we have said necessarily means that least-squares finite element approximations of equations (1)-(7) do not satisfy conditions (32) and (33) as well. All we have concluded so far is that we cannot, using standard estimation techniques, derive these estimates. Indeed, this is the motivation for the computational study of Section 4.
3.5. The connection of conditions (8) and (9) with second-order elliptic problems Consider the problem Ap = g
in t~,
(34)
where A denotes the Laplacian operator and p=P
onF.
(35)
Now, let u=gradp
in fL
(36)
divu=g
in f~.
(37)
Then we have that
Thus, equations (36) and (37) are a first-order system equivalent to equation (34). However, it can be shown that this first-order system is not elliptic and that, in general, least-squares finite
556
PAVELB. BOCHEVand MAX D. GUNZBURGER
element methods for equations (36) and (37) are not stable [23]. However, if one considers system (37), curlta + g r a d p = u ,
c u r l u = 0 and divta = 0
in fL
(38)
then it can be shown that this system is elliptic. Clearly, i f p is a solution of system (37) and (38) with condition (35), then p is also a solution of system (34) with condition (35). Note that the principal part of system (37) and (38), i.e., the differentiated terms, is identical to those of system (1)-(4), and that the boundary condition (35) is the same as equation (8). We can then infer the optimal accuracy of least-squares finite element methods for appropriately defined first-order systems arising from second-order elliptic equations. 4. A C O M P U T A T I O N A L S T U D Y OF A C C U R A C Y The accuracy of least-squares finite element approximations for the Navier-Stokes equations is essentially the same as that for the linear Stokes equations. Moreover, the accuracy is independent of dimension. Thus, here we restrict attention to the Stokes equations in 2-D. We also take for our domain the unit square f~ = {0 ~
Ou
8v
~ + ~yy=
t~v 0u - -Ox 8y
-
•
=
g~ i n t'l,
(39)
g2 in fL
(40)
Oo~ 8p
--7 ~!) + ~x = f i in f~
(41)
and
ap - - 0--x~ +
~yy=f2 in
f~,
(42)
where g~, g2, fi and f2 are given functions. We consider the two sets of boundary conditions: u = U and v = V on F and p(0,0) = P0; L2 Error
for U;
(BCI)
L2 E r r o r
BCI/BC2
for V; BCI/BC2
o
0.01 0.005 0.00 0.0005
°'°°°1I
0.0001 2
3
L2 E r r o r
5
7
for W;
I0.
. 2
15.
BCI/BC2
0.01 0.001
0.01
.
.
.
3
l
. 5
L2 E r r o r 0.1
0.1
.
7
10.
15.
for P; BCI/BC2
. i
0.0001 0.001
0.00001 2
3
5
7
i0.
15.
2
3
Fig. 1. L 2 errors vs number of grid intervals in each direction for example (44)'. ..... , condition (BC2).
5
7
I0.
15.
, condition (BCI);
Accuracy of least-squares methods for the Navier-Stokes equations
HI Error for U; BCI/BC2
HI Error for V; BCI/BC2
0.5
0.5
0.2
0.2 0.i 0.05
0.1 0.05
557
0.02 0.01
0.02 0.01
2
3
5
7
I0.
15.
2
HI Error for W; BCI/BC2 2 1 0.5
3
5
7
i0.
15.
HI Error for P; BCI/BC2
°o15
O.
-.. "'..
0.2 0.i
.
"-,..
0.01
""-,,,.,.
0.005 0.001 2
3
5
7
i0.
15.
2
3
5
Fig. 2. H' errors vs number of grid intervals in each direction for example (44): . . . . . , condition (BC2).
7
i0.
15.
, condition (BCI);
and p = P a n d unl + vn2 = W on F,
(BC2)
where U, V, P a n d W are given functions defined o n F, P0 is a given n u m b e r a n d n,, n2 denote the c o m p o n e n t s of the u n i t outer n o r m a l . We will define the various data functions by choosing a n exact solution (u, v, o~,p) a n d then substituting into the above equations. We will measure the differences (L 2 error)
L2 Error for U; BCI/BC2 0.001 0.0005
0.001 0.000!
0.0001 0.00005
0.0003 0.00005
0.00001
Error
for
V;
2
3
5
7
BCI/BC2
0.00001
2
3
5
7
i0. 15.
L2 Error for W; BCI/BC2 0
L2
.
0
i0. 15.
L2 Error for P; BCI/BC2
1
o.oi ! 0.001
0.001
0.0001
0.0001
0.00001 2
,
,
,
3
5
7
"'.'"'""~,, i0. 15.
0.00001
,
2
g
,
.......
"> 16
Fig. 3. i 2 errors vs number of grid intervals in each direction for example (45) with s = 2.5: - condition (BC1); . . . . . , condition (BC2).
.".
15.
PAVEL B. BOCHEV and MAX D. GUNZBURGER
558
0.05 0.02 0.01 0.005 0.002 0.001:
HI Error for U; BCI/BC2
HI Error for V; BCI/BC2 0.02
00 . 0 1 0.005
i
0.002
0.001
2
3
5
HI Error
7
i0.
2
15.
5
7
i0.
0.1 0.05
1 .
0.021 o.ol I
' ........... .....
o.oo51
.... . .
00:0012[
2
3
5
7
i0,
15,
HI Error for P; BCI/BC2
for W; BCI/BC2
0.2 0.I 0.05 0.02 0.01 0.005 0.002 0.001
3
15.
.... .
2
3
5
7
.....)
I0.
15.
Fig. 4. H I errors vs number of grid intervals in each direction for example (45) with s = 2.5: - condition (BCI), . . . . . , condition (BC2).
and (H l error)
[¢
- -
~h[1 =
~X (¢
-
-
~h)
(~
..[_
_
~h)
dfl
,
w h e r e ~ c o u l d be a n y o f u, v, ~ , o r p. T h r o u g h o u t we use picccwise q u a d r a t i c finite e l e m e n t spaces; thus, for sufficiently s m o o t h s o l u t i o n s a n d for the d o m a i n f l = {0 ~< x ~< 1, 0 ~< y ~< 1}, we expect t h a t for the b o u n d a r y c o n d i t i o n s (BC2) we have t h a t
ll~
-
-
~hll0 = O(h3) a n d [~
-
~hh = O(h2),
L2 Error for V; BCI/BC2
L2 Error for U; BCI/BC2 0.00005
0
o.oooo~ 5. I0
0.0000 5. i0
-6 i. i0-7 5, i0
i. 105. i0 2
3
5
7
i0. 15,
L2 Error for W; BCI/BC2
.
0
2
0.00001 _61
"'....
0
0
3
0
Fig. 5.
3
L 2 errors
5
7
i0, 15.
5
5
7
10.
15,
2
0
1
........................
0.000I0-6 011.
i. i0 2
0
L2 Error for P; BCI/BC2
0 . 0.0001
0.0001
(43)
3
5
7
i0.
vs number of grid intervals ineachdirecfion ~rexampM ~5) w i t h s = 2 . 0 h - - , condition (BCI); . . . . . ,condition(BC2).
15.
559
Accuracy of least-squares methods for the Navier-Stokes equations
HI E r r o r
H1
for U; BCI/BC2
Error
for V;
BCI/BC2
5
i0.
0.001
0.001 0.0005
0.0005
0.0002
0.0002
0.0001
0.000
O. 00005
~,.. 2
3
5
H1 E r r o r
7
i0.
0.00005
15.
2
HI E r r o r
for W; BCI/BC2
0.005
0.005
0.001 0.0005
0.001 0.0005
0.0001 0.00005
"--
2
3
5
7
i0.
3
7
15.
for P; B C I / B C 2
0.0001 0.00005
"x
15.
2
3
5
7
i0.
15.
Fig. 6. H ~ errors vs number of grid intervals in each direction for example (45) with s = 2.01: - condition (BC1); . . . . . , condition (BC2).
w h e r e a g a i n ~ c o u l d be a n y o f u, v, o~ o r p. T h i s is c o n f i r m e d by the c o m p u t a t i o n s t h a t f o l l o w . O f special i n t e r e s t to us h e r e a r e the c o m p u t a t i o n a l l y d e t e r m i n e d rates o f c o n v e r g e n c e for the b o u n d a r y c o n d i t i o n s ( B C 1). T h e first e x a m p l e we p r e s e n t has the s m o o t h e x a c t s o l u t i o n u=
-- c o s n x s i n r~y + l - - y 3,
~o
L2 Error
av
-
v = s i n g x c o s Try + l -- x 3,
~u
~x
~ y , a n d p = sin y cos x + x y 2.
for U;
BCI/BC2
0.005
(44)
L2 Error
for V;
BCI/BC2
0.005
0.001 0.0005
00-001
.0005 k
0.0001 0.00005
0.0001 0.00005 2
3
L2 Error
5
7
for W;
i0.
15.
2
BCI/BC2
3
5
7
i0.
15.
L2 Error
for P; BCI/BC2
2
5
0.i 0.01
0.01
0.001
0.001
0.0001
\
2
3
5
7
i0.
15.
0.0001 3
7
i0.
Fig. 7. L 2 errors vs number of grid intervals in each direction for example (45) with s = 1.5: - condition (BCI); . . . . . , condition (BC2).
15.
560
PAVEL B. BOCHEV and MAX D. GUNZBURGER
HI E r r o r
for U;
BCI/BC2
0.05
HI E r r o r 0
0.02
.
0
5
for V;
BCI/BC2
~
0.02
\
0.01
0.01~
"~" '4
,
0.005
\/•
"4
0.0051 i
2
3
5
HI E r r o r
7
i0.
for W;
15.
BCI/BC2
0.2
--
0.1 0.05
3
5
7
10.
15.
HI E r r o r
for
P; B C I / B C 2
2
5
7
0.2 0.1 0.05
0.02
0.02 0.01 o
2
~ ~
0.01
.oo5!
0.005 2
3
5
7
i0.
15.
3
I0.
15.
Fig. 8. H ~ errors vs number of grid intervals in each direction for example (45) with s = 1.5: - - - , condition (BC1); ..... , condition (BC2). Figures 1 and 2 show log-log plots of the L 2 and H t errors vs the number of grid intervals in each direction for a uniform grid spacing. (Similar results have been obtained for nonuniform grid spacings.) F r o m these figures one can confirm that for the boundary conditions (BC2) one does indeed obtain the best approximation convergence rates of solution (44) (As is usually the case, computed L 2 rates are usually less reliable than their H t counterparts.) In fact, the asymptotic slopes of each of the dotted curves of Fig. 1 are approximately 3, while of Fig. 2 are approximately 2. Surprisingly, the slopes of the solid curves are asymptotically the same as those of the dotted curves; thus, for the boundary conditions (BC1) we seemingly obtain the best approximation rates of convergence of solution (44). Experiments with other smooth solutions yield similar results. L2 E r r o r
for U;
BCI/BC2
0.01 i 0. 005 ~
L2 E r r o r
for V;
BCI/BC2
2
5
i0.
0.02 0.01 0.005 0.002 0.001 0.0005
0.002 0.001 0.0005
0.0002
0.0002 2
3
L2 E r r o r
5
7
for W;
i0.
15.
BCI/BC2
0.i 0.05
3
7
15.
L2 E r r o r
for P; B C I / B C 2
2
5
0.1
0.01 0.005
0.01
0.001 0.00051
0.001 2
3
5
7
i0.
15.
3
7
i0.
15.
Fig. 9. L 2 errors vs number o f grid intervals in each direction for example (45) with s = 1,00h - - condition (BC1); . . . . . , condition (BC2).
561
A c c u r a c y o f l e a s t - s q u a r e s m e t h o d s for t h e N a v i e r - S t o k e s e q u a t i o n s
HI Error
for U;
BCI/BC2
HI E r r o r
0.2 0.15
0.2 0.15
0.i 0.07 0.05
0.i 0.07 0.05 x "'^~
0.03
0.02
0.02 3
HI
5
Error
7
i0.
for W;
BCI/BC2
^
0.03 2
for V;
15.
~,!,,
2
BCI/BC2
3
5
HI Error
0.5;
10.
7
V 15.
for P; BCI/BC2
.
0.2
0.
0.i
0.
0.05
0.0 2
3
5
7
i0.
i
15.
7
i0.
is.
Fig. 10. H t e r r o r s v s n u m b e r o f g r i d i n t e r v a l s in e a c h d i r e c t i o n f o r e x a m p l e (45) w i t h s = 1 . 0 0 h - c o n d i t i o n ( B C I ) ; . . . . . , c o n d i t i o n (BC2).
The next set o f examples uses the exact solution u = v = to = p = [(x - a) 2 + (y - b)2]s/2.
(45)
We use a = b = 0.1234. The exponent s m a y be chosen to adjust the smoothness o f the exact solution, and thus the rates o f convergence o f the best approximations out o f the finite element space used. Thus, except for some exceptional integer values, we have that I[~ - (hll0 =
O(h '+')
and
I~ - ~hh =
O(ht),
where t = min(s, 3).
(46)
and where again ~ denotes any o f u, v, co or p, and (h denotes the corresponding best approximation. Figures 3-10 are l o g - l o g plots o f L ~ and H j errors vs the n u m b e r o f grid intervals in each direction for a uniform grid spacing and for different values o f s. According to equations (46), the Table 1. Rates of convergence of the H ~ and L 2 errors in the best approximation (b.a.) and in the least-squares finite element solution (u, v, to, p )
Exact solution
Boundary condition type
(44)
BCI BC2
(45), s = 2.5
BC2 BC2
(45), s = 2.02
BC1 BC2
(45),s = 2.5
BCI BC2
(45),s = 2.001
BC2 BC2
CAF 22-4/~-K
Rates of convergence Error type
b.a.
u
v
~o
p
Ht L2 H~ L2 H~ L2 HI L2 H~ L2 Hi L2 H~ L2 H~ L: HI L2 H~ Lz
2 3 2 3 2 3 2 3 2 3 2 3 2.5 2.5 2.5 2.5 1.001 2.001 1.002 2.002
1.97 3.54 2.98 3.10 2.93 3.00 1.93 2.96 1.98 2.82 1.98 2.93 1.41 2.39 1.39 2.39 0.99 1.87 0.98 1.92
1.96 3.61 2.98 3.11 1.93 2.98 1.93 2.97 2.98 2.81 1.98 2.95 1.41 2.42 1.39 2.42 0.99 1.87 0.98 1.94
2.02 3.44 1.98 3.04 1.88 2.88 1.90 2.89 1.86 2.70 1.84 2.71 1,82 2.52 1.43 2.52 1.34 1.98 0.97 2.06
2.09 3.21 2.00 3.00 1.89 2.92 1.94 2.94 1.85 1.98 1.76 2.87 1.82 2.36 1.41 2.36 1.37 1.81 0.97 1.89
562
PAVELB. BOCHEVand MAX D. GUNZBURGER
best approximation for the cases s = 2.5 and 2.001 have L 2 and H ~ errors of O ( h 3) and O(h2), respectively. Figures 3-6 show that the least-squares finite element solutions for both boundary conditions (BC1) and (BC2) also seem to converge at approximately this rate for these values of s. Similar conclusions can be drawn from Figs 7 and 8, for which the best approximations and both finite element solutions have L 2 and H ~ errors close to O (h 25) and O (h ~5), respectively. Likewise, from Figs 9 and 10 we have that finite element solutions have L 2 and H ~errors close to O ( h 2) and O(h~), respectively; these rates are again those of the best approximation for s - 1.001~ The above conclusions drawn from Figs 1-10 can be quantized by computing the slope of a least-squares straight-line fit to the various curves in the figure. The results are given in Table 1~ The first column of Table 1 gives the exact solution used to obtain the remaining columns. The second column gives the error type, i.e. either H j or L ~, to which the entries of the remaining columns correspond. The third column gives the rate of convergence of the best approximations to the exact solutions out of the quadratic finite element space. Specifically, that column gives the value o f t for the H ~error and t + 1 for the L 2 error as determined from equations (46). The fourth column gives the exponent fl in the error formula ]iu - uhll = O(h'), where II • II denotes either the L 2 n o r m or the H ~ semi-norm, u is the exact x - c o m p o n e n t of the velocity and u h is its least-squares finite element approximation. The exponent fl is determined computationally from the data that was used to produce the plots for the errors in the approximation to u in Figs 1-10. Specifically, fl was determined as the slope of the best least-squares straight-line fit of the curves in those figures. The remaining three columns give similar data for v, o and p. Note that the behavior of the L ~ errors is much more erratic than that of the H ~ errors; this is usually the case. However, we see that the errors for condition (BCI) are not all that different from those for condition (BC2) and we may conclude that the errors in the former case are, at the least, nearly optimal. 5. C O N C L U D I N G
REMARKS
We have studied the accuracy of a least-squares finite element method for the Navier-Stokes equations based on a velocity-vorticity-pressure formulation. In particular, we have focused on velocity boundary conditions. In this case, standard techniques for error estimation fail so that a computational study of accuracy is called for. In spite of this failure, the computational experiments reported here indicate that the accuracy is, at the least, nearly optimal, i.e. that the rate of convergence as the grid size is reduced is seemingly the same as that for the best approximations out of the finite element spaces. When the least-squares finite element discretization algorithm is coupled with a Newton linearization with continuation, desirable discrete algebraic properties result. Thus, the overall method seems to provide a good combination of accuracy and efficiency. There remains substantial issues to study connected with the least-squares finite element method for incompressible flows. These include practical implementation issues such as the use of iterative linear solvers and theoretical issues such as the derivation of rigorous error estimates. We will address these issues in a forthcoming paper. Acknowledgements--This work was supported in part by the Air Force Office of Scientific Research under Grant No.
AFOSR-90-0179. M. D. Gunzburger was also supported by the Institute for Computational Methods in Propulsion at the NASA Lewis Research Center. REFERENCES I. V. Girault and P.-A. Raviart, Finite Element Methods for Navier-Stokes Equations. Springer-Verlag, Berlin (1986). 2. M. Gunzburger, Finite Element Methods for Viscous Incompressible Flows. Academic Press, Boston, MA (1989). 3. M. Gunzburger and R. Nicolaides (Eds), Incompressible Computational Fluid Dynamics: Trends and Advances. Cambridge Univ. Press, Cambridge (1993). 4. L. Franca, T. Hughes and R. Stenberg, Stabilized finite elements. In Incompressible Computational Fluid Dynamics: Trends and Advances (Edited by M. Gunzburger and R. Nicolaides). Cambridge Univ. Press, Cambridge (in press). 5. C. Chang, Least squares finite-element method for incompressible flow in 3-D (1993). 6. C. Chang, Piecewise linear approach to the Stokes equations in 3-D (in press).
Accuracy of least-squares methods for the Navier-Stokes equations
563
7. C. Chang and B.-N. Jiang, An error analysis of least-squares finite element methods of velocity-pressure-vorticity formulation for the Stokes problem. Comput. Meth. Appl. Mech. Engng 84, 247 (1990). 8. B.-N. Jiang, A least-squares finite element method for the incompressible Navie~Stokes equations. Int. J. Numer. Meth. Fluids 14, 843 (1992). 9. B.-N. Jiang and C. Chang, Least-squares finite elements for the Stokes problem. Comput. Meth. Appl. Mech. Engng 78, 297 (1990). 10. B.-N. Jiang, T. Lin and L. Povinelli, Large-scale computation of incompressible viscous flow by least-squares finite element method (in press). I I. B.-N. Jiang and L. Povinelli, Least-squares finite element method for fluid dynamics. Comput. Meth. Appl. Mech. Engng 81, 13 (1990). 12. B.-N. Jiang and V. Sonnad, Least-squares solution of incompressible Navie~Stokes equations with the p-version of finite elements (in press). 13. D. Lefebvre, J. Peraire and K. Morgan, Least squares finite element solution of compressible and incompressible flows (in press). 14. L. Tang and T. Tsang, A least-squares finite element for time-dependent incompressible flows with thermal convection (in press). 15. H. Keller, Numerical Methods in Bifurcation Problems. Springer-Verlag, Berlin (1987). 16. H. Keller, Global homotopies and Newton methods. In Recent Advances in Numerical Analysis (Edited by C. de Boor and G. Golub), pp. 73-79. Academic Press, New York (1978). 17. W. Rheinboldt, Solution fields of nonlinear equations and continuation methods. SIAM Jl Numer. Analysis 17, 221 (1980). 18. M. Gunzburger and J. Peterson, Predictor and steplength selection in continuation methods for the Navier-Stokes equations. Computers Math. Applic. 22, 72 (1991). 19. F. Brezzi, J. Rappaz and P.-A. Raviart, Finite-dimensional approximation of nonlinear problems, Part I: branches of nonsingular solutions. Numer. Math. 36, 1 (1986). 20. M. Crouzeix and J. Rappaz, On Numerical Approximation in Bifurcation Theory. Springer-Verlag, Berlin (1990). 21. S. Agmon, A. Douglis and L. Niremberg, Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions I, II. Commun Pure Appl. Math. 12, 623 (1959); 17, 35 (1964). 22. W. Wedland, Elliptic Systems in the Plane. Pittman, London (1979). 23. G. Fix, M. Gunzburger and R. Nicolaides, On finite element methods of the least squares type. Computers Math. Applic. 5, 87 (1979).