- Email: [email protected]
A generalized galerkin method for steady convection-diffusion problems with application to quadratic shape function elements
COMPUTER METHODS NORTH-HOLLAND
IN APPLIED
MECHANICS
AND ENGINEERING
48 (1985) 25-43
A GENERALIZED GALERKIN METHOD FOR STEADY CONVECTION-DIFFUSION PROBLEMS WITH APPLICATION TO QUADRATIC SHAPE FUNCTION ELEMENTS
Jean DONEA Applied Mechanics Division, Joint Research Centre-Ispra Establishment, Communities, 27020 Ispra, Varese, Italy
Ted BELYTSCHKO
Commission of the European
and Patrick SMOLINSKI
Department of Civil Engineering, Northwestern University, Evanston, IL 60201, U.S.A.
Revised
Received manuscript
12 July 1983 received 14 July 1984
A generalization of the standard Galerkin finite element method is considered to enable it to deal successfully with steady convection-diffusion problems. The proposed method employs a generalized governing equation which is obtained by subtracting from the original differential equation the scalar product of its gradient by a vector of free parameters associated with each of the coordinate directions. The generalized equation is successively discretized by the standard Bubnov-Galerkin finite element method. The effectiveness of the method is illustrated for the case of quadratic local interpolations in one and two space dimensions.
1. Introduction
Traditionally, the steady convection-diffusion equation has been one of the most difficult equations to solve by numerical methods. The difficulties arise from the unsymmetric nature of the convection operator, which tends to engender physically unrealistic oscillations in the numerical solutions when the transport terms are dominant. This problem arises in both Galerkin/finite element and finite difference formulations with reasonable grid size. For example, when second-order central difference approximations to the convective term are generated, the discrete equations at odd and even nodes tend to become decoupled, leading to an almost singular system of global equations. Several schemes have been devised to effectively deal with these difficulties. The remedy with difference methods has long been to use upwind differencing for the convective term. This was accomplished by shifting the first derivative term in the upwind direction [l, 21. For the simple case of the one-dimensional, constant-coefficient, convection-diffusion equation, an upwind scheme could be devised [l] which gives exact nodal answers. However, problems of accuracy still remained in practical situations through excessive crosswind diffusion and with variable flow fields and source terms. The upwind concept was first introduced in a finite element context by Christie et al. [3], who achieved an upwind character in the finite element equations by skewing the weight functions appropriately. Hughes [4] has achieved similar 00457825/85/$3.30
@ 1985, Elsevier
Science
Publishers
B.V. (North-Holland)
26
T. Belytschko et al., Generalized Galerkin method for convection-difision
problems
results for low-order elements by shifting the numerical integration points in the evaluation of the transport terms. Computationally simple amplification schemes for achieving upwind-type finite element approximations have been proposed by Belytschko and Kennedy [5] and analyzed by Belytschko and Eldib [6]. These early finite element procedures of upwind-type were also subject to the same criticisms as upwind finite differences when applied to multi-dimensional situations or in problems with variable source terms. More recently, however, methods were developed which do not suffer these deficiencies and represent a computation of convection-dominated significant breakthrough in the finite element phenomena. These are the so-called consistent Petrov-Galerkin methods in which a suitably modified weight function is applied to all terms in the equations. Of particular relevance in this context are the papers of Hughes and coworkers (see [7-9]), and of Morton and coworkers (see [lo-121). In the present paper, we propose an alternative method for deriving accurate finite element schemes for the steady convection-diffusion equation which has enabled us to treat higherorder shape function. The proposed approach is inspired by the Taylor-Galerkin method recently developed by Donea [13] for first-order hyperbolic equations. It consists of constructing a generalized governing equation for the convection-diffusion problem which includes all the necessary ingredients to derive upwind-type finite element schemes in a consistent fashion. This generalized equation is then discretized by application of the standard Bubnov-Galerkin process. The proposed procedure is introduced in Section 2 for the simple case of the linear convection-diffusion equation in one space dimension. The resulting scheme for piecewise quadratic local interpolations is presented in Section 3 and shown to be exact under appropriate circumstances. A series of comparison problems is presented in Section 4. Section 5 is devoted to an analysis of the convergence properties of the parabolic element in comparison with those of the linear element, indicating a definitely superior convergence rate for the quadratic shape function element. Finally, the extension of the proposed method to deal with multi-dimensional convection-diffusion problems is discussed in Section 6 and illustrated in Section 7 for the case of biquadratic shape function elements.
2. Generalized
governing
equation for the convection-diffusion
For simplicity, consider first the linear, constant-coefficient one space dimension,
a@
-Udx+Kz
a2@+s=o, OSXSL,
problem
convection-diffusion
equation in
(2.1)
with the following boundary conditions: Q(O) = g 7
(2.2a) (2.2b)
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
27
Here u is the convective velocity, K is a positive diffusion coefficient and s a given source. As is well-known, for large values of the non-dimensional Peclet number uL/K, the application of the Galerkin method to (2.1) often leads to wildly oscillatory solutions. In effect, central difference approximations to the convective term are generated and the resulting system of discrete equations becomes almost singular. This may readily be seen by considering a uniform mesh of linear elements with nodes numbered sequentially. For K = 0 in (2.1) the Galerkin equation at an interior node I in a uniform mesh of elements with length h is - & (@I+, - @I-,> + SI + b(s,-I-
2s, + s,+1) = 0 .
(2.3)
It can immediately be noted that the equations at odd-numbered nodes are fully decoupled from the equations at even-numbered nodes, and vice-versa. This is one explanation for why oscillatory solutions are obtained when solving (2.1) for a high value of the Peclet number. The easiest way to couple the discrete equations (2.3) and thereby suppress the oscillations, is to express the first derivative 8D/%x at a point located halfway between the nodes, i.e. to write
a@ -= dX
@,h
@J_,
(2.4)
’
where J = I if u > 0 and J = I + 1 if u < 0. This is precisely what would be achieved by a full-upwind diff erencing of the convective term. Equation (2.4) still represents a central difference approximation, but at a shifted point, and for consistency, the source term and the diffusion term in (2.1) should also be evaluated at this shifted point. The above elementary considerations indicate that the construction of consistent finite element schemes of an upwind-type requires a coherent shift of all terms in the convectiondiffusion equation. The shift will here be achieved by adding a scalar multiple of the derivatives of (2.1) to (2.4) since the partial derivatives of the governing differential equation with respect to the space variables appear to represent a useful ingredient for the formulation of consistent finite element schemes of upwind-type. With this in mind, and returning to the one-dimensional convection-diffusion equation (2.1) we suggest replacing it by a generalized equation of the form -Ug+K
d2#
U-g--K~-dX
a3cP as
(2.5)
where, for simplicity, u and K have been assumed to be constant; the results can also be generalized for variable u and K. Equation (2.5) is obtained by subtracting from the original differential equation its derivative with respect to the space variable, x, multiplied by a function h(x). Equation (2.5) may be rewritten as
-u$(@-*y-)+K$(@-hy+s-A~=O, ax dX
(2.6)
28
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
and in this form, the generalized equation appears to be ideally suited for achieving the coherent shift of all terms required in a consistent finite element formulation of the upwindtYPee To obtain a Galerkin finite element choose a weight (test) function
approximation
to the generalized
equation
(2.3, we
w(x) = WIN,(X)
(2.7)
7
where upper case subscripts denote nodal values and NI are appropriate shape functions. The element equations are then constructed by representing the approximate solution @ by shape functions in a form identical to (2.7). Here, we choose the same shape functions used for the weight (test) function, i.e.,
where QI are the nodal values of the solution. following equation for problem (2.5):
The Galerkin
a*@ prp -u~+(~+Au)--&~~+s-/I$
method
dx=O,
then provides
all I,
the
(2.9)
where @ is required to satisfy (2.2a) and N,(O) = 0. At this point, it is important to note that, due to the presence of the third derivative in (2.9), problems do arise unless local interpolants with C’ continuity are introduced. These are, of course, very inconvenient in practice and assumptions will therefore be made that enable us to employ finite elements with Co continuity for the spatial discretization. To this aim, the third derivative term in (2.9) is first rewritten as L
I
NIA~ $
dx = ,:~(N,~K~)dx-~o=~(hK~)dX.
(2.10)
0
Applying the Green-Gauss theorem and requiring that the parameter A vanish at all the element boundaries, the first term in the right-hand side of (2.10) may be dropped, We then assume that the second term does not contain interelement contributions and rewrite (2.10) as L N~K
$+$dX.
(2.11)
By further assuming that the generalized flux (K + hu)d@/dx is continuous across the interelement boundaries, the integral in (2.9) may be represented as a sum of element contributions and the generalized Galerkin equation (2.9) may be transformed into
+
NIfh%)sdx}+N&)f=O,
where f is the prescribed flux defined in (2.2b).
(2.12)
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
29
To arrive at (2.12) we have used the Green-Gauss theorem and required the free parameter A to vanish at the boundaries x = 0 and x = L of the domain. Equation (2.12) represents the basis for deriving consistent one-dimensional finite element schemes of an upwind-type for the convection-diffusion equation (2.1). We note that the modified weight function, N1 + AaN1/Jx, acts on both the convective term and the source term. The diffusion term has also been modified, in the sense that, in addition to the standard Galerkin term, an additional term involving A#@/&’ now comes into play. Finally, the usual natural boundary condition on the diffusive flux is maintained. REMARK 2.1. The present generalized Galerkin formulation (i.e. (2.12)) is similar to the Petrov-Galerkin formulation of Hughes and Brooks [S]. In particular, the assumptions made in the treatment of the diffusion term to permit the use of Co finite elements are identical. The principal merit of the approach presented herein is that it automatically leads to a modified weight function which acts on both the convective term and the source term. It will also become apparent from the development presented in Section 6 that in multi-dimensional situations, the present approach naturally introduces the tensorial structure of the ‘upwind’ modification which was advocated by Hughes and Brooks [7,8] in their streamline-upwind Petrov-Galerkin formulation to avoid spurious crosswind diffusion effects. REMARK 2.2. When piecewise linear elements are employed for the spatial discretization, &D/&* is zero within an element, and therefore the ‘upwind’ modification does not affect the discretization of the diffusion term. This is, however, not the case for higher-order elements, as will be seen in Section 3, where a piecewise quadratic approximation will be developed. REMARK 2.3. The parameter A is a free parameter which determines the amount of ‘upwinding’. Its value should be determined case by case so as to maximize the accuracy. For the case of the linear, constant-coefficient convection-diffusion equation without a source term, the optimal choice for A in connection with a uniform mesh of piecewise linear elements is known to be [l] A =&h,
(2.13)
where a
=
coth21.-2 2
Y
(2.14)
and y = uhlu .
(2.15)
With the above choice, exact nodal answers are obtained for all values of the Peclet number y.
3. One-dimensional piecewise quadratic approximation As a first illustration of the generalized Galerkin formulation introduced in Section 2, let us examine the solution of the convection-diffusion equation (2.1) with quadratic shape function elements. We consider a generic element with the nodes 1 and 3 and a midpoint node 2 and a local coordinate system 5 in which the element goes from -1 to +l. The shape functions are then
30
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
N=G2-5),
N2 = 1 - ,$‘,
N3 = i(5’
+ 5) ,
(3.1)
with 5 = x/h, where h is the node-to-node distance. With the above parabolic shape functions, a2@/ax2 is constant within an element and the generalized Galerkin equation (2.12) indicates that, by contrast with the linear element, the upwind modification will now affect the weighting of the diffusion term. Consider a mesh of uniform elements, with nodes numbered in sequence and oddnumbered nodes between elements, even-numbered nodes at element midpoints. Two types of nodal equations are then derived from (2.12): at odd (interelement) nodes 21+ 1, the equations, in the absence of source term (s = 0), are -u
@ [(
21+2 -
2
@21
2h
@2,+3-
)- ( @2I+2
-
>I
@21-l
4h y+,
+
Qj21
@21+3-
2@21+1+
2
+
3hK
@
21+2 -
@21
h3
[(
@21+3-
)- (
@21-l
2h3
)I 7 =
@2,-l
4h2
0
(3.2)
while at even (midpoint) nodes 21 the equations are -U
@ ( 2h 21+1-
a-1
) + (K + *q21+1
- 2;f;-l+ @22-l) = 0 .
(3.3)
The last term in (3.2) represents the upwind modification of the diffusion term. No such modification appears in the midpoint equations (3.3), because JN2/%5 = -25 is an odd function of 5, so
aNzhK a2@
I-ax
%dx=O.
(3.4)
e
It remains now to determine the optimal value of the free parameter h. Actually two different optimal values of A are found: one for the mid-side node equations and another one for the corner-node equations. Since the structure of the mid-side node equations (3.3) is exactly that obtained with piecewise linear finite elements, it follows that the optimal A for such equations is that given by (2.13). The optimal value of A to be introduced in the corner-node equations can be determined by requiring that the truncation error introduced by the discrete equations (3.2) be identically zero. Adapting the procedure used by Heinrich and Envia 1141for linear elements, the optimal value of A for the corner-node equations is found to be given by A=ph,
(3.5)
where 4 2-coshy--+tanhz+y
1 sinh y (3.6)
P= 4 tanh $- sinh y - t sinh y tanh 5
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
Fig. 1. Values
of u and p as a function
of Peclet
31
number.
and y is defined in (2.15). The values of (Yin (2.14) and p in (3.6) are plotted in Fig. 1. With such values, exact nodal answers are obtained for all Peclet numbers in the particular case of the one-dimensional, linear advection-diff usion equation.
4. Numerical examples In this section, the performance of the generalized Galerkin schemes derived in Section 3 for the parabolic element in one space dimension will be illustrated in two test problems. The first test problem consists of solving (2.1) without a source term (s = 0) in the interval [0, l] with essential boundary conditions given by Q(O) = 1
and
The physical properties u = 100 )
@(l) = 0. are taken as
K=l,
(4.2)
corresponding to an element Peclet number (Pe = 2&/K) of 5 for a uniform parabolic elements. The exact solution for this boundary value problem is
CD(X) = 1-
(1
(4.1)
-
e"")/(l
- e”/“) ,
mesh of 20
(4.3)
and Fig. 2 indicates that the optimal quadratic finite element approximations derived in Section 3 effectively produce exact nodal answers. The solution was then repeated using the standard Galerkin formulation (A = 0) and the results displayed on Fig. 2 indicate an oscillatory behaviour of the numerical solution in the vicinity of the outflow boundary. As a second test problem, we have solved (2.1) with K = 0 (pure convection) in the presence of a source term. Leonard [15] has shown that some upwind formulations give very poor results when a source term is present. The test case used in [15] to illustrate the difficulty is shown here in Fig. 3. A bilinear spatial distribution of the source (si(x)) is assumed and the inflow boundary condition is Q(O) = 0. The problem was solved with a uniform mesh of 8
32
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
,900
.920
.940
.960
.980
i.
X Fig. 2. Solution
of the one-dimensional
numerical
example.
parabolic elements and Fig. 3 indicates that exact nodal results were obtained. was then repeated with a quadratic source term in the form
s*(x)= l-$x+&x2,
The problem
(4.4)
and, as indicated in Fig. 3, exact nodal results were again obtained.
5. Convergence properties of quadratic approximation In order to determine the effectiveness of the generalized Galerkin approach, it is necessary to examine the convergence properties of the method. In this study the quadratic approximation, as discussed in Section 3, and the linear approximation are compared. The linear element is equivalent to many of the current upwind methods, thus providing a standard by which the performance of the quadratic element may be measured. The convergence study was performed with K = 0.002 in (2.1) over the interval [0,11.The boundary conditions were G(O) = 0
)
Q(1) = 1 )
(5.1)
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
+
EXACT
33
2
X Fig. 3. Solutions
of a one-dimensional
problem
with the presence
of a source
term.
and the source term is s(x) = 4Ux3 - 12Kx* .
(5.2)
The velocity u was varied to obtain the element solution for this problem is CD(x)= x4.
Peclet numbers given in Fig. 4. The exact
(5.3)
Shown in Fig. 4 are the L2 errors for various Peclet numbers and various mesh sizes with the corresponding rates of convergence summarized in Table 1. The error E is calculated by
E* =
I
o= (cl& - QQ2
dx
,
(5.4)
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
34
LINEAR
QUADRATIC
-1
-2
-3 -2
W
W
5
g-4 -1
-3
-5
-4 -2
IO
-1.500
-J3-
GAM=lO
-I+
GAU=25
+
GAU=50
-1.000
LOG
-.
0
-6 -2
-_I 10
Fig. 4. L2 rates of convergence
CAM=10 GAU=25
+
CAM=50
-1.000
-1.500
LOG
h
-f3-
-&-
-.
JO
h
for linear and quadratic shape function elements.
where QE, aa are the exact and approximate solutions, respectively. For the linear element, the rate of convergence is reduced markedly when advection is introduced into the equation. This phenomenon is also exhibited by the quadratic approximation, but to a lesser extent than in the linear case. It appears that the rate of convergence of the quadratic example approaches the value of two as the problem becomes one of pure advection.
6. Extension to multi-dimensional problems In this section, the generalized Galerkin formulation introduced dimensional problems will be extended to deal with multi-dimensional Table 1 LI rates of convergence for linear and quadratic element at various Peclet numbers, y y = ah/K
0 10 25 50
Linear element 1.998 0.950 0.940 0.937
Quadratic element 2.997 2.222 2.161 2.148
in Section 2 for onesituations.
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
35
Consider the following problems in two or three dimensions:
@=g
onr,;
KZL/r dn
on rN .
Here o is the convective velocity, K is a positive diffusion coefficient and s a given source. TD represents the Dirichlet part of the boundary r of 0 and rN = r - Tr, the Neumann part. Proceeding as in Section 2 for the one-dimensional case, we associate with (6.1) a generalized governing equation which is obtained by subtracting from (6.1) the scalar product of its gradient by a vector
A = (Ax,A,, AL)
(6.3)
of free parameters associated with the coordinate The resulting generalized equation reads
directions.
-u.V@+A*V(~*V~)+V*(KV@)-A.V(V*(KV@))+S-A*Vs=O.
To obtain a Galerkin finite element choose a weight function
and represent
the approximate
approximation
to the generalized
(6.4)
equation
(6.4) we
solution by
CD(x)= @INI(x) .
(6.6)
In (6.5) and (6.6) the shape functions N1 are assumed to be continuous across interelement boundaries and the weight function w is assumed to vanish on the Dirichlet boundary r,,. The Galerkin method then provides the following equation for problem (6.4):
N,h dr,
= 0,
(6.7)
where fil = N1 + A . VNI
(6.8)
are the generalized, discontinuous weighting functions. In deriving (6.7) we have made assumptions equivalent to those used in Section 2 to treat the third spatial derivatives arising from the generalized diffusion term in (6.4). In addition, we have required that the normal component of the vector A of free parameters vanish on the boundary r of the integration domain.
T. Belytschko et al., Generalized Galerkin method for convection-difision
36
problems
Denoting by (u, v,, w) the Cartesian components of the convective velocity ~1,one notes from the first term in (6.7) that the generalized Galerkin method introduces an added diffusivity t? which processes the following tensorial structure:
It follows that, if A,, A, and A, are each evaluated according to the value of the velocity component along the associated coordinate direction, then the added diffusivity will act only in the direction of the flow. For example, if the flow field is such that uf0,
VfO
and
w = 0,
(6.10)
then A, = 0, and the diffusivity matrix (6.9) reduces to
(6.11)
which clearly manifests the absence of crosswind diffusion. It is therefore concluded that the present generalized Galerkin formulation possesses the same attributes as the streamline-upwind Petrov-Galerkin formulation of Hughes and Brooks [7,8]. The proposed method will now be illustrated for the case of biquadratic shape function elements in two dimensions.
7. Two-dimensional quadratic approximation Let us examine the solution of (6.1) with biquadratic shape function elements in two dimensions. We consider a generic 9-node element with nodes numbered as indicated in Fig. 5.
0 8
l 9
8
1
5 -
2
0
Fig. 5. Normalized
9-node
..
0
quadrilateral
element.
T. Belytschko
et al., Generalized Galerkin method for convection-diffusion problems
31
In terms of the normalized coordinates (5, 77) the shape functions for this element are given in Table 2 [16]. The generalized weighting functions in (6.8) are here in the form
(7.1) and, according to the analysis made in Section 3 for the quadratic element in one dimension, we associate two optimal values with each of the free parameters A, and A,: A”,and h’y are for ‘corner-node’ type shape functions and these parameters are evaluated by (3.5)-(3.6) with y defined by the velocity component and the mesh size in the relevant coordinate direction. A? and A? are associated with ‘midside-node’ type shape functions and evaluated according to (2.12)-(2.13). Under these conditions, the generalized weighting functions for the 9-node biquadratic elements are
(7.2)
To assess the effectiveness of the above generalized series of numerical calculations were performed.
biquadratic
finite element schemes, a
EXAMPLE 7.1. Advection skew to the mesh [7-151. The problem is depicted in Fig. 6. The flow is unidirectional, constant (]]u]]= 1) and skewed to the mesh. The diffusivity coefficient K was Table 2 Z
NI
ftx - lh(v
- 1)
a(5 + lh(77 - 1) :4%+ lb-/(.11 + 1) !5(5- lh(11 + 1) !(l - 5’)(11- 1) M5 + 110- 7’) :u - if2)q(q+ 1)
:stc- Ml
- v2) (1 - 5’)(1- T2)
T. Belytschko et al., Generalized Galerkin method for convection-difision
38
\
problems
/
@=l Fig. 6. Advection
skew of mesh: problem statement.
taken to be zero, resulting in a pure advection problem. A five-by-five mesh of equal-sized biquadratic elements was employed. The upwind boundary condition is discontinuous, as shown in Fig. 6, and downwind homogeneous natural boundary conditions are considered. A comparison of downwind boundary profiles of CDis shown in Fig. 7 for flow directions given by 8 = 21.8”, 45.0”, and 69.4”. The results clearly indicate that, except for 8 = 69.4”, the present formulation produces less oscillations than the Galerkin formulation, without introducing spurious cross diffusion effects. EXAMPLE 7.2. Axi-symmetric example. To study the effect of this upwinding technique in another multi-dimensional situation, an exact solution for the axi-symmetric advection-diffusion equation S~ra6, was constructed and compared to a two-dimensional finite element solution. dimensional mesh is shown in Fig. 8. The material properties and solution are K =
1
x lo-6,
The two-
(7.4)
u, = 50/r, G(r) = 1.259 x 10-3r4 + 1.704 x 10-3r3 ,
(7.3)
(7.5) 5 S r s 6,
(7.6)
39
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
2 1.8
DEGREES
1.200 (
LEGEND -
EXACT
-f3-
PRESENT
-r+
GALERKIN
.800 .600
FORMULATION
X 45.0
DEGREES
,800
.600
,200
.ooo
Fig. 7. (a) Downwind
+I--
PRESENT
tr
GALERKIN
boundary
FORMUIATION
profile of @ for 8 = 21.8”. (b)
so the velocity satisfies continuity specified at the inside and outside the stated conditions is
Downwind boundary
profile
of @ for 8 =
for incompressible flow. Dirichlet boundary conditions radii. The source term satisfying the governing equation
s(r) = 0.252r2 + 0.256r. Shown
in Fig. 9 are the results
45.0”.
are for
(7.7) of the standard
and upwinded
solutions.
Oscillations
are
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
40
69.4
1.200
DEGREES
BOO -.600--
.400--
.200--
.000--
-.200-l .ooo
LEGEND -
EXACT
-t9-
PRESENT
d-
GALERKIN
FORMULATION
300
1.000
1.500
2.1
X Fig. 7(c). Downwind
apparent in the standard formulation to the exact solution.
boundary
profile of @ for @= 69.4”.
while the upwinded form provides a good approximation
8. Conclusions A new method has been presented
5
for the derivation of the streamline upwind formulation
6
Fig. 8. Mesh used in axi-symmetric
example.
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
WITHOUT
41
UPWINDING
2Qoo
-200
-600
LEGEND
-800
-10001 5.000
5.200
5.400
5.600
5.600
-E3-
X=Y
-e-
X,Y AXIS
+-
EXACT
6.’
R
Fig. 9(a). The solution
of the axi-symmetric
example
using the standard
Galerkin
method.
proposed by Hughes and coworkers for the steady convection-diffusion equation. Essentially, the proposed method applies the conventional Galerkin procedure to a generalized governing equation which is obtained by subtracting from the original differential equation the scalar product of its gradient and a vector of free parameters associated with each of the coordinate directions. The principal merit of the proposed approach is that it automatically leads to a modified weight function which acts on all terms in the governing equation, and that, in multi-dimensional situations, it naturally introduces an upwind modification possessing the tensorial structure required to avoid crosswind diffusion. The effectiveness of the method has been illustrated for the case of quadratic local interpolations in one and two space dimensions. It was also shown that the quadratic shape function element exhibits a definitely superior convergence rate than the linear element over the whole range of Peclet numbers. In particular, while the convergence of the linear element drops from order h* to order h as the velocity increases, the convergence of the quadratic element only drops from order h3 to order h’.
42
T. Belytschko et al., Generalized Galerkin method for convection-difision
UPWINDED
problems
SOLUTION
LEGEND
5.000
5.200
5.400
5.600
5.800
6.
-a-
X=Y
-a-
X,Y AXIS
+t-
EXACT
0
R Fig. 9(b). The generalized
Galerkin
solution
of the axi-symmetric
example
References [l] D. Allen and R. Southwell, Relaxation methods applied to determine the motion, in two dimensions, of a viscous fluid past a fixed cylinder, Quart. J. Mech. Appl. Math. VIII (1955) 129-145. [2] D.B. Spalding, A novel finite difference formulation for differential expressions involving both first and second derivatives, Internat. J. Numer. Meths. Engrg. 4 (1972) 551-559. [3] I. Christie, D.F. Griffiths, A.R. Mitchell and O.C. Zienkiewicz, Finite element methods for second order differential equations with significant first derivatives, Internat. J. Numer. Meths. Engrg. 10 (1976) 1389-1396. [4] T.J.R. Hughes, A simple scheme for developing ‘upwind’ finite elements, Internat. J. Numer. Meths. Engrg. 12 (1978) 1359-1365. [5] T. Belytschko and J.M. Kennedy, Computer models for subassembly simulation, Nuclear Engrg. Design 43 (1978) 17-38. [6] T. Belytschko and I. Eldib, Analysis of a finite element upwind scheme, in: T.J.R. Hughes, ed., Finite Element Methods for Convection Dominated Flows, AMD Vol. 34 (ASME, New York, 1979). [7] T.J.R. Hughes and A.N. Brooks, A multi-dimensional upwind scheme with no crosswind diffusion, in: T.J.R. Hughes, ed., Finite Element Methods for Convection Dominated Flows, AMD Vol. 34 (ASME, New York, 1979). [8] T.J.R. Hughes and A.N. Brooks, A theoretical framework for Petrov-Galerkin methods with discontinuous weighting functions: application to the streamline-upwind procedure, in: R.H. Gallagher, ed., Finite Elements in Fluids, Vol. 4 (Wiley, New York, 1982).
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
43
191 T.J.R. Hughes, T.E. Tezduyar and A.N. Brooks, A Petrov-Galerkin finite element formulation for systems of conservation laws with special reference to the compressible Euler equations, Proc. IMA Conf. on Numerical Methods in Fluid Dynamics, University of Reading, England (Academic Press, New York, 1982). [lo] J.W. Barrett and K.W. Morton, Optimal finite element solutions to diffusion-convection problems in one dimension, Internat. J. Numer. Meths. Engrg. 15 (1980) 1457-1474. [ll] J.W. Barrett and K.W. Morton, Optimal Petrov-Galerkin methods through approximate symmetrization, IMA J. Numer. Anal. 1 (1981) 439-468. [12] K.W. Morton, Generalized Galerkin methods for steady and unsteady problems, Proc. IMA Conf. on Numerical Methods in Fluid Dynamics, University of Reading, England 1982 (Academic Press, New York, 1982). [13] J. Donea, A Taylor-Galerkin method for convective transport problems, Internat. J. Numer. Meths. Engrg. 20 (1984) 101-119. [14] J.C. Heinrich and E. Envia, Finite element techniques in transport phenomena, Preprint. [IS] B.P. Leonard, A survey of finite differences of opinion on numerical muddling of the incomprehensible Dominated defective confusion equation, in: T.J.R. Hughes, ed., Finite Element Methods for Convection Flow, AMD Vol. 34 (ASME, New York, 1979). [16] O.C. Zienkiewicz, The Finite Element Method (McGraw-Hill, London, 3rd ed., 1977).
IN APPLIED
MECHANICS
AND ENGINEERING
48 (1985) 25-43
A GENERALIZED GALERKIN METHOD FOR STEADY CONVECTION-DIFFUSION PROBLEMS WITH APPLICATION TO QUADRATIC SHAPE FUNCTION ELEMENTS
Jean DONEA Applied Mechanics Division, Joint Research Centre-Ispra Establishment, Communities, 27020 Ispra, Varese, Italy
Ted BELYTSCHKO
Commission of the European
and Patrick SMOLINSKI
Department of Civil Engineering, Northwestern University, Evanston, IL 60201, U.S.A.
Revised
Received manuscript
12 July 1983 received 14 July 1984
A generalization of the standard Galerkin finite element method is considered to enable it to deal successfully with steady convection-diffusion problems. The proposed method employs a generalized governing equation which is obtained by subtracting from the original differential equation the scalar product of its gradient by a vector of free parameters associated with each of the coordinate directions. The generalized equation is successively discretized by the standard Bubnov-Galerkin finite element method. The effectiveness of the method is illustrated for the case of quadratic local interpolations in one and two space dimensions.
1. Introduction
Traditionally, the steady convection-diffusion equation has been one of the most difficult equations to solve by numerical methods. The difficulties arise from the unsymmetric nature of the convection operator, which tends to engender physically unrealistic oscillations in the numerical solutions when the transport terms are dominant. This problem arises in both Galerkin/finite element and finite difference formulations with reasonable grid size. For example, when second-order central difference approximations to the convective term are generated, the discrete equations at odd and even nodes tend to become decoupled, leading to an almost singular system of global equations. Several schemes have been devised to effectively deal with these difficulties. The remedy with difference methods has long been to use upwind differencing for the convective term. This was accomplished by shifting the first derivative term in the upwind direction [l, 21. For the simple case of the one-dimensional, constant-coefficient, convection-diffusion equation, an upwind scheme could be devised [l] which gives exact nodal answers. However, problems of accuracy still remained in practical situations through excessive crosswind diffusion and with variable flow fields and source terms. The upwind concept was first introduced in a finite element context by Christie et al. [3], who achieved an upwind character in the finite element equations by skewing the weight functions appropriately. Hughes [4] has achieved similar 00457825/85/$3.30
@ 1985, Elsevier
Science
Publishers
B.V. (North-Holland)
26
T. Belytschko et al., Generalized Galerkin method for convection-difision
problems
results for low-order elements by shifting the numerical integration points in the evaluation of the transport terms. Computationally simple amplification schemes for achieving upwind-type finite element approximations have been proposed by Belytschko and Kennedy [5] and analyzed by Belytschko and Eldib [6]. These early finite element procedures of upwind-type were also subject to the same criticisms as upwind finite differences when applied to multi-dimensional situations or in problems with variable source terms. More recently, however, methods were developed which do not suffer these deficiencies and represent a computation of convection-dominated significant breakthrough in the finite element phenomena. These are the so-called consistent Petrov-Galerkin methods in which a suitably modified weight function is applied to all terms in the equations. Of particular relevance in this context are the papers of Hughes and coworkers (see [7-9]), and of Morton and coworkers (see [lo-121). In the present paper, we propose an alternative method for deriving accurate finite element schemes for the steady convection-diffusion equation which has enabled us to treat higherorder shape function. The proposed approach is inspired by the Taylor-Galerkin method recently developed by Donea [13] for first-order hyperbolic equations. It consists of constructing a generalized governing equation for the convection-diffusion problem which includes all the necessary ingredients to derive upwind-type finite element schemes in a consistent fashion. This generalized equation is then discretized by application of the standard Bubnov-Galerkin process. The proposed procedure is introduced in Section 2 for the simple case of the linear convection-diffusion equation in one space dimension. The resulting scheme for piecewise quadratic local interpolations is presented in Section 3 and shown to be exact under appropriate circumstances. A series of comparison problems is presented in Section 4. Section 5 is devoted to an analysis of the convergence properties of the parabolic element in comparison with those of the linear element, indicating a definitely superior convergence rate for the quadratic shape function element. Finally, the extension of the proposed method to deal with multi-dimensional convection-diffusion problems is discussed in Section 6 and illustrated in Section 7 for the case of biquadratic shape function elements.
2. Generalized
governing
equation for the convection-diffusion
For simplicity, consider first the linear, constant-coefficient one space dimension,
a@
-Udx+Kz
a2@+s=o, OSXSL,
problem
convection-diffusion
equation in
(2.1)
with the following boundary conditions: Q(O) = g 7
(2.2a) (2.2b)
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
27
Here u is the convective velocity, K is a positive diffusion coefficient and s a given source. As is well-known, for large values of the non-dimensional Peclet number uL/K, the application of the Galerkin method to (2.1) often leads to wildly oscillatory solutions. In effect, central difference approximations to the convective term are generated and the resulting system of discrete equations becomes almost singular. This may readily be seen by considering a uniform mesh of linear elements with nodes numbered sequentially. For K = 0 in (2.1) the Galerkin equation at an interior node I in a uniform mesh of elements with length h is - & (@I+, - @I-,> + SI + b(s,-I-
2s, + s,+1) = 0 .
(2.3)
It can immediately be noted that the equations at odd-numbered nodes are fully decoupled from the equations at even-numbered nodes, and vice-versa. This is one explanation for why oscillatory solutions are obtained when solving (2.1) for a high value of the Peclet number. The easiest way to couple the discrete equations (2.3) and thereby suppress the oscillations, is to express the first derivative 8D/%x at a point located halfway between the nodes, i.e. to write
a@ -= dX
@,h
@J_,
(2.4)
’
where J = I if u > 0 and J = I + 1 if u < 0. This is precisely what would be achieved by a full-upwind diff erencing of the convective term. Equation (2.4) still represents a central difference approximation, but at a shifted point, and for consistency, the source term and the diffusion term in (2.1) should also be evaluated at this shifted point. The above elementary considerations indicate that the construction of consistent finite element schemes of an upwind-type requires a coherent shift of all terms in the convectiondiffusion equation. The shift will here be achieved by adding a scalar multiple of the derivatives of (2.1) to (2.4) since the partial derivatives of the governing differential equation with respect to the space variables appear to represent a useful ingredient for the formulation of consistent finite element schemes of upwind-type. With this in mind, and returning to the one-dimensional convection-diffusion equation (2.1) we suggest replacing it by a generalized equation of the form -Ug+K
d2#
U-g--K~-dX
a3cP as
(2.5)
where, for simplicity, u and K have been assumed to be constant; the results can also be generalized for variable u and K. Equation (2.5) is obtained by subtracting from the original differential equation its derivative with respect to the space variable, x, multiplied by a function h(x). Equation (2.5) may be rewritten as
-u$(@-*y-)+K$(@-hy+s-A~=O, ax dX
(2.6)
28
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
and in this form, the generalized equation appears to be ideally suited for achieving the coherent shift of all terms required in a consistent finite element formulation of the upwindtYPee To obtain a Galerkin finite element choose a weight (test) function
approximation
to the generalized
equation
(2.3, we
w(x) = WIN,(X)
(2.7)
7
where upper case subscripts denote nodal values and NI are appropriate shape functions. The element equations are then constructed by representing the approximate solution @ by shape functions in a form identical to (2.7). Here, we choose the same shape functions used for the weight (test) function, i.e.,
where QI are the nodal values of the solution. following equation for problem (2.5):
The Galerkin
a*@ prp -u~+(~+Au)--&~~+s-/I$
method
dx=O,
then provides
all I,
the
(2.9)
where @ is required to satisfy (2.2a) and N,(O) = 0. At this point, it is important to note that, due to the presence of the third derivative in (2.9), problems do arise unless local interpolants with C’ continuity are introduced. These are, of course, very inconvenient in practice and assumptions will therefore be made that enable us to employ finite elements with Co continuity for the spatial discretization. To this aim, the third derivative term in (2.9) is first rewritten as L
I
NIA~ $
dx = ,:~(N,~K~)dx-~o=~(hK~)dX.
(2.10)
0
Applying the Green-Gauss theorem and requiring that the parameter A vanish at all the element boundaries, the first term in the right-hand side of (2.10) may be dropped, We then assume that the second term does not contain interelement contributions and rewrite (2.10) as L N~K
$+$dX.
(2.11)
By further assuming that the generalized flux (K + hu)d@/dx is continuous across the interelement boundaries, the integral in (2.9) may be represented as a sum of element contributions and the generalized Galerkin equation (2.9) may be transformed into
+
NIfh%)sdx}+N&)f=O,
where f is the prescribed flux defined in (2.2b).
(2.12)
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
29
To arrive at (2.12) we have used the Green-Gauss theorem and required the free parameter A to vanish at the boundaries x = 0 and x = L of the domain. Equation (2.12) represents the basis for deriving consistent one-dimensional finite element schemes of an upwind-type for the convection-diffusion equation (2.1). We note that the modified weight function, N1 + AaN1/Jx, acts on both the convective term and the source term. The diffusion term has also been modified, in the sense that, in addition to the standard Galerkin term, an additional term involving A#@/&’ now comes into play. Finally, the usual natural boundary condition on the diffusive flux is maintained. REMARK 2.1. The present generalized Galerkin formulation (i.e. (2.12)) is similar to the Petrov-Galerkin formulation of Hughes and Brooks [S]. In particular, the assumptions made in the treatment of the diffusion term to permit the use of Co finite elements are identical. The principal merit of the approach presented herein is that it automatically leads to a modified weight function which acts on both the convective term and the source term. It will also become apparent from the development presented in Section 6 that in multi-dimensional situations, the present approach naturally introduces the tensorial structure of the ‘upwind’ modification which was advocated by Hughes and Brooks [7,8] in their streamline-upwind Petrov-Galerkin formulation to avoid spurious crosswind diffusion effects. REMARK 2.2. When piecewise linear elements are employed for the spatial discretization, &D/&* is zero within an element, and therefore the ‘upwind’ modification does not affect the discretization of the diffusion term. This is, however, not the case for higher-order elements, as will be seen in Section 3, where a piecewise quadratic approximation will be developed. REMARK 2.3. The parameter A is a free parameter which determines the amount of ‘upwinding’. Its value should be determined case by case so as to maximize the accuracy. For the case of the linear, constant-coefficient convection-diffusion equation without a source term, the optimal choice for A in connection with a uniform mesh of piecewise linear elements is known to be [l] A =&h,
(2.13)
where a
=
coth21.-2 2
Y
(2.14)
and y = uhlu .
(2.15)
With the above choice, exact nodal answers are obtained for all values of the Peclet number y.
3. One-dimensional piecewise quadratic approximation As a first illustration of the generalized Galerkin formulation introduced in Section 2, let us examine the solution of the convection-diffusion equation (2.1) with quadratic shape function elements. We consider a generic element with the nodes 1 and 3 and a midpoint node 2 and a local coordinate system 5 in which the element goes from -1 to +l. The shape functions are then
30
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
N=G2-5),
N2 = 1 - ,$‘,
N3 = i(5’
+ 5) ,
(3.1)
with 5 = x/h, where h is the node-to-node distance. With the above parabolic shape functions, a2@/ax2 is constant within an element and the generalized Galerkin equation (2.12) indicates that, by contrast with the linear element, the upwind modification will now affect the weighting of the diffusion term. Consider a mesh of uniform elements, with nodes numbered in sequence and oddnumbered nodes between elements, even-numbered nodes at element midpoints. Two types of nodal equations are then derived from (2.12): at odd (interelement) nodes 21+ 1, the equations, in the absence of source term (s = 0), are -u
@ [(
21+2 -
2
@21
2h
@2,+3-
)- ( @2I+2
-
>I
@21-l
4h y+,
+
Qj21
@21+3-
2@21+1+
2
+
3hK
@
21+2 -
@21
h3
[(
@21+3-
)- (
@21-l
2h3
)I 7 =
@2,-l
4h2
0
(3.2)
while at even (midpoint) nodes 21 the equations are -U
@ ( 2h 21+1-
a-1
) + (K + *q21+1
- 2;f;-l+ @22-l) = 0 .
(3.3)
The last term in (3.2) represents the upwind modification of the diffusion term. No such modification appears in the midpoint equations (3.3), because JN2/%5 = -25 is an odd function of 5, so
aNzhK a2@
I-ax
%dx=O.
(3.4)
e
It remains now to determine the optimal value of the free parameter h. Actually two different optimal values of A are found: one for the mid-side node equations and another one for the corner-node equations. Since the structure of the mid-side node equations (3.3) is exactly that obtained with piecewise linear finite elements, it follows that the optimal A for such equations is that given by (2.13). The optimal value of A to be introduced in the corner-node equations can be determined by requiring that the truncation error introduced by the discrete equations (3.2) be identically zero. Adapting the procedure used by Heinrich and Envia 1141for linear elements, the optimal value of A for the corner-node equations is found to be given by A=ph,
(3.5)
where 4 2-coshy--+tanhz+y
1 sinh y (3.6)
P= 4 tanh $- sinh y - t sinh y tanh 5
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
Fig. 1. Values
of u and p as a function
of Peclet
31
number.
and y is defined in (2.15). The values of (Yin (2.14) and p in (3.6) are plotted in Fig. 1. With such values, exact nodal answers are obtained for all Peclet numbers in the particular case of the one-dimensional, linear advection-diff usion equation.
4. Numerical examples In this section, the performance of the generalized Galerkin schemes derived in Section 3 for the parabolic element in one space dimension will be illustrated in two test problems. The first test problem consists of solving (2.1) without a source term (s = 0) in the interval [0, l] with essential boundary conditions given by Q(O) = 1
and
The physical properties u = 100 )
@(l) = 0. are taken as
K=l,
(4.2)
corresponding to an element Peclet number (Pe = 2&/K) of 5 for a uniform parabolic elements. The exact solution for this boundary value problem is
CD(X) = 1-
(1
(4.1)
-
e"")/(l
- e”/“) ,
mesh of 20
(4.3)
and Fig. 2 indicates that the optimal quadratic finite element approximations derived in Section 3 effectively produce exact nodal answers. The solution was then repeated using the standard Galerkin formulation (A = 0) and the results displayed on Fig. 2 indicate an oscillatory behaviour of the numerical solution in the vicinity of the outflow boundary. As a second test problem, we have solved (2.1) with K = 0 (pure convection) in the presence of a source term. Leonard [15] has shown that some upwind formulations give very poor results when a source term is present. The test case used in [15] to illustrate the difficulty is shown here in Fig. 3. A bilinear spatial distribution of the source (si(x)) is assumed and the inflow boundary condition is Q(O) = 0. The problem was solved with a uniform mesh of 8
32
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
,900
.920
.940
.960
.980
i.
X Fig. 2. Solution
of the one-dimensional
numerical
example.
parabolic elements and Fig. 3 indicates that exact nodal results were obtained. was then repeated with a quadratic source term in the form
s*(x)= l-$x+&x2,
The problem
(4.4)
and, as indicated in Fig. 3, exact nodal results were again obtained.
5. Convergence properties of quadratic approximation In order to determine the effectiveness of the generalized Galerkin approach, it is necessary to examine the convergence properties of the method. In this study the quadratic approximation, as discussed in Section 3, and the linear approximation are compared. The linear element is equivalent to many of the current upwind methods, thus providing a standard by which the performance of the quadratic element may be measured. The convergence study was performed with K = 0.002 in (2.1) over the interval [0,11.The boundary conditions were G(O) = 0
)
Q(1) = 1 )
(5.1)
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
+
EXACT
33
2
X Fig. 3. Solutions
of a one-dimensional
problem
with the presence
of a source
term.
and the source term is s(x) = 4Ux3 - 12Kx* .
(5.2)
The velocity u was varied to obtain the element solution for this problem is CD(x)= x4.
Peclet numbers given in Fig. 4. The exact
(5.3)
Shown in Fig. 4 are the L2 errors for various Peclet numbers and various mesh sizes with the corresponding rates of convergence summarized in Table 1. The error E is calculated by
E* =
I
o= (cl& - QQ2
dx
,
(5.4)
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
34
LINEAR
QUADRATIC
-1
-2
-3 -2
W
W
5
g-4 -1
-3
-5
-4 -2
IO
-1.500
-J3-
GAM=lO
-I+
GAU=25
+
GAU=50
-1.000
LOG
-.
0
-6 -2
-_I 10
Fig. 4. L2 rates of convergence
CAM=10 GAU=25
+
CAM=50
-1.000
-1.500
LOG
h
-f3-
-&-
-.
JO
h
for linear and quadratic shape function elements.
where QE, aa are the exact and approximate solutions, respectively. For the linear element, the rate of convergence is reduced markedly when advection is introduced into the equation. This phenomenon is also exhibited by the quadratic approximation, but to a lesser extent than in the linear case. It appears that the rate of convergence of the quadratic example approaches the value of two as the problem becomes one of pure advection.
6. Extension to multi-dimensional problems In this section, the generalized Galerkin formulation introduced dimensional problems will be extended to deal with multi-dimensional Table 1 LI rates of convergence for linear and quadratic element at various Peclet numbers, y y = ah/K
0 10 25 50
Linear element 1.998 0.950 0.940 0.937
Quadratic element 2.997 2.222 2.161 2.148
in Section 2 for onesituations.
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
35
Consider the following problems in two or three dimensions:
@=g
onr,;
KZL/r dn
on rN .
Here o is the convective velocity, K is a positive diffusion coefficient and s a given source. TD represents the Dirichlet part of the boundary r of 0 and rN = r - Tr, the Neumann part. Proceeding as in Section 2 for the one-dimensional case, we associate with (6.1) a generalized governing equation which is obtained by subtracting from (6.1) the scalar product of its gradient by a vector
A = (Ax,A,, AL)
(6.3)
of free parameters associated with the coordinate The resulting generalized equation reads
directions.
-u.V@+A*V(~*V~)+V*(KV@)-A.V(V*(KV@))+S-A*Vs=O.
To obtain a Galerkin finite element choose a weight function
and represent
the approximate
approximation
to the generalized
(6.4)
equation
(6.4) we
solution by
CD(x)= @INI(x) .
(6.6)
In (6.5) and (6.6) the shape functions N1 are assumed to be continuous across interelement boundaries and the weight function w is assumed to vanish on the Dirichlet boundary r,,. The Galerkin method then provides the following equation for problem (6.4):
N,h dr,
= 0,
(6.7)
where fil = N1 + A . VNI
(6.8)
are the generalized, discontinuous weighting functions. In deriving (6.7) we have made assumptions equivalent to those used in Section 2 to treat the third spatial derivatives arising from the generalized diffusion term in (6.4). In addition, we have required that the normal component of the vector A of free parameters vanish on the boundary r of the integration domain.
T. Belytschko et al., Generalized Galerkin method for convection-difision
36
problems
Denoting by (u, v,, w) the Cartesian components of the convective velocity ~1,one notes from the first term in (6.7) that the generalized Galerkin method introduces an added diffusivity t? which processes the following tensorial structure:
It follows that, if A,, A, and A, are each evaluated according to the value of the velocity component along the associated coordinate direction, then the added diffusivity will act only in the direction of the flow. For example, if the flow field is such that uf0,
VfO
and
w = 0,
(6.10)
then A, = 0, and the diffusivity matrix (6.9) reduces to
(6.11)
which clearly manifests the absence of crosswind diffusion. It is therefore concluded that the present generalized Galerkin formulation possesses the same attributes as the streamline-upwind Petrov-Galerkin formulation of Hughes and Brooks [7,8]. The proposed method will now be illustrated for the case of biquadratic shape function elements in two dimensions.
7. Two-dimensional quadratic approximation Let us examine the solution of (6.1) with biquadratic shape function elements in two dimensions. We consider a generic 9-node element with nodes numbered as indicated in Fig. 5.
0 8
l 9
8
1
5 -
2
0
Fig. 5. Normalized
9-node
..
0
quadrilateral
element.
T. Belytschko
et al., Generalized Galerkin method for convection-diffusion problems
31
In terms of the normalized coordinates (5, 77) the shape functions for this element are given in Table 2 [16]. The generalized weighting functions in (6.8) are here in the form
(7.1) and, according to the analysis made in Section 3 for the quadratic element in one dimension, we associate two optimal values with each of the free parameters A, and A,: A”,and h’y are for ‘corner-node’ type shape functions and these parameters are evaluated by (3.5)-(3.6) with y defined by the velocity component and the mesh size in the relevant coordinate direction. A? and A? are associated with ‘midside-node’ type shape functions and evaluated according to (2.12)-(2.13). Under these conditions, the generalized weighting functions for the 9-node biquadratic elements are
(7.2)
To assess the effectiveness of the above generalized series of numerical calculations were performed.
biquadratic
finite element schemes, a
EXAMPLE 7.1. Advection skew to the mesh [7-151. The problem is depicted in Fig. 6. The flow is unidirectional, constant (]]u]]= 1) and skewed to the mesh. The diffusivity coefficient K was Table 2 Z
NI
ftx - lh(v
- 1)
a(5 + lh(77 - 1) :4%+ lb-/(.11 + 1) !5(5- lh(11 + 1) !(l - 5’)(11- 1) M5 + 110- 7’) :u - if2)q(q+ 1)
:stc- Ml
- v2) (1 - 5’)(1- T2)
T. Belytschko et al., Generalized Galerkin method for convection-difision
38
\
problems
/
@=l Fig. 6. Advection
skew of mesh: problem statement.
taken to be zero, resulting in a pure advection problem. A five-by-five mesh of equal-sized biquadratic elements was employed. The upwind boundary condition is discontinuous, as shown in Fig. 6, and downwind homogeneous natural boundary conditions are considered. A comparison of downwind boundary profiles of CDis shown in Fig. 7 for flow directions given by 8 = 21.8”, 45.0”, and 69.4”. The results clearly indicate that, except for 8 = 69.4”, the present formulation produces less oscillations than the Galerkin formulation, without introducing spurious cross diffusion effects. EXAMPLE 7.2. Axi-symmetric example. To study the effect of this upwinding technique in another multi-dimensional situation, an exact solution for the axi-symmetric advection-diffusion equation S~ra6, was constructed and compared to a two-dimensional finite element solution. dimensional mesh is shown in Fig. 8. The material properties and solution are K =
1
x lo-6,
The two-
(7.4)
u, = 50/r, G(r) = 1.259 x 10-3r4 + 1.704 x 10-3r3 ,
(7.3)
(7.5) 5 S r s 6,
(7.6)
39
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
2 1.8
DEGREES
1.200 (
LEGEND -
EXACT
-f3-
PRESENT
-r+
GALERKIN
.800 .600
FORMULATION
X 45.0
DEGREES
,800
.600
,200
.ooo
Fig. 7. (a) Downwind
+I--
PRESENT
tr
GALERKIN
boundary
FORMUIATION
profile of @ for 8 = 21.8”. (b)
so the velocity satisfies continuity specified at the inside and outside the stated conditions is
Downwind boundary
profile
of @ for 8 =
for incompressible flow. Dirichlet boundary conditions radii. The source term satisfying the governing equation
s(r) = 0.252r2 + 0.256r. Shown
in Fig. 9 are the results
45.0”.
are for
(7.7) of the standard
and upwinded
solutions.
Oscillations
are
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
40
69.4
1.200
DEGREES
BOO -.600--
.400--
.200--
.000--
-.200-l .ooo
LEGEND -
EXACT
-t9-
PRESENT
d-
GALERKIN
FORMULATION
300
1.000
1.500
2.1
X Fig. 7(c). Downwind
apparent in the standard formulation to the exact solution.
boundary
profile of @ for @= 69.4”.
while the upwinded form provides a good approximation
8. Conclusions A new method has been presented
5
for the derivation of the streamline upwind formulation
6
Fig. 8. Mesh used in axi-symmetric
example.
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
WITHOUT
41
UPWINDING
2Qoo
-200
-600
LEGEND
-800
-10001 5.000
5.200
5.400
5.600
5.600
-E3-
X=Y
-e-
X,Y AXIS
+-
EXACT
6.’
R
Fig. 9(a). The solution
of the axi-symmetric
example
using the standard
Galerkin
method.
proposed by Hughes and coworkers for the steady convection-diffusion equation. Essentially, the proposed method applies the conventional Galerkin procedure to a generalized governing equation which is obtained by subtracting from the original differential equation the scalar product of its gradient and a vector of free parameters associated with each of the coordinate directions. The principal merit of the proposed approach is that it automatically leads to a modified weight function which acts on all terms in the governing equation, and that, in multi-dimensional situations, it naturally introduces an upwind modification possessing the tensorial structure required to avoid crosswind diffusion. The effectiveness of the method has been illustrated for the case of quadratic local interpolations in one and two space dimensions. It was also shown that the quadratic shape function element exhibits a definitely superior convergence rate than the linear element over the whole range of Peclet numbers. In particular, while the convergence of the linear element drops from order h* to order h as the velocity increases, the convergence of the quadratic element only drops from order h3 to order h’.
42
T. Belytschko et al., Generalized Galerkin method for convection-difision
UPWINDED
problems
SOLUTION
LEGEND
5.000
5.200
5.400
5.600
5.800
6.
-a-
X=Y
-a-
X,Y AXIS
+t-
EXACT
0
R Fig. 9(b). The generalized
Galerkin
solution
of the axi-symmetric
example
References [l] D. Allen and R. Southwell, Relaxation methods applied to determine the motion, in two dimensions, of a viscous fluid past a fixed cylinder, Quart. J. Mech. Appl. Math. VIII (1955) 129-145. [2] D.B. Spalding, A novel finite difference formulation for differential expressions involving both first and second derivatives, Internat. J. Numer. Meths. Engrg. 4 (1972) 551-559. [3] I. Christie, D.F. Griffiths, A.R. Mitchell and O.C. Zienkiewicz, Finite element methods for second order differential equations with significant first derivatives, Internat. J. Numer. Meths. Engrg. 10 (1976) 1389-1396. [4] T.J.R. Hughes, A simple scheme for developing ‘upwind’ finite elements, Internat. J. Numer. Meths. Engrg. 12 (1978) 1359-1365. [5] T. Belytschko and J.M. Kennedy, Computer models for subassembly simulation, Nuclear Engrg. Design 43 (1978) 17-38. [6] T. Belytschko and I. Eldib, Analysis of a finite element upwind scheme, in: T.J.R. Hughes, ed., Finite Element Methods for Convection Dominated Flows, AMD Vol. 34 (ASME, New York, 1979). [7] T.J.R. Hughes and A.N. Brooks, A multi-dimensional upwind scheme with no crosswind diffusion, in: T.J.R. Hughes, ed., Finite Element Methods for Convection Dominated Flows, AMD Vol. 34 (ASME, New York, 1979). [8] T.J.R. Hughes and A.N. Brooks, A theoretical framework for Petrov-Galerkin methods with discontinuous weighting functions: application to the streamline-upwind procedure, in: R.H. Gallagher, ed., Finite Elements in Fluids, Vol. 4 (Wiley, New York, 1982).
T. Belytschko et al., Generalized Galerkin method for convection-diffusion problems
43
191 T.J.R. Hughes, T.E. Tezduyar and A.N. Brooks, A Petrov-Galerkin finite element formulation for systems of conservation laws with special reference to the compressible Euler equations, Proc. IMA Conf. on Numerical Methods in Fluid Dynamics, University of Reading, England (Academic Press, New York, 1982). [lo] J.W. Barrett and K.W. Morton, Optimal finite element solutions to diffusion-convection problems in one dimension, Internat. J. Numer. Meths. Engrg. 15 (1980) 1457-1474. [ll] J.W. Barrett and K.W. Morton, Optimal Petrov-Galerkin methods through approximate symmetrization, IMA J. Numer. Anal. 1 (1981) 439-468. [12] K.W. Morton, Generalized Galerkin methods for steady and unsteady problems, Proc. IMA Conf. on Numerical Methods in Fluid Dynamics, University of Reading, England 1982 (Academic Press, New York, 1982). [13] J. Donea, A Taylor-Galerkin method for convective transport problems, Internat. J. Numer. Meths. Engrg. 20 (1984) 101-119. [14] J.C. Heinrich and E. Envia, Finite element techniques in transport phenomena, Preprint. [IS] B.P. Leonard, A survey of finite differences of opinion on numerical muddling of the incomprehensible Dominated defective confusion equation, in: T.J.R. Hughes, ed., Finite Element Methods for Convection Flow, AMD Vol. 34 (ASME, New York, 1979). [16] O.C. Zienkiewicz, The Finite Element Method (McGraw-Hill, London, 3rd ed., 1977).