Stabilized finite element methods for steady advection—diffusion with production

Computer methods in applied mechanics and engineering ELSEVIER

Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191

Stabilized finite element methods for steady advection-diffusion with production Isaac Harari* Department of Solid Mechanics Materials and Structures, TeI-Aviv University, Ramat Aviv 69978, Israel

Thomas J.R. Hughes Division of Applied Mechanics, Durand Building, Stanford University, Stanford, CA 94305-4040, USA

Received 13 October 1993

Abstract

Finite element methods for solving equations of steady advection-diffusion with production are constructed, employing a linear source to model production. The Galerkin method requires relatively fine meshes to retain an acceptable degree of accuracy for wide ranges of values of the physical coefficients. Numerous criteria for improved accuracy via Galerkin/least-squares are proposed and examined in the entire parameter space. Of these, several provide the lowest error in nodal amplification factors, each in a certain range of ratios of physical coefficients. Guidelines for required mesh resolutions are presented for Galerkin and Galerkin/least-squares, showing that in large parts of the parameter space substantial savings may be obtained by employing Galerkin/least-squares. Galerkin/least-squares/gradient least-squares, a new variation of the Galerkin method, is designed to provide exact nodal amplification factors, offering accurate solutions at any resolution. Numerical tests validate these conclusions.

1. Introduction

A considerable a m o u n t of effort is being invested in developing cost-effective c o m p u t a t i o n a l m e t h o d o l o g i e s for problems of fluid mechanics. Numerical tools that can deal with realistic p r o b l e m s are emerging, as c o m p u t a t i o n a l techniques b e c o m e more sophisticated and operate on faster and m o r e powerful platforms. Experience in developing numerical m e t h o d s in a variety of fields indicates that the successful application of these m e t h o d s to large-scale problems is often p r e c e d e d by a d e e p analytical u n d e r s t a n d i n g of their p e r f o r m a n c e in simplified settings. This is instrumental in designing such m e t h o d s and validating their p e r f o r m a n c e on configurations of practical interest. Performing such an analysis, taking care not to preclude multi-dimensional generalizations, and designing finite e l e m e n t m e t h o d s for a d v e c t i o n - d i f f u s i o n with production is the purpose of this work. F u r t h e r m o r e , in this study we are initiating the d e v e l o p m e n t of finite element m e t h o d s for general s e c o n d - o r d e r elliptic partial differential equations. T h e extension of this investigation to absorbing media, thereby completing the t r e a t m e n t of the general case, will be presented in later work. H a v i n g reached a high degree of mathematical and algorithmic sophistication, finite element m e t h o d s , which are based on variational formulations, have b e c o m e the numerical technique of choice for n u m e r o u s classes of b o u n d a r y value problems. These m e t h o d s are also emerging as strong challengers to e n t r e n c h e d traditional approximate solution m e t h o d s in m a n y other fields, including a large class of problems in fluid dynamics (see, e.g. [1, 2] and references therein). Finite elements have a

* Corresponding author. 0045-7825/94/$07.00 ~ 1994 Elsevier Science B.V. All rights reserved SSDI 0045-7825(93~E0172-5

166

1. Harari, T.J.R. Hughes / Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191

rich mathematical background which can aid in method design and be used to prove convergence of numerical solutions to the exact solution with mesh refinement. However, accurately modeling problems by Galerkin finite elements may become prohibitively expensive in certain ranges of their parameters, and we explore in depth here methods of achieving this goal in the presence of production by less costly means. Modeling production as a linear source, the problem of steady advection-diffusion with production is characterized by two parameters that depend on its physical coefficients. As either of these parameters vanishes the governing equation tends to a limiting case, the advection-diffusion equation in one instance, and the Helmholtz and singular diffusion equations in the other. When both parameters vanish, the equation reduces to the Poisson equation. In this case finite elements based on the Galerkin method possess several advantageous attributes. Superconvergence is exhibited in one-dimensional models, e.g. linear elements yield solutions that nodally interpolate exact results, leading to minimal errors (measured in the H ~ semi-norm, related to error in flux), even on coarse meshes. Together with stability of the method, associated with symmetric, positive-definite coefficient matrices of the ensuing linear algebra problem, this indicates good performance on general, multi-dimensional configurations, which is guaranteed by proven optimal-rate convergence. We would like to endow finite element methods for advection-diffusion with a linear source with as many of these beneficial properties of the Galerkin method for the Poisson equation as possible. Advective-diffusive equations were examined as linear, steady-state models for the Navier-Stokes equations in the development of finite elements for problems of fluid dynamics. The first-order derivatives in this equation engender skew-symmetric coefficient matrices with potentially destabilizing effects. Numerical pathologies immediately became apparent (in a manner analogous to solutions obtained by finite difference methods that employ central difference schemes). Mathematical analysis led to the introduction of additional quantities into the formulation in order to demonstrate convergence. This enabled the design of stable finite element methods retaining the higher-order accuracy of the Galerkin method, at relatively low added computational cost. Such ideas were developed by Hughes and Brooks [3] and originally referred to as 'streamline-upwind/Petrov-Galerkin' (SUPG). In a more recent development, the concept of 'Galerkin/least-squares' (GLS) arose as a generalization of these ideas. This methodology is obtained by appending terms in least-squares form to the standard Galerkin formulation. The added terms contain residuals of the Euler-Lagrange equations of the boundary-value problem usually evaluated over element interiors, thereby preserving the consistency inherent in the Galerkin method (an important ingredient in obtaining improved convergence rates with higher-order interpolation) as well as respecting regularity requirements on the functions employed. The notion of Galerkin/least-squares crystallized in [4] and in the application of ideas of this sort to abstract mixed problems by Franca and Hughes (see [5] and references therein) and their colleagues, and was the key to recognizing these techniques as part of a general framework. In fluid mechanics, Galerkin/least-squares methods are identical to SUPG for hyperbolic cases, but the analysis in the presence of diffusion is simpler [4]. Galerkin/least-squares finite element methods have been successfully employed in a wide variety of applications (for an overview see, e.g. [6]). In the other limiting case, the governing equations contain an undifferentiated term in addition to second derivatives. The sign of this term in singular diffusion problems is such that ellipticity of the Laplace operator is preserved. However, solutions to these problems may contain sharp boundary layers, as in heat conduction with temperature-dependent strong sources due to chemical reactions, diffusion problems in semiconductors and elastic materials on elastic supports. Until recently, the performance of finite element methods for these problems was not analyzed in depth. Close examination of singular diffusion problems indicates that the standard Galerkin method suffers from shortcomings for these problems as well, potentially losing the capability to characterize derivatives of the solution (loss of stability in the H ~ sense, while retaining L 2 stability). The 'Galerkin/gradient least-squares' (GGLS) method proposed by Franca and Dutra do Carmo [7], in which the least-squares terms contain residuals of the gradient of the governing differential equation, successfully remedies this difficulty. This method is also intended for modeling complex boundary layer phenomena, such as arise in the analysis of thin structures. The Helmholtz equation is similar to singular diffusion equations but

1. Harari, T.J.R. Hughes / Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191

167

the sign of the undifferentiated term may also be reversed, thereby describing time-harmonic, propagating phenomena. In this case there is potential loss of ellipticity of the linear operator. Galerkin/least-squares technology has been exploited to relax wave-resolution requirements in these problems [8]. Results obtained by the Galerkin/least-squares and Galerkin/gradient least-squares methods are very similar, in some cases identical, for this class of problems [9]. In studying the advection-diffusion equation with a linear source ,,~e restrict ourselves to scalar, single field equations, but previous work on the limiting cases points to ways in Which the treatment may be be generalized. A simple change of variables shows that the constant-coefficient equation has exact solutions which are products of exponential functions (similar to solutions of advection-diffusion, but with halved rates) and solutions to a Helmholtz equation with a modified coefficient of the undifferentiated term. Solutions to the Helmholtz equation describe propagating or evanescent waves, depending on sign of the undifferentiated term. In summary, two types of physical behavior may thus be observed, either exponential variation at a composite rate or propagation that is modulated exponentially. Therefore, solutions may vary widely in nature depending on values of the coefficients, and possess very rapidly varying values in domain interiors as well as close to boundaries. In Section 2, a boundary-value problem for advection-diffusion with a source is presenied in its strong and weak forms. Finite element approximations, both the classical Galerkin method and several least-squares variants are presented in Section 3. A new combination of GLS with GGLS, 'Galerkin/ least-squares/gradient least-squares,' is proposed here. Experience from the analysis of Galerkin finite elements and design of Galerkin/least-squares and Galerkin/gradient least-squares methods for the above limiting cases indicates that it is beneficial to consider a model problem, typically in one dimension, with a known exact solution. Such a problem is proposed in Section 4. An analysis of numerical solutions obtained by employing linear finite elements in Section 5 indicates that relatively fine meshes may be required by the Galerkin method to retain an acceptable degree of accuracy. Numerous criteria for improved accuracy by employing Galerkin/least-squares are proposed and examined in the entire parameter space. Of these, several are found to provide the lowest error in nodal amplification factors, each in a certain range of ratios of physical coefficients. Guidelines for required mesh resolutions are presented for Galerkin and Galerkin/least-squares, showing that in large parts of the parameter space, substantial savings may be obtained by employing Galerkin/least-squares. The new technique, GLSGLS, is designed so that nodal values in the one-dimensional model are exact, leading to highly accurate computed solutions. Numerical tests in Section 6 bear out these observations. Closing remarks are presented in Section 7. Galerkin/least-squares methods other than the ones that offer the lowest error in amplification factors are presented for reference in Appendix A.

2. Boundary-value problem The domain J2 C I~d is an open set, where d = 2 or 3 is the number of space dimensions, with piecewise smooth boundary F. The outward unit vector normal to F is denoted by n = (n t, n 2. . . . . na). We assume that F admits the partition F =Fg LI Fh, where Fg t3 Fh = 0. We wish to study the following linear steady-state problem of advection-diffusion with a linear source as a model for transport equations of turbulent quantities such as kinetic energy or velocity scale (see, e.g. [10] and references therein): Find 4, : S')---~l~, such that -~b

=f

in O ,

(1)

4' = g

on Fg,

(2)

,¢v6 .n = h

o n r;,,

(3)

~ b := V" (~rV~b) - V " (u~b) + kZq~.

(4)

where

V. is the divergence operator and V is the gradient operator; K is the diffusivity tensor, assumed

168

I. Harari, T.J.R. Hughes / Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191

symmetric and positive definite, u is the given flow velocity and k 2 indicates the strength of the linear source (which can vary with the turbulence model); and f : O---~R, g : Fg---~R and h : Fh----~~ are the prescribed data.

R E M A R K 1. In an isotropic medium, the diffusivity is K = K1 where K is a positive scalar. Furthermore, If the medium is homogeneous, i.e. K is constant then v-(,¢v4,) --

K,a~,

(5)

where A is the Laplace operator.

R E M A R K 2. If the velocity vector is solenoidal, i.e.V, u = O, then v.(u4,) = u.v6.

(6)

R E M A R K 3. The coefficient of the undifferentiated term in (1) with an opposite sign is referred to as the absorption coefficient. In wave phenomena, this is related to evanescence. R E M A R K 4. The advection-diffusion equation and the Helmholtz equation may be considered as two limiting cases of Eq. (1), with k = 0 in the former and u = 0 in the latter. Both these cases have been treated in the finite element literature, see, e.g. [3, 8]. REMARK variables

5. In the isotropic, homogeneous case (k actually need not be constant) the change of U'X

4, = qJ exp(--~--K )

(7)

transforms Eq. (1) to a standard Helmholtz equation with a modified coefficient of the undifferentiated term [ l l , p. 231]. Solutions are therefore products of the exponential functions in (7) (similar to solutions of advection-diffusion, but with halved rates) and solutions to a Helmholtz equation. Solutions to the Helmholtz equation describe propagating or evanescent waves, depending on sign of the undifferentiated term. Two types of physical behavior may thus be observed, either exponential variation at a composite rate or propagation that is modulated exponentially. Therefore, solutions may vary widely in nature depending on values of the coefficients, and possess very rapidly varying values in domain interiors as well as close to boundaries. The Helmholtz equation is well understood since it has been treated in depth in the mathematical physics literature, e.g. [11-13]. Existing finite element methods for this equation such as [7-9] may thus be exploited to solve the boundary-value problem for ¢, from which th is readily obtained. Alternatively, the change of variables may be employed as a tool in the analysis and design of finite element methods that are applied directly to (1). Let L2(g2 ) denote the space of square-integrable functions on O, and H~(/2) denotes the space of functions in L2(O ) with generalized derivatives also in L2(O ). The variational form of the boundaryvalue problem is stated in terms of St, the set of trial solutions, and ~, the space of weighting functions. For the Neumann problem both spaces are equal to H~(O). Otherwise, Dirichlet boundary conditions on Fg and their homogeneous counterparts must be satisfied by functions in 90 and 7/', respectively. The objective of the variational problem is to find 4, E 90 such that Vw E

A ( w , ~) = L(w) ,

(8)

where

A ( w , 4)):= (Tw, K74~) - (Vw, urb ) - (w, k24)), L ( w ) : = (w, f ) + ( w , h ) r ,

