Finite Element Formulations for Hyperbolic Systems with Particular Emphasis on the Compressible Euler Equations Citation Tezduyar, Tayfun Ersin (1982) Finite Element Formulations for Hyperbolic Systems with Particular Emphasis on the Compressible Euler Equations. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/994z-vj07. https://resolver.caltech.edu/CaltechETD:etd-09112006-130749 Abstract A Petrov-Galerkin finite element formulation for first-order hyperbolic systems is developed generalizing the streamline approach which has been successfully applied previously to convection-diffusion and incompressible Navier-Stokes equations. The formulation is shown to possess desirable stability and accuracy properties. The algorithm is applied to the Euler equations in conservation-law form and is shown to be effective in all cases studied, including ones with discontinuous solutions. Accurate and crisp representation of shock fronts in transonic problems is achieved. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: (Mechanical Engineering) Degree Grantor: California Institute of Technology Division: Engineering and Applied Science Major Option: Mechanical Engineering Thesis Availability: Public (worldwide access) Research Advisor(s): Wu, Theodore Yao-tsu Thesis Committee: Sabersky, Rolf H. (chair) Corngold, Noel Robert Jennings, Paul C. Hughes, Thomas J. R. Wu, Theodore Yao-tsu Defense Date: 10 May 1982 Funders: Funding Agency Grant Number Caltech UNSPECIFIED Stanford University UNSPECIFIED NASA Ames Research Center UNSPECIFIED Record Number: CaltechETD:etd-09112006-130749 Persistent URL: https://resolver.caltech.edu/CaltechETD:etd-09112006-130749 DOI: 10.7907/994z-vj07 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 3460 Collection: CaltechTHESIS Deposited By: Imported from ETD-db Deposited On: 27 Sep 2006 Last Modified: 24 Apr 2025 22:36 Thesis Files Preview PDF - Final Version See Usage Policy. 4MB Repository Staff Only: item control page

