THE SIMPLIFIED MOVING FINITE ELEMENT METHOD
APPLIED TO COMPRESSIBLE FLOW
by
James R. Valentine
B.S., University of New Hampshire, 1974
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado at Denver in partial fulfillment
of the requirements for the degree of
Master of Science
Department of Mechanical Engineering
1988
This 'thesis for the Master of Science degree by
James R. Valentine
has been approved for the
Department of
Mechanical Engineering
by
Date
iii
Valentine, James R. (M.S., Mechanical Engineering)
The Simplified Moving Finite Element Method Applied to
Compressible Flow
Thesis directed by Associate Professor John A. Trapp
The Simplified Moving Finite Element (SMFE) method is a
moving mesh finite element scheme developed, by John Dukowicz in
1984. A moving mesh scheme is useful in resolving steep
gradients since nodes tend to move into and thus be
concentrated in the gradients. Here, the SMFE method id
applied to Burgers Equation and to a compressible flow problem.
The compressible flow example is a shock tube, and the analytic
solutition to shock tube flow is developed and compared to the
SMFE solution.
The form and content of this abstract are approved,
recommend its publication..
Signed
iv
ACKNOWLEDGEMENTS
I would like to express my thanks to John Trapp for the
opportunity to be able to write this thesis and for the great
amount of help he gave in the project. I would also like to
thank Ken Ortega and William Clohessy for serving on my thesis
committee, Faye Waitman for help with the use of the word
processor and the UCD Mechanical Engineering Department for
making the research possible.
CONTENTS
CHAPTER
I. INTRODUCTION................................ 1
II. ANANLYTIC SOLUTION FOR SHOCK TUBE FLOW...... 3
Governing Equations For OneDimensional Flow... 3
Shocks...................................... 7
Reflection of a Shock........................... 10
Contact Surfaces............................ 12
Method of Characteristics................... 13
Expansion Waves............................. 17
Shock Tube Flow............................. 20
III. NUMERICAL METHODS........................... 24
Conventional numerical methods.............. 24
The simplified moving finite element method.... 26
Penalty Functions........................... 32
The SMFE Method Applied to Burgers Equation.... 36
The SMFE Method Applied to Compressible Flow... 40
IV. SMFE CODE DEVELOPMENT....................... 48
V. RESULTS..................................... 53
Burgers Equation............................ 53
Compressible Flow.....'......................... 58
VI
Penalties....................................... 67
Other Comments.................................. 67
BIBLIOGRAPHY
70
Vll
FIGURES
Figure
2.1 Expansion wave in the shock tube.................. 17
2.2 Density profile in the shock tube................. 20
2.3 Shock tube flow in the xt plane.................. 21
5.1 Burgers equation 1................................ 54
'5.2 Burgers equation II............................... 56
5.3 Burgers equation III.............................. 57
5.4 Density........................................... 61
5.5 Pressure........................................ 62
5.6 Velocity.......................................... 63
5.7 Reflected shock, density.......................... 64
5.8 Reflected shock, pressure......................... 65
5.9 Reflected shock, velocity......................... 66
CHAPTER I
INTRODUCTION
This thesis presents a model for the flow of a gas in a
shock tube using an adaptive mesh finite element method.
A shock tube consists of a tube of gas separated into
two compartments by a membrane. One compartment contains gas
at a high pressure (the driver gas) and the other contains gas
at a lower pressure (the driven gas). The membrane is broken
and the high pressure gas moves into the low pressure region.
An expansion wave propagates into the high pressure region and
a shock wave and a contact surface propagate into the low
pressure region.
At the contact surface and shock there are very steep
gradients in fluid properties and to accurately resolve these
gradients a very fine mesh is needed. In the adaptive mesh
method used here node locations are adjusted during the
solution process so the mesh is concentrated in areas where a
fine mesh is needed; i.e., where steep gradients exist. In
areas where gradients are not so steep and little resolution is
required, the nodes are spaced farther apart. Thus fewer nodes
are required overall than would be in a conventional fixed mesh
method and fewer calculations are required to follow the flow.
2
In chapter II, relationships between properties across
a shock wave, a contact surface and an expansion wave will be
derived and from these the analytic solution for flow in a
shock tube will be developed. The adaptive mesh finite element
solution can then be compared to the analytic solution.
Chapter III explains the adaptive mesh method and Chapter IV
describes the code development. The final chapter gives the
results of the method, both for a simple case using Burgers
equation and for the more complicated compressible flow shock
tube example.
3
CHAPTER II
ANALYTIC SOLUTION FOR FLOW IN A SHOCK TUBE
GOVERNING EQUATIONS FOR ONEDIMENSIONAL FLOW
The basic governing equations for fluid flow that will
be used here can be derived by applying three physical
principles to a control volume in the flow:
1. Mass must be conserved.
2. The time rate of change of linear momentum of
fluid in the control volume is equal to the net
force exerted on the fluid.
3. Energy must be conserved.
Omitting the details of the derivations, (the
derivations in this section are standard, the details can be
found in reference [1]) the following governing equations in
integral form result for inviscid, adiabatic flow with no body
conservation of mass:
^dV+
r 3t
pVdS=0
forces:
4
conservation of linear momentum:
3(pV)
r 3t
dV+
(pV)VdS=
conservation of energy:
PdS
9(pE)
r 3t
dV+
(pE)VdS=
s
PV.dS
s
2
where E=total energy/mass=i+^
i=specific internal energy
u=velocity.
Assuming onedimensional, steady, uniform, inviscid,
adiabatic flow with no body forces, and considering a very thin
control volume with ends normal to the flow direction, the
integral equations above reduce to the following form where the
subscripts 1 and 2 refer the properties of the fluid at each
end of the control volume.
continuity: P1U=P2U2 (II1a)
2 2
momentum: P +pi u =P2+P2U2 (111b)
energy:
P1U1 [hi
2 2
^1 1 r U2 "I
+J = p2u2Lh2+2~J
(II1c)
or, if pi ui 5^0,
2 2
U1 u2
energy: h+=h+ (II1 d)
p
where h=specific enthalpy=i+ .
The second, simpler form of the energy equation (II1d) will be
used here except in the case of a contact surface.
5
In addition to these equations we have two equations of
state, P=P(p,T) and i=i(p,T). For a perfect gas, the
expressions for the equations of state are
P=pRT (II2)
and
i=cvT (or h=cpT). (II~3)
The onedimensional flow equations above can be used to
analyze flow through a stationary normal shock wave (flow is
perpendicular to a normal shock, otherwise the shock is called
oblique).
In a shock tube, however, the shock isn't stationary,
but the equations can be easily rewritten to apply to the case
of a moving shock by using a coordinate system that moves with
the shock. Let U, ui and U2 be the velocities of the shock,
the fluid in front of the shock and the fluid behind the shock
with respect to a fixed coordinate system as before, and u] and
U2 be the corresponding velocities in the primed coordinate
system moving with the shock. Obviously, the thermodynamic
properties i, p, T and h are the same in both coordinate
systems. Letting the thickness of the control volume become
infinitesimal, the onedimensional governing equations (II1a)
(II1d) with respect to the moving coordinate system can be
written
continuity: puj=p2U2
6
momentum: P+pUi'2=P2+P2u2^
>2 .2
u U2
energy: h+=h2+2~
But,
ui' =u] U (IIMa)
and
u2=u2_U (114t>)
so, for a moving shock,
continuity: p (u U) = P2(U2U) (II5a)
2 2
momentum: P +pi (ui U) =P2+P2(u2U) (II5b)
(utU)2 (u2"U)2
energy: h^ + =h2+^ (II5c)
The equations above are valid for onedimensional, uniform,
inviscid, adiabatic flow with no body forces.
The propagation speed of small sound disturbances is
also important in the analysis of compressible flow. The sound
speed a can be found to be (this is a standard derivation and
can be found in reference [1])
or, for a perfect gas
a=/TP7p=/rRT . (I16)
The mach number is then defined as M= .
a
(I17)
7
SHOCKS
In the following analysis of shocks, relationships will
be developed for the upstream and downstream properties of a
stationary shock then applied to the more general case of a
moving shock.
In general, in the analysis of stationary shocks, the
upstream conditions pi P h (or ii), u and T are known.
There are then three governing equations and two state
equations to solve for the five downstream properties pÂ£, P2
h2 (or ij), U2 and T2. For a perfect gas these are collected
and recorded below:
P1 U1 = P2U2 (II1a)
P"l +P1 U1 =P2+P2u2 (II1b)
2 2
U1 u2 hl+2~=h2+2~ ("II 1 d)
P=pRT , (H2)
h=cpT . (II3)
The equations above can now be used to develop
relations between the upstream and downstream fluid properties
across a shock for a perfect gas.
Equation (II3) can be written h=CpT=[^ry]T
and along with equations (II6) and (II7) substituted into the
energy equation (II2d) to obtain
YRi M1YRT1 yR, M2YRT2
[^]T,
2
2
8
or
t2 HirlHr
(II8)
which relates upstream and downstream temperature across the
shock. Using equation (I16), this can be rewritten as
a2
a . r T1 i2
(II9)
L.
Next, rewriting the continuity equation (II1a), we have
p2 u Mi a^ Mi/yRTi
PI u2 M2a2 ^2yrRT2 M2/fJ '
n[^hV'
or
P2 U1 M1
P1 u2 M2
(1110)
The linear momentum equation (II1b) can be rewritten as
'1
1 +
2i
P1U1
=P2
1 +
2*i
P2U2
p2
and, using equation (II6) to eliminate p/P, this results in
P1
Then, using the definition of the Mach number, (II7), we get
' YU?' =P2 r 2*i Yu2
2 2
. a. L a2
P2 1 +YM
Pl 1+YM2
(II11)
The upstream and downstream mach numbers can also be
related. From the equation of state for a perfect gas we get
p2_P2t2
P1 P1T1
9
Substituting into this equations (I18), (1110) and (1111)
results in a polynomial of M2 in terms of M] and Y. The
solutions of this polynomial are
(1112)
The first of these solutions is trivial. The second solution,
equation (1112), expresses the downstream mach number M2 in
terms of the upstream mach number Mi and Y for a shock of
finite strength. Substituting equation (1112) into equations
(II8), (II9), (1110) and (1111) yields
2 2
T2 h, (2YM1(Y1))((Y1)M1+2)
Ti h.
(Y+1)2M?
ai'
,(2YM1(Y1))((Y1)M1+2)l1/2
(Y+1)2M?
p2 uj (Y+1 )M]
p1 u2 2+(Y1)M2 ,
?f=1+7TT(M'1).
(1113)
(II14)
(1115)
(1116)
Equations (II12)(II16) can be used to solve for the
downstream properties explicitly if the upstream properties are
known.
10
An additional relationship will be needed for analysis
of flow in a shock tube. Using the definition of the Mach
number, equation (II7), we can write
M _f2 _U1 a2u2 U1 ~u2
1 a ^a aa2~ aq
Now, substituting equation (1112) for M2 and equation
CII1 4) for a2/a and simplifying yields
uiu2 2(Mf1)
(H17)
a! M^Y1) .
Equations (II8)(II17) can also be applied to a
moving shock. As noted earlier, a moving shock can be modeled
as a stationary shock by using a primed coordinate system that
moves with the shock so, u} and U2 will be used in place of u
and U2 in the stationary shock equations and using equations
(IlJla), (Il^b) and (II7), M and M2 can be written
u} uiU
M1
M2=
a1 a1
U2 U2_U
a2 a2
(1118)
(1119)
As before, u is the velocity in the stationary reference frame,
u* the velocity in the moving frame and U the propagation speed
of the shock in the stationary frame.
REFLECTION OF A SHOCK
Now, the shock relations developed in the previous
section will be used to analyze what happens when the shock
hits the solid end wall of the shock tube. When a shock hits a
11
rigid surface normal to its propagation direction, the velocity
of the particles induced by the shock immediately behind the
shock must be stopped, so the boundary condition requires a
reflected shock coming from the wall. The strength of the
reflected shock must be enough to stop the motion of the
particles induced by the incident shock, so as it propagates
into the particle motion induced by the incident shock, it
leaves motionless fluid behind it. In the analysis here, it
will be assumed that the incident shock is propagating into a
motionless gas, as is the case in a shock tube. So, the fluid
velocity behind the reflected shock and in front of the
incident shock are both zero, and the reflected shock
propagates into the conditions behind the incident shock.
Applying equation (1117) to the reflected shock and noting
that ur2=0, ari=ai2 and ur1=ui2 results in
or,
UpiUp2 Url~ur2 Uj2 2(M^i1)
ar1 ~ ar1 _ai2Mpi(Y+1)
ui2 rai22(Mr1 ~1)
ai1=Lai1 Mr1(Y+1)
(1120)
where Mr is given by equation (1117) since we are considering
a moving shock. Now, applying equation (1117) to the incident
shock and noting that Uii=0 we get
t 2
ui1~ui2 Ui1~ui2 ui2_2(Mi11)
ai1 ai1 ai1 Mn (Y+1)
(II21)
12
where once again, M^i is given by equation (1118). Equating
equations (1120) and (1121) and rearranging results in
which can be solved numerically for Mr if Mii is known (which
will be the case for the shock tube). Once Mrl has been found,
the fluid properties behind the refleceted shock can be found
using equations (II13)(II19) since the properties in front
of the reflected shock are known.
CONTACT SURFACE
A contact surface is a discontinuity across which
pressure and velocity are constant, but other properties are
discontinuous. In a shock tube, a contact surface forms at the
interface between the high and low pressure gases and moves
with the fluid surrounding it. With no viscosity, the contact
surface retains its identity as a discontinuity. Since the
contact surface moves with the fluid surrounding it, no fluid
flows through it so u =U2 and u]=U2=0.
and moving with the contact surface, we have the governing
equations from earlier,
(Mj11)2(Y+1)2 /2
,(2YMil(Y1)((Y1)Mi1+2).
(1122)
Considering a control volume around the contact surface
continuity: Piu=P2U2
(II1 a)
(II1b)
13
>2 r 2
i r i f p u2 i
energy: p<ui [hj J = p2U2l.h2+2 J . (H1c)
Considered together, the continuity and momentum equations show
that P =?2 across the contact surface. No relationship between
the upstream and downstream density and enthalpy can be
determined from the governing equations above.
METHOD OF CHARACTERISTICS
The method of characteristics is a mathematical
technique for generating the solution of a system of hyperbolic
partial differential equations from a known initial state. The
governing equations for onedimensional compressible flow make
up such a system. The partial differential equations are .
manipulated algebraically, resulting in a system of total
differential equations, the compatibility equations, which are
valid along certain paths, the characteristic curves which are
specified by characteristic equations.
In the onedimensional flow problem, the solution
consists of the velocity and thermodynamic state of the gas at
any location and time. The thermodynamic state of the gas is
specified by any two independent thermodynamic properties,
giving a total of three unknown dependent variables, which in
this case will be chosen to be p, p and u. The three dependent
variables are functions of the independent variables x and t.
The three governing equations are thus sufficient to provide
14
solutions for p, p and u. Any other thermodynamic properties
in the equations must be expressible as functions of p and p
through the state equation(s). The method of characteristics
must provide three total differential equations if a complete
solution is to be found. The three total differential
equations will be the compatibility equations, each of which
will be valid along a characteristic curve whose slope is
specified by a characteristic equation. The characteristic and
compatibility equations will be functions of the dependent
variables, so the dependent variables must be continuous in the
region in which the method of characteristics is to be used,
although their first partial derivatives may be discontinuous
(it can be shown that these discontinuities, if present, will
propagate along characteristics [10]). If discontinuities in
the variables themselves exist, the compatibility equations and
the slopes of the characteristic equations will also be
discontinuous. In the shock tube problem, such discontinuities
exist at the shock and contact surface so the method of
characteristics can be applied to the regions between the
discontinuities, but not to the whole shock tube at once.
To apply the method of characteristics to one
dimensional, inviscid, adiabatic flow, we begin with the
differential form of the governing equations for one
dimensional, inviscid, adiabatic flow with no body forces where
15
the energy equation has been written in terms of entropy
([1],p.184),
continuity: pt+(pu)x=0
momentum: put+puux+px=0
energy: st+usx=0
where the subscripts refer to partial derivatives. Note that
there are four unknowns in the equations above. To eliminate
entropy, the energy equation can be replaced by (pt+upx)
c2(pt+upx)=0 since the flow is isentropic (inviscid, adiabatic
flow is isentropic). Omitting the details of the derivation
(which can be found in Zucrow and Hoffman) three compatibility
equations result, each of which is valid along a corresponding
characteristic curve in the xt plane. They are
du+=0 along the curve with slope
pa & a+u
du^=0 along the curve with slope
pa 6 au
and ds=0 along the curve with slope jj .
The characteristic with slope 1/(a+u) is called the C+
characteristic and the one with slope 1/(au) the C_
characteristic. It can be seen that sound waves propagate
along the C+ and C characteristics and the third
characteristic describes pathlines, the paths of the fluid
particles in the xt plane. Each point in the xt plane will
have one of each of these characteristic curves passing through
it and integrating the corresponding compatibility equations
16
along the characteristic curve from a previous time to the new
time will give three equations for the fluid properties at the
new time. In the shock tube, initially all particles on the
same side of the membrane have the same entropy and since the
flow is adiabatic and inviscid (except in the immediate
vicinity of the shock) the entropy of each particle remains
constant unless it passes through the shock (all particles
passing through the shock experience the same change in
entropy). Therefore the third compatibility equation, ds=0
along a curve with slope 1/u, will carry the same information
along all characteristics throughout the same region of flow
and thus can be neglected. So to analyze the inviscid flow in
a shock tube only the first two chararcteristics are needed.
Integrating the compatibility equations along the C+ and C_
characteristics and using thermodynamic relations for
isentropic changes in a perfect gas gives
2a
J+=u+^7Y=constant along C+ characteristic,
2a
J=u;pi=COnstant along C characteristic,
where J+ and J_ are called the Riemann invariants,
each point in the xt plane,
a=jj(J+J)
u=^( J++J)
where J+ is the value from a previous time on the C+
characteristic and J is the value from a previous time on the
(1123)
(1124)
So, for
(1125)
. (1126)
C_ characteristic.
17
EXPANSION WAVES
The breaking of the membrane in a shock tube creates a
low pressure disturbance that propagates into the region of
high pressure gas as shown below in the xt plane.
figure 2.1: Expansion wave in the shock tube
For convenience, regions 3 and 4 are labeled as they will be
used later in finding the analytic solustion to the shock tube
flow. Region 4 contains the undisturbed high pressure gas and
region 3 the gas behind the contact surface as the contact
surface propagates to the right. The characteristic C_ with
slope 1 /(un~an)=1 /aij represents the head of the expansion wave
18
and C_2 with slope 1/(11333) the tail of the wave. Since all
C+ characteristics begin in the undisturbed region
2a 4
J+=Ujj+^3y along all C+ characteristics. Therefore, J+
has the same value at all points in the expansion, so both J+
and J_ must be constant along any C_ characteristic. This
means that u and a are constant along a C_ characteristic, so
these characteristics must be straight lines.
To find p, T, u and a at any point in the expansion
wave, consider the C+ and C_ characteristics passing through
the point as shown in the previous figure. Since J+ is the
same throughout the wave,
2an 2a
J+=U4+^j=u+^Tj=constant ,
or,
a (Y1)(uiju)
=1 +
an 2aij
(1127)
for any point in the expansion wave. Substituting equation
(II6) for a and aij results in
T
tT
1 +
(Y1 ) (uiju)
2aij
(1128)
For an isentropic thermodynamic process from an initial
state 1 to a final state 2,
Y
^l_r piiY_rtyi
?2~ P2 ~LT2J
and considering the expansion wave to be isentropic (adiabatic
and inviscid), we can write for a perfect gas
pji lpiij lV
19
Combining this with equation (1128) results in
and
P
pT
P__
P4"
f (TD(uhu)
1+
2Y_
Y1
1 +
2aij
(Y1) (ui,u)
2a ij
Y1
(II29)
(II30)
which are valid at any point in the expansion wave.
The expressions above relate pressure, density and
temperature at any point in the expansion wave to the known
properties in region 4 and the velocity at the point in the
wave. To find the velocity, consider the C_ characteristic
through the point with slope 1/(ua). As noted earlier, this
characteristic is a straight line, so every point (x,t) on it
must satisfy ua=(xxg)/(ttg) where xg and tg are the initial
location and time of the disturbance. The expression for a/aij
in equation (1127) can be substituted for a, resulting in
"Y+1
xxg
tVar
(Y1)
Ui
(II3D
which is valid throughout the expansion wave.
Now the properties at any location in the expansion
wave can be found in terms of the known properties in region 4.
By considering points along the C_2 characteristic at the tail
of the expansion wave in the previous figure, the properties in
region 3 can be found in terms of those in region 4.
20
SHOCK TUBE FLOW
As noted earlier, a shock tube consists of a tube
divided into two sections by a membrane. One section contains
high pressure gas at rest and the other section low pressure
gas at rest. The fluid properties in both of these regions are
known. When the membrane is broken, an expansion wave
propagates into the high pressure region and a contact surface
and a shock wave into the low pressure region. At some time
after the membrane is broken, the density profile will appear
as in figure 2.2 below.
figure 2.2: Density profile in the shock tube
Each region is labeled as shown in the figure for convenience,
with the shock separating regions 1 and 2, the contact surface
21
separating regions 2 and 3 and the expansion wave separating
regions 3 and 4. xq is the membrane location. In the xt
plane this appears as in figure 2.3.
figure 2.3: Shock tube flow in the xt plane
The C+ and C characteristics in each region are parallel, so
the fluid properties throughtout these regions will be
constant.
The expressions previously developed for properties
across a shock and an expansion wave are sufficient to
determine the unknown properties in regions 2 and 3. Repeating
the shock relations from earlier, we have
p2 2Y 2
^I^CM,1) (1116)
uru2 2(M?1)
Tj =Mt (y1 )
(II17)
22
Rewriting equation (1116) above for Mf results in
M=
pi,
Y+1
ni/:
2Y
+1 (H32)
Now, substituting equation (1132) into equation (1117) and
noting that for this particular case u=0, we have
_ra1 ir^2_. I
u2L y J L p ^ J
2Y
^(Y+1) + (Y1)
I'M
(1133)
which is an expression for U2 developed by analyzing the shock.
Remembering from the section on contact surfaces that
U2=U3 across the contact surface, we will now develop an
expression for U3 by analyzing the expansion wave. Repeating
from the analysis of the expansion wave, we have
P
1 +
(Y1) (uifu)
Y1
2aif
(1129)
for any point in the wave. Applying equation (1129) above at
the tail of the expansion wave where P=P3 and u=U3> noting that
ui=0 and rewriting results in
u3=
2ai
Y1
,p2
Pn
Y1
2Y
(113^)
Equating equations (1133) and (113^) (since U2=U3), noting
that P2=P3 and rearranging gives
P PJ.
P1~P1
raf ?2 .
i1'
1
/2Y(2Y+(Y+1)(P2/Pl1))
(H35)
23
This expression can be solved numerically for P2/P1, so
P2=P3 can be found. Then, from the shock and expansion wave
relations, the rest of the properties in regions 2 and 3 can be
found.
CHAPTER III
NUMERICAL METHODS
CONVENTIONAL NUMERICAL METHODS
When a very steep gradient or discontinuity such as a
shock is encountered in a numerical model of a flow, large
oscillations in the calculated flow properties may develop
behind the shock unless special provisions are taken to prevent
this. The oscillations are not instabilities; their amplitudes
will remain constant and will not increase with time ([7],
p.328). In the case of a centered finite difference
approximation, the oscillations are the actual solution of the
finite difference equations ([8], pp.1615).
Shockfitting or shockpatching was an early technique
to deal with discontinuities in the flow. In using this
method, the flow on either side of the shock is treated as
continuous and the changes in flow properties across the shock
serve as additional boundary conditions for the inviscid flow
on either side. At the shock, relations between upstream and
downstream properties and information from the inviscid flow
are used to solve for the properties and location of the shock.
25
A problem with this technique is that there is no way to
account for shocks that may arise during the solution process.
A newer technique for dealing with discontinuities was
first published by Richtmyer and von Neumann in 1950. In this
method, an artificial viscosity term is introduced into the
numerical equations and acts as a dissipative force to smear
the shock out over 3 to 5 nodes, smoothing the oscillations in
much the same manner as viscosity acts in a real fluid.
Richtmyer and Mortons method was to add a quadratic velocity
gradient term to the pressure term in the conservation of
linear momentum and energy equations. The additional term then
appears as a second derivative of velocity in the momentum
equation and a first derivative of velocity in the energy
equation, which is the same form as the physical viscosity term
in both equations. Since the additional term is quadratic, it
acts like a viscosity with a coefficient term that varies with
the velocity gradient. A constant viscosity coefficient for
this term would tend to make a very strong shock thinner than a
weaker one, but if. the viscosity increases with the strength of
the shock as in the case here, all shocks will be of about the
same thickness. Artificial viscosity can also be added
implicitly to finite difference approximations of the governing
equations by a scheme where the leading term in the trucation
error is a second derivative of velocity.
26
Since this method was first introduced, variations of
it have been developed, but the principle remains the same. In
the finite element code that will be developed here, a very
small artificial viscosity will be added by including a second
derivative of density in the continuity equation, a second
derivative of mass flux in the momentum equation and a second
derivative of energy in the energy equation. Thus all
parameters are smoothed slightly. Other methods, such as
Godunov's method, have also been developed for dealing with
discontinuities in fluid flow.
THE SIMPLIFIED MOVING FINITE ELEMENT METHOD
As discussed in the previous section, artificial
viscosity smears a shock, spreading it out over several nodes
and dissipating any oscillations that may form behind it. An
alternate technique for getting several nodes in the shock
region is to have a mesh that moves into or moves with the
shock. This is done with a moving mesh method.
In a conventional fixed mesh finite element method, the
approximate solution of a system of partial differential
equations is found at each mesh point by minimizing the
residual or error of the approximate solution over all mesh
points. In the minimization process, only the amplitudes of
the approximate solution at each mesh point are variable. In
27
the moving finite element method (MFE) introduced by Miller and
Miller ([5], [6]), the nodes or mesh points are allowed to
move. For this case, the approximate solution is a function of
both the amplitude and location of each node which leads to two
sets of coupled equations when minimizing the residual, one for
the amplitudes and one for the node locations.
The simplified finite element (SMFE) method, which will
be used here, was introduced by John Dukowicz [2], and provides
a simpler method for moving the mesh. The MFE method solves a
PDE of form qt=L(q) by minimizing [q^LCq)1^ while the SMFE
method minimizes QSqxL(q)2 with respect to the mesh
velocity where Q=qt+Sqx is the total time derivative of q with
respect to the mesh motion and S is the mesh velocity. The
equations resulting from^the SMFE minimization problem are
simpler than those resulting from the MFE method. As an
introduction to the SMFE method, consider the governing
equations for fluid flow in conservation form ([7], p. 301)
continuity: p^+mx=0 , (III1a)
momentum: m^+(m2/p+p)x=0 , (IIIlb)
energy: et+((e+p)m/p)x=0 . (III1c)
where subscripts indicate partial derivatives, m=pu=mass flow
rate and e=total energy/volume = pE = p(i+u^/2). Equations
(III1a)(lil1c) above can be rewritten
continuity: pt+upx+pux=0 ,
(III2a)
28
momentum: m^+umx+mux+px=0 , (III2b)
energy: et+uex+eux+(pu)x=0 . (III2c)
These equations are of the form
Qt+uQx+LU)=0 (HI3)
which will be the general form of the PDE used in the SMFE
method. For a mesh moving at velocity S (S^ will in general be
different at each node), the total time derivative of the fluid
property q following the mesh motion will be
<5=^=qt+Sqx (IIIU)
The general form of the SMFE method PDE in equation (II13) can
then be written
<)+(uS)qx+L(q)=0 (III5)
or, <5=(Su)qxL(q) (III6)
Each of the variables q, <5, u and S in equations (III5)
and (III6) above can be approximated in the k interval between
nodes x^ and x^+i by
i + 1
qk(x)= l qj
j=i J
i + 1
4k(x)= l (III7b)
j = i J
i + 1
uk(x)= l Uj4>k(x) (III7c)
j=i J
i+1
Sk(x)= l Sj<(,k(x)
j=1 J
(III7d)
29
where qj, <5j Uj and Sj are the values of the parameters at the
j node at location xj and Â£ and Â£+i are linear shape
functions defined over the interval [xj,xi+i] by
*Â£(x)
x~xi+1
xixi+1
*
k
i+1
xx^
xi+1"xi
(III8a)
(III8b)
and are zero everywhere else. The residual or error of the
finite element approximation over interval k is given by
Rk= I [ <5j + (UjSj ) 3] (frj (x) +L (qk
(III9)
89
where r1 as evaluated using (III7a), (III8a) and (III8b) is
OX
Bq^ qi+iqj
3x xi+1Xi
(III10)
(5 can be approximated by a numerical equation for (III6),
4Qi(S,q,x)(Siui)[S]1[L(q)]i (III11)
ft rv
where [r^]. and [L(q) ]; are approximations of ^ and
dX 1 dX
L(q) at node i, so each term in equation (III9) can be
evaluated. Mqk) is evaluated separately since it may contain
second derivative viscosity terms. Since the first derivatives
of qk are in general discontinuous at each node, the second
derivatives will not exist and will have to be evaluated
through mollification. In mollification, the function is
smoothed or rounded over an interval
point where the first derivative is discontinuous. Then both
the first and second derivatives exist and are continuous and a
30
terra containing the second derivative can be evaluated in the
limit as 5x>0.
Now we want to minimize the square of the norm of the
residual with respect to the mesh motion. Defining a
functional I,
1= R 2=
fxN _ _
R2dx=Â£
x k
fxi+1
(Rk)2dx
CIII13)
the minimization gives us for N nodes the equations
31
^=0 i=1,2,...,N
(III14.)
From the form of the residual equation (III9), it can be seen
that Sj^ is found only in Rk and Rk1, so
31 3
3Si 3St
fX,
(Rk1 )2,jx+
rxi+i
xii
(Rk)2dx
=0 i=1 ,2......N
or,
Rl',!isrdx<
Xi1 1
>xi+1. 2Rk
Rk7dx=0 i=1 ,2....N
(III15)
Also, from (II16) and (III9), it can be seen that Rk will
contain linear (S^u^) and (S^+iUi+i) terms and Rk1 linear
(Si_Ui) and (S^u^) terms. Therefore, (III15) will result
in N equations, the i^h of which will contain linear
(S^_ 1 Ui 1), (SjUi) and (Si+iUi +  ) terms, so equation
(II115) can be written in matrix form
[A][Su>[B] (III16)
where the coefficient matrix [A] is a tridiagonal NxN matrix.
The variable u will always be an unknown in the problem; in the
31
simplest case of one equation and one unknown, q=u, and in the
case of several unknowns (q's), u can be evaluated in terms of
one or more q's. So, equation (III16) gives us N equations to
solve for the N unknowns, Si,...,Sn. The new mesh locations
can be found since
(III17)
dt
can be discretized. Then, the numerical approximations for q,
equation (III11) can be used to advance q to the new time
level. (In the actual code, equations (III11), (III16) and
(III17) are solved simultaneously.)
When the solution of a system of n PDE's is desired, as
in the case of compressible flow, the residual of the
approximation of each PDE is formed, then the squares of the
norms of the residuals are summed as
I=WTtIt+WT2l2+.+WTnIn (III18)
where WT WT2, ... WTn are weighting terms. As in the case
of a single PDE, the mesh motion equation is determined by
minimizing the derivative of I with respect to the mesh
velocity.
It should be noted that in certain instances the matrix
[A] in the mesh motion equation [A][Su]=[B] may become
singular. This will occur when three points line up, that is
when 9qk V3x=3qk/3x. When executing the program, the effect
of [A] approaching singularity is that the mesh velocity at
32
a node becomes very large and that node may pass an adjacent
node, tangling the mesh. To avoid this problem, so called
penalty functions or regularization terms are introduced into
the minimization equation.
PENALTY FUNCTIONS
The penalty functions are added to the minimization
equation and are of a form so that they are insignificant
unless the internodal distance approaches a specified minimum
distance, at which point they increase in magnitude and push
the nodes apart, preventing mesh tangling. The penalties
consist of a "viscosity" term, e, and a "spring" term, S, both
of which are functions of the internodal distance and property
gradients. The penalties are added to I in equation (III18)
and the sum is minimized with respect to the mesh motion S^, so
the new minimization problem with penalties included becomes
91 ,9C&ASS)2
+ r/ =o
9Si 9Si
(III19)
where AS=SjSi_ As the minimum internodal distance is
approached, e acts as an viscosity force to slow the motion of
the nodes and S as a spring force to keep the nodes apart.
Expanding the derivative of the penalty function term in the
minimization equation gives
UJ. _2 _2
gS.Eisi1+(ei+Â£i+l)Siei+lSi+l=eiSxÂ£i+iS:[ + i
91 _2
h
33
or,
31 2 2 2 2
gs.2^si1llii )+2(Ei+ei+i )(Siui)2'ei + 1 (Si + iui + 1 )
=2Â£iSi2ei+iSi+i 2ef+i (ujuii )2ef(uiUi+i ) (III20)
where ej and are the values of e" and S at node i. The 2
multiplying each term can be dropped, since it acts only as a
scaling factor for the penalty terms. The unpenalized
minimization equation is [A][Su]=[B], so the penalty functions
can be added to the unpenalized problem simply by adding the
appropriate terms in (III20) to matrices [A] and [B].
To see the effect of the penalties on mesh motion,
consider and to be very small for a smoothly moving mesh,
so the minimization equation in this case is unaffected by the
penalty terms. The penalty terms at nodes i1, i and i+1 are
i1 : EiiSi_2+(ei+ei)Si_iefsi=eiSiieiSi (III21a)
i; e'iSii + (Ei + ei+i )Siei+]Si+1=eiSiei + lSi+] (III21 b)
_2 _2_2 _2
i+1: Â£i+lSi+(Ei+i+Â£i+2)si+1Â£i+2si+2=ei+1Si+1~Â£i+2si+2
(III21c)
Now, suppose the property gradients between nodes i1 and i
and nodes i and i+1 begin to approach the same value. As [A]
becomes singular, the magnitude of becomes very large,
causing node i to approach an adjacent node. Suppose that node
i approaches node i1 and that this causes to become very
large. Assuming that the rest of the mesh is moving smoothly,
that is, the velocities of all the rest of the nodes are about
the same, and that e and S at all other nodes are small and can
be neglected, the penalty equations above reduce to
(III22b)
(III22a)
i + 1: 0=0
(III22C)
Equations (III22a) and (III22b) above are identical and can
be rearranged to form
If Sj_ is very small, Sj/ei is also small, so a large has
acted as a "viscous" force, slowing the motion of node i and
causing nodes i and i+1 to move together. As the magnitude of
Sj increases, Sj/Efi beomes larger and beomes larger than
S^_i. So, a large % acts as a "spring" force, pushing nodes
i1 and i apart. If S^ had been such that node i approached
node i+1, ei+T would become large and would thus tend to make
nodes i and i+1 move together. Sj+i would act as a spring
force tending to push the nodes apart.
picture of the actual penalization process and neglects all
terms other than the penalty terms in the minimization
equation, it gives an idea of how the penalty functions affect
mesh motion and thus can help choose an appropriate form for
the penalty functions. It should be noted that the penalty
(II123)
Although the analysis above is a greatly simplified
terms can also keep the mesh motion regular before the minimum
35
internodal distance is approached. In the tests of the SMFE
method (discussed in chapter V), sometimes the initial
conditions cause singularities in the coeffecient matrix [A],
but the mesh moves smoothly even though the nodes are far
apart.
The penalty functions used in the SMFE code were taken
from an MFE package, DYLA [3]. They are
Ei =
1+CivlVql2<
c2 V
1 +
C3V+VqI2. . .(AXik)
n1
k2
n2+k3
(II124)
and
1+CsIvql2+
C2S
C3s+Vq .
1 +
k1
AXi~k
mi
ki
m2
J L(Axi~ki) J ,
(HI25)
where Axi=Xi~Xi_i,
and
Vq=
qjQil
XfXi1 .
(III26)
For the case of a system of n PDE's,
Vq 2=WT1Vqi  2+WT2{Vq2{2+...+WTN J Vqn J2
where WT ,... ,WTn are the same weights used to weight the
functionals in (III8). The remaining terms are all constants
with ki the minimum desired separation distance between nodes.
By changing the constants in the penalty functions, the
behavior of the penalty functions can be altered to fit a
specific problem. The last two bracketed terms in each of the
36
penalty functions above increase as Ax^ approaches the minimum
separation distance, k For Axi>>k, these terms are very
small. The choice of constants in these terms determine how
fast the penalty functions grow as the minimum separation
distance is approached. The first bracketed term in each of
the penalty functions acts as a gradient dependent penalty.
The C3 and C^y terms act to push nodes i and i1 apart when
there is a large property gradient ahead. This action pushes
nodes out of steep areas and into the transition areas at the
edge of steep gradients. Similarly, the C2V and C2S terms tend
to pull nodes into steep gradients. By adjusting constants in
the penalty functions, the relative influences of the spring
force and viscous force can also be adjusted.
THE SMFE METHOD APPLIED TO BURGERS EQUATION
A first test of the SMFE code was done for a single
PDE, Burgers equation,
ut+uux=Juxx (HI27)
where 1/R is a constant viscosity coefficient (R=Reynold's
number). The solution of Burgers equation developes a steep
shock very quickly. From the earlier general form of the SMFE
method PDE (III3),
qt+uqx+L(q)=0,
we have for this particular case q=u and L(q)=(1/R)uxx. This
37
results in the matrix equation (III16) for moving the mesh,
[A][Su]=[B],
and also discretized equations for
fl=(Su)ux+(1/R)uxx (III28)
dXi
*i_dt~~S
(III17)
In this example, a finite difference approximation is used for
(HI28),
3ui 1r32u
where
Caun ui+1~ui~1
L3xjixi+1xii ,
32ut 2 r3uk 9uk"1
T1 =
L J .
[i
9x2Ji (xi+1Xji)L9x 9x
and, as noted earlier in (III10),
]
9uk ui+1ui
(HI29)
(III30a)
(III30b)
9x xj+ixi .
Substituting (III29) into the residual, (III9) and expanding
the summation results in for interval k,
R,,(Sl.ul)[[H]f^(X)
1
+
R
(III31 a)
38
and for interval k1 ,
Rk1(Sl.1ui1)
[Â£].
9uk"1
+(Si~ui)
9x
9ui 9uk1
^(x)
ri 
L 9xJi 9x
(x)
1
+
R
[^]. ] tx) + [Jj] 1 (x)
92u
9xz i i1 L9x1i'l
We also have the derivative terms
92u
9x2
(I1131b)
9Rk~1
9Sj
9Rk
r 9u"] 9uk_1
L Jv I
9xJi 9x
(x)
r 9u 9u^
9Si [L9xJi*9x
<(>k(x)
(III32a)
(III32b)
(II31 a), (HI31 b), (III32a) and (III32b) can now be
substituted into the minimization equation (III15). To
evaluate the terms in (III15), it can be seen that the
following integrals of the shape functions are needed:
[ 1 + [k(x) ]2dx= [ l+t<))k+1 (x) ]2dx=X'1 + ""Xl (III33a)
Y.* Yj J
ri + 1
fX 1 , xi+1_xi
j <{>k(x)<))k+1 (x)dx=g
(III33b)
fXi+1
,1,, \ 92U . 1 r du1'
^(x)9Fdx=^
1r3uk 9uk_1
2L 9x 9x
]
(III33c)
and
fX*+1 k / y92u, 1j 9uk+1 9uki
j ^i + 1(x)9FdX=2t979F] .
J v .
(III33d)
The first two integrals, (III33a) and (III33b), are evaluated
in a straight forward manner using (III8a) and (III8b), but
39
the last two integrals are not quite as simple because of the
second derivative terms. Since 3uk/3x is discontinuous at the
endpoints of the interval 32uk/3x2 will not be defined, but
rounding u so that the first derivative discontinuities
disappear, we can use the approximation
32uk 1 r3uk 3uk~1 i
3x2 Xi+1Xi 3x 3x
in (III33c) and
32uk 1 r3uk+1 3uk
3x2 "xi + iXi* 3x 3x
in (III33d) and then perform the integration.
Substituting (IIX31) (III33) above into the
minimization equation (III15) results in
^i, i1 (Sii Ui_i )+Ai f i (Si~Ui)+Ai j i + i (S^ + i Ui+ )=Bj_
(HI34)
where
I X 1 X H** 1 r 3uI 3uk ^ r3uj 3uk^
l_ 6 L3xJi1 3x L3xJi 3x
(III35a)
'i.i
Xi1) r3uj 3uk_1 2 (xi+1xi) r 3u1 3uk
3 L3xJi 3x * ' 3 L3xJi 3x
(xi+iXj)
r 3u", 3uk
L3xJi~3x
r 3u 1 3uk
L3xJi+1 3x
(III35c)
Bi=^
1
=2R
3uk 3uk_1
3x 3x
1
r 3uI _3u^ 3uk~1
l3xJi 3x 3x
6R(xi'xil)
k(xi+1"xi)
r SuI _3uk_1
L rw J i
rilHi +2r
L3xzJi1 ^L3xzJi
32u
3xJi 3x
r1 +2[1
13xJi 3x '3x2ji+1 L3x2ji
(III35d)
40
and we now have 3 sets of N equations, (11134), (III29) and a
discretized form of (III17) to solve for the 3N unknowns, S^,
Ui and xi, for i=1,...,N.
THE SMFE METHOD APPLIED TO COMPRESSIBLE FLOW
The compressible flow equations in conservation form,
(III2a)(III2c), with an artificial viscosity term added to
each become
continuity: Pt+uPx+Pux_v1 Pxx=<3 (III~36a)
momentum: mt+umx+mux+pxv2nixx=0 (III36b)
energy: et+uex+eux+(pu)x\;3exx=0 (III36c)
where m=pu, e=total energy/mass and the vj are small constant
viscosity coefficients. All of the equations above are of the
form qt+uqx+L(q)=0. A state equation, P=P(p,m,e), is also
needed. Here we will assume a perfect gas, so the state
equation is
POf1) [e^] (III37)
where, for this problem, Y=1.4.
From (III18), we have
IWT) IContinuity+WT2lmomentum+WT3lenergy s0 we want to
calculate
31 3*1 312 3*3
' s7=WTl3Sr+WT23s7+WT33s7=0
(11138)
where each term 3In/3S^, n=1,2,3, is evaluated separately.
41
We will now evaluate each of the three terms in (III38).
The residual for each of the three PDE's is as in (III9),
i+1 9(q
Rk?= l [(qn)j + (uj~sj)~9xn~'k, q2k, q3k) n1,2,3.
j=i
(III39)
When considering the continuity equation, (III36a), we have
q=p and
L1(p,mJe) = p^vi^ f (III40a)
for the momentum equation, (III36b), q2=m and
L2(p,m,e)=m^+g^v2gj2 ^ (III40b)
and for the energy equation, (III36c), q3=e and
L3(p,m,e) = (e+P)^+u^\)3^f (III40c)
In the residual of (III39), we will use the same form of
finite difference approximation for Q as was used for Burgers
equation in (III29),
4r(sruj>[ff]j[L<
where [L(q)]j is a finite difference approximation for L(q) at
location xj. Derivative terms in [L(q)]j will be approximated
using (III30a) and (III30b). So, substituting (III41) into
the residual, (III39), and expanding the summation results in
for interval k,
Rk=
(Siui)
r 3qi 3qk
L3xJi_3x
"[L(q)]i
4>k(x)
(Si +  ui+1)
r 9qi 9qk
L3xJi+1~3x
[L(q)]i+1 ^i+1(x)+L(qk)
(III42a)
and for interval k1,
Rk"1 =
(Si1_ui1)
(SiUi)
[git,
r3qi 9qk 1
L9xJi 9x
It]*?
[L(q)].Uk1(x)+L(qk1)i
(III42b)
Since L(q) is not a function of S, we get the derivative terras
9Rk
9Sj=
r
L9xJi 9x
*{
(IIIJ3a)
and
9Rk~1
9S?
9qk1
[Ml _
L9xJi 9x
<()k1 (x) # (III43b)
Substituting (III42a), (III42b), (III43a) and (III43b) into
the minimization equation (II115), results in
0=
fxi + 1
J.
si^4[Â£]^]'lLt,,1*kx>
(Sl*1Uil) + Cx)*LCqk>
rMi _MH
L9xJi 9x
<})^(x)dx
+ J
fxi
xi1
(Si1_ui1)
(Siui) r 9q1 3qk1 L9xJi 9x a [L(q)]i (x)+L(qk~1 )
M J
r 9qi 9qk~T
L9xJi 9x
(J>k"1 (x)dx
(11144)
This gives us the contribution of each PDE to the total
minimization equation (III38). As noted earlier, the
minimization equation (III14) (and thus also (III38)) can be
written in the form of (III16), [A][Su]=[B]. Row i of the
coefficient matrix [A] will consist of the coefficients of
(Si_1Uii), (Si~ui) and (Sj+iUi+i) in (11144). So
Ai,i1J
f*i
xi~1
3q i 3gk~1
r1
L3xJi1
3x
r 3q  3qk ^
L3xJi 3x
^i1 ^i""^ ^x
or, after performing the integration,
Ai,i1=
Similarly,
(XiXi1)
[1 ~
L3xJi1
3qk
1
3x
j~3q j Sq^
3xJi 3x
Xi1
r3Â£ 3qk~1
l3xj 3x
Yi
[xi+1
+J
Ai
or, after integrating,
(xixi1)
rMi
L3xJi 3x
i.
r3q 3qk 1
l3xJi 3x
kkdx
2 (xi+ixi)
+ 
(III45a)
r 3qi 3qk
L3xJi 3x
(III45b)
and
A i+i = fXl+1 [f1 llf^l ^ld)k'1d>k'1dx
li + 1 I L3xJi + 1 3x L3xJi 3x ^i+ri
or, after integrating,
All of the remaining terms in (11144) go into [B], so
fxi+1
Bi="J.
[L(q)]i<>J.(x) + [L(q)]i+i
r 9qi 9q^
L9xJi 9x
*5(x)dx
rXi
j.
:ii
[L(q) ]i_i (x) + [L(q) (x)+L(qk 1 )
iBr
3q
k1
9xJi 9x
or, after integrating,
(x)dx
Bi=
r 3q i MU
L9xJi 9x
r j (xi + lXi)  (Xi + 1Xi)
[L(q)]i l +[L(q)]i+1 l
r 3qi
L9xJi 9x
r,, (xi~Xii) , (XiXiT)
LL(q) Ji13+lL(q)Jig
J Xj^
L *Jxl_1
(III45d)
(III45a), (III45b), (III45c) and (III45d) when evaluated
for each PDE and weighted appropriately give the contribution
of that PDE to the.matrices [A] and [B].
Obviously, the integrals involving L(qk) and L(qk1) in
(III45d) will have to be evaluated separately for each PDE.
Making use of the integrals of the shape functions (III33a)
(III33d) along with
fxi + 1 fxi+1. (xi+xi)
, k+1(x)dx= <>k(x)dx=
J V i * V i
(11147)
( we can evaluate the remaining integrals in (III45d) For the
45
continuity equation, where q=p and L(q) is as in (III40a),
over interval k1 from x^i and we get
 L(qk'1 )k1 (x)d x= [ p~v1 (x)dx
xi1
Xi1
fXi 3uk_1
J.
:i1
3^[PiI^lj (x) + Pi4>k1 ](j)k1 (x)dx
3uk1
3x
[Pi1
] vi 1(x)dx
' xi1
(Xj[Xil) (XiXjl)
+P1
]
^1r3p^ 3pk~11
2 3x 3x ,
and similarly, for interval k from x^to Xi+i,
(III47a)
For the momentum equation, where q=m and L(q) is as in (III
40b), over interval k1 from x^_i and x^ we get
rXi
j
L(q)4>k_1 (x)dx= [m^+^v2^]4)k 1(x)dx
Xi1 ^xi1
3uk1 fxi
=^j tmilil] (x)+mi4>k"1 (x) ]k_1 (x)dx
v . 
xi1 .
rXi
agifV1 <*><*,*{1 p*{i(x)qx
3x J ri
xi1
3uk~1
3x
[mii
xi1
(XiXi_l) (Xi"Xi_l)
+mi5
3pk1 (xjXj1) vgramk 3mk~1i
L J_2 L3x 3x 1
3x
(III48a)
46
and similarly for interval k from to Xi+i,
fxi+T1, s, 3ukr (xi + 1Xi) (xi+1xi)l
j L(q)k(x)dx= [mi+mi+1g]
For the energy equation, where q=e and L(q) is as in (III40c),
over interval k1 from x^_i and xj we get
J 1 Mq)+^1(x)dx.[ 1
xi1 xi1
=^TH 1 [i1 *11 (x)+eii_1 (x)]^1 (x)^x
xi1
+^TH 1 [pi1^ii (x)+pif_1 (x)]^1 (x)dx
xi1
+^TH 1 [uii<>Jl] (x)+Ui<)k"1 (x)]k_1 (x)dx
xi1
V3J 1 Â§F^1(x)dx
xi1
3uk1
[ei
(XiXi) (XiXil)
3x LJi1
.3uk~1
6 "+ei 3 J
(xiXi1) ^ (xiXiT)
6 +Pl 3
(xixi1) (xixi1)
3x
3Pk1
3x
[Pi
i1
[ui1
+Ui
V3r3eJl 3ek~1 
'2 L3x 3x J
(III49a)
and similarly for interval k from to xi+i ,
(xi+1~xi) (xi+1xi)
i+1 6
^3r3e^ 3ek~1 i
2 L3x 9x J
(IIIi9b)
Now, the contribution of the functional I for each PDE to [A]
and [B] can be evaluated separately with (III45)(III49),
then all contributions weighted appropriately and summed along
with the penalty terms to give the matrix equation
[A][Su]=[B]. We also have the finite difference approximation
(III41) for <5 for each variable q=p,m,e; the equation of state
(III37) for P; the relation u=m/p; and a discretized form of
S=dx/dt, giving seven equations at each mesh point to solve for
the seven unknowns p, m, e, P, S, x and u.
CHAPTER IV
SMFE CODE DEVELOPMENT
For a system of n PDE's, the SMFE method outlined in
the previous chapter results in n+1 equations to be solved at
each node, one for the mesh motion and one for each PDE. At
node i, for the mesh motion we have (1113*0,
^i, i1 (^ii ui1 )+^i, i (Sjuj )+Aj_, x + 1 (Sj+i ~Ui+i )=Bi
where Si=ii=dxj/dt as in (III17),
and for for each of the n PDE's we have a numerical
approximation of the form of CIII11),
9 Q *
<51 i(SiUi) [g7]i[Li (Q1 ,Q2. .Qn) ]i
9 Q
^ni = (Si~Ui) ]^ [L>n(Ql Q2 Qn)
where derivatives are evaluated as in (III30a) and (III30b).
The first test of the SMFE code was Burgers equation
where n=1, q=u and L(q)=(1/R)uxx. Initially, explicit
difference, equations were written for (11111) and (III17)
First the elements of [A] and [B] where calculated using xn and
un, then the system [A][Su]=[B] was solved for (Su)n+* where
the superscript n refers to the current time level and n+1 to
the next time level. Uj[ was advanced with the explicit finite
difference equation
n+1 r._ .n+1 (aunn r92uin, n
. . pn+1 n+1 .n+1
and, since =ui +(SjUi)
n+1
.n+1 n
AtSi +Xi
However,it was found that a very small maximum time
step (At=101*) was required for stability. This could be due
to the equations being stiff or a stability requirement of the
form Ax/At^u. This stability criterion requires that features
of the profile move no more than Ax in one time step At. When
a steady state profile was reached in the first test of Burgers
equation, the smallest value for Ax was about .0006 and the
front moved at a velocity of about .5, meaning the maximum
stable time step would be about At=.0012. The actual maximum
time step that could be used, however, was about half of this.
A stiff ODE has a component in its solution of the form e_k^
where k>>0. The solution to the first test of Burgers equation
is u(x,t)=1/[l+eR(*5x.25t)] for R=1000 [2], so there is a
stiff term in the solution.
To speed the solution process, an implicit ODE solver,
DASSL, was used. DASSL solves a system of m ODE's of form
f'l(k,yi,...,yui,7',...,
: (iv1)
^ra^yi >ym*7l > >^m)"0
given initial conditions
yi(to)=YOi ?i(to)=?Oi i=1,2,...,m
dv
where y=y(t) and .
50
The n+1 equations at each of N nodes resulting from the SMFE
method make up a system of the form of (IV1) where m=(n+1)N.
DASSL requires the user to write a main program (here
called SMFE) which calls the package of subroutines DASSL,
which in turn calls two user supplied subroutines, RES and JAC.
JAC calculates the elements of an mxm matrix of partial
derivatives the i,j element of which is given by
PDij =3fi/3yj+c3fwhere fi is as given by (IV1) and c is
supplied by DASSL. A DASSL option used here, however, allows
DASSL to approximate the elements of the matrix, so JAC in this
case is just a dummy routine. Subroutine RES computes the
elements of an array DELTA of the errors of approximation of
each equation. Given a finite difference approximation
Fi(t,yi,...,yN5h >?n)=0 of the 0DE
fi(t,y1 ,... ,yN.?i, .> then
DELTAi=Fi(t,y ,... ,yN,^i.Â£n) DASSL initially chooses an
approximate time increment At, approximates yn+1 by the
derivative of an interpolating polynomial of y at the last k
time steps (for a kth order method), then solves for yn+1 using
Newtons method. If this solution passes an error test with
tolerances prescribed by the user in the main program, the time
step is accepted. If the error test is not passed, the time
step is decreased and/or the order of the method is increased
51
and the process repeated: All program parameters including
initial conditions are initialized in the main program.
The largest part of the SMFE code is the subroutine
RES. For Burgers equation, the elements of the 1 xm array for
Yi are set to yi=xi, Y2=ui .,Y2n1=xn y2n=Un> where xj
and uj are the location and velocity of node j. Since Y2n1=xn
and i=dx/dt=S, then ?2n1=Sn the mesh velocity of node j. For
the compressible flow problem, yin3=xn> y4n2=Pn> y4n1=mn and
yjjn=en where xn, pn, mn and en are the values at node n. In
addition, u=m/p and, assuming ideal gas behavior,
P=(Y1)e1/2(m2/p), so u and P can be found in terms of the
y's.
At the beginning of subroutine RES, the parameters x,
S, u, p, j5, m, iR, e, Â§ and P (or x and u in the case of Burgers
equation) are assigned their appropriate values from the
elements of the arrays [y] and [Â£] at the next tentative time
level as provided by DASSL, then the array of errors of
approximation, DELTA, is computed and returned to DASSL. The
major part of this process is calculating the coefficients of
the mesh motion equation [A][Su]=[B], where for node i
DELTAn;L_g=Aj., i 1 (Si1 ui1 )+Ai i (Sj^Ui )+Aj. f i+1 (S^+i Ui+ )Bi
The other elements of DELTA for node i are
DELTAiji_2=Pi(SiUi) [jj^]^ + [LCOptinuity( P >m>e)
DELTAnj._i =lfli(Sj.iu) [ j + [kmomentum( P >me) ]^
52
r\ .
DELTAii=5i(Siui) [j^]^+[Lenergy(P >ra>e)
In the manner outlined above, the solution was advanced in
time
53
CHAPTER V
RESULTS
BURGERS EQUATION
Three test runs of Burgers equation were made. The
first and simplest (since no penalty terms were needed) was the
same one done by Dukowicz in his introduction of the SMFE
method [2]. The initial profile was a steep gradient between
x=.1 and x=,2 given by uj=l/2 [1+cos( 10ir(xiX2)) ] for 2 S i ^ 19
where X2=0.1. Boundary conditions were ui=1 at x[=0 and U20=0
at X2o=10 and a viscosity coefficient of 1/R=.001. The front
quickly steepened and formed a steady state profile of
thickness about 0.0015 as it propagated. The profile at the
initial time and three later times are shown in figure 1 where
+ shows the node locations and the solid line interpolates the
SMFE solution between the mesh points. As can be seen in the
figure, the nodes initially placed in the front moved with the
front, providing excellent resolution of the very steep
gradient with only 20 nodes and no oscillations developed.
Dukowicz did not publish numerical data, but these profiles
appear similar to his plots.
54
TIME = .000000
TIME = 4.000000
TIME = 8.000000
TIME = 12.000000
figure 5.1
Burgers equation I
55
The next two tests of Burgers equation were identical
to tests done by Gelinas and Doss using the MFE method in with
a package called DYLA [3]. In both of these cases, penalty
functions were required. Again a very small viscosity
coefficient, 1/R=.0001, was used. In the first of these tests,
a sinusoidal initial profile of 21 evenly spaced nodes was
used, with the initial condition u=sin(2iTx) +.5sin(irx) and
boundary conditions u=0 at x=0 and U2i=0 at X2i=1. A steep
front developed from this profile and propagated to the right
as shown in figure 2. Without penalties, a point was reached
as the front steepened where DASSL could not find a solution.
The penalty terms used here were the same ones used by Gelinas
and Doss: C y=C] s=C2V=c2S=0 C3V=C3S=1 k1 =5*10_5 (the minimun
node separation distance), ni=mi=2, n2=m2=1 k2=103, k^=0 and
kij=10^ (as in equations (III24), (III25)).
The last test for Burgers equation used a trapezoidal
initial profile with very steep sides as in figure 3* Here
penalty functions were needed because the initial profile
produced singularities in the [A] matrix. The penalty terms
were the same as in the case above. As the profile moved,
excellent resolution of the steep front and corners was shown
with only 31 nodes as in figure 3 (some internodal distances
were very small, so not all nodes are clearly visible in the
figure).
56
figure 5.2
Burgers equation II
57
TIME  .000000
TIME = .200000
TIME = .400000
TIME = 1.000000
figure 5.3
Burgers equation III
58
In both of the last two tests, the numerical data
compared very closely with the MFE solutions of Gelinas and
Doss. The mesh motion was very similar for both methods and
profiles of the solution appeared identical. In all three of
the tests for the SMFE method, steep gradients were resolved
very well by a relatively small number of nodes and no
discernable oscillations developed behind the gradients in any
of the cases. In the last two tests, the mesh points moved
into corners and steep fronts, a desirable characteristic for a
mesh motion scheme. Having mesh points within a steep gradient
prevents oscillations from occuring behind the gradient and
mesh points in corners are needed to provide good resolution of
such features.
COMPRESSIBLE FLOW
The shock tube test of the SMFE method was similar to a
test done with the MFE method by Gelinas and Doss [31.
Boundary coniditions were u=0, p=1 and e=2.5 at x^=0 and
UN=0> pn=125 and ejj=.25 at xjj=1 where N is the number of
nodes. The initial gradient at the membrane was linear, had a
thickness of .001, was located between x=.5 and x=.501 and had
40 evenly spaced nodes within it. The artificial viscosity
coefficients used were v =V2=V3=.0001. The same penalty terms
as used with Burgers equation were initially tried here, but
were found not to work well, so several tests were made with
59
various changes. The basic terms were C2V=C2S=0 C3v=C3s=1 >
kl=5x107 (the minimun node separation distance), n=mi=2,
n2=m2=1 k2=10"3, k3=0 and kn=106 (which are similar to those
for Burgers equation) while Cy, Cs and k3 were varied.
Repeating the expressions for the penalty functions from
earlier,
En =
1+C1vlVql2+
C2V .
Si=
1 +C1S jVqj2+
C3V+lVq'
C2S
1 +
k1
1 n,r
Ax^ki
n1
k2
(Ax^.ki )
h2+k3]
(III24)
C3S+IVql
,.jiL_pr ,kt. j
tolKlJ LcAxtkT) 2J
(HI25)
and it can be seen that Cv and Cg should increase the
"viscous" and "spring" forces in a large gradient, while k3
should provide a small viscosity term that is always present.
By setting certain of these terms to zero and giving the others
small values, the effect of each term or combination of terms
could be studied. There seemed to be little correlation
between the expected and actual results. At time = 0.1, all
choices of penalties caused 56 nodes to be present in the
shock and M5 in the contact surface. Possibly the magnitudes
of Civ, CS and k3 were not varied enough, so the mesh motion
wasn't affected to an appreciable degree. Usually the contact
surface and shock were well defined, but occassionally
oscillations developed behind them. In this respect, the best
60
results were obtained when nodes tended to cluster around
corners, but often it wasnt clear why some penalty terms
caused this to happen. The greatest problem was mesh tangling
and irregular motion in the vicinity of the expansion wave.
Of the combinations of penalty terms tried here, the
best results were obtained for Ciy=10"6, [<3=0, Cs=0 and the
rest of the penalty terms as noted earlier. Figures 4, 5 and 6
show profiles of the density, pressure and velocity. +
indicates the node locations and the solid line shows the
analytic solution as computed by the methods of chapter II.
The steady state profile showed the thickness of the contact
surface to be 0.015 and the thickness of the shock to be 0.001.
As can be seen in the figures, the mesh moved with the contact
surface and shock, providing good resolution of both, and no
oscillations of appreciable magnitude developed behind the
gradients. The SMFE solution also stayed very close to the
analytic solution. Some smearing of the contact surface and
expansion wave can be seen due to the artificial viscosity
terms, but since the shock normally tends to steepen with time,
it did not show this effect. When the initial gradient was
located farther to the right, the shock reached the end of the
tube and was reflected, as shown in figures 7, 8 and 9.
However, the penalty terms had to be changed slightly (k3 from
++++++++++
niiumniMMiitMHU
+
TIME = .000000 TIME = .040000
TIME = .100000
TIME = .160000
figure 5.4
Density
62
TIME = .000000
TIME = .040000
TIME = .100000
TIME = .160000
figure 5.5
Pressure
63
TIME = .000000
TIME = .040000
TIME = .100000
TIME = .160000
figure 5.6
Velocity
64
TIME = .000000
TIME = .100000
TIME = .120000
TIME = .140000
figure 5.7
Reflected shock, density
65
TIME = .000000
TIME = .100000
TIME = .120000
TIME  .140000
figure 5.8
Reflected shock, pressure
66
TIME  .000000
TIME = .100000
TIME* .120000
TIME = .140000
figure 5.9
Reflected shock, velocity
67
0 to 105) since the penalties used earlier would not allow the
shock to be reflected (some resolution of the expansion wave is
sacrificed by using the new penalties).
PENALTIES
The penalty terms are a very important part of the SMFE
method. Only in one case, the first test of Burgers equation,
were they unnecessary. The choice of penalty terms is also
very important, but it appears that to obtain good results
trial and error is about as effect in choosing the penalties as
choosing terms in a logical manner according to their
anticipated effect. A good choice of penalties must be made so
that mesh tangling is avoided and so that the mesh is allowed
to move in manner so as to prevent oscillations. Overall, it
seems that the penalty method used here to keep the mesh
regular was somewhat awkward.
OTHER COMMENTS
Overall, the SMFE code seems to work well in giving
excellent resolution of very steep gradients and sharp corners.
In the relatively simple case of Burgers equation, the method
was fairly easy to use. In the more complicated case of
compressible flow, however, some problems were encountered.
The penalty terms had to be chosen carefully and often
68
manipulated to get good output. Also, for Burgers equation,
little computational time was needed (on the UCD Prime 9950
system), but the more complicated compressible flow problem
took a much longer time, 215 hours (the choice of penalties
could greatly influence the run time). This may have been due
in part to the general nature of DASSL. Although DASSL was
very convenient and easy to use, an implicit ODE solver more
specific to this problem probably would have been faster. As
noted earlier, DASSL approximated the elements of the matrix of
partial derivatives in the subroutine JAC, the elements of
which are given by PDjj = 9fj/3yj+c3fwhere f^ is as given
by (IV1). The user had the option of writing the routine to
calculate these elements or letting DASSL approximate them, but
when f^ was one of the 80 equations in the system [A][Su]=[B],
he elements of [PD] were complicated and very difficult for the
user to compute. In the compressible flow case, this matrix
was a sparse 320x320 banded matrix, but DASSL would not treat
it as a banded matrix unless the routine was written by the
user. By not treating the matrix as banded, quite a bit of
unnecessary work was being done. The length of time needed for
a complicated onedimensional problem (such as with
compressible flow) to be run on the UCD Prime system, could
make it difficult to extend the method to two dimensions using
DASSL and the same computer. Overall, however, the main
69
drawback to the SMFE method was the necessary manipulation of
the penalty terms to get good output.
70
BIBLIOGRAPHY
1. Anderson, J.D., Modern Compressible Flow With Historical
Perspective, McGrawHill Book Company, New York, New York,
1982.
2. Dukowicz, J.K., "A Simplified Adaptive Mesh Technique
Derived from the Moving Finite Element Method", Journal of
Computational Physics, 56, 2, November, 1984.
3. Gelinas, R.J. and Doss, S.K., User Instruction Manual for
DYLA A Moving Finite Element Code in 1D, Science
Applications, Inc., 1811 Santa Rita Road, Pleasanton, CA,,
94566, 1981.
4. Gelinas, R.J., Doss, S.K. and Miller, K., "The Moving
Element Method: Applications to General Partial
Differential Equations with Multiple Large Gradients",
Journal of Computational Physics, 40, 202249 (1981).
5. Miller, K. and Miller, R.J., "Moving Finite Elements. I",
SIAM Journal of Numerical Analysis, 18, No. 6, December,
1981.
6. Miller, K., "Moving Finite Elements. II", SIAM Journal of
Numerical Analysis, 18, No. 6, December, 1981.
7. Richtmeyer, R.D. and Morton, K.W., Difference Methods for
Initial Value Problems, Second Edition, Interscience
Publishers, John Wiley and Sons, New York, New York, 1967.
8. Roache, P.J., Computational Fluid Dynamics, Hermosa
Publishers, P.0. Box 8172, Alburquerque, N.M., 87108,
1972.
9. Whitham, G.B., Linear and NonLinear Waves, Interscience
Publishers, John Wiley and Sons, New York, New York, 1974.
10. Zucrow, M.J. and Hoffman, J.D., Gas Dynamics, Vols. I
and II, John Wiley and Sons, Inc., New.York, New York,
1976, 1977.