(9) (10)

and (. , . ) is the L2(J'2) inner product. Subscripts on inner products denote domains of integration other

1. Harari. T.J.R. Hughes I Cornput. Methods Appl. Mech. Engrg. 115 (199,1) 165-191

169

than O. The operator A(., .) is a bilinear form, but clearly neither symmetric nor an inner product, and L(.) is linear. Both the skew-symmetric second term and the negative third term in (9) have potentially destabilizing effects. These observations presage possible difficulties for numerical methods that are based directly on (8).

3. Discrete formulations

Consider a partition of the computational domain into finite elements such that ~ is the union of element interiors. Let behc be, 7 r h c ~ be sets of finite element functions consisting of continuous piecewise polynomials of degree 1. We consider several discrete formulations. 3.1. Galerkin This is the classical method: find ~bh E beh such that V w h E 7rh A(wh, cbh)=L(w"),

(11)

where h is the mesh parameter. The coefficient matrix is made up of a linear combination of standard finite element symmetric stiffness and mass matrices, and a skew-symmetric matrix. The Galerkin formulation is consistent in the sense that ~b, the exact solution of the boundary-value problem (1)-(3), satisfies (11). This property is an important ingredient in obtaining improved convergence rates with higher-order interpolation. The Galerkin method is an approximation of (8) and, as observed, stability may be degraded. This is demonstrated by considering A ( w h, w h) which yields in the isotropic case

m(w",

w") =

'-'Vwhll 2 - ( V w h,

uw

-Ilkwhll

2,

(12)

where I111 is the L2(O ) norm. The second term vanishes and does not contribute to stability for problems with homogeneous Dirichlet boundary conditions and solenoidal velocity fields (this follows from integration by parts). The stability of the Galerkin method for relatively small K or large k is thus quite problematic and we anticipate that it may not always lead to good numerical solutions. We therefore consider variants of the Galerkin method obtained by the addition of terms that contain least-squares forms of residuals. These modifications allow extra freedom in the design of the methods in order to improve their performance, while maintaining the consistency that is inherent in the Galerkin method. 3.2. Galerkin ~least-squares The variational form (8) is adjusted to include a residual of the governing differential equation (1) (which, in the continuous case, is identically zero). The resulting approximation, termed Galerkin/ least-squares (GLS), is AGLS(wh , ~bh) = L(wh) - (*~wh , ~'f)h,

(13)

where A GLs(Wh, ~bh) := A ( w h ' ~bh) + (.LPwh, r~Lpcbh)f~,

(14)

and ~- is a yet undetermined parameter (with dimensions of time) that enables the design of the method to satisfy certain criteria. 3.3. Galerkin / gradient least-squares Here the modification to the original form is in terms of the gradient of the differential equation ACGLs(W h, 6 h) = L ( w h) - (T&C'wh' rvT£)~,

(15)

I. Harari, T.J.R. Hughes / Cornput. Methods Appl. Mech. Engrg. 115 (1994) 165-191

170

where A GGt.s(W", ~bh): = A(w", ~") + (V~w", rvVL~'~bh)fi.

(16)

In this case the parameter r v is likely to take on different values from r in Eq. (13).

3.4. Galerkin ~least-squares/gradient least-squares Combining the above two ideas yields A GLSOLs(W",q~") = L(w") - (L~'w", rf)~ - (~7~wh, rvVf) ~ ,

(17)

where A OLSGLs(W",~b"):= A(w h, dp") + (~w", r~4)")h + (V~w h, rrV~qSh)h •

(18)

This concept, not considered previously elsewhere, gives more freedom in designing the numerical method, while retaining the consistency of all the above methods.

R E M A R K 1. The subscript ,O indicates that the least-squares integrals are evaluated over element interiors, and therefore lack of global smoothness of the finite element functions creates no problems. Furthermore, r and r v should be thought of as local quantities, to be defined consistently on the element level in terms of the element size, highlighting the applicability of these methods to general unstructured meshes.

R E M A R K 2. The standard Galerkin method can be obtained as a special case of any of the least-squares variants by setting ~-= 0 and r v = 0.

R E M A R K 3. Consider the implementation of linear finite elements. All derivatives of order higher than one in the least-squares terms vanish. It was shown in [9] for the Heimholtz and singular diffusion equations that the Galerkin/gradient least-squares method is quite similar to Galerkin/least-squares. Indeed both methods employing a uniform mesh of linear elements may produce identical solutions for any d-dimensional Dirichlet problem with a uniform source distribution. However, for the advectiondiffusion problem, the Galerkin/gradient least-squares method does not modify the Galerkin equation, and thus cannot remedy its shortcomings.

R E M A R K 4. Employing linear elements, the coefficient matrix of each of the least-squares formulations is still made up of a linear combination of standard finite element stiffness and mass matrices, and a skew-symmetric matrix, adding virtually no complexity beyond the Galerkin method. The Galerkin/ least-squares variation to Galerkin is achieved by simply modifying the scalar coefficient multiplying the stiffness contribution by a factor of 1 + zlul'-/K (in the isotropic case), where lul-" :-- u - u, and the scalar coefficient multiplying the mass contribution by a factor of 1 - rk 2. The load terms are adjusted accordingly to maintain consistency. These operations are performed on the element level. The matrix equations of the other least-squares formulations may be obtained in similar fashion.

4. Model problem Experience from previous analysis of Galerkin finite elements and design of Galerkin/least-squares and Galerkin/gradient least-squares methods for a variety of applications indicates that it is beneficial to consider a model problem, typical in one dimension, with a known exact solution. In the homogeneous case the governing differential equation on the open unit interval is K6xx - u~b.,. + k'-~b = O,

(19)

where the physical coefficients are assumed constant and an inferior comma denotes differentiation.

1. Harari. T.J.R. Hughes / Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191

171

This is our model problem. The exact solution may be decomposed by the aforementioned change of variables (7) 4~ = ~0exp

,

(20)

where the solution to the resulting Helmholtz equation is a linear combination of trigonometric or exponential functions {

,flu 2

k 2 "~

g,=expk,-+V-4-7-: Z x,)

(21)

describing propagating or evanescent waves, respectively (depending on the sign of u 2 - 4Kk2), such that boundary conditions are satisfied. Two types of physical behavior may thus be observed, either exponential variation at a composite rate or propagation that is modulated exponentially, i.e. ~b = exp(/3x),

(22)

where 1

/3 = ~-K (u + Vu 2 - 4Kk2).

(23)

Solutions may vary widely in nature depending on values of the coefficients, possessing rapidly varying values in domain interiors as well as near boundaries. For example, solutions to the Dirichlet problem on the unit interval with th(0)= go and ~b(1)= gt are

,/,=

u X/ sin g0 exp \2-K-K

+g,

(

)

U

~

-- 2 (a-x) 4K ( V kz

;

UZ ) ) / s i n V k2 4,.,x ;

uz

(24)

Figures 1 and 2 show some of the variety of possible solutions to this problem. For k = 0 , the weU-known outflow layer of width O(K/]U]) is evident. A comparison of the two figures when k z > 0 (note the difference in scales) indicates that a small modification in boundary values can result in drastic changes in the solution. Values of the exact solution at uniformly-spaced points xa = A h may be parameterized by two non-dimensional parameters, that contain the physical coefficients and the mesh size uh % : - 2K '

a~:-

(25)

(kh) 2 12K

(26)

The first of these is called the element Peclet number. These parameters represent the relative weights of the first-order and undifferentiated terms. When either approaches zero the problem tends to one of the limiting cases. Solutions at the discrete points are thus linear combinations of (exp/3h) a where /3h = a,, +- V a ~ - 12ak,

(27)

emphasizing the two parts of the solution that emanate from the change of variables (20), one part that corresponds to advection-diffusion and the other to Helmholtz. The two types of physical behavior are readily expressed in terms of the numerical parameters exponential variation at a composite rate,

a k / a ] ~< 1/12,

exponentially modulated propagation,

% / a ~ > 1 / 12,

since a k / a ] = K k 2 / ( 3 u 2) is independent of the mesh.

I. Harari, T.J.R. Hughes / Cornput. Methods Appl. Mech. Engrg. 115 (1994) 165-191

172

~2/~ = 0

1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.I 0

'

I

'

I

'

I

'

I

~ /

ul~ = o

2---

,--.

.

.

/

/

" .

.

0.2

0.4

;.

..'./I.

~

. / _/". . /'/5"! /.'/

0.6

1.0

0.8

k2/~ = 4

1.2

'

I

'

I

'

I

'

I

'

/.;~

0.8 1.0 I 0.6 0.4 0.2 0 0

0.2

0.4

0.6

0.8

1.0

k2/K = 16 1.5

I

I

I

I

~llIl

l, 0

0.5

~,,,,~,'~

... •.

-0.5 -1.O -1.5

0

0.2

0.4

0.6

0.8

1.0

Fig. I. Exact solutions to the model problem on a unit interval with &, = 0 and g, = 1 (bold lines denote propagating solutions).

REMARK.

Definitions of the non-dimensional p a r a m e t e r s m a y be generalized by k n o w n p r o c e d u r e s to multi-dimensional configurations [14], e.g. a,,:-

lulh

2K '

(28)

and to systems [15]. Solutions to the Dirichlet p r o b l e m at N + I points x A are 4~(xa) = (go e x p ( A a , ) sin((N - A ) ~ 1 2 a k - a ~ ) + g~ exp((A - N)a,,) s i n ( Z ~ 1 2 a k - a~))/sin(N~/12ak -- a ~ ) .

(29)

REMARK.

At the transition b e t w e e n the two types of physical behavior, 12a k = a] solutions are a linear c o m b i n a t i o n of e x p ( u / 2 K ) and x e x p ( u / 2 K ) . In this case solutions to the Dirichlet p r o b l e m at the discrete points are

6(XA) = (g0(N - A) e x p ( A a , ) + gtA exp((A - N ) a , , ) ) / N .

(30)

1. Harari, T.J.R. Hughes / Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191

173

k2/~ = o 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.I 0

. . . . .

0

'

0.2

0.4

0.6

0.8

I

I

1.0

k-~/,c = 4 12

I

I

1°I

.I"

8

\

/ o/

6

/ w- °

/"

4

~,*

\

, f"

'\ .,'"

,

,,

" .....

'~

"', ",.

,,"

~-

2 ...L_ 0.2

0 0

0.4

0.6

0.8

1.0

k2/~ = 16 15

'

I

'

I

'

I

;

I

'

10 5 "~

0 -5 -10 -15 0.2

0.4

0.6

0.8

1.0

x

Fig. 2. Exact solutions to the model problem on a unit interval with go = 1 and g~ = 0 (bold lines denote propagating solutions).

T h e t w o e x a c t n o d a l a m p l i f i c a t i o n f a c t o r s , e x p flh, c h a r a c t e r i z e d i s c r e t e p o i n t s o n v a l u e s at n e i g h b o r i n g p o i n t s

~b(xa ) -

the dependence

of solutions

at

exp( 2a,,)gO(Xm_ l ) + ~b(xa+,) 2 exp(a,,)cs (~b(xA - t ) + &(XA +1 )) c o s h a,, + (~b(xA _, ) - 4'(XA +1 )) s i n h a,,

2cs

,

(32)

where

cs:=cos

V1 2 a ~ - a ] = c o s h

Va , 2, - 1 2 a , .

(33)

5. Linear finite elements W e c o n s i d e r t h e d o m a i n to b e d i s c r e t i z e d by a u n i f o r m m e s h o f l i n e a r e l e m e n t s o f l e n g t h h. V a l u e s o f t h e finite e l e m e n t s o l u t i o n at n o d e A a r e d e n o t e d ~ba : = C~h(XA). N o d a l a m p l i f i c a t i o n f a c t o r s o f t h e n u m e r i c a l s o l u t i o n a r e e x a m i n e d in o r d e r to c h a r a c t e r i z e its p e r f o r m a n c e .

174 5.1.

1. Harari. T.J.R. Hughes / Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191 Galerkin

The Ath equation of the tridiagonal system is ( - ~ + 1 + 2~bA - 4'.~- 1) + °~,,(~bA., -- ¢A-1) -- 2 % ( ~ + 1 + 4 ~ + qbA_1) = 0.

(34)

The terms multiplied by % are 'mass-type' terms and emanate from the L_, part of the operator A ( . , .) in Eq. (9). The terms multiplied by % emanate from the skew-symmetric part of the operator. The remaining terms are 'stiffness-type' and emanate from the H 1 portion of the operator. One may thus anticipate degradation of numerical solutions with increasing values of % and %. The non-dimensional parameters % and % are element-based quantities. The results presented here are thus applicable to general unstructured meshes, even though the analysis itself is performed on a uniform mesh. The resulting stencil is REMARK.

(1 - % + 2 % ) ~ + 1 - 2(1 - 4 % ) ~ + (1 + % + 2 a k ) ~ _ , = 0.

(35)

Seeking solutions in which ~ , ~ - p a we obtain a quadratic polynomial in p. The two solutions to this equation are the Galerkin nodal amplification factors 1 - 4c~k + ~c~] - 12ak(1 -- ~k) P =

(36)

1 + 2 % -- %

We refer to the larger in absolute value of the two solutions as PlIn the two limiting cases, the amplification factors have special properties simplifying the analysis. In advection-diffusion % = 0, there is an error only in Pt since P2 = 1 which is exact. For the Helmholtz equation % = 0 , the product of the amplification factors is exact PIP2 = 1 and so one need only investigate errors in their ratio which is related to the phase of propagation. In the general case examined here, there are errors in both factors (or in both the product and the ratio) considerably encumbering the analysis and subsequent method design. Two gross numerical pathologies are immediately apparent, one in each region of physical behavior. When exact solutions vary exponentially, %~
(37)

i,e. spurious oscillations occur, beginning at mesh sizes of 3u +- ~ / 9 u 2 - 2 4 K k 2

h =

2k 2

(38)

(In the advection-diffusion limit of k - + 0 , oscillations start at h = 2K/u.) Pl becomes negative at the smaller h and P2 at the larger one (at the other critical points the numerator and denominator of each p switch signs simultaneously). Note that (37) applies to problems in the range a k ~ < a ] / 8 , slightly overlapping the region in which physical solutions are propagating. The nodal amplification factors should be complex to represent propagation when ~k>C~]/12. However, the Galerkin amplification factors lose their imaginary component when c~2, _ 12C~k(1 _ O~k)> 0.

