- Email: [email protected]
Numerical stability in wave propagation problems
compums
& stmtures,
Vol. 1, pp. 495-510. paRanloa Plwa 1971. Printed in orut
NUMERICS
Britain
STABILITY IN WAVE PROPAGATION PROBLEMS YUKSELUCKAN~
and ALFRERO H.-S. ANG~
University of Illinois, Urbana, Illinois 61801, U.S.A. Abstract-Stability characteristics for some finite-difference schemes obtained on the basis of a mathematically consistent lumped-parameter model and a suitable numerical integrator for wave propagation in spherically symmetric and axisymmetric ~nt~uo~ media are investigated. The numerical integrator used is Newmark’s /l-method. It is shown that the difference schemes for large radii which are identical to those for one-dimensional and plane strain wave motions are conditionally stable for values of B between 0 and l/6 and unconditionally stable for fi=1/4 and l/2. For points near the center or axis of symmetry, it is found that for all values of @the difference schemes are unstable. However, conditional stability may be obtained if one uses a dissipative ~h~isrn in the form of linear artificial viscosity; the effect of this dissipative mechanism on the stability characteristics of the considered difference schemes is also examined in detail. Conditionally stable difference schemes are also found by introducing appropriate dependent variable transformations.
Cd velocity of dilatational waves I imaginary constant J?’ t/At = index for time level r/Ar=index for distance of a generic point from the center of symmetry z/At=index for z-ordinate of a discrete point art&ial pseudo-viscous stresses in r- and z-directions, respectively distance of a point from the center of symmetry time independent variable displacement of a mass point in r-direction velocity of a mass point in r-direction acceleration of a mass point in r-dir&ion ; a factor in Newmark% integrator coefficient of linear artificial viscosity Ar, A& : tlnite space increments in r-, 8- and z-directions, respectively At time increment P mass density at, a,, ur direct stresses in r- and z-directions, respectively t eigenvalues of the amplification matrix n P 4 9% 4s r t u u ..
XNTRODUC3ION A WIDELY used method of solution for wave propagation problems in solid media is to discretize the continuous problem using a mathematical or physical model and to integrate the resulting system of differential-difference equations by the aid of a suitable numerical
t Present address: Middle East Technical University, Department of Computer Science, Ankara, Turkey; formerly research assistant in Civil Engineering, Univ. of Illinois, Urbana, Illinois. $ Professor of Civil Engineering. 495
496
YUKSELUCKANWK~
ALFREDOH.-S.ANG
integrator. An important matter associated with this approach is the question of convergence of the numerical solution to the continuous solution of the problem, which, in turn, is closely related to the stability properties of the finite-difference scheme through Lax’s Equivalence Theorem [4]. Unless the stability of a difference scheme ~orres~nd~ng to a properly posed initial value problem is ascertained, the resulting numerical solution may be far from representing the true solution. It is the principal objective of this study to analyze the stability characteristics of certain finite-difference schemes commonly used in wave propagation problems in higher dimensional solid media. Practical stability criteria are derived and presented in a form suitable for direct applications. The difference schemes are based on the assumptions of linearly elastic material behavior and small displacement; however, this does not unduly restrict the applicability of the derived stability criteria. For example, for a class of inelastic solids where the plastic velocity of propagation is smaller than the velocity of the elastic dilatational waves, the results remain applicable. Wave propagation in spherically symmetric and axisymmetric media are considered. As one-dimensional and plane strain wave propagation problems are, respectively, the limiting cases of the spherically symmetric and axisymmetric wave motions for points infinitely distant from the center of symmetry, they naturally fall within the scope of the present study. For all problems, the space variables are discretized using a discrete lumpedparameter model proposed by Ang [If; with this model, the resulting differential-di~erence equations invariably turn out to be the central finite-difference analogues of the corresponding continuous formulations, and hence are mathematically consistent. Time-wise discretization is accomplished with the aid of Newmark’s P-integrator [3]. The simple cases of the one-dimensional and plane strain wave propagation problems are of constant-coefficient type, for which the Fourier transform method [4] is suitable in deriving necessary conditions for stability. With spherically symmetric and axisymmetric wave motions, which involve variable-coefficient systems, the Fourier transform method yields only local stability conditions. The center of symmetry in a spherically symmetric medium and the points on the axis of symmetry in an axisymmetric medium are points of discontinuity at which the equations of motion are not valid. In order to render the displacements finite at such points of discontinuity, it is assumed that the displacements (and hence, the velocities and accelerations) of the points at the center or axis of symmetry vanish for all time. Furthermore, for the purpose of stability analysis, the boundary conditions, if any, are assumed to be linear and homogeneous. It is found that although the discrete formulations for one-dimensional and plane strain wave motions are marginally stable, the difference schemes for wave propagation in spherically symmetric and axisymmetric media are always unstable. Two different methods are suggested to remedy this situation: (1) The first approach, due to von Neumann and Richtmyer 121,uses a damping mechanism in the form of linear artificial viscosity in the original formulations. It is found that the modified difference schemes become conditionally stable, although the solution based on the damped system is invariably slightly different from that of the original problem. The amount of damping that produces a maximum stability range is investigated herein. (2) In the second approach, the original difference scheme is replaced with an equivalent formulation through appropriate transformation of the original variables, such that the resulting scheme becomes conditionally stable.
NumericalStability in Wave Propagation Probkms DI!XREXX
FORMULA’ITONS
AND STABILITY
497
ANALYSES
Spherically symmetric wave motion A displacement type approach for wave propagation in a sphe~~~y symmetric elastic continuum formulated on the basis of a mathematically consistent lard-p~ameter model and Newmark’s integrator yields the following difference equation [6] :
+ (I- 2BMZ
-u;,+:)+j?#+j
r&J;+2
+(1-2/3)U79+$J,
(1)
where u:zu(pAr,
nAt),
(2)
is the displacement in the r-direction, (r, 6, cp) being the orthogonal spherical coordinate system; /I is a numerical constant in Newmark’s integrator which can assume values of 0, l/12, l/8, l/6, l/4, or l/2; At is the uniform time increment, Ar is the uniform space mesh length, cd represents the velocity of propagation of dilatational waves in elastic continua, and &L--- i-l ‘Ar 2’
i=l,2,...,2P-t1
is)
is an index denoting the distance from the center of symmetry, and n = t/At indicates the time level at which computations are being performed. The domain of integration is given as
O
(4)
where B_2K2sm2kAr z 2+
K2 2K2. 2----sm-, P P
kAr 2
(5)
in which
and k is the frequency of the Fourier components of the error term, and 4 represents the eigenvalues of the amplification matrix.
YUKSELUCKANand ALEREDO H.-S. ANG
498
For infinitely large radii, asp+ co, the difference equation, equation (l), reduces to that corresponding to one-dimensional wave propagation, and B, as given by equation (5), becomes real-valued. The solution of the characteristic equation yields the following necessary condition for stability : 1 Ar -J(l-48)
c,bt<
if P<$,
Unconditionally
stable if p 2 $ ,
(7)
which ensures the satisfaction of the von Neumann necessary condition on the eigenvalues, Irl< 1. For points near the center of symmetry, the analysis of local stability is done by solving the quadratic equation (4) and computing the moduli 151for fixed values of p, In this case B is a complex quantity; consequently, the moduli It;] may be computed using the complex arithmetic feature of a digital computer. Figures 1 and 2 illustrate typical plots of the variations of I<]vs. Kf or various values ofp and /?, in which only the maximum moduli are shown. It is concluded that; (1) for all finite values of p, and all geometrically possible values of fi, the von Neumann condition, \
o-4
06
0.8
I.0
I.2
I.4
I.6
I.8
2.0
Value of K FIG. 1.
Spherically symmetric wave propagation. Effect ofp on
I 02
I 04
08
I
I
I
0.8
I.0
I.2
I I.4
stability. B=O.
I I.6
I I.8
Volue of K FIG. 2. Spherically symmetric wave propagation.
Effect of j3 on stability. p=l.
I 2.0
499
Numerical Stability in Wave Propagation Problems Axisymmetric
wave motion
On the same basis of formulation as above, wave pro~gation in ~sy~et~c media in cylindrical coordinates yields the foIlowing difference equations [6]:
elastic
nt2 ~~i-1,jtlffbr_+~,j-l~+~1~28~~U~++~,jf~~~~++~,j-~ n+ 1
-Ui-1,j+1+U;_+,‘ji-1 fUI-i,j-t)]P
)+B(Ut;+i,
j+l-@+I,
j-l-"l-l,
j-i-1
(9)
YUKSELUCKAN and ALFREW H.-S. ANG
500
in which, ~7,jZ u(pAr, qbz, nAl) , WY,j= w(pAr, qdz, nbt) ,
W)
are the displacement components in the I- and z-directions, respectively, with (r, 8, z) being the orthogonal cylindrical coordinate system; Ar and AZ are the uniform space mesh lengths in r- and z-directions, respectively, and cd and c, are the velocities of propagation of the dilatational and shear waves in an elastic medium, respectively. p and q are indices denoting the position of a generic point and defined as i-l
r
p’L\r=27 i-l
z
i=1,2 ).‘.)
2P+1,
j=l,2,...,2Q+l.
q=rz=25
(11)
The domain of integration is defined by O<~S PAr, O
and,
(12) where
&-$& 2
12c,2AtZ. k,Ar -sm, pAr2 2
(13
u2At2 . k,Ar 2Lsm2 c pAr2
- c~)Af2sin~ +$cdz pArAz
1 .
2
Numerical Stability in Wave Propagation Problems
501
For infinitely large radii; i.e. P-P co, the difference system, equations (8) and (9), reduce to those for plane strain wave motions, and B, and B, of equations (13) and (14) become real-valued. For this case, it can be proved [6] that the von Neumann condition is satisfied if 1 AZ -&2-8/I)
*<
for B
Unconditionally stable for /I 2 $,
(15)
where AZ< Ar. Equation (15), therefore, is a necessary condition for stability for plane strain wave propagation. For points at finite distances from the axis of symmetry, the difference scheme is of variable-coefficient type, and the analysis of local stability requires the determination of Iti from the characteristic equations, equation (12). The evaluation of the moduli is more complicated in this case; it can be shown [6] that ItI depends on the Poisson ratio, V, of the material as well. In Figs. 3 and 4, typical variations of the maximum moduli against K, [see equation (611for various values ofp, /?, and v are shown. The following conclusions can be drawn from these results: (1) For all values of p, B and v, the difference scheme for axisymmetric wave motion is locally unstable; (2) as p, or as /3 increases, the stability properties of the difference scheme improve; and (3) as v increases, the stability characteristics deteriorate, and the worst situation occurs for v = 0.50 (an incompressible medium).
I
050
I
I
1
I
I
02
04
0%
0.8
I.0
I I.2
I I.4
I I.8
I I.6
2.0
valye of K Fro. 3. Axisymmetric wave propagation. E&d of p on stability. p=O, v=O-00.
0’50’
I 0.2
I
o-4
I 06
I 08
I
I
I
B = 1/2f
IO 1.2 I.4 1.6 Value of K FIG. 4. Axisymmetric wave propagation. E%ct of fi on stability. p=l,
1.8
2.0
v=O-OO.
502
YUKBL UCKAN and ALFREDO H.-S. ANG
REMEDIES FOR NUMERICAL INSTABILITY Effect of linear artificial viscosity on stability
If a damping mechanism is appropriately introduced in an otherwise instable difference scheme, it may be possible to obtain a conditionally stable system. The linear artificial viscosity suggested by von Neumann and Richtmyer [2] is such a mechanism. In spherically symmetric wave motion, the damping term is introduced through the equation of motion in terms of stresses:
.. U!$ ++ - 06)+ h_ ar-PU,
(16)
where p is the mass density of the medium, a, and creare direct stresses, and qr is an artificial pseudo-viscous stress in the r-direction defined by qr = pl’c,ArE Zr’
(17)
in which I? is a constant referred to as the ‘coefficient of artificial viscosity’. The first two terms on the left hand-side of equation (16) belong to the equation of motion in terms of stresses from which the difference system, equation (I), is derived. aq,/& is discretized using the backward difference formula for the first-order time derivative? ti and the central difference formula for the second-order space derivative in conformity with the way the original scheme is discretized. The resulting difference scheme is a four level difference equation, whose characteristic equation is :
(18) where a, b, c and dare complex quantities given by a=1 +/3B+j?A, b= -2+(1-2/3)B+(l-3/3)A,
c= 1 +j?B-(1-3/Q/t, d= -PA,
(19)
A = 2 J(2)IYK sin”?,
(20)
in which
and B is the same quantity defined by equation (5). Equation (18) is a cubic with complex coefficients; this can be solved by first reducing it to an equivalent system of algebraic equations with real coefficients and isolating one of the roots using the generalized NewtonRaphson procedure, and then solving the factorized quadratic equation directly. The results are exemplified in Fig. 5 in which only the maximum moduli are plotted vs. K. t It has been shown [6] that the forward difference approximation difference scheme, and therefore, should not be used.
for icn yields a highly unstable
503
Numerical Stability in Wave Propagation Problems
0.2
04
06
06
I.0 Vclue
I.2
of
I.4
I.6
I.6
i 0
K
FIG. 5. Spherically symmetric wave propagation. Effect of r on stability. p=l,
8=1/S.
These numerical results indicate that for some values of r, the coefficient of artificial viscosity, the originally unstable difference scheme becomes conditionally stable. For small values of p, larger values of IY are required to obtain stability; and as p increases, the beneficial effect of r becomes more pronounced. The maximum value of the parameter K at which Ill I 1, is defined as the ‘critical value of K’, and is denoted as K,,.K,, is a measure of the stability range of that particular scheme, and for KS K,, the difference scheme is stable. In Figs. 7 through 10, KC, is plotted against r for various values of p and /I. It is seen that generally Kc, attains some maximum at a particular value of r. Furthermore, it is observed that as /3 increases, Kc, corresponding to a given value of I also tends to increase. In axisymmetric wave propagation, damping is introduced similarly through the equations of motion in terms of stresses as follows:
a2 +h
ar
0,
-a,+
dz+-r
a_ ..
-$-W
(21) in which qP and qr are pseudo-viscous stresses
ati
qr= pWrdr, and
(22) The discretized forms of equation (21), without the terms aq,/ar and aq,/dz, are given by equations (8) and (9). In order to include the artificial viscosity in the equations of motion, the quantities aq,/iJr and 8qZ/az are again discretized using backward difference formulas for
504
YUKSELUCKAN and ALFREDO H.-S. ANG
the first-order time derivatives, ti and 3, and central finite-difference approximations for the second-order space derivatives. The characteristic equations of the amplification matrix for this four-level difference scheme are,
azt3 -f-b,? -tc,4;+d2=0,
(23)
where ai=1-t-j&+jL41,
a,=l.+/3B2+j?A,,
b,= -2+(1-2j)B,f(l-3/W,,
h,= -2+(l-2j)B,i-(I-3P)A,,
c,=l+jB,--(l-3/+4,,
f’z= I +~&-(I
d1=-PA,,
ci2= -BA2.
-3@A,. (24)
Here, B, and B2 are the same quantities as those given by equations (13) and (14), and A, and A, are defined by A, E 2J(2)rK
. &,Ar sin , 2
AZ= 2J(2)lX~sm
. zkzAz .
(25)
2
The numerical solution of the complex cubits of equation (23) can be obtained through a procedure similar to that used in equation (18). A typical plot of It1 vs. K is shown in Fig. 6. Again, it can be concluded that for some values of I?, the difference scheme becomes stable; the range of conditional stability for a given value of I? increases with p and 8. It should be observed, however, that for an incompressible medium fv=OSO) the damping mechanism suggested herein does not produce a stable scheme.
05
06 Value
FIG. 6.
07 of
08
09
K
Axisymmetric wave propagation. Effect of r on stability. p--l, jI=O, v =~OGO
Numerical Stability in Wave Propagation
0.5 Value
od
0.2
0.3
0.4 Value
FIG.
9. Spherically symmetric wave propagation.
Relation between
05
06
I
0.7
Kc,
00
and r. p= 1.
I
0.9
of r
FIG.8. Spherically symmetric wave propagation.
Value
0.7
r
of
FIG. 7. Spherically symmetric wave propagation.
00 1
06
505
Problems
Relation between
Kc,
and r. p=5.
of r
Relation between Kc, and r. p=lO.
YUKSELUCKAN and ALFREDO H.-S. ANG
506
Value
of r
FIG. 10. Sphericallysymmetric wave propagation.
Figures 1l-14
Relation between Kcr and r. p=15.
show the variations of K,, vs. IY for various values of p, /I, and v. In
general, K,, reaches a maximum for some value of F, beyond which no improvement the stability range is noted;
it is this value of l? that is recommended
tions.
Value of r FIG. 11. Axisymmetric
wave propagation.
in
for practical applica-
Relation between Ker and F’. p= 1, v=O-@40.
507
Numeric4 Stability in Wave F’ropagatitm Problems
FIG. 12. Axisymmetric wave propagation. Relation between Ker and r. p=S, ~=0-040.
FIG. 13. Axisymmetric wave propagation. Relation between IJ& and r. p=lO, v=O-040.
Fro. 14. Msymmttric
wave Propagation. Relation between Ker and I’. pplS,
~=O-(h$o.
YIJKSEL UCKAN and ALFREDO H.-S. ANG
508
Transformations of dependent variables
Stability of numerical schemes may be assured also by changing the original problem through certain transformations of the dependent variables. The simple transformation,
UE--,
u r
(26)
changes the difference equation, equation (l), for wave motion in a spherically symmetric continuum into q+2_2q+‘++
(Ed!!> *r
2[8(Li;t’f-2U;+2+U1_+Z)+(l-2B)(U;=:
-2u,““+u;-+:)+j?(u:+,-2u~+u;_,)]-~ +(1-2/3)u;+‘+/W:],
P2 cdAt z 2CBU;+2 ( > (27)
in which the difference approximation corresponding to the first-order space derivative is absent. Consequently, the parameter B which appears in the characteristic equation of the amplification matrix for the difference equation of equation (27) becomes real-valued. For local stability, the following necessary condition is derived:
Unconditionally
stable, for /I z&.
G-m
The difference scheme defined by equation (27) is therefore, locally conditionally stable, and hence, is preferable to equation (1) as far as stability is concerned. Another form of transformation is,
u2 L . i3r0 r
(29)
It can be demonstrated that equation (29) reduces the continuous equation of motion for the spherically symmetric problem (30) into the simple wave equation (31) thus, the conditional stability is the same as that given in equation (7).
Numerical Stability in Wave Propagationproblems
509
For axisymmetric wave propagation, no simple transformation similar to equation (26) was found. However, it can be shown that the continuous equations of motion for axisymmetric wave propagation [6] are reduced to the following uncoupled equations
Lip2x+ l y+ a9 c,”
ar2
(32)
4rZ 2
if the following displacement potentials are used in the original differential equations of motions :
(33)
The discrete equivalents of equation (32) constitute
Unconditionally
a locally stable scheme provided
stable, for j32 ) ,
(34)
where ArSAz. Although the displacement potentials 4 and Y yield a conditionally stable difference scheme, the formulation in terms of these displacement potentials will, unavoidably, create additional problems with regard to the side conditions in practical problems. CONCLUSIONS
Discrete formulation of wave propagation problems can be accomplished conveniently on the basis of a lumped-parameter model and Newmark’s generalized integrator. The associated calculational stability conditions necessary to assure meaningful solutions are investigated. It is shown that the constant-coefficient difference schemes for one-dimensional and plane strain wave motions are conditionally stable for 8-z l/4, and unconditionally stable for /I2 l/4. However, the investigation of the local stability characteristics of the variable-coefficient schemes for wave motion in spherically symmetric and axisymmetric elastic media proved that instability is to be expected. In order to correct the numerical instability, a damping mechanism is suggested in the original difference schemes. It is observed that the linear artificial viscosity terms render the previously instable difference schemes conditionally stable at the expense of slight inaccuracy. For relatively small values of the coefficient of artificial viscosity, the deviation of a numerical solution from the true solution is negligible and suitable for practical purposes.
510
Yurcse~ UCKAN and ALFRED~H.-S. ANY
Conditionally stable difference schemes obtained with the same discretization procedure can be found also through certain transformations of the original dependent variables. In particular, the difference schemes formulated in terms of displacement potentials are conditionally stable, However, additional difficulties pertaining to side conditions in most practical problems may limit the applicability of this latter formulation. Acknowledgement-The study described herein was conducted as part of a research program on the numerical calculation of earthquake motions in the Department of Civil Engineering, University of Illinois at Urbana, Illinois. The program was supported by the National Science Foundation under research grant No. GK1858.
REFERENCES [l] A. H. -S. ANG, Numerical approach for wave motion in nonlinear solid media, Proc. Volume, Conjizrence on Matrix Methocis in Structural Mechanics, Wright-Patterson Air Force Base, Ohio, (1965). [2] J. VON NEUMANNand R. D. RICHTMYER,A method for the numerical calculation of hydrodynamic shocks, J. appl. Phy. 21,232 (1950). [3] N. M. NEWMARK, Method of computation for structural dynamics, J. Eng. Mech. Div. AXE, 85, No. EM3, Proc. Paper 2094 (1959). [4] R. D. RXXTMYERand K. W. MORTON, Difference Methods for Initial-Value Problems, 2nd Ed., Interscience, New York, (1967). [S] A. H. M. SAMEHand A, H. -S. ANG, Numerical analysis of axisymmetric wave propagation in elasticplastic layered media, Civil Engineering Studies, Structural Research Series, No. 335, University of Illinois, Urbana, Ill., May 1968. [6] Y. UCKAN and A. H. -S. ANG, A study of the numerical analysis of wave propagation in solid media. Civil Engineering Studies, Structural Research Series, No. 371, University of Illinois, Urbana, (1970). (Received 2 August 1971)
& stmtures,
Vol. 1, pp. 495-510. paRanloa Plwa 1971. Printed in orut
NUMERICS
Britain
STABILITY IN WAVE PROPAGATION PROBLEMS YUKSELUCKAN~
and ALFRERO H.-S. ANG~
University of Illinois, Urbana, Illinois 61801, U.S.A. Abstract-Stability characteristics for some finite-difference schemes obtained on the basis of a mathematically consistent lumped-parameter model and a suitable numerical integrator for wave propagation in spherically symmetric and axisymmetric ~nt~uo~ media are investigated. The numerical integrator used is Newmark’s /l-method. It is shown that the difference schemes for large radii which are identical to those for one-dimensional and plane strain wave motions are conditionally stable for values of B between 0 and l/6 and unconditionally stable for fi=1/4 and l/2. For points near the center or axis of symmetry, it is found that for all values of @the difference schemes are unstable. However, conditional stability may be obtained if one uses a dissipative ~h~isrn in the form of linear artificial viscosity; the effect of this dissipative mechanism on the stability characteristics of the considered difference schemes is also examined in detail. Conditionally stable difference schemes are also found by introducing appropriate dependent variable transformations.
Cd velocity of dilatational waves I imaginary constant J?’ t/At = index for time level r/Ar=index for distance of a generic point from the center of symmetry z/At=index for z-ordinate of a discrete point art&ial pseudo-viscous stresses in r- and z-directions, respectively distance of a point from the center of symmetry time independent variable displacement of a mass point in r-direction velocity of a mass point in r-direction acceleration of a mass point in r-dir&ion ; a factor in Newmark% integrator coefficient of linear artificial viscosity Ar, A& : tlnite space increments in r-, 8- and z-directions, respectively At time increment P mass density at, a,, ur direct stresses in r- and z-directions, respectively t eigenvalues of the amplification matrix n P 4 9% 4s r t u u ..
XNTRODUC3ION A WIDELY used method of solution for wave propagation problems in solid media is to discretize the continuous problem using a mathematical or physical model and to integrate the resulting system of differential-difference equations by the aid of a suitable numerical
t Present address: Middle East Technical University, Department of Computer Science, Ankara, Turkey; formerly research assistant in Civil Engineering, Univ. of Illinois, Urbana, Illinois. $ Professor of Civil Engineering. 495
496
YUKSELUCKANWK~
ALFREDOH.-S.ANG
integrator. An important matter associated with this approach is the question of convergence of the numerical solution to the continuous solution of the problem, which, in turn, is closely related to the stability properties of the finite-difference scheme through Lax’s Equivalence Theorem [4]. Unless the stability of a difference scheme ~orres~nd~ng to a properly posed initial value problem is ascertained, the resulting numerical solution may be far from representing the true solution. It is the principal objective of this study to analyze the stability characteristics of certain finite-difference schemes commonly used in wave propagation problems in higher dimensional solid media. Practical stability criteria are derived and presented in a form suitable for direct applications. The difference schemes are based on the assumptions of linearly elastic material behavior and small displacement; however, this does not unduly restrict the applicability of the derived stability criteria. For example, for a class of inelastic solids where the plastic velocity of propagation is smaller than the velocity of the elastic dilatational waves, the results remain applicable. Wave propagation in spherically symmetric and axisymmetric media are considered. As one-dimensional and plane strain wave propagation problems are, respectively, the limiting cases of the spherically symmetric and axisymmetric wave motions for points infinitely distant from the center of symmetry, they naturally fall within the scope of the present study. For all problems, the space variables are discretized using a discrete lumpedparameter model proposed by Ang [If; with this model, the resulting differential-di~erence equations invariably turn out to be the central finite-difference analogues of the corresponding continuous formulations, and hence are mathematically consistent. Time-wise discretization is accomplished with the aid of Newmark’s P-integrator [3]. The simple cases of the one-dimensional and plane strain wave propagation problems are of constant-coefficient type, for which the Fourier transform method [4] is suitable in deriving necessary conditions for stability. With spherically symmetric and axisymmetric wave motions, which involve variable-coefficient systems, the Fourier transform method yields only local stability conditions. The center of symmetry in a spherically symmetric medium and the points on the axis of symmetry in an axisymmetric medium are points of discontinuity at which the equations of motion are not valid. In order to render the displacements finite at such points of discontinuity, it is assumed that the displacements (and hence, the velocities and accelerations) of the points at the center or axis of symmetry vanish for all time. Furthermore, for the purpose of stability analysis, the boundary conditions, if any, are assumed to be linear and homogeneous. It is found that although the discrete formulations for one-dimensional and plane strain wave motions are marginally stable, the difference schemes for wave propagation in spherically symmetric and axisymmetric media are always unstable. Two different methods are suggested to remedy this situation: (1) The first approach, due to von Neumann and Richtmyer 121,uses a damping mechanism in the form of linear artificial viscosity in the original formulations. It is found that the modified difference schemes become conditionally stable, although the solution based on the damped system is invariably slightly different from that of the original problem. The amount of damping that produces a maximum stability range is investigated herein. (2) In the second approach, the original difference scheme is replaced with an equivalent formulation through appropriate transformation of the original variables, such that the resulting scheme becomes conditionally stable.
NumericalStability in Wave Propagation Probkms DI!XREXX
FORMULA’ITONS
AND STABILITY
497
ANALYSES
Spherically symmetric wave motion A displacement type approach for wave propagation in a sphe~~~y symmetric elastic continuum formulated on the basis of a mathematically consistent lard-p~ameter model and Newmark’s integrator yields the following difference equation [6] :
+ (I- 2BMZ
-u;,+:)+j?#+j
r&J;+2
+(1-2/3)U79+$J,
(1)
where u:zu(pAr,
nAt),
(2)
is the displacement in the r-direction, (r, 6, cp) being the orthogonal spherical coordinate system; /I is a numerical constant in Newmark’s integrator which can assume values of 0, l/12, l/8, l/6, l/4, or l/2; At is the uniform time increment, Ar is the uniform space mesh length, cd represents the velocity of propagation of dilatational waves in elastic continua, and &L--- i-l ‘Ar 2’
i=l,2,...,2P-t1
is)
is an index denoting the distance from the center of symmetry, and n = t/At indicates the time level at which computations are being performed. The domain of integration is given as
O
(4)
where B_2K2sm2kAr z 2+
K2 2K2. 2----sm-, P P
kAr 2
(5)
in which
and k is the frequency of the Fourier components of the error term, and 4 represents the eigenvalues of the amplification matrix.
YUKSELUCKANand ALEREDO H.-S. ANG
498
For infinitely large radii, asp+ co, the difference equation, equation (l), reduces to that corresponding to one-dimensional wave propagation, and B, as given by equation (5), becomes real-valued. The solution of the characteristic equation yields the following necessary condition for stability : 1 Ar -J(l-48)
c,bt<
if P<$,
Unconditionally
stable if p 2 $ ,
(7)
which ensures the satisfaction of the von Neumann necessary condition on the eigenvalues, Irl< 1. For points near the center of symmetry, the analysis of local stability is done by solving the quadratic equation (4) and computing the moduli 151for fixed values of p, In this case B is a complex quantity; consequently, the moduli It;] may be computed using the complex arithmetic feature of a digital computer. Figures 1 and 2 illustrate typical plots of the variations of I<]vs. Kf or various values ofp and /?, in which only the maximum moduli are shown. It is concluded that; (1) for all finite values of p, and all geometrically possible values of fi, the von Neumann condition, \
o-4
06
0.8
I.0
I.2
I.4
I.6
I.8
2.0
Value of K FIG. 1.
Spherically symmetric wave propagation. Effect ofp on
I 02
I 04
08
I
I
I
0.8
I.0
I.2
I I.4
stability. B=O.
I I.6
I I.8
Volue of K FIG. 2. Spherically symmetric wave propagation.
Effect of j3 on stability. p=l.
I 2.0
499
Numerical Stability in Wave Propagation Problems Axisymmetric
wave motion
On the same basis of formulation as above, wave pro~gation in ~sy~et~c media in cylindrical coordinates yields the foIlowing difference equations [6]:
elastic
nt2 ~~i-1,jtlffbr_+~,j-l~+~1~28~~U~++~,jf~~~~++~,j-~ n+ 1
-Ui-1,j+1+U;_+,‘ji-1 fUI-i,j-t)]P
)+B(Ut;+i,
j+l-@+I,
j-l-"l-l,
j-i-1
(9)
YUKSELUCKAN and ALFREW H.-S. ANG
500
in which, ~7,jZ u(pAr, qbz, nAl) , WY,j= w(pAr, qdz, nbt) ,
W)
are the displacement components in the I- and z-directions, respectively, with (r, 8, z) being the orthogonal cylindrical coordinate system; Ar and AZ are the uniform space mesh lengths in r- and z-directions, respectively, and cd and c, are the velocities of propagation of the dilatational and shear waves in an elastic medium, respectively. p and q are indices denoting the position of a generic point and defined as i-l
r
p’L\r=27 i-l
z
i=1,2 ).‘.)
2P+1,
j=l,2,...,2Q+l.
q=rz=25
(11)
The domain of integration is defined by O<~S PAr, O
and,
(12) where
&-$& 2
12c,2AtZ. k,Ar -sm, pAr2 2
(13
u2At2 . k,Ar 2Lsm2 c pAr2
- c~)Af2sin~ +$cdz pArAz
1 .
2
Numerical Stability in Wave Propagation Problems
501
For infinitely large radii; i.e. P-P co, the difference system, equations (8) and (9), reduce to those for plane strain wave motions, and B, and B, of equations (13) and (14) become real-valued. For this case, it can be proved [6] that the von Neumann condition is satisfied if 1 AZ -&2-8/I)
*<
for B
Unconditionally stable for /I 2 $,
(15)
where AZ< Ar. Equation (15), therefore, is a necessary condition for stability for plane strain wave propagation. For points at finite distances from the axis of symmetry, the difference scheme is of variable-coefficient type, and the analysis of local stability requires the determination of Iti from the characteristic equations, equation (12). The evaluation of the moduli is more complicated in this case; it can be shown [6] that ItI depends on the Poisson ratio, V, of the material as well. In Figs. 3 and 4, typical variations of the maximum moduli against K, [see equation (611for various values ofp, /?, and v are shown. The following conclusions can be drawn from these results: (1) For all values of p, B and v, the difference scheme for axisymmetric wave motion is locally unstable; (2) as p, or as /3 increases, the stability properties of the difference scheme improve; and (3) as v increases, the stability characteristics deteriorate, and the worst situation occurs for v = 0.50 (an incompressible medium).
I
050
I
I
1
I
I
02
04
0%
0.8
I.0
I I.2
I I.4
I I.8
I I.6
2.0
valye of K Fro. 3. Axisymmetric wave propagation. E&d of p on stability. p=O, v=O-00.
0’50’
I 0.2
I
o-4
I 06
I 08
I
I
I
B = 1/2f
IO 1.2 I.4 1.6 Value of K FIG. 4. Axisymmetric wave propagation. E%ct of fi on stability. p=l,
1.8
2.0
v=O-OO.
502
YUKBL UCKAN and ALFREDO H.-S. ANG
REMEDIES FOR NUMERICAL INSTABILITY Effect of linear artificial viscosity on stability
If a damping mechanism is appropriately introduced in an otherwise instable difference scheme, it may be possible to obtain a conditionally stable system. The linear artificial viscosity suggested by von Neumann and Richtmyer [2] is such a mechanism. In spherically symmetric wave motion, the damping term is introduced through the equation of motion in terms of stresses:
.. U!$ ++ - 06)+ h_ ar-PU,
(16)
where p is the mass density of the medium, a, and creare direct stresses, and qr is an artificial pseudo-viscous stress in the r-direction defined by qr = pl’c,ArE Zr’
(17)
in which I? is a constant referred to as the ‘coefficient of artificial viscosity’. The first two terms on the left hand-side of equation (16) belong to the equation of motion in terms of stresses from which the difference system, equation (I), is derived. aq,/& is discretized using the backward difference formula for the first-order time derivative? ti and the central difference formula for the second-order space derivative in conformity with the way the original scheme is discretized. The resulting difference scheme is a four level difference equation, whose characteristic equation is :
(18) where a, b, c and dare complex quantities given by a=1 +/3B+j?A, b= -2+(1-2/3)B+(l-3/3)A,
c= 1 +j?B-(1-3/Q/t, d= -PA,
(19)
A = 2 J(2)IYK sin”?,
(20)
in which
and B is the same quantity defined by equation (5). Equation (18) is a cubic with complex coefficients; this can be solved by first reducing it to an equivalent system of algebraic equations with real coefficients and isolating one of the roots using the generalized NewtonRaphson procedure, and then solving the factorized quadratic equation directly. The results are exemplified in Fig. 5 in which only the maximum moduli are plotted vs. K. t It has been shown [6] that the forward difference approximation difference scheme, and therefore, should not be used.
for icn yields a highly unstable
503
Numerical Stability in Wave Propagation Problems
0.2
04
06
06
I.0 Vclue
I.2
of
I.4
I.6
I.6
i 0
K
FIG. 5. Spherically symmetric wave propagation. Effect of r on stability. p=l,
8=1/S.
These numerical results indicate that for some values of r, the coefficient of artificial viscosity, the originally unstable difference scheme becomes conditionally stable. For small values of p, larger values of IY are required to obtain stability; and as p increases, the beneficial effect of r becomes more pronounced. The maximum value of the parameter K at which Ill I 1, is defined as the ‘critical value of K’, and is denoted as K,,.K,, is a measure of the stability range of that particular scheme, and for KS K,, the difference scheme is stable. In Figs. 7 through 10, KC, is plotted against r for various values of p and /I. It is seen that generally Kc, attains some maximum at a particular value of r. Furthermore, it is observed that as /3 increases, Kc, corresponding to a given value of I also tends to increase. In axisymmetric wave propagation, damping is introduced similarly through the equations of motion in terms of stresses as follows:
a2 +h
ar
0,
-a,+
dz+-r
a_ ..
-$-W
(21) in which qP and qr are pseudo-viscous stresses
ati
qr= pWrdr, and
(22) The discretized forms of equation (21), without the terms aq,/ar and aq,/dz, are given by equations (8) and (9). In order to include the artificial viscosity in the equations of motion, the quantities aq,/iJr and 8qZ/az are again discretized using backward difference formulas for
504
YUKSELUCKAN and ALFREDO H.-S. ANG
the first-order time derivatives, ti and 3, and central finite-difference approximations for the second-order space derivatives. The characteristic equations of the amplification matrix for this four-level difference scheme are,
azt3 -f-b,? -tc,4;+d2=0,
(23)
where ai=1-t-j&+jL41,
a,=l.+/3B2+j?A,,
b,= -2+(1-2j)B,f(l-3/W,,
h,= -2+(l-2j)B,i-(I-3P)A,,
c,=l+jB,--(l-3/+4,,
f’z= I +~&-(I
d1=-PA,,
ci2= -BA2.
-3@A,. (24)
Here, B, and B2 are the same quantities as those given by equations (13) and (14), and A, and A, are defined by A, E 2J(2)rK
. &,Ar sin , 2
AZ= 2J(2)lX~sm
. zkzAz .
(25)
2
The numerical solution of the complex cubits of equation (23) can be obtained through a procedure similar to that used in equation (18). A typical plot of It1 vs. K is shown in Fig. 6. Again, it can be concluded that for some values of I?, the difference scheme becomes stable; the range of conditional stability for a given value of I? increases with p and 8. It should be observed, however, that for an incompressible medium fv=OSO) the damping mechanism suggested herein does not produce a stable scheme.
05
06 Value
FIG. 6.
07 of
08
09
K
Axisymmetric wave propagation. Effect of r on stability. p--l, jI=O, v =~OGO
Numerical Stability in Wave Propagation
0.5 Value
od
0.2
0.3
0.4 Value
FIG.
9. Spherically symmetric wave propagation.
Relation between
05
06
I
0.7
Kc,
00
and r. p= 1.
I
0.9
of r
FIG.8. Spherically symmetric wave propagation.
Value
0.7
r
of
FIG. 7. Spherically symmetric wave propagation.
00 1
06
505
Problems
Relation between
Kc,
and r. p=5.
of r
Relation between Kc, and r. p=lO.
YUKSELUCKAN and ALFREDO H.-S. ANG
506
Value
of r
FIG. 10. Sphericallysymmetric wave propagation.
Figures 1l-14
Relation between Kcr and r. p=15.
show the variations of K,, vs. IY for various values of p, /I, and v. In
general, K,, reaches a maximum for some value of F, beyond which no improvement the stability range is noted;
it is this value of l? that is recommended
tions.
Value of r FIG. 11. Axisymmetric
wave propagation.
in
for practical applica-
Relation between Ker and F’. p= 1, v=O-@40.
507
Numeric4 Stability in Wave F’ropagatitm Problems
FIG. 12. Axisymmetric wave propagation. Relation between Ker and r. p=S, ~=0-040.
FIG. 13. Axisymmetric wave propagation. Relation between IJ& and r. p=lO, v=O-040.
Fro. 14. Msymmttric
wave Propagation. Relation between Ker and I’. pplS,
~=O-(h$o.
YIJKSEL UCKAN and ALFREDO H.-S. ANG
508
Transformations of dependent variables
Stability of numerical schemes may be assured also by changing the original problem through certain transformations of the dependent variables. The simple transformation,
UE--,
u r
(26)
changes the difference equation, equation (l), for wave motion in a spherically symmetric continuum into q+2_2q+‘++
(Ed!!> *r
2[8(Li;t’f-2U;+2+U1_+Z)+(l-2B)(U;=:
-2u,““+u;-+:)+j?(u:+,-2u~+u;_,)]-~ +(1-2/3)u;+‘+/W:],
P2 cdAt z 2CBU;+2 ( > (27)
in which the difference approximation corresponding to the first-order space derivative is absent. Consequently, the parameter B which appears in the characteristic equation of the amplification matrix for the difference equation of equation (27) becomes real-valued. For local stability, the following necessary condition is derived:
Unconditionally
stable, for /I z&.
G-m
The difference scheme defined by equation (27) is therefore, locally conditionally stable, and hence, is preferable to equation (1) as far as stability is concerned. Another form of transformation is,
u2 L . i3r0 r
(29)
It can be demonstrated that equation (29) reduces the continuous equation of motion for the spherically symmetric problem (30) into the simple wave equation (31) thus, the conditional stability is the same as that given in equation (7).
Numerical Stability in Wave Propagationproblems
509
For axisymmetric wave propagation, no simple transformation similar to equation (26) was found. However, it can be shown that the continuous equations of motion for axisymmetric wave propagation [6] are reduced to the following uncoupled equations
Lip2x+ l y+ a9 c,”
ar2
(32)
4rZ 2
if the following displacement potentials are used in the original differential equations of motions :
(33)
The discrete equivalents of equation (32) constitute
Unconditionally
a locally stable scheme provided
stable, for j32 ) ,
(34)
where ArSAz. Although the displacement potentials 4 and Y yield a conditionally stable difference scheme, the formulation in terms of these displacement potentials will, unavoidably, create additional problems with regard to the side conditions in practical problems. CONCLUSIONS
Discrete formulation of wave propagation problems can be accomplished conveniently on the basis of a lumped-parameter model and Newmark’s generalized integrator. The associated calculational stability conditions necessary to assure meaningful solutions are investigated. It is shown that the constant-coefficient difference schemes for one-dimensional and plane strain wave motions are conditionally stable for 8-z l/4, and unconditionally stable for /I2 l/4. However, the investigation of the local stability characteristics of the variable-coefficient schemes for wave motion in spherically symmetric and axisymmetric elastic media proved that instability is to be expected. In order to correct the numerical instability, a damping mechanism is suggested in the original difference schemes. It is observed that the linear artificial viscosity terms render the previously instable difference schemes conditionally stable at the expense of slight inaccuracy. For relatively small values of the coefficient of artificial viscosity, the deviation of a numerical solution from the true solution is negligible and suitable for practical purposes.
510
Yurcse~ UCKAN and ALFRED~H.-S. ANY
Conditionally stable difference schemes obtained with the same discretization procedure can be found also through certain transformations of the original dependent variables. In particular, the difference schemes formulated in terms of displacement potentials are conditionally stable, However, additional difficulties pertaining to side conditions in most practical problems may limit the applicability of this latter formulation. Acknowledgement-The study described herein was conducted as part of a research program on the numerical calculation of earthquake motions in the Department of Civil Engineering, University of Illinois at Urbana, Illinois. The program was supported by the National Science Foundation under research grant No. GK1858.
REFERENCES [l] A. H. -S. ANG, Numerical approach for wave motion in nonlinear solid media, Proc. Volume, Conjizrence on Matrix Methocis in Structural Mechanics, Wright-Patterson Air Force Base, Ohio, (1965). [2] J. VON NEUMANNand R. D. RICHTMYER,A method for the numerical calculation of hydrodynamic shocks, J. appl. Phy. 21,232 (1950). [3] N. M. NEWMARK, Method of computation for structural dynamics, J. Eng. Mech. Div. AXE, 85, No. EM3, Proc. Paper 2094 (1959). [4] R. D. RXXTMYERand K. W. MORTON, Difference Methods for Initial-Value Problems, 2nd Ed., Interscience, New York, (1967). [S] A. H. M. SAMEHand A, H. -S. ANG, Numerical analysis of axisymmetric wave propagation in elasticplastic layered media, Civil Engineering Studies, Structural Research Series, No. 335, University of Illinois, Urbana, Ill., May 1968. [6] Y. UCKAN and A. H. -S. ANG, A study of the numerical analysis of wave propagation in solid media. Civil Engineering Studies, Structural Research Series, No. 371, University of Illinois, Urbana, (1970). (Received 2 August 1971)