Adaptive node movement in finite-difference solutions of partial differential equations with applications to gasdynamics. Part 2

7

ADAPTIVE NODE MOVEMENT

IN FINITE-DIFFERENCE SOLUTIONS

OF PARTlAL DIFFERENTlAL EQUATIONS

WITH APPLICATIONS TO GASDYNAMICS

February, 1990

(Part 2)

by

Donaid Frederick Hawken

UTIAS Report No. 333

CN ISSN 0082-5255

rI\l.

'2..

..

ADAPTIVE NODE MOVEMENT

IN FINITE-DIFFERENCE SOLUTIONS

OF PARTlAL DIFFERENTIAL EQUATIONS

WITH APPLICATIONS TO GASDYNAMICS

(Part 2)

by

Donaid Frederick Hawken

Subrnitted December,1989

February,1990

© Donaid Frederick Hawken 1990

UTIAS Report No. 333

CN ISSN 0082-5255

ACKNOWLEDGEMENTS

I wish to thank Professors J. J. Gottlieb and J. S. Dansen for their encouragement in the writing of this thesis. My gratitude also goes to my wife. Jill. for all her help in proofreading the various drafts.

Financial assistance fr om the Natural Sciences and Engineering Research Council of Canada was very helpful and is gratefully acknowledged.

SUMMARY

This thesis presents an in-depth study of two adaptive numerical methods of solving partial differential equations. A review of the adaptive-method literature is included and has been published in the open literature. The methods reviewed involve the adaptive movement of nodes. so as to obtain a low level of solution truncation error while miniDlizing the number of nodes used in the calculation. Such methods are applicable to the solution of nonstationary flow problems that contain moving regions of rapid change in the flow variables. surrounded by regions of relatively smooth variation. Flows with shock waves. contact surfaces. slip streams. phase-change

interfaces. and boundary layers can be modelled with great precision by these methods. Significant economies of execution can be attained if nodes are moved so that they remain concentrated in regions of rapid variation of the flow variables.

One of the reviewed methods. the ~Ioving Finite Element (MFE) method originally introduced by Milier and Milier. has been implemented and is compared to a second method. called the Moving Finite Difference (MFD) method. which was developed for this thesis. Both methods solve

simultaneously for the flow variables and node locations at each time-step. These adaptive methods move the nodes so as to minimize an 'error' measure

that contains a function of the time derivatives of the solution. The error measure is manipulated to obtain a matrix equation for node

velocities. Both methods make use of penalty functions to prevent node crossing. The penalty functions result in extra terms in the matrix equation that promote node repuls ion and that become large when node separation beeomes small.

Extensive work applying the MFE and MFD methods to one-dimensional gasdynamics problems has been performed. The test problems consist of Burgers' equation. ideal viscid planar flow of gas within a shock-tube. analogous spherically symmetrie flow of gas resulting from an implosion/ explosion proeess. and propagation of a shock wave or rarefaction wave in a duet with an area change. The MFD method typically executes with at least twice the speed of the ~WE method. and the solution exhibits eomparable or greater quality. A stability analysis of the solution of Burgers' equation was performed. The results of the analysis imply an increase in stability if adaptive methods compute node-velocities that keep the nodes moving in step with the regions of rapid change in the flow variables.

...

CONfENTS Page ACKNOWLEDGEMENTS _ii SUMMARY

. .

.

. .

. .

.

_iii TABLE OF CONfENTS _iv

LIST OF FREQUENTLY USED SYMBOLS _ix

LIST OF FREQUENTLY USED ABBREVIATIONS _x

1.0 INTRODUCTION . . . ₁

2.0 REVIEW OF ADAPTIVE METHODS . ₃

2.1 Introduction ₃

2.2 Overview of Adaptive Methods ₅