(39)

It is convenient to examine the behavior of numerical solutions in the V'-%kO~,,plane (see Fig. 3). For each given physical problem, h changes in this plane along a straight line emanating from the origin. Specifically, the abscissa and the ordinate correspond to the Helmholtz and advection-diffusion problems, respectively. Spurious oscillations in Galerkin solutions occur above the dashed curve % = 1 + 2%. The transition between propagating and evanescent exact solutions is marked in this plot by the bold straight line % = ~ (propagation occurs below this line). The numerical analog follows

1. Harari, T.J.R. Hughes / Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191 6

I

'

I

'

I

175

'

!

advection5 l- diffusion

| l g

physical

a

• spurious ,f" jr." 2 [- { oscillations f J . ' / [ .~'~'U"---""~ 1~ 0 I."

0

~

~ ,

numerical

' SUofl]~it~ntn~/ t ran:iti0n :,

,

0.5

I

1.0

~.~oltz

1.5

Fig. 3. Behavior of exact and Galerkin solutions in the parameter space. from (39), namely the solid curve a,, = ~/12ak(1 - ak) in the lower part of the plot. Clearly, Galerkin solutions properly represent physical propagation, even qualitatively, only in part of the parameter space. Degradation of computed solutions is to be expected when waves are insufficiently resolved by a mesh. We consider four elements per wave to be the minimum number sufficient to resolve propagating phenomena. This is shown by the dotted line a,, = ~/12otk - ,rr2/4. The portion of the plane below this curve represents meshes that cannot be expected to provide meaningful numerical solutions by any method. Beyond the two gross pathologies evident in the Galerkin representation, we would like to investigate its quantitative performance more closely. The nodal amplification factors are complex and thus somewhat cumbersome to analyze. In order to facilitate examining the accuracy of Galerkin solutions (of exponential rates in one case or phase and amplitude in the other), it is convenient to decompose the numerical amplification factors in a manner similar to the aforementioned change of variables of the continuous solution (20) which leads to the composite exponent (27),

~ +2ak+a,,1-4ot_k+--j_._a_~--l_2ffk(__l]-- ~k) P =

+ 2o~k

~/(1 + 2ak) 2 - a~

a,,

V~,p2~-expa,,

~-*'

(40)

= exp _+~/a~ - lZa~

The product of the amplification factors thus corresponds to the exponential solution related to advection-diffusion and their ratio, the propagating or evanescent wave emanating from the modified Helmholtz problem. We may thus employ two real quantities to characterize the quality of numerical solutions: the modulation error which is related to the product of amplification factors, and the phase error emanating from their ratio. (Note that in propagation, when the amplification factors should be complex, p, = P2, i.e. IPt/P2I = lexp = 1 and there is only an error in phase related to the ratio of amplification factors.) Figure 4 shows the relative error in the product of the Galerkin amplification factors. When the relative error is less than - 1 , below the bold line (i.e. PtP2 < 0 ) , there are oscillations in the Galerkin solution. For a~ = 0, we have the well-known behavior of solutions to the advection-diffusion equation, with the onset of oscillations at a,, = 1. There is severe degradation of the Galerkin representation of the modulation well before oscillations occur. Up to a k = a ~/ 12 this figure represents the error in one of the two exponential functions that compose the solution. Beyond this value, the curves show the error in exponential modulation of the propagating wave. The onset of spurious oscillations occurs at

~/a,~ 12akl

176

1. Harari, T.J.R. Hughes 1.0

Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191

/

''/)'

