BIBLIOTHEEK TU Delft P 1955 6037
CONSIDERATIONS ON DIFFERENCE
APPROXIMATIONS
FOR FLOW PROBLEMS
PROEFSCHRIFT
TER VERKRIJGING VAN DE GRAAD VAN DOCTOR IN DE TECHNISCHE WETENSCHAPPEN AAN DE TECHNISCHE HOGESCHOOL DELFT, OP GEZAG VAN DE RECTOR MAGNIFICUS IR. H. R. VAN NAUTA LEMKE, HOOG-LERAAR IN DE AFDELING DER ELEKTROTECHNIEK, VOOR EEN COMMISSIE UIT DE SENAAT TE VERDEDIGEN
OP WOENSDAG 3 MAART 1971 TE 16 UUR DOOR
ARIE CORNELIS VLIEGENTHART
DOCTORANDUS WIS- EN NATUURKUNDE GEBOREN TE DORDRECHT
DRUKKERIJ HOOGLAND DELFT — 1971
Dit proefschrift is goedgekeurd door de promotor
PROF. DR. E. VAN SPIEGEL.
Aan mijn ouders
Aan Lienke
CONTENTS
INTRODUCTION T References 16 DISSIPATIVE DIFFERENCE SCHEMES FOR SHALLOW WATER EQUATIONS 19
1 Introduction 19 2 Preliminary considerations 19
3 Numerical damping and phase shift 21
h Finite-difference schemes 22
5 Shallow water flow 23
6 Numerical results 28
7 A local instability 2 9
8 Conclusions 31 References 32 THE SHUMAN FILTERING OPERATOR A N D T H E NUKERICAL COMPUTATION OF SHOCK
WAVES 33 1 Introduction 33
2 One-dimensional considerations 33 3 A two-dimensional example 35
References 1(0
ON FINITE-DIFFERENCE METHODS FOR T H E KORTEWEG-DE V R I E S EQUATION Ul
1 Introduction h^
2 General remarks h}
3 Exact solutions h2
h A centered finite-difference scheme it3
5 Dissipative difference schemes ^9
Appendix 31
References 58 SAMENVATTING 6l
INTRODUCTION
Calculations of solutions of partial differential equations which are
based on the replacement of derivatives by finite-difference approxi-mations have in fact been known for many years. For hyperbolic problems
Massau [1] (l899) introduced the so-called characteristic methods. In
addition, the work by Courant, Friedrichs and Lewy [2] (1928) has been
a basic contribution to the numerical solution of hyperbolic partial
differential equations. Courant, Friedrichs and Lewy were only interested in difference equations as a tool for proving the existence of solutions
of partial differential equations. Their ideas, i.e. the famous stability
condition, the energy inequalities, and the leap frog difference method
turned out to be basic when during the Second World War numerical solutions
of partial differential equations had to be provided for various kinds of problems. At present, industry demands more and more numerical solutions
to increasingly complicated problems, which demand has been stimulated by
the development of high-speed computers. This development of fast computers
has led to the solution of a wide variety of complex time-dependent flow
problems by using numerical techniques to solve the equations of motion for compressible and incompressible fluid dynamics. Various methods are available
for solving certain classes of problems, but each of these methods has certain
limitations. Needless to mention, one demands that a difference approximation
be consistenc, convergent, and stable. In practice, however, additional
re-quirements are made concerning the necessary human programming time, the machine computing time, and the computer memory capacity. Another requirement,
based on practical grounds, is that the method be flexible. For instance,
suppose that one has to design conditions which assure that a flow is
shock-less. For this purpose a number of calculations will be made, and, in some af
these shocks will certainly be present. It is necessary then that the shocks are well represented in order to get a good impression of what is happening.
In fact, one should require here that the method be capable of handling not
only shocks but also contact discontinuities and discontinuities in the
deri-vatives of the flow variables (rarefaction heads).
The two fundamental representations for the equations of fluid dynamics
are the Eulerian and the Lagrangian. These two representations of the
partial differential equations are formally equivalent, but this may not
be true for the solutions which are obtained from the appropriate dif-ference approximations. Each of the two representations has its advantages
and disadvantages. In the Lagrangian form the independent space variables
are referred to a reference frame which is fixed in the fluid and thus
moves with it. Consequently, a network as applied for the numerical
cal-culations may be distorted from its original form, which possibly leads to serious errors as well as severe computational instabilities. Generally
speaking, in the case of strong fluid distortions it will be doubtful
whether the numerical results represent the true solution to the problem.
To restore the original form of the meshes a new network can be introduced
if necessary, though this process is difficult due to interpolation pro-blems. Kolsky [3], Blair et al. [It], and Schultz [5] obtained good results
by doing this. Generally, the calculations in the Lagrangian representation
become too cumbersome if more than one space coordinate is involved. The
major advantage of the Lagrangian approach is the possibility to represent
interfaces between fluids of differing properties. In the case of one space coordinate, i.e. plane flow, flow with cylindrical symmetry, and flow with
spherical symmetry, the Lagrangian representation will often be preferred.
In the Eulerian form the reference frame is fixed relative to the observer
and the fluid moves through this spatial coordinate system with time. Here extreme fluid distortions can be handled, but material interfaces become
hard to identify as time goes on. Until about I96O, i.e. before the
develop-ment of the Lax-Wendroff second-order scheme [6] the only stable difference
schemes for the Eulerian representation were less accurate than the Lagrangian
ones. Therefore, the latter were preferred. The large amount of publications over the last ten years concerning approximations for the Eulerian formulations
shows the opposite to be true.
To eliminate the disadvantages of each of both representations, mixed
formulations have been suggested by Frank and Lazarus [7], and Noh [8]. In
this respect the Particle-in-Cell method for compressible fluids should also
be mentioned. This method has been developed by Harlow et al. [9].
Of the numerous methods for the numerical solution of flow problems,
we first discuss the so-called characteristic methods. Here the system of
equations is first expressed in cheiracteristic coordinates. Integration then takes place along these coordinates. In the case of one or two first-order
(one second-order) equations in two independent variables the favourable
properties of the characteristic methods are generally recognized.
Difficulties arise, however, in the treatment of interfaces between different
fluids, contact discontinuities, shock waves, and free surfaces. The positions of such boundaries in the solution are often not known in advance. Additional
information is generally required at these boundaries. In the case of a shock
this is given by the Rankine-Hugoniot conditions which represent the physical
consideration that an element of mass that crosses the shock must conserve its
mass, momentum, and energy (see also [lO] ).
A special problem arises in the characteristic methods if the computed values
of the dependent variables are to be found at some fixed time. For instance,
with two independent variables, two-dimensional interpolation is required.
This difficulty is avoided by Hartree's method [11] where a fixed network is
chosen in advance, and where the calculations progress over the constant time steps by means of linear interpolation. A similar method with simpler but less
accurate equations has been introduced by Courant, Isaacson and Rees [12].
Thomas [13] and Ansorge [lU] used Adam's procedures for the integration of
ordinary differential equations as an integral part of their characteristic
method. In the case of more than two independent variables good results were obtained by among others Bruhn and Haack [15], Tornig [l6],and Daubert and
Graffe [17].
Another class of numerical methods, in itself certainly not less
important, is not based on the reduction of the equations to characteristic
form. We first mention the well-known method by von Neumann and Richtmyer
[18]. To be able to treat shock waves, a pseudo-viscosity term was added to the pressure here. The unnecessary smearing out of rarefaction waves is
prevented if the viscosity term is replaced by the term as introduced by
Rosenbluth (unpublished work, Los Alamos). The von Neumann-Richtmyer method
which has been introduced for the Lagrangian representation and plane geometry
was also used in cylindrical and spherical symmetry by Erode [19, 20]. Another choice of the pseudo-viscosity term has been proposed by Landshoff [21], whilst
Ludford, Polachek and Seeger [22] have applied a viscosity term which is based
on the true viscosity.
Filler and Ludloff [23] first expressed the equations of viscous,
heat-conducting gas flow in divergence form and then constructed appropriate implicit and explicit difference schemes. Longley [2lt] similarly expressed
the equations in divergence form with the addition of explicit viscosity,
heat conduction, and mass diffusion terms. Various difference schemes with
the Landshoff viscosity term and the Lagrangian equations in conservative
form have been considered by Fromm [25]. Schultz [26] presented tensor artificial viscosity in the case of higher dimensions.
The Lelevier scheme [27] for the Eulerian equations of gas flow is
accurate to the first order. The same is true for the Lax scheme [28] which
has been devised for the equations in conservative form. This scheme has the
effect of introducing dissipation without the addition of explicit dissipative
terms. For cylindrical geometry Payne [29] has used the Lax scheme, though in this case the equations cannot be put into conservative form. For spherical
symmetry similar work has been done by Roberts [30].
A very successful method is the second-order accurate Lax-Wendroff scheme [6]
for both the Eulerian and Lagrangian equations in conservative form. The
oscillations which generally appear immediately behind the shock are typical
for this method. To reduce these oscillations. Lax and Wendroff added an artificial term to their difference scheme. It should be noted that this
differs, in principle, from the addition of the pseudo-viscosity term in the
von Neumann-Richtmyer method. Here dissipative effects were introduced at the
level of the partial differential equations. The Lax-Wendroff scheme is usually
complicated because of the explicit appearance of the Jacobian matrices of the system of partial differential equations. The two-step modification (see [31])
which also is accurate to the second order may therefore often be preferable.
The methods as obtained by Strang [32], Gary [33], and Gourlay and Morris [3'*]
should be considered as modifications and extensions of the original Lax-Wendroff
ideas.
The method which is called "shock fitting" has been applied by Thomas [I3],
Gary [35], and Quarles [36]. Here the shock is treated as an internal boundary. The necessary conditions at the shock and the motion of the shock front itself
are determined by means of the so-called Rankine-Hugoniot conditions. Difficulties
arise, however, if shocks develop spontaneously within the fluid, since it is not
known where shock fitting is to be applied. Shock fitting in two dimensions is
described by Richtmyer and Morton
This method has been constructed for the one-dimensional gas flow equations
in the Lagrangian representation.
The Particle-in-Cell method, which has been developed at Los Alamos by
Harlow and his associates [9], is not based on differencing the partial
differential equations, but on certain elementary physical models. This method
has been constructed for the calculation of compressible fluids and is a
combination of the Eulerian and Lagrangian method, as has already been mentioned. Though the Particle-in-Cell method is not very accurate, it has
turned out to be surprisingly versatile for, among others, multifluid problems
with large fluid distortions. Both memory storage requirements and computing
time can be reduced here by a different transport calculation which does not
require the use of particles. This has been realized in the Fluid-in-Cell method as introduced by Gentry et al. [39], which is most suitable for
one-material problems involving large fluid distortions. Another method which
has some features similar to the Particle-in-Cell and Fluid-in-Cell procedures
is the Particle-and-Force method as described by Daly et al. [Uo]. This method
has also been designed for compressible flow problems. However, a big and fast computer is needed since the calculations are time-consuming and require a
large amount of storage.
A similar method for time-dependent viscous incompressible flow of fluid with
a free surface in the Marker-and-Cell method introduced by Welch et al. [I4I].
Also for this method a big and fast computer should be available.
Finally we mention the Alternating Direction Implicit method as developed
by, among others, Mitchell and Fairweather [1*2], and D'YaJconov [1+3] and the Locally One-Dimensional method introduced by various Russian authors, in
particular D'Yakonov, Marchuk, Samarskii and Yanenko (see e.g. [I4I4]). These
methods may sometimes be favourable for reasons of economy.
As has already been remarked, a difference approximation should be consistent.
For instance, let us consider an initial-value problem for which the two-level
difference scheme IT = C(At)U has been so constructed that a unique vector-valued function U at the time t=(n+l)At is calculated from any sufficiently
smooth function Ü at t=nAt, by applying the difference operator C(At). The
order of accuracy of the difference operator is then defined as the largest
real number p for which
||u(x,t+At) - C(At)u(x,t)| I = 0 [(At)^"^^] , as At * 0,
Here u(x,t) and u(x,t+At) denote the genuine solutions of the partial
differential equations at the times t and t+At respectively. The norm
in this definition is the L operator norm.
The difference operator is called consistent if the order of accuracy
is positive or, as p will generally be an integer, if the approximation is
at least accurate to the first order. We point out that the refinement of
the mesh has to be set in advance, i.e. in the case of p space variables
x,,x„,...,x the functions g.(At) with Ax. = g.(At) (i=1,2 p) must
I <r p 1 1 1
be known beforehand. This phenomenon has, of course, no counterpart in
ordinary differential equations. The next examples may be illustrative
(see also Ames [1*5]). We discuss three difference schemes for the parabolic
(1) partial differential equation — - —— = 0
at 3x
u°*'
J"S:;
J = XU^^^ + ( 1 - 2 A ) U ^ + X U ^ _ ^ , - 2 ( U 1 / X ) U " * ^ + U " ! i = - { U > i - 2 ( 1 - U ^ - l = 2X ( U " , - U ° - ^ - U"^^ + U" , J J+1 J J J - 1 (2) } , (3) J J-1 (Crank-Nicholson (DU Fort-Frankel)where X = At/(Ax) , and U. = U(jAx, nAt); j and n are integers.
We substitute any exact solution u(x,t) of equation (l) into the three
difference schemes. By means of Taylor's series and using equation (l) itself, we then get for
u(t+At) - C(At)u(t)
At
(M
in the point (jAx, nAt) respectively,
2 1* 1.,/3 u vn ^ 1,, x2 ,3 u vn
. A t ( - ^ ) . . - ( A X ) ( - ^ ) .
^ ( A x ) ^ ( ^ )' - ^ ( A t ) ^ ( ^ - f ^
^ )'+... ,
^2 3x^ J ^ St^ 2 3^23^2 J
1,„. >2,a-^u ,n _^ 1,. ,2,3 u,n ,At ,2 ,3^u >n ^
^ 3t3 J ^2 3^1* J Ax 3^2 J
These expressions are called the local truncation errors and only the terms
which constitute the principal parts of these errors have been given
explicitly here. We notice that the truncation errors of the schemes (2)
and (3) tend to zero as At, Ax -*• 0, and these schemes are consistent
indeed. Scheme ('t), however, only is consistent with equation (1 ) if 4t goes faster to zero than Ax, i.e. if, for instance. At = constant (Ax)
holds good. If on the other hand we have At = vAx, where v is constant,
then scheme (h) is consistent with the hyperbolic partial differential equation
2 3 u 3 u . 3u
3t 3x'^ 3t
Let the solutions of the partial differential equations and of the
difference equations be denoted again by u(x,t) and U(x,t). The difference
method is called convergent if
||u(x,t) - U(x,t)|| ->• 0, as At -»• 0, holds good.
Here, as before, the space increments Ax. will tend to zero according to
the definite functions g.(At).
Suppose that the calculations are carried from t=0 to t=T. The difference
operator C(At) is called (Lax-Richtmyer) stable if for some T > 0 the set
n
[C(At)] for 0 < At < T
0 4 nAt 4 T
is uniformly bounded.
According to Lax's Equivalence Theorem (see also [1*1]), for a certain
class of problems convergence and stability are equivalent if the consistency
condition is satisfied.
Restricting ourselves for the moment to linear problems where Fourier analysis is applicable, the stability requirement involves the uniform boundedness of the
amplification matrices (see also [37] )• If the difference operator is
approxi-mated with m th order accuracy, then the same must be true for the amplification
matrix according to the theory. We illustrate this for the case of%the
one-dimensional hyperbolic partial differential equation
Suppose that we are dealing with the initial-value problem for this
equation and that the initial condition is such that the solution can be represented by a Fourier series,
u(x,t) = I , , u(k) exp[i(kx-lt) ] .
Substitution of the k th term of this series into equation (5) gives
for the exact amplification factor g of equation (5) the relationship
g(k,At) = e , with 1 = ak.
Let us consider the next difference schemes for equation (5):
(Lax)
U"*1=U" - (p-p)(u" - u" J , J J-p J-p J-p-1
with p 4 y 4 p +1 (Lelevier)
u"+i - u" = -fv(u"^^ - u"^J + u" - u" ) ,
J J ^ J+1 J - 1 J+1 j - 1 '( C r a n k - N i c h o l s o n )
u f ^ = u " - iy [ ( i - y ) ö " . i + 2 y ö " - ( i + y ) u " _ , ] .
J J J J J
(Lax-Wendroff)
where y = aAt/Ax and, as before, U. s U(jAx, nAt).
The amplification factors of these difference schemes turn out
to be in the same order
cos 5 - i ysin 5,
[ c o s ( p 5 ) - i sin(pC ) ] [ 1 - ( y - p ) ( 1 - c o s 5+ i s i n C ) ] ,
( 1 - i i p s i n 5 ) / ( l + li M s i n 5 ) ,
1 - i U s i n 5 - 2y s i n ^ ( 5 / 2 ) , where S = kAx.
For small values of 5 we readily obtain for the first two amplification factors
g(k,At) = l-iy? + 0(5^),
and for the last two amplification factors we get
g(k,At) = 1-iy5 - i ^'^ + 0(5^).
Consequently, we have here modulo terms of second and third order
respectively in 5 •
g(k,At) = e-"^« ,
which is the exact amplification factor of equation (5), since
lAt = akAt = yC .
These results are in accordance with the fact that the schemes (6)
and (7) are accurate to the first order, whilst the schemes (8) and (9)
are accurate to the second order. In our opinion, one could devise
difference schemes by simply basing the construction on requirements
concerning the amplification factor (matrix).
In hyperbolic problems the wave components in any solution should
be retained with unchanged amplitude. In some nonlinear cases, however,
a slight numerical dissipation is indispensable, which makes the difference scheme irreversible, i.e. it cannot be used to compute solutions both
forward and backward in time. A difference scheme of this kind has been
introduced by Friedrichs [1*6]. In the literature the term "overstability"
is sometimes used in this connection as such a scheme is in fact too
stable from a practical, computing point of view, which in certain cases may lead to a very poorly converging solution. The numerical dissipation
is generally introduced in a scheme by one-sided differences. The three
schemes (6), (7) and (9) are all irreversible. For the construction of a
nondissipative difference scheme it is sufficient that all differences
are correctly centered around the same point, just as in the case of scheme (8) both space and time differences are centered around the point
(jAx, (n+5)At) in the x-t plane.
There is another feature which is due to the discretization, but which generally cannot be eliminated completely by a special choice of
This is the phase error of the Fourier components in the solution, present
in both explicit and implicit methods, and which may be quite disadvantageous
for high wave-number components. For both amplitude and phase errors the paper of Roberts and Weiss [1*7] is enlightening. In this paper the phase
distortion of a second-order scheme is considerably reduced by modifying the
scheme in such a way that the scheme becomes accurate to fourth order in
space, but remains accurate to second order in time. A successful method
reducing the phase errors has been introduced by Fromm [1*8]. This method involves a linear combination of forward and backward time steps.
REFERENCES
[l] J. Massau, Memoire sur 1'integration graphique des equations aux
derivees partielles, F. Meyer-van Loo, Gent, Belgium (1899).
[2] R. Courant, K. Friedrichs and H. Lewy, Math. Ann., 100, 32 (1928).
[3] H. Kolsky, Report LA-1867, Los Alamos Scient. Lab. (1955).
[1*] A. Blair, W.C. Metropolis, J. von Neumann, A.H. Taub and M. Tsingow, Maths. Comput., 13, ll*5 (1959).
[5] W.D. Schultz, Methods in Comp. Phys., Vol. 3, p. 1, Acad. Press, New York (196I*).
[6] P.D. Lax and B. Wendroff, Comm. Pure Appl. Math., 13, 217 (196O).
[7] R.M. Frank and R.B. Lazarus, Methods in Comp. Phys., Vol. 3, p. 1*7, Acad. Press, New York (I96I*).
[8] W.F. Noh, Methods in Comp. Phys., Vol. 3, p. 117, Acad. Press, New York (1961*).
[9] F.H. Harlow, Methods in Comp. Phys., Vol. 3, p. 319, Acad. Press, New York (1961*).
[10] A. Jeffrey and T. Taniuti, Non-linear wave propagation, Acad. Press, New York (1961*).
[11] D.R. Hartree, Numerical Analysis, 2nd Edit., Oxford Univ. Press, London and New York (1958).
[12] R. Courant, E. Isaacson and M. Rees, Comm. Pure Appl. Math., 5, 2l*3 (1952).
[13] L.H. Thomas, Comm. Pure Appl. Math., 7, 195
[Il*] R. Ansorge, Numer. Math., 5, 1*1*3 (1963).
[15] G. Bruhn and W. Haack, Zamp, 9b, 173 (1958).
[16] W. Törnig, Numer. Math., 5, 353 (1963).
[17] A. Daubert and O. Graffe, La Houille Blanche, nr. 8, 81*7 (1967).
[18] J. von Neumann and R.D. Richtmyer, J. Appl. Phys., 21, 232 (1950).
[19] H.L. Brode, J. Appl. Phys., 26, 766 (1955).
[20] H.L. Brode, Research Memorandum RM - I82I* - AEC, Rand Corp. (1956).
[21] R. Landshoff, Report LA - 1930, Los Alamos Scient. Lab. (1955).
[22] G. Ludford, H. Polachek and R.J. Seeger, J. Appl. Phys., 2l*, 1*90 (1953).
[23] L. Filler and H.F. Ludloff, Maths. Comput., 15, 26l (I961).
[2l*] H.J. Longley, Report LA - 2379, Los Alamos Scient. Lab. (196O).
[25] J.E. Fromm, Report LA - 2535, Los Alamos Scient. Lab. (1961).
[26] W.D. Schultz, J. Math. Phys., 5, 133 (196I*).
[27] R. Lelevier, private communication (1953).
[28] P.D. Lax, Comm. Pure Appl. Math., 7, 159 (I95l*).
[29] R.B. Payne, J. Fluid Mech., 2, I85 (1957).
[30] L. Roberts, J. Math. Phys., 36, 329 (1958).
[31] R.D. Richtmyer, NCAR Techn. Notes 63-2, Colorado (1963).
[32] W.G. Strang, Num. Math., 6, 37 (196I*).
[33] J. Gary, Maths. Comput., 18, 1 (196I*).
[3I*] A.R. Gourlay and J. Morris, Maths. Comput., 2 2 , 28 (1968).
[35] J. Gary, Report NYO - 9603, Courant Inst. Math. Sci., New York Univ. (1962).
[36] D. Quarles, Ph.D. Thesis, New York Univ. (1961*).
[37] R.D. Richtmyer and K.W. Morton, Difference methods for initial-value problems, 2nd Edit., Wiley (Interscience), New York (1907).
[38] S.K. Godunov, Mat. Sb. , 1*7, 271 (1959).
[1*0] B.J. Daly, F.H. Harlow, J.E. Welch and E.E. Sanmann, Report LA - 311*1*, Los Alamos Scient. Lab. (1965).
[1*1] J.E. Welch, F.H. Harlow, J.P. Shannon and B.J. Daly,
Report LA - 3l*25, Los Alamos Scient. Lab. (1965).
[1*2] A.R. Mitchell and G. Fairweather, Numer. Math., 6, 285 (1961*).
[1*3] Ye.G. D'Yakonov, Z. Vycisl. Mat. i Mat. Fiz., 3, 385 (1963).
[1+1*] A.A. Samarskii, Z. Vycisl. Mat. i Mat. Fiz., h, 21 (196I*). [1*5] W.F. Ames, Numerical methods for partial differential equations,
Nelson, London
(1969)-[1+6] K.O. Friedrichs, Comm. Pure Appl. Math., 7, 3l*5 (l95l*).
[1*7] K. Roberts and N. Weiss, Maths. Comput., 20, 272 (I966).
[1*8] J.E. Fromm, J. Comp. Phys., 3, 176 (1968).
Dissipative Difference Schemes for Shallow Water Equations
S U M M A R Y
Three dissipative finite-difference schemes are discussed for the numerical calculation of discontinuous shallow water flow. The shallow water equations have been derived on assumptions which are not acceptable in the case of disconti-nuous flow. However, they may give satisfactory results if only weak jumps are present.
It will be shown that a coarse network may affect the velocity of propagation of the computed bore.
1. Introduction
The nonlinear shallow water equations derived by Stoker [1 ] are frequently used for continuous solutions to problems concerning shallow fluid flow.
In this paper, we discuss three dissipative finite-difference schemes for approximating dis-continuous shallow fluid flow governed by the shallow water equations.
The schemes are tested for stability, numerical damping and phase shift. The Courant-Friedrichs-Lewy stability condition appears, though the schemes are dissipative if this condi-tion is slightly strengthened.
In [2] Leendertse had to deal with numerical damping as a falsification of the solutions of his calculations of long-period water waves. This damping is necessary if we are to obtain a stable solution in our calculations.
As in [2] the numerical phase shift influences the solution, but is minimized by a refinement of the superimposed grid.
By a nonlinear combination of the shallow water equations, we obtain a second system of equations. The two systems are equivalent in the case of smooth solutions, but yield different solutions if a discontinuity occurs in the flow. We shall arrive at a system which is physically acceptable for our computations.
The results obtained from the jump conditions at the discontinuity in the How are favourably compared with the numerical solutions. Results from laboratory experiments carried out by Cavaillé [3] give an indication of the usefulness of the shallow water equations as a computa-tional model for discontinuous shallow fluid flow.
2. Preliminary Considerations
We consider the initial-value problem for the class of equations of the form
",+ƒ.. = 0 , O^t^T, and u{x,0) = (p(x), (2.1) where u(x, t) is a p-component vector and/(u) is a known vector-valued function of u.
Suppose that equation (2.1) can be transformed into
u, + A{u)u, = 0, (2.2) where A (u) denotes the Jacobian matrix of/with respect to u. Let the quasi-linear system (2.2)
be hyperbolic, i.e., the matrix A has p real and different eigenvalues a, for all values of the ar-gument u.
We approximate the system (2.2) by means of explicit difference schemes of the form
v(x,t + At) = Sj,v{x,t), (2.3) in which At represents the time increment and S^^, a finite sum of translation operators with
matrix coefficients,
(J)
The operator T denotes translations by amounts Ax = At/X in the x-direction, A being constant independent of Ax and At, T'v(x, t) = v{x+jAx, t).
Definition. The differential equations (2.2) are approximated with m-th order accuracy by means of the difference scheme (2.3) if for all sufficiently smooth, genuine solutions u(x, t) of (2.2),
\\u{x,t + At)-S^Mx,t)\\ = 0{At"'*^) as At^O holds uniformly in t.
The norm appearing in this definition denotes the Lj norm. We assume that the calculations will be performed from t = 0 to f = T. For this time range we require the difference scheme (2.3) to be stable, ignoring the influence of the boundary conditions.
Definition. The difference scheme (2.3) is (Lax-Richtmyer) stable if for some a >0,a constant M exists such that
WSUèM
for all n. At satisfying 0 < At< a, and 0 ^ nAt ^ T.
In the nonlinear case we consider the first variation of the operator Sj„ and then replacing the variable coefficients by their value at some given point, we obtain the "localized" operator. For stability of Sj^ in the nonlinear case we require the associated "localized" operator to be stable at every point. Suppose that u(x, 0) has period q in x, then the local amplification matrix G of S^jj will be defined by
G(x,zlt,0 = I C , e ' J ^
(j)
where I, = kAx; k is dual to x in the Fourier transform and ranges over the values 2iirlq, r being any integer.
According to Richtmyer and Morton [4] the aforementioned stability condition involves the uniform boundedness of the matrices G" and yields the von Neumann necessary condition for the eigenvalues 3^ of G:
\gXx,At,(,)\^\ + 0(At), 0 < 4 t < a , v = l , 2 , . . . , p . (2.4) In [4] the following is proved:
Theorem. If the amplification matrix G has a complete set of linearly independent eigenvectors, then condition (2.4) is sufficient as well as necessary for stability if a constant x exists such that zl ^T >0, where A^ is the Gram determinant of the normalized eigenvectors.
In the next sections we shall discuss dissipative difference schemes.
Definition. The difference scheme (2.3) is dissipative of order 2r, where r is a positive integer, if there is a positive constant 5, such that
\gXx,At,i)\^l-dm'', v = l , 2 , . . . , p
for all X and all (, and At satisfying \i\^n, and At < a, where a is some positive constant number.
3. Numerical Damping and Phase Shift
The definition of the amplification matrix G in the preceding section presupposes a component u(/c)e''* of the Fourier series for i^(x, 0). This component changes for the genuine solution u(x, () of system (2.2) with constant matrix A into
ü(/c)e'"'^-"', (3.1) l/k = a,, (3.2) where a, denotes an eigenvalue of A.
Relation (3.2) can be found by inserting (3.1) into system (2.2). We assume that k and / are real. If the system (2.2) is replaced by the difference scheme (2.3), the frequency / is generally chang-ed into a different, and possibly complex frequency which we denote by /'.
Consequently, we have in case of equation (2.3) for the amplitude and velocity of the Fourier component, respectively,
«(fe)e'"'<''" and Re(/')/fe •
The concept of a complex propagation factor T has been introduced by Leendertse [2],
gi(*x-I'r)
We follow the propagation factor over a time interval in which the component with frequency / propagates over its wavelength. Then, for t = 27r//,
T(x, At, i) = exp[27ii{l - (I'/l)}] . (3.3) The modulus and argument of T represent measures for the numerical damping and phase
shift respectively.
Using (3.2) and (3.3), we may derive
\T(x, At, 0\ = exp[27tlm(r)/(ka,)] (3.4) and
arg[T(x, At, i)] = - 27r[Re(/')/()taJ- 1] . (3.5) The velocity ratio Q will be given by Q=Re(/')/(Jlca,). (3.6)
We substitute the component ü(k)é""~'''^ into equation (2.3) with constant coefficients; then it can be shown that
[G(x, At, £,)-e-"'"I'\ü(k) = 0 ,
where I is the unit matrix, and G is the amplification matrix of the difference scheme again. Thus we find a relation for l':
g,(x. At, i) = &-""• ; (3.7) here g^ denotes an eigenvalue of G again.
Remembering the definition of a dissipative difference scheme in section 2, it is clear from (3.7), that for such a scheme we have Im(/') < 0, for all ^ T^O. Hence dissipation goes together with damping of the Fourier components.
4. Finite-Difference Schemes
We shall use the notations:
u(jA.x, nAt) = u"j = u",
Du" = u ; ; + i - t / ; _ , ,
D,u'' = u"j^,-Wj, D,u" = u ; - u ; , , ,
Mu" =(u%,+u';.^)/2
M,u" = (u%,+u'J);2,
M,u" = (u"j + u"j.,)/2.
Hereafter we shall denote also the tiumericdl solution of (2.2) by u. The following three dissipative difference schemes for the system (2 2) will be discussed.
Scheme A u"+' = Mu"-A/2D/", Scheme B «"*' = u''-A/2Df" + /.^/2[(M,A")(DJ")-(M2A-)(DJ")] , Scheme C tt"+' = Mu"-A/2D/" u"*2 = u"-;.D/"*' .
Scheme A is a first-order-accuratc difference scheme given by Lax [5].
In [6] Lax and Wendroff introduced scheme B This scheme has second-order accuracy. Scheme C is the two-step Lax-Wendroff procedure found by Richtmyer [7]. If the coetficients are constant, it changes over to scheme B, though with mesh spacings 2Ax and 2At.
The local amplification matrix G of scheme A, as derived in [8], is
ƒ cos i-iXA sin 4 , (4.1) in case of scheme B matrix G is
/ - a / l s i n ( ^ - / l M ' ( l - c o s i J ) , (4.2) whilst in case of scheme C matrix G is
/-i/l/4 sin2c-/lM'^(l-cos2(J). (4.3) Denote Xa^ by n; then it is clear that the eigenvalues g of the amplification matrices (4.1),
(4.2) and (4.3) are, in the same order,
cos ^-ifi sin i , (4.4) l-i^ism^-H^\-cosi), (4.5) 1-i> sin 2i;-/^^(l-cos 2c). (4.6)
The amplification matrices of the schemes above are all three polynomials in the matrix A, hence their eigenvectors are the same as those of A, constituting a complete set of linearly in-dependent eigenvectors, not depending on X and ^. Therefore, according to the theorem of section 2, the von Neumann condition is sufficient as well as necessary for stability. It is easily established that applying this condition to the eigenvalues (4.4), (4.5) and (4.6) yields for the three schemes the well-known necessary stability condition by Courant, Friedrichs and Lewy
[9]:
IMI g 1 • (4.7)
From (4.4) we obtain for the eigenvalues g„ of matrix (4.1):
\qJ^=\-(l-U^)smU- (4.8) According to the definition in section 2 we conclude that scheme A is dissipative of second
order, if \n\< 1.
From (4.5) we obtain in case of scheme B
l9vl' = l - 4 / / ^ ( l - / / ^ ) s i n * K - (4.9) whilst (4 6) yields for scheme C
ISJ^ = 1-4^^(1-^^)sin* <J. (4.10) Relations (4.9) and (4.10) show that the schemes B and C are dissipative of fourth order,
provided that |/j| < 1.
Relations (3.4) and (3.6) together with relation (3.7) and each of the expressions (4.4) through (4.6) can be combined into the following relations for
scheme A
\T\ = (cos^ ^ + n^ sin' {)"•"«' (4.11) Q =arctan(/itan(S)/(/j^), (4.12) scheme B
I ri = [ 1 - 4/i'(l - /i') sin* i <ï]"'"'<' (4.13)
scheme C
i n = [ l - 4 / / ' ( l - / i ' ) s i n * (?]"<'"<' (4.15) II sin 2£,
, 1 - 2 / i ' s i n ' O
The phase shift is related to the velocity ratio by
arg(T) = - 2 7 t ( Q - l ) . (4.17) Figures 1 through 6 show sets of curves of relations (4.11) through (4.17) for different values
of n; the wavelength is denoted by L, i.e., L/Ax = 2n/^.
In section 3 it has been argued that dissipation goes together with damping of the Fourier components. One can see from Figures 1, 3 and 5 that especially the short-wavelength compo-nents are attenuated (see also [4]).
Figures 2, 4 and 6 show clearly that the schemes A, B and C are dispersive. Again especially the short-wave components are affected.
In section 6 we shall apply scheme B to the equations of the shaflow water theory while providing at once a confirmation of the results above.
5. ShaUow Water Flow
The following experiments have been carried out by Cavaillé [3].
A basin with a horizontal, rectangular and smooth bottom is divided into two parts by a slide. The basin is filled with water; we denote the depths of the fluid on either side of the slide by hg and hi, with ho >hi.
The slide is rapidly pulled up and after a short time interval a discontinuity in the elevation of the fluid (bore) is observed, propagating with side 1 as the front side. (Hereafter side 0 will be called the back of the bore).
In the other direction a rarefaction wave propagates, resulting in a decreasing elevation.
Q = a r c t a n ( , T" 2 J I (^nO • (4.16)
u- 09- 08-ar 06- 05-04. 03- 02-av
oo-?
1
/—^'-\ /
\ y
l\\,,
-^^
1.0 / / / M / " ^ 0.7/ 0 3 /^
^
/ ^
^ /
// /
/ / Q l /^
2 3 4 5 7 10 15 20 50 100 • L A x Fig 1 Scheme A : Numericol damping os given by (4.11).r - 46S0
50 100
rig,2. Schemt A: Vtlocrty rotio ond ptW5« shrft os given
by (t.12)and(i:.17).
2 3 4 5 7 10 15 20 50 100
Fig 3. Scheme BrNunnericat domping os given by(t.l3).
org(T) Wtgj. 72
288
2 3 4 5 7 10 IS 20 50 100
rig. 4. Scheme B: Velocity rotio and phase shift os given by(AU]ond(417).
i 5 7 ld 1^ 20 5b 100 Fig 5 Scheme C- Numerical domping as given by(&-15).
28S
360
Fig.6 SchemeC:Vtlocitynitloand phase shift os given by(tl6)an(l(t.17).
TABLE 1 depths at 1 = 0 ho 1 1 1 I Al 0.696 0.348 0.174 0.043 laboratory results
«
0.96 0.93 090 098 numerical results system (5.4) 9 0.96 0.91 0.87 0.84 ht 0.84 0.63 0.51 0.39":
0.17 041 0,57 0 75 system (5.5) q /ij 096 0.94 096 1.07 0.84 0.63 0,48 0.29 <••! 0.17 0.42 061 092Dimensionless values of the velocity of propagation of the bore, denoted q, measured for different ratios of fig and hi, have been hsted in Table 1. For dimensionless variables see below in this section.
We shall compare the experimental results by Cavaillé with solutions of the shallow water theory (see Stoker [1]), determined by the equation of conservation of mass
M w = o.
„.„
and the equation of motion in the ;c-direction
dv dv dh „ , ,
Vt^"Tx^'ê-x = '- (")
The space and time coordinates are denoted by x and t; h denotes the depth of the fluid, V the velocity in the (horizontal) x-direction and g the acceleration due to gravity.
Multiplying equation (5.1) by i; and (5.2) by h, and adding the equations together, yields the equation of conservation of momentum
| ( H . | ^ ( / ^ . - ^ ) + . ^ | = 0. (5.3)
This equation can also be derived by integration of (5.2) with respect to the vertical spacecoordinate and making use of (5.1).
By combining equations (5.1) and (5.2) we also get
| [ i ( / , r ^ + S / , ^ ) ] + ±l^v(hv'+gh')-\+\gj-^(h'v) = 0,
expressing the conservation of energy, i.e. the sum of kinetic and potential energy. We introduce dimensionless variables:
X = x/ho v = v/^/gho I = ty/gjhn g = glg=\ : fi = hlho
ho denotes again the depth on the upstream side before the slide is pulled up.
Using the dimensionless variables and leaving out the bars, we may write the system which consists of the equations (5.1) and (5.2) in the notation (2.1)
rH + ^fi?\) = 0. (5.4)
Similarly equations (5.1) and (5.3) yield the system
— I \ + -:r\ •>„ , , , ) = 0 ;(» is called the discharge, (B = fti'. (5.5) St\(p) dx\(p^lh + \h^)
It will be clear that the systems (5.4) and (5.5) are equivalent in the case of a continuous solu-tion. Consequently, the eigenvalues of the matrix A are for both systems
a„ = v±^, (5.6) representing the speeds of progressive and retrogressive waves.
In the case of continuous solutions one would prefer system (5.4) for calculations; then, it follows from above that momentum and energy are conserved.
We shall show that in the case of a discontinuous solution the two systems are not equivalent at all.
6. Numerical Results
The calculations described below have been carried out with difference scheme B. Suppose the slide in the basin (see section 5) is pulled up at t = 0; thus, we take as initial con-ditions v = 0, ho=l, and, for example, hi=0.5.
Figure 7 shows a computed bore extended over about three intervals Ax. In section 4 we
F1g.7 Initiol depths and bore computed by opptying scheme G to system
|5.5);ix.a2,At.a04. exact solution
distance covered
.1/Ax 1 2 3 4 i 6 7 8 9 10 Fig 8 Distance covered by computed bore versus V^^iht-t , h,«0.5.
system j S i l system (5.5)
described a falsification of the phases of the Fourier components for the numerical solution in the case of too large meshes Ax (see also Figure 4). Though a linear analysis is not appropriate here, it may be possible that a too coarse network imposed on the x, t-space, causes a retarda-tion of the computed bore.
To show this numerically, we carried out the calculations with difierent values of Ax, whilst At is such that the ratio X = At/Ax keeps the constant value 0.2 in order to satisfy the stability condition
X(\v\+^h)èl,
which is obtained from relations (4.7) and (5.6).
The distance covered by the computed before versus l/Ax is shown in Figure 8 for the time intervals Og t^ 5 and 5 g t ^ 10.
The x-value for which the depth of the fluid equals the average value of the depths on both sides of the bore will also represent the position of the bore; it is determined by linear interpola-tion.
Figure 8 shows that the distance covered during the time interval 0 ^ t g 5 keeps increasing
28
au-4.8 IS I*a
for the considered values of Ax, whilst the distance covered during the time interval 5 g I g 10 is approximately constant for Ax < 0.25.
This can be explained by the following argument.
At the instant t = 0 a genuine discontinuity is present. The computed bore presents itself not just as a discontinuity, however, but as a strong gradient of depth and velocity of the fluid (see
also Figure 7). Thus, it may be expected that the Fourier series for the computed bore contains less short-wave components than the Fourier series for the discontinuity at ( = 0, and from sec-tion 4 it is clear that the retardasec-tion of the bore occurs especially at the start of the computasec-tion.
In the case of a stronger bore we obtained the same results: the graphs in Figure 9 show about
601 58 56+ I I ( 42 40 distance covtrtd _s«t«io "(Swf 5«t>10 -l/Ax 1 2 3 4 5 6 7
Fig,9 Similar to figure fl, except that h, =0,02.
10
the same feature as those in Figure 8, though in the case of system (5.5) the bore velocity is significantly larger than in the case of system (5.4).
We compared our numerical results with the experimental results obtained by Cavaillé. To minimize the retardation described above, we used
Ax = 0.1 and At = 0.02 for 0 < t < 5 , and Ax = 0.2 and At = 0.04 for t^5 .
Values of the computed bore velocity q, together with values of depth and velocity of the fluid behind the bore, /ij and V2, have been listed in Table 1.
In addition, analytic calculations were carried out using the generalized Rankine-Hugoniot relations (see the next section) applied to the systems (5.4) and (5.5). The results agreed well with the numerical resuhs (see [8]).
7. A Local Instability
In the case of system (5.5) and for /ij =0.174 and h, =0.043 an instability was observed. This occurred at the position of the jump, and directly after t = 0.
Figure 10 shows the depth/i as a function of x for various numbers of time-steps. For t >0.8 the depth becomes even negative locally, whilst the velocity increases more and more.
depth h
FigtO. System(5,S),scheme B: ALocol instobility;h,c1, h,.0,174, Ax.a25i&t.Q0S
By usmg different mesh sizes, we may conclude that the instability is independent of the fineness of the net and also of the ratio X = At/Ax. This local and nonlinear instability has been found by Burstein [10] during numerical calculation of supersonic flow in three independent variables. Our argument, as does Burstein's, is based on a local and partial absence of numerical damping.
For the eigenvalue t' — .Jh of the matrix A we have on the front side of the bore v-^= - ^ < 0 .
The generalized Rankine-Hugoniot relation (see also Lax [11]) for equation (2.1) takes the form
<?[«] = [ ƒ ] • (7.1) where [ ] denotes the jump across the discontinuity, and q is again the velocity of propagation
of the discontinuity.
We apply relation (7.1) to system (5.4) and, taking into account that the velocity i; on the front-side of the bore equals zero, it can be shown that behind the bore we have
V2 = (h2-hi)
y hi + h^'
and thus also
,-^h = (h2-hi)^j^-^^ - ^ 2 = .
where N^ is definite-positive.
The right member of this relation, and thereby the eigenvalue v-yjh is positive for h2/h, > 4.56. (In fact also for /12//11 < 0.44, but this is not taken into consideration here).
Similarly, applying relation (7.1) to system (5.5) yields with Uj = 0 :
and behind the bore,
„ rh-(h h\ M^^ fh - h\-h\h2-^hihi+h\
V-Jh-(h2-h,)^^j^-^h2- . where Nj is definite-positive.
Therefore, in this case, the eigenvalue 11 — ^/1 is positive for ftz/^i >3.21. It follows that for both systems (5.4) and (5.5), and for a strong jump, i.e.. a large ratio /i2/'ii, the sign of one of the eigenvalues of the matrix A changes across the jump. We may expect that this eigenvalue will be zero somewhere in the region of the computed bore, and considering relation (4.9) leads to the conclusion that here the difference scheme is not dissipative: no damping of the short-wave components is performed by the eigenvalue v — yfh.
Hence dissipation is absent just where it is necessary for a stable solution. The remedy to prevent the instability in our calculations is to give, as initial condition, the jump in h over two intervals Ax, instead of over one interval zlx (see Figure 11).
I I
!^ I I " r ! Si-,*' IN
J L
Fig,n Initial conditions to prevent nonlinear instability.
8. Conclusions
The equivalence of the systems (5.4) and (5.5) in the case of smooth solutions has been mentioned in section 5. From the results in section 6 it follows, however, that the discontinuous solutions of the two systems are different.
Lax [11] pointed out that in hydrodynamics the system consisting of the equations of con-servation of mass, momentum, and energy on one hand, and the system of the equations of conservation of mass, momentum, and entropy on the other hand, are equivalent in cjse of smooth solutions, but are not equivalent in the discontinuous case.
We consider in the shallow fluid flow a column at the position of the jump, as shown in Figure 12.
(2)
L
^ .^x
Rg.12. Cross-section view with column at the position of the Jump
For this column the conservation of momentum is expressed by
hAv. -•])'-'i-hAv2-9)v2+ighl- ïghl ==0. (8.1) The first two terms on the left side of (8.1) repiesent the change in momentum of the column
per unit time, whilst the force acting on the end sections of the column is represented by the last two terms. Equation (8.1) can also be derived by applying the generalized Rankine-Hugo-niot relation (7.1) to equation (5.3).
On the contrary, application of (7.1) to equation (5.2) yields a relation which cannot be mier-preted as a conservation law.
As it can also be established that the Rankine-Hugoniot relation of equation (5.1) expresses the conservation of mass across the jump, we may conclude that the conservation of mass and momentum is exactly satisfied by the jump relations of system (5.5).
Remembering the theory of weak solutions (see Lax [11]), it can be said that the physically correct jump relations belong to the piecewise continuous weak solutions of (5.5).
Therefore we prefer (5.5) for the calculation of a discontinuous solution in shallow fluid theory.
In [12] the problem of calculating a bore from the shallow water equations has been solved by using the von Neumann-Richtmyer method (see [4]); here however, use has been made of system (5.4), which is not acceptable from the physical point of view.
The solutions of system (5.5) are, as can be seen in Table 1, not very satisfying, especially if the jump grows stronger. This can be explained by the argument that the shallow water equa-tions are not appropriate to this problem: they have been derived on the assumption that the pressure is given as in hydrostatics (see Stoker [1]). This assumption, however, is completely incorrect if a discontinuity is present in the flow.
Nevertheless, as Table 1 shows, a weak bore is governed rather well by the shallow water equations.
Acknowledgement
The author wishes to express his gratitude to Prof. dr. E. van Spiegel. His helpful comments and advice made this work possible.
R E F E R E N C E S
[1] J. J. Stoker; Water waves, Interscience Publishers, New York, 1957.
[2] J. J. Leendertse; Aspects of a computational model for long-period water-wave propagation. Doctoral Thesis, Delft University of Technology, 1967.
[3] J. Cavaillé; Contribution a l'étude de l'écoulement variable accompagnant la vidange brusque d"une retenue,
Publ. Ministère de I'Air, No. 410, 1965.
[4] R. D. Richtmyer and K. W. Morton; Difference methods for initial-value problems, 2nd edition, Interscience Pub-lishers, New York, 1967.
[5] P. D. Lax; Weak solutions of nonlinear hyperbolic equations and their numerical computation, Comm. Pure Appl.
Math., Vol. 7, pp. 159-193, 1954.
[6] P. D. Lax and B. Wendroff; Systems of conservation laws, Comm. Pure Appl. Math. Vol. 13, pp. 217-237, 1960. [7] R. D. Richtmyer; A survey of difference methods for non-steady fluid dynamics, N.C.A.R. Tech. Note 63-2,1963. [8] A. C. Vliegenthart; Berekeningen van discontinuïteiten in ondiep water, Report NA-4, Math. Inst., Delft
Univer-sity of Technology, 1968.
[9] R. Courant, K. O. Friedrichs and H. Lewy; On the partial difference equations of math, physics, N.Y.U. Courant
Inst. Math. Sci. Res. Dept., N.Y.O.-16S9, 1956.
[10] S. Z. Burstein; Finite-difference calculations for hydrodynamic flows containing discontinuities, J. Comp.
Physics, Vol. 1, pp. 198-222, 1967.
[11] P. D. Lax; Hyperbolic systems of conservation laws, Comm. Pure Appl. Math., Vol. 10, pp. 537-566, 1957. [12] A. Preissmann and J. A. Cunge; Calcul du mascaret sur machine électronique. La Houille Blanche, No. 5,1961.
The Shuman FUtering Operator and the Numerical Computation of
Shock Waves
S U M M A R Y
Shuman's method of filtering short-wave components is discussed. In certain cases, this method can be used to prevent nonlinear instabilities in numerical calculations. As illustrative example, the computation of a detached shock in front of a blunt body is considered.
1. Introduction
For numerical work in weather prediction Shuman [1] devised a method of filtering short-wave components.
The application of dissipative finite-difference schemes for the numerical computation of flowfields containing shock waves, yields a damping of short Fourier components. However, the damping is absent for selective components of the solution. These encountered oscillations, which have no physical meaning, can be suppressed by applying Shuman's method. In certain cases this even appears to remove nonlinear instabilities as with the computation of a detached shock in front of a blunt body. This problem has been solved numerically by Bohachevsky and Rubin [2], Burstein [3] and Lapidus [4].
Bohachevsky and Rubin used Lax's scheme which has first-order accuracy and generally does not show nonlinear instabilities. However, it yields a shock which is spread out over a large number of mesh spacings. The scheme obtained by Lax and Wendroff has been applied by Buistein. To prevent nonlinear instabilities, an artificial viscosity term is added here.
Lapidus described the smoothing of short components by means of a finite-difference smooth-ing operator which is not an integral part of the difference approximation.
The complicacy of the second-order accurate methods used by Burstein and Lapidus makes them laborious to program. Our objective was to devise a method which has the simplicity of Lax's method, but yields a shock which is considerably narrower.
2. One-Dimensional Considerations
Consider the system of equations of the form
":+fx=0, - 0 0 < X < -1-00 , t g O , (2.1)
where u is a p-component vector function of x and t, and ƒ (u) is a known vector function of u. After differentiation we obtain the quasi-linear system u,-{-A(u)u^ = 0, with .4 = grad/. The system will be hyperbolic, i.e., the matrix A has p real and different eigenvalues for all values of u.
The numerical integration of system (2.1) will be carried out by means of explicit difference schemes. To this end an orthogonal mesh system with constant mesh spacings Ax and At is imposed on the half plane —oo<x< -l-oo, t^O.
A difference approximation to system (2.1) which has been introduced by Richtmyer [5] is
«r'=iK--.+«i- .)-mf'j.. -f'j-1) ,, ,>
where u(jAx, lAt) = u'j, and X = At/Ax.
The provisional values at the time t = (l-\-1) At are obtained from the first step, whilst the second step yields the final values at t = (l-i-2)At. It is weU known that this scheme is consistent with system (2.1) and has a truncation error which is 0[ (zlt)^] in the smooth part of the solution. We suppose that /. = 0(1).
If H(X, t) is periodic in x, then substitution of the component fl(„e"'* of the Fourier series for u(x, r) yields in the case of a constant matrix A
fl!,V = G(Jr,9)fl!„,
with
G(At,q)= I-iXAsin(2^)-X^A^\-cos(2i)} . (2.3) For the derivation of relation (2.3) see [6].
G is called the amplification matrix, I is the unit matrix, whilst ^ = qAx.
If the eigenvalues of matrix A are denoted by a„ then the eigenvalues g„ of matrix G are gAAt,q)= l-i,i sm(2i)-n' {\-cos(2i)} ,
where /x = /la„ v = l, 2,..., p.
In [7] we obtained the necessary and sufficient stability condition for scheme (2.2):
and derived
\gj'=l-4n^(l-n')sm*i. (2.4)
From the last relation we concluded that Richtmyer's scheme (2.2) is dissipative of fourth order, provided that \n\< 1. For definitions of stability and dissipation see [8].
Scheme (2.2) has no dissipation in three cases (see (2.4)); (1) l/'l = i ;
(2)/' = 0 ; (3) sin <^ = 0 .
If the wavelength is denoted by L, then we may conclude from the relation L/Ax = 2n/^, that case (3) is present if, for instance, L = 2^x, representing the length of the longest wave which is not damped by our difference approximation.
I>0
A
"\i
t - s V j\A/vy
...I
0 S 10 15 20 »• XFigure 1. Shock calculation with equation u, + u,=0 and Richtmyer's scheme; Ax=O.M, J( = 0.02.
Figure 1 shows the initial values and solutions at the times f = 5 and r= 10 of the equation M,-l-u^=0, which is a special form of system (2.1). The amplified short-wave pulses can clearly be observed.
These oscillations can be damped by the following smoothing element described by Shuman [1]:
Substitution of the Fourier component fl(,|e''* into (2.5) yields
ö;i;'=i(i+cos<ï)fl!,V. (2.6) The multiplicative factor of (2.6) is real, so that the phases of the components are unaffected.
However damping occurs for the selected components with about cos <^= — 1, or, for instance, again L = 2Ax.
Because of the linearity of (2.5) the factor j(l -(-cos {) is independent of ƒ
It is clear that the mean value of a field of infinite extent is not affected by (2.5) so that the conservation (see [8]) is still present.
By means of a Taylor expansion about x=jAx we obtain from (2.5): «• = «+i(4x)^ii„+0[(/lx)*].
Thus the element (2.5) yields dissipation to our difference approximation.
t.Q t-S \ t-IO \
Figure 2. Similar to Fig. 1, except that smoothing has been applied. . exact solution.
Figure 2 is similar to Figure 1, except that the operator (2.5) was applied after every two steps of scheme (2.2).
In the nonlinear case not all oscillations may be damped out completely, as the example in the next section shows.
The three-step scheme consisting of the combined equations (2.2) and (2.5) has first-order accuracy.
3. A Two-Dimensional Example
3.1. Preliminary
The two-dimensional version of Richtmyer's scheme for the hyperbolic system of partial differential equations w,-\-f^-i-h^ = 0, is
Ax'
At At
„,l + 2 _ „ l _ lfl+1 ft+l \ 11,' + ' _ ! • ' + ' ^
' • ' j . * - '*;.* 7 7 u j +1 .t - y J - 1 .tl - - r - 1 " j.k +1 " j.t - 1 J Ax Ay- (3.1)
Here the notation w'ji, = w(jAx, kAy, I At) has been used.
Suppose that the matrices A and B denote the Jacobians of/(w) and h(w) respectively. Substitution of the Fourier component *,,,) e'"* *'"'' into difference scheme (3.1) with constant matrices A and B yields
<t,^ = G(zlf,q,r) *{,.„, where
G(At,q,r)= ƒ -1'(cos a-I-cos P) — (A sin a-l-B sin ^ ) - 2 | — (/t sin a-(-B sin/?)[^; Ax I Ax
a = qAx, P = rAy .
G is called the amplification matrix again.
We apply scheme (3.1) to the two-dimensional equations of a polytropic gas, dp dm dn
T^ + ^ + ^ = 0 dt ox dy
dm d f3—v m^ ,, , / ^ , n^ \] d /mn\
6n d (mn \ 8 ( 3 - y n^ ,, . / ^ , m^W
dE 8 {i—y m , , ,, mE] 8 (1—y n , , ,, nE) „ ,,.,, p is the density, £ is the energy, i.e., the sum of internal and kinetic energy, and y is the ratio of specific heats, m and n denote pu and pt; respectively, where u and v are the horizontal and vertical flow speed. The pressure p has been eliminated by means of the relation
p = ( 7 - l ) { £ - ( m ^ + n^)/(2p)}.
The four equations (3.2) represent the conservation of mass, momentum (two components) and energy, in that order. By means of a similarity transformation Richtmeyer and Morton [8] derived in the case Ax = Ay for the eigenvalues g„ of matrix G
I^J^ = (1 -2n^,)^ + fi^(cos a-l-cos P)^ (v = 1, 2, 3, 4) where
H, = (At/Ax) Tsin^o + s m ^ I'.'.'^A (3-3)
and u' = {u sin a-i-1' sin ^)/.ysin^a-l-sin^)S; c is the local sound speed. The condition (At/Ax)-(c-\- s/u'' + v^)< 1/V2 appears to be necessary for stability.
From relation (3.3) it is clear that the difference scheme (3.1) has no dissipation if at least one of the n„ equals zero and in certain cases the difference scheme may be unstable then. This will occur with the calculation of a detached shock before a rectangular body moving with constant supersonic speed along its axis of symmetry. This stationary problem will be solved by means of a converging process using the non-stationary equations (3.2).
For this purpose we introduce dimensionless variables through the free-stream values of the density and sound speed, denoted by p ^ and c„ in the upstream region before the shock:
P = P/Pco. m = m/(p„cj, n = n/(p^c^), E = E/(p^ci).
Substitution of these new variables together with x=x/L, y=y/L and l = tc^/L, where L
denotes some representative length, into the equations (3.2) does not change these equations. Below we shall omit the bars.
Because of the symmetry of the problem only the upper half of the body with its surrounding flow is taken into account.
We shaU use constant values of Ax, Ay and At, whilst Ax = Ay. 3.2. Initial Conditions
Figure 3 shows a discontinuity that was initially prescribed to facilitate the process at the early stages. At the downstream side of the discontinuity the values of the variables p. m, n and £ are calculated by using the Rankine-Hugoniot conditions. For the lower part, where the dis-continuity is perpendicular to the axis of symmetry, these conditions are given by [ / ] = ö , where [ ] denotes the jump across the discontinuity.
For the upper part the conditions are given by [ ƒ ] — [A] = 0. The dimensionless free-stream values are
Px = 1. "loo = -Moo. "oo = 0, £a, = Wi + :;x3p;;
M is the local Mach number.
Therefore, together with the geometry of the body, only the values of M^ and y are to be specified to determine the problem uniquely.
3.3. Boundary Conditions
Except for the left boundary, the boundary conditions are imposed by means of rows of virtual points (see Fig. 4).
Figure 3. Initial discontinuity before the body.
Figure 4. Boundaries with lines of virtual points.
At the axis of symmetry (y = 0) we take as conditions at every time t=(/-i- \)At:
Pj.- • Pj.i " j . i ^ j . - i - ^J.'
Similar conditions are prescribed on the upper boundary of the body, and, interchanging the conditions for m and n, on the front side of the body.
On the upstream boundary of the flow region (x = 0) the variables will have the free-stream values.
On the downstream boundary and the upper boundary of the flow region we impose the conditions, representing the continuity of the first derivatives, in directions which are situated between the directions of the characteristics in the steady case (see Fig. 6). In this way, the influence of the condition which has been specified in some point of the boundary is directed into the domain of dependence of this point. On the downstream boundary {x = JAx) we thus have at every time t = (l-t-\)At:
PJ^ 2pj.k-pj~
and the same condition for m, n and £.
On the upper boundary (y = KAy) we have ^
Pj+ i.K+ \ — ^Pi.K~ Pj- l.K-1 >
and the same condition for m, n and £.
The upper boundary will be taken so far from the body that the lower characteristics pro^ ceeding from this boundary behind the shock do not intersect or touch the sonic line (see Fig. 6). In the stationary case, which is approximated here, this region is subsonic and a disturbance in any of its points would influence the whole region.
3.4. Numerical Results
The computation performed by applying scheme (3.1) to the system (3.2), together with the initial and boundary conditions described above, yields unstable results at the early stages. We observed an unlimited growth of short-wave components near the axis of symmetry: within the shock and at the body, and also at the top of the body near the comer.
Around the stagnation point we have uvO and t'aO, and thus M'*0, whilst near the top of the body, i.e., at the lower part of the sonic line, u'iic, and thus u'-^c.
It follows that in both cases relation (3.3) yields \g^\x\ and no damping is present locally. This is especially true for the wave components for which is valid cos a-Hcos /? = 0.
The following finite-difference smoothing operator appears to stabilize the solution
Figure 5. Densities a l o n g the c o o r d i n a t e ƒ = 30 at the times t = 0,50 At, 100 At 700 At. M „ = 4 , 3 , y
At = 0.2.
--l,4;Ax = Ay=l,
''J.k i(w'jXl, + w'jll, + w
+ 2
j.k + 1 ' V k - (3.4)
Substitution of the Fourier component »i'(,.r)e'''''*'^'' into (3.4) yields K*') = i(cos a-l-cos li)w\l,] ,
with, as before, a = q /dx and ^ = r Ay. Again the multiplicative factor is real and damping occurs for the components with cos a-l-cos ^ = 0 .
Taylor expansions of equation (3.4) about the point (jAx, kAy) yields K = w + k(Ax)'w,, + i(Ay)'w„.-\-Ol{Ax)\ (Ay)*] ,
and thus dissipation is introduced again.
Scheme (3.4) is applied after the two steps of equation (3.1) at every time t = (l-i-2)At and a stable and convergent process is obtained (see Fig. 5).
The positions of shock and sonic line, together with the characteristics in two points of the upper and right boundaries are shown in Figure 6.
Densities along some lines y = constant have been plotted in Figure 7. Because of the smooth-ing element, the accuracy of our approximation is of first order, i.e. we have to deal with the
Figure 6. Shock, sonic line and tharacteristics in two boundary points (t = 900 At).
Figure 7. Density profiles along lines parallel to the axis of symmetry {f = 900^().
accuracy of Lax's scheme. This scheme consists of the first step of scheme (3.1), and applying Lax's scheme to this problem, which needs approximately the same computation time as (3.1) and (3.4) together, yields a shock which is spread out over about ten meshes, as against three or four meshes in the case of our approximation.
The smoothing described by Lapidus [4] does not completely stabilize our calculations. R E F E R E N C E S
[1] F. G. Shuman, Numerical methods in weather prediction: II. Smoothing and filtering. Monthly Weather Review, 85, 11 (1957)357-361.
[2] I. O. Bohachevsky and E. L. Rubin, A direct method for computation of nonequilibrium flows with detached shock waves, AIAA Journal, 4, 4 (1966) 600-607,
[3] S. Z. Burstein, Finite-difference calculations for hydrodynamic flows containing discontinuities, J. Comp. Physics,
2 (1967) 198-222.
[4] A. Lapidus, A detached shock calculation by second-order finite differences, J. Comp. Physics, 2 (1967) 154-177. [5] R. D. Richtmyer, A survey of difference methods for non-steady fluid dynamics, NCAR Tech. Note 63-2, (1963). [6] A. C. Vliegenthart, Berekeningen van discontinuïteiten in ondiep water. Report NA-4, Math. Inst., Delft University
of Technology (1968).
[7] A, C. Vliegenthart, Dissipative difference schemes for shallow water equations, J. Eng. Maths., 3, 2 (1969). [8] R. D. Richtmyer and K. W. Morton, Difference methods for initial-value problems, second edition, Interscience
Publishers, New York, (1967).