2.3 Summary of Var ious Adapt i ve }{ethods • ₈ 2.3.1 Methods Based on Integral Minimization

or Equidistribution of Error Measure . 8 2.3.2 Methods Based on Attraction and Repulsion

Pseudoforces Between Nodes • . . . . . 30 2.4 C10sing Discussion of the Adaptive Methods

Surveyed . . . . . . . . . . 41 3.0 OPTIMAL SOLUTION OF STIFF PARTlAL DIFFERENTlAL

EQUATIONS • . • • • . . . . • . • . . • • . • • ₄₅ 3.1 Imp1icit Versus Explicit Methods App1ied to Stiff

Prob1ems • • . • • • • . • 45

3.2 Characterization of StUf Ordinary Differentia1 Equation Solvers

. · · ·

·

·

.

.

46 3.3 Implic i t Runge-Kutta Methods

. .

. .

_·

_·

.

.

₄₉ 3.4 Imp1icit Mu1ti-Step Backward-Difference Method

of Gear

.

.

.

.

_{· ·}

_{· . ·}

_·

·

· ·

50 3.5 Comparison of the Runge-Kutta and Gear

Backward-Difference Methods

_·

51

3.6 Modified Newton Iteration for the Gear Method

and lts Descendants

_{· · ·}

. · · ·

₅₂

3.7 _{Details of Imp1ementing the Method of Gear for} Adaptive Methods

_·

.

. . . . ·

·

· ·

55 iv

3.7.1 The P + 1 Value Normal Form of the Gear Equations

3.7.2 Time-Step Size and Order-Control Formulae 3.7.3 The Variable-Coefficient Gear Method • • • • 3.7.4 The Fixed-Leading-Coefficient Gear Method 4.0 THE MODIFIED MOVING FINITE-ELEMENT METHOD

4.1 Calculation of Left-Hand Side Terms . 4.2 Calculation of Right-Hand Side Terms

4.2.1 Terms Containing a Source Function •

4.2.2 Terms Containing First-Order Spatial Derivatives • • • • • . . • • 4.2.3 Terms Containing a Second-Order Spatial

Derivative . • . • • . . . • • • 4.3 Symbolic Input of Partial Differential Equations 4.4 Means of Improving MFE Performance

4.4.1 Overemphasis of Viscosity 4.4.2 Horizontal Emphasis

4.5 Analytical Estimates of the Reduced Finite-Element Integrals . . • . • • • • • • •

4.5.1 Integral Containing a Constant •

4.5.2 Integrals Containing One Component of the Solution . . • . • • • • • 4.5.3 Integrals Containing Two Components of

the Solution • . . • • • . • . • 4.5.4 Integrals Containing Three Components

of the Solution • • • • • 4.5.5 Removal of Singularities from the

Integration Formulae • • • • • • • • . 5.0 THE MOVING FINlTE-DIFFERENCE METHOD

5.1 Use of an Error Measure Based on Time Derivatives

56 57 59 61 64 68 70 71 72 74 78 79 79 80 80 81 81 82 86 89 93 of the Solution • • . • • • • . 94 5.2 Failure of

MFE

Method and Simple MFD Method

on Stationary Problems 100

5.3 Improvement of the Method by Use of a

Solution-Based Error Measure • • • • . • . . • . • . 100 5.4 Adaptive Method Based on Minimizing the Integral

of an Error Measure • • • . • • • . • . . • 105 5.5 The Ultimate Form of the MFD Equation for

Node Velocities • ₁₁₃

5.6 Calculation of Right-Hand Sides of PDE's . . ₁₁₄ 5.6.1 Terms Containing a Source Function • ₁₁₅ 5.6.2 Terms Containing First-Order Spatial

Derivatives . • • . . . . • • 115 5.6.3 Terms Containing a Second-Order Spatial

Derivative . . • • • . . • . . 116 5.6.4 Rejected Differencing Methods

_·

_·

·

·

117

6.0 NODE-MOVEMENT PENALTY FUNCTIONS

.

.

.

.

₁₁₉

6.1 Old Penalty Functions

.

_{· ·}

·

·

120 6.2 New Penalty Functions

. . .

·

· · ·

121 6.3 Gradient-Weighted New Penalty Functions

.

. ·

_·

_{· ·}

122 6.4 Amplitude-Change-Weighted New Penalty Functions 123 6.5 Penalty Functions Derived Using the Calculus

of Variations . . . • • . • . • . . 124 7.0 DERIVATION OF ONE-DIMENSIONAL TEST PDE'S ₁₂₆ 8.0 BURGERS' EQUATION TEST PROBLEM •

.

.

.

.

. . . .

₁₃₁

8.1 Boundary Conditions • . ₁₃₁

8.2 Stability Analysis with Stationary Norles ₁₃₂ 8.3 Simple Stability Analysis with Moving Nodes . ₁₃₈ 8.4 Effect of an Adaptive Equation on Stability

Analysis ₁₄₀

8.5 Effect of the MFD Method on Stability • ₁₄₃

8.6 Results of MFD Calculations ₁₄₄

8.6.1 Effect of Variation of Right-Hand-Side

Parameters • •

. . .

₁₄₅

8.6.2 The Effects of Varying Other Parameters 145

8.6.3 Burgers' Equation in Conservation-Law Form. 146

8.7 Results of MFE Calculations . 146

9.0 SHOCK-TUBE PROBLEM • 148

9.1 Contact-Surface Profile and Thickness 149 9.2 Shock-Wave Profile and Thickness 150

9.3 Initial-Condition Profiles 152

9.4 Initial Placement of Nodes 152

9.5 Boundary Conditions . . • • • 153

9.6 Results of MFD Calculations • 155

9.7 Results of MFE Calculations • 159

10.0 SPHERICAL EXPLOSION/I~WLOSION PROBLEM . 161 10.1 Description of the Explosion/Implosion Process 161

10.2 Results of MFD Calculations 162

10.3 Results of MFE Calculations 168

11 .0 SHOCK WAVE IN A DUCT WITH AN AREA ENLARGEMENl' 169 11.1 Physics of Interaction of Shocks and Area

Enlargements • • . • • • • • •• 170

11.2 Initial Conditions 171

11.3 Description of RCM Solution 173

11.4 Results of MFD Calculations 174

12.0 RAREFACTION WAVE IN A DUCT WITH AN AREA REDUCTION . 176 12.1 Physics of Interaction of Rarefaction Waves

with Area Reduct ions . • • . • • • • 177 12.2 Initial-Conditions .

12.3 Description of RCM Results • 12.4 Results of MFD Calculations 12.5 Results of MFE Calculations 13.0 FINAL DISCUSSION AND CONCLUSION • •

13.1 Comparison of the Performances of the NFO and MFE Methods • • . • • • • • . • • • •

178 178 179 180 181 182

13.2 Recommendations for Further Work . 13.3 Contributions

REFERENCES • • • • FIGURES AND TABLES .

. . .

viii 184 185 187 199

a Pe Pr Re t At

x

T y p

LIST OF FREQUENfLY USED SnmOLS

speed of sound

undisturbed speed of sound

vector having the principle gasdynamical variables of the problem to be solved as components

component m of

A

value of component m of

A

at node j

~

_

~-l

j-l -~

value of

A

at node j

dot product of

AAj

with itself Pec let number

Prandtl number Reynolds number time

size of time step distance

spatial coordinate at node j

xj+l - Xj-l temperature

undisturbed temperature ratio of specific heats density

undisturbed density

PDE

ODE

~E

mro

CFL

LIST OF FREQUENTLY USED ABBREVIATIONS

Partial Differential Equation

Ordinary Differential Equation

Moving Finite Element

Moving Finite Difference

t

6.0 NODE-MOVEMENT PENALTY FUNCTIONS

These were introduced by Miller and Miller [64,65], and further developed in Gelinas, Doss, and Milier [67], Djomehri [73], and Milier [74] as a means of preventing nodes from crossing each other in regions of large solution variation and also to prevent singularity of the matrices

associated with the Aloving Finite Element (MFE) method first introduced by Milier and Milier. The penalty functions have been adopted, in modified

form, for use in the Moving Finite Difference (MFD) method described. in section

S.O.

In both the MFE and the lUFD methods there is a tendency for the nodes to concentrate excessively in the interiors of shock waves since the nodes tend to follow the characteristics, which congregate in a shock

(and cross for inviscid problems). Extra terms are inserted into the

MFE

and AWD equations used for calculating the node velocities. These terms become dominant when the nodes are very close together and in addition prevent degeneracy of the system of equations in regions of zero (or, in

the case of the lfFE method, constant) slope. The extra terms are the re8ult of Dlinimizing the following measure of relative node velocity and node position,

N

I

2

[

p ei[ii - ii-l] si]2 (6.1)

i = 1 by requiring that al p 0 aij (6.2)

for all nodes j. In this measure, ei and si are positive functions ofAX i and are defined in such a way that they become very large when AX i is very smallor approaches some specified Dlinimum value. Under such conditions the change in node velocity from node to node will approach a value dictated by the ratio of si to ei.

The result of the minimization,

= (6.3)

is added to the AfFE or

MFD

equation for calculating the velocity of each node j. Consideration of cases where the terms in equation (3) are dominant will clarify the operation of the penalty functions. IfAXj is excessively small, then the terms containing ej and sj dominate the equations for calculating the velocity of node j - 1 and the velocity of node j. One obtains

.j X [SECTION 6.0] . j-l

X

+ (6.4) 119

thus resulting in a tendency for the nodes to separate again. In the case where AXj - AXj+1 and where both are excessively small. one obtains

' j

X

2

(6.5)

In other words. the velocity of node j tends to become the average of the veloeities of the two adjacent nodes plus an additional amount that reduces motion towards the closest node. It should be noted that there is no

limitation on the absolute mot ion of the nodes. only on their relative mot ion. When the penalty function terms dominate. there is a tendency for the nodes to become equidistributed. Hence. for example. in the interior of a shock the nodes are stringently equispaced. It is desirabie to extend the region of equidistribution to include the entire shock so that node separation will change significantly only in regions of small slope. Variation of node separation in regions of large solution variation can

lead to increased discretization error. In the ideal case. a cloud of equidistributed nodes will follow the moving shock. In practice. using current forms of the penalty functions. node separation does start to

increase near the top and bottom of the shock where the solution has a corner. This can lead to flattening of the corners (MFE) or the

development of small oscillations or kinks in the corners (MFD). Judicious selection of penalty function parameters can reduce th is effect.

Note. the geometrical factor aj is not included in the penalty functions so that they will be large enough to prevent node crossing near the origin in problems with cylindrical and spherical symmetry. This is the reason why the common factor of aj is retained in equátion (5.50). the

MFD equation for node velocities.

The actual forms of ej and sj have varied considerably. The most usefuI forms are discussed below.

6.1 OId Penalty Functions

The 'oId' node-movement penalty funct ions were first introduced by )filler [65] and are discussed here for historical purposes. They take the forms