t 0.5

I

0

/

/ /

/

/ /

i

./

1/12--i

i~y.~-;=___

~3

"~

+'

X\ -0.5

1/6-- '

\,

xx "'.

1 .....

~,

\ N

\,

-1.0

-1.5 0.5

1.0

1.5

2.0

2.5

3,0

ol u

Fig. 4. Modulation e r r o r - - r e l a t i v e error in the product of Galerkin amplification factors. Curves below the bold line c o r r e s p o n d to oscillating solutions.

higher values of a,, as the ratio a k / a ~ (which is independent of the mesh) increases. Beyond a k = a ~ / 8 , no oscillations occur. Furthermore, the accuracy of the Galerkin representation of the modulation improves. There is no modulation error at any a k in the Helmholtz limit a,, = O. The error in the ratio of the Galerkin amplification factors is examined separately for evanescent and propagating solutions. The former is shown is Fig. 5, which is similar to Fig. 4 limited to the range 0 <~ a k / a ~ ~< 1/12. Increasing a,, always leads to spurious oscillations in this range, indicated by relative errors of less than - 1 (below the bold line). For advection-diffusion a k = 0, the error in the ratio is identical to that of the product since P2 = 1. As the ratio a k / a ~ increases, there is less improvement in the accuracy here than there was in the modulation error. Furthermore, at the transition ct~, = 12a k, there is significant degradation in accuracy. This observation presages difficulties for the numerical methods in the vicinity of the physical transition, apparently due to the change in the nature of solutions. For propagating solutions, the amplification factors are complex conjugates so the magnitude of their

1.0

'

'

'1

11

I

~,/'

I

ok/o~ = 0 - 1/24---1/13----1/12---

o.5 / / / . !

0

-0.5

-1.0

-1.5

L

I 0.5

,

I 1.0

,IX

I, 1.5

I/, 2.0

i 2.5

, 3.0

