- Email: [email protected]
A reliable solver for nonlinear biharmonic equations
Computers& FluidsVol. 18, No. 4, pp. 405~20, 1990 Printed in GreatBritain.All rightsreserved
0045-7930/90 $3.00+ 0.00 Copyright© 1990PergamonPressple
A RELIABLE SOLVER FOR NONLINEAR BIHARMONIC
EQUATIONS ANDREAS LIPPKEI and HELMUTWAGNER2 ~Fachbereich Mathematik, Technische Universit/it Berlin and 2Dr Wagner Software GmbH, Berlin, Federal Republic of Germany (Received 29 January 1990; in revised form 22 May 1990)
Al~traet--Efficient and reliable numerical techniques of high-order accuracy are presented for solving problems for nonlinearperturbed biharmonicequations. The method is widelyapplicable,e.g. to problems of elasticity and in fluid mechanics. Here it is used to obtain accurate solutions for the driven cavity applying fourth order approximation on all convectiveterms. Solutions are calculated up to Reynolds number 30000.
I. INTRODUCTION We shall describe a numerical procedure for solving nonlinear biharmonic equations. This procedure compares very favorably with existing methods with respect to flexibility, accuracy and reliability, and should be a very efficient method in elasticity or in fluid mechanics. Briefly, the numerical method is this: the fourth order partial differential equation is discretized using central differences with three unknowns per grid point on an orthogonal net. This yields an approximation with truncation error expansion proceeding in powers of the mesh-width squared. The difference molecule involves only three grid lines in each direction for the three unknowns, the function and both the first derivatives. All but the fourth derivatives can even be discretized with fourth order error expansions. The nonlinear difference equations are solved by a sequence of Newton and Chord iterations using a modified LU decomposition. This special decomposition reduces the time and storage requirements of the procedure and allows calculations on a personal-computer. Numerical methods for nonlinear biharmonic equations can be tested on plane incompressible two dimensional flows, especially cavity flows. We applied the method to a plane cavity flow using fourth order approximations for the convective terms of the Navier-Stokes equations and a transformation to concentrate the grid lines in the boundary layers. This yields a highly accurate approximation even at high Reynolds numbers. 2. FORMULATION OF THE PROBLEM Let us consider the Dirichlet problem for the nonlinear biharmonic equation in an open set G c R 2 with a piecewise Lipschitz continuous boundary ~G. A2~k(~, r / ) + L[~b(~, r/)] = 0 ,
(~, r/) E G,
(1)
~k(~, ~/) = g,,
(¢, ~/) ~ dG,
(2)
~bn(~, r/) = g2,
(~, r/) ~ dG.
(3)
Here ft, denotes the normal derivatives of ~k with respect to the exterior normal. Let us assume that the nonlinearity L has a continuous and bounded Fr6chet derivative. In linear elasticity, ~b(~, ~7) can represent the Airy stress function or, as in the theory of thin plates, the vertical displacement due to an external force. In fluid mechanics, the equations describe the stream function formulation of the stationary Navier-Stokes equations. For the sake of simplicity we assume throughout the paper that the domain G is mapped onto a rectangle [ax, bx] x lay, by] by a simple transformation (x, y)~-*[¢ (x), r/(y)] with ~, ~/e C4(G). This assumption is not essential for the principles of the proposed numerical procedure but it reduces the notational efforts without losing the significant properties. CArls/4--~
405
406
ANDREAS LIPPKE and HELMUT WAGNER
Using ~x"= ~ x '
and
Oy..=~yy,
we derive the transformed representation of the biharmonic operator, ~f~2:=a4.odx--a2.2~x~y+ 4 2 2 a0,4 ~4 +
3 v u ~ av.~,Cgxdy, v+,a=0
(4)
and rewrite eqn (1) using the coordinates x, y in the mathematical domain and the notation/S for the transformed nonlinearity L and achieve
/%27s[x(~), y(r#)] + £Ts[x(¢), y(r#)] = O, (x, y) ~ (ax, bx) x (ay, by).
(5)
The coefficients av.~are given by the derivatives of the transformation functions ~ and q, mapping the physical domain onto the mathematical one, e.g.
(Ox'V'
a4.0 : t ~ ) • We assume that the transformed boundary conditions, e.g. on the "west" boundary can be written as a . ~ ¢ ( a x , y ) = gW(y) x -
X.[aW(y)O~qs(a~ ' Y) + b ~ ( Y ) O z ¢ ( a x , y ) = r~(Y)
(6)
with a ~ ( y ) 2 + b ' ( y ) 2 4= O. This boundary condition includes second order boundary values and is slightly more general than the one listed in eqn (3). These assumptions include not only the plane Navier-Stokes equations in stream function formulation with outflow boundary conditions but also, e.g. the axially symmetric Navier-Stokes equations. 3. D I S C R E T I Z A T I O N
For the discretization of the biharmonic operator in eqn (5) we introduce an uniform mesh with mesh sizes hx,= bx - ax N-1 '
and
hy := by - ay M-I
and grid points xi,=a,<+(i - 1)hx, y/=ay+(j
i = I(1)N,
- 1)hy, j = I(1)M.
We use the notation ~s;,s for the unknowns of the discrete system and ~s(xi, Yi) for the values of the exact solution, respectively. The centered second-order approximation of the biharmonic operator on a rectangular grid involves 13 unknowns arranged in the 13-points molecule and leads to the well known block five diagonal matrix. Although quite a number of authors investigated numerical methods for the solution of this system of linear equations [1-4] the results are not satisfactory for nonlinear problems, i.e. the Navier-Stokes equations [5]. Some of the numerical problems arise from the discretization of the third and fourth derivatives, c3~s (xi,
1
yy) = ~ (~st_ 2,j - 4 ~ i - l , j "-4-6~1i, j
h2~ 6 - 41/1i+ l,j + I/#i+ 2,7) - ~ O x ~ l ( x i , Yy) + O ( h 4 )
which has to be changed into a non symmetric discretization in the vicinity of the boundaries using an approximation of the boundary conditions, which usually leads to unsatisfactory approximations of the normal derivative near the boundaries [6]. On the other hand, solvers designed for the block five diagonal coefficient matrix can not include the nonlinear terms without losing their efficiency. Explicit treatment of the nonlinearity, e.g. in a
A reliable solver for nonlinear biharmonic equations
407
biharmonic driver scheme [5] leads to convergence problems. Therefore the standard approximation of the biharmonic operator yields a coefficient matrix and hence a numerical algorithm which is not suitable for nonlinear problems. In order to overcome these problems we propose mixed approximations using the first partial derivatives ~,x and qJY as additional unknowns. This yields, e.g. a second order three point approximation 12 6 ~ h2 O4x~l(Xi, Yj) = -~ ($i-- ,.j -- 2~li, j + ~1i+ 1,j) + ~x (I]/i + l,y -- l//x- l,j) -- " ~ 06x~l(Xi, Yj) + O(h ~) of the fourth derivative with respect to x. These approximations share some properties with the well known Hermitian approximation used, e.g. in fluid mechanics by Schkalle et al. [7], and Thiele and Wagner [8] but they avoid the difficulties arising from the use of the mixed second derivative. In order to accomplish the linear system with a sufficient number of equations for the unknowns ~ and ~Y, respectively, we use the trivial identities 'i +
I~i+l,j--~li-l,j
=
fx
I
i = 2(I)N--
i - I
f Yj + 1
~¢i,j + l
--
~fi, j --1 :
((~y ~k ) (X,, y ) d Yj
1,
(rxO)(x, yj)dx, j=
I(I)M
i=
I(I)N,
dy,
- I
j = 2(1)M - I
and approximate the integrals by Simpson's rule to achieve the following fourth order approximation: i = 2 ( I ) N - I,
Pi,(')j(~l).__~li, Y 3 (~kij_, - tpi, j+,), -- Yj_ 1 + 4~b{j. + ~,,j+, + ,-ny
(3)
(7)
I(1)M
~=l(l)N,
3
x
, j=
Pi, j(~k ):=~kX-l,j + 40~,j + O i+l,j2l--hxx(~¢i-l, j --
(8)
¢ i + l,j), j=
2(I)M - 1
This leads to a complete approximation of the biharmonic operator which is of second order with three unknowns for each grid point. The problem mentioned above of the inaccurate approximation of the boundary conditions is avoided by this approximation. A numerical procedure designed to solve the system of linear equations efficiently will be outlined in the subsequent chapters. It will prove to be much more flexible than the standard solvers for the block five diagonal system; the superiority will be demonstrated with a driven cavity calculation. 4. SOLUTION OF THE LINEAR SYSTEM In this section we describe an efficient solution procedure for linear systems arising from mixed approximations. We assume that the system of linear equation can essentially be described by three equations, noted by p~l), p(2) and pO). p(2) reflects the approximation of the biharmonic operator n!2!,= 1~ I,J
1
~
t.j,i,
± j,i,j ,t,~
,.ij ,I,y
i = 2(I)N--
1,
ak, ltlti+k,j+ l T Uk, IW i + k , j + l -~ t'k, IW i + k , j + l '
k , l ffi - 1
y ffi 2 ( I ) M - 1
and the first boundary condition in eqn (6) gives i= 2(I)N-
I,
p!Z! t,J = g(xi, yj), j = 2(1)M - 1
We introduce a more suitable matrix formulation using the abbreviations
~,i,=(O,., ..... x 0 f , = (q,,.',,,
~,i,,,,)T, ~O,.M), • , 0 , ,y M)
T,
-]
)
i = I(1)N.
(9)
408
ANDREAS LIPPKE a n d HELMUT WAGNER
The coefficients o f 0~.j, 0 x~.j and ~,{j from eqn (9) form the matrices A~.k, A 7.k and A y~.,, respectively, 60,k
A~'=
0
0
k,-I
i,2 ak,o
0
"..
"..
"..
0
•
0
a i'M-I k,-1
i,M-1 ak,O
i,M-I all
0
6o,k
0
i,2 ak, l
0
......
"
0
0
i,2 b,._l
i,2 bk,o
i,2 bk.l
0
0
"'.
"'.
"'.
0
•
0
bi, M - [
i,M- 1 bk,o
bk,i , Mt - I
0
0
0
A~,=
......
ai,2
. . . . . .
k,-1
0
......
0
0
ci,2 k,-I
i,2 Ck,O
0
"'.
"'.
"'.
•
0
C i,M- I k,- I
i,M- I Ck,O
0
. . . . . .
0
i,2 Ck, l
0
......
"
0 i,M Ck, l
0
I
0
where 60,k denotes the K r o n n e c k e r symbol• The additional eqns (7) and (8) are accomplished with the discretized boundary condition from eqn (6):
p~)(qJ),=~., ~ky, + ~i.2 ~OY2+ ai.l ~Oi., + ag.2~J,.2 = r~,
i = I(1)N,
pi,M(~l)(I) '----~i,M-I~Ii,
i
(3)
1--~
P i , j (I]j ) - -
x
"JF~i,M~CYM'JI-~Ii, M_I~Ii,M_I-}-ai, M~li,M=r~, ,
,
I,j~t1 l,j + ~2 jl]J~ j + a l , j l f l l , j -1- a2,j~J2,j -~x
x
~
e
p~?, (~,),=Gu_ ,.j~O u _ ,,j + ~u.jqJ u.j + aN_ t.flJN-,.j + aN, jl]iN, j = r j , We introduce
C,,v =
~,J
~t,2
0
...
0
1
4
1
0
:
0
".
"'.
"'.
0
:
0
1
4
1
~;,M-I
~,U
0
C,,e =
......
a;,t
a;,2
0
...
0
3~by
0
- 3/hy
0
:
0
"..
"..
"'.
0
:
0
3/hy
0
- 3/hy
a~,M- I
?t;,M
0
......
I(1)N,
(lO)
j = l(l)M,
r 7,
~
=
j = l(1)M.
(ll)
A reliable solver for nonlinear biharmonic equations and Et,0==diag(~t,i,
•
•
,
409
al,u),
E [ 0 , = d i a g ( E u , • • , ~'~,M), Et.i.'=diag(t~2.1,. E~' 1-'=diag(/72,1,. EN,0==diag(a~¢,l,
a2,u), •, ~'~,M), .,
..,
a~,M),
E~v,0,=diag(bN, l, .., #,,,M), E~c,-l:=diag(aN-1,1
.....
aN-
I.M),
ETv,_~==diag(EN_ u . . . . . EN_ Lu), 3
E_t==-E,:=ffxxlM,
Eo==0,
Ex_fi=E~==IM,
E~:=4lu.
which serve to change eqns (7), (8), (10) and (11) into an appropriate matrix formulation. Rewriting the right hand sides o f the b o u n d a r y conditions eqn (10), eqn (11) yields
g~==(rW(yl ) , . . . , rW(yu)) r, gxN==(re(yl ) . . . . . re(yu)) r, gY'.=(r'(x,), 0 . . . . . O, r"(xi)) r,
i = I(1)N.
If we d r o p the suffix i for the m o m e n t and define the corrections d-'=tp - co(°)'
dX==~x _ ~ o ) , dy..=$y - ~ky(0). the complete system o f 3 M N linear equations written in terms o f correction reads d I = 0~ 1
x x Y Y Ai,,di+k + Ai,kdi+k + Ai,kdi+k
= r~ 2),
i = 2(1)N -
1,
k= -1
dN=0, (12)
Ei.odl + E u dE + E~,od~ + E~,I d~ = O, 1
x x _ ~., Ekdi+k + Ekdi+k -- "i-(3), i = 2(1)N - 1,
k=-I
E~c-ldtc-l+ENdN+ E ~N-t d ~N-t + E } d ~ = O, Ci, ldYl
+
Ci.2di,2 = 0 ,
i = I(1)N.
with the right hand sides r(2) i,i = 0 r(2)
~.3 = g i . j
- - n (2) h b
1-" i.j x'r
(°)'tJ,
j--2(1)M-1
r(2) i , M ~_ 0
i = I(1)N.
r O) -- 0 i,I --
(3) =.
ri, j
r (3) i,M __ -
-
,., (3) [ d , (0)"~
- - F i , j~,~ ,
0
j,
j = 2(1)M -- 1
(13)
410
ANDREAS L1PPKE and HELMUT WAGNER
In order to obtain the special structure of eqn (13) we use an initial guess ~(0) and ~k~(°), which satisfies the boundary conditions. The initial values ~y(0) are then calculated from equation
Ci,1 ~//y(O)= g ~ -
Ci,2~/(0), i = I(1)N.
This choice of an initial solution yields a great reduction in numerical costs due to the zero right hand side in eqns 02). To be more precise, we write the system in block tridiagonal form: ~,2
~,3
0
...
0
~,,
~,2
~,3
0
:
0
"'.
"'.
"'.
0
:
0
TN_,:
TN-1,2
T~_,,3
TN, 1
TN, 2
0
......
RI
Ol
=
:
•
I i
R.
D N ..,
using the block tridiagonal matrices
TI,I :=
~,1 : =
Ti,3'=
Cl, l
CI, 2
0
o
L,
o
0
El,o
E~,o
T|,2:=
0
0
0
hYi_l
Ai , l
A xi,
0
E_I
E ~- I
0
0
0
Ti,2'=
I
x
AYl
Ai. I
Ai, l
0
El
E~
CMA
Cu,2
0
o
G
o
0
EM.O
E,~,o
r~ 2)
and
D:=
TM, I:=
Ei,j E~j
Cia
Ci,2
0
A~ilo
Aio,
A xi,O
0
Eo
E~
0
0
0
0
0
0
EM , - l
E~M , - I
0 r~ 3)
4 d~
The standard LU decomposition of the coefficient matrix yields
L,=
Bl
0
0
...
0
T2a
B2
0
0
:
0
"'.
"'.
"'.
0
:
0
T,N_I,I
BN_I
0
TNj
B.
0
......
J
0
and the vectors
Ril =
0
°°° 1
0
0
TM.:=
0
(14)
411
A reliable solver for nonlinear biharmonic equations
and
r
z
3M
0
u:=
0
~ 0
w,
0
...
0
13&f
w2
0
i
**.
.a.
.a.
0
0
0
...
...
’
L&f WM.-,
0
L4
,
The non diagonal matrices are given by the following scheme [see eqn (14)]
4 = T1.2 i = 2(l)?/
Wi_,=B,?J-,,,, Bi= T,,z- Ti,, WI-15
(15)
i = 2( 1)N.
with the matrices
In order to solve eqn (15) we have to decompose the matrices Bi into
with the components e:,, = 41.1C,’ Q:,, = q;,2 -
Q:J Ci.2
Q :,3= &
(16)
Q:,, = d,, C,’ Qi,2 = CSi.2- Qi.1Ci,z)(Qi,2)-’ Q:., = cr:,>- QLQ:J Rewriting eqn (15) we obtain 0
A{l 0
0
0
(17)
A;l
Ai,l
E,
E;
and a simple calculation using Vi Wi = (LB)-‘Ti,3,
1
-Q:,&
0
0
Ai,,
A;,
E, - Q:,A,,
ET- Q:,4,
]
yields the solution of eqn (17) which very simply is due to the special structure of Ti,3.
412
ANDREAS LIPPKE a n d H E L M U T W A G N E R
We find W3,1 = -(Q~,3)-'Q~,zA! I, W 2 1,
i -1 =(Q2,2)
y i (A,,i-Qz3W3,,)
W3,z = (Q~,3)-I(E2 - Q~,2Ai,, ) W2.2 = (Q~2,2)-2(Ai, t - Q~,3 w3,2) __
i
--2
__
i
--2
W3.3-(Q3,3)
x
i
x
( E 2 - Q3,2A,,2) x
i
Wz3 - (Q2,2) (Ai, l -- Q2,3 W3,3) W2,1 =
- - C i~21 El, 2 W2,2
WL2 = - C i31 Ci.2 W2.2
W2,3 = - G 3 1 G.2 W2.3.
(19)
Using this decomposition eqn (14) could be solved using the additional vectors Zi B2 Z2 = R2,
i = 2(1)N
BiZi=R,-T~,2Z~_~,
(20)
DN= ZN, D,=Z~-I,V~Di+2
i=N-
1(-1)1
(21)
However, rearranging the backward sweep eqn (21) yields
Di= Z i - -
Bi-ITi,3Di+l
,¢~
Bi(Z i -- D,) = T,.3D,+ 2
(22)
and hence ON=
Bi V =
ZN~
T/,3Di+
i = N -- 1(-- 1)1
i,
Oi = Z i - V
(23)
with one additional vector V. By this modification we avoid the storage of large matrices W~. An estimate for the numerical complexity of the algorithm is obtained by calculation of the significant number of operations (multiplications and divisions) necessary to solve the linear system. Let To denote the number of operations to decompose a three band matrix and Ts the number of operations to solve a linear system with a decomposed three band matrix. In the case of non sparse systems we use An and As, respectively. The decompositions using eqn (16) and the solution step for eqn (15) with the matrices B~, i = I(1)N require Din= T n + 2 A o + 2 M T s + M A s + M 3 + 6M 2 and S~ = 33M z operations, respectively. The evaluations of the matrices W~using eqns (18), eqn (19) yield Diw = 6 M A s + 3 M T s + 3 M 3 + 18M z operations which add up to the "cost" of a complete LU decomposition of N
DT=
i=l
N
Z
i=2
+
operations. The complexity of the solution step is then S r = (2N - 1)(2(M 2 + A s ) + Ts).
Using the well known "costs" To = 2M, Ts = 3M, An = (M 3 - M ) / 3 and A s = M 2 yields Dr - 11M3N + 7 0 M Z N
(24)
Sr - 12MZN.
(25)
and
A reliable solver for nonlinear biharmonicequations
413
5. A NONLINEAR EXAMPLE Numerical methods for the Navier-Stokes equations are frequently tested on driven cavity flows. Until recently the available computational results were not very accurate, even for moderate Reynolds numbers. Therefore it seems worthwhile to test our numerical procedure with a driven cavity calculation and compare the results with those of Schreiber and Keller [9]. Thus we consider the plane steady laminar flow of an incompressible viscous fluid. In terms of the stream function the Navier-Stokes equations become
A25 + Re(O¢qJO, A ¢ - O,~bO¢A$) = 0.
(26)
The details are omitted here, we refer to the introductory book of Roache [10] or the article of Cebeci et al. [11] for a more detailed description of the various formulations and their numerical treatment. We just want to point out that the velocity components are represented in terms of the stream function ~ by
u(x,y)=~3y~(X,y)
and
v(x,y)=-~x~b(x,y
).
Therefore the velocities are easily and very accurately implemented in our numerical procedure. Furthermore, we can achieve a fourth order accurate difference approximation of the convection terms of the Navier-Stokes equation which leads to very accurate results even at high Reynolds numbers, e.g. we use 7.5
1.5
h~ (~3~b),j-I-
¢
O(h~)
for the discretization of the third derivative or
1 ;~] (~,,,j- ~,~,,,) - ~ (4~,~,j + 2~, L) + o(h~) 6
for the second derivative, Let D2, ,,,, D[j and D~j denote the difference formulas for the operators A2, a,A and a~A, respectively, at the point E~,t/j, then the local discretization errors zq.,=D?,,j
,,J"l'~" -
A25 (~l, t/j),
z ~c".,=D. , , j ,,J,'b -- A~3,$ (~i, t/j),
Tct--D¢ i , j . - - i j , - ,b - A a ¢ ¢ ( ¢ ~ , t/j), are easily calculated using eqn (26):
D~jd/ + Re[(O¢d/ ),,jD[y~k - (O, ql ),,jS,¢.j~b] = z ° j + (O¢$ )i,jRezC5 - (O,~b ),,jRezC~ =:%,j. If c denotes the absolute maximum of the velocity components in the whole discretization area c,=max {1(0¢~,),jI, I(d,¢),J I}, i,j
we obtain the following upper bound of the local discretization error
I .1
1 ,51 + Re c(lffTl+ Iff l)
(27)
with
= O(max(h¢, h,) 2) [%C]l + Iz,,c}l = O(max(h¢, hn)4) Equation (27) demonstrates the importance of an accurate approximation of the convective terms in the Navier-Stokes equations due to the magnification of the error by the Reynolds number. Our procedure uses fourth order difference molecules for all convective terms which leads to the desired very accurate results especially for high Reynolds numbers.
414
ANDREA$ LIPPKE a n d HELMUT WAGNER
In order to reduce the local discretization errors we used a coordinate transformation. It is well known that coordinate transformations are introduced to map an irregularly shaped domain onto a suitable mathematical domain to avoid problems related to non equispaced difference approximations. Many methods are available for this task, e.g. Thompson et al. [12] and Warsi [13, 14]. On the other side such transformations can successfully be used to redistribute the coordinate lines with the aim of a better approximation and hence smaller local errors in the boundary layers [6, 15-17]. We use a simple orthogonal transformation: { ( x ) = axX + (1 -- ~x)X3(3 -- 2x),
r/(y) = % y + (1 - ~y)y3(3 - 2y),
(ctx,%) ff (0, 11 x (0, 11.
This leads to more complex differential equation, for example, the transformed fourth derivative reads
040
=
1
:
(o~) ' ° 4 0 - ~ ( 6o x a~~){
,~3,h+ ~"
(a:¢): 15 ( ~
~(axe) -
4
0~~ ) (?~0 2 2
3
+ 10ax~ax~
(axe)6
2
,.(axe)
')(a~
3
4
"
ax~ ~0 '
7 ( O - - ~ J x~'
But these efforts will reduce the error significantly and are hence worthwhile. We employ a Newton-Chord method directly on the nonlinear differential equations to achieve a sequence of linear problems for the correction d. This leads to the following iterative scheme (k = 0(1). • .) A2d(~,l) + Re((3,,~,(k'°)t~yAd(k'° -- dy ~k(k'°)8x A d (ka)
+ &xd(k't)t~yh~) (k,O) - ~y d(k'O~ A~ (~,o))
t =0(1)K
= - AE~(ka) -- R e (t~ ~, (k'°ayA~, (k.t)
_ a, ~,(k,')~xa~,(k,')) 0(ka+ I),=~(k.0 + d(kJ)
0(~+,,o) = 0(~,,o In this context, the coefficient matrix is a very large sparse matrix, so a full Newton step will be many times more expansive than a chord step. To be more precise we recall eqns (24) and (25). The ratio of the number of operations for the decomposition to the number of operations for the solution of the linear system is Dr . llM3N Sr -
+ 70M2N_ M +6,
12M2N
hence, independent of N. On the other hand, the chord method is at best linearly convergent, while the Newton iterates converge quadratically to the solution if the differential equation has an isolated solution. Either method can fail to converge while the other method does converge. Therefore we determine the better method by estimating the correction per "cost" of computation. This idea was introduced by Schreiber and Keller [9] and we adapted this to our algorithm. We used a relative precision of 4 digits at every grid point for the following sequence of Reynolds numbers: R e = 0, 200, 500, 1000, 2000, 5000, 10,000. Our numerical experiments indicate that if the above mentioned Newton-Chord method is combined with a linear extrapolation step to predict the initial guess for the next Reynolds number, one Newton step followed by three or four chord steps lead to the final numerical solution.
A reliable solver for nonlinear biharmonic equations 6. N U M E R I C A L
415
RESULTS
In order to conform to the results of Schreiber and Keller [9] we will give the stream function and vorticity values for all vortices. We want to point out, that due to the compression of grid points in the corners of the cavity three corner vorticies are found in both lower corners, although only 121 x 121 grid points are used in the discretization. This is roughly half the number of unknowns used by Schreiber and Keller [9] and only a quarter of the unknowns used by Ghia et al. [18]. Thus we could reduce the required storage to ~ 400,000 words of central memory. A streamline plot at R e = 1000 (Fig. 1) and partial view of the underlying grid with some streamlines at Re = 15,000 (Fig. 2) are presented to show the high resolution in the boundary layers.
Fig. 1. Re = 1000.
Fig. 2. Re = 15,000, lower left corner.
416
ANDREAS LIPPKE a n d HELMUT WAGNER Table 1. Stream function values
A D G J M P
0.110000 0.060000 0.010000 0.000100 - 0.000001 -0.001000
B E H . K N Q
0.100000 0.040000 0.002000 0.000010 -0.000010 -0.002000
Fig. 3. Re = 5000.
Fig. 4. Re = 10,000.
C F I L O R
0.080000 0.020000 0.001000 0.000001 -0.000100 - 0.010000
A reliable solver for nonlinear biharmonic equations
417
Our results agree well with those of Ghia et al. [18] and Benjamin and Denny [19] at medium Reynolds numbers as well with those of Schreiber and Keller [9] in the whole range up to Re = 10,000. Nevertheless the tables include Reynolds numbers up to 30,000 which are not found in literature. On the subsequent pages streamline plots for various Reynolds numbers are given. They show the development of a central, almost circular vortex, with bottom vortices that do not shrink when Re grows. At Re = 2000 the third secondary vortex, near the upstream top corner, is present in the data. We choose identical streamline values (Table 1) in all plots to achieve comparable results. When the Reynolds number is increased to and beyond Re = 10,000 further vorticies appear in
Fig. 5. Re = 15,000.
Fig. 6. Re = 20,000.
ANDREAS LIPPKE a n d HELMUT WAGNER
418
Table 2 30x30
1 5 x 15
15x30
3 0 x 15
296 s 18.3 s
21.5 s 3.02 s
43.8 s 6.26 s
145 s 8.92 s 24.2 s 2.42 s
PC 80386/25 M H z
Newton Chord
SUN Sparcl
Newton Chord
50.4 s 5.66 s
3.85 s 0.883 s
7.93 s 1.65 s
Cray XMP
Newton Chord
17.9 s 0.983 s
1.31 s 0.155 s
2.63 s 0.320 s
8.82 s 0.475 s
N u m b e r of operations
eqn (24) eqn (25)
7.9~o + 5 4.1~o + 4
1.61o + 6 8.1no + 4
5.41o + 6 1.61o + 5
1.1 io + 7 3.21o + 5
Table 3
(a) Re
x
y
tp
co
0.1001 0.1188 0.1214 0.1206 0.1198 0.1193 0.1192 0.1192
0.325~o + 01 0.206~o + 01 0.192~o+ 01 0.187to+ 01 0.185to+ 01 0.183Lo+ 01 0.1821o+ 01 0.181~o+01
Primary vortex 0 1000 5000 10,000 15,000 20,000 25,000 30,000
0.500 0.466 0.489 0.489 0.489 0.477 0.477 0.477
0.233 0.433 0.466 0.477 0.477 0.489 0.489 0.477
(b) x
Re
y
¢
x
Secondary vortex lower left 0 1000 5000 10,000 15,000 20,000 25,000 30,000
0.034 0.133 0.184 0.203 0.222 0.231 0.241 0.241
0.966 0.891 0.926 0.939 0.950 0.956 0.966 0.966
x
y
-0.22751o -0.172510 -0.30921o -0.3260to -0.3088to -0.29071o -0.2710~o -0.25851o
y
~,
Tertiary vortex lower left -
05 02 02 02 02 02 02 02
0.003 0.009 0.020 0.061 0.068 0.061 0.055 0.055
0.997 0.994 0.980 0.932 0.912 0.905 0.891 0.891
+0.50801o +0,357910 +0.144010 + 0.1372to + 0.3385to +0.44131o-+0.5118~o + 0.56351o -
10 07 05 03 03 03 03 03
(c)
Re
~
x
Secondary vortex lower right 0 1000 5000 10,000 15,000 20,000 25,000 30,000
0.966 0.919 0.926 0.939 0.945 0.950 0.956 0.961
0.966 0.919 0.859 0.833 0.825 0.806 0.788 0.769
x
y
-0.227510 -0.23921o --0.1433~o -- 0.17591o -0.19081o - 0.19191o -0.1916~o -0.18921o
y
Tertiary vortex lower right ---
05 03 02 02 02 02 02 02
0.997 0.994 0.991 0.980 0.950 0.926 0.919 0.905
0.997 0.994 0.991 0.976 0.950 0.939 0.939 0.939
+0.508010 + 0.55131o -+ 0.7053to -+0.201 l ho -+0.46651o +0.14861o +0.2320~o +0.295910 -
10 08 07 05 04 03 03 03
(d)
Re
ff
x
Secondary vortex upper right 5000 10,000 15,000 20,000 25,000 30,000
0.939 0.932 0.926 0.919 0.919 0.919
y
0
Tertiary vortex upper right
0.087 0.087 0.087 0.087 0.087 0.081
- 0.14431o - 02 - 0.2587~o - 02 -0.31981o - 02 -0.35811o-02 -0.38651o - 02 -0.40961o - 02
0.988 0.976 0.971 0.971
y
t,O
x
0.166 0.175 0.175 0.184
+0.1215no+0.61621o +0.10951o +0.15171o-
04 04 03 03
(e)
Re
x
Quatery vortex lower left 10,000 15,000 20,000 25,000 30,000
0.003 0.006 0.006 0.009 0.012
0.994 0.994 0.994 0.991 0.988
- 0.2822to - 08 - 0 . 1 1 3 3 1 o - 07 -0.24151o - 07 -0.66811o-07 -0.217110-06
y
~b
Quatery vortex lower right 0.997 0.994 0.994 0.994
0.997 0.997 0.994 0.994
- 0 . 1 2 9 0 1 o - 08 -0.3473~o - 08 -0.718710-08 - 0 . 1 6 8 5 1 o - 07
the lower corners although they are not found in the streamline plots due to their small values. These vortices are nevertheless not spurious because the relative discretization error o f the calculation is smaller than 10 -3 at all grid points. At Re = 12,000 a small tertiary vortex appears within the secondary vortex near the upstream top corner (Fig. 5).
A reliable solver for nonlinear biharmonic equations
419
Fig. 7. Re = 30,000. A l l calculations were c a r r i e d o u t using the m a i n f r a m e a n d the P C version o f the p r o g r a m . W e used the C r a y - X M P at the Z I B Berlin a n d v a r i o u s p e r s o n a l c o m p u t e r s . T h e efficiency o f the a l g o r i t h m is d e m o n s t r a t e d b y a test p r o b l e m . W e c h o o s e a driven cavity with a small n u m b e r o f g r i d p o i n t s at a R e y n o l d s n u m b e r o f 10. It t o o k 4 i t e r a t i o n s (one N e w t o n - a n d three C h o r d - S t e p s ) to achieve a n a c c u r a c y o f 7 digits using o n l y the trivial initial a p p r o x i m a t i o n . All calculations are carried o u t with a n identical p r o g r a m written in C a n d all three m e n t i o n e d c o m p u t e r s are r u n n i n g a U N I X o p e r a t i n g system. H e n c e a full sequence d e s c r i b e d in T a b l e 2 ( R e = 0, 200, 500, 1000, 2000, 5000, 10,000) with a 30 x 30 grid takes a p p r o x i m a t e l y 120 C P U seconds o n the C r a y - X M P a n d 420 C P U seconds on a Sun S p a r c l c o m p u t e r (Table 3). REFERENCES 1. M. Amara and P. Destuynder, A numerical method for the biharmonic problem. Int. J. numerical Meth. Engng 17, 1515-1523 (1981). 2. L. Bauer and E. L. Reiss, Block five diagonal matrices and the fast numerical solution of the biharmonic equation. Math. Comput. 26, 311-326 (1972). 3. D. Braess and P. Peisker, On the numerical solution of the biharmonic equation and the role of squaring matrices for preconditioning. IMA J. numerical Analy. 6, 393-404 0986). 4. P. Bj6rstad, Fast numerical solution of the biharmonic Dirichlet problem on rectangles. S I A M J. numerical .4naly. 20, 59-71 (1983). 5. P. J. Roache and M. A. Ellis. The BID-method for the steady-state Navier-Stokes equations. Comput. Fluids 3, 305-320 (1975). 6. A. Lippke, Fortsetzungsmethoden und scknelle elliptische Liiser bei ziihen StdJmungen, Diss, TU Berlin (1985). 7. Y. D. Schkalle. F. Thiele and H. Wagner, A numerical method to solve the steady state Navier-Stokes equations for natural convection in enclosures. Recent Contributions to Fluid Mechanics (Edited by Haase), pp. 222-234. Springer, Berlin (1982). 8. F. Thiele, H. Wagner and J. Wepler, Solution of the Navier-Stokes equations using the finite difference method of hermitian type. 4. G A M M on Num. Meth. in Fluid Mech., Paris, pp. 326-336. Vieweg, Braunschweig (1981). 9. R. Schreiber and H. B. Keller, Driven cavity flows by efficient numerical techniques. J. comput. Phys. 49, 310-333 (1983). 10. P. J. Roache, Computational Fluid Dynamics. Hermosa (1982). l 1. T. Cebeci, R. S. Hirsh and H. B. Keller et al., Studies of numerical methods for the plane Navier-Stokes equations. Comput. Meth appL Mech. Engng 27, 13-44 (1981). 12. J. F. Thompson et aL, Numerical grid generation. Symp. on Numer. Gener. of Curv. Coord. Sys., Nashville, Tennessee (1982). 13. Z. U. A. Warsi, Numerical generation of orthogonal and non-orthogonal coordinates in two dimensional simply- and doubly-connected regions, yon Karman Institut for Fluid Dynamics, TN 151, VKI (1984). 14. Z. U. A. Warsi, Basic differential models for coordinate generation. Symp. on Numer Gener. of Curv. Coord. Sys., Nashville, Tennessee 0982).
420
ANDREASLWPKE and HELMUTWAGNER
15. Carcaillet et al., Generation of solution-adaptive computational grids using optimisation. Comput. Meth. appl. Mech. Engng 57, 279-295 (1986). 16. de Rivas, On the use of nonuniform grids in finite-difference equations. J. Comput. Phys. 10, 32-45 (1972). 17. H. Wagner, Ein Differenzenverfahren fiir die Navier-Stokes-Gleichungen in mehrfach zusammenhiingenden Gebieten. Diess. TU Berlin (1985). 18. U. Ghia, K. N. Ghia and C. T. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method. J. comput. Phys. 48, 387-411 (1982). 19. A. S. Benjamin and V. E. Denny, On the convergence of numerical solutions for 2-D flows in a cavity at high Re. J. Comput. Phys. 12, 348 (1973).
0045-7930/90 $3.00+ 0.00 Copyright© 1990PergamonPressple
A RELIABLE SOLVER FOR NONLINEAR BIHARMONIC
EQUATIONS ANDREAS LIPPKEI and HELMUTWAGNER2 ~Fachbereich Mathematik, Technische Universit/it Berlin and 2Dr Wagner Software GmbH, Berlin, Federal Republic of Germany (Received 29 January 1990; in revised form 22 May 1990)
Al~traet--Efficient and reliable numerical techniques of high-order accuracy are presented for solving problems for nonlinearperturbed biharmonicequations. The method is widelyapplicable,e.g. to problems of elasticity and in fluid mechanics. Here it is used to obtain accurate solutions for the driven cavity applying fourth order approximation on all convectiveterms. Solutions are calculated up to Reynolds number 30000.
I. INTRODUCTION We shall describe a numerical procedure for solving nonlinear biharmonic equations. This procedure compares very favorably with existing methods with respect to flexibility, accuracy and reliability, and should be a very efficient method in elasticity or in fluid mechanics. Briefly, the numerical method is this: the fourth order partial differential equation is discretized using central differences with three unknowns per grid point on an orthogonal net. This yields an approximation with truncation error expansion proceeding in powers of the mesh-width squared. The difference molecule involves only three grid lines in each direction for the three unknowns, the function and both the first derivatives. All but the fourth derivatives can even be discretized with fourth order error expansions. The nonlinear difference equations are solved by a sequence of Newton and Chord iterations using a modified LU decomposition. This special decomposition reduces the time and storage requirements of the procedure and allows calculations on a personal-computer. Numerical methods for nonlinear biharmonic equations can be tested on plane incompressible two dimensional flows, especially cavity flows. We applied the method to a plane cavity flow using fourth order approximations for the convective terms of the Navier-Stokes equations and a transformation to concentrate the grid lines in the boundary layers. This yields a highly accurate approximation even at high Reynolds numbers. 2. FORMULATION OF THE PROBLEM Let us consider the Dirichlet problem for the nonlinear biharmonic equation in an open set G c R 2 with a piecewise Lipschitz continuous boundary ~G. A2~k(~, r / ) + L[~b(~, r/)] = 0 ,
(~, r/) E G,
(1)
~k(~, ~/) = g,,
(¢, ~/) ~ dG,
(2)
~bn(~, r/) = g2,
(~, r/) ~ dG.
(3)
Here ft, denotes the normal derivatives of ~k with respect to the exterior normal. Let us assume that the nonlinearity L has a continuous and bounded Fr6chet derivative. In linear elasticity, ~b(~, ~7) can represent the Airy stress function or, as in the theory of thin plates, the vertical displacement due to an external force. In fluid mechanics, the equations describe the stream function formulation of the stationary Navier-Stokes equations. For the sake of simplicity we assume throughout the paper that the domain G is mapped onto a rectangle [ax, bx] x lay, by] by a simple transformation (x, y)~-*[¢ (x), r/(y)] with ~, ~/e C4(G). This assumption is not essential for the principles of the proposed numerical procedure but it reduces the notational efforts without losing the significant properties. CArls/4--~
405
406
ANDREAS LIPPKE and HELMUT WAGNER
Using ~x"= ~ x '
and
Oy..=~yy,
we derive the transformed representation of the biharmonic operator, ~f~2:=a4.odx--a2.2~x~y+ 4 2 2 a0,4 ~4 +
3 v u ~ av.~,Cgxdy, v+,a=0
(4)
and rewrite eqn (1) using the coordinates x, y in the mathematical domain and the notation/S for the transformed nonlinearity L and achieve
/%27s[x(~), y(r#)] + £Ts[x(¢), y(r#)] = O, (x, y) ~ (ax, bx) x (ay, by).
(5)
The coefficients av.~are given by the derivatives of the transformation functions ~ and q, mapping the physical domain onto the mathematical one, e.g.
(Ox'V'
a4.0 : t ~ ) • We assume that the transformed boundary conditions, e.g. on the "west" boundary can be written as a . ~ ¢ ( a x , y ) = gW(y) x -
X.[aW(y)O~qs(a~ ' Y) + b ~ ( Y ) O z ¢ ( a x , y ) = r~(Y)
(6)
with a ~ ( y ) 2 + b ' ( y ) 2 4= O. This boundary condition includes second order boundary values and is slightly more general than the one listed in eqn (3). These assumptions include not only the plane Navier-Stokes equations in stream function formulation with outflow boundary conditions but also, e.g. the axially symmetric Navier-Stokes equations. 3. D I S C R E T I Z A T I O N
For the discretization of the biharmonic operator in eqn (5) we introduce an uniform mesh with mesh sizes hx,= bx - ax N-1 '
and
hy := by - ay M-I
and grid points xi,=a,<+(i - 1)hx, y/=ay+(j
i = I(1)N,
- 1)hy, j = I(1)M.
We use the notation ~s;,s for the unknowns of the discrete system and ~s(xi, Yi) for the values of the exact solution, respectively. The centered second-order approximation of the biharmonic operator on a rectangular grid involves 13 unknowns arranged in the 13-points molecule and leads to the well known block five diagonal matrix. Although quite a number of authors investigated numerical methods for the solution of this system of linear equations [1-4] the results are not satisfactory for nonlinear problems, i.e. the Navier-Stokes equations [5]. Some of the numerical problems arise from the discretization of the third and fourth derivatives, c3~s (xi,
1
yy) = ~ (~st_ 2,j - 4 ~ i - l , j "-4-6~1i, j
h2~ 6 - 41/1i+ l,j + I/#i+ 2,7) - ~ O x ~ l ( x i , Yy) + O ( h 4 )
which has to be changed into a non symmetric discretization in the vicinity of the boundaries using an approximation of the boundary conditions, which usually leads to unsatisfactory approximations of the normal derivative near the boundaries [6]. On the other hand, solvers designed for the block five diagonal coefficient matrix can not include the nonlinear terms without losing their efficiency. Explicit treatment of the nonlinearity, e.g. in a
A reliable solver for nonlinear biharmonic equations
407
biharmonic driver scheme [5] leads to convergence problems. Therefore the standard approximation of the biharmonic operator yields a coefficient matrix and hence a numerical algorithm which is not suitable for nonlinear problems. In order to overcome these problems we propose mixed approximations using the first partial derivatives ~,x and qJY as additional unknowns. This yields, e.g. a second order three point approximation 12 6 ~ h2 O4x~l(Xi, Yj) = -~ ($i-- ,.j -- 2~li, j + ~1i+ 1,j) + ~x (I]/i + l,y -- l//x- l,j) -- " ~ 06x~l(Xi, Yj) + O(h ~) of the fourth derivative with respect to x. These approximations share some properties with the well known Hermitian approximation used, e.g. in fluid mechanics by Schkalle et al. [7], and Thiele and Wagner [8] but they avoid the difficulties arising from the use of the mixed second derivative. In order to accomplish the linear system with a sufficient number of equations for the unknowns ~ and ~Y, respectively, we use the trivial identities 'i +
I~i+l,j--~li-l,j
=
fx
I
i = 2(I)N--
i - I
f Yj + 1
~¢i,j + l
--
~fi, j --1 :
((~y ~k ) (X,, y ) d Yj
1,
(rxO)(x, yj)dx, j=
I(I)M
i=
I(I)N,
dy,
- I
j = 2(1)M - I
and approximate the integrals by Simpson's rule to achieve the following fourth order approximation: i = 2 ( I ) N - I,
Pi,(')j(~l).__~li, Y 3 (~kij_, - tpi, j+,), -- Yj_ 1 + 4~b{j. + ~,,j+, + ,-ny
(3)
(7)
I(1)M
~=l(l)N,
3
x
, j=
Pi, j(~k ):=~kX-l,j + 40~,j + O i+l,j2l--hxx(~¢i-l, j --
(8)
¢ i + l,j), j=
2(I)M - 1
This leads to a complete approximation of the biharmonic operator which is of second order with three unknowns for each grid point. The problem mentioned above of the inaccurate approximation of the boundary conditions is avoided by this approximation. A numerical procedure designed to solve the system of linear equations efficiently will be outlined in the subsequent chapters. It will prove to be much more flexible than the standard solvers for the block five diagonal system; the superiority will be demonstrated with a driven cavity calculation. 4. SOLUTION OF THE LINEAR SYSTEM In this section we describe an efficient solution procedure for linear systems arising from mixed approximations. We assume that the system of linear equation can essentially be described by three equations, noted by p~l), p(2) and pO). p(2) reflects the approximation of the biharmonic operator n!2!,= 1~ I,J
1
~
t.j,i,
± j,i,j ,t,~
,.ij ,I,y
i = 2(I)N--
1,
ak, ltlti+k,j+ l T Uk, IW i + k , j + l -~ t'k, IW i + k , j + l '
k , l ffi - 1
y ffi 2 ( I ) M - 1
and the first boundary condition in eqn (6) gives i= 2(I)N-
I,
p!Z! t,J = g(xi, yj), j = 2(1)M - 1
We introduce a more suitable matrix formulation using the abbreviations
~,i,=(O,., ..... x 0 f , = (q,,.',,,
~,i,,,,)T, ~O,.M), • , 0 , ,y M)
T,
-]
)
i = I(1)N.
(9)
408
ANDREAS LIPPKE a n d HELMUT WAGNER
The coefficients o f 0~.j, 0 x~.j and ~,{j from eqn (9) form the matrices A~.k, A 7.k and A y~.,, respectively, 60,k
A~'=
0
0
k,-I
i,2 ak,o
0
"..
"..
"..
0
•
0
a i'M-I k,-1
i,M-1 ak,O
i,M-I all
0
6o,k
0
i,2 ak, l
0
......
"
0
0
i,2 b,._l
i,2 bk,o
i,2 bk.l
0
0
"'.
"'.
"'.
0
•
0
bi, M - [
i,M- 1 bk,o
bk,i , Mt - I
0
0
0
A~,=
......
ai,2
. . . . . .
k,-1
0
......
0
0
ci,2 k,-I
i,2 Ck,O
0
"'.
"'.
"'.
•
0
C i,M- I k,- I
i,M- I Ck,O
0
. . . . . .
0
i,2 Ck, l
0
......
"
0 i,M Ck, l
0
I
0
where 60,k denotes the K r o n n e c k e r symbol• The additional eqns (7) and (8) are accomplished with the discretized boundary condition from eqn (6):
p~)(qJ),=~., ~ky, + ~i.2 ~OY2+ ai.l ~Oi., + ag.2~J,.2 = r~,
i = I(1)N,
pi,M(~l)(I) '----~i,M-I~Ii,
i
(3)
1--~
P i , j (I]j ) - -
x
"JF~i,M~CYM'JI-~Ii, M_I~Ii,M_I-}-ai, M~li,M=r~, ,
,
I,j~t1 l,j + ~2 jl]J~ j + a l , j l f l l , j -1- a2,j~J2,j -~x
x
~
e
p~?, (~,),=Gu_ ,.j~O u _ ,,j + ~u.jqJ u.j + aN_ t.flJN-,.j + aN, jl]iN, j = r j , We introduce
C,,v =
~,J
~t,2
0
...
0
1
4
1
0
:
0
".
"'.
"'.
0
:
0
1
4
1
~;,M-I
~,U
0
C,,e =
......
a;,t
a;,2
0
...
0
3~by
0
- 3/hy
0
:
0
"..
"..
"'.
0
:
0
3/hy
0
- 3/hy
a~,M- I
?t;,M
0
......
I(1)N,
(lO)
j = l(l)M,
r 7,
~
=
j = l(1)M.
(ll)
A reliable solver for nonlinear biharmonic equations and Et,0==diag(~t,i,
•
•
,
409
al,u),
E [ 0 , = d i a g ( E u , • • , ~'~,M), Et.i.'=diag(t~2.1,. E~' 1-'=diag(/72,1,. EN,0==diag(a~¢,l,
a2,u), •, ~'~,M), .,
..,
a~,M),
E~v,0,=diag(bN, l, .., #,,,M), E~c,-l:=diag(aN-1,1
.....
aN-
I.M),
ETv,_~==diag(EN_ u . . . . . EN_ Lu), 3
E_t==-E,:=ffxxlM,
Eo==0,
Ex_fi=E~==IM,
E~:=4lu.
which serve to change eqns (7), (8), (10) and (11) into an appropriate matrix formulation. Rewriting the right hand sides o f the b o u n d a r y conditions eqn (10), eqn (11) yields
g~==(rW(yl ) , . . . , rW(yu)) r, gxN==(re(yl ) . . . . . re(yu)) r, gY'.=(r'(x,), 0 . . . . . O, r"(xi)) r,
i = I(1)N.
If we d r o p the suffix i for the m o m e n t and define the corrections d-'=tp - co(°)'
dX==~x _ ~ o ) , dy..=$y - ~ky(0). the complete system o f 3 M N linear equations written in terms o f correction reads d I = 0~ 1
x x Y Y Ai,,di+k + Ai,kdi+k + Ai,kdi+k
= r~ 2),
i = 2(1)N -
1,
k= -1
dN=0, (12)
Ei.odl + E u dE + E~,od~ + E~,I d~ = O, 1
x x _ ~., Ekdi+k + Ekdi+k -- "i-(3), i = 2(1)N - 1,
k=-I
E~c-ldtc-l+ENdN+ E ~N-t d ~N-t + E } d ~ = O, Ci, ldYl
+
Ci.2di,2 = 0 ,
i = I(1)N.
with the right hand sides r(2) i,i = 0 r(2)
~.3 = g i . j
- - n (2) h b
1-" i.j x'r
(°)'tJ,
j--2(1)M-1
r(2) i , M ~_ 0
i = I(1)N.
r O) -- 0 i,I --
(3) =.
ri, j
r (3) i,M __ -
-
,., (3) [ d , (0)"~
- - F i , j~,~ ,
0
j,
j = 2(1)M -- 1
(13)
410
ANDREAS L1PPKE and HELMUT WAGNER
In order to obtain the special structure of eqn (13) we use an initial guess ~(0) and ~k~(°), which satisfies the boundary conditions. The initial values ~y(0) are then calculated from equation
Ci,1 ~//y(O)= g ~ -
Ci,2~/(0), i = I(1)N.
This choice of an initial solution yields a great reduction in numerical costs due to the zero right hand side in eqns 02). To be more precise, we write the system in block tridiagonal form: ~,2
~,3
0
...
0
~,,
~,2
~,3
0
:
0
"'.
"'.
"'.
0
:
0
TN_,:
TN-1,2
T~_,,3
TN, 1
TN, 2
0
......
RI
Ol
=
:
•
I i
R.
D N ..,
using the block tridiagonal matrices
TI,I :=
~,1 : =
Ti,3'=
Cl, l
CI, 2
0
o
L,
o
0
El,o
E~,o
T|,2:=
0
0
0
hYi_l
Ai , l
A xi,
0
E_I
E ~- I
0
0
0
Ti,2'=
I
x
AYl
Ai. I
Ai, l
0
El
E~
CMA
Cu,2
0
o
G
o
0
EM.O
E,~,o
r~ 2)
and
D:=
TM, I:=
Ei,j E~j
Cia
Ci,2
0
A~ilo
Aio,
A xi,O
0
Eo
E~
0
0
0
0
0
0
EM , - l
E~M , - I
0 r~ 3)
4 d~
The standard LU decomposition of the coefficient matrix yields
L,=
Bl
0
0
...
0
T2a
B2
0
0
:
0
"'.
"'.
"'.
0
:
0
T,N_I,I
BN_I
0
TNj
B.
0
......
J
0
and the vectors
Ril =
0
°°° 1
0
0
TM.:=
0
(14)
411
A reliable solver for nonlinear biharmonic equations
and
r
z
3M
0
u:=
0
~ 0
w,
0
...
0
13&f
w2
0
i
**.
.a.
.a.
0
0
0
...
...
’
L&f WM.-,
0
L4
,
The non diagonal matrices are given by the following scheme [see eqn (14)]
4 = T1.2 i = 2(l)?/
Wi_,=B,?J-,,,, Bi= T,,z- Ti,, WI-15
(15)
i = 2( 1)N.
with the matrices
In order to solve eqn (15) we have to decompose the matrices Bi into
with the components e:,, = 41.1C,’ Q:,, = q;,2 -
Q:J Ci.2
Q :,3= &
(16)
Q:,, = d,, C,’ Qi,2 = CSi.2- Qi.1Ci,z)(Qi,2)-’ Q:., = cr:,>- QLQ:J Rewriting eqn (15) we obtain 0
A{l 0
0
0
(17)
A;l
Ai,l
E,
E;
and a simple calculation using Vi Wi = (LB)-‘Ti,3,
1
-Q:,&
0
0
Ai,,
A;,
E, - Q:,A,,
ET- Q:,4,
]
yields the solution of eqn (17) which very simply is due to the special structure of Ti,3.
412
ANDREAS LIPPKE a n d H E L M U T W A G N E R
We find W3,1 = -(Q~,3)-'Q~,zA! I, W 2 1,
i -1 =(Q2,2)
y i (A,,i-Qz3W3,,)
W3,z = (Q~,3)-I(E2 - Q~,2Ai,, ) W2.2 = (Q~2,2)-2(Ai, t - Q~,3 w3,2) __
i
--2
__
i
--2
W3.3-(Q3,3)
x
i
x
( E 2 - Q3,2A,,2) x
i
Wz3 - (Q2,2) (Ai, l -- Q2,3 W3,3) W2,1 =
- - C i~21 El, 2 W2,2
WL2 = - C i31 Ci.2 W2.2
W2,3 = - G 3 1 G.2 W2.3.
(19)
Using this decomposition eqn (14) could be solved using the additional vectors Zi B2 Z2 = R2,
i = 2(1)N
BiZi=R,-T~,2Z~_~,
(20)
DN= ZN, D,=Z~-I,V~Di+2
i=N-
1(-1)1
(21)
However, rearranging the backward sweep eqn (21) yields
Di= Z i - -
Bi-ITi,3Di+l
,¢~
Bi(Z i -- D,) = T,.3D,+ 2
(22)
and hence ON=
Bi V =
ZN~
T/,3Di+
i = N -- 1(-- 1)1
i,
Oi = Z i - V
(23)
with one additional vector V. By this modification we avoid the storage of large matrices W~. An estimate for the numerical complexity of the algorithm is obtained by calculation of the significant number of operations (multiplications and divisions) necessary to solve the linear system. Let To denote the number of operations to decompose a three band matrix and Ts the number of operations to solve a linear system with a decomposed three band matrix. In the case of non sparse systems we use An and As, respectively. The decompositions using eqn (16) and the solution step for eqn (15) with the matrices B~, i = I(1)N require Din= T n + 2 A o + 2 M T s + M A s + M 3 + 6M 2 and S~ = 33M z operations, respectively. The evaluations of the matrices W~using eqns (18), eqn (19) yield Diw = 6 M A s + 3 M T s + 3 M 3 + 18M z operations which add up to the "cost" of a complete LU decomposition of N
DT=
i=l
N
Z
i=2
+
operations. The complexity of the solution step is then S r = (2N - 1)(2(M 2 + A s ) + Ts).
Using the well known "costs" To = 2M, Ts = 3M, An = (M 3 - M ) / 3 and A s = M 2 yields Dr - 11M3N + 7 0 M Z N
(24)
Sr - 12MZN.
(25)
and
A reliable solver for nonlinear biharmonicequations
413
5. A NONLINEAR EXAMPLE Numerical methods for the Navier-Stokes equations are frequently tested on driven cavity flows. Until recently the available computational results were not very accurate, even for moderate Reynolds numbers. Therefore it seems worthwhile to test our numerical procedure with a driven cavity calculation and compare the results with those of Schreiber and Keller [9]. Thus we consider the plane steady laminar flow of an incompressible viscous fluid. In terms of the stream function the Navier-Stokes equations become
A25 + Re(O¢qJO, A ¢ - O,~bO¢A$) = 0.
(26)
The details are omitted here, we refer to the introductory book of Roache [10] or the article of Cebeci et al. [11] for a more detailed description of the various formulations and their numerical treatment. We just want to point out that the velocity components are represented in terms of the stream function ~ by
u(x,y)=~3y~(X,y)
and
v(x,y)=-~x~b(x,y
).
Therefore the velocities are easily and very accurately implemented in our numerical procedure. Furthermore, we can achieve a fourth order accurate difference approximation of the convection terms of the Navier-Stokes equation which leads to very accurate results even at high Reynolds numbers, e.g. we use 7.5
1.5
h~ (~3~b),j-I-
¢
O(h~)
for the discretization of the third derivative or
1 ;~] (~,,,j- ~,~,,,) - ~ (4~,~,j + 2~, L) + o(h~) 6
for the second derivative, Let D2, ,,,, D[j and D~j denote the difference formulas for the operators A2, a,A and a~A, respectively, at the point E~,t/j, then the local discretization errors zq.,=D?,,j
,,J"l'~" -
A25 (~l, t/j),
z ~c".,=D. , , j ,,J,'b -- A~3,$ (~i, t/j),
Tct--D¢ i , j . - - i j , - ,b - A a ¢ ¢ ( ¢ ~ , t/j), are easily calculated using eqn (26):
D~jd/ + Re[(O¢d/ ),,jD[y~k - (O, ql ),,jS,¢.j~b] = z ° j + (O¢$ )i,jRezC5 - (O,~b ),,jRezC~ =:%,j. If c denotes the absolute maximum of the velocity components in the whole discretization area c,=max {1(0¢~,),jI, I(d,¢),J I}, i,j
we obtain the following upper bound of the local discretization error
I .1
1 ,51 + Re c(lffTl+ Iff l)
(27)
with
= O(max(h¢, h,) 2) [%C]l + Iz,,c}l = O(max(h¢, hn)4) Equation (27) demonstrates the importance of an accurate approximation of the convective terms in the Navier-Stokes equations due to the magnification of the error by the Reynolds number. Our procedure uses fourth order difference molecules for all convective terms which leads to the desired very accurate results especially for high Reynolds numbers.
414
ANDREA$ LIPPKE a n d HELMUT WAGNER
In order to reduce the local discretization errors we used a coordinate transformation. It is well known that coordinate transformations are introduced to map an irregularly shaped domain onto a suitable mathematical domain to avoid problems related to non equispaced difference approximations. Many methods are available for this task, e.g. Thompson et al. [12] and Warsi [13, 14]. On the other side such transformations can successfully be used to redistribute the coordinate lines with the aim of a better approximation and hence smaller local errors in the boundary layers [6, 15-17]. We use a simple orthogonal transformation: { ( x ) = axX + (1 -- ~x)X3(3 -- 2x),
r/(y) = % y + (1 - ~y)y3(3 - 2y),
(ctx,%) ff (0, 11 x (0, 11.
This leads to more complex differential equation, for example, the transformed fourth derivative reads
040
=
1
:
(o~) ' ° 4 0 - ~ ( 6o x a~~){
,~3,h+ ~"
(a:¢): 15 ( ~
~(axe) -
4
0~~ ) (?~0 2 2
3
+ 10ax~ax~
(axe)6
2
,.(axe)
')(a~
3
4
"
ax~ ~0 '
7 ( O - - ~ J x~'
But these efforts will reduce the error significantly and are hence worthwhile. We employ a Newton-Chord method directly on the nonlinear differential equations to achieve a sequence of linear problems for the correction d. This leads to the following iterative scheme (k = 0(1). • .) A2d(~,l) + Re((3,,~,(k'°)t~yAd(k'° -- dy ~k(k'°)8x A d (ka)
+ &xd(k't)t~yh~) (k,O) - ~y d(k'O~ A~ (~,o))
t =0(1)K
= - AE~(ka) -- R e (t~ ~, (k'°ayA~, (k.t)
_ a, ~,(k,')~xa~,(k,')) 0(ka+ I),=~(k.0 + d(kJ)
0(~+,,o) = 0(~,,o In this context, the coefficient matrix is a very large sparse matrix, so a full Newton step will be many times more expansive than a chord step. To be more precise we recall eqns (24) and (25). The ratio of the number of operations for the decomposition to the number of operations for the solution of the linear system is Dr . llM3N Sr -
+ 70M2N_ M +6,
12M2N
hence, independent of N. On the other hand, the chord method is at best linearly convergent, while the Newton iterates converge quadratically to the solution if the differential equation has an isolated solution. Either method can fail to converge while the other method does converge. Therefore we determine the better method by estimating the correction per "cost" of computation. This idea was introduced by Schreiber and Keller [9] and we adapted this to our algorithm. We used a relative precision of 4 digits at every grid point for the following sequence of Reynolds numbers: R e = 0, 200, 500, 1000, 2000, 5000, 10,000. Our numerical experiments indicate that if the above mentioned Newton-Chord method is combined with a linear extrapolation step to predict the initial guess for the next Reynolds number, one Newton step followed by three or four chord steps lead to the final numerical solution.
A reliable solver for nonlinear biharmonic equations 6. N U M E R I C A L
415
RESULTS
In order to conform to the results of Schreiber and Keller [9] we will give the stream function and vorticity values for all vortices. We want to point out, that due to the compression of grid points in the corners of the cavity three corner vorticies are found in both lower corners, although only 121 x 121 grid points are used in the discretization. This is roughly half the number of unknowns used by Schreiber and Keller [9] and only a quarter of the unknowns used by Ghia et al. [18]. Thus we could reduce the required storage to ~ 400,000 words of central memory. A streamline plot at R e = 1000 (Fig. 1) and partial view of the underlying grid with some streamlines at Re = 15,000 (Fig. 2) are presented to show the high resolution in the boundary layers.
Fig. 1. Re = 1000.
Fig. 2. Re = 15,000, lower left corner.
416
ANDREAS LIPPKE a n d HELMUT WAGNER Table 1. Stream function values
A D G J M P
0.110000 0.060000 0.010000 0.000100 - 0.000001 -0.001000
B E H . K N Q
0.100000 0.040000 0.002000 0.000010 -0.000010 -0.002000
Fig. 3. Re = 5000.
Fig. 4. Re = 10,000.
C F I L O R
0.080000 0.020000 0.001000 0.000001 -0.000100 - 0.010000
A reliable solver for nonlinear biharmonic equations
417
Our results agree well with those of Ghia et al. [18] and Benjamin and Denny [19] at medium Reynolds numbers as well with those of Schreiber and Keller [9] in the whole range up to Re = 10,000. Nevertheless the tables include Reynolds numbers up to 30,000 which are not found in literature. On the subsequent pages streamline plots for various Reynolds numbers are given. They show the development of a central, almost circular vortex, with bottom vortices that do not shrink when Re grows. At Re = 2000 the third secondary vortex, near the upstream top corner, is present in the data. We choose identical streamline values (Table 1) in all plots to achieve comparable results. When the Reynolds number is increased to and beyond Re = 10,000 further vorticies appear in
Fig. 5. Re = 15,000.
Fig. 6. Re = 20,000.
ANDREAS LIPPKE a n d HELMUT WAGNER
418
Table 2 30x30
1 5 x 15
15x30
3 0 x 15
296 s 18.3 s
21.5 s 3.02 s
43.8 s 6.26 s
145 s 8.92 s 24.2 s 2.42 s
PC 80386/25 M H z
Newton Chord
SUN Sparcl
Newton Chord
50.4 s 5.66 s
3.85 s 0.883 s
7.93 s 1.65 s
Cray XMP
Newton Chord
17.9 s 0.983 s
1.31 s 0.155 s
2.63 s 0.320 s
8.82 s 0.475 s
N u m b e r of operations
eqn (24) eqn (25)
7.9~o + 5 4.1~o + 4
1.61o + 6 8.1no + 4
5.41o + 6 1.61o + 5
1.1 io + 7 3.21o + 5
Table 3
(a) Re
x
y
tp
co
0.1001 0.1188 0.1214 0.1206 0.1198 0.1193 0.1192 0.1192
0.325~o + 01 0.206~o + 01 0.192~o+ 01 0.187to+ 01 0.185to+ 01 0.183Lo+ 01 0.1821o+ 01 0.181~o+01
Primary vortex 0 1000 5000 10,000 15,000 20,000 25,000 30,000
0.500 0.466 0.489 0.489 0.489 0.477 0.477 0.477
0.233 0.433 0.466 0.477 0.477 0.489 0.489 0.477
(b) x
Re
y
¢
x
Secondary vortex lower left 0 1000 5000 10,000 15,000 20,000 25,000 30,000
0.034 0.133 0.184 0.203 0.222 0.231 0.241 0.241
0.966 0.891 0.926 0.939 0.950 0.956 0.966 0.966
x
y
-0.22751o -0.172510 -0.30921o -0.3260to -0.3088to -0.29071o -0.2710~o -0.25851o
y
~,
Tertiary vortex lower left -
05 02 02 02 02 02 02 02
0.003 0.009 0.020 0.061 0.068 0.061 0.055 0.055
0.997 0.994 0.980 0.932 0.912 0.905 0.891 0.891
+0.50801o +0,357910 +0.144010 + 0.1372to + 0.3385to +0.44131o-+0.5118~o + 0.56351o -
10 07 05 03 03 03 03 03
(c)
Re
~
x
Secondary vortex lower right 0 1000 5000 10,000 15,000 20,000 25,000 30,000
0.966 0.919 0.926 0.939 0.945 0.950 0.956 0.961
0.966 0.919 0.859 0.833 0.825 0.806 0.788 0.769
x
y
-0.227510 -0.23921o --0.1433~o -- 0.17591o -0.19081o - 0.19191o -0.1916~o -0.18921o
y
Tertiary vortex lower right ---
05 03 02 02 02 02 02 02
0.997 0.994 0.991 0.980 0.950 0.926 0.919 0.905
0.997 0.994 0.991 0.976 0.950 0.939 0.939 0.939
+0.508010 + 0.55131o -+ 0.7053to -+0.201 l ho -+0.46651o +0.14861o +0.2320~o +0.295910 -
10 08 07 05 04 03 03 03
(d)
Re
ff
x
Secondary vortex upper right 5000 10,000 15,000 20,000 25,000 30,000
0.939 0.932 0.926 0.919 0.919 0.919
y
0
Tertiary vortex upper right
0.087 0.087 0.087 0.087 0.087 0.081
- 0.14431o - 02 - 0.2587~o - 02 -0.31981o - 02 -0.35811o-02 -0.38651o - 02 -0.40961o - 02
0.988 0.976 0.971 0.971
y
t,O
x
0.166 0.175 0.175 0.184
+0.1215no+0.61621o +0.10951o +0.15171o-
04 04 03 03
(e)
Re
x
Quatery vortex lower left 10,000 15,000 20,000 25,000 30,000
0.003 0.006 0.006 0.009 0.012
0.994 0.994 0.994 0.991 0.988
- 0.2822to - 08 - 0 . 1 1 3 3 1 o - 07 -0.24151o - 07 -0.66811o-07 -0.217110-06
y
~b
Quatery vortex lower right 0.997 0.994 0.994 0.994
0.997 0.997 0.994 0.994
- 0 . 1 2 9 0 1 o - 08 -0.3473~o - 08 -0.718710-08 - 0 . 1 6 8 5 1 o - 07
the lower corners although they are not found in the streamline plots due to their small values. These vortices are nevertheless not spurious because the relative discretization error o f the calculation is smaller than 10 -3 at all grid points. At Re = 12,000 a small tertiary vortex appears within the secondary vortex near the upstream top corner (Fig. 5).
A reliable solver for nonlinear biharmonic equations
419
Fig. 7. Re = 30,000. A l l calculations were c a r r i e d o u t using the m a i n f r a m e a n d the P C version o f the p r o g r a m . W e used the C r a y - X M P at the Z I B Berlin a n d v a r i o u s p e r s o n a l c o m p u t e r s . T h e efficiency o f the a l g o r i t h m is d e m o n s t r a t e d b y a test p r o b l e m . W e c h o o s e a driven cavity with a small n u m b e r o f g r i d p o i n t s at a R e y n o l d s n u m b e r o f 10. It t o o k 4 i t e r a t i o n s (one N e w t o n - a n d three C h o r d - S t e p s ) to achieve a n a c c u r a c y o f 7 digits using o n l y the trivial initial a p p r o x i m a t i o n . All calculations are carried o u t with a n identical p r o g r a m written in C a n d all three m e n t i o n e d c o m p u t e r s are r u n n i n g a U N I X o p e r a t i n g system. H e n c e a full sequence d e s c r i b e d in T a b l e 2 ( R e = 0, 200, 500, 1000, 2000, 5000, 10,000) with a 30 x 30 grid takes a p p r o x i m a t e l y 120 C P U seconds o n the C r a y - X M P a n d 420 C P U seconds on a Sun S p a r c l c o m p u t e r (Table 3). REFERENCES 1. M. Amara and P. Destuynder, A numerical method for the biharmonic problem. Int. J. numerical Meth. Engng 17, 1515-1523 (1981). 2. L. Bauer and E. L. Reiss, Block five diagonal matrices and the fast numerical solution of the biharmonic equation. Math. Comput. 26, 311-326 (1972). 3. D. Braess and P. Peisker, On the numerical solution of the biharmonic equation and the role of squaring matrices for preconditioning. IMA J. numerical Analy. 6, 393-404 0986). 4. P. Bj6rstad, Fast numerical solution of the biharmonic Dirichlet problem on rectangles. S I A M J. numerical .4naly. 20, 59-71 (1983). 5. P. J. Roache and M. A. Ellis. The BID-method for the steady-state Navier-Stokes equations. Comput. Fluids 3, 305-320 (1975). 6. A. Lippke, Fortsetzungsmethoden und scknelle elliptische Liiser bei ziihen StdJmungen, Diss, TU Berlin (1985). 7. Y. D. Schkalle. F. Thiele and H. Wagner, A numerical method to solve the steady state Navier-Stokes equations for natural convection in enclosures. Recent Contributions to Fluid Mechanics (Edited by Haase), pp. 222-234. Springer, Berlin (1982). 8. F. Thiele, H. Wagner and J. Wepler, Solution of the Navier-Stokes equations using the finite difference method of hermitian type. 4. G A M M on Num. Meth. in Fluid Mech., Paris, pp. 326-336. Vieweg, Braunschweig (1981). 9. R. Schreiber and H. B. Keller, Driven cavity flows by efficient numerical techniques. J. comput. Phys. 49, 310-333 (1983). 10. P. J. Roache, Computational Fluid Dynamics. Hermosa (1982). l 1. T. Cebeci, R. S. Hirsh and H. B. Keller et al., Studies of numerical methods for the plane Navier-Stokes equations. Comput. Meth appL Mech. Engng 27, 13-44 (1981). 12. J. F. Thompson et aL, Numerical grid generation. Symp. on Numer. Gener. of Curv. Coord. Sys., Nashville, Tennessee (1982). 13. Z. U. A. Warsi, Numerical generation of orthogonal and non-orthogonal coordinates in two dimensional simply- and doubly-connected regions, yon Karman Institut for Fluid Dynamics, TN 151, VKI (1984). 14. Z. U. A. Warsi, Basic differential models for coordinate generation. Symp. on Numer Gener. of Curv. Coord. Sys., Nashville, Tennessee 0982).
420
ANDREASLWPKE and HELMUTWAGNER
15. Carcaillet et al., Generation of solution-adaptive computational grids using optimisation. Comput. Meth. appl. Mech. Engng 57, 279-295 (1986). 16. de Rivas, On the use of nonuniform grids in finite-difference equations. J. Comput. Phys. 10, 32-45 (1972). 17. H. Wagner, Ein Differenzenverfahren fiir die Navier-Stokes-Gleichungen in mehrfach zusammenhiingenden Gebieten. Diess. TU Berlin (1985). 18. U. Ghia, K. N. Ghia and C. T. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method. J. comput. Phys. 48, 387-411 (1982). 19. A. S. Benjamin and V. E. Denny, On the convergence of numerical solutions for 2-D flows in a cavity at high Re. J. Comput. Phys. 12, 348 (1973).