(6.6)

and

+ (6.7)

where Co is the specified m1nlmum node separation and is usually set to several times smaller than the expected shock width. The value of C1/C3 is chosen to be several times smaller than the expected shock velocity so as

•

,

to limit the relative velocity of repuls ion of adjacent closely spaced nodes within a shock. The value of C3/C4 is chosen to be several times larger than the expected shock width so as to cause a sharp reduction in the ratio of sj to ej when the node separation approaches and exceeds the expected shock width, thus reducing the repulsion of adjacent nodes. Another guideline in the selection of C4 is that the square of C4 should not dominate the other coefficients of node velocity in regions away from the shock.

The value of C2 may be specified so as to provide a small velocity of attraction between adjacent nodes when the square of the node separation exceeds Cl/C2' thus preventing excessive node separation, but in practice C2 is usually set to zero. When C2 is given a nonzero value, the value is selected to make the product of C2/C4 and the average node separation somewhat larger than the expected shock velocity. Therefore, the nodes in advance of a shock can equally space themselves and withdraw faster than the advance of the shock.

Cl is adjusted to make the terms eJeJ and ejsj large enough to prevent the crossing of nodes within a shock but small enough to allow proper adaptivity. The value of Cl can of ten be increased above the minimum value required to prevent node crossing. causing the nunmer of nodes actually within a shock to be reduced. Reduction of the number of nodes within a shock of ten results in much larger time-steps even if an implicit ODE solution method is being employed. An optimal value for Cl is most easily arrived at by trial and error. The values of C2' C3' and C4 are set in relation to the value of Cl in accordance to the rules given above.

6.2 New Penalty Functions

The 'old' node-movement penalty functions are fairly complicated to specify. A simpler set of penalty functions was introduced in Djomehri

[73] and Miller [74] for use with the MFE method.

A PDE with second-order spatial derivatives will tend to be more stable than a similar PDE without these diffusion effects. This is because the presence of viscosity broadens shock waves and opposes the tendency of characteristics to cross within the shocks. Examination of equation (4.49) reveals that the diffusion effects resulting from eval~ati~n of (~j.~Axx) show up as terms that are roughly proportional to ~[AAJ/AXJ]2, where ~ is a diffusion coefficient and Axx is the second-order spatial derivative of the solution. These terms tend to prevent node crossing and hence have an effect complementary to that of the ejsj penalty terms. Thus. a reasonable functional form for ejsj introduced by Djomehri [73] is

(6.8)

where C4 = 2 and C2 is chosen to be 20 or 30 times smaller than ~ times the square of the allowable spatial truncation error. Since eJsJ varies with respect to node separation in the same manner as typical diffusion terms when node separation is not close in value to CO, the diffusion terms will

dominate ejsj unless AAj is somewhat smaller than the allowable spatial truncation error. If node separation approaches the minimum allowed value CO, then the terms in ejsj become dominant and tend to discourage node crossing.

Examination of equation (4.28) reveals that evaluation of (~j.~j)

results in MFE matrix terms proportional to [AAj]2/AXj. Thus. areasonabie

form for ejej introduced by Djomehri [73] is

(6.9)

where

Cs

= C4 - 1 and where C3 is chosen to be 20 or 30 times ~a~ger than the square of the allowable spatial truncation error. Since eJeJ varies with respect to node separation in the same manner as the (~j.~j) terms when node ~e~aration is ,not close in value to CO, the (~j.~j) terms will dominate eJeJ unless AAJ is reduced to the order of magnitude of the

allowable spatial truncation error. The value of C3 may be made larger in order to make the nodes move more smoothly together. In some cases. Co can be made very smallor zero without harmful effect.

In practice. for a constant value of C3' increasing C2 will increase

the tendency of the nodes to spread out and. when carried to an extreme.

will remove them fr om areas where they are most needed. Too small a value of C2 may result in failure to drive enough nodes out of the interior of shock waves so as to properly resolve the corners in the solution at the top and bottom of the shocks. The use of nonzero overemphasis of viscosity

reduces sensitivity to the value of C2.

On

average. t~e velocity of

repulsion of two adjacent nodes will vary as C2/[C3[AXJ - CO]].

Given a fixed ratio of C2 to C3' if C3 is increased there will be an

increased tendency towards equidistribution of nodes which may reduce beneficially the spatial separation of nodes near the top and bottom of

shock waves. However. when such equidistribution is carried to an extreme. it also increases the spatial separation of nodes unduly in the interior of the shock waves. A balance must be found by trial and error.

The same penalty functions can be used in

MFD

calculations. if C4 is set to 3. to match the variation of the diffusive

MFD

terms with respect to

node separation in the tridiagonal form of the

MFD

method. Similar values of C2 and C3 can be used for some problems whether they are solved using the

MFD

or the MFE method.

6.3 Gradient-Weighted New Penalty Functions

It is desirabie to make the penalty functions sensitive to solution properties as weIl as to node separation. Gradient weighted penalty

functions were originally introduced by Gelinas and Doss [69]. The present formulation is based on that of Djomehri [73] and may be written

where CG is a small positive constant (CG «1). Equation (3) may rewritten in the form

xj+1ej+1ej+1 + Xj-1 e jej ej+1ej+1 + ejej ej+1sj+1 - ejsj ej+1ej+1 + ejej (6.10) (6.11)

Consider the case where node J 1S at the tail of a right-running shock

wave. As.sume the following: AXj

>

AXj+1, the slope is ste~per in the direction of node j + 1 than in that of node j - 1, and AXJ is very small

so that equation

(3)

dominates the

MFE

or ~WD equation for the velocity of node j. As a consequence of the first part of equation (11), the velocity of node j will tendto be a weighted average of the velocities of the adjacent nodes with an emphasis on the velocity at node j + 1 (the closest node). The gradï"ent we ight ing increases this emphas is and at the same time

reduces the negative component of velocity introduced by the second part of equation (11). Node j will not move as much towards node j - 1 as it might have without the use of gradient weighting. In the interior of a shock where the gradients are constant and even larger, the gradient weighting

redu~es even more the magnitude of the second part of equation (11), thus tending to freeze the relative positions of the nodes unless they come extremely close together. The net effect is to stabilize the relative mot ion of the nodes in areas of large gradient. reducing the bleeding of nodes from the interior of shocks over many time-steps.

6.4 Amplitude-Change-Weighted New Penalty Functions

When the

MFE

penalty functions were applied to ~um problems, it was discovered that in s·ome cases the node-to-nod.e change in

Am

became excessive in shock interiors, even if gradient weighting was used. In orde-r to prevent this, the following modificat ion to equat ion (8) was

introduced:

(6.12)

The term ejsj is normally positive, but if IIAAjl12 (the value of the summation) exceed$ I/CA' which can occur in a shock wave if too many nodes are. driven eut, then the term reverses sign. This has the net effect of changing node repulsion into node attraction. The value of IIAAjl12 will reach an equilibrium value determined by the value of CA'