Fig. 5. Dispension e r r o r - - r e l a t i v e error in the ratio of Galerkin amplification factors for exponentially varying solutions. C u r v e s below the bold line c o r r e s p o n d to oscillating solutions.

I. Harari, T.J.R. Hughes / Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191

1.6 ~~

, a~/a~

1.5

=

0 -

177

,

/\

1.4 ......

f"

1.3

~~

1.2 1.1 1.0

L.-L_. . . . . . . .

0.9 0.8

0

0.2

0.4

0.6

1" 0.8

1.0

4-aZ

Fig. 6. Phase error in the Galerkin solution for propagating solutions. (Bold lines denote physical and numerical transition to propagation on the left and right, respectively; solid curve indicates limit of sufficient resolution.) ratio equals unity which is the exact value. However, there is a phase error in the Galerkin representation (see Fig. 6). Curves of phase error are restricted by the bold lines denoting physical transition to propagation on the left and the numerical transition on the right. In the Helmholtz limit % = 0, there is the known phase lag in the Galerkin solution [9]. However, as the transition is approached with increasing values of a 2u / a k, the phase error becomes a lead and, in general, the error becomes quite severe. The solid curve indicates the limit of sufficient resolution. G o o d numerical performance is not expected to the right of this line, but improved accuracy is desirable to its left. All in all, Figs. 4 - 6 give an indication of the ranges of parameters at which the Galerkin method performs well, and what improvements might be called for. Nodal values of the Galerkin solution for the Dirichlet problem are A

N-A

C~A = (gO(PiP2) (Pt

N-A

A N N ) + g t ( P ,A -p2))l(p, -p,_)

,

(41)

cf. (29) for the exact counterpart. The dependence of Galerkin solutions on neighboring nodal values, cf. (32) for the exact case P t P 2 ~ A - t + ~A+I

¢~A =

/91 +

P2

((])A- 1 q- ~bA*I)(1 + 2%) + (~bA_, -- &A+,)% 2(1 - 4 % )

(42) (43)

where

P~+P2-

Pl + P2 Pl P2

2(1 -4%)

(44)

1+2%-%'

2( 1 - 4ork ) 1 + 2a k + % "

(45)

These two quantities, p, + P2 and (p~ + R2)/(PlP2), are real even when p is complex and they may also be employed to examine the performance of numerical methods since the nodal amplification factors are uniquely determined by these terms, namely

178

1. Harari. T.J.R. Hughes I Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191

1(

P =2" 5.2.

PI +p2---

V

,

(Or + & ) - - 4 ( p

'

+

].

(46)

Galerkin ~least-squares

The Galerkin/least-squares algebraic equations for a uniform mesh of linear elements may be obtained from the Galerkin equations by replacing % and % by % "= 1 + 4 a ] r K / h

(47)

2 '

1 - 12akrK/h 2

8~ := % 1

+

~

_ ,

4a?,rK/h-

'

(48)

respectively. The preceding analysis may then be repeated except that the numerical solution now depends on 8,, and 8k, which can be determined to satisfy certain design criteria by the selection of r in order to improve the performance of numerical methods. The Galerkin method is clearly a special case of GLS obtained by setting r = 0. The performance of finite element methods is generally assessed on the basis of the value of the error measured in an integral norm. In many previous applications, Galerkin/least-squares methods employing linear elements were designed to provide exact nodal amplification factors (e.g. [3, 9]) which is a property of the Galerkin method for the Laplace operator in one dimension. This criterion assures a high degree of accuracy and, in particular, nodal exactness is attained for Dirichlet problems with homogeneous equations. In one dimension, linear finite element functions that nodally interpolate the exact solution minimize the gradient of the error, typically related to error in flux, in the L 2 norm. Nodal exactness is thus a bona fide measure of the performance of linear finite elements. For the advection-diffusion equation, one of the amplification factors (p,) is a (1, 1) Pad6 approximation of the exact exponential, and the other factor equals the exact constant. The G a l e r k i n / least-squares method is easily designed to correct the first factor, thereby yielding exact amplification factors. For the Helmholtz equation, the amplification factors have the property of being complex conjugates of unit magnitude, just like the exact factors. In this case, there is only a phase error, which Galerkin/least-squares is designed to correct in order to provide exact amplification factors. For the more general case that is being examined here, this is not possible. There is an error in both exponential rates when solutions vary exponentially, and an error in magnitude and phase for propagating solutions. Both these errors cannot be corrected simultaneously by the single design parameter r in GLS. This is easy to demonstrate by examining the magnitude and phase of the Galerkin/least-squares solution. For an exact magnitude, we require Pl P2 = exp 2% ,

(49)

or, in terms of the non-dimensional parameters, 0/u

1 + 28 k - tanh % .

(50)

Clearly, this requirement is identically satisfied in the Helmholtz limit. The phase is related to the ratio of amplification factors. An exact phase requires P , / P 2 = exp 2~/a~ - 12a k

(51)

and in terms of the Galerkin/least-squares parameters, 1 - 48 k ~ ( 1 + 2ak)-

-2

= ¢s.

- - O/it

Solving Eqs. (50) and (52) for the non-dimensional parameters yields

(52)

I. Harari, T.J.R. Hughes / Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191 1

cosh

179

a,, - cs

c~k - 2 2 cosh a,, + cs '

(53)

3 sinh a,, a,, - 2 cosh a,, + cs "

(54)

In both the Helmhoitz and advection-diffusion limits, one of these equations is trivially satisfied in each case (both sides vanish) and the other determines the known design of ~-. For example, in advectiondiffusion a k = 0 which leads to ~k = 0, satisfied by definition, and c~,,= tanh a,,, the nodally exact solution. But in the general case both equations need to be satisfied independently, which is not possible with a single parameter. Therefore, other design criteria are needed, each leading to a Galerkin/least-squares method with certain properties. Evaluation of the numerical error is based on the amplification factors, which are related directly to the nodal values. However, in designing the Galerkin/least-squares methods, we also consider the magnitude and phase, since these are always real. As in the analysis of the Galerkin method, these quantities are related to the product and ratio of the amplification factors. Another pair of real quantities that is employed in the method design is composed of the terms pt + P2 and ( P t + P z ) / ( P t P 2 ) which determine the dependence of nodal values of a computed solution on neighboring values. 5.2.1.

D e s i g n criteria

Numerous criteria for improved accuracy via Galerkin/least-squares were examined in the entire parameter space. Of these, several are recommended for use in practice as each was found to provide the lowest error in nodal amplification factors in a certain range of physical coefficients. The recommended guidelines are listed in the following. Additional criteria are listed in Appendix A for reference. These may prove useful in special applications. 1. C r i t i c a l : The critical method guarantees no spurious oscillations in the region of exponential variation a~ > 12a k by designing z so that 1 + 2c~k - ~,,/> 0.

(55)

This is obtained by a Galerkin/least-squares method with '0, 1 + 2a k

rK h "- -

1 + 2a k - a , , / > 0 , (56) - a,,

4(6a~ -a~,) '

1 + 2a k

-a,,

< 0

(see Fig. 7)). We note that in the second case, the denominator vanishes for certain combinations of the non-dimensional parameters. This need not be a concern since the accuracy of the amplification factors deteriorates to an extent that the solution is rendered meaningless long before this region of singularity is approached. 2. E x a c t p : The Galerkin/least-squares method can be designed so that o n e o r t h e o t h e r of the amplification factors is exact p = exp/3h ,

(57)

where we select Pt for a,, > 0 and P2 otherwise. Since p satisfies the quadratic equation that follows from (35), in terms of the Galerkin/least-squares parameters, we require (1 + 2c~k) cosh(/3h) - ~,, sinh(/3h) - (1 - acid) = 0 .

(58)

This is obtained by zK h2

-

(1 + 2ak) cosh(/3h) - a,, sinh(flh) - (1 - 4ak)

(59)

4((6a~ - a~,) cosh(flh) + (12a~ + a~))

Although amplification factors should be complex when a~. > a~,/12, solutions to the boundary-value problem ( 1 ) - ( 3 ) as it is posed are real. However, the Galerkin/least-squares parameter in (59) is

180

I. Harari, T.J.R. Hughes / Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191 0.10

I

I

I

I

I

ok/o~ = o - -

0.09

1/48---1/18 . . . . 1/12------

0.08 0.07 0.06

% "--

