Nonlinear hydrodynamic shock propagation analysis by the finite element method

Ial. 1. &?ng sci. Vol. IS. pp. 351-356 @ Persamon Press Ltd.. 1980. Printed in Greal Britain

NONLINEAR HYDRODYNAMIC SHOCK PROPAGATION ANALYSIS BY THE FINITE ELEMENT METHOD JOHN E. JACKSON apartment of Mechanical En~nee~ng, Clemson U~v~rsity, Clemson, SC 29631,U.S.A. and THOMAS L. COST Department of Aerospace Engineering, Mechanical Engineering and Engineering Mechanics, The University of Alabama, Tuscaloosa, AL 35486,U.S.A. Abstract-The method of weighted residuals is usedin injunction with the finite element method to solve the nonlinear hydrodynamic field equations. The algebraic equations which result from the finite element approximations were solved using the Standard Newton-Raphson method, the Modified Newton-Raphson method, the Sell-Correcting Incremental Method, and the Method of Successive Substitutions. Also, the Gear 2-Step and 3-Step temporal operators were employed along with the Park temporal operator and a first order backward difference method to obtain the transient response. Solutions for l-dimensional shook and rarefaction waves obtained using the method agree with theoreticalJy predicted values and with values obtained using a finite difference method. Of the four temporal operators investigated, none was determined to be superior. However, based upon compu~tional cost and solution accuracy, the Self-Correcting Incremental Method and the Modified Newton-Raphson Method with Jacobian reformation after every three iterations both proved superior to the other methods studied for solution of the nonlinear equations.

INTRODUCTION UNTIL RECENTLY, numerical

formulation of fItrid dynamics problems has been fi~te-d~erence oriented. However, the finite element method has now begun to be applied to fluid flow problems through the use of variational methods and the method of weighted residuals. The finite element method offers a simpler means of handling irregular boundaries than do standard finite difference techniques and suggests a straightforward extension to problems involving combined fluid and structural systems. The occurrence of shocks in a fluid causes severe numerical difficulties. Since a shock embodies, in general, a discontinuity of fluid density, velocity, pressure and entropy, the partial differential equations do not apply across the shock. However, in a numerical solution procedure, such as a finite difference or finite element method, the shock discontinuity is replaced with a wave front of finite, though small, thickness over which the fluid parameters change rapidly but are stilI continuous. To smear out the shock discontinuity, Von Neumann and Richtmyer[l] added an “artificial viscosity” term to the finite difference formulation of the hydrodynamic equations. Taylor et al.[2] proposed a finite element formulation for an ideal fluid in which wave motion was to be obtained by adding to the system of equations a similar artificial viscosity term. Lax[3,4] and Lax and Wendroff 151 demonstrated a finite difference technique for handling fluid shock propagation problems which required no explicit in~oduction of additions terms into the differential equations describing the fluid system. The fundamental concept of the LaxWendroff scheme was to express all equations in conservation law form and then to solve these stepwise in time using a Taylor series expansion. Chung[6] proposed solution of fluid flow problems by a combination of the finite element technique and the Lax-Wendroff method. Chung and Chioul71 have demonstrated the use of the finite element technique for spatial discret~ation in conjunction with an implicit temporal operator in solving a boundary value problem of a perfect gas behind a normal shock wave. In this study, the fluid field equations are obtained from conservation principles and then specialized for the case of a perfect gas. Next, spatial discretization achieved by means of finite element techniques and Galerkin’s method gives a system of nonlinear algebraic equations. The numerical sohrtion methods are demonstrated by several I-Dimensions examples. Solutions given by several implicit temporal operators are compared. An evaluation of costs and accuracy 351

JOHN E.JACKSONandTHOMASL.COST

352

for various nonlinear equation solution algorithms is performed. ,A discussion of results follows, with comments on application of these methods to multidimensional problems and to fluidstructure systems. GOVERNINGEQUATIONS