6.5 Penal1y Functions Derived Using the Calculus Qf Variations

In section 5.4 the calculus of variations was nsed to derive equation (5.39). which incorporated smoothness terms. and equation (5.42). which

included an alternative form of the smoothness terms. Equation (5.46) was

subsequently developed without including the smoothness terms. If the first form of the smoothness terms had been used. then

L [ a-I] - X + 2X - X + [

·j+l .j ·j-l

- -

[l)Xiïl]a+jJ--(6.13)

would have been added to equation (5.46). Because of the similarity of equation (5.46) and equation (5.50) (with Q' ~ 0). the smoothness terms above could be used as alternative penalty-function terms for the MFD method. However. since the signs of the coefficients of the node velocities would not always enhance the diagonal dominance of the MFD

matrix and since there is no strong repulsion of extremely close nodes. the second form of the smoothness terms wonld be more suitable.

If the second form of the smoothness terms had been used then. in the limit of small

p.

= +

+

(l'L

[Axjïa+~-T (6.14)

would have been added to equation (5.46). The left-hand side of the expression can be constructed from the new-penalty-fnnction ejej terms if

Cs

=

a +

p

=

2 and if C3

=

L[a +

fl -

1]

=

1.. The right-hand side varies in the same way as the right-hand side of equation (5.50) with respect to node

separation and can be constructed from the new-penalty-function ejsj terms if C2

=

n'L and if C4

=

a +

p -

1

=

1. In contrast to the first form of the smoothness terms. the signs of the coefficients of the node velocities always enhance the diagonal dominance of the MFD matrix. and the right-hand side provides a fairly strong inflnence against node crossing.

Eqnation (14) is a useful alternative to the new penalty functions. defined by equations (8) and (9). Note. from the above paragraph. that C4

= 1, in contrast to the new penalty fnnctions where C4 = 3 to match the variation of the diffnsive

MFD

terms with node separation. Therefore. if equation (14) is used. sj/ej is not as large as it would be if the new penalty functions were employed, and thus node repuls ion within a shock

to equation (14) but with O' = 0, and so they are even less effective in preventing node crossing. Owing to the frequent regridding used by Petzold, this is not a serious problem in her work.

7.0 DERIVATION OF ONE-DIMENSIONAL TEST PDE'S

The set of partial differential equations describing the time evolution of a calorically and thermally perfect gas, in the case where there is no internal heat generation, no long-range electrical or magnetic

interactions, and no body forces can be written in the form

= continuity equation

av

p-=

at

Navier-Stokes equations

ae

p

-at

= \lP + 'V

[[~

+ I1r]\l

.y]

3 P\7'Y +

\7 .

D ..

\7T]

+ I

equation of conservation of internal energy

(7.1)

\l.UI1'VXY]

(7.2)

(7.3)

where t is time,

\7

is the vector-valued spatial gradient operator, p is the density , Y is the gas velocity, and T is the temperature. P = pRT is the pressure and e = RT/[y-l] is the internal energy per unit mass in the gas; R is the gas constant and y is the ratio of the specific heats.

Within the second-order spatial-derivative terms, À is the coefficient of

thermal conductivity, 11 is the ordinary coefficient of viscosity, and I1r is the coefficient of (rotational) bulk viscosity. See Anderson, Tannehill, and Pletcher [15], Schreier [128), and Hughes and Gaylord [129] for a more detailed treatment. The dissipation function I represents the rate at which mechanical energy is expended into heat energy in the process of deformation of aviscid fluid, and it may be written

I = + + +

+ (7.4)

where the tij are components of the deformation rate tensor (i.e.,

derivatives of each component of velocity with respect to each coordinate direction) defined, for example, in Hughes and Gaylord [129). The

relaxation phenomena quantified by 11, I1r' and À are sources of dissipation, which acts to increase the thickness of transitions.

If a gas is disturbed, the translational degrees of freedom generally relax af ter four or five molecular collisions, and the rotational degrees of freedom af ter a few more, while vibrational degrees of freedom or

excited states of atoms may require many more collisions in order to relax. The ordinary coefficient of viscosity, 11, is a quantity that describes the effect of relaxation of mot ion in the translational degrees of freedom. If an excited internal degree of freedom has a relaxation time that is smaller than the effective time constant of wave mot ion, then a bulk viscosity can

be defined. which accurately describes the effect of the internal degree of freedom. It happens that the rotational relaxation time of diatomic gas molecules is of the same order as the shock transition time. If a

~Iaxwellian distribution of molecules is assumed. then the bulk viscosity may be approximated (see Bradley [130]) by

't -fJ.

8

(7.5)

where 't is measured in units of number of intermolecular collisions and has

a value of about 5.3 for air. Some of the test problems solved for this thesis were inspired by those solved by Honma and Glass [131] using the random-choice method. Af ter Honma and Glass. ~r will be set to 2~/3 end ~ end À are assumed to be constant in what follows since the problems to be solved involve relatively weak pressure transitions.

In the case of one-dimensional planar flow (Q

=

0). cylindrically

symmetric flow (Q = 1). or spherically symmetric flow (Q

=

2). equations

(1) through (3) can be reduced to ap 1

a

x <lp V =

--

-

(7.6) at xQ ax

av

_~v RapT +

2-

_

~[ xQ a2~V

] _

Q2~V

---at

ax p ax x<lp ax ax x 2p (7.7) and aT _~T [y - 1]T axQv + y - 1 ~[ xQ aÀT ]

at

ax xQ ax RX<lp ax ax +

2p[Y~Q

[

[::f+

Q[if

1

Rp (7.8)

If ~ and À are set to zero. a set of equations is obtained that may be solved by the method of characteristics (Moe). If equations (6) to (8) are transformed into alocal reference frame that moves with the velocity of the characteristics. then the time derivative on the left-hand side of the equations gives rise to a term equal to the characteristic velocity times the slope of the variabIe on the left-hand side. This extra term cancels out with the inviscid right-hand side of the equations so that p. V. and T are constant in the local time frame. The Moe is not applicable to the viscid form of the equations but some numerical techniques based on the Moe such as the random-choice method of Glimm [1] are able to handle the extra terms by means of operator splitting. Some adaptive methods. such as that of Kansa. Morgan. and Morris [42] or that of Ghia. Ghia. and Shin [43]. try to move the nodes along the characteristics so that the first two terms are more or less cancelled out of the

PDE's.

This reduces or removes the effect of the truncation error in these first-order spatial-derivative terms on the usually small. second-order spatial-spatial-derivative terms.

By use of simple manipulations. the set of planar equations can also be written in the form used in Honma and Glass [131].

ap

=

at apv

=

at apv _ QPV ax x apv2 + P ax a[E + P]V ax

~V2

a22JlV +

-x ax2

~[E

+ P]V

x

+ +

_-

Q

_"

_{- -}

a2JlV x ax + (7.9) (7.10) Q aÄ.T + JlV2 x ax (7.1.1)

where E

=

p[e + v2/2] is the total energy per unit volume. In the case of planar flow (Q

=

0) equations (9) to (11) are in conservation-law form. That is. eitber the coefficients of the derivative terms are constant. or,

if they are variabie. their derivatives appear nowhere in the equation. This conservation-law form is. in facto an expression. in different ia!

format. of the conservation of some quantity (i.e .• mass. momentum. energy) and of ten gives improved results in nonadaptive finite-difference or

finite-element solutions of flows containing shock waves. The equations can be put in near conservation-law form for nonplanar flows. Defining the momentum density.

M

= pV. one can also rewrite the equations so that the

right-hand sides are expressed only in terms of variables contained on the left-hand sides. Tbe results are

ap = at aM

=

at aE = at + 1 ax<4I x Q ax 1 ax<\t2/ p

-

-

_

.

_

-xQ ax

~

-

~[

xQ[ x Q ax + y - 1 aM2/p - 2E 2 ax Q

2J.lM

EM y - 1 M3

]

]

y -p 2 _p2

~

_![

x~[ Ä.~ ~

+ _{[ Jl} -xQax ax R p (7.12) (7.13) y -1

~

] ] Ä. 2R (7.14) Only equation (13) is not in conservation-law form. because of the last term. This last term contains a factor of Jl and should only be important for small values of

X

where. however. M/p = V will be smal! or zero to match typical boundary conditions. Tbe value of conservation law

form is somewhat diminished in adaptive methods since. when the equations are transformed to the frame in which the nodes are stationary. each left-hand side gives rise to an extra term equal to the product of the node velocity and the slope of the variabIe on the left-hand side. The transformed equations are no longer in conservation-law form for this reason and also because. in the case of nonplanar problems. the spatial derivatives have coefficients l/XQ that are not constant if the nodes are

moving. One means of avoiding the latter problem is to recognize that the coefficients l/XQ can be expressed as spatial derivatives of functions of X using 1 X aln [X]

ax

and 1 x2 = (7.1S)

This also has the bene fit of avoiding division by zero at the or1g1n of the nonplanar coordinate system when the system of PDE's is discretized.

The one-dimensional equations in conservative or nonconservative form will be used in section 9.0 for a shock-tube problem (Q = 0) and used in section 10.0 for a spherical explosion/implosion problem (Q

=

2). In the first case. a diaphragm between regions of differing pressure is suddenly removed at time t

=

O. In the second case. a pressurized spherical region starts to expand due to the removal of a spherical diaphragm at time t

=

O. In both cases. a shock wave moves into the region of lower pressure and a rarefaction wave moves into the region of higher pressure. A contact surface (surface separating regions with the same pressures but different pressure histories and therefore different final temperatures. densities. and energies) initially located at the position of the diaphragm moves at the gas velocity into the region of lower pressure.

Solution of other equations that can be derived-from the form of equations (6) through (14) will be discussed in sections 8.0. 11.0. and 12.0.

In order to avoid numerical difficu1ties associated with disparate magnitudes of the solution components. it is desirabIe to normalize the set of equations. This can be accomplished by replacing P by P/P1' V by V/al' T by T/T1' E by E/P1[a1]2. X by X/LN' and t by a1t/LN. Pl' al. and Tl are the density. speed of sound. and temperature in the low-pressure region before removal of the diaphragm if. for instance. a shock-tube problem is being solved. LN is some length characteristic of the gas flow. such as the length of the region of high pressure. or the length of an area change. or. if the behaviour of the gasdynamical variables within transitions is being investigated. the mean-free path between collisions of gas molecules. Note that al

=

JrRT1 because the gas is calorically and thermally perfect.

Only the normalized variables will be used in sections 8.0 through 12.0. and the PDE's will be in the corresponding normalized form. In a slight abuse of notation. the same names will be used for the normalized variables as for the non-normalized variables since there is no danger of confusion.

..

Note that the coefficients of the second-order terms in the

normalized equations will be proportional to the inverse of the Reynolds number of the flow. given in this case by

Re = (7.16)

In some problems dealing with steady or quasi-steady flow. the velocity at infinity would replace al in the normalization. leading to the more

familiar definition of the Reynolds number. but in this case al better characterizes the problems to be solved.

Recall that. for a perfect gas.

= yR /.l

y - 1 Pr

(7.17)

where Cp is the specific heat at constant pressure and Pr is the Prandtl number. which has a value close to unity. Note that Cp/Pr = Ä//.l is approximately constant for most gases.

The product of the Reynolds and Prandtl numbers. given in this case by

P1 a 1LN yR

Pe =

--

-

--

(7.18)

Ä y - 1

is of ten referred to as the Peclet number. NOde spacing must be reduced as the Reynolds and Peclet numbers increase in size. to prevent spurious

oscillations near steep portions of the numerical solution •

S.O BURGERS' EQUATION TEST PROBLEM

The equation

au

at

= ~u+

ax

2

ax

2

(S.l)

introduced by Burgers [132], models many of the features of the systems of one-dimensional gas-dynamics equations introduced in section

7.0.

The model equation has a time dependent term, a nonlinear advection term, and a diffusive or dissipative term with 1

»

6

>

O. Many of the stability considerations that apply to the more complex systems of equations are more clearly illustrated by a treatment of Burgers' equation.

Consider a finite-difference calculation utilizing stationary nodes with uniform separation AI. Tbe nodes are numbered fr om the left starting at i

=

0 and ending at i

=

N. Application of central differences to

equation (1) yields the ODE

aui

where Ui+l _ Ui-l Ui 2AI + 6 Ui+1 2 2Ui + Ui-1 [AX]2 (S.2) (S.3) at any internal node i. The size of AI must be chosen small enough that the truncation error in the advection term does not dominate the diffusive term, and hence AX must be reduced as 6 is reduced.

8.1 Boundary Conditions

A boundary condition is specified at each end of the domain. The value of

U may be specified at node zero, in which case EO is given the value of the corresponding time derivative (usually zero). In other cases reflection symmetry or antisymmetry is specified about the boundary so that

EO = 6 ~-~ (S.4)

[AX]2 or

EO = _~ UI - ~ (S.S)

A fourth alternative is to employ periodic boundary conditions, in which case

u

O _{= UN} and UO UI - UN- 1 AX

&

UI - 2UO _{+ UN-1} + 2 [AX]2

The boundary at node N is treated in a similar way.

8.2 Stability Analysis with Stationary Nodes

(8.6)

(8.7)

This section and sections 8.3 through

8.S

present a progressive development of an advanced stability analysis of adaptive methods as applied to

Burgers' equation. Equation (2) may be written in vector form

a!!

=

at

(8.8)

where

!!

and

li

are column vectors with components Ui and Ei, respectively. As detailed in section 3.2, it is the eigenvalues of the Jacobian of

li

that determine the stability of the solution obtained by a particular ODE

differencing method. Assume that the solution at each node has an error, si, resulting perhaps from inaccuracies in the initial conditions. Replace Ui by Ui +si in equation (2) and subtract the original equation from the

·newequation. Assume that si is very small and neglect its square. The result, at each internal node,

=

at 2AX Ui_{+1 -}

_u

_{i -}1 si 2AX

,

+ 2 [AX]2 (8.9)

may be thought of as an ODE for the propagation of error in the initial conditions at node i. By inspection it can be inferred that the error will grow most quickly in regions of large U or in regions of large negative slope of U.

Examination of equation (9) reveals that the coefficients of si+1, si, and si-1 are the derivatives of Ei with respect to Ui+1, Ui, and Ui-I, respectively. The coefficients thus form the i th row of [J], the Jacobian matrix of

li

with respect to~. Equation (9) can be expressed in vector

form as

(8.10)

where ~ is a column vector with components ei. [J] is easily seen to be tri-diagonale The N + 1 eigenvalues of the Jacobian are the roots of the determinant of [J] - À[I]. where [I] is the identity matrix. In general. the resultant expressions for À are too complex to be easily interpreted.

It is possible to obtain a bound on the eigenvalues of [J] by

applying the Gerschgorin Circle Theorem discussed. for example. in Burden and Faires [133]. A circle is inscribed on the complex plane for each diagonal coefficient in an n by n matrix. Each circle has a centre given by the value of a diagonal coefficient and a radius given by the sum of the absolute values of the off-diagonal coefficients in the same row. The Gerschgorin Circle Theorem states that the n eigenvalues of the n by n matrix must be contained within the regiones) enclosed by the inscribed circles. Moreover. if k of the circles do not intersect the remaining n - k circles. then precisely k eigenvalues are contained within the

regiones) enclosed by the k circles. Applying the theorem to [J]. one can show that the circles are defined by

I

Ö Ui+l - Ui-1

1

2 Z + - -- + = [AX]2 2AX

I

ö

Uil

[AX]2 +

~

(8.11)

'for each interior node. Here Z = ZRE +jZIM is a complex number with j being the square root of minus one. IfAX is sufficiently small, then ö/[AX]2 is larger than the magnitude of Ui/AX. so that the absolute value signs on the right-hand side may be removed. Simple manipulation yields

[ ZRE + _[AX]2 Ö +

Ui+l _ Ui-l ]2 2AX

+ = (8.12)

The specific manner in which the boundary conditions are handled affects ~ and hence the coefficients and eigenvalues of [J]. If the value of U is constant at a boundary, the corresponding circle has zero radius with centre Z

=

O. If reflection symmetry is specified at node zero, then the circle is defined by

(

ZRE +

- -

Ö

t

+ [ZIM]2 =

---

ö2

[AX]2 [AX]4

(8.13)

A circle defined by

[

ZRE +

- -

ö +

Ul~!JO

t

+ [ZIM]2

₌

---

ö2

[AX]2 [AX]4

(8.14)

is the result if reflection antisymmetry is specified at node zero. The other boundary is treated in a similar manner. If periodic boundary conditions are employed an expression very like that of equation (12) is obtained.

"

In other words, the eigenvalues of [1] are within the area covered by N - 1 circles of radius &2/[AX]4 having centres at

z

=

---

& Ui+1 - Ui-1 (8.15)

[AX]2 2AX

with i equal to 1 through N - I, or are within the two regions defined by the coefficients in the first or last rowof [1]. Tbe circles will cross the imaginary axis in regions of negative slope; hence, it is possible for some eigenvalue(s) of large magnitude to be close to the negative real axis unless &/[AX]2 greatly exceeds the magnitude of the slope. Existence of eigenvalues of large magnitude near the real axis can adversely affect the stability of the numerical solution. Tbus, it is important to understand the effect on the eigenvalues of [1] produced by large magnitudes of Ui/AX and by large negative values of the slope .

It is useful to examine the behaviour of

ei+1 - e i - 1

C---=

at

2AX + & ei+1 - 2ei+ 2 [AX]2 (8.16) which is closely related to equation (9). Bere C is defined to be the largest value of U in tbe domain and C' is the most negative value of the slope of U in the domain. Tbe corresponding lacobian will be designated as

[.T'].

If periodic boundary conditions are employed, then a procedure closely related to the Fourier stability analysis of PDE's, described in Anderson, Tannehill, and Pletcher US], can used to determine the

eigenvalues of [1']. Since the differential equation is linear in ei, one can easily obtain the eigenvalues by expanding ei in a complex Fourier series and considering each Fourier component separately. That is,

(8.17) where em, the amplitude of the mth Fourier component, is a function of time only. Km = nm/NAX where m = 0, I, 2 .•. N. Tbe period of the fundamenta! frequency (i.e., m = 1) is assumed to be twice the length of the domain or 2NAX. Rence ei ±l

=

_{eiexp[±jPm] where Pm} =

'mAX,

and equation (16) is transformed into exp[jPm] - exp[-jPm]

C

~m 2AX +

=

---~em [AX]2

at

2 C'e m (8.18) [SECTION 8.2] 134

for eaeh of the N + 1 values of Pm. Thus, '" m = _ö_[l -_[AX]2 cosP m] C' .C . A J-S1n"m AX (8.19)

is the mth eigenvalue of [J']. Incidentally, the eorresponding eigenvector

of [J'] is obtained by evaluating equation (17) for each component of ~.

As detailed in section 3.2, the product of At and any eigenvalue with negative real part must be within the absolute-stability region of the numerical method used to solve the system of ODE's. If the explicit Euler method were used, then the region of -absolute stability is a cirele of unit radius about "'AT equal to negative one. In other words, using equation

(19), there is a requirement that

[

-

C'At _{+ 1}

r

+

[CAt]2 _AX sin Pm 2

i

1 (8.20) If C' is set to zero, then detailed analysis yields the stability condition

I

C~:12

_{a A}

i

I

_ö[AX]2 At

I

_i

₁

This is in fact the stability condition for the linearized Burgers' equat ion

au

at

=

Cau

ax

+ ö

a

2

u

2

ax

2

,

treated. for example, in Anderson, Tannehill, and Pleteher [15].

(8.21)

(8.22)

If C' is less than zero, the condition of equation (20) is violated in the vicinity of cosPm equal to one. Some of the eigenvalues "'m in the vicinity of eosPm equal to one are such that "'mÄt is outside the region of absolute stability and has negative real part. Recall that applieation of the Gerschgorin Circle Theorem to the original Jacobian [J] results in circles that cross the imaginary axis in regions of negative slope. Benee, that analysis also indicates that some of the eigenvalues of [J] multiplied by At may have negative real parts outside the reg ion of absolute stability of the explicit Euler methode Therefore, the solution wil10ften be

decaying while the error grows. Thus, negative slopes will have serious effects on the stability of the solution if the explicit Euler method is used.

Application of the Gersehgorin Circle Theorem to [J'] yields a bound on the eigenvalues given by a circle that, for ö/[AX]2

2

IC/AXI, is defined by

..

& 2 &2

[ ZRE + + C' ] + [ZIM]2 =

---[AX]2 [AX]4

(8.23)

IfZ is replaced by À,m' the Ie ft-hand side becomes

& 2 2

[

_{- - -}

]

cos2~m

+ _[

~

_]

sin2~m

[AX]2

(8.24)

If &/[AX]2

»

IC/AXI. then the actual eigenvalues are weIl within the bounds given by the Gerschgorin Circle Theorem except near sin~m equal to

zero where they approach but do not exceed the bounds. Clearly. the

Gerschgorin Circle Theorem provides only limited information on stability. For reasons stated in section

3.S.

the implicit multi-step methods of Gear are used to solve the system of ODE's. The Gear method of first order is the implicit Euler method for which the region of absolute stability is the exterior of a circle on the complex plane with unit radius and centre )"At equal to one. Thus. for Burgers' equation. there is a requirement that

[

-

& At [1 -

cos~m]

[AX]2

C'At

1

r

+ [CAt]2 _AX sin 2 ~m

1

1

(8.2S)

If C' is zero. it is easily seen that the condition is always satisfied. The first-order Gear method is unconditionally stabie when applied to the

linearized Burgers' equation. The size of At will be limited by accuracy requirements alone. If C' is negative. the condition is violated in the vicinity of cos~m equal to onei however. the eigenvalues that violate the condition will always have positive real parts. Dence. the solution will of ten be growing faster than the error and the violation of the condition will not have serious consequences. The solution is 'relatively stabie' rather than absolutely stabie. The permissible value of At must be smaller for the full Burgers' equation than for the linearized Burgers' equation to ensure that the error grows more slowly than the solution.

Similar differences in the stability of the solution of the linearized and full Burgers' equations occur for the higher-order Gear methods. Note that the eigenvalues of [J] may not project as far out of the regions of absolute stability as the eigenvalues of [J'] because a pessimistic definition of C and C' was used to define [J']. Dowever.one can conclude that the solution of the full Burgers' equation by a given method is somewhat less stabie than the solution of the linearized Burgers' equat ion.

All the Gear methods of use are absolutely stabie for eigenvalues in a region about the entire negative real axis. but this region gets closer to the real axis as the order of the method is increased. It is of benefit to make &/[AX]2 much larger in magnitude than C/AX. Then. if C' is

negative. the real part of the eigenvalue will be positive only when cos~m

is very close to one, in which case the imaginary part of the eigenvalue will be smalle More specifically, if C' is negative the real part of Àm is

zero for

cos~m 1 + (8.26)

If Ici/AI is much smaller than &/[AX]2 and one makes the reasonable

assumption that Ic'l is at most the same order of magnitude as Ici/AI, then

~m is much smaller than one. Therefore, the imaginary part of Àm is

approximately cJ2Ic'I/& and is thus much smaller than &/[AX]2. Recall that the same conclusion about the desirability of &/[AI]2 being large relative to Ici/AI follows from application of the Gerschgorin Circle Theorem to the original Jacobian [J]. Alternative methods of differencing can also reduce the influence of C' relative to the first term in the eigenvalue.

Conservation-law differencing of the advection term in equation (1) yields

[Ui+l]2 _ [Ui-l]2 4AI

+ & Ui+l - lUi + Ui-l 2 [AI]2

The associated simplified Jacobian [J'] has eigenvalues given by

À

=

m C'cos~m .C . A J-slnpm AI (8.27) (8.28)

Comparison with equation (19) shows that the second term in the expression for the eigenvalue does not have as severe an effect on the stability as in the central-difference case for most values of cos~m.

If it is known that Ui

L

0 everywhere on the domain, then backward differencing of the advection term in equation (1) yields

Ui - Ui-l

[2Ui - Ui-l] _ _ _ _ ₊ & Ui+l - '2Ui + Ui -1 2 [AI]2

The associated simplified Jacobian [J'] has eigenvalues given by

À

=

m -&-[1 -_[AX]2

cos~m]

j[~

+ 2C'

]sin~m

C' [3 - 2cos~m]

(8.29)

(8.30)

In this case the influence of the second term is partially cancelled or reversed by the third term for most values of cos~m since C is positive.

•

The imaginary part of the eigenvalue is slightly reduced in magnitude when C' is negative.

If the sign of Ui varies. then upwind differences can be used instead with a similar influence on the stability. Central differences are

retained for the second-order spatial derivative since diffusion is a bi-directional process.

8.3 Simple Stability Analysis with Moving Nodes

The stability conditions are somewhat more complex when an adaptive method is used. If equation (1) is transformed to the frame (ç.t). where the nodes are stationary and have unit separation (Aç

=

1). and if the spatial derivatives are approximated by central differences in the physical frame. the result is

ü

i - Xi Ui+1 - Ui-1 Ui Ui+1 - Ui-1 Xi+1 - Xi-1 Xi+1 - Xi-1

&

Xi-1[

Ui+1 - Ui Ui

_-

Ui-1

]

(8.31)

+

Xi+1

_-

Xi+1 - Xi Xi

_-

Xi-1

ü

i and Xi are the time derivatives. taken at fixed ç. of the solution and spatial coordinate. respectively, at node i.

As a first approximation, assume the node positions are analytical functions of time alone and do not depend on the values of U at each node. Therefore, the solution at each node has an error, si, but there is no error in node location. Replace Ui by Ui + ai in the equations above, subtract the original equation, and manipulate the viscous term to obtain. af ter discarding terms nonlinear in ai,

=

ei Ui+1 - Ui-1 Xi+1 _ Xi-1

+ &/2 [ 1

Xi+1 - Xi-1 Xi+1 - Xi

+ &/2 [ 1

Xi+1 - Xi-1 Xi+1 - Xi

. •. ei+1 - si-1 [Ul - Xl] _ _ _ _ _ Xi+1 - Xi-1 + 1 ][Si+1 Xi - Xi-1 (8.32) In order to simplify the analysis, one can investigate the behaviour of the closely related equation

= [SECTION 8.3] ~[ai+1 - a i - 1 ] 2 138 + ![a i +1 - 2a i + a i - 1 ] 2 (8.33)

where C' is the most negative value of [Ui +1 - Ui-l)/[Xi+l - Xi-I) in the domain. At the node with the given value of C'. G is the value of

Ui -

ii

2

- - - -

-

-Xi+l - Xi - 1

1 ] (8.34)

Xi - Xi-l • and H is the value of

+

_1_]

Xi - Xi-l

(8.35)

The eigenvalues of the corresponding Jacobian [J') are given by

C' _jGsinJlm (8.36)

Time-integration of Burgers' equation by means of the Gear methods will be most stabie when H is much larger than Ic'l or IGI so that

eigenvalues near the imaginary axis will have small magnitudes. One requirement for stability is that the adaptive method uses small node separations in regions of negative slope so that H will greatly exceed

Ic' I. Nodes should also be concentrated in regions of large variation of the slope to ensure accuracy. Equation (36) should be compared to equation

(19). A solution in which the nodes are moved adaptively wil1 not be more stabie than a solution with fixed nodes. unless the node velocities

provided by the adaptive method result in a value of IGI that is smaller than the value of IC/All. It is desirable. but not necessary. that G be zero.

The value of G will be zero if the node velocities are computed using

~2

_{2L Xi+1 - Xi} 1

(8.37)

If the second derivative in Burgers' equation had been approximated by central differences in the numerical frame. the last two terms in equation

(32) would have been

+ 2& [ei+l - 2e i + ei-I) [Xi+l - Xi-l)2

and

2& _{[Xi+l - 2Xi + Xi-I) [e i +1 - ai-I)} [X i +1 - Xi-l)3

The alternative differencing would result in a slight change in the definition of Hand G. and the node velocities would be given by

=

+ 2& [Xi+1 - 2Xi + Xi-I) [Xi+l - Xi-l)2

(8.38)

(8.39)

This last equation is equivalent to the equation for node-movement of Ghia. Ghia. and Shin [43]. which was reviewed in section 2.3.1.

Neither vers ion of the equation for node-movement guarantees stability since there is no direct mechanism for increasing the

concentration of nodes in regions of negative slope to ensure that H is much larger than

Ic'l.

There is also no mechanism to increase node concentration in regions of large variation of positive slope so as to increase solution accuracy in those areas. Note that. for either form of the equation for node-movement. the nodes tend to move at the

characteristic velocity Ui and would cross in regions of negative slope but for the influence of the second term. The second term has an effect

equivalent to the right-hand sides of the penalty functions described in section 6.2. Mutual repulsion of excessively close nodes is strongest for the first form of the equation for node movement.

The repuls ion behaviour of the equations for node movement can be increased by multiplying the second terms by a factor greater than unity. Ghia. Ghia. and Shin [43] take this approach in their adaptive methode Note that. as a consequence. G is no longer zero.

Any adaptive method that moves the nodes so as to approximately satisfy the first or second equation for node movement can drastically reduce the magnitude of those eigenvalues of [J] that are close to the imaginaryaxis. As a consequence. the Gear methods can use larger time-steps than is possible with nonadaptive methods without inducing

instability. as long as nodes remain sufficiently concentrated in regions of rapid variation of the solution.

8.4 Effect of an Adaptive Eguation on Stability Analysis

The analysis given above is not complete since it was assumed that the adaptive algorithm was exact and that therefore there was no error in the node locations. The effect on the eigenvalues of the Jacobian due to error in the node locations can be illustrated by using the first equation for node movement.

The discretized Burgers equation is reduced to

=

&/2 [ 1

Xi+l - Xi-l Xi+l - Xi +

1 ][Ui+1 Xi - Xi-1

(8.41) if the nodes move as specified in equation (37). Let the error in Xi be designated by yi. The system of ODE's for the propagation of error in equations (41) and (37) has rows defined by

..

&/2 Xi-l[

1 1

][e i +1 - 2e i + e i - 1]

e1 ₊

Xi+1 Xi+1 - Xi Xi - Xi-1

6 yi+1 _ 2y i + yi-1[ 1 1 _rUi 4 Xi+1 - Xi-1 [Xi+1 - Xi]2 Bi _ Xi-1]2

aç2

and

~

yi+l -

y~-I[

1 + 4 Xi+l - Xi-l [Xi+l - Xi]2

& yi+l - yi-l [ _1 _ _ 2 [Xi+l - Xi-I] 2 Xi+l - Xi

1 r2Ui

+ Xi - Xi-l aç2 (8.42)

1

+ 1 _ ][yi+l -2y i + yi-l] [Xi - Xi-l]2 1

-[Xi+l - Xi]2 ____ 1 _ _ _ ][yi+1 _ yi-l] [Xi - Xi-1]2 (8.43)

It is convenient to investigate the related system of equations with rows defined by

." B " 1 " "1 K' "1 " "1 K " 1 "1 el = _[e1+ - _2e1 + e1 - ] - _[y1+ - 2yl + yl- ] __ [y1+ _ yl- ]

and

2 2 2

si + ~[yi+l -2y i + yi-1]

2

+ ~[yi+1 _ yi-l]

2

Here.

K

is the value of

Xi+1~~2Xi-1[

1 1 _y2Ui

+

[Xi+1 - Xi]2 [Xi _ Xi-1]2 aç2

: Xi-1 ]2[

1 1 _y2Ui

+ +

[Xi+1 Xi+1 - Xi Xi - Xi-1 a"ç2

(8.44)

(8.45)

(8.46)

with the largest magnitude.

K'. L,

and

L'

are evaluated at the same node using

K'

= &/2 [ 1 1

r

2Ui (8.47) Xi+l - Xi-l [Xi+l - Xi]2 [Xi - Xi-1]2 3ç2

L

~[

1 [X i _ 1 X i -1] 2 ] (8.48) = + [Xi+l - Xi]2

•

•

and

L'

=

(8.49)

H is also evaluated at the same node using equation (35).

By applying periodic boundary conditions and expanding ei and yi in a Fourier series, one obtains for each Fourier component the system of

equations

~m

₌

+ (8.50)

(8.51) where ~m

=

nm/~l with m

=

0, 1, 2 ••• N. The eigenvalues are given by

+

[8 + L][1 - cos~m] - jL'sin~m

2

J [

[H - L][1 -

cos~m]

+

jL'sin~m

]2 + 4K'[1 -

cos~m]

-

j4Ksin~m

2

(8.52) H, L, and L' are roughly proportional to &/[AI]2, where AI is the average of AIi+1 and AIi at the node with maximum magnitude of K. K and K' are roughly proportional to &/[AX]3[a2u/a

ç

2]. For sufficiently small AX or for a sufficiently small numerical second derivative of U, Hand L will always greatly exceed

IK

and L' will always greatly exceed

JK'.

The tot al number of nodes that must be used depends on how weIl the adaptive method minimizes the numerical second derivative of U. Given a sufficient number of nodes, the imaginary part of eigenvalues close to the imaginary axis will be small, as is required for stability. the eigenvalues for a particular Fourier component are

H[1 - cos~m] (8.53)

and

L[1 - cos~m] + jL'sin~m (8.54)

If the nodes are evenly separated in the reg ion of maximum K, L equals H and L' is zero. The eigenvalues are then the same as obtained by the simple analysis used previously •

In general, an improved adaptive algorithm would monitor the

numerical second derivative of U and related quantities and increase node concentration where required, while approximately satisfying equation (37).

8.5 Effect of the MFD Method on Stability

A simple MFD equation for node movement with wl

=

1. including penalty function terms, with Q'

=

Co

=

0 may be written

rii

C 3

-[AXi ]2

AXi+l

.. 1 Ui+l - Ui-l rüi+l - 2Ûi + U1- ] _ _ _ _ _

Xi+l _ Xi-l

n

.

1 Ui+1 - Ui-1

+ _[Ui+1 - 2Ui + U1- ] _ _ _ _ _

2 Xi+1 - Xi-l

+

W2~[Xi+l

- 2Xi + Xi-l][Ui+1

2 Xi₊1

Ui-lJ2 • Xi-l where

n

= 0/[1 + OAt/4].

(8.56)

A detailed stability analysis of the MFD equations quickly becomes unwieldy. An heuristic analysis will be given instead. If the time derivative of U at node i is eliminated from equation (56) by using equation (31). the result. af ter rearrangement. is

4

X1

.

.

-

Ui +

~[

1 Xi+l - Xi

1

Xi - Xi-l W2[Xi+l - 2X i + Xi -1_][U,]2

+ W2-[X1+ - 2Xi + Xi-1] [U,]2 U . 1 2 +

[Ü

i +1 +

Ü

i -1_]U' +

~[Ui+l

- 2Ui _{+ Ui-1]U'} 2 + ___ 1 _ _ ][Ui+1 Xi - Xi-1

rii

C2 C2 C3 - +

-[AXi]2 AXi+l !Xi

2Ui + Ui-1 lU'

(8.57)