0.05 0.04 /

0.03

./"

0.02 O.O1

I

0

I

I

0.5

1.0

1.5

,

2.0

I 2.5

, 3.0

~u

Fig. 7. Critical Galerkin/least-squares method.

required to take on complex values in propagation (when/3h becomes complex) due to the presence of the hyperbolic sine (except at the Helmholtz limit). In other words, employing this criterion in the propagation region leads to a method with the undesirable property of requiring a complex parameter although the problem has real exact solutions. We therefore restrict the applicability of this criterion to exponentially varying solutions ak~
I

I

I

I

0.09 0.08 0.07

~k/~2,

. " ~~ _ _ ~ --

~

~

"~

~

I

0 1/48-=

1/18

........

~

-

-

- -

-

-

-

1 / 1 2 - - -

0.06 0.05 •****-,

'****-.

0.04 0.03 0.02 0.01 0

I

0.5

~

I

1.0

,

1

1.5

L

I

2.0

~

I

2.5

Ot u

Fig. 8. Galerkin/least-squares method for exact p~.

J

3.0

1. Harari, T.J.R. Hughes / Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191 dp.4 = ( g o ( N -

181 (60)

A)p a + g,Apa-N)/N,

cf. (30). This is achieved with G L S by requiring that the phase (or the ratio of amplification factors) be exact, employing -cs=(1 + 2t~k)(a 2. - 6a~) + (1 - 4t~a)(a2, + 12a 2)

rK

h2

2 2 4(cs 2 (a?,9 - 6 a 2 ) 2 - (aT,9 + 12ak) ) •t /

2~

2,,

"~

2

'~

2

"~

csvotutcs ta'u - 6ot2) 2 - (ot2u + 12ak)" ) + 3 6 a k ( 3 a k + a~,) + 4(csE(a~, -- 6 a 2 ) 2 -- (a~ + 12a2) 2)

(61)

This equation holds for a k <~ [ ( X / 3 - 1)/6]a2, (which includes the transition). For a k > [ ( V 3 - 1)/6]a2, the sign on the second term is reversed. This m e t h o d reduces to the well-known G a l e r k i n / l e a s t - s q u a r e s m e t h o d that provides nodal exactness in the Helmholtz and advection-diffusion limits. For a d v e c t i o n diffusion with a source, the applicability of this m e t h o d is limited to the vicinity of transition (which is indeed problematic) and it is somewhat costly for computation. 4. Offset: In this case the relative errors in (Pl + Oz) and (t91 + P 2 ) / ( P l P 2 ) , the coefficients in Eq. (42), P,P2

2cs exp a,,

2cs -

-

Pl + P2 exp a,,

-

1 =

1

Pt + P2

,

(62)

are offset by employing zK ((l + 2 a k ) C O s h a , , - a , s i n h a , , ) c s - ( 1 - 4 a k) h" 4((6a~ - a ] ) c s cosh a,, + (12akz + a ] ) )

(63)

(see Fig. 9). Again, this m e t h o d is nodally exact in the Helmholtz and advection-diffusion limits. In general, it is similar to the Phase m e t h o d , but simpler and hence less costly. 5.2.2. R e c o m m e n d a t i o n and range o f applicability T h e numerical methods are r e c o m m e n d e d on the basis of errors in nodal amplification factors. T h e measure of the relative error is defined as

e"p := l P t / e x p fllh -

112 +

(64)

Ip2/exp/32h - 1[ 2 .

I

I

I

I

I

0,04 0.02 0 -0.02 -0.04 -0.06 -0.08 -0.I0

-_gf

10 B - -

50 - - I00 . . . . . .

v~ i

I

I

0.5

1.0

,

I

1.5

~

I

2.0

~

I

2.5

,

3.0

Otk

Fig. 9. Galerkin/least-squares method for offsetting relative errors of coefficients in (42).

182

1. Harari, T.J.R. Hughes / Cornput. Methods Appl. Mech. Engrg. 115 (1994) 165-191

Figure 10 shows several cases in the region of exponential variation. The Exact p~ GLS method (59) provides the best results close to the advection-diffusion limit, and the Phase method (61) and the Offset method (63) virtually coincide. (Recall that at the advection-diffusion limit, ep = 0 for all these Galerkin/least-squares methods.) Close to the transition the Phase method is the best while the other two GLS methods are very similar. The Phase method is more costly to employ, however, and thus the Exact p~ method is recommended in the entire region up to propagation. In the region of propagation itself the Exact p~ method was seen to call for a complex parameter so it is not recommended here. As a matter of fact in the intermediate region the improvement of any of the Galerkin/least-squares methods over Galerkin is not of great consequence. However, as the effects of propagation become stronger, the GLS methods may offer markedly better results (see Fig. 11). In these cases, results for the Phase method are not plotted since its performance is quite similar to that of the Offset method, yet more costly to compute. At the Helmholtz limit, the Offset method leads to e. = 0 . Of the methods proposed, the one that provides the lowest error in nodal amplification factors in each range of ratios of the physical coefficients is thus recommended, namely

ok = ~/48

'~¢"

u u Galerkin - -

20 18 16 14 12 10

I

G LS/Offset -- -- --

0

1

,

_ ~

.I

/

/ /

7/ ,

/

/

/

4 0 -/

/

/

6

E-5

' /I

/

G L S / E x a c t Pl -- -GLS/Ph~e -----

8

/l

,

,

2

,

3

Ou

ok = og,/18 20 18 16 14 12 "L~"

I

/,

10

/

/i

E-5

I

0.2

--

0.4

20 18 16 14 12 10

/ /

//

/// I

i

0.6

0.8

I

]

[

i

1.0

/I

1,2

/'# //

//

8

E-5

/

/

8 6 4

~'~"

,

I

'

6 4 2 0

/. I

0

0.I

0.9

0.3

0.4

0.5

0.6

~u

Fig. 10. Relative error in nodal amplification factors in the range of exponential variation.

1. Harari, T.J.R. Hughes

I

Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191 Ot k ~

%~"

E-5

I

2°I 18 16 14 12 10 8 6 4 2 0

I

183

O t u2

'

I

'

I

Galerkin ~

' /I

I

/i II

GLS/Offset - -- --

, 0

I 0.03

0.06

0.09

0.12

0.15

0.18

Ot u

a~ = 10a~

%~"

20 18 10 14 12 10 8 6 4

I

'

I

'

Ii

I

..,

l I I I / f

i

I

0.02

0.04

I

I

~

l

0.06

0.08

0.10

~t u

ak : 100 ~

%~"

E-5

20 18 16 14 12 10 8 6 4 2 0

I~ 1

~

/

I

I

/ / / /

0

F_ 0.01

--I-- ~ ~ ' ~ 0.02

,/

/

/

/

/

jJ I

0.03

i

0.04

Otu Fig. 11. R e l a t i v e e r r o r m n o d a l amplification factors in the r a n g e o f p r o p a g a t i o n .

(59),

]rK/h2=O, L(63),

ak/a~,
Galerkin/least-squares with (63) may also be used for increased accuracy near the transition in the range 1/12<-Ctk/a~, <~1/11. The required mesh resolution is determined by the acceptable error in nodal amplification factorS: The guideline that is often employed by numerical analysts solving wave propagation problems by the Galerkin method is 10-20 linear elements per wave. This corresponds to a relative error of 0.1-1% in the amplification factor. Designers of integration algorithms for transient analysis allow errors of 1 / 4 - 1 / 2 % . Based on this experience, the methods for the advection-diffusion problem with a linear source considered here are selected to keep the error below 1/2%. We thus require e~ < 5 x 10-5

(65)

1. Harari, T.J.R. Hughes / Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191

184

4

\

I

\

2

I

Galerkin --

--

\ \

1

GLS --- --

\ 0

g

I

\

3

"\

-1 -2 -3 -4 -5 -6 -6

-4

-2

0

2

4

6

I g( ,/ ~) F i g . 12. R e q u i r e d

mesh resolution.

allowing the mesh resolution presented in Fig. 12. The abscissa here represents the physics of the problems and the ordinate corresponds to the numerics. The bold line marks the transition to propagation. Galerkin/least-squares offers significant relaxation in mesh resolution requirements, particularly approaching the two limiting cases. In the intermediate region, the Galerkin method is more cost effective. Note the logarithmic scale of the plot so that even a small distance between the curves corresponds to a large difference in required resolution. 5. 3. Galerkin ~least-squares/gradient least-squares