The governing fluid equations may be obtained from the principles of conservation of mass, conservation of momentum, and balance of energy. From Fungi81 these are

(2) and

II $ f 5 (EQ)) I, dV =

v

((uijS1.j- qi,i)d V,

(3)

J

where E is total energy per unit volume, p is mass density and the ui are velocity components. The Xi are the body force ~om~nents, q&irepresents heat input and oij is the stress tensor. Neglecting body forces, viscosity, and heat input terms, these laws may be reduced to PJ +

m/j = 0

(4) .=O,

mi,f+(mjUj)j+(v-l)(E-~mjUj)

(5)

.I

E,t + (EUj),j + (7 -

l)(EUj -$mtuIui),j =0,

(6)

and Mi-pUi =O.

(7)

Equations (4-6) are the continuity, momentum, and energy equations, respectively. The mj are components of linear momentum, y is the ratio of specific heats and the perfect gas equation of state

has been assumed. FINITEELEMENTFORMULATION

Numerical solution of the nonlinear hydr~yn~ic domain in the manner

u; = C,u”

equations begins by discretizing the time

+ c2,

(9)

where n implies time t = t, and the forms of Cr and CZ depend upon the particuliar temporal operator chosen. Next, the finite element method requires that continuous variables be approximated by expressions of the form N

u” =

z hru;,

where N is the number of nodal points per element.

(10)

Nonlinear hydrodynamic shock propagation analysis by the finite element method

353

To solve the nonlinear equations resulting from the above process, linearize by assuming

u!+’ = l,lk + AUi

(11)

and dropping higher order terms. Application of Galerkin’s method to the linearized equations gives

I ~Y47+‘,A4i)hndV, V

(12)

where $i represents the degrees of freedom at each nodal point i. After combining the elemental results of eqn (12) the global finite element equations may be written

Here, k is an iteration cycle and J$ is the Jacobian evaluated at the kth iteration. Equation (13) is solved at each successive time step. Equation (13) lends itself to several solution schemes for sets of nonlinear equations. The Standard Newton-Raphson method requires reformulation of the Jacobian after each iteration cycle, whereas the Modified Newton-Raphson method requires reformulation only after several iteration cycles. Iteration continues until the delta quantities vanish. If only one iteration per time step is performed, the solution is a Self-Correcting Incremental Method. Iteration on the right-hand side of eqn (13) without reformulation of the Jacobian results in what has been termed the Method of Successive Substitutions[9, lo]. EXAMPLE PROBLEMS