Finite Element Formulat ions for Hyperbol ic Systems with Particula r Emphasis on tne Compres sible Euler Equation s Thesis by Tayfun Ersin Tezduyar In Partial Fulfillme nt of the Requirem ents for the Degree of Doctor of Philosoph y Californi a Institute of Technolo gy Pasadena , Californi a 1982 (Submitt ed May 10, 1982) ACKNOWLEDGEMENTS I wish to thank my advisor Thomas J.R. Hughes, for his guidanc e and patienc e through out my graduat e studies. I am especial ly grateful for his assistan ce and encoura gement during the prepara tion of this thesis. This researc h was made possible by the suppor t of the Californ ia Institut e of Technol ogy, the Stanfor d Univers ity and Compu tational Fluid Dynami cs branch of the NASA Ames Researc h Center. I would like to thank H. Lomax for guidanc e and encoura gement and J. Barton, T. Holst and J. Steger for valuabl e suggest ions. iii Abstract A Petrov-Galerkin finite element formulation for first-order hyperbolic systems is developed generalizing the streamline approach which has been successfully applied previously to convection-diffusion and incompressible Navier-Stokes equations. The formulation is shown to possess desirable stability and accuracy properties. The algorithm is applied to the Euler equations in conservation-law form and is shown to be effective in all cases studied, including ones with discontinuous solutions. Accurate and crisp representation of shock fronts in tran- sonic problems is achieved. iv TABLE OF CONTENTS Acnowledgements ii Abstract iii Chapter 1 Introduction 1 Chapter 2 One-dimensional Hyperbolic Systems 5 2.1 Initial/Boundary-value Problem 5 2.2 Weighted Residual Formulation 7 2.3 Semi-discrete Equations 10 2.4 Transient Algorithms 12 2.5 Selection of 16 2.6 Finite Difference Equations Chapter 3 20 Stability and Accuracy Analysis of Algorithms for the One-dimensional Linear Hyperbolic Equation 22 3.1 Development of the Tools for the Analysis 22 3.1.l Introduction 22 3.1.2 Finite Element Solution 23 Unified Analysis of Algorithms 32 3.2.l Description and Classification of Algorithms 33 3.2.2 Implicit Algorithms 35 3.2.3 Explicit Algorithms 37 3.2.4 Stability and Accuracy Studies 40 3.2 Comparison of Algorithms 41 3.3.1 Implicit Algorithms 41 3.3.2 Explicit 1-pass Algorithms 44 3.3.3 Explicit 2-pass Algorithms 46 Numerical Applications in One Dimension 48 4.1 Introduction 48 4.2 Numerical Applications in the Linear Transient Case 48 3.3 Chapter 4 T v 4.3 Chapter 6 Chapter 7 Appendix I Appendix II Propagation of a Small Disturbance in a Gas 48 4.2.2 Initial/Boundary-value Problem 50 4.2.3 Finite Element Solutions 51 Numerical Applications in the Nonlinear Transient Case 58 4.3.l Barotropic Compressible Flow 58 4.3.2 Initial-value Problems 59 4.3.3 Finite Element Solutions 60 Numerical Applications in the Nonlinear Steady Case 67 4.4.1 Isothe.tn1.al Flow in a Nozzle 67 4.4.2 Boundary-value Problem 69 4.4.3 Finite Element Solutions 70 Multi-dimensional Hyperbolic Systems 83 5.1 The Initial/Boundary-value Problem 83 5.2 Weighted Residual Fo.tn1.ulation 84 5.3 Semi-discrete Equations 85 5.4 Transient Algorithms 86 5.5 Selection of T. 86 Numerical Applications in Two Dimensions 88 6.1 Introduction 88 6.2 Computational Domain and Algorithmic Features 92 6.3 Subcritical Case 94 6.4 Supercritical Cuse 97 4.4 Chapter 5 4.2.l l Conclusions 101 The Euler Equations 104 I.l General Principles 104 I.2 Three-dimensional Case 106 I.3 Two-dimensional Case 111 I.4 One-dimensional Case 114 Stability and Accuracy Analysis of Algorithms for the One-dimensional Linear Parabolic Equation 117 vi II.l II.2 II.3 Appendix III Development of the Tools for the Analysis 117 II.1.1 Introduction 117 II.1.2 Finite Element Solution 118 Unified Analysis of Algorithms 120 II.2.1 Introduction 120 II.2.2 Implicit Algorithms 121 II.2.3 Explicit Algorithms 123 Summary Stability and Accuracy Analysis of Algorithms for the One-dimensional Linear Second-order Hyperbolic Equation 129 III.l Development of the Tools for the Analysis 129 III.1.1 Introduction 129 III.1.2 Finite Element Solution 130 III.2 Unified Analysis of Algorithms 137 III.2.1 Introduction 137 III.2.2 Implicit and Explicit 1-pass Algorithms 139 Explicit 2-pass Algorithms 140 III.2.3 III. 3 References 124 Summary 141 144 l CHAPTER l Introduction Analysis of inviscid, compressible fluid flows, especially ones with discontinuities, has been an interesting and challenging part of the research done in the field of computational fluid dynamics (see [L6]). Numerous workers in this field have employed finite difference techniques. The following represents a brief sampling of some recent works. Ballhaus et al. ~B4] used implicit approximate fac~orization schemes to solve the transonic small disturbance equation. Holst and Ballhaus [H3] applied approxi- :rnate factorization schemes to the full potential equation in conservation form. Holst and Brown [H4] utilized solution adaptive grids for the full potential equation in conservation form. Warming and Beam [W3] used approxi- mate factorization schemes to solve the Euler equations in conservation law form. Steger [S3], and Steger and Warming [s2] applied flux vector split- ting ideas to the solution of the Euler equations in conservation form. Finite difference schemes of the above type are mostly limited to problems with simple geometries. Finite element methods, on the other hand, can easily handle arbitrary geometries. In the finite element method, the problem domain is discretized into sub-domains (elements) , and, via a weighted residual formulation, the governing differential equation system is translated into a system of ordinary differential equations. 2 In a weighted residual formulation, selecting the weighting functions from the same class that the interpolation functions are selected from, leads to a {Bubnov) Galerkin formulation. When applied to differential equation systems with symmetric operators (e.g. diffusion equations, most structural and solid mechanics problems) Galerkin formulations produce solutions with a "best approximation" property. That is, the error is mini- mized with respect to a certain norm. For systems with non-symmetric operators (e.g. first-order hyperbolic systems), however, the Galerkin formulation does not possess a best approximation property. This, in some cases, may result in solutions with spurious node-to-node oscillations. In fact, this problem is not limited to Galerkin finite element formulations. It also arises for finite difference schemes when non-symmetric operators are approximated centrally. Instead of using weighting functions which lead to a Galerkin formulation, one can employ a Petrov-Galerkin formulation by modifying those weighting functions according to an optimal rule. The basic idea is to minimize the spurious oscillations without introducing excessive diffusion to the solution. An optimal streamline upwind/Petrov-Galerkin formulation for convection dominated flows was recently developed by Hughes and Brooks {see [B7, BS, Hl2, Hl4, Hl5]) and was successfully applied to the solution of advection-diffusion and incompressible Navier-Stokes equations. In this work we present a Petrov- Galerkin algorithm which is a generalization of the streamline upwind/PetrovGalerkin algorithm to hyperbolic systems. The weighting functions (which would normally lead to a Galerkin formulation) are perturbed by the product 3 of the coefficient matrix of the hyperbolic system, the gradient of the weighting function and a time parameter. The alternatives of transposing or not transposing the coefficient matrix in the weighting function, and the selection of the time parameter are among the subjects discussed here. The algorithm presented, under some very special conditions, reduces to the Lax-Wendroff scheme, which is known as a shock capturing algorithm. By incor- porating the coefficient matrix of the hyperbolic system into the weighting function we automatically inject the eigenvalue/eigenvector information of the system into our finite element formulation. In chapter 2 we briefly review the properties of one-dimensional hyperbolic systems and introduce the Petrov-Galerkin algorithms. of the weighting functions is discussed in detail. The selection We also investigate under what circumstances the weighted residual formulation of a system can be duced to that of uncoupled single degree-of-freedom equations. re- The proce- dure of finite element discretization, and the transient algorithm used for solving the semi-discrete equation are described. Further, for a special case, we write the finite difference equations for the Petrov-Galerkin formulation. In chapter 3 we perform a detailed stability and accuracy analysis of algorithms for the linear one-dimensional hyperbolic equation. rithms ~f Several algo- interest are studied and compared. Chapter 4 reports nu~erical resul~s in one space dimension. Several linear and nonlinear, steady and transient problems are solved using various techniques. Special emphasis is placed on problems with discontinuous solu- tions (shocks) . 4 In chapter 5 we introduce the multi-dimensional versions of the PetrovGalerkin algorithms. Chapter 6 covers numerical applications in two space dimensions. A biconvex thin airfoil problem is solved for subsonic and transonic cases. Several algorithms are tested. In chapter 7, we draw our conclusions\and make suggestions for future research. Appendix I reviews the properties of the compressible Euler equations. In appendix II, a stability and accuracy analysis of algorithms for the one-dimensional, linear parabolic equation is performed. In appendix III a similar analysis for the one-dimensional linear second-order hyperbolic equation is performed. The methods used in appendices II and III are essentially the same as that used in chapter 3. 5 CHAPTER 2 One-dimensional Hyperbolic Systems Initial/Boundary-value Problem 2.1 D = ]0,L[ Let Q [O, denote the open interval of length 0 and IT t E [O, T] x E and is f = {O, L}, that is, Spatial and temporal coordinateE are denoted by L. the points D The boundary of denote its closure. L] and let L , respectively. , Consider the following systei-a of u -,t + AU ~-,x m partial differential equations: + G = (2.1.l) 0 where u U(x, t) (2.1.2) A = A(U, x, t) (2.1.3) G = G(U, x, t) (2.1.4) and a comma denotes partial differentiation. We are concerned with the case in which (2.1.l) is hyperbolic, that is when S A has real eigenvalues and there exists a transformation matrix such that (2 .1.5) where A is the diagonal matrix of eigenvalues of A (2.1.l) is called a balance law if there exists a vector ~ = $(U, x, t) (2.1.6) 6 such that (2.1. 7) A balance law in which G = is said to be a conservation law. 0 In the linear case (2.1.8) A (x, t) A and G and /(x, B and A In the constant-coefficient case + t)U B(x, (2.1.9) t) are independent of x t Classical references for the study of hyperbolic systems are Courant - Hilbert [c3] and Courant - Friedrichs [C2]. Consideration of the eigenstructure of enables the specification For a general treatment of this topic of appropriate boundary conditions. see Yee [Yl]. A For the present purposes, it suffices to assume that the boundary conditions take the abstract form d u where g(t) = d is a boundary operator and g (2.1.10) is a prescribed function. The initial/boundary-value problem for (2.1.l) consists of finding a function U which satisfies (2.1.1), the boundary conditions (2.1.10), and the following initial condition: U(x, where ~O is a given function of ~o(x) 0) x E ~ . (2.l.11) 7 2.2 Weighted Residual Formulation Consider a discretization of 2 , ... , net , where Q into element subdomains is the number of elements. ne.Q, n u = De e l I We assume ne (2.2.1) e=l (2.2.2) All functions considered in the finite element formulation will be smooth on the element interiors (i.e. Qe's). Two classes of functions are impor- tant in the developments which follow. The classes are distinguished by their continuity properties across the element boundaries. Functions of the first class are assumed to be continuous across element boundaries. These functions are denoted by c 0 = 0 c (Q} and may be recognized as containing the standard finite element interpolations. Functions of the second class are allowed to be discontinuous across element boundaries and are denoted by c- 1 = c- 1 (Q). A weighted residual formulation of (2.1.l) is given by 0 where = 1 W• ( Q - U -,t + AU + G) dQ --,x W is a weighting function and cases we assume U interpolations. denotes the dot product. is approximated by standard, c0 , In all finite element The weighting functions may be selected from a different set of functions than the trial solutions. Petrov-Galerkin formulation (see e.g. Hl5, M2, Rl, Wl]). (2.2.3) Thus (2.2.3) gives rise to a [Bl, B7, B8, Cl, Dl, Hl, H8, Hl2, Hl4 8 An important class of Petrov-Galerkin methods, which is emphasized in the se~uel, = W where and W is defined by the following expression for W + (2.2.4} T W - -,x W is a member of the same class of functions as the trial solutions T is either TA or TAT where T is a parameter which is chosen to optimize accuracy according to some criterion. This class of methods represents a generalization to hyperbolic systems of the streamline-upwind/ Petrov-Galerkin formulation which has been successfully applied heretofore to the advection-diffion and incompressible Navier-Stokes equations [B7, BS, Hl2, Hl4, HlS]. Both choices of T have interesting consequences. assume the linear, constant-coefficient case in which For example, G 0 Choose Then (2.2.3} reduces to the canonical form 0 .J"<~ + T~~,x) • <~,t + ~~.x)dQ (2.2.5) Q where w == STW and u = s- 1 u . of uncoupled scalar equations. Thus (2.2.5) is equivalent to a system Scalar equations of this form are extensively analyzed in Chapter 3. Furthermore, the choice T = TAT leads to difference equations which, under special circumstances, have essential features in corrunon with the well-known Lax-Wendroff method [R2]. Under the circumstances which led to (2.2.5), choosing T = TA does not result in the canonical form (2.2.5) unless the weighted residual formulation is generalized to 9 O = ~(~ + T~~,x) • Q(~,t + ~~.x)dQ (2.2.6) Q where (2.2.7) = w = s w (2.2.8) u = s u (2.2.9) Computational experiences with generalizations of formulations of this type proved''cumbersome and unreliable in the nonlinear regime when compared with (2.2.3) and (2.2.4), and thus were abandoned. Despite the fact that T =TA does not canonically reduce (2.2.3), it leads to another optimality property which will be described subsequently (see §2. 5) . If T is taken to be zero then we have the usual Galerkin method which possesses central-difference like character. 10 2.3 Semi-discrete Equations Spatial discretization of the weighted residual equation (2.2.3) via finite elements leads to the following semi-discrete system of ordinary differential equations: . + Cv = F Mv where M = M(v, t) (2.3.1) is the generalized "mass" matrix, generalized convection matrix, F == F(v , t) the vector of (unknown) nodal values of U , C == C(v, t) is the force vector, is the v is· and a superposed dot denotes time differentiation. The initial-value problem for (2.3.1) consists of finding a function v(t) v satisfying (2.3.l) and the initial condition v (0) where :::: ~o (2.3.2) is determined from (2.1.11). ~o The arrays in (2.3.l) are assembled from element contributions: (2.3.3) e m e ~ab [me J f (Na I + Na,x TT)N -ab = ~"Ge (2.3.4) - b cill (2.3.5) ne,Q, c A (ce) (2.3.6) e=l (2.3.7) 11 e ~ab (2.3.8) F (2.3.9) fe {fe} (2.3.10) ~a fe -a where = ~ n n en e •e e e + c (m -ab ~b) -a b 2°b 2: b=l Na G dQ (2.3.11) /).. represents the finite element assembly operator; (local) element node numbers; 1 ~a , b < n - where en N function associated with node m x m e 2b a ; is the I en and b are is the number is the element shape of nodes for the element under consideration; a n a identity matrix; and is a vector which contains the boundary condition data emanating from (2.1.10). The dimensions of the nodal arrays and the dimension of and e ~ab and e ~ab are m x m , are The reader is reminded that n np !I<~· t) = L A=l where A v ~A NA (x)~ (2. 3. 12) (t) A is the global nodal index, n np is the total number of nodes and refers to the components of X associated with node A. 12 Transient Algorithms 2.4 We consider first a family of one-step.implicit methods defined by + a M -n+Y -n+y + v ~n+l c -n+Y -n v ~n+y (2.4.1) F -n+Y lit a -n+a (2.4.2) where M = ~(Y:n+y , tn+y) (2.4.3) c -n+y = ~(~n+Y , tn+yl (2. 4.4) F - ~(Yn+y , tn+y> (2.4.5) -n+y -n+y a n+y + y~n+l (2.4.6) ~n = (1 - Y)v + YYn+l (2.4.7) a = (1 - a)a + aa -n+l (2.4.8) + ytn+l (2.4.9) -n+a t n+y In the above, and - Y)a v -n+Y a (1 ~n ~n (1 - y)t 6t n is the time step, n is the step number, and y are parameters which determine stability and accuracy properties. The starting value, ~ 0 , ~o c:i:o may be determined from So Yo (2.4.10) 13 where ~o = ~(~O' 0) (2.4.11) ~o = ~(~O' 0) (2.4.12) ~(~O' 0) (2.4.13) ~o y If =1 the above algorithm reduces to the generalized trapezoidal y method, whereas if =a , i t reduces to the generalized midpoint method. These methods have been contrasted in [H7, Hl6]. A general family of predictor/multi-corre ctor algorithms, based on the preceding implicit methods, is implemented as follows: 1. i 2. v (0) -n+l (2.4.14) is the iteration counter) i = 0 = v-n + b.t(l - a)a-n (2.4.15) (predictor phase) (0) (2.4.16) 3. ~n+l 4. R = F (i) -n+y 5. M* 6a = R 6. ~n+l :::: ~n+l = 0 (i+l) ( i) M(i) -n+y ~n+y (i) M* (i) + (i) ~n+y ~n+y (residual force) (2.4.17) ( 2. 4 .18) is the "effective mass") (2.4.19) 6a (corrector phase) 7. (i+l) ~n+l = (i) ::'.'n+l + a Lit Lia (2.4.20) If additional iterations are to be performed, i + 1 , and calculations resume with step 4. i is replaced by Either a fixed number of 14 iterations may be performed or iterating may be continued until a convergence condition. at step n + 1 (i+l) '.:n+l = ~n+l ) . R satisfies When the iterative phase is completed the solution (i+l) is defined by the last iterates (i.e. v = v and -n+l -n+l At this point n is replaced by n + 1 and cal:ulations for the next time step may begin. The properties of the algorithm are strongly influenced by the choice of the effective mass. There are various possibilities. For example, a fully implicit procedure may be defined by taking M* M(i) + a6t c(il + a6t H(i) -n+Y -n+y -n+y (2.4.21) where M(i) -n+Y c (i) -n+Y H (i) -n+Y ( ( i) - --n+Y tn+y> (2.4.22) C(v(i) - -n+Y , tn+Y) (2.4.23) H(v(i) - -n+y I tn+Y) (2.4.24) == M v = = I ne.Q, H = A (!:e) (2.4.25) e=l and = (2.4.26) = (2.4.27) has dimensions m x m . In general, this definition of leads to a non-symmetric band-profile matrix. M* 15 An explicit algorithm may be constructed by taking M* to be "lumped" (i.e. diagonal) : M* = (2.4.28) -Md.iag There are several schemes for obtaining suitable -Md.J.ag In the present work we assume that the diagonal element array is defined by e m a I a if [Zl] =b (2.4.29) 0 if a 'I b where (2.4.30) = (2.4.31) In the present work we confine our attention to the implicit and explicit schemes defined above. Stability and accuracy analyses are presented in Chapter 3. However, there are\other possibilities. Implicit-explicit finite element mesh partitions [H9, HlO, Hll, Hl3] may prove useful, for example. Additionally, to obtain the stability properties of implicit methods, while eliminating the equation-solving burden imposed by (2.4.21), approximate factorization schemes may be employed. We are presently experimenting with element-by-element factorizations which are very convenient from an implementational standpoint [ Hl 7]. 16 2.5 Selection of T T Two expressions for have been employed in the numerical calculations. One is based upon spatial discretization and the other upon temporal discretization. They are given as follows: spatial criterion In this case we assume = F a h/p T where F A , p A. (A) 's and the 1 T is the element length, p that is = max 1 < i < m are the eigenvalues of IA. (A) I 1 A . - Note that (2.5.1) (2.5.2) is a local - specification of values of h is a non-dimensional parameter, is the spectral radius of (2.5.l) A T in that i t depends upon the element lengths and eigen- which vary from point to point. Rationale for this form of is provided by the following exar.iples: Examples 1. Consider the scalar model equation u, t where A + A u ,x = o (2.5.3) Raymond and Garder [Rl] have shown that is assumed constant. if F a (2.5.4) 17 then the semi-discrete equations achieve fourth-order phase accuracy. Consider the steady analog of (2.5.3) regularized by a diffusion 2. term, A. u ,x As s + 0 , = E: u,xx (2.5.5) the choice (2.5.6) leads to nodally exact solutions. The general case for the advection- diffusion equation is described in Hughes-Brooks [B7, BB, Hl2, Hl4, Hl5]. Remarks 1. (2.5.4) and (2.5.6), need to The preceding optimality conditions, be altered for higher-order elements. For example, in the case of three- node quadratic elements (2.5.6) should be changed to F a = 114 [Nl] Throughout this work only low-order elements are employed. 2. A weighted residual formulation of a linear, constant-coefficient, hyperbolic system with G = 0 can be given in which each uncoupled scalar equation is treated optimally, viz. 0 = f ~ where W == ST w and (W + F a h sgn ~ ~,x> u = s-l u cf. • (2.2.5). (~,t + ~ ~,x>cill (2.5.7) Unfortunately, there does not appear to be a nonlinear or multidimensional generalization of (2.5.7) and consequently, we have not explored the subject further. 18 temporal criterion In this case we assume T = F a 6t (2.5.8) Note that (2.5.8) is a global specification in that all elements. ~t is the same for Rationale for this choice is provided by the following examples. Examples 1. Assume H = 0 . If T = TA , symmetry of the implicit operator and F = l , M* (see Eq. 2.4.21). from the definitions of the element contributions to me + a~t e ~ab -ab !Na = J ~~ I + a~t Qe + (a~tl 2 then (2.5.8) leads to This can be seen M* (Na,x Nb Qe ~Na,x Nb,x ATA dQ = if AT + N N a b,x ~)dQ )T (~a + a6t Ce -ba (2.5.9) The obvious advantage in this case is the decreased storage and factorzation costs. Symmetric element arrays are also advantageous in implicit- explicit finite element mesh partitions [Hl3]. This choice also leads to an optimality condition in that for a specified residual, R , 6a is optimal with r:spect to the norm defined by ting this is to say that the increment of the synunetric bilinear form which generates U M* . the increment Another way of put- is optimized with respect to M* . This concept of optimality is related to the following optimal steady formulation 19 a(W, U) where a(• , •) A 2. and G (2.5.10) is a symmetric, bilinear form defined by = and -f (~ + T A W-,x ) • G dft st == fst T - l (W - + T A W -,x ) • (p + T A U -,x are assumed to be independent of (2.5.11) that is A = A(x) (2.5.12) G = G (x) (2.5.13) F = 1 , The choices u , ) a.It leads to an explicit Lax-Wendroff type method. and M* = M . - -diag I We shall explore this point further subsequently. Remark The factor, F , for nonlinear effects. in (2.5.1) and (2.5.8) has been included to account It has been our experience that a value of F greater than one needs to be employed to adequately handle shock-wave phenomena. 20 2.6 Finite Difference Equations The finite difference equations for the preceding algorithms are needed subsequently for the stability and accuracy analyses and are of interest in their own right. In explicating the finite difference equations for an inter- nal node we have made the following assumptions: (i) linear elements are employed; and (ii) h , over each element. Furthermore, for notational clarity we have dropped the superscript i A and H are constant; and (iii) n + y. and subscript G varies linearly The equations are as follows: implicit case .+ (~<~ + atlt H)D r + Mtl h T T A (- TT + ailt A - ex.lit TT ~)Dl D 2 )tl~(j) = - Dl)~(j) '( - -h I Dr + TT 2 ~ 2 + + ( - A D1 - h h· T ) ( - 2 I Dr + T Dl ~ (j) (2.6.1) where D r and v (j) - 2(r ~(j-l) + (1 - 2r)~(j) + r ~(j+l)) (2.6.2) 1 (2.6.3) Dl ~ (j) = 2<- ~(j-1) + ~ (j+l)) D2 v (j) = - ~ (j) number j 1 2 -(j-1) -(v - 2 v (j) + "'.:1:(j+l)) = v(x.) is the subvector of , The value of - -J etc. rule employed. r v (2.6.4) which is associated with node is determined by the element quadrature The following are the most important cases: 21 r rule 1/4 1 - point Gauss 1/6 2 - point Gauss (exact) 0 trapezoidal explicit case In the explicit case, the right-hand side of (2.6.1) is the same, but the left-hand side simplifies to (2.6.5) Remarks 1. It is interesting to observe that even though upwind influence has been introduced via the weighting function defined by (2.2.4), the resulting difference equations are centered about node 2. G = 2 , Assume a = l T = TA T , j T = Llt/2 pass (i.e. one-iteration) case, then from (2.4.15), and the explicit one- (2.4.16), (2.4.20) and (2.6.1) we get ~n+l (j) Eq. = ( ~ - L\t h A)2 D2)v-n (j) ( L\t h - (2.6.5) defines the Lax-Wendroff method. (2.6.5) 22 CHAPTER 3 Stability and Accuracy Analysis of Algorithms for the One-dimensional Linear Hyperbolic Equation 3.1 3.1.l Development of the Tools for the Analysis Introduction Model Problem: Convection Equation One-dimensional convection of a function U(x,t) is governed by the following hyperbolic equation: u ,t + \ u ,x = o where \ is the convection velocity. (3.1.1) An initial condition of the following form is assumed: U(x,O) = e where i = (-1) ~ and k ikx is the wave number. is an element of the set of Fourier functions. (3.1.2) We note that the function 8 ikx This is a complete set and any piecewise regular function in the range of (0,2TI) can be represented as an expansion in this set. Exact Solution The exact solution of (3.1.1) and forward. (3.1.2), for constant A, is straight~ Assuming a solution of the form: U(x,t) = X(x)T(t) (3.1.3) 23 leads to: X(x) = e ikx T (t) = e ~ We define the damping coefficient, Vt - i \) (3.1.4) , (3.1.5) A. k (3.1.6) and the frequency, w , as the com- v ponents of the complex parameter (3~1.7) \) Thus: ~ = w 3.1.2 0 (3.1.8) A. k (3.1.9) Finite Element Solution Semi-Discrete Equation The Petrov-Galerkin formulations described in chapter 2 leads to the following semi-discrete equation: . (3.1.10) Mv+cv=O where v ( 3 .1.11) = {v.} J (3 .1.12) Here j stands for an interior node, j that node. and (2 . 3 . 6 ) . The matrices M and C , and x. J is the coordinate of were previously defined by (2.3.3) 24 For the purpose of analysis, we assume constant spacing h . A and constant mesh Further, we assume that the finite element solution, also has a separable form: =XT (3 .1.13) T T (t) (3.1.14) x ~x .} (3'.1.15) v where - J (3.1.16) Consequently, we get the following semi-discrete form: . MX T + C X T =0 (3.1.17) The spatial component of the numerical solution is determined by the nodal interpolation of the imposed initial condition (3.1.2): X . -_ e ikxJ· -_ eiqj J (3.1.18) where the dimensionless wave number, (3.1.19) q = k h ' is a measure of the spatial refinement of the numerical method. The jth equation of the system of equations of (3.1.17) is I(% ~r - T A e1)• I + (A ~l + TA x.J- 1 2 2 h e2H X. J X.J+ 1 0 (3.1.20) 25 where T is the parameter that appeared in section 2.2, and are the stencils for node j ~r , ~l ~ , 2 corresponding to the assembly of the element level matrices: - . o-r (3.1.21) ~l (3.1.22) +l f Na ' ~ - Nb, t; d t: (3.1.23) -1 These stencils are directly related to the difference operators given by (2.6.2)-(2.6.4). Depending on the numerical integration technique used, they can assume several fonns. D -r = For example, 2[1/6, 4/6, 1/6] e1 = [-112, a, 1121 Q2 = (exact integration) (3.1.24) (exact integration) (3.1.25) (3.1.26) [-1/2, 1, -1/2] Further, we define the following array: E = [ e -iq , 1, e+iq] (3.l.27) Then [X. l' X., X. J JJ ]+ 1 = [e iq (j-1) = e iqj E , e iqj e iq(j+l)l J (3.1.28) 26 Substituting (3.1.28) into (3.1.20) leads to (h D 2 ~r (3.1.29) 0 Now define the complex scalars M and C corresponding to M and c respectively: M 1 2TA - = -(D 2 ~r h D )E ~l (3.1.30) ~ (3.1.31) i'.ltC where 6t is the time step of the time integration algorithm. We define the following non-dimensional parameters. Cllt = C2T = is called the Courant number. time, we can view 2T . c6t If we set All.t (3.1.32) h 2TA (3.1.33) h Considering that T has units of as an algorithmic "Courant number" based on c2 T = O , we obtain the usual Galerkin formulation. With the definitions of (3.1.30) and (3.1.31), (3.1.29) reduces to the following ordinary differential equation: . + CT M T 0 (3.1.34) Transient Algorithm Transient algorithms were described in section 2.4. We adopt the repre- sentation: y . [TI T] (3.1.35) 27 and let Y _n Given denote the approximation of v _n ~ to calculate I th . -time step. :n+l· y ( 0) -n+l P n we go through a predictor phase and an iterative phase In the predictor phase we calculate where at the Y y (0) -n+l by the following operation: p y (3 .1.36) - _n is the predictor matrix defined by l (1-a) llt p == 0 (3.l.37) _J The iteratiqe phase starts with the zeroth-iteration value, Y(O) -n+l ' and continues according to the recurrence rule below: Given solve the following system for ML .6.T<i> n+l + CL .6.T(i) n+l .6.T (i) n+l = a.6.t .6.T (i) n+l • (i+l) Tn+l = • (i) Tn+l + • ( i) .6.Tn+l (3 .1.40) T(i+l) n+l = T(i) n+l + .6.T (i) n+l (3 .1.41) The superscripts (·3 .1. 38). and yCi+l) . ,;,n+l · c L and = +CR T(i)) - (MR T(i) n+l n+l (3.l.38) (3.1.39) R refer to the left and right-hand sides of We reserve the option of having different evaluations for M on different sides, so that we can accomodate all of the algorithms described in chapter 2. 28 The recurrenc e rule defined above can be expressed as: v (i+l) J J y (i) = ~n+l (3 .1.42) - -n+l is the iteration matrix emanating from the recurrenc e rule: 1-a fltCR M = J (3.1.43) 1 where M Combining the predicto r (3 .1.44) and iterative phases, we have y A y -n+l (3 .1.45} - -vn (3 .1.46) Here S is the number of iteration s., Exploitin g the ·fact that det P = 0 , that is det JS det p det A = 0 (3.1.47) we can write A21 A22 All Al2 = J.1 (3.1.48) Then, from (3.1.45): T n J.1 Tn . T n+l J.1 Tn+l (3 .1.49) 29 Substituting (3.1.49) into (3.1.45): T (tr A)T n+l - ( 3 .1. 50) n Numerical Frequency and Damping Coefficients For the temporal component, we assume a solution of the form: T n where v = e v tit n (3.l.51) is the numerical counterpart of the exact v defined by (3.1.6). From (3.1.50) and (3.1.51): e Now, we need to calculate tr A :::: tr(J - tit \) = tr A tr A s P) s + (1 - a).6t s = Jll J21 are components of where (3.l.52) JS (3.l.53) can be calculated by way of the Cayley-Hamilton theorem: (3.1.54) The coefficients of the matrix a and b \ 1 and \ 2 , J (3.l.55) a b are functions of the eigenvalues, :::: The eigenvalues are given by (3.l.56) 30 .A (3.1.57) l (MR+ a 6t ·cR)/M (3.1.58) l + ~ ( ( 1 - R) S - 1) (3.l.59) l = = 1 - From (3.1.53)-(3.1.58): tr A = where (3.l.60) (3.1.61) One can note that, from (3.1.30)-(3.1.33), the terms expressed in terms_ of the dimensionless parameters M= l 2 q M and 6tC can be I c6t and c2T (D (3 .l. 62) -r (3.l.63) From (3.1.52) and (3.1.59)-(3.1.63), we can express the complex parameter v6t in terms of the dimensionless parameters q I c6t and c2T , V6t = ln(tr A) For comparison, we need to express less parameter set. v6t that is, ( 3 .1. 64) in terms of the same dimension- From (3.1.6): v6t = - iAk 6t = - i A.6t h kh (3.1.65) By means of {3.1.64) and (3.1.65) we can study the stability and accuracy of a wide variety of algorithms for solving the convection problem. 31 We define the analytical and numerical amplification factors Z and z bv z = e z- = e \)f:. t vt:.t ( 3 .1. 66) (3.1.67) 32 3.2 Unified Analysis of Algorithms Introduction In section 3.1 we developed the tools for stability and accuracy analysis of a large class of algorithms. type. The algorithms can be of implicit or explicit Unification of algorithms into one class enables us to perform one general stability and accuracy analysis for the entire class and then study individual cases within the framework of this analysis. For the purpose of analysis, we classify the algorithms we study into two groups: the implicit types are the base algorithms and the explicit types are the ones derived from the base algorithms. To each implicit algorithm defined, we can (at least in principle) associate an explicit algorithm. This concept was described in section 2.4. Further, we unify all the algorithms we study into one general PetrovGalerkin class. Two main parameters c2 T and r (defined in section 2.6) determine the particular algorithm in this class. We obtain closed form expressions for the modulus frequency (jzj = e-~ 6 t) and the (w 6t). These expressions are simple for the implicit and explicit 1-pass algorithms, and somewhat more complex for the explicit 2-pass aJgorithms. Unconditional-stability proofs can easily be made for the implicit algorithms; stability limits can be determined for the explicit 1-pass alga~ithrns with the same ease. Expressions for the exact and numerical frequencies are not in easily comparable forms. One can expand these expressions in q or c6 t and compare them in series form; this requires extra caution and patience in the algebraic bookkeeping. Alternately, one can, with a general-purpose program, compute 33 w/w the ratio for the desired ranges of c 2 T , C~t and q . The latter approach is adapted herein. 3.2.1 Description and Classification of Algorithms In the following sections we briefly describe the (implicit) algorithms and place them in the general class. (Bubnov-) Galerkin Algorithms In this group we study two Galerkin (c T = 0) algorithms; we call these 2 GC and GL GC is a galerkin algorithm· with consistent (exactly integrated) mass. GL, on the other hand, has a lumped mass (integrated with nodal trapezoidal rule). The stencils D _r and !?1 are !?r = 2· ( r, 1 ~l = The parameter GC and to r , O for [ -~, - 2r, r] (3.2.1) 0, ~ (3.2.2) which was described in section 2.6, is set to 1/6 for GL . The form of the stencil D -1 corresponds to exact integration and is equivalent to central differences. Petrov-Galerkin Algorithms We study four Petrov-Galerkin algorithms. They are PG(Pade), PG(C T = 1), 2 PG(C2T = 2//15) and PG(C2T = c~t). PG(Pade) has weighting function element shape functions: W constructed from the following 34 N1 ( ~> = o (3.2.3) N (~> = i (3.2.4) 2 where ~ (3.2.5) E [-1, +l) is the usual isoparametric coordinate. In this case we need to replace 2[~, ~, o] and [-1, +l, o], (3.2.7) respectively. This method corresponds to the Pade finite difference approximation (see [ W3]) • The other methods, that is [Hl5) , 1) PG ( C2T = 21.fls ) [ Rl J and PG(C2-r = Ctit) employ (3.1.24), rules. (described in chapter 2, section 5) (3.1.25) and (3.1.26), corresponding to exact integration The name of each algorithm implies the way c T is chosen. 2 35 .Implicit AlgorLthms 3.2.2 From (3.1.62)-(3.1.63), the scalars M and M l - 2r (1 - cos q) - Lite become C2T sin q 2 • .1. (3.2.8) Lite (3.2.9) The stencils of the Bade approximation, (3.2.6) and (3.2.7), lead to, l M 2 (1 + e -iq) (3.2.10) L1tc The same expressions can also be obtained by setting Remark: and (3.2.11) c2T = 1 in r = 1/4 Therefore, for the purpose of analysis, (3.2.8)-(3.2.9). we can include the PG(Pade) algorithm in the two-parameter ce 2 T, r) family of Galerkin/Petrov-Galerkin algorithms. For the implicit class, (J.1.52), (3.1.59) and (3.1.67) reduce to the following expression for the numerical amplification factor: ~ z where Q 1 Q (3.2.12) is given by (3.1.61). Frequency Analysis The following expression for tan(w 6t) w (;.t is found: 2 cl:it W(G + c2T V/2)/ (H +a c6t W) (H - (1 - a)Ci::.t W) J (3.2.13) 36 where v 1 - cos q (3.2.14) w sin q (3 .2 .15) G 1 = = H v (3.2.16) - C2T W/2 (3.2.17) 2 r Eq. (3.2.13) is a general expression for any combination of a . For example, for the trapezoidal rule ta.'1. (w l.lt) = r , c2T and (a = ~) the expression reduces to (C T = 0), 2 then (3.2.18) becomes: 2 cl.lt W(G + c2T V/2)/ Further, if we have a Galerkin algorithm = c6 t W G/[G 2 - tan(w ~t) (~ CL'it W) 2 ] (3.2.19) 2 [ 1 - (~ cl.lt w) J (3.2.20) Lumping the mass term reduces this expression to: = tan(w L'it) c.6t w / Modulus/Stabi lity Analysis The following inequality needs to be satisfied for the stability of an algorithm: Izl 2 < (3.2.21) i This inequality translates to: - ( 1 - 4 r) c 2 T v - (2 a - 1) ( c~ T v + ( 2 - v) ) < o (3.2.22) 37 As we can easily observe, for a ~ ~ , are unconditionally stable, provided that GL(C2T = 0, r = 0) have no modulus error, that is, !z! = 1 . GC(C2T 3.2.3 = 0, r = 1/6) , all the algorithms considered r 2_ 1/4 . and In particular, PG(Pade) (C2T = 1, r = 1/4) Explicit Algorithms For explicit algorithms, the scalar quantities MR MR , 6tCR and M are: = G+ i H (3.2.23) 6tCR = c6tc2T v + i c6t w M :::: 1 (3. 2. 24) . (3.2.25) Explicit 1-Pass Algorithms For 1-pass algorithms, the numerical amplification factor takes the form Z = 1 - 6tcR (3.2.26) that is (3.2.27) For the frequency analysis, we get a relatively simple expression: tan (W 6t) Clearly, for 1-pass algorithms, (3.2.28) r and a have no effect. For modulus/stability analysis, the following inequality is considered: < l (3 .2. 29) 38 Utilizing the identity w2 = 2v - v 2 I this condition can be written as: > 0 (3.2.30) We can further wri· te th·is inequality · in two ot'!er forms: The first form , + (1 - v2 + (3.2.31) provides us with the stability limits on the algorithmic parameter The second form, < (3.2.32) provides the stability limit on the Courant number once we select the parameter c2 T . Consider the following cases: If c2 T = 0 , then, for the stability condition we get v - 2 > O . This implies unconditional unstability. If we set c2 T = 1 , By setting c 2 T = C~t , then, the stability condition is c 6t ~ 1 . we get the same stability condition: c6 t < 1 . Explicit 2-Pass Algorithms For the 2-pass algorithms, the numerical amplification factor is ~ z (3.2.33) 39 ~ With this form of Z , the expressions for condition become rather complicated. we ~an, tan(w 6.t) and. the stability Defining the variables = (3. 2. 34) = (3.2.35) without any further algebraic elaboration, write the expressions for tan(w 6.t), that is, tan(w L\t) = AlB2 + A2Bl l + A B - A2B2 1 1 (3.2.36) and for the stability condition (3.2.37) 40 3.2.4 Stability and Accuracy Studies As the algebra of the analysis gets lengthier, especially for the 2-pass algorithms, we find it more convenient to conduct the analysis numerically for the desired ranges of the parameters involved. For the algorithms previously considered we utilized a general purpose program to determine the stability and accuracy pro2erties. algorithm, the parameters w/w and r and c 2T were given. The quantities were computed and plotted for the values of 0 • 8 , 1. 0 and q E ]0, 'IT [ For each c6t ~/w 0.2, 0.4, 0.6, • For future reference, we describe the following concepts: Unit CFL condition [M3]: An algorithm satisfies the unit CFL condition if it produces nodally exact solutions for Order of accuracy: algorithm as q O , The behaviors of the ~/w and w/w curves for an reveal the order of accuracy of that algorithm. ther of these curves has a finite slope as order accurate. =1 c t 6 If ei- q + 0 , then the algorithm is first If both curves have slopes approaching zero as q + 0 , then the algorithm is at least second-order accurate. The quantity ~/w is called the algorithmic damping ratio (see [H2]) and is related to the logarithmic decrement 0 where T = = 2Tr(~/W) (3.2.38) is the numerical prediction for the dependent variable. The dimensionless wave number Scaling 6 [H2] via the following expression q by 21T , we get: q is a measure of spatial refinement. 41 h .L q* 2TI (3.2.46) ;\ which represents the number of elements per wave length. to limit study to the range q E Jo, TI[ For example It is reasonable = TI/2 q translates to having 4 elements for one full wave form. a = ~ For the entire graphical analysis 3.3 Comparison of Algorithms 3.3.1 Implicit Algorithm Fig. 3.1 shows the frequency ratio and PG(Pade) algorithms. accurate. term GC CPL condition. (w/w) for finite GL which vanishes for For finite for the implicit q , GC decreases; the opposite is true for algorithms. PG(c T = c t) 2 6 and like PG(C2T = 1) and PG(Pade) GL become more accurate as behaves like GC as as satisfies the unit c t 6 PG(Pade). , PG(C2T c t 6 -+ 6t -+ 2//15 . (~/w) = 2//15) They are .all second-order accurate. PG(C2T = 2/115) GL· This is due to the q GL. Fig. 3.2 shows the algorithmic damping ratio ratio for the implicit GC, They have no modulus error, and are second-order is more accurate than 2 r(l - cos q) (trapezoidal) . O , and and the frequency PG(C2T = c6t) As one might expect, like PG(C T = 1) 2 as c6 t -+ 1 42 GC IMPLICIT fAJ Q N c.0.t 1.0-- .a-·- "' .6--.4 ------ .2 ......... 0 I- -o a: - • 'l:: >-- e.J zw. II> ::::lO 0 w a:: r.... Q ci o.o Q o.s 1.0 1.s 2.5 2.0 3.0 DIMENSIONLESS WAVE NUMBER GL IMPLICIT fBl "4 0 -o ~..,;,_ ___ 0::: >- C.J cs"! E ::::l Cl 0 Cl ci-+-----......~----------~--......~--.....~--~ o.o o.s 1.0 1.5 2.0 2.5 3.0 DIMENSIONLESS WAVE NUMBER Cl PGCPADEJ IMPLICIT CCJ N "' ...;- 0 -o I-..: a: a:: >-- c..J z.,, ~ c;;0 w a:: t... Cl Q 0.0 0.5 1.0 1.5 2.0 2.5 3.0 DIMENSIONLESS WAVE NUMBER Figure 3.1 Frequency ratios for implicit Galerkin and Petrov-Galerkin (Pade) algorithms 43 PGCC2T•ll IMPLICIT ~ CAI d I .: //j' // I/ , ,,,-· /.,.:.·~· /.-~/ ,#;.efJP' c~, 0 ~ lQ 0:: • Q C!I z ~~. 8 ,,,,.,..-· "' 0 - c• i- cc- 1.0- 0:: .&--- t.J .8-·- ~ !: ig a::: • >< u z.,.• .4 - .2 ......... o'i' t!) ::>a 8 e: i~ c? o.o o.s 1.0 1.5 2.0 Oit'IENSJONl..ESS MFIVE NUMBER 2.s . d J.O PGIC2T•2/SORT!l5ll IMPLICIT ~ CCl PG CC2T• 1l . IMPLICIT (8) Q .... Q a.a o.s 1.0 1.s 2.0 DIMENSIONLESS WAVE NUMBER rDJ 3.o 2.s PGCC2T-2/SORT(15ll IMPLICIT .... "' _. 0 i- cc - • 0:: >- ~ ... t..l • ::::> Cl 0 t..l 0:: to.. .. ci o.s 1.0 1.5 2.0 D111ENSIONl..ESS WAVE NUMBER 2.5 o.o J.O o.s 1.0 1.5 2.0 DIMENSIONLESS WAVE NUMBER PGCC2T•COTJ IMPLICIT PG!C2T•CO!l IMPLICIT [f") - 3.0 2.5 Q .,, 0 ./ / I- - Cl• --~ ·· .. cc0:: >- ~ ... w. ::>a w a::: 0 to.. Q c:i 0.5 1.0 1.5 2.0 DIMENSIONLESS WAVE NUMBER Figure 3.2 2.5 3.0 o.o 0.5 1.0 l.S 2.0 DIMENSIONLESS WAVE NUMBER Algorithmic damping and frequency ratios for Petrov-Galerkin algorithms 2.5 3.0 44 3.3.2 Explicit 1-Pass Algorithms Fig. 3.3 shows the results for the explicit 1-pass and GL are uncondition ally unstable; not shown.) PG algorithms (GC PG(Pade) and PG(C 2 T are equivalent because, for 1-pass algorithms, there is no dependence on In this group, only PG(C T = CDt) 2 is second-orde r accurate. 1) r This superior- ity with respect to order-of-ac curacy can easily be seen from (3.2.29), which can be written as: l -zl2 = 1 - CDt V(2(C2T - CDt) Clearly, the departure of However, if c 2 T = CDt , 1 -zl 2 2 + CDt V(l - C2T)) from unity is first-order in ( 3. 3 .1) CDt and then = and thus, the departure of and V . lz-12 (3.3.2) from unity is now second-orde r in cut V . The algorithms PG(Pade), PG(C 2T = 1) and PG(C 2 T : The stability limits are c6t < 2/115 for PG(C 2 T = 2/liS) CDt) satisfy the unit CFL condition. for the other PG algorithms. and CDt < 1 45 i PG!PRDE,C2T•ll El CRJ Cl PGCPAOE:,C2T•l l E:l CBI ,,; C.c.1 1.0 .8 .6 .4 --·- ----.2 ·········· "' 0 ,/ ,::;::_. -c:i """ cc -. ·-·--......:::.:::.=::_:--·····-·-···········... a: ,.. ~ u ~i;& i.. •. ··..... a.a a.s 1.0 1.s ~a 2.5 3.0 .. a.a DIMENSIONL.ESS WAVE NUMBER PGCC2T·2/SORT!l5ll El !Cl ... Cl ~IQ \ "· Cl t?-1---~~--......----......--...,.----~----t i ', ··... z.,, w. ::i Cl 0 w a: !: !I 4? !5 C!I "· .. ....... 2.5 o.s 2.a 1.5 1.0 DIMENSIONLESS WAVE NUMBER (OJ l.O PGCC2T·2/SORTl15ll El Q "' a: Cl. ~ - 0 """ -. cc a: -c:i ,.. ~- u ·--~~:~ O"! ::i Cl 0 w e: o.5 1.0 1.s 2.s l.o 3.D . ·----···:;::;: ··....~.~.::~::!': Cl a.a . i PG!C2T•CDTl E:l !El a.s 1.a 1.s 2.a OIMENSIONLESS WAVE: NUMBER DIMENSIONLESS WAVE NUMBER (f") Cl 2.s 3.a PG!C2T•COTJ E:l ,,; "' ..... a: 0 - Cl --~- a: ~~;?:::;;;::;-.;~~ ,..u z.,, w. ::i Cl 0 w a: i.. . ------- // Cl a.s 1.0 1.5 2.a DIMENSIONLESS WAVE NUMBER Figure 3.3 2.5 3.0 a.a 2.s 2.0 1.5 1.0 o.s OIMENSIONLE:SS WAVE NUMBER Algorithmic damping and frequency ratios for explicit 1-pass algorithms l.O 46 Explicit 2-Pass Algorithms 3.3.3 Fig. 3.4 shows the algorithmic damping ratio and the frequency ratio for explicit 2-pass not shown.) PG(Pade) for algorithms. (GC and c~t = GL are unconditionally unstable; In this group, all the algorithms are second-order accurate. satisfies the unit CFL condition. C~t .::_ l around PG , except for 0.8. PG(C 2 T = 2/115) All the algorithms are stable which exceeds the stability limit 47 PGIPAOEl E2 ~ fAl PGfPAOEl E2 (91 C6t 0 ..... 1.0-- ~ .4-- .a-·.6--- ~~0 .2 .......... . a.. 5§ 2 00 u -- --- .. :..-·:;::::;~ ~ "!:: lS .--· ------ ~c;i C.ll ~!<i 0 o+-...........,..............,._..........,.............,..............,.............~ o.o 0.5 1.0 1.S 2.0 2.5 3.0 DIMENSIONLESS WAVE NUMBER c;i-1-~~--~.....~~...---..--.....-~~ o.o a,5 1.0 1.s 2.0 2.5 3.a DIMENSIONLESS WAVE NUMBER PGCC2T•IJ E2 i CCJ PGCC2T·l l E2 0 -........ · u ~ = lS ~c;i C.ll iE !<i "' 0-1-~~...............,,.....................................................__. ~-i-~~....-~......,....~.....~~-i--~~~--! o.o o.5 1.0 2.0 1.s 2.s DIMENSIONLESS WAVE NUMBER 3.o a.a o.s 1.0 1.5 2.0 3.o 2.5 DIMENSIONLESS WAVE NUMBER PGCC2T•2/SORTCl5ll E2 PG!C2T•2/SORT!1Sll E2 (fl C> N.,..~.,.....,.....,.....,.....,.....,........,..,.......,..,.....,.......,..,.......,.....,. 0 / . I /II / . I /I I / .I ....../ ,// ; / 0 _, ..><-:::::-./ / o.o --- .\ -.&E. . "-·1-------~Qiii~~~ ........ ...... _// ,.... ·· ...~ u ~ ~~ ..... ~ "' 0-+-~~ .....~-...................~..................,_.....~ a.s 1.0 1.s a.o 2.s 2.0 DIMENSIONLESS WAVE NUMBER 0.5 l.O l.5 2.0 2.5 3.0 OIMENS!ONLESS WAVE NUMBER PGCC2T-COTl E2 PGCC2T•CO!l E2 ~ •. · ..... 0 0-1---~......~--..--~--~--.......--~........--~ 0.0 0,5 1.0 1.5 2.0 DIMENSIONLESS WAVE NUMBER Figure 3.4 2.5 3.0 o.o a.s t.a 1.s 2.0 DIMENSIONLESS WAVE NUMBER 2.s Algorithmic damping and frequency ratios for explicit 2-pass algorithms 3.0 48 CHAPTER 4 Numerical Applications in One Dimension 4.1 Introduction To test the capabilities of our finite element schemes and to estimate the effects of several algorithmic parameters involved, we experimented with problems from various classes of one-dimensional hyperbolic systems. We started with a linear transient problem which had a discontinuous solution. Then we studied nonlinear transient problems with continuous and discontinuous solutions; these illustrated the concepts of shock stability and admissibility. Further, we experimented with a set of nonlinear steady problems; these resulted in continuous or discontinuous solutions depending on the way the boundary conditions were specified. 4.2 4.2.1 Numerical Applications in the Linear Transient Case Propagation of a Small Disturbance in a Gas Barotropic Compressible Flow Equations The Euler equations in one dimension are given in appendix I. The mass and momentum conservation equations can be uncoupled from the energy conservation equation by assuming that the flow is barotropic, that is: p = p(p) (pressure) (4.2.1) Then, the barotropic flow equations can be written as a system of conservation equations with conservation variables and flux vector defined as: 49 (4.2.2) (4.2.3) where p and u are density and velocity, respectively. Small Disturbance Equations We assume that the departures of and uo are very small. stated as (see Taking 1 p-' (p) p .. (Po) - l u p .. (p) and u from constant values p 0 to be zero, this assumption can also be [w4]): £ Po where p << 1 (4.2.4) << 1 (4.2.5) << l (4.2.6) is the derivative of p(p) with respect to its argument. With these assumptions, we get the following linearized version of the original conservation equation system: P,t Po u,t + (4.2.7) 0 = 0 (4.2.8) 50 Further, by introducing the parameter 2 c c ' P ... <Po> (4.2.9) PIP 0 - 1 (4.2.10) - u/c (4.2.11) = and the scaling P* = u* :::;: + A we get u~,t u -,x = 0 (4.2.12) where u [:: 1 : :;: (4.2.13) and A Clearly, the eigenvalues of A 0 :::;: [ -c -c] (4.2.14) 0 are: (4.2.15) 4.2.2 Initial/Boundary-value Problem The equation system of (4.2.12), together with the following initial/ boundary-value data was studied numerically by Hughes in [H6]: p (x, O) u (x, O) = = 0 1 l x E ] 0, 10[ (4.2.16) 51 p(O, t) = u(lO, t) = t > 0 ( 4. 2 .17) (We drop the asterisks for notational simplicity.) The parameter c was chosen to be unity. 4.2.3 Finite Element Solutions Naming Convention for the T=T b If we choose to define the operator T =T ~ , and T T = T ~T according to the criterion then the resulting formulation would have a second-order term containing the product ATA . Therefore, we name this criterion the form." If we choose the criterion T = T AT second-order term would contain the product criterion the 2 "A -form". "ATA- on the other hand, then the 2 A . Therefore we name this Thus, in the study of this and all the othe!:' problems, we adopt the fo].lowing naming convention: T = T = T A - (4.2.18) - (4.2.19) This naming convention will also be used for the multi-dimensional cases. Algorithmic Features In this problem, both ATA and A2-forms result in the same formula- ·tion due to the symmetry of the opera tor A We can also set T = 0 and obtain the usual Galerkin technique. The finite element mesh contains 20 elements with uniform mesh spacing 52 of .5 . The quadrature rules chosen provide exact integration of the vectors and matrices involved. The transient algorithm parameters respectively. Y and a were set to be 1.0 and 0.5 Implicit, explicit 1-pass (designated by El), and explicit 1-pass (designated by E2) algorithms were tested. Two different time steps, 0.50 and 0.25, were used; these time steps correspond to Courant numbers 1.0 and 0.5, respectively. For both time steps, we used the temporal criterion for given by (2.5.8). for T , The parameter F was set to unity. as T The spatial criterion (2.5.1), would result in the same formulation for both time steps, pro- vided that we set F = 0.5 when ~t = 0.25 Results Fig. 4.1 shows the implicit Galerkin solution for Courant number 1.0 . As can be seen, this technique produces spurious oscillations. Fig. 4.2 shows the solutions produced by the explicit 1-pass Petrov-Galerkin algorithm with a = 1.0 and a = 0.5 . One would expect the solutions to be very similar because the explicit 1-pass algorithrn is independent of Slight differen~es occur a because the initial start-up conditions are handled differently by each algorithm. Both algorithms satisfy the unit CFL condition defined in chapter 3 in this problem, a = 0.5 up conditions. produced results inferior to a = 1.0 yet, due to start- Therefore, for the explicit 1-pass algorithm, when it comes to a choice between a= 1.0 and a= 0.5, the fo:nner value is preferred. Fig. 4.3 shows the solutions (for Courant number 1.0) produced by Petrov- EXACT - - GALERKIN ............. ::J N 10.0 15.0 s.o I Figure 4 .. 1 T 10.0 20.0 N I 0.0 5.0 s.o lI l T 10.0 T 10.0 15.0 15.0 Small disturbance propagation in a gas: Galerkin implicit, ne\l, = 20, 6t = . 5, Courant number = 1.0 . 15.0 0 ....z I 0 ::::> - r-d a: I Xe i 0 - I I 0 i ....... (0) 0 <'i ' 0.0 0 I T ...i I o.o .I · .. · 20.0 0 r-d a: Xo I LO 0 r-l s.o C> r-.iI 0 C> I ::::> f-oO a: ~ ~J - ~ c-,1 (8) 0 I 0 (Cl C> ::J N a: 0.0 1 !q~ f-oO I 0 :1 (A) 20. - 20.0 Ul w 0 '. ::::> (\J 5.0 ' EXACT - - - 0 0.0 5.0 I Figure 4.2 15.0 20.0 r-lI o.o 5.0 5.0 I'\ 10.0 ... .. '. \\ r .. . T \ T 10.0 I ''" explicit 15.0 15.0 \ -------. Small disturbance propagation in a gas: global T 1-pass, ne!/, = 20, L1t = . 5, Courant number= 1.0 . T 10.0 I : ; ~J f- 0 a: o.o -I Xo I - 0 0 (0) 0 r-3 0 I I 20.0 r-lI ., 15.0 rJ, - , - 0 ' - o I I I .Y : T 10.0 a• 1.0 ......... a =0.5 -------- .... ::> I 0 I- 0 a: o.o :1 o r-iI 0 ' C! I Xo Vl (C) ::::> (\J ~ ~~ a: . ..:.- , 0 r-o d r-l (B) o i-d a: J: .. I .: _ _ __ If) C! <"! 9X o ·:/ (A) 20.0 20.0 U1 .!» 55 Galerkin implicit, explicit I-pass and, explicit 2-pass algorithms. The explicit 2-pass results are not distinguishable from the implicit results. This implies that, for this problem, as the number of passes increase, the explicit algorithm converges quite rapidly to the implicit one. Fig. 4.4 shows the results for Courant number 0.5 produced by the same set of algorithms. The implicit and explicit 2-pass algorithms are, again, .indistinguishable. The explicit 1-pass algorithm produces slightly greater oscillations. ::::::> C'IJ a: .. 20.0 ::> a:: I 0 f--0 Xo I - 0 ~ N (D) o ' 0.0 Figure 4.3 s.o • \\......... \ 5.0 \ .....- I : \ .\ \\ "' T 10.0 '~ '. T 10.0 \ ../ - 15.0 15.0 I \ \\ "i\ -,;_ ,_./ Small disturbance propagation in a gas: global T , ne~ = 20, 6t = .5, Courant number = 1.0 (Il and E2 are virtually indistinguishabl e) . T 0.0 ' 0.0 .-S 0 I 0 NI 15.0 - 20.0 ::> a:: 0 10.0 • 15.0 • Xe f--0 lJ) N 5.0 T 10.0 1 1: 0 0 f--0 Xo lO I ~ r-l : <'iI r---- - - - - - , 0.0 5.0 0 E2 - - - ' I EXACT ..;- 0 Ii l1 El (C) o ::::::> C'IJ cc f--0 Xo ~ '\ ~ ~t7'' :; ,..,j (8) 0 ,..,j (A) o _ _ _ _ _ __:__ _ _ _ __ 2C.0 20.0 - Ul ()'\ s \..'Z...... '" K" '.::) a: Xe 0 0 Figure 4.4 .....; ····' \\,.,-_, I\\ •• Small disturbance propagation in a gas: global T , ne~ = 20, 6t = 2.5, Courant number= 0.5 (Il and E2 are virturally indistinguishabl e) . T ''\ ·· .• , T ·--·· -~ Li ., 20.0 ~-+------.------~ ----------' 0.0 5.0 10.0 15.0 20.0 T 0 20.0 • N 0 ~-+----------.----------.---------. . . . ---------4 \'. ,~ ~---------------------~ --------~ o.o 5.0 10.0 15.0 0 I 0 i-a· I 0 (D) 5 I 15.0 ·--:-:\ 20.0 ----·:·\ X~\. • ...... 0 '--""'<c,__.,.,., ....:" ..........-----------~::: a: .... m I 0 lO.O T 15.0 --- E2 10.0 -------- EI \ c I 5.0 .. -.~' ) II:' ...r'~. 7· ,: 5.0 == EXA~~ ---. N (8) 0 0 0.0 ~I 0 N 0 ~ 0.0 0 I 0 ~ 0 x I LJ) (Cl g a: Xe t-- 0 oI ·-- (A)o~ ~~~~~~~~~~~N -..J Ul 58 4.3 4.3.l Numerical Applications in the Nonlinear Transient Case Barotropic Compressible Flow Barotropic compressible flow was defined in section 4.2. The differen- tial equations need to be satisfied everywhere except a.t the s!1ock front where the Rankine-Hugoniot conditions (see [Ll, L2, L3]), = have to be satisfied. and [ Here s (4.3.1) is the propagation speed of the shock ] is the jump operator. That is, for any variable fro~t Q : (4.3.2) where the superscripts "-" "+" and the shock front, respectively. refer to the left and the right of For barotropic flow these conditions are: s[p] s[pu] = = [p u] (4.3.3) [Pu 2 + p] (4.3.4) For shock profile to be stable, the entropy condition also must be satisfied (see [Ll, L2, L3]). inequalities in terms of s The entropy condition is given in the form of and the eigenvalues of the jacobian matrix A . In the present case, the eigenvalues are = (4.3.5) 59 4.3.2 Initial-value Problems We considered two initial-value problems, both studied by Hughes in [H6], with the following equation of state: 1 27 = p(p) (4.3.6} The first problem has the following set of initial data: where H(x} p{x, 0) = 1 + 2H(- x) (4.3.7) u(x, O} = f H(- x) (4.3.8) is the Heaviside step function. satisfy the jump relations. This initial data does not Therefore, the initial shock profile splits into a stable shock which propagates to the right and a simple wave which propagates to the left. The second problem has the following set of initial data: p(x, O) = 1 + 2H(+ x) p {x, 0) = 3 H(+ x) 2 (4.3.9) (4.3.10) This is the mirror image of the previous data with respect to the point x =0 . This initial data does not satisfy the entropy condition and there- fore represents an unstable shock. ing to the right. The result is a rarefaction wave travel- 60 4.3.3 Finite Element Solutions Algorithmic Features 2 A -forms were tested on these problems. and . We also tried the Galerkin algorithm for one case only. The finite element mesh contains 40 elements with a uniform element length of 1.0. The transient algorithm parameters respectively. y and a were set to 1.0 and 0.5, The time steps were taken to be 0.6 and 0.3 corresponding to Courant numbers (based on the maximum eigenvalue) 1.0 and 0.5, respectively. The parameter by (2.5.8). T was chosen according to the temporal criterion given The parameter F was usually taken to be one, however we tested cases where it was greater than one. Results Fig. 4.5 shows how the Galerkin algorithm performed for We used an implicit 3-iteration scheme. ~t = 0.6. The location of the shock (that is the shock speed) is in agreement with the exact solution; but there are spurious oscillations behind the shock. The Petrov-Galerkin algorithms, in general, performed quite well. The conunon discrepancies between the numerical and exact solutions were, with varying magnitudes, overshoots at the shock fronts and oscillations behind the simple wave. We tested implicit schemes with I, 2 and 3 iterations (designated by II, I2 and I3) and explicit schemes with 1, 2 and 3 passes (designated by El, E2 and E3). We found that at least 3 iterations were needed to get the correct shock structu~es when an implicit scheme is used. This observation is in agreement with the findings of Baker [B2, B3]. 61 (A) 0 ..,...,..----------------------------------------------------0 . •' C4 ' "' 0 "'':i .- ~ a a::: 0 ·EXACT GALERKIN 0 . a.a ' 10.0 I I 20.0 30.0 x 40.0 ( 8) 0 ~...,...------------------------------------------- 0 ...:- - ..... . ... .. ·.. ,-· 0 0 ,. .. ·.. ,.,. ... ...... · ....... "· \! l ·- 1· 0 T-+--------------,....,------------...-,--------------,----------.....J 0.0 Figure 4.5 10.0 20.0 x 30.0 40.0 Barotropic compressible flow, stable shock propagation: Galerkin implicit 3-pass, net = 40, 6t = .6, nts = 25, Courant number (based on maximum spectral radius) ~ 1.0 . 62 The element level "mass" matrices and the vectors corresponding to them were integrated exactly. For the integration of all the other matrices we tested 1, 2 and 3 point Gaussian quadrature rules. Fig. 4.6 shows the com- parison of the integration rules for ~t = 0.6 2 A -form temporally. with and chosen T We conducted the comparison tests on implicit 3-iteration (I3) and explicit 1-pass (El) algorithms. We observe that the results change One can therefore use 1-point qua~ only slightly with the quadrature rule. drature for economy reasons without much decrease in accuracy. However, in the following problems we used the 3-point quadrature rule. Fig. 4.7 shows the results for with the exact solution. ~t = 0.6; all are in close agreement All have, with comparable magnitudes, overshoots at the shock front and oscillations behind the simple wave. explicit 2 and 3-pass results are very similar. For the We observe that ATA-forrn the explicit 1-pass algorithm became unstable; for this form we also tested the implicit 3-iteration scheme with F = 2 . This slightly reduced the oscil- lations and the magnitude of the overshoot to the left of the shock front. Fig. 4.8 shows the results for ~t = 0.3. The solutions are in agree- ment with the exact solution with slightly more oscillations behind the shock front compared to the ~t = 0.6 case. For the A2 -form, there was virtually no difference between explicit 2 and 3-pass algoritluns. We also tested the ~ explicit 2-pass algoritlun with F = 2 . This algorithm reduced the oscilla- tions and the magnitude of the overshoot to the left of the shock front. Fig. 4.9 shows the results for the unstable shock problem. the 2 A -form with the temporal choice of T We used the time step was 0. 6. We tested the implicit 3-iteration and explicit 2-pass algoritluns. Both solu- tions were in close agreement with the exact solution. algorithm The El 0 a::: I 10.0 Figure 4.6 30.0 10.0 ::J ? ~ d (D) IS 0.0 l- I -,--10.0 I - -- x 20.0 30.0 . ·1 I I JI 11 ft ¥ I \ ···Ll 30.0 I ~.. ·~ I - x I 20.0 ,t ........ oe::.7'c·_ ../ 10.0 I '· " t\""., \ • ;:·::::;;;-;c.,.;-- ;;;.;;~~ .,.-- ..... ~~~'!; I I El 0.0 13 iO.O 10.0 Barotro pic compres sible flow, stable shock propaga tion: global T , A2 -form, ne£ = 40, 6t = .6, nts = 25, Courant number (based on maximum spectra l radius) = 1.0. Compari son of differe nt element quadrat ure rules. x 20.0 I a.o ~I I 11 I I ,, ..................... . c- ... 10.0 I ~1 30.0 I -v~ I --- ~r x I 20.0 I f\\ ::J I I I I c:i CD (B) in \\ ~! I I I I I 10.0 2- PT ---------3-PT - - - I-PT ............ EXACT \I '''.) ?\i n 13 El dia.a (C) "J ~ 0 n (A) "I w (j) 0::: Cl (C) Cl 0::: (A) (31 .. (4) f • I ,E3 x . ... ' ~- \' ... " ... =i (0) :::::> (B) 12) (I) ,., (3) ,., (3)1_2)_ _ (I) x -L ' '·' t~ L~ ~- -L~l~''~· v •"l....,...,,......-t -~-~:. GLOBAL ATA x -- ---- GLOBAL A2 , F=I Barotrop ic compress ible flow, stable shock propagati on: nei = 40, 6t = .6, nts = 25, Courant number (based on maximum spectral radius) = 1.0. Compariso n of various SU/PG schemes. . ' .. !31 F•l,E2 /·~ ·"\ . f • I , 13 .. .. ... - .. 12) F•2,I3 (I) (2) '.ll GLOBAL ATA .... x .' .··\1 I I\ J. .- D:. -· -- - (4} E3 .. '-..·-- 131 E2 121 EI (II 13 --...., .. Figure 4.7 (4) I (4) (3) (2) (I) GLOBAL A2 , F =I ""' (j'\ a:: 0 (C) a:: 0 (A) (3) .. (4) 12 l 13 l . :.~ x I " F l: ' .' f ·' 1· ~ ·. I . .. I\ ":· . .1 t: ~. ~ ~.. ' ' ~I ,l ~• => (0) => (ffi C4l 131 I (41 I • x #" x , I' .. ..t lIL'.: \'v .·-..l :: ! 'l 1\\. -~ • GLOBAL A1 A, F =I (l_l__ (3) - - (2) Ill GLOBAL A2 Barotrop ic compress ible flow, stable shockypro pagation: ne£ = 40, 6t = .3, nts = 50, Courant number (based on maximum spectral radius) = 0.5. Compariso n of various SU/PG schemes. .. ··-· ........... ..... ·. E! E 2 -E3 (I) I 3 12) ( 11 x .. ( . .. ,.:·~ GLOBAL ATA, F = I (41 F =2.E2 ( 11 F = I, 13 (2) F = l,E I l3lF=l,E2 ........... .. ---...... ... .•, ~ Figure 4.8 (4) I I (4) (31 (2) !II GLOBAL A2 O'I Ul 66 (A) 0 0 :.,, 0 N 0 ........ 0 0 EXACT . !3 ··········· E1 ----- 0 7'4------------------------------------------------------1 40.0 10.0 30.0 0.0 20.0 x (B) 0 ~------------------------------------------------------- 0 N 0 0 Cl-+-------.__,,,.._-0 . ______________________ ______ ';""_,_ _.. 0.0 Figure 4.9 10.0 20.0 x ~~~~--~~--~~~-' 30.0 40.C 2 Barotropic compressible flow, rarefaction wave: global T A -form, neR. = 40, 6t = . 6, nts = 25. Courant number (based on maximum spectral radius) = 1.0 . . 67 produced slightly more oscillations behind the simple wave. All algorithms produced good results for this case. 4.4 4.4.l Numerical Applications in the Nonlinear Steady Case Isothermal Flow in a Nozzle We consider the one-dimensional isothermal flow \in a nozzle with cross- sectional area varying along the axis. The governing balance law equations, provided by Lomax et al. [L7], possesses the following conservation variable, flux and source vectors: = (4.4.1) A (4.4.2) = (4.4.3) u G where the acoustic speed c 2 is constant and the cross-sectional area is A(x) A (4.4.4) The jacobian matrices are: 0 A a.-wau 1 = (4.4.5) -u2 + c2 2u 68 0 0 a~;a~ = (4.~.6) 2 -c A, /A x The eigenvalues of A 0 are ~ u + c ,\1,2 Assuming that A(x) is a continuous function of (4.4.7) x , the Rankine- Hugoniot conditions for steady flow reduce to: [p u] = [pu 2 O 2 + pc ] (4.4.8) = 0 (4. 4. 9) 69 4.4.2 Boundary~value Problem We studied steady flows suggested by Lomax et al. [L7]. The cross- sectional area and the acoustical speed were given by A(x) 1. c 0 + (x - 2 • 5) 12.5 2 0 < x = 1.0 < 5. (4.4.16) (4.4.17) The problems considered were: 1. Subsonic inflow - subsonic outflow with no shock. 2. Subsonic inflow - supersonic outflow with no shock. 3. Subsonic inflow - subsonic outflow with shock. The exact solutions, which can be obtained by the integration of the square of the Mach number (in this case u 2 ), were provided by Lomax et al. [L7]. 70 4.4.3 Finite Element Solutions For the boundary conditions of these problems, we set the values of the conservation variables as given by the exact solution. The number of variables to be specified at each boundary depends on the nature of the flow at that boundary. For supersonic inflow, two variables are set; for sub- sonic inflow or outflow, one variable; and for supersonic outflow, no variable is specified. Considering all the combinations of possible boundary conditions in terms of conserv~cion variables, we have the following cases for each problem: For subsonic inflow - subsonic outflow problems: UlUl: ul at the inflow/u1 at the outflow UlU2: ul at the inflow/u 2 at the outflow U2Ul: u2 at the inflow/u at the outflow 1 U2U2: u2 at the inflow/u 2 at the outflow For the subsonic inflow - supersonic outflow problem: Ul: u U2: u2 1 at the inflow at the outflow Transient Introduction of the Source Term In these problems, we introduced the source term into the equation system in a transient fashion. the term ally. A,x That is, instead of having the full value of right at the beginning, we let it reach its full value gradu- This is done by taking A ,x as a linear function of time during 71 an initial time interval at the end of which The numerical A,x where nti interval. reaches its full value. can be expressed as follows: (A, A A ,x x) NUMERICAL ,x = (4.4.18) denotes the number of time steps marking the end of the transition For the problems solved, nti = 10. Algorithmic Features Both ATA 2 A -forms and were employed. We also tested the Galerkin algorithm. The finite element mesh has 40 elements with uniform element length of 0.125. The element level "mass" matrices were integrated exactly; all the other matrices and vectors were integrated by the 3-point Gaussian quadrature rule. We set the transient-algorithm parameters Y and a to unity and employed implicit schemes with 2 iterations. The parameter T was chosen according to both spatial and temporal criteria given by (2.5.1) and (2.5.8), respectively. The time step for each problem was usually chosen to be ten times the estimated critical time step for that problem. step ~tCR We define the critical time as the time step for which the Courant number, based on the maximum spectral radius, is unity. That is = h/max p (~) (4.4.19) 72 p(~) where is the spectral radius of A. While full convergence to the steady state solutions was attained in about 100 steps, the 50-step solutions were close enough to the steady state solutions for practical purposes. Results for the Subsonic Inflow-Subsonic Outflow Problem With No Shock For this problem the time step was chosen to be 0.735. We attempted to solve the problem with all possible boundary-condition types. Of the four types tried, only one, UlU2, failed to give the expected solution. For each boundary condition type, we tested the usual Galerkin algorithm, and 2 A -forms, the latter-two with temporal choice of For the type UlUl we also tested the In all cases F = l 2 A -form with spatial choice of T • T . For each boundary-condition type, there was no difference between the solutions produced by different algorithms. However, the solutions differed slightly from one boundary condition type to another. agreement with the exact solution was very close. For all types, the Fig. 4.10 shows the results. Results for the Subsonic Inflow-Supersonic Outflow Problem With No Shock The time step was chosen to be 0.460. We solved the problem with both boundary-condition types Ul and U2. each type, we tested the Galerkin algorithm and latter two with temporal choice of T . The results are shown in Fig. 4.11. T A A , and In all cases F = 1 For 2 A -forms, the For each boundary-condition type, there was no difference between the solutions produced by different algorithms. However, the solutions differed slightly from one boundary-condition type to another. For all types, the numerical solutions were in close agreement with 73 (A) 0 N,....------------------------------------------. Lf) -~- a~ a::: - EXACT UIUI U2UI . I.I'> 0 U2U2 .......... --- --- --- 0 0 o.o 1.0 2.0 x 3.0 4.0 .S.O ( 8) 0 N..,...-----------------------------------------. t.f') 0 . . I./) 0 0 0-f----------.,----------..,....--------.....--~--~------~~~~ 0.0 Figure 4.10 1.0 2.0 x 3.0 4.0 5.0 Steady nozzle flow, subsonic inflow-subsonic outflow, with no shock: nei = 40, 6t = 0.735. Comparison of different boundary conditions. (All methods give essentially the same results.) 74 (A) 0 N-r-------------------------------------------- ao 0:: EXACT . Lf) 0 UI ·············· U2 ------- 0 o-r---------r---------,.-------~r------------------...J a.a 1.0 2.0 x 3.0 4.a 5.0 ( 8) 0 N-r--------------------------------------. If) 0 l./') 0 0 . o-;---------...,.----------r----~----.--------~--------~~ a.a 1.0 2.0 3.0 4.0 5.0 x gure 4.11 Steady nozzle flow, subsonic inflow-supersonic outflow, with no shock: nei = 40, 6t = 0.460. Comparison of different boundary conditions. (All methods give essentially the same results.) 75 the exact solution. Results for the Subsonic Inflow-Subsonic Outflow Problem With Shock Unless specified otherwise, the time step was taken to be 0.5. We also attempted to solve this problem with all possible boundary condition types. Types UlUl and U2Ul gave the expected solutions. We first describe the solutions obtained for type UlUl: Fig. 4.12 shows the results for the choice of and 2 A - forms with the temporal The solutions are in close agreement with the exact solution T . everywhere except at the shock front where the shock front is not very crisp and shifted to the left by half an element length. Fig. 4.13 shows the results for the F The parameter assumes values 1, 2, 5 and 10. close agreement with the exact solution. near the shock front for low across one element only. element length. F . For with spatial choice of T The solutions are in very There are very slight oscillations F = 1 and F = 2 the shock front is The error in the shock location is about half an As F. increases, the shock front becomes smeared. Fig. 4.14 shows the results for the T . 2 A -form T A A-form The results are similar to that of Fig. 4.13. a. The shock fronts are slightly less crisp. b. For with spatial choice of The only differences are: ! we observe oscillations behind the shock front. F = 1 , It is interesting to note that the oscillatjons are located in the region between the shock front and the point where the flow velocity reaches the speed of sound. UlU2: 2 A -forms, solutions with no shock. with temporal choice of T , produced smooth 76 (A) (2) ( 1) (2) x ( 8) (2) x ·ure 4.12 Steady nozzle flow, subsonic inflow-subsonic outflow, with shock: boundary conditions UlUl, global T, ne£ = 40, 6t = 0.500. 77 (A ) (4 ) (1) 0 0::: I I I ' ( 1) F = 1 (2) F= 2 ( 3) F = 5 (4 ) F = 10 x ( 8) 4) ~igure 4.1 3 x Ste ad y no zz le flo w, su bs on ic inf low -su bs on ic ou tfl ow , wi th sh bo un da ry co nd iti on s oc k: Ul Ul , lo ca l T , A2 , ne z = 40 , ~t = 0.5 00 . 78 (A ) (4) , ,, ( 1) I 0 I a::: (1) F = 1 (2 ) F = 2 (3 ) F = 5 ( 4) F = 10 x ( 8} (4) Lgu re 4.1 4 x Ste ady no zzl e flo w, sub son ic inf low -su bso nic ou tflo w, wi th sho ck: bou nda ry co nd itio ns UlU l, loc al T , ATA, ne i = 40 , 6t = 0.5 00 . 79 U2Ul: The time step was taken to be 10, 20, and 40 times the estimated critical time step. The results for 2 A forms and temporal choice of are very similar and are shown in Figs. 4.15 and 4.16 respectively. T The shock fronts are shifted to the right about one element length. U2U2: Both ric ATA solutions. and 2 A with terrq:>oral choice of T , produced smooth syrnmet- For this case, it was only the Galerkin algorithm which sensed the shock and located it almost at the exact location, but with severe oscillations. Fig. 4.17 shows t~e resnlt produced by the GaJerkin algorithm. Remark We observed that the location of the shock front was shifted about half an element length to the left for the boundary condition type UlUl and about one element length to the right for the boundary condition type U2Ul. This implies that the location of the shock front is dependent to some extent on the type of boundary conditions specified. In particular, for two boundary condition types (UlU2 and U2U2) the exact solution was not obtained. Proper specification of the boundary conditions in problems of this type is a very important subject which has attracted several researchers {see [B6, Ml, Yl]}, but does not yet seem to be fully understood. 80 (A) (1) ,----- , ,' ,' '' ' ,, ' '' ,'' ,, '' '' 0 a::: ( 1) At = 10 At CR (2 ) At = 20 AtCR (3 } Lit = 40 AtCR x ( 8) .'. '' ' I I ( 1) (3) Pi gu re 4.1 5 x St ea dy no zz le flo w, su bs on ic in flo w- su bs on ic ou tfl ow , wi th bo un da ry co nd iti on s sh oc k: U2 Ul , gl ob al T , A~.., , ne i = 40 . 81 (A ) ( 1) ,,, .. "'" .. ,, 0 a::: ,, I , , . I I I I I (1 > at = 10 atc R , I I I (2) At = 20 AtCR I I (3 ) At = 40~tCR x ( 8) ( 1) (3) Fig ure 4.1 6 x Ste ady no zzl e flo w, sub son ic inf low -su bso nic ou tflo w, wit h sho ck: bou nda ry con dit ion s U2U l, glo bal T , ATA, nei = 40. 82 (A) 0 c'1------------------------------------------------------------..... :'\. \..•.../ ....... \_ ..... .....·... __ ·-...... ·"\.. ·......... ·... ..... ·. ···•··.... "····· .. ···.........\ ... !-. ,/\, v ·............. 00 er:: /\ .... .......,: . I.I) 0 EXACT GALERKIN ··············· . 0 o...,______...___________________________.,..__________,....________.... o.o l. 0 2.0 x 3.0 4.0 5.0 ( B) 0 c'1111------------------------------------------------------------- I./) 0 => - . . LI') 0 . ·O o...,-----------.------------------------------------------------1 o.o 1.0 4.0 2.0 3.0 5.0 Figure 4.17 x Steady nozzle flow, subsonic inflow-subsonic outflow, with shock: boundary conditions U2U2, ne£ = 40, 6t = 0.500. 83 CHAP TER 5 Mul ti-di men sion al Hyp erbo lic Syst ems In this chap ter the pres enta tion of Cha pter 2 is gen eral ized to the mul ti-di men sion al case . 5.1 Initi al/B oun dary -val ue Prob lem be an open regi on of lRnsd Let spac e dime nsio ns. The boun dary of dina tes are deno ted be n wher e is the num ber of is deno ted by I' Spa tial coor - x = {x.} . - J. Con side r the follo wing syst em of U t -, + A. U -J -,J. m + par tial diff eren tial equa tion s: G = 0 (5.1 .1) wher e U A. -J G = u (x, t) ~j (!}I :-~ (U, (5.1 .2) t) X f x, (5.1 .3) t) (5.1 .4) oU/C lx. J (5.1 .5) nsd A. U . -J -,J Eq. 2:. -J -,J j=l A. U . (sum mati on conv entio n) (5.1 .1) is the mul ti-di men sion al anal og of (2.1 .1). (5.1 .6) 84 Eq. (5.1.1) is said to be hyperbolic if for each there exists a transformation matrix S k {k.} l E such that (5.1. 7) where A is a real, diagonal matrix. Eq. (5.1.l) is called a balance law if there exist vectors A. = ~J l If, in addition to (5.1.8), we have that G 0 , < n < j .7. such that ~-J (5.1.8) sd (5.1.l) is called a conservation law. (5.1.1) is called a symmetric hyperbolic system if A. ~J 1 < j = The initial condition for U (x, 0) < n - sd (5.1.8) (5.1.l) is x Est = (5.1.9) and boundary conditions are assumed to take the abstract form (2.l.10). 5.2 Weighted Residual Formulation In the present case the weighted residual formulation is given by 0 f ~ . (~' st t + ~j ~I j + G) d~ (S.2.1) 85 where W is typically assumed to have the following form: w = w + T. T. w. T. -J.. (5.2.2) - , J.. and ·_i ]. -A. (no sum) (5.2.3) (no sum) (5.2.4) ]. or T. -l. Eqs. = T. l. A. T - l. (5.2.1) and (5.2.2) are the multi-dimensional analogs of (2.2.3) and (2.2.4)' respectively. 5.3 Semi-discrete Equations The semi-discrete equations of Section 2.3 remain in force except for the definitions of the element arrays (cf. fined as follows e ~ab = e ~ab and e Sab which need to be rede- (2.3.5) and (2.3.8), resp.): f T aQ T.)Nb (Na I+ Na,i. i f (N (5.3.l) Qe e ~ab = a I+ N a, i T . T.)A.Nb . dQ 'J -J.. -J (5.3.2) 86 5.4 Tra nsie nt Alg orith ms The tran sien t algo rith ms of Sec tion 2.4 also per tain to the pres ent case . he ~ab The only chan ge nec essa ry is to the def init ion of the elem ent arra y whic h now take s the form (cf. (2.4 .27) ): (5.4 .l) 5.5 Sele ctio n of T. l spa tial crit eria we con side r two mul ti-d ime nsio nal gen eral izat ion s of the loca l criteri on, (2.5 .l): T. l = F Cl h/p l < i < nsd (5.5 .1) 1 < i < n (5.5 .2) and T. l whe re p. l F Cl h./p , l (no sum) l is the spe ctra l radi us of (J . l A. sd tha t is t ~1 = max 1 ~j < m IA.J (A.) I ~1 (5.5 .3) = II p II (p. p.) ~ (5.5 .4) and p h = h. l l h. p./p l (5.5 .5) l 2 J JV x. ~ l l IJ (5.5 .6) 87 Eq. (5.5 .6) hold s for isop aram etric map ping s. V , is take n in term s of the natu ral Cart unit nsd- cube . nsd For exam ple, The grad ient ope rato r, esia n coor dina tes of the bi- (5.5 .6) yiel ds the follo wing form ulas : = l h = = 2 h. = l 2 I dX/dt_; I (5.5 .7) (5.5 .8) If othe r type s of fini te elem ents are emp loye d, (5.5 .6) need s to be suit ably mod ified . temp oral crit erio n The gen eral izat ion of the glob al crit erio T. l = F n, (2.5 .8), is a L'it 1 < i < nsd (5.5 .9) The exam ples whic h follo w (2.5 .8) conc erni ng sym metr ic imp licit oper ators , incr eme ntal opti mal ity, and LaxWen drof f type meth ods, may be gene rali zed to the mul ti-di men sion case in stra ight forw ard fash ion. Rem ark The form ulas pres ente d in this sect ion of ones used prev ious ly repr esen t esth etic impr ovem ents [B7, BS, Hl2, Hl4, Hl5] . 88 CHAPTER 6 Numerical Applications in Two Dimensions Introduction 6.1 Problem Geometry and Governinq Equations We consider the problem of a thin biconvex airfoil placed in a uniform flow field. The axis of the parabolic arc Fig. 6.l(A) shows the configuration, the unifo:nn flow (non-lifting case). b where length. is aligned with the direction of denotes the ratio of the maximum airfoil thickness to the cord The notatibns used for flow variables are defined in appendix I. The subscript " 00 " refers to the free stream. Since we know that the solution will be symmetric with respect to the x l - axis, we need only consider the half plane domain. as our problem The parabolic arc bounding the airfoil in this-plane is described by the following expression: (6.1.l) The governing equations are the Euler equations described in section I.l and I.3. The ratio of the specific heats is y = (6.1.2) 1.4 (The reader should not be confused by the use of y for the transient algo- rithm parameter of chapter 2 and the present usage.) The eigenvalues of the coefficient matrices obtained from (I.3.8) by setting respectively. ~l {k , k ) = (1, 0) 2 1 and and ~ 2 can be (k , k ) = 0.1 2 1 89 (A) BOUNDARY VALUE PROBLEM I> Pr:IJ I> U1 rIJ I> U2cu I> e'° (B) b x, .. COMPUTATIONAL DOMAIN u2 -- 0 f X2 L, fJ} u =u : p= 1 L1 L2 1 e = eCl) u 2 =-4 bx 1u1r:IJ +0.5 u =0 2 Figure 6.1 x, I Boundary value problem and computational domain for thin parabolic are airfoil. 1 90 Boundary Conditions The free stream parameters are taken to be (6.1.3) l. The value of will be set according to the following formula, which depends on the specified value of the free stream mach number M00 (6.1.4) :::: M;, y(y - 1)/2 + 1 Along the x 1 - axis, outside the airfoil, we impose the following condition u 2 :::: 0 (6.1.5) On the surface of the airfoil, the velocity vector must be perpendicular to the surface normal vector. This restriction can be expressed as: (6.1.6) Assuming that the airfoil thickness is small enough, such that the uniform flow field is perturbed only slightly, u 1 can be approximated in (6.1.6) I by its free stream value (see [L4]): (6.1. 7) From (6.1.1) and (6.1.7), the boundary condition on the surface of the airfoil can then be expressed as == (6.1.8) 91 Nature of the Solution Once we set the free stream parameters as given by (6.1.3) the nature of the solution depends on the free stream Mach number and the airfoil thickness ratio. We fix the thickness ratio to be 0.10 b (6.1.9) and study the problem for two different values of the free stream Mach number. The subcritical value ~ = 0.5 results in a symmetric subsonic solution, while the supercritical value a shock around x _ = .3 . 1 M00 = 0.84 gives a solution with 92 6.2 Computational Domain and Algorithmic Features Finite Element Mesh and Boundary Conditions The computational domain is shown in Fig. 6.l(B). finite element meshes with different overall L = L = 3.5, 1 2 meshes with mesh has SN 4N the medium and fine and the coarse mesh with elements in the x -airection; in the 4N The number of nodal points are: L 1 = L 2 = 2.0. 1 elements each upstream and downstream. (16N + 1) x (4N + 1), total, and For the coarse mesh, Each x -direction, there are 2 elements across the airfoil and on the airfoil. size~: We utilize three N = 1 (SN + 1) for the medium mesh N = 2 , N = 4 . and for the fine mesh The meshes are shown in Fig. 6.2; they are symmetric with respect to the x 2 - axis. and e At the upper boundary, we impose the condition u At the left boundary we set p , 2 = 0 , which can physically be interpreted as a channel wall. Along the conditions of (6.1.5) and (6.1.8). Imposing the boundary condition of (6.1.8) along the x 1 - axis x to their free stream values. 1 - axis we take the boundary instead of on the airfoil surf ace is a standard thin- airfoil approximation. Other Algorithmic Features The integrations of all the element level vectors and matrices were performed by using 2 x ~ Gaussian quadrature rule. We used the temporal definition of the parameter the 2 A -form. T. l with For the subcritical case we also employed the The transient algorithm parameters y and a F =1 and T A A-form. were both set to unity. 93 COARSE MESH MEDIUM MESH FINE MESH 1111111111 11111111111 Figure 6.2 Finite element meshes. Coarse mesh, elements; fine mesh, 1016 elements. 64 elements; medium mesh, 256 94 Implicit methods with one itercttion were employed. The time step for each problem was chosen to be ten times an estimated critical time step for that problem. Similar to the definition of (4.4.19), the critical time step is defined here as: = min (h: IP(~.)) e,j J where the subscript j the element number. Since estimation of we use the free stream value of ~tCR (no sum) is the space dimension and the superscript p(~j) {6.2.l) J e is is unknown prior to execution, for the p (A.) -J The steady boundary condition of (6.1.8) was implemented in the same way as for the nozzle problems of chapter 4. That is, during an initial time period of certain length, the thickness ratio was taken as a linear function of time, and at the end of this time period (4 time steps) it reached its steady-state value. 6.3 Subcritical Case We compare our results, at free stream Mach number .5, to the analytical solution of the Cauchy-Riemann equations of the small-perturbation problem. Expressing the velocity vector as a sum of its free-stream value and a perturbation (uloo, 0) (ui, u;): (6 .3 .l} the analytical solution (see [L5]) can be written as: 95 ~ Su~ - iu 2 where z = xl + i8 i x~ .!. =· Ula:, b *( = (- 1) ~ and B (1 - = 1 z - z ln( z +- u2 (6.3.2) M )~ (6.3.3) 00 One can observe from (6.3.2) that, along the condition for 0.5)) 0.5 x - axis, the boundary 1 has been satisfied, while, for u 1 , the following expression is obtained: xl + o.51 = The pressure coefficient x1 - ) a.sf (6.3.4) defined as. Cp = (6.3.5) can be obtained by the following equation (see [L4]): (6.3.6) If the second order terms are neglected, then we get: (6.3.7) = Finite Element Solutions Fig. 6.3 shows the analytical and finite element (medium and coarse mesh) solutions for subcritical flow at M00 = .5 . The time steps for the regular and coarse meshes were set to 0.23 and 0.46, respectively. state solutions. The solutions a~ the er,rl of 30 steps were taken as steady For the medium mesh, conv.3rgence was achieved in about c "'d 0 COARSE MESH x 0.1 0.2 o.3 0.1 o.5 :::> Figure 6.3. GLOBAL AT A .... ..,--...- - - -......... ;..::.:..:.:::.;.::;.:.·.:.·:·:·;;::,;.;:.:·.:.::;':I·•--"~.~--.-.~---------- -O.S-0.4 -0.3 -0 . 2 -0.1 x 0.0 O.I 0.2 0.3 0.4 0.5 d_,___-?"----.~--.--~....-~-....~.......----.~--..---....----1 d d l'l d n d Steady state solutions for 10% thickness airfoil, subcritical flow, M00 = .50: Lit (medium mesh) = 0.23, Lit (coarse mesh) = 0.46. ' -o.5-0.1 -o.3 -0.2 -0.1 x ~·~·•·6.:•·•·:;a.·1.·•· .......- -...., ....,"""li:'w ........... _ .. .,.., •..,,.;.·:.;,."'""'"'.,•tt""• • "' c)"T"-------------------------------------------------. c o.o •••-~-a••.,. ...... ••.~-• {D) 2 c d+-~-T~-....~---~--..,....--....-~.......~--,r----...---....-~-t -O.S-0.4 -0.3 -0.2 -O.I 0.0 0.1 0.2 0.3 0.4 0.5 - d d l'l d n 0 ATA .... d - GLOBAL x .....~~..-~~..-~~r-~--.,.....~~~~--.~~-...~~-t -O.S -0.4 -0.3 -0.2 -O. t 0.0 0.1 0.2 0.3 O.t 0.5 :::> GLOBAL A d~~~~~·--~~~~~~~~~~~--~~~- "' (8) "'? u· ,o CL o ( c) ANALYTICAL 2 A MEDIUM MESH GLOBAL 1-+~~"T""~~ 0 ? Lil u· ,o CL c "'d (A) \.!) (Jl 97 20 steps. Figs. 6.3(A), (B) and 6.3(C), The values of A A-forms are shown in respectively. (D) and T and The finite element solutions for u 1 (U) are plotted along the airfoil. The numerical solutions are in close agreement with the analytical solution except for the variable u 1 (U) T A A-f orrn when the is used. discrepancy is the main empirical reason we have for favoring the When it comes to other variables, such as p and were observed between the solutions produced by e , A 2 and This 2 A -form. similar discrepancies T A A-forms. Remarks 1. The analytical solution predicts infinite, and thus discontinuous, values at the leading and trailing edges. 2. The computational boundaries did not seem to notably influence the solution near the airfoil. 3. It is interesting to observe that the results for the coarsest mesh are in close agreement with the analytical solution. 4. The Galerkin algorithm, on the other hand, produced highly oscilla- tory results which fed back into the operators of the problem and caused the results to diverge. 6.4 Supercritical Case We compare our solutions, at free stream Mach number 0.84, to the finite- difference solution of Barton [BS] who used the flux vector splitting scheme of Steger [S2, S3]. The finite difference grid has 97 x 33 points with 75 points along the airfoil. The boundary condition on the airfoil was taken to be the same as 98 ours, as given by (6.1.8). However, at the outer boundaries of the compu- tational domain, free stream boundary conditions were imposed in Barton's calculations. Finite Element Solutions Fig. 6.4 shows the 800-iteration solution of Barton together with the finite element solutions for coarse, medium and fine meshes. The time steps were set to 0.40 for the coarse mesh and to 0.20 for the medium and fine meshes. For the coarse mesh, the 60-step solution, and for the medium and fi~e meshes, the 120-step solutions, were taken as the steady-states. We observe that the medium and fine meshes produced very similar results and they are in agreement with the Barton solution, except for a 3-4% shift in the location of the shock front. It is known that the way the boundary conditions are imposed can change the solution considerably [BS, HS, Sl]. Because we imposed different boundary conditions at the outer boundaries of the computational domain, and especially, because we imposed a "wall" boundary condition, rather than a free stream boundary condition at the upper boundary, we are not surprised that there is some difference between the results of Barton and ours. This may be explained as follows: It is known that as the free stream Mach number gets higher, the shock front moves downstream. When we check our Mach number at the points along the upper computational boundary, we see that it is not equal to the free stream Mach number 0.84, but higher (about 0.85-0.86). with the downstream shift of the shock front we obtain. This is consistent . I • I ....-; (\] 0 ! en 0 I "'1"' . Uc). (L 0 0 ~ 0 en ~/ .r .7·> .... ..,.· ~- \ " ~' \ I \\ I~ \ I \ x 0.0 0. 1 0.2 0.3 0.4 Steady state solution s for 10% thlck~ess airfoil , supercr itical flow, M = .84: L1t (medium and fine 00 mesh) = 0.20, L1t(coar se mesh) = 0.40. -o. 5 -o. 4 -o . 3 -o. 2 -o. l Figure 6.4 'II I, 0. c ~------..--~.----....-~-----....~--~--.. .--------.---~~--t Barton (Steger Flux Vector Splittin 9) Fine Mesh Medium Mesh Coarse Mesh ~'?!.'.~·· _.~~--~...-~~---.-----.----...-~--~....- . .-l ~ l.D l.D 100 The coarse mesh (which is really much too coarse for this type of problem) did the best that could be expected. The peak point of the shock is about at the same location found by the medium and fine meshes. The coarse mesh, however, has only two elements between the peak and the trailing edge, and this is not enough for a distinct representation of the pressure profile to t~e right of the shock front. The medium and fine finite element meshes, compared with Barton's differ- ence grid, are also very coarse. Thus we can conclude that, even with rela- tively crude meshes, the finite element algorithm performed very well on this problem. 101 CHAPTER 7 Conclusions In this work we presented a Petrov-Galerkin finite element algorithm for first-order hyperbolic equation systems. The immediate purpose was to solve fluid dynamics problems governed by the conservation-law form of the compressible Euler equations. Finite element algorithms, inherently, can be easily applied to problems with arbitrary geometries and boundary conditions. This permits us to utilize the algorithms developed for complicated domains that finite difference algorithms would, normally, have difficulty with. The Petrov-Galerkin algorithm presented here is a generalization to hyperbolic systems of the streamline upwind/Petrov-Galerkin algorithm developed by Hughes and Brooks [B7, B8, Hl2, Hl4, Hl5]. We conducted an extensive stability and accuracy analysis on a linear model problem and observed that the algorithms suggested have desirable properties. Compared to the usual Galerkin algorithms they minimize spurious oscillations without loss of accuracy. We tested our algorithms on several problems with governing equations in conservation law form. Particular attention was paid to the cases with shocks. In one dimension, numerical experiments were made on linear transient, nonlinear transient, and nonlinear steady problems. We compared the numerical solutions to analytical ones and observed that they were generally in very close agreement. The algorithms handled shock fronts very satisfactorily; 102 the shock fronts were, for the most part, very crisp, with minimal spurious oscillations. The usual Galerkin algorithm on the other hand, was ineffec- tive for problems with discontinuous solutions, while for problems with smooth solutions it performed satisfactorily. In two space dimensions, we tested the algorithms on a thin biconvex airfoil problem. For both subsonic and transonic cases, the algorithms proved to be successful. For the transonic problem, the location and magni- tude of the shock front was in good agreement with Barton (Steger flux-vector splitting) solution (BS, S2, S3]. The subsonic case results were in close agreement with a linear analytical solution. For the subsonic case, we For the transonic obtained a good solution even with a very coarse mesh. case, the coarse mesh solution was not as good as the medium and fine mesh solutions, yet, was qualitatively satisfactory. Overall, the finite element algorithms suggested here performed very well for problems with smooth and discontinuous solution. of the time parameter, l The optimal selection which appears as a factor in the perturbation part of the weighting functions, needs further investigation. This needs to be pur- $Ued from the standpoint of nonlinearities and shocks which are, of course, prime concerns in solving the compressible Euler equations. We believe that, with the recent advances in the development of PetrovGalerkin algorithms, the finite element method has now become a viable alternative in computational fluid dynamics. However, the efficiency of finite element algorithms still needs to be improved, especially with respect to decreasing storage requirements. Recently, an "element-by-element" approch to the finite element formulation has been proposed by Hughes, Levit and Winget 103 (see [Hl7]). The preliminary results seem to be promising, particularly for problems with symmetric operators. Eventually, with the help of such new concepts, the finite element method can be expected to become an economically competitive and powerful analysis tool in the fluid dynamics. fi~ld of computational 104 APPENDIX I The Euler Equations I.l General Principles The compressible Navier-Stokes equations, with no source terms, can be written as a system of conservation equations: U -,t where U +:7.. (I.1.1) =O -J,J is the vector of conservation variables and the vectors, which are, in general, functions of U :!7. 's J and its spatial are flux derivative~;: = (I.1.2) If we neglect dissipative effects (i.e. conduction, viscosity, etc.) then the flux vectors are functions of 87, -J = U only: ::¥, (U) (I.1.3) -J - and (I.1.1) is called the Euler equations (or the inviscid gas dynamics equa-illions) . Let us write the Euler equations in a quasi-linear form: !:!,t + A·U -J-, J· = 0 (I.1.4) where A. -J (I.1.5) 105 Consider the linear combination of coefficient matrices: = A where k is a real vector. be a unit vector. each k (I.1.6) k.A. J-J Without loss of generality, we take k to The equation system of (I.1.4) is hyperbolic if for there exists a non-singular transformation matrix S such that the similarity transformation § diagonalizes.the matrix A. -1 A s Here = A (I. l . 7) A is a real diagonal matrix. It turns out that, for the Euler equations, this similarity transformation also symmetrizes the individual coefficient matrices simultaneously. This, in general, cannot be expected for all hyperbolic systems. To determine the transformation matrix S we need to go through the usual procedure of finding the eigenvalues and eigenvectors of the matrix (see [Tl, W2, W3]). We go through this procedure in two phases. A The first phase consists of a transformation into the primitive variables form: U,. -,t + A":u". -J-,J = 0 (I.1.8) where (I.1.9) Q = au;au~ (I.1.10) and A" ::::: k.A": J-J (I.l.11) 106 II Here II refers to the frame of primitive variables. In this frame, the coefficient matrices have simpler forms, thus, it is easier to find the eigenvalues and eigenvectors. In the second phase, the operators A ... R and -1 R ' which diagonalize are constructed and another transformation is performed: u~ = o (I.1.12) (Q R) -lA. (Q R) (I.l.13) A"".u~. + -,t -J-,J where = A"." -J = -J - - au ...;au"' = R - - (I.1.14) and (I.1.15) = Here 11 ~ " is a diagonal matrix. The of (I.1.7) is, then, equal to the product QR. refers to the frame in which transformation matrix S A"' Further, in this last frame, the individual coefficient matrices are syrnmetric: (A~)T -J = J!.~ -J (I.1.16) In the following sections, we define all the arrays involved in three, two and one space dimensions. I.2 Three-dimensional Case The conservation variables vector and flux vectors are (I.2.1) 107 r -, µ.P u~pu J l + ojl P + oj2 P u.pu J 3 + oj3 P u. (pe + p) l = ~ ~ -J ujpu 2 J p , where c5.1). u and p (I.2.2) _,,, are density, velocity and pressure respectively; and The total energy per unit mass, is the Kronecker delta. sum of the internal and kinetic energies per unit mass. state relates the pressure to the other variables. = p i = i where ~ 1~1 is the An equation of That is: (I.2.3) p(p, i,} e - e , 2 is the internal energy per unit mass. (I.2.4) If we have an ideal gas then the equation of state becomes: p Here y = {y - 1) p is the ratio of the specific heats. The coefficients matrices are: i (I.2.5) 10 8 0 0 ' + : u. J + (I .• 2 .6 ) u. J -2 oj 3 ru /2 - uj u3 6j lu 3 - 6j 3Y ul I <Yu2 - ye )u oj 2 u 3 - cj 3yu 2 -~-j oj le : - Y uj ul oj 3 u + 3 - oj 3 :Yu3 Oj 3Y u. J I ! oj2E: - Yu.ju2 oj3E: - yu ju 3 yu j 109 where y = E y - 1 (I.2.7) -2 ye - yu /2 (I.2.8) Matrices for the transformation to the frame of primitive variables are: 1 ul Q l == p (I.2.9) p u2 p 3 u2 u /2 Pu pu2 1 PU3 l/y and -, 1 - u /p 1 Q -1 = l/P - u 3 /p 2 Y"u 12 (I.2.10) l/p - u2/p l/p - yu l - Yu 2 - Yu 3 y where a blank slot indicates a zero term. The primitive variables vector is: u,. (I.2 .11) 110 The coefficient matrices in the frame of primitive variables are defined by u. J 0 oj3P oj2P oj1P u. oj 1 /P J ,. A. u. J ::::: ~J oj2/P u. J 0 2 oj2pc 2 oj3/p oj3pc 2 u. J is the acoustic speed: c where ojlpc (I.2.12) 2 yp/p (I.2.13) /.. 5 = /.. 1 - c (I.2.14) c = The eigenvalues are: k·U· J J /..4 ::::: 1..1 + c ,. The matrices for diagonalizing R ~ = r R -1 ::::: (see [ W2, W3)) are ~ 2 k3 p/ (/2 c) 3 k2 kl/12 - kl//2 k2//2 k2//2 0 k3//2 - k3!12 0 pc//2 pc//2 p/ (/2 c) ..., kl k 0 - k k3 0 - k - k2 kl 0 0 kl 0 k3 - k 2 - k /c 1 k2 - k3 0 kl k 2 /c k3 k2 - k 0 k /c 3 0 kl//2 k2//2 k3//2 1/(/2 pc) 0 - k1/l2 k2//2 k3//2 1/(/2 pc) 1 1 (I. 2 .15) 2 2 2 (I.2.16) 111 The symmetrized coefficient matrices are [W2, W3] A. ~J = •l. J u. J oj2k3 - <5j3k2 oj3k'1 - oj1k'3 oj1k"2 - oj2k"1 oj2k3 -- 6j3k2 oj3:kl. - ojlk3 6 oj2k'J - oj3k'2 oj2k'3 - oj3k"2 0 ji'1 - 0 j1k'3 oj3k1 - ojlk'J oj1k"2 - oj2k"1 oj1k"2 - oj2k"1 jl2 - 6 j2k1 (I.2.17) where k. J I.3 = ck .//2 (I.2.18) J Two-dimensional Case The conservation variables vector and flux vectors are defined by u (I.3.1) 112 u.p J ujpu Ul'. = "'] 1 + oj1 p (I.3.2) ujpu 2 + oj 2 p uj(Pe + p) The coefficient matrices are defined by: 0 0 0 - 2 jlYu /2 - ujul ojlul - oj1Yu1 6j2ul - oj 1:Yu 2 oj1Y (I.3.3) + u. J . ---·- ···--- ·- -· . -·-- - 2 oj2yu /2 - uju 2 6jlu2 - oj2Yu1 6j2 yu2 oj2u2 6j2y + u. J --·--- c:Yu 2 -· -··· - ye)u.J OjlE - yujul oj 2 E - yu.u J 2 yu. J The matrices for transformation to the frame of primitive variables are: 1 ul p Q u2 2 u /2 (I.3.4) p pul u2 l/Y 113 and l -1 l/p - u /P l = Q (I.3.5) l/p - u2/P - 2 yu /2 - yu - yu l y 2 'I'he priltlitive variables vector is u .. (r.3.6) The coefficient matrices in the frame of primitive variables are defined by u. J .... A· ~J t\ 1 P o.)2 p 0 u. 6j 1/P J = u. oj 2/P J 0 ojlpc 2 oj 2 pc (I.3.7) 2 u. J The eigenvalues are: Al = A2 = k.u. J J A4 = A.l - c A3 = Al + c The matrices for diagonalizing A (I.3.8) are [W2, W3]: l 0 p/ {/2 c) p/{/2 c} 0 k2 kl//2 - kl/12 k2//2 k2//2 pc.;//2 pc/v'2 R 0 - k 0 0 l (I.3.9} 114 and R-1 1 0 0 0 k2 - k = 1/c 1 0 2 0 kl/12 k2//2 1/ (/2 Pc 0 - kl/12 - k2//2 1/ (/2 pc) (I.3.10) The symmetrized coefficient matrices are defined by A7° ~J = r,, ,... 0 v ~. 0 J I.4 0 u. 0 oj1k2 - oj2k"1 0 oj1:k2 J oj1k2 - oj2'k1 8 j1 k2 u. + c k. J - oj2k'1 - 8 j2kl (I.3.11) 0 J u. 0 J - c k.J One-dimensional Case: The conservation variables vector and flux vector are 1 u p = u (I.4 .1) e up F = upu + p u (pe The coefficient matrix is + p) (I.4.2) 115 l 0 2 3)u /2 (y - A (y 3)u ye - 3yu 2 /2 - fru 2 - ye) u 0 (I.4.3) yu Matrices for transformation to the frame of primitive variables are p (I.4.4) p u and r Q -1 1 - u/p l/p - - yu yu 2 /2 (I.4.5) y The primitive variables vector is p u (I.4.6) p The coefficient matrix in the frame of primitive variables is ~ A = u p 0 0 u l/P 0 pc A3 = Al - c 2 (I.4.7) u The eigenvalues are: Al = u A2 Al + c (I.4.8) 116 The mat rice s for dia gon aliz ing -R A~ are [W2, W3] : 1 p/(1 2 c) p/( /2 c 0 1/1 2 - 1/12 0 pc/ /2 pc/1 2 1 0 0 1/1 2 1/ c/2 pc) 0 - 1/1 2 1/ c/2 pc) (I.4 .9) and R -1 = The diag ona lize d coe ffic ien t mat - l/c 2 (I.4 .10 ) .cix is \ A = (I.4 .11) 117 APPENDIX II Stability and Accuracy Analysis of Algorithms for the One-dimensional Linear Parabolic Equation Development of the Tools for the Analysis II.l II.1.1 Introduction Model Problem: Diffusion Equation One-dimensional diffusion of a function U(x,t) is governed by the following parabolic equation: u,t where K K U,XX is the diffusion coefficient. = (II.1.1) 0 An initial condition of the follow- ing form is assumed: U (x, O) (II.1.2) The analysis, arguments, parameters, etc., defined in this appendix are mostly the same, or very similar to, those of chapter 3; we will just state state the ones which are different. Exact Solution The exact solution is obtained by separation of variables, as was done in chapter 3. (II.1.2). The spatial component of the solution is given by The damping coefficient and the frequency of the temporal com- ponent are: - (l;, w) = V = - (K 2 k , 0) (II.1.3) 118 II.1.2 Finite Element Solution Spatial Discretization For the weighted residual formulation of the problem we use the usual weighting function which results in the (Bubnov-) Galerkin formulation. is w= w That (II.1.4) The resulting semi-discrete equation Mv + Cv = (II.1.5) 0 is obtained by following a procedure which is very similar to the one we followed in chapter 2. The element level matrices e m and c e are defined differently: (II.1.6) = J KN (II. l . 7) a,x Nb,x dD Qe Assuming constant K and mesh spacing h , and going through the same arguments that we went through in the corresponding parts of chapter 3, the jth equation of the system of equations can be written as: (II.1.8) thus: .b_ D E T + K 2 D E T 2 -r - h -2 = 0 (II.1.9) 119 We define the scalars M and tit c C corresponding to M = = CK and C as: !_ D E 2 -r - (II.1.10) 2 K tit h The non-dimensional parameter M 2 E>2 (II .1.11) !E is defined as: = We note that K/h has units of velocity; thus (II.1.12) CK can be regarded as a "diffusion Courant number". The problem is reduced to solving the following ordinary differential equation . MT + C T = 0 (II.1.13) Numerical Frequency and Damping Coefficient Going through the same steps we went through in chapter 3, we obtain e vti.t = tr A = 1 + i ((1 - R) S - 1) (II. l.14) where R = Q LltC M = (MR + a L\t CR)/M (II.1.15) ~,t cR/M (II .l .16) !22 ~ (II.1.17) ML + a Llt CL (II.1.18) CK 120 Then ln(tr A) For comparison, we need to express v6t (II.1.19) in terms of the same dimen- sionless parameters; from (II.1.3): = - l2 cK q 2 v6t (II.1.20) We define the analytical and numerical amplification factors, z z and as: z = e z = e v6t (II.l.21) v6t (II.l.22) II.2 Unified Analysis of Algorithms II.2.1 Introduction When we study the stability and accuracy of algorithms for the diffusion equation, we classify the algorithms considered as implicit and explicit types. This time we do not use Petrov-Galerkin algorithms, therefore, the type of implicit algorithm will solely depend on the parameter r (defined in section 2.6) which specifies the way mass terms are integrated. In the diffusion problem, we observe from (II.l.21) that the exact a.mplification factor has no imaginary part and O < Z < l plification factor has no imaginary part either. make surA that the numerical solution is stable. The numerical am- However, we also have to 121 Graphically, we compare the exact and numerical amplification factors and inspect the ratio of the numerical damping coefficient to the exact damping coefficient for several algorithms. Obtaining closed form expressions for Z and the damping ratio is relatively easy, yet, these expressions are not in a simple form. expand the damping ratio in 6t and h , ~/~ One can and observe the temporal and spatial accuracies of the algorithms. The Galerkin algorithms use the integration stencils: .... \ 2(r, l-2r, _, D = _r (- 1 , 1 , e2 Then the scalars M and 2 - 12 (II.2.1) (II.2.2) become: 6tC = M 1 = 6tc - 2rV (II.2.3) CK v (II.2.4) - cos q (II.2.5) when v II.2.2 = 1 Implicit Algorithms For implicit algorithms, the amplification factor from (II.1.14) is z 1 - Q (II.2 .6) The stability limits are determined by the condition: 121 < l (II.2.7) 122 Z< 1 The inequality dictates that 1 + 2aC > K (II.2.8) 4r This condition is automatically satisfied as long as inequality Z> - l 1 r < -4 long as 1 < -4 The second requires that: cK V(2a - 1) Since r 4rv - 2 > (II.2.9) from the first inequality, this condition is satisfied as 1 a~2· We summarize the unconditional stability conditions for the implicit algorithms as: r < l. -4 a>! and - (II.2.10) 2 For the accuracy analysis, we can employ asymptotic expansions of the damping ratio in while, for expansion in (f) h = =0 we first set h -~/~ We expand h . and 6t l + CK k 2 2 (K + 0 (6tJ) k 6t) 6t = 0 . 6t after setting h = 0 These expansions are: (~- - a) 6t) + in 2 ( t + a.2 - a.) ..., (f) 6t 0 l + (k h) 2 + (r - l~) 4 O(h ) Observe that from (II.2.11), accuracy is time. (II. 2 .11) (II.2.12) a= ~ is the condition for second-order From (II.2.12), we see that second-order accuracy in space is automatic, and that the value r = 1 ;12 results in fourth-order accuracy. 123 II.2.3 Explicit Algorithms Explicit 1-pass Algorithm For this algorithm, the amplification factor is: z = 1- Q The inequality Z< 1 V = 2 ; when is automatically satisfied since Z> - 1 second inequality then, (II. 2 .13) dictates that CK < l cK v -> 2 • CK V > O . The The worse case occurs is the stability limit. For the accuracy analysis, the expansions of the damping ratio in and h ~t are: h 0 (II.2 .14} = (II.2.15) Thus the 1-pass algorithm is first-order accurate in time and secondorder accurate in space. Explicit 2-pass Algorithms The amplification factor for these algorithms is Z = l + Q(R - 2) The inequality Z< l (II. 2 .16) requires that CK V(- 1 - 2rV +Ci.CK V) < 0 this leads to < 1 + 4r 2Cl (II. 2 .17) 124 The inequality -z > - l dictates that CK V(l + 2rV - Cl.CK V) - 2 < 0 ; one can show that this is satisfied as long as a> l -2 and r l 2. 4 • Combin- ing these inequalities with the inequality of (II.2.17), we can state that the maximum value CK The expansions in WI a = ~ can assume is 2 ,. corresponding to h = 0 6t and l + (K k h 2 and r =~ . are: 6t) . (} - C) 1 - Ct ) + (K k 2 6t) 2 ( 3 (II .2 .18) WI ~t = =0 l + (k h) 2 (r - {-2 ) (II.2.19) a =~ Thus, with in the irnplici t case' second-order accuracy is achieved in time; and, as r = ~ can raise the spatial accuracy from second to fourth-order. s ununary II . 3 Asymptotic expansions, stability limits determined analytically, and graphical representations of Z (figure II.l) and - l;/E;, (figure II. 2) are utilized for the stability and accuracy analysis of the algorithms considered. Stability We observe from Fig. II.l that for and q e]o, ~[ , Cf. = 2l , CK = 0.2, 0.4, 0.6, 0.8, 1.0 all the algorithms considered had amplification factors which remained in the interval [-1, +l]; that is, all were stable. This can also be deduced from the stability limits determined analytically. For the implicit algorithms (GC and GL) we have stability as long as 125 !Al '=! EXACT SOLUTION Cl GC !MPUCIT CBI -"'~~::~~······· 0.'~··· a::o"' 0 .... ........ ., .. ~ a:.. Cl -...a: - ...• a:.. .... a: .4 .2 ·····--·. ~ ...J Cl o.o o.s ~ 1.5 1.0 2.0 DIMENSIONl.ESS WRVE NUMBER !Cl ···.......... ... ' ........ '=! I 2.s l.O o.o o.s Cl 1.0 1.s 2.0 Dl MENS IONLESS HAVE NUMBER ..... rn• GL IMPLICIT 25.,;"' "' !!§.,; ff lE "- Cl ..... -.... t::~ c: "'! ad a: c.J -'j' - ...J ~ ~ ' o.o ~ o.s 1.0 1.5 2.0 DIMENSIONLESS WAVE NUMBER 2.s ~~ ::S. ··-... - --...... ~<~·-..... '-... Ill !5 ci lE "5.,; ... ... ................ ··.. . ...... Cl -·---·- ···-- , ., .... .,'".,""·'"....."' .......... ' ,~·- ·::--,. . --·-····~.-........... -..........__ cc """-..., ."" .......___ ......___ 2.0 2.5 '- .............. ·- I a.o a.s 1.0 1.5 0 IMENS IONLESS HAVE NUMBER 3.0 GL E2 CF"l ";;:;~~~;=~..:: "' ct:.,; 0 ..... c.J a:.. Cl z· Oe -....a: c.J c.J t: "? :i 'i' t: "! '.:j'i' ~ ............................ a: ' ...... .... ~ '=! Cl ·:·:~~ - ..- ..---· Cl I 3.0 GC E2 !El ---·~. ·· z· 0 Cl 5 ~ "'."";"' ""' -.........--, '·0·........ ..... a:.. Cl 3.0 2.s GC CANO GLl El .,:~~:::·::~~~- - ..... ··--....__ .... _·- . t: U'! -<;> ---- ~ I ··· ..... c.J .6-- -'i' ~ z· o"' --·- .a c.J ··· .., ' a: a:.. Cl CK 1.0 ·· ..... "" ""···...····... ""'""'·'--........... ......... "··· ....... ....... -.........._ ....0 c.J ao . '"'~···... ·. "' ct:.,; ~ '=! I o.o o.s 1.0 1.s 2.0 DIMENSIONLESS WAVE NUMBER Figure II.l 2.s 3.0 Cl I o.o o.s 1.0 1.5 2.0 2.5 DIMENSIONLESS WAVE NUMBER Amplification factors for heat equation algorithms. 3.0 126 l r < - l 2 a> - and -4 - both the consistent mass With rithm (GC) and the lumped-mass (r = 0) al go- algorithm (GL) are therefore stable. For the explicit 1-pass algorithm (El) , we condition is c K< l this clearly. We reca11 that (independent of c K r < l and fou~d that the stability Fig. II.l (D) also shows a) ; z is the condition for not to drop below - 1. O; in Fig. II.l (D) we observe that the curve for just touches the stability limit at = q CK = l 'TT. For the explicit 2-pass algorithms (E2) , the stability requirements C are K < (1 + 4r) rithm with ~ r 4l . (r = The consistent mass satisfies the restrictions on the values of as a restriction on rithm with l a.:_ 2 , a = 21 cK The lumped mass we get has a restriction !) al go- 6 r and ,. (r = 0) algo- cK< 1.0 . Accuracy For the accuracy of the algorithms, Fig. II.2 provides information for l a=2, = 0.2, 0.4, 0.6, 0.8, 1.0 and points away from q = O asymptotic expansions on the other hand provide information needed neighborhood of implicit underdamps. damps. in the q = O . We observe from Fig. II.2 that the GL The Similarly, GC GC - E2 implicit algorithm overdarnps wmile overdamps while GL- E2 under- The explicit 1-pass algorithm, on thP. other hand, overdamps or under- damps depending on the value of damping occurs, is about CK CK . The crossover point above which over- 0.4 . The asymptotic expansions show that all the algorithms considered have second-order spatial accuracy for l r -:f 12 For r 1 12 ' = - the implicit 127 0 ~ "! ill :i:i o~ -""" oE! -a: """ Cic cc 1.0- z.., .8-·.6--- 0::: C!l "' ii: c:i l:: a: 0::: C!l "' Z" 4--- a: c:i .2 ············ a: l:: 0 05! 0 a.oo a.25 a.so a.15 1.00 1•.:is DIMENSIONLESS WAVE NUMBER :ii .; a.DO I.SO GC !ANO" GLJ El !Cl -·------./ ------ ---- - / / / I 1.511 GC E:2 S! I iq _.,..,. o~..L.~-------------· -- ...•...---····----············· o~ +----!Cl!!!!!!-!§':::;:::::::::-.: ·-·--·--·--··--w·-------·~··-• -""" o.so 0.75 1.00 1.25 DIMENSIONLESS WAVE NUMBER 0.25 !Ol Ii! iq GL IMPLICIT !Bl GC IMPLICIT !Al a: 0::: a:~leQ lf c:;r 0 a.oo a.25 a.so a.15 1.00 1• .:is OIMENSIONL.tSS HAVE NUMBER !El o.oo t.50 a.25 a.so o.7s 1.00 1.25 DIMENSIONLESS WAVE NUMBER GL E:2 ~ ......~~~~--~--~~~--~~~~~ o.oo Figure II.2 0.25 a.so 0.75 1.00 1.25 1.50 DIMENSIONLESS WAVE NUMBER Damping Ratios for heat equation algorithms. I.SO 128 and explicit 2-pass algorithms attain fourth-order accuracy. considered have first-order temporal accuracy for . .+ -l G, T 2 • The algorithms For l a. = 2 , implicit and explicit 2-pass algorithms attain second-order accuracy. the 129 APPENDIX II I Stability and Accuracy Analysis of Algorithms for One-dimensional Linear Second-order Hyperbolic Equation III.l Development of the Tools for the Analysis III.1.1 Introduction Model Problem: Wave Equation One-dimensional, undamped wave motion is governed by the following second-order hyperbolic equation for U(x, t) : K U u,tt l where (K)~ ,xx = 0 (III.1.1) is the wave-propagation speed. The initial conditions associated with the problem are assumed to be of the form: u (x, . U (x, 0) = (III.1.2) 0) = (III.1.3) We note that these could be viewed as Fourier components of a general set of initial conditions for U and U • Exact Solution For constant propagation speed, assuming a solution of the form U (x, t) X(x)T(t) leads to the following spatial and temporal components for the function (III.1.4) 130 U(x, t ) : where the constants A 1 T (t) = and A 2 e = X(x) ikx (III.1.5) (III. l . 6) depend on the constants u0 and v0 and (III.l. 7) ~ We define the damping coefficient (- ~' W) and the frequency w as = V = i/K k (III. l . 8) The exact solution is seen to have no damping. III.1.2 Finite Element Solution Spatial Discretization For the weighted residual formulation of the problem, we use the usual weighting functions (W = W) which lead to the (Bubnov-) Galerkin fol'.'l!lulation. The resulting semi-discrete equation .. + M v C v = (III. l . 9) 0 is obtained by following a procedure very similar to the one we followed in chapter 2 and appendix II. The element level matrices Me and Ce are de- fined as: == (III.1.10) (III.1.11) 131 Going over a sequence of steps as we did in chapter 3, we end up with the following ordinary differential equation .. M T Here the scalars M and C + C T = = l 2 Pr E are M 2K lit 2 lit c (III.1.13) 2 For the purpose of analysis, we assume constant CK cK K and h . The is defined as = Here, (III.1.14) !?2 E h2 dimensionless parameter (III. l.12) 0 = (III.l.15) is a Courant number based on the propagation speed /K . Time Integration The ordinary differential equation of (III.1.12) can be solved by a family of time integration schemes described in [Hll]. We adopt the representation y and let Y [T, denote the approximation to -n Given = Y to calculate T, T] Y (III.1.16) at the nth time step. we go through a predictor phase and an iterative phase -n v :'.:.n+l · In the predictor phase, we calculate the zeroth iteration value, Y(O) -n+l ' by the following operation: 132 p y (III.1.17) - -n Here is the predictor matrix defined as~ P p where tit = Lit 1 tit (1 - 0 1 (1 - a.) Lit 0 0 0 2 2 (III.1.18) S are the Newmark parameters which' a , is the time step and 213) control the stability and accuracy of the algorithm. The iterativ.e phase starts with the zeroth-iteration value and continues according to the recurrence rule below: Given y (i) -n+l solve the following system for .. ( i) ML tiT n+l + CL Li ( i) = Tn+l = a. tit L'lTn+1 tiT(i) n+l == B tit 2 6Tn+l .. (i+l) n+l = .. ( i) T n+l T •(i+l) n+l T (i+l) n+l The superscripts L (III.l.21) .. (i) L'lTn+l T (i) n+l and R + (III.1.22) (III.l.23) + T (III.l.19} (III.1.20) .. ( i) + -n+l (i) ~(i)) - (MR;rn+l + C n+l .. ( i) L'lT (il n+l y (i+l) /::,.T (i) n+l (III.l.24) refer to the left and right-hand sides of 133 (III.l.19). M and C We reserve the option of having different evaluations for on different sides, as in chapter 3. The recurrence rule defined above can be expressed as: (i) = J y-n+l (i+l) -n+l y J (III.1.25) is the iteration matrix: M J = 1 - B ilt 2 cR 0 .6.t CR 1 CR 0 Cl (III.1.26) M ....J where (III.l.27) M Combining the predictor and iterati11e phases: !n+l A Here S :::: AY (III.1.28) - _n = (III.1.29) is the number of iterations performed in the iterative phase. Exploiting the fact that det P = 0 , det A det JS det P we determine that at any time step, n+l , of and . Tn+1 . Tn+l 0 (III.1.30) can be expressed in terms That is: .. T n+l (III.1.31) 134 We further observe from the structure of the matrix T(i+l) n+l and .. (i+l) Tn+l (III.l.12). and P This observation suggests that we extract the submatrices from the matrices = J [ 1 p - = rl P and J M and rewrite This is, of course, a con- T n+ 1 . term in the ordinary differential equation T sequence of the absence of the of • ( i) have no dependence on that the iterates J M - l 6.t 0 0 respectively: S tit 2CR - c - S 6t 2MR] - R R (l - 26) 0 M - M L'.lt2 ] 2 (III.l.32) . (III.1.33) (III.l.28) in a different form: (i) Tn+l T(i+li n+l == .. (i+l) T n+l (III. l. 34) J .. ( i) T n+l y = J8 p- -n . T Eq. J n+l .. == Tn + (1 - a)6.t Tn + a 6.t Tn+l (III.1.35) (III .l.36) (III.l.35) reveals that the relation of (III.l.31) degenerates to the form: (III.1.37) with µ being 135 µ where = (III.l. 38) are the components of and -S J By adopting the equivalents form of (III.l.35)-(III.l.36) , we need only to calculate the Sth power of a 2 x 2 matrix instead of a 3 x 3 matrix. Besides, the only difference between the matrix of chapter 3 is that the term S6t 2 here and the matrix J replaces the tenn in mind, by analogy, we write the eigenvalues of J a6t . J With this from (3 .1. 57) -( 3. l. 58') : l (III.1.39) = (III.l.40) By way of a procedure similar to the one in chapter 3; -s J == 1 (III.1.41) where b == 1 - (1 - R) 5 (III.1.42) R (III .l .43) R It is interesting to note that the only difference betTveen the matrix J and its Sth power is the insertion of a accounts for the number of iterations. µ From b term which (III.l.38): (III. I. 44) 136 and from (III. 1.35) -(III.l.37) {::::} ~ ~·{ :: } A* with - (III.1.45) given as -r p A* == r- tit [ ll.tµ[ (1 - Cl) + arp] (1 + a)r i1] (III.l. 46) where µ == ll.t2µ 2 l + (.0.t /2) (1 - 213) µ p = and -r == (M - B b lit 2cR)/M (III.1.47) Numerical Frequency We are mainly interested in the invariants of the matrix A* which are: tr A* = 2 - E(l + 2CI.) (III.l.48) det A* = 1 + E(l - 2a) (III.l.49) where E = b 2 (III.l. 50) ~ The eigenvalues and z2 A* of the matrix are found by solving the quadratic equation z2 - tr A*Z + det A* 0 (III.1.51) 137 Assuming a solution of the form: T (III.1.52) = n one can show that the eigenvalues of and v2D.t A* are related to v1 6t by (III. l . 53) e The eigenvalues can be used for comparing the numerical v to the v . exact z1 We now have parameters q and same parameters. z2 and CK . calculated in terms of two dimensionless We also need to express v in terms of the From (III.1.8): V~t = (0, 6t./K k) = (O, CK q) (III.1.54) while = III.2 (III. l . 55) Unified Analysis of Algorithms III.2.1 Introduction Stability and accuracy analyses of algorithms for the wave equation are made for the implicit and explicit cases. are not introduced. Petrov-Galerkin algorithms The type of implicit algorithm depends only on the r term (defined in section 2.6). The Galerkin algorithms utilize the integration stencils 2 [ r, 1 - 2r, r J D ~r = 1 1 [-2,l,-2] (III.2.1) (III.2.2) 138 The scalars M and 2 6t c then become: 1 - 2rV = M 2 (III. 2. 3) V (III.2 .4) 1 - cos q (III.2.5) = 2c K where v For consistency with the exact solution, we need to have Z Z This condition is satisfied if 2 = (III.2.6) a = 1 2 For stability, we require that only when det A* = 1 . = l 2 v ~t lz1 1 , lz 2 l .:'.:_ 1 ; The eigenvalues z-1 and but this is possible z- 2 are given by the expression: = tr A* 2 tr ~* 2 = (III. 2. 7) where Clearly < 0 1 - E (III. 2. 8) E(E - 2) (III.2.9) if and only i f 0 < E < 2 (III. 2 .10) 139 < 0 This is the stability condition, because, if jz 1 Otherwise, we would have either III.2.2 1 or lz 2 j then greater than unity. Implicit and Explicit 1-pass Algorithms For an implicit or explicit 1-pass algorithm, the term E is 2 = E cK v (III. 2 .11) (1 - 2rV) + 2Sc 2 v K In this expression, if we set r = B= 0 , explicit 1-pass algorithm (Warning: the algorithm in general). we get the value of E for the This applies to (II.2.11) only, but not From the inequality of III.2.10 we determine the stability limits. One can show that E > 0, .if and only if c2 It is quite clear that if r 4r - 1 > K < -l , -4 (III.2.12) 48 then this condition is satisfied uncon- ditionally. E < 2 The inequality implies the following condition: (III.2.13) If B -> .!. 4 this inequality dictates that: I c2 K Clearly, if r 1 < - - 4 > 4r - 1 4B - 1 (III. 2 .14) this condition is satisfied unconditionally. If not, I then we need to satisfy both (III.2.12) and (III.2.14). that (III.2.14) implies (III.2.12). One can show Therefore (III.2.14) is the key condition. 140 If 1 B< 4 , we cannot satisfy (III.2.13) unless this happens, the restriction on CK is given as: c2 < 1 - 4r 1 4S K For the implicit algorithm with ditional stability. B= 0 , and B = 41 Provided that we attain uncon- For the explicit 1-pass algorithm with r =O and the condition of III.2.15 becomes < III.2.3 . (III.2 .15) = _!_6 r 1 r ::_ 4 1 . (III.2 :16) Explicit 2-pass Algorithms For 2-pass al-gori thms (III.2.17) One can show that E > 0 if and only if < For E < 2 , l B= 4 ' (III.2 .18) 46 the following condition must be met: B > With 4r + l both and r (1 + 4r) 16 =0 2 (III.2.19) satisfy this condition. Then the inequality of (III.2.18) implies that < 5 3 < l for . for (III.2.20) r =0 (III.2.21) 141 III. 3 Summary According to the stability guidelines of the previous sections, with 1 a = 2 and 1 B= 4 ' all the algorithms considered are stable (with no modulus CKE [O, l] error) for the ranges of and q E [O, ~] The accuracy infor- mation is provided by the graphs of Fig. III.l. The implicit, lumped mass algorithm (GL) becomes less accurate as increases. CK The implicit, consistent mass algorithm (GC), on the other hand, maximizes its accuracy around CK = 0.6 The explicit 1-pass algorithm (El) satisfies the unit CFL condition defined in section 3.2.4. This is as we would expect, because if A* tr -"'-2 = = 1 - E = 1 - E(E - 2) v cK = 1 (III.3.1) ::::: V (V - (III.3.2) 2) and therefore the eigenvalues are ~\,2 Jv<v - 2) :::: (1 - = cos q + i sin q V) + - +iq e- (III.3.3) This is the same expression as for the exact solution. El increases as C K The accuracy of increases within the stability regime. The explicit 2-pass (E2) version of GC is more accurate than the 142 GC IMPLICIT GL IMPLICIT c,~ r.o-- .a-·.s--- a.a J.a 1.0 l.S :z.a 2.5 DIMENSIONLESS HAVE NUMBER O.S a.a GC !ANO Gl.l El a.a 1.0 1.s 2.a DIMENSIONLESS WAVE NUMBER 2.s 3.a GC E2 ~ (Q) ~.,...~~~~~~~~--~~~--~~--- a.s 1.0 1.s i.a :z.s DIMENSIONLE:SS WAVE NUMBER ~ !El o.o Figure III.l o.s GL o.s o.o 3.a o.s 1.0 1.5 2.0 2.5 OIMENSlONLE:SS WAVE NUMBER E:.2 1.0 1.s 2.0 OIHENSIONLE:SS WAVE NUMBER 2.5 3.0 Frequency ratios for wave equation algorithms. 3.0 143 implicit version until about CK = 0.8 . At this point they almost have the same accuracy, and after that the implicit version is more accurate. GC-E2 maximizes its accuracy between c K = 0.6 and 0.8 . 144 REFERENCES Bl. A.J. Baker, Research on Numerical Algoritluns for the Three-Dimensional Navier-Stokes Equations, I. Accuracy, Convergence and Efficiency, Technical Report AFFDL-TR-79-3141, Wright-Patterson Air Force Base, Ohio, December 1979. B2. A.J. Baker and M.O. Soliman, "On the Utility of Finite Element Theory for Computational Fluid Dynamics,", Proceedings of the 5th AIAA Computational Fluid Dynamics Conference, June 1981, Palo Alto, California. B3. A.J. Baker, Private Communication, 1981 B4. W.F. Ballhaus, A. Jameson and J. AJbert, "Implicit Approximate-Factorization Schemes for Steady Transonic Flow Problems", AIAA Journal, Vol. 16, No. 6, 1978, pp. 573-579. BS. J. Barton, Private Communications, 1982. B6. R.M. Beam, R.F. Warming, and H.C. Yee, "Stability Analysis of Numerical Boundary Conditions and Implicit Difference Approximations for Hyperbolic Equations", Proceedings of Numerical Boundary Condition Procedures Symposium, October 19-20, 1981, Ames Research Center, NASA, Moffet Field, California. B7. A. Brooks and T.J.R. Hughes, "Streamline Upwind/Petrov-Galerkin Methods for Advection Dominated Flows", Proceedings of the Third International Conference on Finite Element Methods in Fluid Flow, Banff, Canada, 1980. B8. A. Brooks and T.J.R. Hughes, "Streamline Upwind/Petrov-Galerkin Formulations for Convection Dominated Flows with Particular Emphasis on the Incompressible Navier-Stokes Equations", to appear in the Proceedings of FENOMECH '81 and Computer Methods in Applied Mechanics and Engineering. Cl. I. Christie, D.F. Griffiths, A.R. Mitchell and o.c. Zienkiewicz, "Finite Element Methods for Second Order Differential Equations with Significant First Derivatives", International Journal for Numerical Methods in Engineering, Vol. 10, 1389-1396, 1976. C2. R. Courant and K.O. Friedrich, Supersonic Flow and Shock Waves, Interscience Publishers, New York, 1948. C3. R. Courant and D. Hilbert, Methods of Mathematical Physics, Interscience Publishers, New York, 1962. \ 145 Dl. J.E. Dendy, "Two Methods of Galerkin Type Achieving Optimum L 2 Rates of Convergence for First Order Hyperbolics", SIAM Journal of Numerical Analysis, Vol. 11, pp. 637-653, 1974. Hl. J.C. Heinrich, P.S. Huyakorn, O.C. Zienkiewicz and A.R. Mitchell, "An 'Upwind' Finite Element Scheme for Two-Dimensional Convective Transport Equation", International Journal for Numerical Methods in Engineering, Vol. 11, pp. 134-143, 1977. H2. H.M. Hilber, T.J.R. Hughes and R.L. Taylor, "Improved Numerical Dissipation for Time Integration Algorithms in Structural Dynamics", Earthquake Engineering and Structural Dynamics, Vol. 5, 283-292, (1977). H3. T.L. Holst and W.F. Ballhaus, "Fast, Conservative Schemes for the Full Potential Equation Applied to Transonic Flows", AIAA Journal, Vol. 17 1 No. 2, February 1979, pp. 145-152. H4. T.L. Holst and D. Brown, "Transonic Airfoil Calculations Using SolutionAdaptive Grids", Proceedings of the 5th AIAA Computational Fluid Dynamics Conference, June 1981, Palo Alto, California. HS. T.L. Holst, Private Communications, 1982. H6. T.J.R. Hughes, A study of the One-Dimensional Theory of Arterial Pulse Propagation, Report No. UC SESM 74-13, Structural Engineering Laboratory, University of California, Berkeley, 1974. H7. T.J.R. Hughes, "Stability of One-Step Methods in Transient Nonlinear Heat Conduction", Transactions of the Fourth International Conference on Structural Mechanics in Reactor Technology, San Francisco, California, August 1977. H8. T.J.R. Hughes and J. Atkinson, "A Variational Basis for 'Upwind' Finite Elements", IUTAM Symposium on Variational Methods in the Mechanics of Solids, Northwestern University, Evanston, Illinois, September, 1978. H9. T.J.R. Hughes and W.K. Liu, "Implicit-Explicit Finite Elements in Transient Analysis: Stability Theory", Journal of Applied Mechanics, Vol. 45, (1978), 371-374. HlO. T.J.R. Hughes and W.K. Liu, "Implicit-Explicit Finite Elements in Transient Analysis: Implementation and Numerical Examples", Journal of Applied Mecha~ics, Vol. 45, (1978), 375-378. Hll. T.J.R. Hughes, K.S. Pister and R.L. Taylor, "Implicit-Explicit Finite Elements in Nonlinear Transient Analysis", Computer Methods in Applied Mechanics and Engineering, 17/18 (1979) 159-182. 146 Hl2. T.J.R. Hughes and A. Brooks, "A Multidimensional Upwind Scheme with no Diffusion", in AMD Vol. 34, Finite Element Methods for Convection Dominated Flows, T.J.R. Hughes (ed.), ASME, New York, 1979. C~osswind Hl3. T.J.R. Hughes, "Implicit-Explicit Finite Element Techniques for Symmetric and Nonsymmetric Systems", Proceedings, First International Conference on Numerical Methods for Non-Linear Problems, Swansea, U.K., 1980. Hl4. T.J.R. Hughes and A. Brooks, "Galerkin/Upwind Finite Element Mesh Partitions in Fluid Mechanics", pp. 103-112, in Boundary and Interior Layers - Computational and Asymptotic Methods, J.J.H. Miller (ed.), Boole Press, Dublin, 1980. Hl5. T.J.R. Hughes and A. Brooks, "A Theoretical Framework for Petrov-Galerkin Methods with Discontinuous Weighting Functions: Application to the Streamline Upwind Procedure", to appear in Finite Elements in Fluids, Vol. 4, R.H. Gallagher (ed.), J. Wiley and Sons. Hl6. T.J.R. Hughes, "Analysis of Transient Algorithms with Particular Reference to Stability Behavior", to appear in Computational Methods in Transient Analysis, T. Belytschko and T.J.R. Hughes (eds.), North-Holland Publishing Company, Amsterdam. Hl7. T.J.R. Hughes, I. Levit and J. Winget, "An Implicit Unconditionally Stable Element-by-Element Algorithm for Heat Conduction Analysis", to appear in Journal of the Engineering Mechanics Division, ASCE. Ll. P.D. Lax, "Weak Solutions of Nonlinear Hyperbolic Equations and Their Numerical Computation", Communications on Pure and Applied Mathematics, Vol. 7, pp. 159-193, 1954. L2. P.D. Lax, "Hyperbolic Systems of Conservation Laws II", Communications on Pure and Applied Mathematics, Vol. 10, pp. 537-566, 1957. L3. P.D. Lax "Nonlinear Hyperbolic Systems of Conservation Laws", Nonlinear Problems, R.E. Langer (ed.), pp. 3-12, Madison, The University of Wisconsin Press, 1963. L4. H.W. Liepmann and A. Roshko, Elements of Gasdynamics, J. Wiley and Sons. New York, 1957. LS. H. Lomax and E.D. Martin, "Fast Direct Numerical Solution of the Nonhomogenous Cauchy-Riemann Equations", Journal of Computational Physics, Vol. 15, No. 1, May 1974, p~. 55-80. L6. H. Lomax, "Some Prospects for the Future of Computational Fluid Dynamics", Proceedings of the 5th AIAA Computational Fluid Dynamics Conference, June 1981, Palo Alto, California. L7. H. Lomax et al., Private Communications, 1981-1982. 147 Ml. G. Moretti, "A Physical Approach to the Numerical Treatment of Boundaries in Gas Dynamics", Procee9ings of Numerical Boundary Condition Procedures Symposium, October 19-20, 1981 Ames Research Center, NASA, Moffet Field, California. M2. K.W. Morton and J.W. Barrett, "Optimal Finite Element Methods for Diffusion-Convection Problems", in Boundary and Interior Layers-Computational and Asymptotic Methods, J.J.H. Miller (ed.) Boole Press, Dublin, pp. 134148, 1980. M3. K.W. Morton and A.K. Parrott, "Generalized Galerkin Methods for FirstOrder Hyperbolic Equations", Journal of Computational Physics, Vol. 36, No. 2, 1980. Nl. s. Nakazawa, Ph.D. Thesis, University of Swansea, to appear. Rl. W.H. Raymond and A. Garder, "Selective Damping in a Galerkin Method for Solving Wave Problems with Variable Grids", Monthly Weather Review, Vol. 104, pp. 1583-1590, 1976. R2. R.D. Richtmyer and K.W. Morton, Difference Methods for Initial-Value Problems, Interscience Publishers, New York, 1967. Sl. J.O. Steger, Private Communications, 1981-1982. S2. J.L. Steger and R.F. Warming, "Flux Vector Splitting of the Inviscid Gasdynamics Equations with Application to Finite Difference Methods", Journal of Computational Physics, Vol. 40, No. 2, April 1981. S3. J.L. Steger, "A Preliminary Study of Relaxation Methods for the Inviscid Conservative Gasdynamics Equations Using Flux Splitting", NASA Contractor Report 3415, March 1981. Tl. E. Turkel, "Symmetrization of the Fluid Dynamic Matrices with Applications", Mathematics of Computation, Vol. 27, No. 124, October 1973, pp. 729-736. Wl. L.B. Wahlbin, "A Dissipative Galerkin Method for the Numerical Solution of First Order Hyperbolic Equations", pp. 147-169 in Mathematical Aspects of Finite Elements in Partial Differential Equations, Carl de Boor, (ed.}, Academic Press, New York, 1974. W2. R.F. Warming, R.M. Beam and B.J. Hyett, "Diagonalization and Simultaneous Symmetrization of the Gas-Dynamic Matrices", Mathematics of Computation, Vol. 29, No. 132, October 1975, pp. 1037-1045. W3. R.F. Warming and R.M. Beam, "On the Costruction and Application of Implicit Factored Schemes for Conservation Laws", SIAM-AMS Proceedings, Vol. 11, 1978. W4. H.F. Weinberger, A First Course in Partial Differential Equations, Ginn and Company, Waltham, Massachusetts, 1965. 148 Yl. H.C. Yee, "Numerical Approximation of Boundary Conditions with Applications to Inviscid Equations of Gas Dynamics", NASA TM 81265, March 1981. Zl. o.c. Zienkiewicz, The Finite Element Method, McGraw-Hill, London, 1977.