The Galerkin/least-squares/gradient least-squares equations for a uniform mesh of linear elements may be obtained from the Galerkin equations by making the replacements O~u O~U <"

2 2 ' 4' 1 + 4a.~'K/h + 144a-kZvK/h

(66)

1 - 12ak~'K/h 2 ' 2 4• ak *--ak 1 + 4 a ~ r K / h z + 144ak'rvK/h

(67)

Since there are two design variables, equations analogous to (53) and (54) can be solved to yield nodal exactness, leading to zK h 2

R

6a k sinh a. - a.(cosh a. - cs) 72a~ sinh a.

(68)

Zv~ _ 6a2k(2a,, cosh a. - 3 sinh a. + csa,,) _ ( h4

2592a ~. sinh a,,

a,, ~2 6a k sinh a.

\ 6ak /

- a.(cosh a,, - cs)

72a ~ sinh a,,

(69)

Note that the second term in (69) contains the right-hand side of (68), somewhat reducing the cost of computing the parameters.

6. Numerical examples The results of several computations are presented in support of the analysis of the numerical methods performed herein, and to validate the design of the proposed methods and their recommended range of

1. Harari, T.J.R. Hughes I Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191

185

applicability. The first set of results is for problems with exponentially varying solutions, at a k = a2,/18, approaching the transition to propagation. Closer to the advection-diffusion limit, the performance of Galerkin/least-squares methods is clearly far superior to that of Galerkin (see Fig. 12), but is less so near the transition. Figure 13 shows results of three cases with increasingly coarser resolution. At this point, the suggested allowable resolution for Galerkin is %, = 0.24 and for GLS, a,, = 0.87. It would appear from the plots that these limits are quite conservative since the Galerkin solution starts to deteriorate only at a,, = 1 and the Galerkin/least-squares solution is accurate in all three cases. The Galerkin/least-squares/gradient least-squares results are accurate at all resolutions as expected. The next set of results is obtained for the same ratio of physical coefficients, but the boundary values are swapped. This changes the results drastically (see Fig. 14). At a,, = 1, there is already considerable degradation in the Galerkin solution, and in GLS at a, = 1.5. It appears that to guarantee accurate solutions in general configurations of practical interest, the allowable limits of resolution should be adhered to. G L S G L S continues to perform well. Results for propagating solutions at a k = 10a~ are presented in Fig. 15. The allowable resolution in this case for Galerkin and GLS is a,, = 0.045 and 0.064, respectively. (Even this seemingly small difference can lead to fewer equations by a factor of three in three-dimensional analysis, which has a

a~ = 0.5

"O-

12 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

I

I

I

'

I

Exact Galerkin -----4::1--- ~ GLS - - o - GLSGLS

- •

-

u

0.2

0.6

0.4

au

1. 0

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

1.0

= 1.0 n

I

I

i

0'9I 0

0.8

0.2

0.4

0.6

'

i

,

__i

0.8 1.0

o~ = 1.5 1.0

I

I

'

I

'

/

0.8 0.6 0.4 0.2 0', -0.2

,

I

0.2

,

I

,

0.4

I

0.6

,

I

0.8

~

1.0

x

Fig. 13. Exponentially varying solutions a~. = a~/18 with go = 0 and gt = 1.

186

1. Harari. T.J.R. Hughes / Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191 a . = 0.5 4.0 3.5 3.0 2.5 2.0 1.5

-~

1.01

0.5 0 0.2

0.4

0.6

a . -40

'

i

I

'

0,8

1.0

1,0

'

I

I

'

35 30 25 20 15 10

5 0.2

0.,I

0.6

0.8

1,0

a~ = 1.5 ,100

i

'

i

'

i

'

i

350 I

300

./ ,~

250 200

.

~

.

150 100

50 0',

~ ,~

~

0.2

Fig.

14.

Exponentially

varying

~

~<"

~

I

0.4

solutions

0.6

ak

0.8

,

1.0

= e~:,/18 with g. = l and g, = 0.

considerable effect on the computational cost.) For a,, = 0.05, there are 12 elements per wave. All methods perform well, although a slight deterioration is already evident in the Galerkin solution. Reducing the resolution to approximately six elements per wave (a,, = 0.10) degrades the Galerkin results further. Again the proposed limits of resolution appear conservative. At a,, = 0.15, there are just under four elements per wave, the minimum number sufficient to resolve the physical phenomenon. Here the Galerkin solution is completely wrong. Some error is evident in GLS as well. Galerkin/leastsquares/gradient least-squares is still nodally accurate as expected, Changing the boundary values (see Fig. 16) again enhances the validity of the proposed limits of resolution. Only the GLSGLS method retains accuracy of solutions as mesh resolution decreases.

7. Conclusions In this work we have developed and analyzed finite element methods for solving problems of steady advection-diffusion with production (1)-(3). If one is interested in solving this problem in homogeneous isotropic media, it is best to perform the transformation via the change of variables (7) prior to

1. Harari. T.J.R. Hughes / Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191

187

a . = 0.05

'

1.0

tExacl

II

I

I

Galerkin ------O-GLS --0--

0.5

m~~l~

. ~,~

~ .

"fl-

-0.5 -1.0

,

0

I

~

•

0.2

,

0.4

I

I

0.6

I

i

0.8

1.0

~ , = 0.10 1.0 0.8 0.6 0.4 0.2 0 -0.9 -0.4 -0.6 -0.8 -1.0

'

0

'

-

S

0.2

'

'

0.4

0.6

'

,.4

0.8

1.0

a~ = 0.15 1.0

'

'

'

'

'

'

.

4

"-

-0.5

/

/

0.5 0

~----

'

/'-../.'>\

,

"W'

'=-"

-1.0 I

0

0.9

~

I

i

0.4

I

0.6

~

I

0.8

v 7

1.0

z

Fig. 15. P r o p a g a t i n g s o l u t i o n s ctk = 1 0 a ~ w i t h go = 0 a n d g~ = 1.

computation. The resulting modified Helmholtz equation engenders no skew-symmetric terms so that standard solvers may be employed. Furthermore, GLS alone provides superconvergence in the one-dimensional case, a high degree of accuracy on general configurations of practical interest, along with proven convergence at optimal rates and a considerable amount of numerical experience. Overall, this is the more accurate and less costly route when available. For designing and analyzing numerical methods that solve (1)-(3) directly, we employ experience from the analysis of Galerkin finite elements and design of Galerkin/least-squares and Galerkin/ gradient least-squares methods, which indicates that it is beneficial to consider a model problem, typically in one dimension, with a known exact solution. Values of exact solutions of a one-dimensional model of advection-diffusion with a source at uniformly spaced points are parameterized by two non-dimensional parameters that contain the physical coefficients and the mesh size. The two nodal amplification factors, characterizing the dependence of exact solutions at discrete points on values at neighboring points, are readily expressed in terms of the two parameters. An analysis of numerical solutions obtained by the Galerkin method employing linear finite elements indicates that relatively fine meshes may be required to retain an acceptable degree of accuracy for certain ranges of values of the physical coefficients. In general, Galerkin representation of propagating

188

I. Harari, T.J.R. Hughes / Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191 a . = 0.05 1.5 1.0 Galerkin -----421---

0.5

GLS - - o ~ ~GLSGLS--~-

0 -0.5

/

-1.0 -1.5 -2.0 0

0.2

0.4

0.6

0.8

1.0

0.8

1.0

a~ = 0.10

"O-

2.0 1.5 1.0 0.5 0 -0.5 -1.0 -1.5 -2.0 -2.5 0.2

0

0.4

0.6

au = 0.15 4

'

~

'

~

-2

'

~ A ~ / q N

.,/'

/

~0 /

~ " .11"

-4

\~.' / 0.2

0.4

0.6

0.8

1.0

Fig. 16. Propagating solutions a~ = lOa~, with g, = 1 and gt = 0 .

solutions contains errors in phase and magnitude. When exact solutions are exponentially varying, errors in exponential rates may eventually lead to spurious oscillations on coarser meshes as one or both numerical amplification factors take on signs opposite to the exact ones. In previous work, the Galerkin/least-squares method was designed to provide exact nodal amplification factors for the two limiting cases. This criterion respected restrictions placed on the method by error analysis and good performance was subsequently demonstrated on a variety of numerical tests. In this case it is not possible to design the method so that both amplification factors are exact. Instead, numerous criteria for improved accuracy are proposed and examined in the entire parameter space. Of these, several are found to provide the lowest error in nodal amplification factors, each in a certain range of ratios of physical coefficients. Defining the method so that the larger of the amplification factors is exact (59) provides the best results in the region in which solutions are exponentially varying. In most of the region that contains propagating solutions, balancing relative errors of the coefficients that determine the dependence of solutions on values at neighboring points (63) provides the best results. In the remaining part of this region, adjacent to the transition zone, results obtained by the Galerkin method are accurate enough so that employing Galerkin/least-squares is not justified. There is a certain deterioration in accuracy in a small region around the transition (where the two exact nodal

1. Harari, T.J.R. Hughes / Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191

189

amplification factors approach each other in value) which may be overcome by employing a G a l e r k i n / least-squares method with (63). Guidelines for required mesh resolutions are presented for Galerkin and Galerkin/least-squares, showing that in large parts of the parameter space substantial savings may be obtained by employing Galerkin/least-squares. G L S G L S , a new variation of the Galerkin method combining GLS with G G L S , is designed so that both nodal amplification factors are exact in the entire parameter space, as in the case of the Galerkin method for the Laplace operator and Galerkin/least-squares for advection-diffusion and Helmholtz. These observations are validated by computation.

Appendix A: Additional design criteria for Galerkin/least-squares Several Galerkin/least-squares methods were proposed in Section 5, each recommended in a different range of the parameters. Of the many methods that were considered, the ones selected are those that provide the lowest error in nodal amplification factors. For reference purposes, the other criteria are listed in the following: (1) For an exact magnitude of the nodal amplification factors, see Eqs. (49) and (50), we require ~-K 1 + 2a k - a,, coth %, h2 4(6a2k--a~)

(70)

An added benefit of satisfying this criterion is that the coefficients in the numerator of (43) also have exact values. This method provides nodal exactness in the advection-diffusion limit. (2) The Galerkin/least-squares method may be designed so that the coefficient of ~bA_~ + ~bA+t in (43) is exact, 1 -- 4or k

CS

1 + 2~ k - cosh %

(71)

namely

1 cosh%,-cs elk = 2 2 cosh a,, + cs

(72)

This method also provides the exact value of the arithmetic mean of the coefficients in (42) Pl P2 + 1 p~ + p ~ -

cosh a,, (73)

cs

Here we require (1

TK

+ 20tk)CS --

(1

- 40~k)

cosh %,

4((6a2k - a])cs + (12a~ + a~) cosh %,)

