- Email: [email protected]
A time-split finite element scheme for the PC
Software 23 (1995) 15-26 8 1995 Else&r Science Limited Printed in Great Britain. All rights reserved 0965-9978/95/%09.50
Advances
0965-9978(95)00022-4
ELSEVIER
in Engineering
A time-split finite element scheme for the PC Darrell W. Pepper Depurtment
of
Mechanical Engineering, University
of Nevada,
Las Vegas, Las Vegas, NV 891.54-4027, (IsA
(Received 14 February 1994; revised version received 27 April 1995; accepted 12 May 1995) A time-split finite element scheme is developed based on one-dimensional, linear basis functions. The transient form of the transport equation is discretized and a recursion equation obtained through assemblage over two adjacent elements. The resulting one-dimensional scheme can be readily applied to multi-dimensional problems, requiring only the solution of a simple tri-diagonal banded matrix (in each direction). Computer program listings are included for solving twodimensional advection and displaying a surface plot of a rotating cosine hill distribution on a PC. Key words: time-split, finite element, PC. recursion relation.
1. INTRODUCTION
mass matrix to be easily calculated. However, the use of mass lumping has a tendency to introduce some error in the transient analysis; for problems which converge to steady state, the steady state results are essentially as accurate as those obtained without lumping. Another technique which yields reasonably good results uses one-point (Gauss) quadrature in lieu of 2 x 2 quadrature for four node quadrilaterals (or 2 x 2 x 2 for eight node hexahedrals). Utilizing mass lumping in an explicit time marching scheme with onepoint quadrature, Gresho et al.’ showed that a modified finite element method could produce results comparable in accuracy to finite differences. While not quite the speed of the finite difference method, the technique was significantly faster than the more conventional finite element method. The modifications made by Gresho et al.’ to the finite element method have the effect of altering some of the terms and matrix operations to behave similarly to the finite volume method.2 Development of a finite element-finite volume method was conducted by Prakash & Patankar’ in an effort to more formally combine the finite element method with the finite volume technique; two-dimensional elements were used. While this technique exhibited the advantages of geometrical flexibility and computational speed, there was some loss of accuracy. The concept of time-splitting can be applied to the finite element method to yield a finite element algorithm with the advantages of the finite difference technique. The method of time-splitting, or method of fractional steps4 creates a series of one-dimensional relations which can be multiplied (tensor product) to yield solutions to multi-dimensional problems. The
The advantages of the finite element method over the finite difference method are well known: geometrical flexibility, computational accuracy and stability, ease in altering boundary conditions without requiring changes in the program, and ability to easily use various order interpolation functions. However, the one disadvantage which clearly becomes evident is the amount of computational speed and storage required for problems involving many nodes. In an implicit difference method, the use of one-sided and central differencing techniques to approximate derivative terms produces a set of symmetric banded matrices consisting of usually three to five diagonal terms. The Thomas algorithm is a simple procedure which quickly solves tri-diagonal matrices without involving the zero value off-diagonal terms. If an explicit time scheme is used, the mass matrix associated with the time derivative term in the finite difference technique consists of only one value; hence, the matrix equivalent equation can be solved without the need for a matrix solver. The mass matrix of the finite element method, on the other hand, is symmetric and banded according to the elemental node numbering, regardless of the type of time differencing employed. Over the last few years, efforts have been made to alter the finite element method in order to reduce computational storage requirements and computing times, without significantly decreasing accuracy. The use of mass lumping is one such method whereby the off-diagonal terms are lumped into a single, main diagonal. This procedure permits the inverse of the 15
16
D. W. Pepper
technique has the advantages of simplicity, i.e. onedimensional element formulation, and computation speed, i.e. tri-diagonal matrix solution common to finite difference methods. The time-split algorithm has been used to solve multi-dimensional problems involving two- and three-dimensional scalar transport and fluid dynamics, and has been shown to produce accurate results at computational speeds equivalent to finite difference methods.5’6
2. THE MULTI-DIMENSIONAL ALGORITHM
with a4
+
L
*
%nx+k
a4
^
a4
-n,,+k,~nz+q=O yy ay
n
$(X,Y, z, 1= 0) = $0
= -WA + by1+ M + kl + FYI+ M&V + [M]At(O{F}“+l
TIME-SPLIT
The three-dimensional transport equation for a scalar quantity (e.g. temperature) can be written in the form
4
where [vJ, [r+,], [Q] denote the advection matrices for the u, w, and w velocity components, and {Fi} is the load vector consisting of known values (FXFJZ). Approximating the time derivative term using the 6 method, eqn (5) becomes
(2)
(3)
where eqns (2) and (3) represent the generalized boundary condition constraints and initial values for eqn (1). The finite element weak statement of eqn (1) is obtained as
+ (1 - e){F}“)
(6)
where 0 5 13I 1 (if B = l/2, one obtains CrankNicolson time averaging). The one-dimensional algorithm yields for the left hand side of the matrix equation, ([Ml + BAt( [w,] + [k,])), which is only tri-diagonal using linear one-dimensional elements, or penta-diagonal using quadratic elements. In contrast, the left hand side (matrix) of eqn (6) is ‘very large’, since it contains the three-dimensional coupling of all derivatives. Formulation of this matrix and subsequent solution of eqn (6) can quickly saturate the available memory of a computer using even a coarse mesh and the conventional finite element procedure. However, a number of methods exist which can be used to more efficiently solve eqn (6), such as Gauss-Seidel iteration, alternatingdirection implicit (ADI), approximate factorization (AF), and operator-splitting. In each instance, the multi-dimensional derivative operators are approximated by a sequence of one-dimensional operations. The tensor matrix product of two vectors produces a matrix of higher rank (or order), e.g. the multiplication of two 2 x 2 matrices produces a 4 x 4 matrix. The product of two vectors, e.g. two basis functions, produces a higher order vector. The basic building block is the one-dimensional mass matrix [Ml. Refining notation for n-dimensions, the onedimensional (x-component) form is written as (7)
=
qNi dS (4 If s Note that eqn (4) is ‘functionally’ identical to a onedimensional finite element relation; the three-dimensional form simply generalizes Ni and the advection and diffusion terms possess scalar components parallel to each component of the x, y, z coordinate system. Using subscripts to denote the x, y, z terms, the matrix equivalent expression for eqn (4) is f ff
QNi dx dy dz -
4 2 =-
IN{ $}+(Lvxl + Lwyl + [w.zl){4i1 + ([kxl+ [kyl+ [kI){A) = {Fi:i)
where Ax E l(‘), which is the element length in the xdirection. Notice that the element order and rank are both equal to 2. The two-dimensional mass matrix is constructed using the tensor product (denoted by @) of the two one-dimensional mass matrices, [M,] and [My], i.e.
(5)
AxAy 36
1 2
2 4 2 1 1 2 4 2 2
1 2 4
(8)
A time-split jinite element scheme for the PC which is analogous to the mass matrix relation for the two-dimensional linear quadrilateral element. The threedimensional mass matrix is similarly formed by multi-
plying ew (8) by NJ. The true value of the tensor product is accurately and efficiently approximates the left eqn (6). Every matrix given in eqn (6) is computable using [M,], [MY], [M,], the tensor product, e.g.
that it side of directly matrix
17
into (12) the left side of eqn (12) becomes
[A,I @M,pl@iA:1= Wx+ 6’At&l @[. I@[-1 = Pxl @Wyl@WLI + 6’ At( [My] @ [M,] @ [Dz] + . . .)
+ (0At)2([~x1 @[Dyl@RI + . .I + (~~~)‘(Px1@pv1@Pzl) = [M + QAtD] + (0At)2[.] + (OAt)3[.]
where Cr= i represents the assemblage over all elements (m). For development simplicity, we can express eqn (6) as [M + 0AtD]{A$}
= -At[D]{$}
+ At[M]{F}
(10)
with the definitions,
(14)
Comparing eqns (6) and (14) the tensor product formation of ([Ml + 6 At[D]) produces a truncation error of order O(At)2 plus higher order terms. Therefore, only for 8 = 0 (explicit integration) do these errors vanish, while for 0 > l/2 an upper limit on useful At exists for a given problem. Conversely, in proceeding to a steady-state solution, this truncation error will vanish since (A4) --+ (0) in the limit of the right side of eqn (12) going to zero. The goal now is to rearrange eqn (12) into an efficient numerical procedure. Two choices exist: ‘approximate factorization’ and ‘time-splitting’.4 The former (AF) replaces eqn (12) by the computational sequence,
(11)
Further, [D] = [D(t + At)] and {F} = O{F(t + At)} + (1 - f3){F(t)}, should the advection velocity and source terms be given as functions of time. Recalling that eqns (10) and (11) are actually formed by assembly over each finite element, eqn (6) can be expressed as
(12) where
Back-substituting these three expressions confirms that eqn (15) has replaced the tensor product 8 with the algebraic product [A,] [Ay] [A,]. The computational advantage gained by eqn (15) is that each numerical solution involves only a banded matrix inversion, e.g. MA lAyIandW are each only tridiagonal using onedimensional linear basis functions. This is a tremendous advantage over the full three-dimensional matrices defined in eqns (6) and (12). Hence, use of relatively smaller time steps (At) to control the truncation error is more than fully compensated by the orders of magnitude speed improvement of the solution. The time-split (TS) numerical solution algorithm is even faster than the AF algorithm. The TS computa-
[AxI= 2 [Axle= 2 [M-x+ QAtD,], e=l
e=l
[A,]= 2 [A,le= 2 WY+ 0AtD,l, e=l e=l
(13)
L&l = 2 [Ale= 2 [M,+ OAtD,l, e=l
e=l
Even though the tensor matrix product is exact, eqn (12) defines an approximation. Substituting eqn (13)
(16)
The relative computational advantage is that eqn (16) involves only one-dimensional matrices (and matrix products) to form both the left and right side terms. It also requires formation and use of the intermediate
18
D. W. Pepper
solution updates {&} with the definitions
with the boundary constraint a4 ^
There is a tendency to associate {$1} and {$2} with {+(t + At/3)} and {$(t + 2At/3)} but this is incorrect due to the right side factorization error. As with the AF solution sequence, the TS algorithm uses only tridiagonal matrices which are very efficient. Only one right side computation is required for the AF solution, but it uses three-dimensional element matrices. Conversely, the TS algorithm requires three right side evaluations, following matrix additions to form {&}, but only one-dimensional element matrix products. Both algorithms require three matrix solutions to advance the overall solution by At. Finally, the TS algorithm introduces an additional second order 0(19At)~ truncation error that may not vanish in the limit {A#} -+ (0). H ence, the TS solution sequence may produce a (bounded) oscillation over a range of { A4) as steady state solution is approached.7
3. A TIME-SPLIT ONE-DIMENSIONAL
RECURSION RELATION ELEMENTS
m$+k,,~n,--q=O
(19)
The one-dimensional
finite element analog of eqn (18) is
{I:> [MxI{$.} 1 +u6 (I;} K
1
-1
+ AX(~) -1
1
-3
{:}{;} “”
1{A) = {.fYI
(20)
where $ = d$i/dt. A recursion relation is established by ‘assembling’ eqn (20) over two adjacent elements, as shown in Fig. 1. Ontheintervalform-ltoi(Z(e)=xi-xi-~ PAX-)
USING
A recursion relation can be derived for the onedimensional time dependent advection-diffusion equation using linear basis functions and global assemblage over two adjacent elements (grid intervals). The basis functions have the appearance of peaked hats, or ‘chapeau’ functions (a common reference by geophysical fluid dynamists). Examples of the use of chapeau function solutions to geophysical problems are reported by Raymond & Garder,’ Long & Hicks,’ Pinder & Gray,” and Long & Pepper. 11J2 Utilizing the TS procedure, and a basic one-dimensional algorithm, multi-dimensional problems can be solved without additional programming complexity. A threedimensional, TS chapeau function procedure is discussed by Pepper et a1.13 To illustrate the procedure, the one-dimensional form of eqn (1) is (18)
If the same procedure is applied on the interval i to i + 1 (AX+ = xi+ 1 - xi) eqn (20) becomes (after global assemblage over the two intervals)
;
I I 2Ax-
Ax-
Ax-
2Ax- + 2Ax,
0
Ax+
-2Uj&l
+i
-Vi-l
- Vi -2Ui
0
3-1 4i +Kz . &+I I
2Uj-1 + Vi
0
Vi-1 - Vi+,
2Ui + Ui+l
-Ui-2Ui+l
Uif2Ui+f
-1 Ax-
- 1
Ax1
- 1
I
0 - 1
-1
Ax- + Ax, 0
-1 AX+
Ax+ (22)
i-
Fig.
1
i
i+l
1. Assemblageover two adjacent one-dimensional
elements.
The equivalent difference-differential
equation for node i
A time-split finite element scheme,for the PC
Equation (23) holds for variable element length and non-uniform velocity, and is tridiagonal, i.e. only the i - 1, i, and i + 1 grid points are utilized. If Ax_ = Ax,, f; = 0, and the velocity is constant, eqn (23) reduces to
(the middle row of the 3 x 3 coefficient matrices) is i [$i- FAX.. + 2$i(Ax4i-l(-~Ji&l
+ AX+) + Ji+rAx+] +4iC”i-l
-2ui)
+$i+1(2ui
-
19
ui+l) I
-I- ut+l)
Notice that the initial value terms are obtained from the mass matrix, and is not ‘lumped’. Expressing the time derivative term (4) by an implicit one-step difference and employing the 0 method, eqn (24) can be rewritten as
(4 2.5 r
; [(qqy - &- 1) + 4(&+’
lj//, 7 0.0 0.5
1.0 X
, , 1.5
+&[H(4i+l
+
2
(b) 2.5
r
-“)(4i+l
-#i-I)“]
-2$j+$b-,)"I=0
(25)
where the superscript n denotes present value at time level n. As before, best results are obtained for 0 = l/2 (Crank-Nicolson). Performing a Taylor series expansion about Xi, eqn (25) (with K = 0) can be transformed to the relation z+%-
2.5 r
-di-l)“+’
+ (1 - u4i+,
84 Cc)
(l
- 4;) + (qy,‘: - &,)]
~3)
UAx4 a’$ - + O(At’, Ax’) 180 a.G
(26)
which indicates that the scheme is fourth order accurate in space and second order accurate in time for uniform element length and constant velocity. The technique easily extends to two- or threedimensions by writing a series of one-dimensional equations. For example, eqn (1) is split into three onedimensional equations such that (27)
Cd) 2.5 ,-
(28)
(29)
Fig. 2. Advection of concentration in one-dimension with CT= 0.4 (dashed lines - analytical solution). (a) Initial distribution; (b) Advection only; (c) Advection with K = 0.01 for diffusion and (d) Advection with K = 0.001 for diffusion.
where superscript ’ and ++ denote interim values. Equations (23) and (27)-(29) have been used by Pepper et al., l3 Long & Pepper,12 and Pepper6 to model three-dimensional transport. Equation (25) is used to solve the advection (K = 0) of a ‘cosine hill’ passive scalar in one-dimension. Figure 2(a) shows the initial distribution at t = 0. Figure 2(b)
20
D. W. Pepper
(a)
progressively more dramatic. The distribution is not significantly altered at t = 150; it should ideally remain the same as at t = 0. Extension of the numerical technique to three dimensions is portrayed in Fig. 4; in this example, a spherical distribution is advected diagonally across a cubical enclosure. Grid volumes are equally incremented, and the velocity vector is uniform. Particles are used to enhance visualization of the distribution, with each particle representing a unit value of the distribution. The particles are plotted within each element volume based on the mass distribution calculated within the element volume from eqns (27)(29). The scatter of particles in Fig. 4(b) is caused by the generation
maximum
of numerical
dispersion.
The initial
value at the center of the spherical distri-
bution is 100 and ideally should remain constant
Fig. 3. Two-dimensional advection of a cosine hill passive scalar. (a) Initial distribution and (b) Distribution after 150 time steps.
shows the analytical solution compared with the solution obtained from eqn (25). The Courant number (cr = u At/Ax) is equal to 0.4. For small Courant numbers, the comparison between the two solutions is nearly identical. The small wake created by the finite element solution is caused by numerical dispersion errors. I2 Raymond & Garder’ and Pepper et ~21.,‘~ discuss ways of eliminating dispersive wakes. As the Courant number increases, more dispersion occurs. The inclusion of diffusion is shown in Fig. 2(c) and 2(d), for K = 0.01 and K = 0.001, respectively. The lateral spreading and damping of the distribution occur as expected because of the inclusion of diffusion (as K increases, more ‘physical damping’ of the solution occurs). The ability of the scheme to solve two-dimensional problems is analyzed by calculating the advection of a two-dimensional cosine hill, shown in Fig. 3. The Courant numbers in the x and y directions are equal to 0.10. The initial distribution is shown in Fig. 3(a) with a subsequent direction distribution shown at t = 150 in Fig. 3(b). In this test, eqn (27) was solved for values in the x direction and then eqn (28) was solved for values in the y direction (analogous to an AD1 technique). The mesh size was 33 x 33. Core requirements were minimal (since only a tridiagonal solution was required at one time, approx. 120KB were needed). This problem was also solved with a conventional finite element program using bilinear isoparametric quadrilateral elements; core requirements were approx. 700KB and running times were about an order of magnitude greater. However, as core size increased, the comparisons became
as it
is advected within the enclosure. The peak value in Fig. 4(b) is 99, with the distribution still remaining essentially spherical. (4
(b)
9 ?:
Fig. 4. Advection of a 3-D spherical domain (each particle represents a unit value). (a) Initial distribution and (b) Distribution after 150 time steps.
21
A time-split jinite element scheme for the PC (a) 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 do you
3
4
5
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 10 15 2 20 41 50 10 41 72 85 15 50 85 100 10 41 72 85 2 20 41 50 0 2 10 15 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 want graphical
6
7
0 0 0 0 0 0 0 0 0 0 0 0 0 0 10 2 41 20 72 41 85 50 72 41 41 20 10 2 0 0 0 0 0 0 0 0 0 0 0 0 output7
8
9
10
11
12
13
14
15
16
17
18
19
20
0 0 0 0 0 0 0 0 2 lo 15 lo 2 0 0 0 0 0 0 0 (y/n)?
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
3D VIEWOFCOSINEHILL - INITIAL DISTRIBUTION
plot is finished - press anykey
Fig. 5. Advection 4. ADVECTION COSINE HILL
of a cosine hill distribution.
OF A TWO-DIMENSIONAL
A simple TS program is included (see Appendix) to illustrate the solution technique for unsteady twodimensional problems using tensor products of onedimensional linear elements. The computer code consists of two separate programs; the first program is called COSINE.BAS and calculates the solid body rotation of a two-dimensional ‘cosine hill’ distribution on a 20 x 20 mesh. The second program is named PLOTHILL.BAS and plots the results from COSINE.BAS on the monitor screen as a three-dimensional view of the cosine hill distribution r (time) = 0 and at 12= 20 time steps. Both programs were written using Quick Basic and compiled on an IBM compatible 486, 66 MHz PC.
I (a) Initial
distribution
and (b) Perspective view
COSINE.BAS solves the two-dimensional time-split advection equation using the one-dimensional recursion relation on a uniform mesh with constant advection velocity. The two-dimensional advection equation is
dc dc de dr+ua.w+w-=o aY
(30)
where c is the scalar field distribution and u and ‘u are the X, y Cartesian velocity components. The time-split algorithm effectively partitions eqn (30) into strictly x and y components,
22
D. W. Pepper (a> 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 do you
3
4
5
0 0 0 0 0 0 0 -0 0 0 1 0 0 -0 -0 -0 0 0 0 -0 0 -0 1 -0 0 -1 -0 0 -0 1 -1 0 -1 1 0 -1 -0 -0 1 -0 -1 1 0 -0 1 -1 0 0 -10 0 0 0 0 0 0 0 -0 0 -0 want graphical
0 0 -0 0 0 -0 -0 0 0 -1 1 1 -0 -1 -1 -0 0 0 1 1
6
7
8
9
10
0 0 0 0 0 0 0 0 0 0 -1 0 1 -0 -0 0 0 -1 -0 1 -0 -0 1 0 -1 -0 0 -0-o 10 0 -0 1 -0 -0 -0 0 -1 1 -0 -1 0 -0 -0 0 1 -1 0 -0 -0 0 -1 0 -0 -0 -2 1 -0 -0 1 -0 2 -2 1 2 0 0 -0 -0 2 -2 2 -0 2 -1 -1 1 -1 -0 1 -2 2 1 -2 -2 -1 -0 1 0 -1-l -1 0 2 0 -3 -0 -1 0 2 2 output? (y/n)?
I
11
12
13
14
15
16
17
18
19
20
0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 1 1 2 2 2 4 5 5 1 4 9 14 14 2 9 21 31 32 3 15 30 59 62 2 18 51 84 92 0 14 47 84 100 -2 4 26 54 71 -3 -2 5 16 27 -2 -2 -5 -5 -1 1 -1 -4 -9 -8 1 1 0 -5 -4 2 4 2 1 -0 0 2 3 2 -2 -0 2 1 2 -1 -2 -3 -2 -0
,” -0 1 1 4 10 23 45 72 82 64 29 1 -8 -7 -1 4 2 -0
0 0 0 0 1 2 5 12 24 38 45 33 10 -6 -9 -4 2 3 1 -0
0 0 -0 -0 0 10 2 4 8 14 15 10 -1 -9 -9 -4 2 3 0 -1
0 0 0 -0 -0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 -0 1
0 1 2 2 0 -5 -9 -9 -3 2 4 1 -1 -0
3DVIEWOFCOSINE HILL- MTER28 TIME STEPS
plat is finished - press anykeyI Fig. 6. Advection of a cosinehill after 20 time steps.(a) Final cosinehill distribution and (b) Perspectiveview. Equation (31) is first solved for the intermediate
vector c+, which is then substituted into eqn (32) to solve for ‘+l, the solution at the new time step, t”l. Equation F31) takes the final solution form n Ax c;:; - cl- 1 it+1 - cy n+l - ci+l ci+l + 4ci 2At + 6 2At 2At
1
+ (Ui-1 - Ui+l)(Cy+l -CT) + (2Ui + Vi+ l)(Cyz: - CY+I)] = 0
(33)
Equation (33) can be rewritten as: Aci-l+Bci+Cci+l=D,
(34)
where \ A = -&2Ui B=4-
+ Uj+l) - 1
gpi+,-
u,-1)
C=&2Ui+
Vi-l)-1
D,=c~+I
1-&(2Ui+U,+l)
(
(35) >
+Cr
+cy-‘_l ( 1+&(ZLi,
> + Ui-1)
)J
)
A time-splitjinite
element scheme for the PC
The corresponding form for eqn (32) is obtained by replacing Ax with Ay, U with V, and i f 1 with j f 1. The cosine hill distribution has been used extensively to test the accuracy of numerical advection schemes. Under the applied pure rotational velocity field, the hill should be advected without distortion, which turns out to be a very difficult test. Essentially all numerical solution procedures can exhibit fairly large degrees of inaccuracy, and the linear finite element technique is no exception. For any algorithm, the level of distortion is directly dependent on the Courant number. Conventional finite difference methods do not solve this deceptively simple appearing problem very well, but the linear finite element method does better than most other algorithms. More sophisticated (and complex) methods do perform better, but they also have their limitations, and cost more to execute. The initial hill distribution is defined as C;,j
= 50 1 +
[
(i - x0)* + (j - YO)* COS
(”
I
r0
i, j=
1,n
(36)
where x0 = 5, y. = 10, and r. = 4. This creates a distribution with a peak value of 100, centered at node i = 5, j = 10 and decreasing to zero at node i = 1 and 9, j = 6 and 14, i.e. the radius of the distribution is 4 cell (x) increments. The rotational velocity distribution is calculated from a rigid body rotation field as: ui,=(j-lO)g
alli
Vi=(lO-i)g
allj J
1>
(37)
The program PLOTHILL.BAS is a graphics program which plots the output from COSINE.BAS. The program can be combined with COSINE.BAS. The most demanding computing time is spent in solving the transport equations in COSINE.BAS. Hence, it is advantageous to compile this part of the program and dump the output to a file. The execution time of PLOTHILL.BAS is relatively fast. PLOTHILL.BAS displays a three-dimensional picture of the initial cosine hill distribution and the final distribution after 20 time steps, as shown in Figs 5(a,b)6(a,b). The final distribution shows a significant amount of information about the accuracy of the finite element method. Small ripples which appear in the distribution are due to the computational dispersion associated with the linear finite element algorithm. This numerical dispersion arises as a consequence of the inability of the method to transport the short wave distributions inherent in any numerical method. Also observe the lack of computational damping whereby the peak value is still 100; however, the distribution is spread over a wider radius of cells. This is due principally to truncation
23
error. We also see a slight phase lag in the solution (although this may not be evident at first glance). Phase lag is the result of a numerical method transporting the distribution (waves) at a slower speed than the true advection velocity; hence, the numerical solution lags behind the actual solution. A more detailed analysis of these types of errors as they affect finite element methods is discussed by Pepper & Baker’ and Long & Pepper.** Values for cij are obtained at each i, j (or x, y) location in COSINE.BAS using the TS algorithm as the solution proceeds in time. The three-dimensional representation of the results is obtained by equating the z coordinate values to the cii values. These values are then plotted as variations in the x, y plane, similar to plotting topographic projections on a flat surface. The program is divided into two major parts. The first part reads the data from the disk file created by COSINE.BAS. The second part plots the distributions at t = 0 and t = 20 on the monitor screen, utilizing a series of subroutines to define the picture area and profiles. In the enclosed example listing, the initial distribution is rotated within a 20 x 20 mesh. The solution is stopped after 20 time steps, whereby the distribution has been rotated approx. 180”. The solution required very little computational time. First program
Load the programs into drive A, B, or C (hard disk preferred). COSINE displays a message and begins calculating the advection, listing each time step up to 20. please wait, solution is proceeding time step= 1 time-step= 2 time-step= 3 time-step= 4 time-step=
20
At the end of 20 time steps, the program terminates and the C > prompt appears on the monitor. Files INIT.DAT and FINAL.DAT are created. Second program
To examine the results of the calculation and view a three-dimensional perspective of the distribution at t = 0 and t = 20, type PLOTHILL. The user inputs INIT.DAT or FINAL.DAT to examine the initial or final solution results, as shown in Figs 5(a,b) and 6(a,b). If ‘y’ is typed to request graphical output, the program displays the final result on the monitor screen. After the graphical display is drawn on the monitor. the program terminates and returns the user to system operation.
D. W. Pepper
Suggested computer program exercise 1. Following the procedure for executing the advection of the two-dimensional cosine hill, rerun the example problem with At = 0.40 (this increases the Courant number) and compare results with those obtained for At = 0.10. 2. Alter the code so that a two-dimensional square block is advected from the lower left corner of the computational domain to the upper right corner. Make alterations to the code COSINE.BAS as follows: u= 1.0 ?I= 1.0 cO(t = 0) = 1.O for i = 2, 3; j = 2, 3
Examine the block after 20 time At = 0.10, and (b) At = 0.50.
steps when (a)
5. CONCLUSIONS There are several ways in which the conventional finite element method can be enhanced to decrease computing time and/or storage requirements. The most commonly used techniques include mass lumping of the time dependent terms, reduced quadrature for numerically evaluating the integral terms, and the method of fractional steps (time-splitting) utilizing one-dimensional tensor matrix products. Application of any of these techniques can lead to significant improvement in computational efficiency over the conventional finite element procedure. Both mass lumping and reduced quadrature can be easily applied to practically any finite element program with only a few changes in those program statements dealing with assemblage and mass matrix formulation. The time-split procedure deals principally with the one-dimensional algorithm, and the subsequent multiplication (tensor matrix products) of the one-dimensional components. A one-dimensional recursion relation can be derived which, from outward appearance, has the look of a centered finite difference technique. However, close examination reveals that the relation is a finite element algorithm which is equivalent in accuracy to the conventional finite element method, but with the speed and reduced storage requirements of the finite difference method.
Application of these alterations is becoming more widespread. Pepper & Baker14 show that over twentyfour different numerical techniques can be derived from the Galerkin weak statement procedure and Taylor series expansions for the basis functions.
REFERENCES 1. Gresho,P. M., Chan, S., Upson,C. & Lee, R., A modified finite element method for solving the time-dependent, incompressibleNavier-Stokes equations. Znt. .Z. Numer. Meth. Fluids, Part I: Theory, 1984,4, 557-98. 2. Patankar, S. V., Numerical Heat Transfer and Fluid Flow. Hemisphere,Washington,DC., 1980. 3. Prakash, C. & Patankar, S. W., A control volume based finite element method for solving the Navier-Stokes equation using equal order variable interpolation. Num. Heat Transfer, 1985,8, 259-80. 4. Yanenko, N. A., The Method of Fractional Steps. Springer, Heidelberg,1971. 5. Pepper, D. W. & Baker, A. J., A simpleone-dimensional finite element algorithm with multidimensionalcapabilities. Num. Heat Transfer, 1979,2, 81-95. 6. Pepper, D. W., Modeling of natural convection in three dimensionsusinga time-splitfinite elementmethod. Num. Heat Transfer, 1987, 11, 31-55. 7. Fletcher, C. A. J., Computational Techniques for Fluid Dynamics. Vol. I, Springer,New York, 1988. 8. Raymond, W. H. & Garder, A., Selectivedampingin a Galerkin methodfor solving wave problemswith variable grids. Mon. Weather Rev., 1976,104, 1583-90. 9. Long, P. E. & Hicks, F. J., Simplepropertiesof Chapeau functions and their application to the solution of the advection equation. NOAA TDL Office note 75-8, Silver Spring, MD, 1975. 10. Pinder, G. F. & Gray, W. G., Finite Element Simulation in Surface and Subsurface Hydrology. Academic Press,New York, 1977. 11. Long, P. E. & Pepper, D. W., A comparison of six numericaladvection schemes for calculatingthe advection of atmosphericpollution, Proc. AMS 3rd Symp. Atmospheric Turbulence Dtjiision and Air Quality, Raleigh, NC, 19-22 Oct. pp. 181-7, 1976. 12. Long, P. E. & Pepper, D. W., An examination of some simplenumericalschemes for calculatingscalaradvection. J. appl Meteor., 1981,20(2), 146-56. 13. Pepper,D. W., Kern, C. D. & Long, P. E., Modeling the dispersionof atmosphericpollution usingcubicsplinesand Chapeaufunctions. Atmos. Environ., 1978,13, 223. 14. Pepper, D. W. & Baker, A. J., Chapter 13: Finite differencesvs finite elements,Handbook of Numerical Heat Transfer, ed. W. Minkowycz et al. J. Wiley and Sons, New York, 1988.
A time-split finite element scheme for the PC APPENDIX QUICK
BASIC PROGRAM COSINE
LISTINGS
25
26
D. W. Pepper
Advances
0965-9978(95)00022-4
ELSEVIER
in Engineering
A time-split finite element scheme for the PC Darrell W. Pepper Depurtment
of
Mechanical Engineering, University
of Nevada,
Las Vegas, Las Vegas, NV 891.54-4027, (IsA
(Received 14 February 1994; revised version received 27 April 1995; accepted 12 May 1995) A time-split finite element scheme is developed based on one-dimensional, linear basis functions. The transient form of the transport equation is discretized and a recursion equation obtained through assemblage over two adjacent elements. The resulting one-dimensional scheme can be readily applied to multi-dimensional problems, requiring only the solution of a simple tri-diagonal banded matrix (in each direction). Computer program listings are included for solving twodimensional advection and displaying a surface plot of a rotating cosine hill distribution on a PC. Key words: time-split, finite element, PC. recursion relation.
1. INTRODUCTION
mass matrix to be easily calculated. However, the use of mass lumping has a tendency to introduce some error in the transient analysis; for problems which converge to steady state, the steady state results are essentially as accurate as those obtained without lumping. Another technique which yields reasonably good results uses one-point (Gauss) quadrature in lieu of 2 x 2 quadrature for four node quadrilaterals (or 2 x 2 x 2 for eight node hexahedrals). Utilizing mass lumping in an explicit time marching scheme with onepoint quadrature, Gresho et al.’ showed that a modified finite element method could produce results comparable in accuracy to finite differences. While not quite the speed of the finite difference method, the technique was significantly faster than the more conventional finite element method. The modifications made by Gresho et al.’ to the finite element method have the effect of altering some of the terms and matrix operations to behave similarly to the finite volume method.2 Development of a finite element-finite volume method was conducted by Prakash & Patankar’ in an effort to more formally combine the finite element method with the finite volume technique; two-dimensional elements were used. While this technique exhibited the advantages of geometrical flexibility and computational speed, there was some loss of accuracy. The concept of time-splitting can be applied to the finite element method to yield a finite element algorithm with the advantages of the finite difference technique. The method of time-splitting, or method of fractional steps4 creates a series of one-dimensional relations which can be multiplied (tensor product) to yield solutions to multi-dimensional problems. The
The advantages of the finite element method over the finite difference method are well known: geometrical flexibility, computational accuracy and stability, ease in altering boundary conditions without requiring changes in the program, and ability to easily use various order interpolation functions. However, the one disadvantage which clearly becomes evident is the amount of computational speed and storage required for problems involving many nodes. In an implicit difference method, the use of one-sided and central differencing techniques to approximate derivative terms produces a set of symmetric banded matrices consisting of usually three to five diagonal terms. The Thomas algorithm is a simple procedure which quickly solves tri-diagonal matrices without involving the zero value off-diagonal terms. If an explicit time scheme is used, the mass matrix associated with the time derivative term in the finite difference technique consists of only one value; hence, the matrix equivalent equation can be solved without the need for a matrix solver. The mass matrix of the finite element method, on the other hand, is symmetric and banded according to the elemental node numbering, regardless of the type of time differencing employed. Over the last few years, efforts have been made to alter the finite element method in order to reduce computational storage requirements and computing times, without significantly decreasing accuracy. The use of mass lumping is one such method whereby the off-diagonal terms are lumped into a single, main diagonal. This procedure permits the inverse of the 15
16
D. W. Pepper
technique has the advantages of simplicity, i.e. onedimensional element formulation, and computation speed, i.e. tri-diagonal matrix solution common to finite difference methods. The time-split algorithm has been used to solve multi-dimensional problems involving two- and three-dimensional scalar transport and fluid dynamics, and has been shown to produce accurate results at computational speeds equivalent to finite difference methods.5’6
2. THE MULTI-DIMENSIONAL ALGORITHM
with a4
+
L
*
%nx+k
a4
^
a4
-n,,+k,~nz+q=O yy ay
n
$(X,Y, z, 1= 0) = $0
= -WA + by1+ M + kl + FYI+ M&V + [M]At(O{F}“+l
TIME-SPLIT
The three-dimensional transport equation for a scalar quantity (e.g. temperature) can be written in the form
4
where [vJ, [r+,], [Q] denote the advection matrices for the u, w, and w velocity components, and {Fi} is the load vector consisting of known values (FXFJZ). Approximating the time derivative term using the 6 method, eqn (5) becomes
(2)
(3)
where eqns (2) and (3) represent the generalized boundary condition constraints and initial values for eqn (1). The finite element weak statement of eqn (1) is obtained as
+ (1 - e){F}“)
(6)
where 0 5 13I 1 (if B = l/2, one obtains CrankNicolson time averaging). The one-dimensional algorithm yields for the left hand side of the matrix equation, ([Ml + BAt( [w,] + [k,])), which is only tri-diagonal using linear one-dimensional elements, or penta-diagonal using quadratic elements. In contrast, the left hand side (matrix) of eqn (6) is ‘very large’, since it contains the three-dimensional coupling of all derivatives. Formulation of this matrix and subsequent solution of eqn (6) can quickly saturate the available memory of a computer using even a coarse mesh and the conventional finite element procedure. However, a number of methods exist which can be used to more efficiently solve eqn (6), such as Gauss-Seidel iteration, alternatingdirection implicit (ADI), approximate factorization (AF), and operator-splitting. In each instance, the multi-dimensional derivative operators are approximated by a sequence of one-dimensional operations. The tensor matrix product of two vectors produces a matrix of higher rank (or order), e.g. the multiplication of two 2 x 2 matrices produces a 4 x 4 matrix. The product of two vectors, e.g. two basis functions, produces a higher order vector. The basic building block is the one-dimensional mass matrix [Ml. Refining notation for n-dimensions, the onedimensional (x-component) form is written as (7)
=
qNi dS (4 If s Note that eqn (4) is ‘functionally’ identical to a onedimensional finite element relation; the three-dimensional form simply generalizes Ni and the advection and diffusion terms possess scalar components parallel to each component of the x, y, z coordinate system. Using subscripts to denote the x, y, z terms, the matrix equivalent expression for eqn (4) is f ff
QNi dx dy dz -
4 2 =-
IN{ $}+(Lvxl + Lwyl + [w.zl){4i1 + ([kxl+ [kyl+ [kI){A) = {Fi:i)
where Ax E l(‘), which is the element length in the xdirection. Notice that the element order and rank are both equal to 2. The two-dimensional mass matrix is constructed using the tensor product (denoted by @) of the two one-dimensional mass matrices, [M,] and [My], i.e.
(5)
AxAy 36
1 2
2 4 2 1 1 2 4 2 2
1 2 4
(8)
A time-split jinite element scheme for the PC which is analogous to the mass matrix relation for the two-dimensional linear quadrilateral element. The threedimensional mass matrix is similarly formed by multi-
plying ew (8) by NJ. The true value of the tensor product is accurately and efficiently approximates the left eqn (6). Every matrix given in eqn (6) is computable using [M,], [MY], [M,], the tensor product, e.g.
that it side of directly matrix
17
into (12) the left side of eqn (12) becomes
[A,I @M,pl@iA:1= Wx+ 6’At&l @[. I@[-1 = Pxl @Wyl@WLI + 6’ At( [My] @ [M,] @ [Dz] + . . .)
+ (0At)2([~x1 @[Dyl@RI + . .I + (~~~)‘(Px1@pv1@Pzl) = [M + QAtD] + (0At)2[.] + (OAt)3[.]
where Cr= i represents the assemblage over all elements (m). For development simplicity, we can express eqn (6) as [M + 0AtD]{A$}
= -At[D]{$}
+ At[M]{F}
(10)
with the definitions,
(14)
Comparing eqns (6) and (14) the tensor product formation of ([Ml + 6 At[D]) produces a truncation error of order O(At)2 plus higher order terms. Therefore, only for 8 = 0 (explicit integration) do these errors vanish, while for 0 > l/2 an upper limit on useful At exists for a given problem. Conversely, in proceeding to a steady-state solution, this truncation error will vanish since (A4) --+ (0) in the limit of the right side of eqn (12) going to zero. The goal now is to rearrange eqn (12) into an efficient numerical procedure. Two choices exist: ‘approximate factorization’ and ‘time-splitting’.4 The former (AF) replaces eqn (12) by the computational sequence,
(11)
Further, [D] = [D(t + At)] and {F} = O{F(t + At)} + (1 - f3){F(t)}, should the advection velocity and source terms be given as functions of time. Recalling that eqns (10) and (11) are actually formed by assembly over each finite element, eqn (6) can be expressed as
(12) where
Back-substituting these three expressions confirms that eqn (15) has replaced the tensor product 8 with the algebraic product [A,] [Ay] [A,]. The computational advantage gained by eqn (15) is that each numerical solution involves only a banded matrix inversion, e.g. MA lAyIandW are each only tridiagonal using onedimensional linear basis functions. This is a tremendous advantage over the full three-dimensional matrices defined in eqns (6) and (12). Hence, use of relatively smaller time steps (At) to control the truncation error is more than fully compensated by the orders of magnitude speed improvement of the solution. The time-split (TS) numerical solution algorithm is even faster than the AF algorithm. The TS computa-
[AxI= 2 [Axle= 2 [M-x+ QAtD,], e=l
e=l
[A,]= 2 [A,le= 2 WY+ 0AtD,l, e=l e=l
(13)
L&l = 2 [Ale= 2 [M,+ OAtD,l, e=l
e=l
Even though the tensor matrix product is exact, eqn (12) defines an approximation. Substituting eqn (13)
(16)
The relative computational advantage is that eqn (16) involves only one-dimensional matrices (and matrix products) to form both the left and right side terms. It also requires formation and use of the intermediate
18
D. W. Pepper
solution updates {&} with the definitions
with the boundary constraint a4 ^
There is a tendency to associate {$1} and {$2} with {+(t + At/3)} and {$(t + 2At/3)} but this is incorrect due to the right side factorization error. As with the AF solution sequence, the TS algorithm uses only tridiagonal matrices which are very efficient. Only one right side computation is required for the AF solution, but it uses three-dimensional element matrices. Conversely, the TS algorithm requires three right side evaluations, following matrix additions to form {&}, but only one-dimensional element matrix products. Both algorithms require three matrix solutions to advance the overall solution by At. Finally, the TS algorithm introduces an additional second order 0(19At)~ truncation error that may not vanish in the limit {A#} -+ (0). H ence, the TS solution sequence may produce a (bounded) oscillation over a range of { A4) as steady state solution is approached.7
3. A TIME-SPLIT ONE-DIMENSIONAL
RECURSION RELATION ELEMENTS
m$+k,,~n,--q=O
(19)
The one-dimensional
finite element analog of eqn (18) is
{I:> [MxI{$.} 1 +u6 (I;} K
1
-1
+ AX(~) -1
1
-3
{:}{;} “”
1{A) = {.fYI
(20)
where $ = d$i/dt. A recursion relation is established by ‘assembling’ eqn (20) over two adjacent elements, as shown in Fig. 1. Ontheintervalform-ltoi(Z(e)=xi-xi-~ PAX-)
USING
A recursion relation can be derived for the onedimensional time dependent advection-diffusion equation using linear basis functions and global assemblage over two adjacent elements (grid intervals). The basis functions have the appearance of peaked hats, or ‘chapeau’ functions (a common reference by geophysical fluid dynamists). Examples of the use of chapeau function solutions to geophysical problems are reported by Raymond & Garder,’ Long & Hicks,’ Pinder & Gray,” and Long & Pepper. 11J2 Utilizing the TS procedure, and a basic one-dimensional algorithm, multi-dimensional problems can be solved without additional programming complexity. A threedimensional, TS chapeau function procedure is discussed by Pepper et a1.13 To illustrate the procedure, the one-dimensional form of eqn (1) is (18)
If the same procedure is applied on the interval i to i + 1 (AX+ = xi+ 1 - xi) eqn (20) becomes (after global assemblage over the two intervals)
;
I I 2Ax-
Ax-
Ax-
2Ax- + 2Ax,
0
Ax+
-2Uj&l
+i
-Vi-l
- Vi -2Ui
0
3-1 4i +Kz . &+I I
2Uj-1 + Vi
0
Vi-1 - Vi+,
2Ui + Ui+l
-Ui-2Ui+l
Uif2Ui+f
-1 Ax-
- 1
Ax1
- 1
I
0 - 1
-1
Ax- + Ax, 0
-1 AX+
Ax+ (22)
i-
Fig.
1
i
i+l
1. Assemblageover two adjacent one-dimensional
elements.
The equivalent difference-differential
equation for node i
A time-split finite element scheme,for the PC
Equation (23) holds for variable element length and non-uniform velocity, and is tridiagonal, i.e. only the i - 1, i, and i + 1 grid points are utilized. If Ax_ = Ax,, f; = 0, and the velocity is constant, eqn (23) reduces to
(the middle row of the 3 x 3 coefficient matrices) is i [$i- FAX.. + 2$i(Ax4i-l(-~Ji&l
+ AX+) + Ji+rAx+] +4iC”i-l
-2ui)
+$i+1(2ui
-
19
ui+l) I
-I- ut+l)
Notice that the initial value terms are obtained from the mass matrix, and is not ‘lumped’. Expressing the time derivative term (4) by an implicit one-step difference and employing the 0 method, eqn (24) can be rewritten as
(4 2.5 r
; [(qqy - &- 1) + 4(&+’
lj//, 7 0.0 0.5
1.0 X
, , 1.5
+&[H(4i+l
+
2
(b) 2.5
r
-“)(4i+l
-#i-I)“]
-2$j+$b-,)"I=0
(25)
where the superscript n denotes present value at time level n. As before, best results are obtained for 0 = l/2 (Crank-Nicolson). Performing a Taylor series expansion about Xi, eqn (25) (with K = 0) can be transformed to the relation z+%-
2.5 r
-di-l)“+’
+ (1 - u4i+,
84 Cc)
(l
- 4;) + (qy,‘: - &,)]
~3)
UAx4 a’$ - + O(At’, Ax’) 180 a.G
(26)
which indicates that the scheme is fourth order accurate in space and second order accurate in time for uniform element length and constant velocity. The technique easily extends to two- or threedimensions by writing a series of one-dimensional equations. For example, eqn (1) is split into three onedimensional equations such that (27)
Cd) 2.5 ,-
(28)
(29)
Fig. 2. Advection of concentration in one-dimension with CT= 0.4 (dashed lines - analytical solution). (a) Initial distribution; (b) Advection only; (c) Advection with K = 0.01 for diffusion and (d) Advection with K = 0.001 for diffusion.
where superscript ’ and ++ denote interim values. Equations (23) and (27)-(29) have been used by Pepper et al., l3 Long & Pepper,12 and Pepper6 to model three-dimensional transport. Equation (25) is used to solve the advection (K = 0) of a ‘cosine hill’ passive scalar in one-dimension. Figure 2(a) shows the initial distribution at t = 0. Figure 2(b)
20
D. W. Pepper
(a)
progressively more dramatic. The distribution is not significantly altered at t = 150; it should ideally remain the same as at t = 0. Extension of the numerical technique to three dimensions is portrayed in Fig. 4; in this example, a spherical distribution is advected diagonally across a cubical enclosure. Grid volumes are equally incremented, and the velocity vector is uniform. Particles are used to enhance visualization of the distribution, with each particle representing a unit value of the distribution. The particles are plotted within each element volume based on the mass distribution calculated within the element volume from eqns (27)(29). The scatter of particles in Fig. 4(b) is caused by the generation
maximum
of numerical
dispersion.
The initial
value at the center of the spherical distri-
bution is 100 and ideally should remain constant
Fig. 3. Two-dimensional advection of a cosine hill passive scalar. (a) Initial distribution and (b) Distribution after 150 time steps.
shows the analytical solution compared with the solution obtained from eqn (25). The Courant number (cr = u At/Ax) is equal to 0.4. For small Courant numbers, the comparison between the two solutions is nearly identical. The small wake created by the finite element solution is caused by numerical dispersion errors. I2 Raymond & Garder’ and Pepper et ~21.,‘~ discuss ways of eliminating dispersive wakes. As the Courant number increases, more dispersion occurs. The inclusion of diffusion is shown in Fig. 2(c) and 2(d), for K = 0.01 and K = 0.001, respectively. The lateral spreading and damping of the distribution occur as expected because of the inclusion of diffusion (as K increases, more ‘physical damping’ of the solution occurs). The ability of the scheme to solve two-dimensional problems is analyzed by calculating the advection of a two-dimensional cosine hill, shown in Fig. 3. The Courant numbers in the x and y directions are equal to 0.10. The initial distribution is shown in Fig. 3(a) with a subsequent direction distribution shown at t = 150 in Fig. 3(b). In this test, eqn (27) was solved for values in the x direction and then eqn (28) was solved for values in the y direction (analogous to an AD1 technique). The mesh size was 33 x 33. Core requirements were minimal (since only a tridiagonal solution was required at one time, approx. 120KB were needed). This problem was also solved with a conventional finite element program using bilinear isoparametric quadrilateral elements; core requirements were approx. 700KB and running times were about an order of magnitude greater. However, as core size increased, the comparisons became
as it
is advected within the enclosure. The peak value in Fig. 4(b) is 99, with the distribution still remaining essentially spherical. (4
(b)
9 ?:
Fig. 4. Advection of a 3-D spherical domain (each particle represents a unit value). (a) Initial distribution and (b) Distribution after 150 time steps.
21
A time-split jinite element scheme for the PC (a) 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 do you
3
4
5
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 10 15 2 20 41 50 10 41 72 85 15 50 85 100 10 41 72 85 2 20 41 50 0 2 10 15 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 want graphical
6
7
0 0 0 0 0 0 0 0 0 0 0 0 0 0 10 2 41 20 72 41 85 50 72 41 41 20 10 2 0 0 0 0 0 0 0 0 0 0 0 0 output7
8
9
10
11
12
13
14
15
16
17
18
19
20
0 0 0 0 0 0 0 0 2 lo 15 lo 2 0 0 0 0 0 0 0 (y/n)?
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
3D VIEWOFCOSINEHILL - INITIAL DISTRIBUTION
plot is finished - press anykey
Fig. 5. Advection 4. ADVECTION COSINE HILL
of a cosine hill distribution.
OF A TWO-DIMENSIONAL
A simple TS program is included (see Appendix) to illustrate the solution technique for unsteady twodimensional problems using tensor products of onedimensional linear elements. The computer code consists of two separate programs; the first program is called COSINE.BAS and calculates the solid body rotation of a two-dimensional ‘cosine hill’ distribution on a 20 x 20 mesh. The second program is named PLOTHILL.BAS and plots the results from COSINE.BAS on the monitor screen as a three-dimensional view of the cosine hill distribution r (time) = 0 and at 12= 20 time steps. Both programs were written using Quick Basic and compiled on an IBM compatible 486, 66 MHz PC.
I (a) Initial
distribution
and (b) Perspective view
COSINE.BAS solves the two-dimensional time-split advection equation using the one-dimensional recursion relation on a uniform mesh with constant advection velocity. The two-dimensional advection equation is
dc dc de dr+ua.w+w-=o aY
(30)
where c is the scalar field distribution and u and ‘u are the X, y Cartesian velocity components. The time-split algorithm effectively partitions eqn (30) into strictly x and y components,
22
D. W. Pepper (a> 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 do you
3
4
5
0 0 0 0 0 0 0 -0 0 0 1 0 0 -0 -0 -0 0 0 0 -0 0 -0 1 -0 0 -1 -0 0 -0 1 -1 0 -1 1 0 -1 -0 -0 1 -0 -1 1 0 -0 1 -1 0 0 -10 0 0 0 0 0 0 0 -0 0 -0 want graphical
0 0 -0 0 0 -0 -0 0 0 -1 1 1 -0 -1 -1 -0 0 0 1 1
6
7
8
9
10
0 0 0 0 0 0 0 0 0 0 -1 0 1 -0 -0 0 0 -1 -0 1 -0 -0 1 0 -1 -0 0 -0-o 10 0 -0 1 -0 -0 -0 0 -1 1 -0 -1 0 -0 -0 0 1 -1 0 -0 -0 0 -1 0 -0 -0 -2 1 -0 -0 1 -0 2 -2 1 2 0 0 -0 -0 2 -2 2 -0 2 -1 -1 1 -1 -0 1 -2 2 1 -2 -2 -1 -0 1 0 -1-l -1 0 2 0 -3 -0 -1 0 2 2 output? (y/n)?
I
11
12
13
14
15
16
17
18
19
20
0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 1 1 2 2 2 4 5 5 1 4 9 14 14 2 9 21 31 32 3 15 30 59 62 2 18 51 84 92 0 14 47 84 100 -2 4 26 54 71 -3 -2 5 16 27 -2 -2 -5 -5 -1 1 -1 -4 -9 -8 1 1 0 -5 -4 2 4 2 1 -0 0 2 3 2 -2 -0 2 1 2 -1 -2 -3 -2 -0
,” -0 1 1 4 10 23 45 72 82 64 29 1 -8 -7 -1 4 2 -0
0 0 0 0 1 2 5 12 24 38 45 33 10 -6 -9 -4 2 3 1 -0
0 0 -0 -0 0 10 2 4 8 14 15 10 -1 -9 -9 -4 2 3 0 -1
0 0 0 -0 -0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 -0 1
0 1 2 2 0 -5 -9 -9 -3 2 4 1 -1 -0
3DVIEWOFCOSINE HILL- MTER28 TIME STEPS
plat is finished - press anykeyI Fig. 6. Advection of a cosinehill after 20 time steps.(a) Final cosinehill distribution and (b) Perspectiveview. Equation (31) is first solved for the intermediate
vector c+, which is then substituted into eqn (32) to solve for ‘+l, the solution at the new time step, t”l. Equation F31) takes the final solution form n Ax c;:; - cl- 1 it+1 - cy n+l - ci+l ci+l + 4ci 2At + 6 2At 2At
1
+ (Ui-1 - Ui+l)(Cy+l -CT) + (2Ui + Vi+ l)(Cyz: - CY+I)] = 0
(33)
Equation (33) can be rewritten as: Aci-l+Bci+Cci+l=D,
(34)
where \ A = -&2Ui B=4-
+ Uj+l) - 1
gpi+,-
u,-1)
C=&2Ui+
Vi-l)-1
D,=c~+I
1-&(2Ui+U,+l)
(
(35) >
+Cr
+cy-‘_l ( 1+&(ZLi,
> + Ui-1)
)J
)
A time-splitjinite
element scheme for the PC
The corresponding form for eqn (32) is obtained by replacing Ax with Ay, U with V, and i f 1 with j f 1. The cosine hill distribution has been used extensively to test the accuracy of numerical advection schemes. Under the applied pure rotational velocity field, the hill should be advected without distortion, which turns out to be a very difficult test. Essentially all numerical solution procedures can exhibit fairly large degrees of inaccuracy, and the linear finite element technique is no exception. For any algorithm, the level of distortion is directly dependent on the Courant number. Conventional finite difference methods do not solve this deceptively simple appearing problem very well, but the linear finite element method does better than most other algorithms. More sophisticated (and complex) methods do perform better, but they also have their limitations, and cost more to execute. The initial hill distribution is defined as C;,j
= 50 1 +
[
(i - x0)* + (j - YO)* COS
(”
I
r0
i, j=
1,n
(36)
where x0 = 5, y. = 10, and r. = 4. This creates a distribution with a peak value of 100, centered at node i = 5, j = 10 and decreasing to zero at node i = 1 and 9, j = 6 and 14, i.e. the radius of the distribution is 4 cell (x) increments. The rotational velocity distribution is calculated from a rigid body rotation field as: ui,=(j-lO)g
alli
Vi=(lO-i)g
allj J
1>
(37)
The program PLOTHILL.BAS is a graphics program which plots the output from COSINE.BAS. The program can be combined with COSINE.BAS. The most demanding computing time is spent in solving the transport equations in COSINE.BAS. Hence, it is advantageous to compile this part of the program and dump the output to a file. The execution time of PLOTHILL.BAS is relatively fast. PLOTHILL.BAS displays a three-dimensional picture of the initial cosine hill distribution and the final distribution after 20 time steps, as shown in Figs 5(a,b)6(a,b). The final distribution shows a significant amount of information about the accuracy of the finite element method. Small ripples which appear in the distribution are due to the computational dispersion associated with the linear finite element algorithm. This numerical dispersion arises as a consequence of the inability of the method to transport the short wave distributions inherent in any numerical method. Also observe the lack of computational damping whereby the peak value is still 100; however, the distribution is spread over a wider radius of cells. This is due principally to truncation
23
error. We also see a slight phase lag in the solution (although this may not be evident at first glance). Phase lag is the result of a numerical method transporting the distribution (waves) at a slower speed than the true advection velocity; hence, the numerical solution lags behind the actual solution. A more detailed analysis of these types of errors as they affect finite element methods is discussed by Pepper & Baker’ and Long & Pepper.** Values for cij are obtained at each i, j (or x, y) location in COSINE.BAS using the TS algorithm as the solution proceeds in time. The three-dimensional representation of the results is obtained by equating the z coordinate values to the cii values. These values are then plotted as variations in the x, y plane, similar to plotting topographic projections on a flat surface. The program is divided into two major parts. The first part reads the data from the disk file created by COSINE.BAS. The second part plots the distributions at t = 0 and t = 20 on the monitor screen, utilizing a series of subroutines to define the picture area and profiles. In the enclosed example listing, the initial distribution is rotated within a 20 x 20 mesh. The solution is stopped after 20 time steps, whereby the distribution has been rotated approx. 180”. The solution required very little computational time. First program
Load the programs into drive A, B, or C (hard disk preferred). COSINE displays a message and begins calculating the advection, listing each time step up to 20. please wait, solution is proceeding time step= 1 time-step= 2 time-step= 3 time-step= 4 time-step=
20
At the end of 20 time steps, the program terminates and the C > prompt appears on the monitor. Files INIT.DAT and FINAL.DAT are created. Second program
To examine the results of the calculation and view a three-dimensional perspective of the distribution at t = 0 and t = 20, type PLOTHILL. The user inputs INIT.DAT or FINAL.DAT to examine the initial or final solution results, as shown in Figs 5(a,b) and 6(a,b). If ‘y’ is typed to request graphical output, the program displays the final result on the monitor screen. After the graphical display is drawn on the monitor. the program terminates and returns the user to system operation.
D. W. Pepper
Suggested computer program exercise 1. Following the procedure for executing the advection of the two-dimensional cosine hill, rerun the example problem with At = 0.40 (this increases the Courant number) and compare results with those obtained for At = 0.10. 2. Alter the code so that a two-dimensional square block is advected from the lower left corner of the computational domain to the upper right corner. Make alterations to the code COSINE.BAS as follows: u= 1.0 ?I= 1.0 cO(t = 0) = 1.O for i = 2, 3; j = 2, 3
Examine the block after 20 time At = 0.10, and (b) At = 0.50.
steps when (a)
5. CONCLUSIONS There are several ways in which the conventional finite element method can be enhanced to decrease computing time and/or storage requirements. The most commonly used techniques include mass lumping of the time dependent terms, reduced quadrature for numerically evaluating the integral terms, and the method of fractional steps (time-splitting) utilizing one-dimensional tensor matrix products. Application of any of these techniques can lead to significant improvement in computational efficiency over the conventional finite element procedure. Both mass lumping and reduced quadrature can be easily applied to practically any finite element program with only a few changes in those program statements dealing with assemblage and mass matrix formulation. The time-split procedure deals principally with the one-dimensional algorithm, and the subsequent multiplication (tensor matrix products) of the one-dimensional components. A one-dimensional recursion relation can be derived which, from outward appearance, has the look of a centered finite difference technique. However, close examination reveals that the relation is a finite element algorithm which is equivalent in accuracy to the conventional finite element method, but with the speed and reduced storage requirements of the finite difference method.
Application of these alterations is becoming more widespread. Pepper & Baker14 show that over twentyfour different numerical techniques can be derived from the Galerkin weak statement procedure and Taylor series expansions for the basis functions.
REFERENCES 1. Gresho,P. M., Chan, S., Upson,C. & Lee, R., A modified finite element method for solving the time-dependent, incompressibleNavier-Stokes equations. Znt. .Z. Numer. Meth. Fluids, Part I: Theory, 1984,4, 557-98. 2. Patankar, S. V., Numerical Heat Transfer and Fluid Flow. Hemisphere,Washington,DC., 1980. 3. Prakash, C. & Patankar, S. W., A control volume based finite element method for solving the Navier-Stokes equation using equal order variable interpolation. Num. Heat Transfer, 1985,8, 259-80. 4. Yanenko, N. A., The Method of Fractional Steps. Springer, Heidelberg,1971. 5. Pepper, D. W. & Baker, A. J., A simpleone-dimensional finite element algorithm with multidimensionalcapabilities. Num. Heat Transfer, 1979,2, 81-95. 6. Pepper, D. W., Modeling of natural convection in three dimensionsusinga time-splitfinite elementmethod. Num. Heat Transfer, 1987, 11, 31-55. 7. Fletcher, C. A. J., Computational Techniques for Fluid Dynamics. Vol. I, Springer,New York, 1988. 8. Raymond, W. H. & Garder, A., Selectivedampingin a Galerkin methodfor solving wave problemswith variable grids. Mon. Weather Rev., 1976,104, 1583-90. 9. Long, P. E. & Hicks, F. J., Simplepropertiesof Chapeau functions and their application to the solution of the advection equation. NOAA TDL Office note 75-8, Silver Spring, MD, 1975. 10. Pinder, G. F. & Gray, W. G., Finite Element Simulation in Surface and Subsurface Hydrology. Academic Press,New York, 1977. 11. Long, P. E. & Pepper, D. W., A comparison of six numericaladvection schemes for calculatingthe advection of atmosphericpollution, Proc. AMS 3rd Symp. Atmospheric Turbulence Dtjiision and Air Quality, Raleigh, NC, 19-22 Oct. pp. 181-7, 1976. 12. Long, P. E. & Pepper, D. W., An examination of some simplenumericalschemes for calculatingscalaradvection. J. appl Meteor., 1981,20(2), 146-56. 13. Pepper,D. W., Kern, C. D. & Long, P. E., Modeling the dispersionof atmosphericpollution usingcubicsplinesand Chapeaufunctions. Atmos. Environ., 1978,13, 223. 14. Pepper, D. W. & Baker, A. J., Chapter 13: Finite differencesvs finite elements,Handbook of Numerical Heat Transfer, ed. W. Minkowycz et al. J. Wiley and Sons, New York, 1988.
A time-split finite element scheme for the PC APPENDIX QUICK
BASIC PROGRAM COSINE
LISTINGS
25
26
D. W. Pepper