The finite element formulation of the nonlinear hydrodynamic equations was applied in the solution of several l-dimensional fluid shock propagation problems. The numerical solutions were initiated by specifying initial values of particle velocity, pressure and mass density at each system nodal point. Data was chosen to enable comparison with the finite difference results of Lax[3]. Linear finite elements were used. Figure 1 illustrates the pressure waves obtained for a 99-element system. Results given by application of four implicit temporal operators are compared with those of Lax and with the theoretical shock wave. The four temporal operators investigated are the Gear 2-Step and 3-Step operators, the method of Park and a first-order backward difference scheme [ 111. Figure 2 compares the shock waves obtained from solution of a 29-element system by

Ax= .DC5&=.DCMl4 LI GEAR 3-STEP 0 FIRST ORDER SACKWARD 0 GEAR Z-STEP

IO

THEORETKXL

SOLUTION

I 46

I

c

O$?O

I

.&I .72 -so .56 LOCATION,x Fig. 1. Comparison of temporal operators-shock wave after 49 time increments.

JOHN E.JACKSONand THOMAS L.COST

354

2

0 STANDARD

NEWTON-RAPHSON

D MODIFIED

NEWTON-

0

t

SELF- CORRECTING

-THEORETICAL

METHOD

RAPHSCN INCREMENTAL

MET

SOCUTICN

LOCATION, x

Fig. 2. Comparison of nonlinear equation solution algorithms-shock

wave after 30 time increments.

means of various nonlinear equation solution algorithms. Solutions resulting from utilization of the Standard and Modified Newton-Raphson Methods and the Self-Correcting Incremental Method are compared for accuracy in Fig. 2 and for computational costs in Table 1. Finally, a finite element solution is obtained for the case of a shock propagating through a fluid medium and reflecting from a rigid wall. Results are shown in Fig. 3. DISCUSSIONOFRESULTS

The numerical solution of the l-dimensional hydrodynamic equations did not converge until the equations were cast in conservation law form. This result is analogous to the findings of Harlow and Amsden[ 121for the finite difference method. Before convergence could be obtained, it was necessary to eliminate all product terms from expressions involving time derivatives. Note that the form of eqns (4-7) meets both this requirement and the requirement of conservation law form. It is evident from the sample problem illustrated by Fig. 1 that the finite element solutions produce good approximations of the shape and position of the theoretical shock front and also compare favorably with Lax’s finite difference formulation. No clear advantage for any of the temporal operators was established.

50-

J=I

pp. .333

U! = 4

b=O

PI =a

Pr=o

---THEORETICAL PRESSURE

“=

NODE

TIME STEP

NUMBER

Fig. 3. Shock wave reflecting from a rigid wall.

Noulinear hydrodynamic shock propagation analysis by the finite element method

355

Table 1. Relative costs of nonlinear equation solution algorithms METHOD

CPU

SECONDS

NEWTON-RAF?iSON

69.3

MODIFIED

Nab-RA~SON,

566

JACOBIAN

REFORMULATED

STANWRD

EVERY

3RD

MODIFIED

ITERATION NEWTON-

JACOBIAN EVERY

6TH

SIXXESSIVE *29-ELEMENT tCOC

RAPHSON,

63 5

REFORMULATE0 ITERATION

SELF-CORRECTING

tN0

‘+

CYBER

INCREMENTAL RESTITUTIONS SYSTEM. 175

t 50

TIME

26 6 --STEP3

COMPUTER

CONVERGENCE

The Standard and Modified Newton-Raphson algorithms yield identical pressure waves. However, the Self-Correcting Incremental Method tends to drift from the true solution as the analysis progresses in time. Note from Table 1 that among the Newton-Raphson procedures investigated, the method requiring reformulation of the Jacobian every third iteration is most efficient. Thus, there would seem to be little purpose in using either the Standard NewtonRaphson Method or methods requiring numerous iterations between Jacobian reformu~tions. The Method of Successive Substitutions possessed serious convergence problems, as has been shown to be the case in other highly nonlinear problems[9]. The Self-Correcting Incremental Method was significantly less expensive than the Newton-Raphson Methods. For the nonlinear equation solution algorithms investigated, the best choice would appear to be either the Self-Correcting Incremental Method or a Modified Newto~Raphson Method with Jacobian reformulation after only a small number of iterations. Choice between these two methods for a given problem will depend upon the relative importance of accuracy and cost. Since the effect of the above procedures on solution costs will vary somewhat, depending upon the simultaneous equation solution algorithms utilized and the specific problem definition, the cost comparisons in Table 1 should not be taken as exact ratios but rather as indicative of general trends. Data for Table 1 was generated using Potter’s method for solution of pi-boded simultaneous equations [ 131. Figure 3 illustrates the ability of the numerical procedures to handle strong shocks. Note the ratio of reflected to incoming pressure of five, whereas for an acoustic wave this ratio would be two1141. The value for the reflected pressure, 40, is precisely the magnitude predicted theoretically. ~thou~ the i~ustrative examples are all l-Dimensions systems, the demonstrated solution procedures are identical for 2dimensional and 3-dimensional problems. Figure 3 suggests an extension of the procedures discussed herein to the solution of systems involving fluid-structure interaction. If, instead of a rigid wall, there is a structural system at the fluid boundary, the interaction problem may be formulated by finite eiement discretization of the structural system and combination of the resulting simult~eous equations with those of the fluid system. Theoretica~y, finite element discreti~tion of both systems would allow solution of problems possessing complicated interface geometries, SUMMARY ANDCONCLUSIONS

The nonlinear hydrodynamic equations for a perfect gas have been derived and then solved nume~~ly by the finite element method. It was found that for convergence reasons these equations shoufd be cast in conservation law form and should contain no product terms in the time derivative expressions. Several l-dimensional examples illustrate the ability of the numerical procedures to solve shock propagation problems. Of the four implicit temporal operators investigated, none was established as superior. However, an investigation of nonlinear equation solution algorithms indicates that of the methods studied, the choice considering both cost and solution accuracy lies between the Self-Correcting Incremental procedure and a Modified Newton-Raphson Method with Jacobian III23 Vol.

18. No. 2-H

356

JOHN E. JACKSON and THOMAS L. COST

reformulation each few, say three, iterations. The procedures developed are applicable to multidimensional systems and to problems involving fluid-structure interaction. The advantage of the finite element formulation over finite difference procedures lies in the ease of handling complex boundary geometries. REFERENCES [l] J. VON NEUMANN and R. R. RICHTMYER, J. Appl. Phys. 21,232 (1950). [2] R. L. TAYLOR,G. DASQUPTA and T. B. KHALiL, FiniteElementMethods in Row Problems (Edited by J. T. Oden, 0. C. Zienkiewic~, R. H. Gallagher and C. Taylor), 305-307.UAH Press, Unive~ity of Alabama in Huntsville, Alabama (1974). 131P. D. LAX, CO~~KRS&ire Appt. Math. 7, I59 (1954). [4] P. D. LAX, CornmansPure Appt. Math. IO, 537 (19%). [5’] P. D. LAX and B. WENDROFF, Communspure uppi. Math. 13,217 (1960). [6] T. J. CHUNG, Finite Element Analysis in Fluid Dynamics. University of Alabama in Huntsville, Huntsville, Alabama (1974). 171T. J. CHUNG and J. N. CHIOU, Analysis of Unsteady Compressible Boundary Layer Flow Via finite Elements. University of Alabama in Huntsville, Huntsville, Alabama (1974). [El Y. C. FUNG, Foundations of Solid Mechanics. Prentice-Hall, Englewood Cliffs, New Jersey (1965). [9] F. A. BROGAN and B. 0. ALMROTH, AIAA J. 9(12), 2321(1971). [lo] R. KAO, Comput. Struct. 4, 1091(1974). [ 111K. C. PARK, Comput. Struct. 7,343 (1977). {I21 F. H. H~LOW and A. A. AMSDEN, ~ajd ~namics: An ~ntrodaeto~ Text. Los tiamos Scientihc Laboratory Report LA-4100(1970). [13] Y. TENE, M. EPSTEIN and I. SHEINMAN, Cornput. Sfrucf. 4, 1099(1974). 1141S. GLASSTONE (editor), The Effects of Nuclear Weapons. LJSAEC,U.S. Government Printing OtTice(l%2). [ISIJ. E. JACKSON, Analysis of Fluid-Structure Interaction Under Hydrodynamic Shock Conditions by the Finite Element Method. Doctoral Dissertation, University of Alabama, Tuscaloosa, Alabama (1977).

APPENDIX The implicit temporal operators applied in the above study may be defined as follows1111. 1. Gear 3-Step operator

2. Gear 2-Step operator 3 “+l__*n+1+_(_4u”+U”-.‘)

il.1 -

1 2At

2At

3. Park operator

1 n+, _ 3 uJ -zG u”+’ +-661 { - 15~”+&“-’ - un-*). 4. First-Order Backward Difference Operator I

U;+I=tU

“+I__

1 At u I)