h2

(74)

which is exact in the Heimhoitz limit. (3) Similarly, the Galerkin/least-squares method may be designed so that the coefficient of ~bA_1 -~b,~÷~ in (43) is exact 4~ k

1 -

-

cs -

%,

(75)

sinh %,

by employing rK

--7 h =

(1 4(a~, + 12a~)

% , c s / s i n h %, -

40¢k)

This method is exact in the advection-diffusion limit. (4) Making the sum of amplification factors pl + P2 exact

(76)

190

I. Harari, T.J.R. Hughes / Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191

1 -4~ k - - exp(a,,)cs 1 + 2C~k-- a,,

(77)

by employing ~'K

e x p ( a , , ) c s ( 1 + 2 ~ k - a,,) - (1 - 4 a k)

h2

k - aT, ) + (12a~ + aT,) )

4(exp(a,)cs(6a-

(78)

yields a method that is exact in both the advection-diffusion and Helmholtz limits. (5) Likewise, the term (Pt + R 2 ) / ( P t P 2 ) may be made exact 1 -

4~ k

cs

1 + 2elk + ~,, - exp a,,

(79)

by employing CS(1 + 2 a k + a,,) -- exp(a,,)(1 -- 4ak)

"r___KK_

(80)

--4(CS(6a2k -- a ~ ) + exp(a,)(12a2k + a ~ ) ) "

h2

This method, too, is exact in both the advection-diffusion and Helmholtz limits. (6) Making the ratio of the non-dimensional parameters c~k and c~,, from Eqs. (53) and (54) exact, ~, --

a,,

cosh -

~,, -

cs

(81)

6 sinh a,,

leads to a method that is identical to Galerkin/least-squares/gradient least-squares with ~'v = 0, i.e. ~- is defined in Eq. (68). (7) A least-squares fit of the coefficients in (43), min

1-4elk

cs

]

+

1--46~

~

/

,

(82)

yields a rather complicated expression for r. This method is exact in both the advection-diffusion and Helmholtz limits. (8) Least-squares fit of the non-dimensional parameters d, and 6,, from Eqs. (53) and (54), ( min

6k

1 cosh_a±-cs'~2 ( 2 2 c o s h a,, + cs / + a "

3sinha,, 2 c o s h a,, + cs

)2 "

(83)

(9) Criteria such as a least-squares fit of the product and the phase, of p~ and P2 or obtaining the exact difference p~ - P 2 require parameters that are too complicated to be defined explicitly.

Acknowledgment

The authors wish to thank Ken Jansen for helpful discussions on turbulence modeling and Greg Hulbert for valuable comments on allowable errors in time-integration schemes. In the course of this work, Isaac Harari was supported in part by the Center for Absorption in Science, Ministry of Immigrant Absorption, State of Israel.

References

[11 C. Johnson. Streamline diffusion methods for problems in fluid mechanics, in: R.H. Gallagher, G.F. Carey, J.T. Oden and O.C. Zienkiewicz, eds., Finite Elements in Fluids--Vol. 6 (Wiley, Chichester, 1986) 251-261. [2] F. Shakib, T.J.R. Hughes and Z. Johan, A new finite element formulation for computational fluid dynamics: X. The compressible Euler and Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 89 (1991) 141-219.

1. Harari, T.J.R. Hughes / Comput. Methods Appl. Mech. Engrg. 115 (1994) 165-191

191

[3] 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) 19-35. [4] T.J.R. Hughes, L.P. Franca and G.M. Hulbert, A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations, Comput. Methods Appl. Mech. Engrg. 73 (1989) 173-189. [5] L.P. Franca and T.J.R. Hughes, Two classes of mixed finite element methods, Comput. Methods Appl. Mech. Engrg. 69 (1988) 89-129. [6] I. Harari and T.J.R. Hughes, What are C and h?: inequalities for the analysis and design of finite element methods, Comput. Methods Appl. Mech. Engrg. 97 (1992) 157-192. [7] L.P. Franca and E.G. Dutra do Carmo, The Galerkin/gradient least-squares method, Comput. Methods Appl. Mech. Engrg. 74 (1989) 41-54. [8] I. Harari and T.J.R. Hughes, Galerkin/least-squares finite element methods for the reduced wave equation with nonreflecting boundary conditions in unbounded domains, Comput. Methods Appl. Mech. Engrg. 98 (1992) 411-454. [9] 1. Harari and T.J.R. Hughes, Finite element methods for the Helmholtz equation in an exterior domain: model problems, Comput. Methods Appl. Mech. Engrg. 87 (1991) 59-96. [10] K. Jansen, Z. Johan and T.J.R. Hughes, Implmentation of a one-equation turbulence model within a stabilized finite element formulation of a symmetric advective-diffusive system, Comput. Methods Appl. Mech. Engrg. 105 (1993) 405-433. [11] I. Stakgold, Boundary Value Problems of Mathematical Physics, Vol. II (Macmillan, New York, 1968). [12] R. Dautray and J.L. Lions, Mathematical Analysis and Numerical Methods for Science and Technology, Vol. 1 (Springer, Berlin, 1990). [13] P.M. Morse and H. Feshback, Methods of Theoretical Physics, Vol. 2 (McGraw-Hill, New York, 1953). [14] T.J.R. Hughes, M. Mallet and A. Mizukami, A new finite element fomulation for computational fluid dynamics: II. Beyond SUPG, Comput. Methods Appl. Engrg. 54 (1986) 341-355. [15] T.J.R. Hughes and M. Mallet, A new finite element formulation for computational fluid dynamics: III. The generalized streamline operator for multidimensional advective-diffusive systems, Comput. Methods Appl. Mech. Engrg. 58 (1986) 305-328.