- Email: [email protected]
Decay of strong shocks in non-linear elasticity
Journal ofSound and Vibration (1985) 103(3), 371-378
DECAY OF STRONG M.
SHOCKS EL~ANOWSKI
IN NON-LINEAR AND
ELASTICITY
M. EPSTEIN
Department of Mechanical Engineering, University of Calgary, Calgary, Alberta, Canada (Received 3 August 1984, and in revised form 16 November 1984)
The method of singular surfaces is used to develop a numerical technique for calculating the decay and growth of the amplitude of shock waves in non-linear elastic materials. No restrictions are imposed on the magnitude of the initial conditions, which results in a coupling between the shock and the induced waves, governed by an infinite system of ordinary differential equations. Numerical examples illustrate the applicability of the technique in highly non-linear situations. 1. INTRODUCTION The literature on the theory of shock waves can be roughly divided into two major categories. The first consists of works devoted to the mathematical theory of hyperbolic
systems of partial differential equations and centers on the problem of existence and uniqueness of weak solutions. Typical of such approach are the works of Lax [ 1,2], Courant and Friedrichs [3] and others. In the second category, typified by the works of Truesdell and Toupin [4] and Chen [5], the point of view is adopted that only the events taking place at the wave front are of relevance. Information about these events can be gathered directly from the field equations, without resorting to a global solution, by means of the singular surface approach, pioneered by Hadamard [6]. Obviously, there is no contradiction between both points of view and often the second is included in the first. The theory of singular surfaces has been applied successfully to the propagation of discontinuities in derivatives of the displacement field of order 2 (acceleration waves) and higher. The reason for this success is that, regardless of the elastic constitutive law of the material, the decay of the amplitude of the wave is governed by an ordinary differential equation (along the bi-characteristics) which can be obtained from the original field equations. In the case of shock waves, however, this is not the case. In trying to obtain a decay equation for the shock, one obtains a certain coupling between its amplitude, the velocity of propagation and the amplitudes of secondary waves, which are all unknown. To avoid this situation, use is often made of the weak wave approximation [7], in which strains and shocks are assumed to be very small. Unfortunately, as is shown in this paper, the predictions of such an approximation are qualitatively wrong and, therefore, cannot be used as a starting point for a solution of cases in which the smallness assumption is not valid. The aim of this paper is to show how the infinite system of differential equations obtained by iterating the compatibility condition of the theory of singular surfaces can be used as a basis for a quantitative approximation to the correct decay of the shock amplitude. As is traditional in the literature [5], only the one-dimensional case is considered, to show the salient features of the method. The analysis is confined to non-linear elasticity, but it could obviously be extended to other kinds of materials. The numerical examples show that, as expected [ 11, only the initial second-order discontinuity plays a major role in determining the behaviour of the shock, at least for relatively short 371 0022460X/85/230371 +08 $03.00/O
1985 Academic Press Inc. (London) Limited
372
M. ELiANOWSKl
AND
M. EPSTEIN
time intervals. The examples have been solved numerically by a simple finite difference technique with the aid of a small computer code. As suggested above, and as illustrated by the examples, the iterated compatibility conditions can be successfully used for the analysis of the shock propagation phenomena. On the other hand it should be made absolutely clear that most of the qualitative information about propagation of shock waves, such as the breakdown of continuous solutions, propagation conditions for a shock (entropy-like conditions), as well as the dependence of the decay or growth of a solution on the character of the initial conditions and the nature of the constitutive law, can be gathered from the mathematical theory of hyperbolic systems of conservation laws [2,8].
2. THE GENERAL ITERATED COMPATIBILITY CONDITION Let R be a region in two-dimensional space-time, and let r be a curve with parameter s dividing R into two parts, denoted by R+ and R-. One says that r is a singular curve for a field IJ if I,!Jitself and/or its derivatives are discontinuous across r but smooth elsewhere in R. Under these conditions, according to Hadamard’s lemma [6], the chain rule of differentiation can be used on either side of r for evaluating derivatives of (1,on r itself. Thus, if $’ and I,!- denote appropriate limits of + at r from RC or R-,respectively, one can write d++/ds = $I; dx/ds + +; dtlds,
d$-/ds
= +, dx/ds+
$,t dt/ds, (2.1,2.2)
and hence one obtains
dL9llds = [$,x1Wds+[~,tl dtlds,
(2.3)
where brackets are used to indicate the jump of a quantity across r, i.e.,
[$I = ++- +-.
Figure
1. Sketch illustrating
(2.4)
the shock region.
In particular, with t chosen as a parameter, equation (2.3) can be rewritten as
(2.5) where [-I can be recognized as the Thomas or displacement V=dx/dt
derivative [5], and where (2.6)
STRONG
SHOCKS
IN
373
ELASTICITY
has the physical meaning of the instantaneous speed of a signal along K Equation (2.5) is called the compatibility condition for I/J.If now (I,is taken to be the nth spatial derivative of a function, IJ = u,,..., = P,
(2.7)
one obtains from equation (2.5) [~‘“‘I=,~,_
v[u(“+i)],
and similarly
[tic”‘] = [F]
- V[ tic”+‘)], (2.8,2.9)
where a superimposed dot denotes partial time differentiation. Combining equations (2.8 1 and (2.9) one obtains finally [;‘“‘I=,~]_ ~[U’“+“]_2vv[u~)]+ V2[U(n+2)]. (2.10) This is the most general form of the iterated compatibility condition that will be needed in this paper. In the above, advantage clearly has been taken of the one-dimensional character of the material space. The three-dimensional case has been treated in reference [9] along different lines. 3. THE GOVERNING
EQUATIONS
FOR A SHOCK
A propagating discontinuity of the velocity field in a continuum is called a shock, provided the displacement field is continuous. In the one-dimensional case the equation of motion in Lagrangian form in the absence of body forces is [l] T,, = pii,
(3.1)
where p is the density in the reference configuration, is related to the Cauchy stress component CTby
u is the displacement
T=(l+uJZa.
field and T (3.2)
In any case, for an elastic material in one dimension, u,, can be considered as a proper measure of strain E and the constitutive law can always be expressed uniquely as a function T= T(F).
(3.3)
On the shock itself the Lagrangian form of the equation of motion has the form [5] [T]+pV[ti]=O.
(3.4)
Let ai c [U(i)] S [ccl-l)]
(3.5)
denote the amplitude of the jump of the ith derivative of u. Then the propagation condition follows from equations (3.4) and (2.8): pV2a, = [ T].
(3.6)
Note that when T is a linear function of the strain e, i.e., T=AE,
(3.7)
equation (3.6) yields (A-pV*)a,
=O,
(3.8)
from which one obtains a definite value for the speed of propagation, V=J(AIP),
(3.9)
374
M. ELiANOWSKl
AND
M. EPSTEIN
which is independent of the amplitude of the shock. Similarly, for the case of the general constitutive law but limited to small strains and amplitudes (“weak shocks” [7]), one again obtains equation (3.8), this time through a Taylor expansion and neglecting terms of order a: or higher. However, in the general case of a strong shock propagating in a non-linear elastic material it is obvious that it is impossible to decouple wave speed from wave amplitude. In order to obtain further information in this situation one is led to consider the jumps of the field equation (3.1) and possibly of its derivatives. Taking the jump of equation (3.1) and using equation (2.10) with n =0 one obtains [T,,] = -p{ CC++ 2 va”,- V2a2},
(3.10)
where the notation (3.5) has been used. Again, in the linear case (3.7), and also for weak shocks, one has for the left-hand side of equation (3.10) [ Txl
= Aaz,
(3.11)
where the homogeneous case (A = const) has been considered for simplicity. Substituting this into equation (3.10) and using the propagation condition (3.8) one finds that the amplitude a, is completely determined by equation (3.10) and, in fact, in the homogeneous case it is a constant. In the general, genuinely non-linear case (“strong shock”), the term involving a, in equation (3.10) will not cancel out, so that equation (3.10) will not provide an explicit decay equation for the amplitude ur. 3ne is forced to consider higher order derivatives of equation (3.1), namely [ 7+fl)] = P[ &I’],
(3.12)
or, upon using equation (2.10), [ T(“+r)] = p{& - Gu,,, -2K?,+1+
V2a,+,},
(3.13)
which shows that the decay of a, is governed by an infinite system of ordinary differential equations involving the simultaneous decay of all of the a,‘~. The only exception to this statement is the Riemann problem: namely, the case in which initially all the a, (n 2 2) vanish. In this case, and only in this case, the equations (3.13) de-couple and in fact all the amplitudes, including a,, remain constant in time. On the other hand, suppose that a,, say, is initially non-zero. Then for the strong-shock case equation (3.10) cannot be satisfied with a, = constant and, therefore, one arrives at the conclusion-consistent with the mathematical theory of hyperbolic systems [3]-that the prediction of the weak shock approximation is qualitatively wrong! Depending on the nature of the constitutive law and of the initial conditions the shock amplitude may grow or decay (cf [5]), even when the shock is actually small. In the next section an algorithm to estimate such growth or decay numerically is proposed.
4. NUMERICAL
SOLUTION
As shown in the previous section, the decay or growth-and also the speed-of a shock are governed by an infinite system of equations: namely, the algebraic equation (3.6) and the sequence of equations (3.13) (note that, with some abuse of notation, equation (3.10) is obtained from equation (3.13) with n = 0). For the sake of definiteness a quadratic constitutive law can be considered for the numerical example: i.e., T(E) = AE + BE’,
(4.1)
STRONG
SHOCKS
IN
375
ELASTICITY
where A and B are material constants (homogeneous case). It is assumed, moreover, that the solid ahead of the shock is in its natural, stress-free, state. Under these conditions it follows that, according to definition (2.4), [T] = Aa, - Ba;,
[T’“‘] = Aa,,,
and
- B 1
aian-i+2.
(4.2,4.3)
i=l
With the use of equation (4.2), the propagation A-pV2=
condition (for a, f 0) reduces to Ba,
(4.4)
and equation (3.13) becomes (4.5)
where equation (4.4) has been utilized to eliminate three terms. Note that for n = 0 equation (4.5) can be rewritten as (2)
=
(4.6)
Wp)4a2,
which is the form of equation (3.10) in the particular introducing the change of variables
case under consideration.
b,=a,/(i-l)!
By (4.7)
equation (4.5) can be simplified to the form ~,+1+2@L+,
-~&,=B(n+l)“~‘bib,,., P
fornal,
,=,
and
(Vb:) = (B/p)b:b,,
n =O.
for (4.8,4.9)
The structure of the above system (upon assuming that V has been eliminated with the help of equation (4.4)) is as follows: 6, = K(b,, M 62=
F2$,,
6,
b,,
b2,
b,)
2”=EUk,h-1,b,, bz,. . . , bn+,).
(4.10)
One can observe, first of all, that despite the presence of second$erivatives this is a first order-system. Indeed, the first equation allows one to express b, in terms of b,, b2, b, and b,, and similarly all second derivatives are implicitly expressed in terms of first derivatives. In discretizing this system by a simple-minded forward finite difference scheme, one obtains an algebraic system of equations of the form
b,,i+,= b1.i+ hG,(b,,i,b2,i) b,,i+,= b2.i+ hGz(bl,z,b.r, h,i, bl,i-t, b,,i+,) bn,i+l= b,i + hGn(b,,i,h,i, . . . bn+l,i,bn-I,(-1,bn-l,i+l)v 3
(4.11)
where h is the time increment and where the second subscript denotes the respective time step. Thus formulated, it is obvious that for a fixed h and a given i, b,,i+, can be obtained in principle (exactly, within the discretization adopted) as a function b,,i+, = H(br.0, b2.0,. . . 9bi+2,0)*
(4.12)
376
M. ELiANOWSKI
AND
M. EPSTEIN
In other words, for a fixed time span, the number of initial conditions, and of equations from the system, involved in finding 6, is equal to the number of time increments plus 1. In the examples this method for finding the decay or growth of b, on the basis of the solution of the complete system of difference equations will be referred to, as method No. 1. The question now arises as to whether all of the initial values listed as arguments of H are of relevance, or whether the higher ones can be systematically neglected without affecting the solution appreciably. One way to effect this neglecting of higher order derivatives consists of truncating the system by arbitrarily assuming b, = 0 for some M in equations (4.10). This will automatically truncate the difference system as well, to yield (4.13)
b,,i+,= K(b,,,, b2.0,. . . bra,,) 7
for all i. Obviously, K will coincide with H for the first M - 1 steps, but will differ for all successive steps, even if the higher derivative jumps are originally zero. It is observed, however, that the tendency of these higher jumps is to grow (or decay) at a higher rate than the growth (or decay) of the lower jumps. This implies that whenever the tendency of b, is to decay, the truncation procedure will yield extremely accurate results for reasonably long periods of time. On the other hand, when the tendency to grow occurs, the influence of the missing terms will be felt sooner. In the examples, this truncation procedure is referred to as Method 2. More sophisticated truncation techniques are conceivable, but will not be discussed here. For either method, a simple computer code has been written in BASIC language. The storage needed is so small that a home VIC-20 computer can handle cases up to n = 40. For higher n and also to attenuate the possible influence of round-off error, however, the program was also translated into FORTRAN and run in double precision arithmetic on a CDC main frame. For the numerical examples p = 1, A = 1 and B = 0.5, and different values are introduced for the initial conditions, as shown on Table 1. These values have been chosen so as to test the method under extreme conditions and do not necessarily correspond to any TABLE
1
Initial conditions Example 1 2 3
aI -1 -0.001
a2
a,N>2
-1 +1 -1
0 0 0
physical situation, apart from the fact that initially Ba,
Convergence check at t = 4 for Method 1, Example 2 At M
V aI a2
0.4 10 1.802724 4.499627 6.237010
0.26 20 1.841950 4.785560 5.064460
0.16 40 1.861603 4.932621 4.863638
0.05 80 1.871554 5m5430 4.791594
STRONG
4.5
-
4.0
-
3.5
SHOCKS
IN
371
ELASTICITY
-
Met hod 2
#=20.
0
~=4, 2.5
At=o.\ A,=o,I
-
Example
Figure
2. Relative
shock (a) amplitude
I
and (b) speed versus time.
378
M.
ELiANOWSKI AND M. EPSTEIN
how the results for Example 2 are affected at t = 4 by the change of the time increment. These results show in fact the converging character of Method 1. On the other hand, convergence of the finite difference scheme of Method 2 was tested by adopting At = 0.1, 0.01 and O+OOl, while to check convergence with respect to truncation, the first example was run with M = 6 and M = 11. Tables 3 and 4 show how the results for Example 1 at t = 4 are affected by these changes. TABLE
3
Truncation-wise convergence check at t =4 Method 2, Example 1 (At = 0.1) M V
4 1.14293
6 1.14297
1 1 1.14297
a1
-0.61258 -0.38412 -0.028059
-0.61278 -0.38337 -0.030335
-0.61278 -0.38337 -0.030355
a2 a3
for
TABLE 4 Diference
scheme
convergence
At
0.1
V a1
1.14293 -0.61258
at t = 4 for Method
1 (M = 4)
2, Example
0.01 1.14404 - 0.61766
0.001 1.14@5 -0.61770
ACKNOWLEDGMENT This work has been supported in part through and Engineering Research Council of Canada.
an operating
grant of the Natural
Sciences
REFERENCES 1. P. D. LAX 1957 Communications on Pure and Applied Mathematics 10, 537-566. Hyperbolic systems of conservation laws, II. 2. P. D. LAX 1973 Conference Board of the Mathematical Sciences 11, SIAM. Hyperbolic systems of conservation laws and the mathematical theory of shock waves. 3. R. COURANT and K. 0. FRIEDRICHS 1948 Supersonic Flow and Shock Wuues. New York: Interscience. 4. C. TRUESDELL and R. TOUPIN 1960 Encyclopedia of Physics III (Editor S. FIugge). Berlin: Springer. The classical field theories. 5. P. J. CHEN 1973 in Encyclopedia ofphysics Via/3, (Editor S. Flugge). Berlin: Springer. Growth and decay of waves in solids. 6. J. HADAMARD 1903 Lecons sur La Propagation Ondes et les Equations de I’Hydrodynamique. Paris: Hermann. 7. H. COHEN and A. B. WHITMAN 1977 Journal ofSound and Vibration 51, 283-302. Waves in elastic rods. 8. J. SMOLLER 1983 Shock Waoes and Reaction-Di$usion Equations. Berlin: Springer. 9. H. COHEN and M. EPSTEIN 1983 Mechanics Research Communication, 10, 37-44. A note on nonlinear elastic wave propagation. (Also 1982 Department of Mechanical Engineering, University of Calgary, Report No. 220.)
DECAY OF STRONG M.
SHOCKS EL~ANOWSKI
IN NON-LINEAR AND
ELASTICITY
M. EPSTEIN
Department of Mechanical Engineering, University of Calgary, Calgary, Alberta, Canada (Received 3 August 1984, and in revised form 16 November 1984)
The method of singular surfaces is used to develop a numerical technique for calculating the decay and growth of the amplitude of shock waves in non-linear elastic materials. No restrictions are imposed on the magnitude of the initial conditions, which results in a coupling between the shock and the induced waves, governed by an infinite system of ordinary differential equations. Numerical examples illustrate the applicability of the technique in highly non-linear situations. 1. INTRODUCTION The literature on the theory of shock waves can be roughly divided into two major categories. The first consists of works devoted to the mathematical theory of hyperbolic
systems of partial differential equations and centers on the problem of existence and uniqueness of weak solutions. Typical of such approach are the works of Lax [ 1,2], Courant and Friedrichs [3] and others. In the second category, typified by the works of Truesdell and Toupin [4] and Chen [5], the point of view is adopted that only the events taking place at the wave front are of relevance. Information about these events can be gathered directly from the field equations, without resorting to a global solution, by means of the singular surface approach, pioneered by Hadamard [6]. Obviously, there is no contradiction between both points of view and often the second is included in the first. The theory of singular surfaces has been applied successfully to the propagation of discontinuities in derivatives of the displacement field of order 2 (acceleration waves) and higher. The reason for this success is that, regardless of the elastic constitutive law of the material, the decay of the amplitude of the wave is governed by an ordinary differential equation (along the bi-characteristics) which can be obtained from the original field equations. In the case of shock waves, however, this is not the case. In trying to obtain a decay equation for the shock, one obtains a certain coupling between its amplitude, the velocity of propagation and the amplitudes of secondary waves, which are all unknown. To avoid this situation, use is often made of the weak wave approximation [7], in which strains and shocks are assumed to be very small. Unfortunately, as is shown in this paper, the predictions of such an approximation are qualitatively wrong and, therefore, cannot be used as a starting point for a solution of cases in which the smallness assumption is not valid. The aim of this paper is to show how the infinite system of differential equations obtained by iterating the compatibility condition of the theory of singular surfaces can be used as a basis for a quantitative approximation to the correct decay of the shock amplitude. As is traditional in the literature [5], only the one-dimensional case is considered, to show the salient features of the method. The analysis is confined to non-linear elasticity, but it could obviously be extended to other kinds of materials. The numerical examples show that, as expected [ 11, only the initial second-order discontinuity plays a major role in determining the behaviour of the shock, at least for relatively short 371 0022460X/85/230371 +08 $03.00/O
1985 Academic Press Inc. (London) Limited
372
M. ELiANOWSKl
AND
M. EPSTEIN
time intervals. The examples have been solved numerically by a simple finite difference technique with the aid of a small computer code. As suggested above, and as illustrated by the examples, the iterated compatibility conditions can be successfully used for the analysis of the shock propagation phenomena. On the other hand it should be made absolutely clear that most of the qualitative information about propagation of shock waves, such as the breakdown of continuous solutions, propagation conditions for a shock (entropy-like conditions), as well as the dependence of the decay or growth of a solution on the character of the initial conditions and the nature of the constitutive law, can be gathered from the mathematical theory of hyperbolic systems of conservation laws [2,8].
2. THE GENERAL ITERATED COMPATIBILITY CONDITION Let R be a region in two-dimensional space-time, and let r be a curve with parameter s dividing R into two parts, denoted by R+ and R-. One says that r is a singular curve for a field IJ if I,!Jitself and/or its derivatives are discontinuous across r but smooth elsewhere in R. Under these conditions, according to Hadamard’s lemma [6], the chain rule of differentiation can be used on either side of r for evaluating derivatives of (1,on r itself. Thus, if $’ and I,!- denote appropriate limits of + at r from RC or R-,respectively, one can write d++/ds = $I; dx/ds + +; dtlds,
d$-/ds
= +, dx/ds+
$,t dt/ds, (2.1,2.2)
and hence one obtains
dL9llds = [$,x1Wds+[~,tl dtlds,
(2.3)
where brackets are used to indicate the jump of a quantity across r, i.e.,
[$I = ++- +-.
Figure
1. Sketch illustrating
(2.4)
the shock region.
In particular, with t chosen as a parameter, equation (2.3) can be rewritten as
(2.5) where [-I can be recognized as the Thomas or displacement V=dx/dt
derivative [5], and where (2.6)
STRONG
SHOCKS
IN
373
ELASTICITY
has the physical meaning of the instantaneous speed of a signal along K Equation (2.5) is called the compatibility condition for I/J.If now (I,is taken to be the nth spatial derivative of a function, IJ = u,,..., = P,
(2.7)
one obtains from equation (2.5) [~‘“‘I=,~,_
v[u(“+i)],
and similarly
[tic”‘] = [F]
- V[ tic”+‘)], (2.8,2.9)
where a superimposed dot denotes partial time differentiation. Combining equations (2.8 1 and (2.9) one obtains finally [;‘“‘I=,~]_ ~[U’“+“]_2vv[u~)]+ V2[U(n+2)]. (2.10) This is the most general form of the iterated compatibility condition that will be needed in this paper. In the above, advantage clearly has been taken of the one-dimensional character of the material space. The three-dimensional case has been treated in reference [9] along different lines. 3. THE GOVERNING
EQUATIONS
FOR A SHOCK
A propagating discontinuity of the velocity field in a continuum is called a shock, provided the displacement field is continuous. In the one-dimensional case the equation of motion in Lagrangian form in the absence of body forces is [l] T,, = pii,
(3.1)
where p is the density in the reference configuration, is related to the Cauchy stress component CTby
u is the displacement
T=(l+uJZa.
field and T (3.2)
In any case, for an elastic material in one dimension, u,, can be considered as a proper measure of strain E and the constitutive law can always be expressed uniquely as a function T= T(F).
(3.3)
On the shock itself the Lagrangian form of the equation of motion has the form [5] [T]+pV[ti]=O.
(3.4)
Let ai c [U(i)] S [ccl-l)]
(3.5)
denote the amplitude of the jump of the ith derivative of u. Then the propagation condition follows from equations (3.4) and (2.8): pV2a, = [ T].
(3.6)
Note that when T is a linear function of the strain e, i.e., T=AE,
(3.7)
equation (3.6) yields (A-pV*)a,
=O,
(3.8)
from which one obtains a definite value for the speed of propagation, V=J(AIP),
(3.9)
374
M. ELiANOWSKl
AND
M. EPSTEIN
which is independent of the amplitude of the shock. Similarly, for the case of the general constitutive law but limited to small strains and amplitudes (“weak shocks” [7]), one again obtains equation (3.8), this time through a Taylor expansion and neglecting terms of order a: or higher. However, in the general case of a strong shock propagating in a non-linear elastic material it is obvious that it is impossible to decouple wave speed from wave amplitude. In order to obtain further information in this situation one is led to consider the jumps of the field equation (3.1) and possibly of its derivatives. Taking the jump of equation (3.1) and using equation (2.10) with n =0 one obtains [T,,] = -p{ CC++ 2 va”,- V2a2},
(3.10)
where the notation (3.5) has been used. Again, in the linear case (3.7), and also for weak shocks, one has for the left-hand side of equation (3.10) [ Txl
= Aaz,
(3.11)
where the homogeneous case (A = const) has been considered for simplicity. Substituting this into equation (3.10) and using the propagation condition (3.8) one finds that the amplitude a, is completely determined by equation (3.10) and, in fact, in the homogeneous case it is a constant. In the general, genuinely non-linear case (“strong shock”), the term involving a, in equation (3.10) will not cancel out, so that equation (3.10) will not provide an explicit decay equation for the amplitude ur. 3ne is forced to consider higher order derivatives of equation (3.1), namely [ 7+fl)] = P[ &I’],
(3.12)
or, upon using equation (2.10), [ T(“+r)] = p{& - Gu,,, -2K?,+1+
V2a,+,},
(3.13)
which shows that the decay of a, is governed by an infinite system of ordinary differential equations involving the simultaneous decay of all of the a,‘~. The only exception to this statement is the Riemann problem: namely, the case in which initially all the a, (n 2 2) vanish. In this case, and only in this case, the equations (3.13) de-couple and in fact all the amplitudes, including a,, remain constant in time. On the other hand, suppose that a,, say, is initially non-zero. Then for the strong-shock case equation (3.10) cannot be satisfied with a, = constant and, therefore, one arrives at the conclusion-consistent with the mathematical theory of hyperbolic systems [3]-that the prediction of the weak shock approximation is qualitatively wrong! Depending on the nature of the constitutive law and of the initial conditions the shock amplitude may grow or decay (cf [5]), even when the shock is actually small. In the next section an algorithm to estimate such growth or decay numerically is proposed.
4. NUMERICAL
SOLUTION
As shown in the previous section, the decay or growth-and also the speed-of a shock are governed by an infinite system of equations: namely, the algebraic equation (3.6) and the sequence of equations (3.13) (note that, with some abuse of notation, equation (3.10) is obtained from equation (3.13) with n = 0). For the sake of definiteness a quadratic constitutive law can be considered for the numerical example: i.e., T(E) = AE + BE’,
(4.1)
STRONG
SHOCKS
IN
375
ELASTICITY
where A and B are material constants (homogeneous case). It is assumed, moreover, that the solid ahead of the shock is in its natural, stress-free, state. Under these conditions it follows that, according to definition (2.4), [T] = Aa, - Ba;,
[T’“‘] = Aa,,,
and
- B 1
aian-i+2.
(4.2,4.3)
i=l
With the use of equation (4.2), the propagation A-pV2=
condition (for a, f 0) reduces to Ba,
(4.4)
and equation (3.13) becomes (4.5)
where equation (4.4) has been utilized to eliminate three terms. Note that for n = 0 equation (4.5) can be rewritten as (2)
=
(4.6)
Wp)4a2,
which is the form of equation (3.10) in the particular introducing the change of variables
case under consideration.
b,=a,/(i-l)!
By (4.7)
equation (4.5) can be simplified to the form ~,+1+2@L+,
-~&,=B(n+l)“~‘bib,,., P
fornal,
,=,
and
(Vb:) = (B/p)b:b,,
n =O.
for (4.8,4.9)
The structure of the above system (upon assuming that V has been eliminated with the help of equation (4.4)) is as follows: 6, = K(b,, M 62=
F2$,,
6,
b,,
b2,
b,)
2”=EUk,h-1,b,, bz,. . . , bn+,).
(4.10)
One can observe, first of all, that despite the presence of second$erivatives this is a first order-system. Indeed, the first equation allows one to express b, in terms of b,, b2, b, and b,, and similarly all second derivatives are implicitly expressed in terms of first derivatives. In discretizing this system by a simple-minded forward finite difference scheme, one obtains an algebraic system of equations of the form
b,,i+,= b1.i+ hG,(b,,i,b2,i) b,,i+,= b2.i+ hGz(bl,z,b.r, h,i, bl,i-t, b,,i+,) bn,i+l= b,i + hGn(b,,i,h,i, . . . bn+l,i,bn-I,(-1,bn-l,i+l)v 3
(4.11)
where h is the time increment and where the second subscript denotes the respective time step. Thus formulated, it is obvious that for a fixed h and a given i, b,,i+, can be obtained in principle (exactly, within the discretization adopted) as a function b,,i+, = H(br.0, b2.0,. . . 9bi+2,0)*
(4.12)
376
M. ELiANOWSKI
AND
M. EPSTEIN
In other words, for a fixed time span, the number of initial conditions, and of equations from the system, involved in finding 6, is equal to the number of time increments plus 1. In the examples this method for finding the decay or growth of b, on the basis of the solution of the complete system of difference equations will be referred to, as method No. 1. The question now arises as to whether all of the initial values listed as arguments of H are of relevance, or whether the higher ones can be systematically neglected without affecting the solution appreciably. One way to effect this neglecting of higher order derivatives consists of truncating the system by arbitrarily assuming b, = 0 for some M in equations (4.10). This will automatically truncate the difference system as well, to yield (4.13)
b,,i+,= K(b,,,, b2.0,. . . bra,,) 7
for all i. Obviously, K will coincide with H for the first M - 1 steps, but will differ for all successive steps, even if the higher derivative jumps are originally zero. It is observed, however, that the tendency of these higher jumps is to grow (or decay) at a higher rate than the growth (or decay) of the lower jumps. This implies that whenever the tendency of b, is to decay, the truncation procedure will yield extremely accurate results for reasonably long periods of time. On the other hand, when the tendency to grow occurs, the influence of the missing terms will be felt sooner. In the examples, this truncation procedure is referred to as Method 2. More sophisticated truncation techniques are conceivable, but will not be discussed here. For either method, a simple computer code has been written in BASIC language. The storage needed is so small that a home VIC-20 computer can handle cases up to n = 40. For higher n and also to attenuate the possible influence of round-off error, however, the program was also translated into FORTRAN and run in double precision arithmetic on a CDC main frame. For the numerical examples p = 1, A = 1 and B = 0.5, and different values are introduced for the initial conditions, as shown on Table 1. These values have been chosen so as to test the method under extreme conditions and do not necessarily correspond to any TABLE
1
Initial conditions Example 1 2 3
aI -1 -0.001
a2
a,N>2
-1 +1 -1
0 0 0
physical situation, apart from the fact that initially Ba,
Convergence check at t = 4 for Method 1, Example 2 At M
V aI a2
0.4 10 1.802724 4.499627 6.237010
0.26 20 1.841950 4.785560 5.064460
0.16 40 1.861603 4.932621 4.863638
0.05 80 1.871554 5m5430 4.791594
STRONG
4.5
-
4.0
-
3.5
SHOCKS
IN
371
ELASTICITY
-
Met hod 2
#=20.
0
~=4, 2.5
At=o.\ A,=o,I
-
Example
Figure
2. Relative
shock (a) amplitude
I
and (b) speed versus time.
378
M.
ELiANOWSKI AND M. EPSTEIN
how the results for Example 2 are affected at t = 4 by the change of the time increment. These results show in fact the converging character of Method 1. On the other hand, convergence of the finite difference scheme of Method 2 was tested by adopting At = 0.1, 0.01 and O+OOl, while to check convergence with respect to truncation, the first example was run with M = 6 and M = 11. Tables 3 and 4 show how the results for Example 1 at t = 4 are affected by these changes. TABLE
3
Truncation-wise convergence check at t =4 Method 2, Example 1 (At = 0.1) M V
4 1.14293
6 1.14297
1 1 1.14297
a1
-0.61258 -0.38412 -0.028059
-0.61278 -0.38337 -0.030335
-0.61278 -0.38337 -0.030355
a2 a3
for
TABLE 4 Diference
scheme
convergence
At
0.1
V a1
1.14293 -0.61258
at t = 4 for Method
1 (M = 4)
2, Example
0.01 1.14404 - 0.61766
0.001 1.14@5 -0.61770
ACKNOWLEDGMENT This work has been supported in part through and Engineering Research Council of Canada.
an operating
grant of the Natural
Sciences
REFERENCES 1. P. D. LAX 1957 Communications on Pure and Applied Mathematics 10, 537-566. Hyperbolic systems of conservation laws, II. 2. P. D. LAX 1973 Conference Board of the Mathematical Sciences 11, SIAM. Hyperbolic systems of conservation laws and the mathematical theory of shock waves. 3. R. COURANT and K. 0. FRIEDRICHS 1948 Supersonic Flow and Shock Wuues. New York: Interscience. 4. C. TRUESDELL and R. TOUPIN 1960 Encyclopedia of Physics III (Editor S. FIugge). Berlin: Springer. The classical field theories. 5. P. J. CHEN 1973 in Encyclopedia ofphysics Via/3, (Editor S. Flugge). Berlin: Springer. Growth and decay of waves in solids. 6. J. HADAMARD 1903 Lecons sur La Propagation Ondes et les Equations de I’Hydrodynamique. Paris: Hermann. 7. H. COHEN and A. B. WHITMAN 1977 Journal ofSound and Vibration 51, 283-302. Waves in elastic rods. 8. J. SMOLLER 1983 Shock Waoes and Reaction-Di$usion Equations. Berlin: Springer. 9. H. COHEN and M. EPSTEIN 1983 Mechanics Research Communication, 10, 37-44. A note on nonlinear elastic wave propagation. (Also 1982 Department of Mechanical Engineering, University of Calgary, Report No. 220.)