- Email: [email protected]
A Riemann solver and a second-order Godunov method for elastic-plastic wave propagation in solids
Int. J. Impact Enong Vol. 13, No. 3, pp. 463-478, 1993 Printed in Great Britain
0734-743X/93 $6.00 + 0.00 © 1993 Pergamon Press Lid
A RIEMANN SOLVER AND A SECOND-ORDER GODUNOV METHOD FOR ELASTIC-PLASTIC WAVE PROPAGATION IN SOLIDS X. LIN a n d J. BALLMANN Lehr-und Forschungsgebiet fiir Mechanik der RWTH Aachen, Templergraben 64, 5100 Aachen, Germany
{Received 5 February 1991; and in revised form 27 March 1992)
Summary By use of three basic paths of elastic-plastic loading, a Riemann solver is established for the one-dimensional combined longitudinal and torsional stress wave problem of the thin-walled elastic-plastic tube. Based on this, a second-order Godunov method is developed for the numerical computation of elastic-plastic waves. Finally the numerical method is applied to some basic problems with known exact solutions and the comparison of results is made.
1. I N T R O D U C T I O N
Many solid materials behave plastically when they are subjected to high impulse loads. This may result in locally plastic wave propagation and propagating discontinuity surfaces like shocks, i.e. jumps of stresses and strains, or material interfaces. The wave structure, which has a great influence on the dynamic yield pattern and fracture processes, can only be correctly preserved by methods taking into account the probability occurring physical discontinuities. Unfortunately, most practical problems of mechanical waves in solids cannot be solved analytically due to the shape of the boundaries, stress hardening effects and complicated loading functions. During the last decades different numerical approaches have been introduced by numerous authors, involving finite element analysis [ 1], boundary element methods [2] and methods of characteristics [3,4]. The advantage of FEM is the easy implementation of nonlinearities of the PDEs and complicated boundaries, but it remains difficult to get the correct behaviour of the wave fronts. BEM is mainly restricted to linear problems. The method of characteristics truly integrates mathematically exact equations which have a direct physical interpretation in wave propagation, but this method has to consider discontinuities as interior boundaries within the solution domain, which are not easily managed for multi-dimensional problems. Therefore, it is necessary to have numerical methods for multi-dimensional wave problems which capture discontinuities and preserve simultaneously the principles governing the wave propagation phenomenon. In order to overcome the described difficulties, we take the idea of Godunov's method [5] in gasdynamics as an example, where a Riemann problem is the starting point of the method of solution. By this method, Nagayama [6] considered the behaviour of one-dimensional solids under high-pressure conditions, Trangenstein and Pember [7] dealt with the longitudinal wave motion in an elastic-plastic bar. In our case, the basic problem is the propagation of combined longitudinal and torsional waves in a thin-walled tube, where simple wave solutions can be used for the analysis of the wave pattern. The problem has been studied, e.g. by Clifton [8], Ting [9-11] and Goel and Malvern [12]. Their fundamental contributions formed the basic understanding for the treatment of waves involving elastic-plastic effects by simple wave analysis. In this paper, we will first present the construction of a Riemann solver for combined longitudinal and torsional stress waves in a thin-walled tube. In order to render the possibility of any arbitrary loading paths in stress space, three elementary paths will be introduced for the elastic-plastic loading history. Then we will deal with a second-order 463
464
X. LIN and J. BALLMANN
Godunov method. Finally, some numerical results will be given and compared with exact solutions. 2. BASIC GOVERNING EQUATIONS Consider a thin-walled tube being composed of an elastic-plastic stress-hardening material, which is subjected to a combined longitudinal and torsional shock loading, see Fig. 1. Following [ 8 - 1 2 ] for the derivation of the system of PDEs to be solved, we start with the local balance of momentum with components in longitudinal and circumferential direction, and the linearized strain-displacement relations partially derived with respect to time: Ow
Ot
Of
-
(1)
Ox'
where
w =
,
f=
.
u, tr, e shall represent the longitudinal particle velocity, stress and strain, and v, r, 7 the circumferential ones; p is the mass density, t is the time, and x is the distance measured along the axis of the tube. If the material of the tube is an isotropic work-hardening one, the incremental formulations of the material law have the following form
de=
+ 02 /I d~ + H~z dr'
dT= Hcrzda + (G+ HO2r2)dz,
(2)
where E is the Young's modulus and G the shear modulus. H is a nonlinear function of the yield stress x. The yield condition is given by
with 0 = v/3 for the yon Mises yield condition and 0 = 2 for the Tresca yield condition. The function H is defined by
where, g(tr) is the slope of the stress-strain curve in a simple tension test, i.e. v(~)
=
dtr --. dE
(5)
In the elastic region, g(Ox) = E, so H = 0, in the plastic region for the linearly isotropic work-hardening material, g(Ox) = EF (const).
FIG. 1. Thin-walled tube subjected to a sudden combined loading Or, r).
Riemann solver and Godunov method for solids
465
Substituting Eqns (2) into Eqn (1), the governing equations can be rewritten as AOf _ Of dt dx'
(6)
where
A =
0
0
p
0
0
0
p
-E + 0---T-
Hez
0
0
Hez
--+HO2z 2 G
0
0
1
He 2
1
For Eqn (6) a simple wave solution can be found, which connects two domains of different homogeneous states. Therein f is a function of x / t = c only. Equation (6) is then rewritten in the form (cA + I) df = 0,
(7)
where c is a free parameter. The necessary condition for a non-trivial solution df :/= 0 to exist, is llcA + Ill -- 0.
(8)
Solving this equation for c yields two solutions. In the elastic region (when H = 0), these are the characteristic wave speeds: Co=x/~,
c 2 = ~ p G-,
(9)
and in the plastic region, the two solutions are the wave speeds c~ and cf, which are correlated with the stresses (a, T) and satisfy 0 ~< cs ~< c2 ~< cf ~ Co.
(10)
The speeds c~ and cf are distinguished as slow and fast. The solutions of Eqn (8) contain the sign + for right-running waves and the sign - for left-running waves. In the right-ward simple wave region: x / t = c ( x > 0), where c may assume the values of Co, cf, c2 or c~, Eqn (7) results in the following compatibility relations:
,
du=---,
de pc
dr=---.
dT pc
(11)
(12)
The compatibility relations for the left-ward wave region (x < 0) can be achieved by simply changing the sign of c in the above two equations. 3. RIEMANN SOLVER Suppose the initial conditions at t = 0 to be two different homogeneous states. We denote the initial states by f_ in x < 0 and f+ in x > 0, respectively. When t > 0, there will arise two families of centered simple waves, one family propagating to the left-hand side, and another to the right-hand side. Depending on the initial values in the two wave regions, the wave speeds in the two wave families may adopt all values occurring in Eqn (10), or
466
X. LIN and J. BALLMANN
t Cs
Os
ca
% l,l c, o
FIG. 2.
A sketch in the (x, t) plane for a Riemann problem of combined stresses (a, ~).
only some of them, i.e. any one, two, three or four of the different kinds of waves can be absent 1"8]. As an example, one possible configuration of simple waves for certain given initial conditions is sketched in Fig. 2. Following a line x = const advancing in time, we can only cross the occurring waves in the inverse sequence as written in Eqn (10). This fact helps us to find a rule for the loading path in the stress space to be used in the Riemann solver, which will be discussed in the next section. Here, we want to get the solution f at the interface: x = 0, t > 0. The method applied here is well-known in gasdynamics from the shock tube problem, where the pressure and velocity are continuous at the contact interface. These conditions can be used to connect the regions of initially high and low pressures by an iterative process, where the velocity usually serves as a control variable, see e.g. [ 13]. Bearing in mind the assumption in continuum mechanics of solids (no dislocation), the stresses a, 3 and the velocities u, v are continuous at the interface. Therefore, we have to develop a method in the (a, 3, u, v) space. We prefer (u, v) as control variables. From Eqn (12), in the left-ward wave region where x / t = - c , the components of f satisfy the relation u=u_+
-jo_ pc
v=v_+
-_ pc
(13)
Generally, c in these equations is a function of (rr, 3), and the integral is dependent on the path in the stress space. So we cannot get a closed-form solution of the integrals. However, when in an iterative process an approximated final state (a, 3) is assumed, then by Eqn (11) together with (8) and (10), there exists one unique path in the stress space to connect this point with the initial state, see Section 4 for details. The integration of Eqn (11) can be carried out numerically in order to find this path for the integration of Eqns (13), which then can be integrated to get the iterated values of (u, v). In the right-ward wave region, we have to proceed in a similar way in order to calculate iterated values (u, v) at the right-hand side of the interface, using the same final values (a, 3). The related equations are u=u+-
,
v=v+-
j ~+ pc
--.
(14)
. pc
Since the interface points x = 0, t > 0 connect the left-ward wave region and the right-ward one, f will satisfy both Eqns (13) and (14). As in gasdynamics, if we draw a curve for rr as a function of u (and in the same way 3 as a function of v), the curve of Eqn (13) extended in the up-right direction, while that of Eqn (14) in the up-left direction. Then, the solution f is located at the crossing point of the two curves, see Fig. 3. To obtain directly the crossing point of two curves is difficult. In a linear tangential approximation the solution point will be approached by the crossing point of two straight tangential lines. Supposing fl = (al, 31, u 1, /)1) T to satisfy Eqn (13), i.e. ul=u_+
--, _
pc
vl=v_+
--, _
(15)
pc
according to Eqns (12), we take the tangential lines through this point in the (a, u)-plane
Riemann solver and Godunov method for solids
467
t7
12.
FIG. 3.
A sketch in the (o, u) plane for Riemann solver.
and in the (3, v)-plane, respectively. Then u-u
l=Al(a-al),
v-v
(16)
x=Bl(3-31),
where, A 1 is the slope of the a ( u ) curve and B 1 is that of the 3(v) curve. A 1 is sometimes different from BI, e.g. in the elastic case, A 1 = 1/(pCo), B1 = 1~(pc2). Therefore, we understand A 1 and BI as the following limiting values AI
= lim
1
A,,~OpC(a 1 + Aa, 31)
,
1
B1 = lira
~opc(al,31
+AQ
.
(17)
Supposing similarly fz to satisfy Eqns (14) in point 2 of Fig. 3, one gets the tangential approximations u-u2=
- A 2 ( c r - a2),
v - v2 = - B 2 ( z
- 32).
(18)
Then the crossing points of the tangential lines in the (~, u) and (3, v) planes can be solved from Eqns (16) and (18) as b = A l a l + /1"20"2 -j- u2 -- ul,
~ = 8131
A1+ A2
d- B 2 3 2 + v 2 - - 01
(19)
B1 + B 2
Substituting (6", ~) into Eqns (13) and (14), two new points (t~x, 01) and (~2, 02) can be attained. If I~1 - t~21 and 101 - 021 are small enough, the work is finished; otherwise, we can continue the iteration with the new points fl and i"2. In our computed examples, the iteration never failed when we simply started in the first step with the interaction of two elastic waves. Figure 3 is a sketch for the. above Riemann solver's procedure. Since the stresses a and 3 have different paths for loading and unloading, every time when (a, 3) is substituted into Eqns (13) and (14), we must find out a loading path, namely, the sound speed c = c(a, 3) along the lines from (a_, 3_) to (al, 31) and from (a+, 3+) to (a2, 32). In the following we discuss this problem. 4. L O A D I N G
PATH
Suppose (tri, T~) to be the initial state at a point before the next loading increment being taken under consideration. It is of course within the local yield elliptic surface xi
Denoting the final state after the loading increment by (a., 3c), the loading change is an elastic one if (ac, 3,) is also within the yield surface x~
+ 3o ~< ~ .
(21)
In this case, the loading amount has been generated by at most two kinds of waves, namely
468
X. LIN and J. BALLMANN "7"
" ' - . o°.~o*'" FIG. 4. Elastic loading path in the (or, r) plane.
.
e
"F .....
......
,
O"
**"
(7
",
......
(a)
(b)
(c)
FIo. 5. Three basic plastic loading paths in (a, z) plane.
c o and c2, and fa°do- = t roe--i - , ,]~r, pC
f~r'dz-ze-ti-
pc 0
, pC
(22)
pc 2
Since Co > c2, the path, see Fig. 4, can be expressed as (o-i, tl)
Co ' (o-c, ti)
c2 ~((re, t¢).
(23)
The problem for an elastic-plastic loading is more complicated because up to four kinds of waves with different sound speeds cs, c2, cf, Co are possible. Furthermore, cs and cr are dependent on the stresses (o-, t). To each of the different occurring waves corresponds a part of the loading path. Next, we will give three basic paths for the elastic-plastic loading case, and show later that all paths can be combined by these three paths. Owing t6 Eqn (1 l) we know that the elastic-plastic paths in the four quarters of (o-, T) are symmetrical to origin. So, without loss of generality, we may restrict the explanation by pictures to the first quadrant (a >/0, t/> 0). Suppose (a0, TO) is a point on the initial yield surface rl, and after loading, the final point ((re, t,) is outside the yield surface r~. Since the different speeds of sound satisfy the relation in Eqn (10), then, if c = c s appears, that will be the sound speed in the last part of the path ending at (o-~, z,). Therefore Eqn ( l l ) can be integrated backward from the point (g~, t~) by putting c = cs. This integration along a c~-path will end at a state named by (o.*, t*), which is placed at the initial yield surface x~ or at the crossing point with the o.-axis (t = 0) outside x~. There exist altogether three possible states for (o.*, t*), which lead to the three basic paths as follows. Path 1: [O'b[t> Io.*l, In this case, (o.*, t*) is on the yield surface. Then the path from (O'b, tO) to (O.*, t*) is expressed as (O.b' tb )
Co ) (O'*, tb)
c2 ) (O'*, r * )
~" ' (O',, t©),
(24)
as in Fig. 5(a). If Io.0l < Io-*l, two different paths named path 2 and path 3 will be possible. But in both cases we must integrate first through (O'b, t0) forward along a cf-path. Path 2: If there is a crossing point (o'k, tk) of the cf-path and the c~-path, as in Fig. 5(b),
Riemann solver and Godunov method for solids
469
the loading path is (ab, %)
' (ak, Zk)
~'
~' ' (ae, Z,).
(25)
P a t h 3: If the q-path crosses the line a = a~ at a value of [TI > I~1, then the c : p a t h is replaced by the line a = a~, which represents an unloading c2-path. The complete loading path is (ao, %)
c,
' (ak, rE)
~'
(26)
, ( a , , zc),
as in Fig. 5(c). Now we turn to discussing the general loading cases. We still denote the initial point by (a~, zl) and the final point by (a~, re). All the paths including plastic loading can be classified by two groups. (i) z~ and r e have the same sign, i.e. T;z~/> 0. Then point b on the initial yield surface can be determined as follows: % = z i,
ab = sign(a~)0x/~ - ~ .
(27)
The sound speed along the path from (a~, zi) to (ab, %) is C0. The paths after reaching (ab, %) can be determined as one of the basic paths discussed previously. (ii) q and % have a different sign, i.e. T : , < 0. In this case, we must first select a point (a,, z,) on the yield surface with z, = z i,
(28)
2 a a = sign(a~)O V/~i2 - %.
Then (ab, %) is determined by reflection of the point (a., r°) with the e-axis, % = --Za,
ab = a,.
(29)
In this group, the point b is not really reached during the wave process, but it may serve as an auxiliary point in order to reap benefits of the three basic paths examined previously. Regarding again the points (ab, %) and (a., ~.), if path 1 is applicable, the true path will be (a.~i)
co
,
(a*,Ti)
.
.
. . , (a*, :)
, ( % ~e).
(30)
If path 2 is applicable, the true path comes out as (al, r,)
Co
,
(%ra)
c,
,
(ak,--~)
.
.
.
.
, ( a k, rk)
, ( a , , x,),
(31)
where the q-path represents the mirror-inverted cf part of the basic path 2 at the a-axis. The continuation to point k is realized by a cz-path. If path 3 is to be selected, the true path is Gi, Ti)
Co
ct
' (aa, ~a)
' (%
--Tk)
cz
'
(32)
( % ~e).
Figures 6(a)-(c) represent the different cases of group (ii) discussed here. e
T
....
T
:,.-~ co c~
7
Cs/""
T
...... T....
a ',
..... t .... !.... "':ie ~
/ C2
',,
.
!
I
(a) FIG. 6.
(b)
(c)
Three true elastic-plastic loading paths for the case Zo~, < 0.
,; c2
e
470
X. LIN and J. BALLMANN
In gasdynamics, the waves to be treated in a Riemann problem are only single waves, either a shock or a rarefaction wave. There are only two possibilities which can be easily decided by comparing the pressure. In comparison with that, in elastic-plastic solids, up to four kinds of waves may occur, and the evaluation of the two wave configuration requires the solution of the accurate loading path connecting the different states. This complicates of course the logical structure of the computer program. Nevertheless, the loading paths described above make the evaluation possible. 5. THE SECOND-ORDER GODUNOV METHOD Suppose the x-axis is divided into a limited number of cells with equal lengths Ax and the cell centre is indicated by j. At the time level t = t", the functions a, 3, u, v, e, ? and in all cells are known and denoted by a~, z~, etc. These quantities are assumed to be constants in every cell. We want to calculate their new values at the time level t "+l Compared with gasdynamics, Eqns (1) have a similar conservation form, but contrary to the finite equations of state for a gas, Eqns (2) possess only an incremental formulation for changes of state• This results in many difficulties in numerical applications. For instance, let us consider the two-step L a x - W e n d r o f f scheme for Eqns (1) w.+½ j+½
1
At
.
~-~ ~ ( W ) "{'- W;+I) + 2----~x(f~j+, --
wj,+1 = w ; +
At
~),
.+, - .j_½,, e+L
(33)
+½ .+½ • • . where, • ~J+2.~ - f ( w , -~+2 ~). No. problems anse wRh . the calculation of the flux ~+2 in gas . . . J dynamtcs, since the equatton of state ts an algebram one. However, it ts dtfficult for ~qns (2), because there may be different yield stresses ~cj in cell-j and ~cj+t in cell-(j + 1), and similarly for the function H. A feasible way to overcome this difficulty is G o d u n o v ' s method [5], in which the flux at the cell interfaces is evaluated directly• Since we have different initial values in the different cells, for t > t" a Riemann problem is to be solved at every cell interface. From our presentation in the last two sections, the flux f can be worked out from the Riemann solution• Then, the Godunov scheme [5] for updating the values for the next time step is expressed as ,+1 At wj = w7 + ~x(fJ+½ - fj_½),
(34)
where the flux fj+½is denoted as the Riemann solution at t = t" + At~2 in the cell interface. Equation (34) was proved as a first-order accurate scheme. Before we present our method, let us see the most simple situation sketched in Fig. 7, where the distribution of the stress a is assumed as a solution of the Riemann problem at
,
L
) X
L
ax/~
,1.
ax/2
j+l
FIG. 7. A sketch for constructing a second-order Godunov method.
Riemann solverand Godunov method for solids
471
t -- t" + At~2. In the first-order Godunov method, fj+½is the value at the interface x = xj+½. It becomes evident from the profile in Fig. 7 that much information in the Riemann solution can be lost in this manner. Glimm [14] used a random number to select the solution point. Van Leer [13] and later Ben-Artzi and Falcovitz [15] achieved a second-order accuracy by reorganizing a linear distribution for the initial values, and adding a factor of Of/Or to the flux f. There are still many other Godunov-type methods, e.g. of Harten et al. [ 16]. Our work differs from that of the above authors in that the construction of the second-order scheme can be formulated directly from the Riemann solution. In our method, the flux fj+½ in Eqn (34) is replaced by its integral average in two half cells, i.e.
L½_±I x''' AxJx, fdx,
(t = t"+½).
(35)
Since the wave speed c is known having solved the Riemann problem, Eqn (35) can be rewritten in a more practical form. Assuming xj+½ as the origin, integrating Eqn (35) by parts, and noting x = - c A t / 2 for x < 0 and x = cat/2 for x > 0, we then have
At(f
fj÷~'
~+½ = ~ ( ~ + ~+') + 2-~x ,, j q
f'J*'c d f ) . cdf + or.'
(36,
The scheme Eqn (34) with the flux Eqn (36) will possess second-order accuracy in domains of smooth solution. In order to prove this, we use Eqn (7) to rewrite Eqn (36) in an approximate form: t - ~). ~+½ = ~1 ( ~ + ~+t) + 2-~x Aj+~(~+ -
(37)
Substituting Eqn (37) into (34) we would arrive at the original Lax-Wendroff scheme. The numerical procedure consists of three steps: the first one is to solve the Riemann problem, and to get the solution for the fan-shaped simple wave regions on both sides of the interface; the second one is to calculate the flux by Eqn (36); the final step is the updating of wj,+1 in Eqn (34). When ejn + l and ?jn + l are attained from Eqn (34), %n + l and Tjn + l can be updated by integration from Eqns (2) ej
- ej =
;
+
O2 J
dtr +
Htrzda+
77+' --),7 =
;
Htrz dr,
+HOaz 2 dz.
(38)
; Solving Eqns (38) for (a~ + 1, rT+ 1) is always difficult since we don't know the loading path in the j-th cell a priori. In order to evaluate the integrals, a simple wave path could be assumed, but the necessary iteration will result in a time-consuming numerical computation. In our proposal, the total increases of strain are divided into several small increments (fie, 6?) with the same factor of proportionality 6e = (e7 +' - e7)3).,
6)' = (?7 +1 - ?7)62,
(39)
where 62 > 0, Z6y = 1. Then the algebraic equations
6e =
+ 02 ]
+ Htrzt~,
6?=H(rz6~ + ( l + Ho2z2)6z,
(40)
are used for solving (&r, 6z), where, in every step, (a, z) in the coefficients are taken as the constants attained from the last step. The final stress will be the summation from all steps
trjn + l = ajn + Era,
~jn + l = ~jn + Erz.
(41)
472
X. LIN and J. BALLMANN
Since in the elastic region, H = 0, the first step for (6a, &r) should be taken to be large so that (a, ~) just reaches the yield surface.
6. TEST EXAMPLES In this section, we are going to present some test examples in order to demonstrate the efficiency of the proposed Riemann solver and the related second-order G o d u n o v method described in the foregoing sections. Since only little numerical results are available for one-dimensional combined stress problems, at first, we refer to the shock-tube problem which is a well-known example in gasdynamics. However, the symbols in gasdynamics are different from those of solid mechanics. We therefore remove this example to the Appendix. In this section, we consider the test problems from solid mechanics, and we will compare our numerical results to the exact analytical solutions. Suppose the tube is initially at rest, but subjected to a homogeneous initial stress distribution (fro, %). At time t = 0, a sudden shock-like load (a, z) is applied on the left side of the tube. The material parameters are chosen as p = 2.8 gm cm -3, Co = 6.39 km s - ~, c2 = 3.15 km s - ~, x o = 0.15 GPa, 0 = 2. The following are the results of the instantaneous distributions of the unknown functions from three test runs taken at time step N = 100. The computations were carried out with Ax = 1 and Courant number coAt/Ax = 1. The solid lines represent the exact solutions which can be attained by the methods in [ 8 - 1 2 ] , and the circles indicate numerical results from our computations. Test 1. (ao, %) = (0, 0.15), (a, r) = ( 1, 0.5). In this case, there are two kinds of waves in the simple wave region. The loading path is (0,0.15)
c, ,(0.32581,0.00096)
c' ,(1, 0.5).
(42)
It should be noted that in the region of the Q-wave r will decrease. An illustration of the loading path in the (a, r) plane and the space distributions of the seven unknown functions from our calculation is drawn in Fig. 8. Test 2. (ao, %) = (0.3, 0), (a, r) = (0.5, 1). In this case there are three kinds of waves in the simple wave region. The loading path is
(0.3,0)
.
. . O) . ,(0.27184,
. ,(0.27184, . 0.06345)
,(0.5, I).
(43)
Here, across the Co-wave a decreases. The results are shown in Fig. 9. Test 3. (a o, %) = (0, - 0 . 1 3 ) , (a, r) = (0.5, 0.5). In this case, there will be four kinds of waves in the simple wave region. The loading path is (0, - 0 . 1 3 )
~° , ( 0 . 1 4 9 6 7 , - 0 . 1 3 ) c,
,(0.29992,0.04733)
~' ,(0.29992, - 0 . 0 4 7 3 3 ) c' ) (0.5,0.5).
(44)
The results are shown in Fig. 10. It can be seen from the results that the calculations by the proposed second-order G o d u n o v method are basically in good agreement with the exact solutions, even for the n o n - m o n o t o n e functions. Of course, there are still some defects for regions with non-smooth solution, e.g. the rear part of the slow-wave, in which the numerical results propagate a little slower than the exact ones, and the strain ~, in Test 1 near the impact point has some defect, too. Still some effort is needed in order to achieve better agreement in such singular points.
Riemann solver and Godunov method for solids
0.55
1.1 0.9-
"o O_
0.45
..... THISWORK EXACT
0.7-
0.35
0.5
o
E
W .....OTHS IR f
0.25
0.3
•~ ~o
I--
-
-0. I
2b
4b
6'o
-0.05
8'o 1~o
2;
o
4b
6;
do "1;o
X (mm)
0.02
m,
0.00
0.00
-0.02 -0.04
-0.02
-0.06 -0.08-0. I0 -0.12-
-0.16-0.18
0
/
I'¢4
2;
.... THISWORK
..... THISWORK EXACT
-0.04
EXACT
-
-
-0.06 -0.08 -0.
4;
I0 -
-0.12
~o" 8~'1;o
2'o
0
4'o
6'o
8'o 1;o
X (mm)
X(mm)
0.095
0.II
0.085-
0.09-
0.0750.065-
c
EXACT
0.15
X (mm)
-0.14
-
0.05
0.I
E
473
..... THI5WORK
_90 0 . 0 5 5 ~n 0.045w 0.0350.025
EXACT
0.07 -
~
THIS WORK
0
E E
--
0.05 -
EXACT
0.03 -
0.015
o.01
0.005 -0.005
2'o
4'o
6'0
-o.01 810
'
) IO0
o
2'0
X(mm)
46
6'o
t~o
X (mm)
0.7-
0.6"5 o_ (...9 v 0 (7 ,7 0 ,.<,
LOAO ING PATH
0.5-
..... THI5 WORK - - EXACT
0.4-
T
'~.=.°.tsb.
0.3-
f
G
0.20.I
0
2'0
4'0
&
IOO
X(mm) FIG. 8.
Comparison of the exact solution with numerical results for combined longitudinal and torsional waves in a thin-walled tube (two waves case).
474
X. LIN and J. BALLMANN
0.55
I.I
L
o.so I
o~ I o.~o~
0.9
.....THIS WORK
\I
.~
__ ~,~c~
0.7..... THIS WORK
O_
°.~I
\
o I-
°.~°I
EXACT
--
0.s0.3 0.1 -0.1
0
20
40
60 X (mm)
80
2'o
,O0
•
6'0
4'o
i
8'o ,oo
X(mm) 0.05
0.005
I
~
-0.005
0.00
r
-0.05
-
-0.015
-0.025
IiOQI
-
~-0.035-
THIS WORK EXACT
-0.
IO-
--
40. IS >
-0.045-
-0.20
-0.055-
-0.25
-0.065-
-0.30
,..,
-0.075
o
A
~o
10
&
"
-0.35
|
4'o
A
,oo
0,07
•
0.04 -
u~ o.
0.03 -
0.20
.....
THIS WORK
o
--
EXACT
E E O (.9
0.02 O.Ol
o. Is
.....
THIS WORK
--
EXACT
O.)O
0.05 0.00
o
2'o
4'o
6'0
•
i
0'o ,oo
-0.05
J
"'
;o
2'o
•
io
'
1.0 LOADING
PRTH
(~,T)I
O.S ~L 0.6 0 C3.
.....
THIS
--
EXACT
HORK
0.4
0
=
-<:
Cs
0.15
/ ..........................
oo ....""
0.2
0.0
o
2'o
;o
6'o
8'0 ,Go
X (mm) F,G. 9.
i
8'o ,DO
X (mrr9
X(mm)
cz
i
,DO
O. 25 -
-
c: o
O.O0
~o
0.30 -
0.05 -
w
~o
X (mm)
X (mm)
0.06
THIS WORK EXACT
.....
E
Comparison of the exact solution with numerical results for combined longitudinal and torsional waves in a thin-walled tube (three waves case).
Riemann solver and Godunov method for solids
0.55
"S
0.5
I ,\
EL
o.351
o
0.25 -
475
.....
THIS WORK
--
EXACT
E
0,15-
"L
0.4 "6
EL I--
0,05
0.3 ""
0.2
WORK
THIS
0.1 0.0 -0.1
-o.os
2'o
4'o
6'o
8'0
•
-0.2
i
,oo
o
0.005 -0.005-0.015-
0.00
,r--------
-
-0.02 -0.04 -
--" THIS WORK - - EXACT
-0.045-
8'0 i~o
X (mm)
J
-0.025-J -0.035-
~o
4'0
~o
X (mrn)
E
-0.06
-
-0.08
-
.... THIS WORK - - EXRCT
v
>
-0.055-
-o. I0-
-0.12-
-0.065-0.075-
-0.14
-0.085
-0.16
o
2'o
4~
6'0 8'o i;o
2'o
4'o
lO0
X(mm)
X(mm)
0. 1 0-~ o 0.12
0 . 0 4 5 - ~ 0.035c
_q
0,025-
C~ LLI
0.015-
0.08-
THIS WORK EXRCT
E
•--- THIS WORK - - EXRCT
0.06-
E
0.040.02-
0.005-0.005
0.00-
o
~o
~o
4~
do
iGo
-0.02
o
~o
4;
X (mm)
6'o
8'0
i~o
X(mm)
0.6 0.5-
LOROINGTPRTH
"s EL L9 or-i 123. 0 v
'
0.4-
(o'.'r)/
THIS WORK EXACT
0.3:""
0.1
|
""
a
....................... .o.'cz I. c~c~ ?
0.2-
o
2'0
4'o
6'0
6'0
i&o
X(mm) Fm. 10. Comparison of the exact solution with numerical results for combined longitudinal and torsional waves in a thin-walled tube (four waves case).
476
X. L]N and J. BALLMANN 7. C O N C L U S I O N S
We can draw the following conclusions. (i) The construction ofa Riemann solver for elastic-plastic solids includes more difficulties than in gasdynamics, because there exist more components of stress and velocity in a sectional plane, which can undergo discontinuous changes. This resulted in complicated loading paths. Three basic elastic-plastic loading paths and a linear tangential approximation method are introduced as useful tools for the numerical analysis. (ii) The second-order Godunov method proposed in this paper can be used to deal with one-dimensional combined tangential and torsional stress waves in thin-walled tubes. The results of numerical tests are in considerable good agreement with the exact solutions. By some proper techniques we will enable this method for applications to more general two-dimensional problems. Therefore, this paper may have a finite significance in the impact dynamics of elastic-plastic solids. Acknowledoements--The authors wish to express their gratitude to the Alexander yon Humboldt Foundation for the research fellowship ofXiao Lin, and to the East China Institute of Technology for providing the opportunity to realize this cooperation.
REFERENCES 1. J. A. AaERSON, J. M. ANDERSON and W. W. KING, Dynamic analysis of cracked structures using singularity finite elements. In Mechanics of Fracture 4, Elastodynamic Crack Problems (edited by G. C. SIH), pp. 249-294. Noordhoff International Publishing, Leyden (1977). 2. H. ANTES, Anwendungen der Methode der Randelemente in der Elastodynamik und der Fluiddynamik. Matheraatische Methoden in der Technik 9 (edited by J. LEHN, H. NEUNZERTand H. WACKER). B. G. Teubner, Stuttgart (1988). 3. R.J. CLIFTON, A difference method for plane problems in dynamic elasticity. Q. appl. Math. 25, 97- l 16 (1967). 4. K.S. KIM, SpannungswellenanGrenzflfJcheninlinearelastischenScheiben. VDIVerlag, Reihe 18, Nr. 91 (1991). 5. S. K. GODUNOV, A finite difference method for the numerical computation of the discontinuous solutions of the equations of fluid dynamics. Mat. Sb. 47, 271-306 (1959). 6. K. NAGAYAMA,Solution of the high-pressure Rieman problem for solids including rigidity effects. J. Phys. Soc. Japan 58, 1631-1638 (1989). 7. J. A. TRANGENSTEINand R. B. PEMRER,The Riemann Problem for longitudinal motion in an elastic-plastic bar. Lawrence Livermore National Laboratory, UCRL-101214 (1989). 8. R. J. CLIftON, An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube. Proc. Fifth U.S. National Congress of Applied Mechanics, University of Minnesota, pp. 465-480 (1966). 9. T. C. T. TING, Elastic-plastic boundaries in the propagation of plane and cylindrical waves of combined stress. Q. appl. Math. 27, 441-449 (1970). 10. T. C. T. TING, Plastic wave propagation in linearly work-hardening materials. J. appl. Mech. 40, 1045-1049 (1973). l 1. T. C. T. Tn~G and NmG NAN, Plane waves due to combined compressive and shear stresses in a half space. J. appl. Mech. 36, 189-197 (1969). 12. R. P. GOEL and L. E. MALVERN,Biaxial plastic simple waves with combined kinematic and isotropic hardening. J. appl. Mech. 37, 1100- 1106 (1970). 13. B. VAN LEER, Towards the ultimate conservative difference scheme, V. A second-order sequel to Godunov's method. J. comp. Phys. 32, 101-136 (1979). 14. G. A. SOD, A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws. J. comp. Phys. 27, 1-31 (1978). 15. M. BEN-ARTZl and J. FALCOVITZ, A second-order Godunov-type scheme for compressible fluid dynamics. J. comp. Phys. 55, 1-32 (1984). 16. A. HARTEN, P. D. LAX and B. VAN LEER, On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Rev. 25, 35-61 (1983).
APPENDIX In order to compare the second-order Godunov method presented in this paper with other methods, we choose a well-known test problem from gasdynamics. The one-dimensional flow of an inviscid and non-heatconducting gas is described by the Euler equations in conservation form Ow
aF
--+--= Ot Ox
O,
(45)
Riemann solver and G o d u n o v method for solids
1.o
477
I.O--
exocL
thLs~;
.... ]-order :)'3 .v ..J .J m C 03 C3
to o
0.5 --
exoc
1 -onden LhLS
0.5.
'~;.. ~--~
work
o.o 0'.2
0'.4
0.0
0~.8 '
0'.6
f
i
0.2
i
0~.6
0.4
0.8
X
X 3.0
k~.%
O t._ D t~ tn m L (3-
.~.r
m
2.5
O3 L
0.5
O C
62
.... I - o r ' d e n • Lh~.s w o r k
•LhLw sork
13.0 0
0.2
2.0
0.6
0.4
,
,
l.S
0.8
0
0.2
X
"
0.6
0'.4
,
0.8
X
FIG. 11. C o m p a r i s o n of the second-order G o d u n o v method with the first-order method and the exact solution for the shock tube problem, C F L = 1, output at N = 40 time step, t = 0.185.
where
w =
,
F =
pu 2 + p
,
\(E + p)u/ with p the density, u the velocity, p the pressure and E the total energy per unit volume. Let e be the internal energy per unit mass, then for a perfect gas with a constant ratio of specific heats 7, P = (7 -
1)pe
= (y - l ) ( E - ½ p u 2 ) .
(46)
7 is a constant. When the initial values wT, F7 in every cell are given, the Riemann problem can be solved at the cell interfaces. The flux at the interfacej + ½ for the second-order G o d u n o v method is taken as an integral average over two half cells (suppose x j + ~ = 0):
lCJ+}
=
Ax
-~/2
F dx,
t = t" +
.
(47)
Since in the Riemann solution, we have determined the distribution of characteristic lines X = - =
(48)
U - - C, U, ;~ 9- C,
t
or shock lines X
= - = u - D, u + D,
(49}
t
in the fan-shaped simple wave region of the (x, 0-plane (where c, D are the sound and shock speed, respectively), we also know the function F(~} for the different ~. Then, Eqn (47) can be rewritten using integration by parts and putting x = ~ A t / 2 At
I F ; '~
Fj+~ = ~(F7 + F~+ i) - 2-~x jF:
~dF.
(50)
Equation (50) would be the same as Eqn (36), if u = 0 in Eqn (48). Having solved the fluxes, they will be used
478
X. LtN and J. BALLMANN
for the updating of w,
w~+! = w~
-~(rj+~pj_;). Ax 2
(5~)
The numerical test problem is taken from Sod [14], where the initialvalues are: u = 0; p ~- p = I for x < 0.5 and g = 0.125, p = 0.I for x > 0.5; ~,= 1.4; ~ = 0.01. Since any second-order method will produce oscillations near shock waves and aecxleration waves of finite amplitudes at the neighbourhood of x = 0.5, a better way is using a hybrid method which suitably combines a first-and a second-order method. O n e stronger point in our paper is that both the first-and the second-order method can be realized by Eqn (51), in which only F is needed to change. In the test example, the fluxes ~i+~ were taken to be first-order accurate in the two interfaces near the head of the acceleration wave at x = 0.5, and one to three interfaces in the shock wave region. Elsewhere, they are all taken as second-order fluxes. With the Riemann solution for the flux, we can use the Courant number CFL = I in the computations. Figure II shows the results for the spatial distribution of p, u, p and e after 40 time steps in comparison to the exact solution and to the results from the first-order G o d u n o v method. W e can see that in regions with continuous solution a significant improvement has been achieved by the second-order G o d u n o v method.
0734-743X/93 $6.00 + 0.00 © 1993 Pergamon Press Lid
A RIEMANN SOLVER AND A SECOND-ORDER GODUNOV METHOD FOR ELASTIC-PLASTIC WAVE PROPAGATION IN SOLIDS X. LIN a n d J. BALLMANN Lehr-und Forschungsgebiet fiir Mechanik der RWTH Aachen, Templergraben 64, 5100 Aachen, Germany
{Received 5 February 1991; and in revised form 27 March 1992)
Summary By use of three basic paths of elastic-plastic loading, a Riemann solver is established for the one-dimensional combined longitudinal and torsional stress wave problem of the thin-walled elastic-plastic tube. Based on this, a second-order Godunov method is developed for the numerical computation of elastic-plastic waves. Finally the numerical method is applied to some basic problems with known exact solutions and the comparison of results is made.
1. I N T R O D U C T I O N
Many solid materials behave plastically when they are subjected to high impulse loads. This may result in locally plastic wave propagation and propagating discontinuity surfaces like shocks, i.e. jumps of stresses and strains, or material interfaces. The wave structure, which has a great influence on the dynamic yield pattern and fracture processes, can only be correctly preserved by methods taking into account the probability occurring physical discontinuities. Unfortunately, most practical problems of mechanical waves in solids cannot be solved analytically due to the shape of the boundaries, stress hardening effects and complicated loading functions. During the last decades different numerical approaches have been introduced by numerous authors, involving finite element analysis [ 1], boundary element methods [2] and methods of characteristics [3,4]. The advantage of FEM is the easy implementation of nonlinearities of the PDEs and complicated boundaries, but it remains difficult to get the correct behaviour of the wave fronts. BEM is mainly restricted to linear problems. The method of characteristics truly integrates mathematically exact equations which have a direct physical interpretation in wave propagation, but this method has to consider discontinuities as interior boundaries within the solution domain, which are not easily managed for multi-dimensional problems. Therefore, it is necessary to have numerical methods for multi-dimensional wave problems which capture discontinuities and preserve simultaneously the principles governing the wave propagation phenomenon. In order to overcome the described difficulties, we take the idea of Godunov's method [5] in gasdynamics as an example, where a Riemann problem is the starting point of the method of solution. By this method, Nagayama [6] considered the behaviour of one-dimensional solids under high-pressure conditions, Trangenstein and Pember [7] dealt with the longitudinal wave motion in an elastic-plastic bar. In our case, the basic problem is the propagation of combined longitudinal and torsional waves in a thin-walled tube, where simple wave solutions can be used for the analysis of the wave pattern. The problem has been studied, e.g. by Clifton [8], Ting [9-11] and Goel and Malvern [12]. Their fundamental contributions formed the basic understanding for the treatment of waves involving elastic-plastic effects by simple wave analysis. In this paper, we will first present the construction of a Riemann solver for combined longitudinal and torsional stress waves in a thin-walled tube. In order to render the possibility of any arbitrary loading paths in stress space, three elementary paths will be introduced for the elastic-plastic loading history. Then we will deal with a second-order 463
464
X. LIN and J. BALLMANN
Godunov method. Finally, some numerical results will be given and compared with exact solutions. 2. BASIC GOVERNING EQUATIONS Consider a thin-walled tube being composed of an elastic-plastic stress-hardening material, which is subjected to a combined longitudinal and torsional shock loading, see Fig. 1. Following [ 8 - 1 2 ] for the derivation of the system of PDEs to be solved, we start with the local balance of momentum with components in longitudinal and circumferential direction, and the linearized strain-displacement relations partially derived with respect to time: Ow
Ot
Of
-
(1)
Ox'
where
w =
,
f=
.
u, tr, e shall represent the longitudinal particle velocity, stress and strain, and v, r, 7 the circumferential ones; p is the mass density, t is the time, and x is the distance measured along the axis of the tube. If the material of the tube is an isotropic work-hardening one, the incremental formulations of the material law have the following form
de=
+ 02 /I d~ + H~z dr'
dT= Hcrzda + (G+ HO2r2)dz,
(2)
where E is the Young's modulus and G the shear modulus. H is a nonlinear function of the yield stress x. The yield condition is given by
with 0 = v/3 for the yon Mises yield condition and 0 = 2 for the Tresca yield condition. The function H is defined by
where, g(tr) is the slope of the stress-strain curve in a simple tension test, i.e. v(~)
=
dtr --. dE
(5)
In the elastic region, g(Ox) = E, so H = 0, in the plastic region for the linearly isotropic work-hardening material, g(Ox) = EF (const).
FIG. 1. Thin-walled tube subjected to a sudden combined loading Or, r).
Riemann solver and Godunov method for solids
465
Substituting Eqns (2) into Eqn (1), the governing equations can be rewritten as AOf _ Of dt dx'
(6)
where
A =
0
0
p
0
0
0
p
-E + 0---T-
Hez
0
0
Hez
--+HO2z 2 G
0
0
1
He 2
1
For Eqn (6) a simple wave solution can be found, which connects two domains of different homogeneous states. Therein f is a function of x / t = c only. Equation (6) is then rewritten in the form (cA + I) df = 0,
(7)
where c is a free parameter. The necessary condition for a non-trivial solution df :/= 0 to exist, is llcA + Ill -- 0.
(8)
Solving this equation for c yields two solutions. In the elastic region (when H = 0), these are the characteristic wave speeds: Co=x/~,
c 2 = ~ p G-,
(9)
and in the plastic region, the two solutions are the wave speeds c~ and cf, which are correlated with the stresses (a, T) and satisfy 0 ~< cs ~< c2 ~< cf ~ Co.
(10)
The speeds c~ and cf are distinguished as slow and fast. The solutions of Eqn (8) contain the sign + for right-running waves and the sign - for left-running waves. In the right-ward simple wave region: x / t = c ( x > 0), where c may assume the values of Co, cf, c2 or c~, Eqn (7) results in the following compatibility relations:
,
du=---,
de pc
dr=---.
dT pc
(11)
(12)
The compatibility relations for the left-ward wave region (x < 0) can be achieved by simply changing the sign of c in the above two equations. 3. RIEMANN SOLVER Suppose the initial conditions at t = 0 to be two different homogeneous states. We denote the initial states by f_ in x < 0 and f+ in x > 0, respectively. When t > 0, there will arise two families of centered simple waves, one family propagating to the left-hand side, and another to the right-hand side. Depending on the initial values in the two wave regions, the wave speeds in the two wave families may adopt all values occurring in Eqn (10), or
466
X. LIN and J. BALLMANN
t Cs
Os
ca
% l,l c, o
FIG. 2.
A sketch in the (x, t) plane for a Riemann problem of combined stresses (a, ~).
only some of them, i.e. any one, two, three or four of the different kinds of waves can be absent 1"8]. As an example, one possible configuration of simple waves for certain given initial conditions is sketched in Fig. 2. Following a line x = const advancing in time, we can only cross the occurring waves in the inverse sequence as written in Eqn (10). This fact helps us to find a rule for the loading path in the stress space to be used in the Riemann solver, which will be discussed in the next section. Here, we want to get the solution f at the interface: x = 0, t > 0. The method applied here is well-known in gasdynamics from the shock tube problem, where the pressure and velocity are continuous at the contact interface. These conditions can be used to connect the regions of initially high and low pressures by an iterative process, where the velocity usually serves as a control variable, see e.g. [ 13]. Bearing in mind the assumption in continuum mechanics of solids (no dislocation), the stresses a, 3 and the velocities u, v are continuous at the interface. Therefore, we have to develop a method in the (a, 3, u, v) space. We prefer (u, v) as control variables. From Eqn (12), in the left-ward wave region where x / t = - c , the components of f satisfy the relation u=u_+
-jo_ pc
v=v_+
-_ pc
(13)
Generally, c in these equations is a function of (rr, 3), and the integral is dependent on the path in the stress space. So we cannot get a closed-form solution of the integrals. However, when in an iterative process an approximated final state (a, 3) is assumed, then by Eqn (11) together with (8) and (10), there exists one unique path in the stress space to connect this point with the initial state, see Section 4 for details. The integration of Eqn (11) can be carried out numerically in order to find this path for the integration of Eqns (13), which then can be integrated to get the iterated values of (u, v). In the right-ward wave region, we have to proceed in a similar way in order to calculate iterated values (u, v) at the right-hand side of the interface, using the same final values (a, 3). The related equations are u=u+-
,
v=v+-
j ~+ pc
--.
(14)
. pc
Since the interface points x = 0, t > 0 connect the left-ward wave region and the right-ward one, f will satisfy both Eqns (13) and (14). As in gasdynamics, if we draw a curve for rr as a function of u (and in the same way 3 as a function of v), the curve of Eqn (13) extended in the up-right direction, while that of Eqn (14) in the up-left direction. Then, the solution f is located at the crossing point of the two curves, see Fig. 3. To obtain directly the crossing point of two curves is difficult. In a linear tangential approximation the solution point will be approached by the crossing point of two straight tangential lines. Supposing fl = (al, 31, u 1, /)1) T to satisfy Eqn (13), i.e. ul=u_+
--, _
pc
vl=v_+
--, _
(15)
pc
according to Eqns (12), we take the tangential lines through this point in the (a, u)-plane
Riemann solver and Godunov method for solids
467
t7
12.
FIG. 3.
A sketch in the (o, u) plane for Riemann solver.
and in the (3, v)-plane, respectively. Then u-u
l=Al(a-al),
v-v
(16)
x=Bl(3-31),
where, A 1 is the slope of the a ( u ) curve and B 1 is that of the 3(v) curve. A 1 is sometimes different from BI, e.g. in the elastic case, A 1 = 1/(pCo), B1 = 1~(pc2). Therefore, we understand A 1 and BI as the following limiting values AI
= lim
1
A,,~OpC(a 1 + Aa, 31)
,
1
B1 = lira
~opc(al,31
+AQ
.
(17)
Supposing similarly fz to satisfy Eqns (14) in point 2 of Fig. 3, one gets the tangential approximations u-u2=
- A 2 ( c r - a2),
v - v2 = - B 2 ( z
- 32).
(18)
Then the crossing points of the tangential lines in the (~, u) and (3, v) planes can be solved from Eqns (16) and (18) as b = A l a l + /1"20"2 -j- u2 -- ul,
~ = 8131
A1+ A2
d- B 2 3 2 + v 2 - - 01
(19)
B1 + B 2
Substituting (6", ~) into Eqns (13) and (14), two new points (t~x, 01) and (~2, 02) can be attained. If I~1 - t~21 and 101 - 021 are small enough, the work is finished; otherwise, we can continue the iteration with the new points fl and i"2. In our computed examples, the iteration never failed when we simply started in the first step with the interaction of two elastic waves. Figure 3 is a sketch for the. above Riemann solver's procedure. Since the stresses a and 3 have different paths for loading and unloading, every time when (a, 3) is substituted into Eqns (13) and (14), we must find out a loading path, namely, the sound speed c = c(a, 3) along the lines from (a_, 3_) to (al, 31) and from (a+, 3+) to (a2, 32). In the following we discuss this problem. 4. L O A D I N G
PATH
Suppose (tri, T~) to be the initial state at a point before the next loading increment being taken under consideration. It is of course within the local yield elliptic surface xi
Denoting the final state after the loading increment by (a., 3c), the loading change is an elastic one if (ac, 3,) is also within the yield surface x~
+ 3o ~< ~ .
(21)
In this case, the loading amount has been generated by at most two kinds of waves, namely
468
X. LIN and J. BALLMANN "7"
" ' - . o°.~o*'" FIG. 4. Elastic loading path in the (or, r) plane.
.
e
"F .....
......
,
O"
**"
(7
",
......
(a)
(b)
(c)
FIo. 5. Three basic plastic loading paths in (a, z) plane.
c o and c2, and fa°do- = t roe--i - , ,]~r, pC
f~r'dz-ze-ti-
pc 0
, pC
(22)
pc 2
Since Co > c2, the path, see Fig. 4, can be expressed as (o-i, tl)
Co ' (o-c, ti)
c2 ~((re, t¢).
(23)
The problem for an elastic-plastic loading is more complicated because up to four kinds of waves with different sound speeds cs, c2, cf, Co are possible. Furthermore, cs and cr are dependent on the stresses (o-, t). To each of the different occurring waves corresponds a part of the loading path. Next, we will give three basic paths for the elastic-plastic loading case, and show later that all paths can be combined by these three paths. Owing t6 Eqn (1 l) we know that the elastic-plastic paths in the four quarters of (o-, T) are symmetrical to origin. So, without loss of generality, we may restrict the explanation by pictures to the first quadrant (a >/0, t/> 0). Suppose (a0, TO) is a point on the initial yield surface rl, and after loading, the final point ((re, t,) is outside the yield surface r~. Since the different speeds of sound satisfy the relation in Eqn (10), then, if c = c s appears, that will be the sound speed in the last part of the path ending at (o-~, z,). Therefore Eqn ( l l ) can be integrated backward from the point (g~, t~) by putting c = cs. This integration along a c~-path will end at a state named by (o.*, t*), which is placed at the initial yield surface x~ or at the crossing point with the o.-axis (t = 0) outside x~. There exist altogether three possible states for (o.*, t*), which lead to the three basic paths as follows. Path 1: [O'b[t> Io.*l, In this case, (o.*, t*) is on the yield surface. Then the path from (O'b, tO) to (O.*, t*) is expressed as (O.b' tb )
Co ) (O'*, tb)
c2 ) (O'*, r * )
~" ' (O',, t©),
(24)
as in Fig. 5(a). If Io.0l < Io-*l, two different paths named path 2 and path 3 will be possible. But in both cases we must integrate first through (O'b, t0) forward along a cf-path. Path 2: If there is a crossing point (o'k, tk) of the cf-path and the c~-path, as in Fig. 5(b),
Riemann solver and Godunov method for solids
469
the loading path is (ab, %)
' (ak, Zk)
~'
~' ' (ae, Z,).
(25)
P a t h 3: If the q-path crosses the line a = a~ at a value of [TI > I~1, then the c : p a t h is replaced by the line a = a~, which represents an unloading c2-path. The complete loading path is (ao, %)
c,
' (ak, rE)
~'
(26)
, ( a , , zc),
as in Fig. 5(c). Now we turn to discussing the general loading cases. We still denote the initial point by (a~, zl) and the final point by (a~, re). All the paths including plastic loading can be classified by two groups. (i) z~ and r e have the same sign, i.e. T;z~/> 0. Then point b on the initial yield surface can be determined as follows: % = z i,
ab = sign(a~)0x/~ - ~ .
(27)
The sound speed along the path from (a~, zi) to (ab, %) is C0. The paths after reaching (ab, %) can be determined as one of the basic paths discussed previously. (ii) q and % have a different sign, i.e. T : , < 0. In this case, we must first select a point (a,, z,) on the yield surface with z, = z i,
(28)
2 a a = sign(a~)O V/~i2 - %.
Then (ab, %) is determined by reflection of the point (a., r°) with the e-axis, % = --Za,
ab = a,.
(29)
In this group, the point b is not really reached during the wave process, but it may serve as an auxiliary point in order to reap benefits of the three basic paths examined previously. Regarding again the points (ab, %) and (a., ~.), if path 1 is applicable, the true path will be (a.~i)
co
,
(a*,Ti)
.
.
. . , (a*, :)
, ( % ~e).
(30)
If path 2 is applicable, the true path comes out as (al, r,)
Co
,
(%ra)
c,
,
(ak,--~)
.
.
.
.
, ( a k, rk)
, ( a , , x,),
(31)
where the q-path represents the mirror-inverted cf part of the basic path 2 at the a-axis. The continuation to point k is realized by a cz-path. If path 3 is to be selected, the true path is Gi, Ti)
Co
ct
' (aa, ~a)
' (%
--Tk)
cz
'
(32)
( % ~e).
Figures 6(a)-(c) represent the different cases of group (ii) discussed here. e
T
....
T
:,.-~ co c~
7
Cs/""
T
...... T....
a ',
..... t .... !.... "':ie ~
/ C2
',,
.
!
I
(a) FIG. 6.
(b)
(c)
Three true elastic-plastic loading paths for the case Zo~, < 0.
,; c2
e
470
X. LIN and J. BALLMANN
In gasdynamics, the waves to be treated in a Riemann problem are only single waves, either a shock or a rarefaction wave. There are only two possibilities which can be easily decided by comparing the pressure. In comparison with that, in elastic-plastic solids, up to four kinds of waves may occur, and the evaluation of the two wave configuration requires the solution of the accurate loading path connecting the different states. This complicates of course the logical structure of the computer program. Nevertheless, the loading paths described above make the evaluation possible. 5. THE SECOND-ORDER GODUNOV METHOD Suppose the x-axis is divided into a limited number of cells with equal lengths Ax and the cell centre is indicated by j. At the time level t = t", the functions a, 3, u, v, e, ? and in all cells are known and denoted by a~, z~, etc. These quantities are assumed to be constants in every cell. We want to calculate their new values at the time level t "+l Compared with gasdynamics, Eqns (1) have a similar conservation form, but contrary to the finite equations of state for a gas, Eqns (2) possess only an incremental formulation for changes of state• This results in many difficulties in numerical applications. For instance, let us consider the two-step L a x - W e n d r o f f scheme for Eqns (1) w.+½ j+½
1
At
.
~-~ ~ ( W ) "{'- W;+I) + 2----~x(f~j+, --
wj,+1 = w ; +
At
~),
.+, - .j_½,, e+L
(33)
+½ .+½ • • . where, • ~J+2.~ - f ( w , -~+2 ~). No. problems anse wRh . the calculation of the flux ~+2 in gas . . . J dynamtcs, since the equatton of state ts an algebram one. However, it ts dtfficult for ~qns (2), because there may be different yield stresses ~cj in cell-j and ~cj+t in cell-(j + 1), and similarly for the function H. A feasible way to overcome this difficulty is G o d u n o v ' s method [5], in which the flux at the cell interfaces is evaluated directly• Since we have different initial values in the different cells, for t > t" a Riemann problem is to be solved at every cell interface. From our presentation in the last two sections, the flux f can be worked out from the Riemann solution• Then, the Godunov scheme [5] for updating the values for the next time step is expressed as ,+1 At wj = w7 + ~x(fJ+½ - fj_½),
(34)
where the flux fj+½is denoted as the Riemann solution at t = t" + At~2 in the cell interface. Equation (34) was proved as a first-order accurate scheme. Before we present our method, let us see the most simple situation sketched in Fig. 7, where the distribution of the stress a is assumed as a solution of the Riemann problem at
,
L
) X
L
ax/~
,1.
ax/2
j+l
FIG. 7. A sketch for constructing a second-order Godunov method.
Riemann solverand Godunov method for solids
471
t -- t" + At~2. In the first-order Godunov method, fj+½is the value at the interface x = xj+½. It becomes evident from the profile in Fig. 7 that much information in the Riemann solution can be lost in this manner. Glimm [14] used a random number to select the solution point. Van Leer [13] and later Ben-Artzi and Falcovitz [15] achieved a second-order accuracy by reorganizing a linear distribution for the initial values, and adding a factor of Of/Or to the flux f. There are still many other Godunov-type methods, e.g. of Harten et al. [ 16]. Our work differs from that of the above authors in that the construction of the second-order scheme can be formulated directly from the Riemann solution. In our method, the flux fj+½ in Eqn (34) is replaced by its integral average in two half cells, i.e.
L½_±I x''' AxJx, fdx,
(t = t"+½).
(35)
Since the wave speed c is known having solved the Riemann problem, Eqn (35) can be rewritten in a more practical form. Assuming xj+½ as the origin, integrating Eqn (35) by parts, and noting x = - c A t / 2 for x < 0 and x = cat/2 for x > 0, we then have
At(f
fj÷~'
~+½ = ~ ( ~ + ~+') + 2-~x ,, j q
f'J*'c d f ) . cdf + or.'
(36,
The scheme Eqn (34) with the flux Eqn (36) will possess second-order accuracy in domains of smooth solution. In order to prove this, we use Eqn (7) to rewrite Eqn (36) in an approximate form: t - ~). ~+½ = ~1 ( ~ + ~+t) + 2-~x Aj+~(~+ -
(37)
Substituting Eqn (37) into (34) we would arrive at the original Lax-Wendroff scheme. The numerical procedure consists of three steps: the first one is to solve the Riemann problem, and to get the solution for the fan-shaped simple wave regions on both sides of the interface; the second one is to calculate the flux by Eqn (36); the final step is the updating of wj,+1 in Eqn (34). When ejn + l and ?jn + l are attained from Eqn (34), %n + l and Tjn + l can be updated by integration from Eqns (2) ej
- ej =
;
+
O2 J
dtr +
Htrzda+
77+' --),7 =
;
Htrz dr,
+HOaz 2 dz.
(38)
; Solving Eqns (38) for (a~ + 1, rT+ 1) is always difficult since we don't know the loading path in the j-th cell a priori. In order to evaluate the integrals, a simple wave path could be assumed, but the necessary iteration will result in a time-consuming numerical computation. In our proposal, the total increases of strain are divided into several small increments (fie, 6?) with the same factor of proportionality 6e = (e7 +' - e7)3).,
6)' = (?7 +1 - ?7)62,
(39)
where 62 > 0, Z6y = 1. Then the algebraic equations
6e =
+ 02 ]
+ Htrzt~,
6?=H(rz6~ + ( l + Ho2z2)6z,
(40)
are used for solving (&r, 6z), where, in every step, (a, z) in the coefficients are taken as the constants attained from the last step. The final stress will be the summation from all steps
trjn + l = ajn + Era,
~jn + l = ~jn + Erz.
(41)
472
X. LIN and J. BALLMANN
Since in the elastic region, H = 0, the first step for (6a, &r) should be taken to be large so that (a, ~) just reaches the yield surface.
6. TEST EXAMPLES In this section, we are going to present some test examples in order to demonstrate the efficiency of the proposed Riemann solver and the related second-order G o d u n o v method described in the foregoing sections. Since only little numerical results are available for one-dimensional combined stress problems, at first, we refer to the shock-tube problem which is a well-known example in gasdynamics. However, the symbols in gasdynamics are different from those of solid mechanics. We therefore remove this example to the Appendix. In this section, we consider the test problems from solid mechanics, and we will compare our numerical results to the exact analytical solutions. Suppose the tube is initially at rest, but subjected to a homogeneous initial stress distribution (fro, %). At time t = 0, a sudden shock-like load (a, z) is applied on the left side of the tube. The material parameters are chosen as p = 2.8 gm cm -3, Co = 6.39 km s - ~, c2 = 3.15 km s - ~, x o = 0.15 GPa, 0 = 2. The following are the results of the instantaneous distributions of the unknown functions from three test runs taken at time step N = 100. The computations were carried out with Ax = 1 and Courant number coAt/Ax = 1. The solid lines represent the exact solutions which can be attained by the methods in [ 8 - 1 2 ] , and the circles indicate numerical results from our computations. Test 1. (ao, %) = (0, 0.15), (a, r) = ( 1, 0.5). In this case, there are two kinds of waves in the simple wave region. The loading path is (0,0.15)
c, ,(0.32581,0.00096)
c' ,(1, 0.5).
(42)
It should be noted that in the region of the Q-wave r will decrease. An illustration of the loading path in the (a, r) plane and the space distributions of the seven unknown functions from our calculation is drawn in Fig. 8. Test 2. (ao, %) = (0.3, 0), (a, r) = (0.5, 1). In this case there are three kinds of waves in the simple wave region. The loading path is
(0.3,0)
.
. . O) . ,(0.27184,
. ,(0.27184, . 0.06345)
,(0.5, I).
(43)
Here, across the Co-wave a decreases. The results are shown in Fig. 9. Test 3. (a o, %) = (0, - 0 . 1 3 ) , (a, r) = (0.5, 0.5). In this case, there will be four kinds of waves in the simple wave region. The loading path is (0, - 0 . 1 3 )
~° , ( 0 . 1 4 9 6 7 , - 0 . 1 3 ) c,
,(0.29992,0.04733)
~' ,(0.29992, - 0 . 0 4 7 3 3 ) c' ) (0.5,0.5).
(44)
The results are shown in Fig. 10. It can be seen from the results that the calculations by the proposed second-order G o d u n o v method are basically in good agreement with the exact solutions, even for the n o n - m o n o t o n e functions. Of course, there are still some defects for regions with non-smooth solution, e.g. the rear part of the slow-wave, in which the numerical results propagate a little slower than the exact ones, and the strain ~, in Test 1 near the impact point has some defect, too. Still some effort is needed in order to achieve better agreement in such singular points.
Riemann solver and Godunov method for solids
0.55
1.1 0.9-
"o O_
0.45
..... THISWORK EXACT
0.7-
0.35
0.5
o
E
W .....OTHS IR f
0.25
0.3
•~ ~o
I--
-
-0. I
2b
4b
6'o
-0.05
8'o 1~o
2;
o
4b
6;
do "1;o
X (mm)
0.02
m,
0.00
0.00
-0.02 -0.04
-0.02
-0.06 -0.08-0. I0 -0.12-
-0.16-0.18
0
/
I'¢4
2;
.... THISWORK
..... THISWORK EXACT
-0.04
EXACT
-
-
-0.06 -0.08 -0.
4;
I0 -
-0.12
~o" 8~'1;o
2'o
0
4'o
6'o
8'o 1;o
X (mm)
X(mm)
0.095
0.II
0.085-
0.09-
0.0750.065-
c
EXACT
0.15
X (mm)
-0.14
-
0.05
0.I
E
473
..... THI5WORK
_90 0 . 0 5 5 ~n 0.045w 0.0350.025
EXACT
0.07 -
~
THIS WORK
0
E E
--
0.05 -
EXACT
0.03 -
0.015
o.01
0.005 -0.005
2'o
4'o
6'0
-o.01 810
'
) IO0
o
2'0
X(mm)
46
6'o
t~o
X (mm)
0.7-
0.6"5 o_ (...9 v 0 (7 ,7 0 ,.<,
LOAO ING PATH
0.5-
..... THI5 WORK - - EXACT
0.4-
T
'~.=.°.tsb.
0.3-
f
G
0.20.I
0
2'0
4'0
&
IOO
X(mm) FIG. 8.
Comparison of the exact solution with numerical results for combined longitudinal and torsional waves in a thin-walled tube (two waves case).
474
X. LIN and J. BALLMANN
0.55
I.I
L
o.so I
o~ I o.~o~
0.9
.....THIS WORK
\I
.~
__ ~,~c~
0.7..... THIS WORK
O_
°.~I
\
o I-
°.~°I
EXACT
--
0.s0.3 0.1 -0.1
0
20
40
60 X (mm)
80
2'o
,O0
•
6'0
4'o
i
8'o ,oo
X(mm) 0.05
0.005
I
~
-0.005
0.00
r
-0.05
-
-0.015
-0.025
IiOQI
-
~-0.035-
THIS WORK EXACT
-0.
IO-
--
40. IS >
-0.045-
-0.20
-0.055-
-0.25
-0.065-
-0.30
,..,
-0.075
o
A
~o
10
&
"
-0.35
|
4'o
A
,oo
0,07
•
0.04 -
u~ o.
0.03 -
0.20
.....
THIS WORK
o
--
EXACT
E E O (.9
0.02 O.Ol
o. Is
.....
THIS WORK
--
EXACT
O.)O
0.05 0.00
o
2'o
4'o
6'0
•
i
0'o ,oo
-0.05
J
"'
;o
2'o
•
io
'
1.0 LOADING
PRTH
(~,T)I
O.S ~L 0.6 0 C3.
.....
THIS
--
EXACT
HORK
0.4
0
=
-<:
Cs
0.15
/ ..........................
oo ....""
0.2
0.0
o
2'o
;o
6'o
8'0 ,Go
X (mm) F,G. 9.
i
8'o ,DO
X (mrr9
X(mm)
cz
i
,DO
O. 25 -
-
c: o
O.O0
~o
0.30 -
0.05 -
w
~o
X (mm)
X (mm)
0.06
THIS WORK EXACT
.....
E
Comparison of the exact solution with numerical results for combined longitudinal and torsional waves in a thin-walled tube (three waves case).
Riemann solver and Godunov method for solids
0.55
"S
0.5
I ,\
EL
o.351
o
0.25 -
475
.....
THIS WORK
--
EXACT
E
0,15-
"L
0.4 "6
EL I--
0,05
0.3 ""
0.2
WORK
THIS
0.1 0.0 -0.1
-o.os
2'o
4'o
6'o
8'0
•
-0.2
i
,oo
o
0.005 -0.005-0.015-
0.00
,r--------
-
-0.02 -0.04 -
--" THIS WORK - - EXACT
-0.045-
8'0 i~o
X (mm)
J
-0.025-J -0.035-
~o
4'0
~o
X (mrn)
E
-0.06
-
-0.08
-
.... THIS WORK - - EXRCT
v
>
-0.055-
-o. I0-
-0.12-
-0.065-0.075-
-0.14
-0.085
-0.16
o
2'o
4~
6'0 8'o i;o
2'o
4'o
lO0
X(mm)
X(mm)
0. 1 0-~ o 0.12
0 . 0 4 5 - ~ 0.035c
_q
0,025-
C~ LLI
0.015-
0.08-
THIS WORK EXRCT
E
•--- THIS WORK - - EXRCT
0.06-
E
0.040.02-
0.005-0.005
0.00-
o
~o
~o
4~
do
iGo
-0.02
o
~o
4;
X (mm)
6'o
8'0
i~o
X(mm)
0.6 0.5-
LOROINGTPRTH
"s EL L9 or-i 123. 0 v
'
0.4-
(o'.'r)/
THIS WORK EXACT
0.3:""
0.1
|
""
a
....................... .o.'cz I. c~c~ ?
0.2-
o
2'0
4'o
6'0
6'0
i&o
X(mm) Fm. 10. Comparison of the exact solution with numerical results for combined longitudinal and torsional waves in a thin-walled tube (four waves case).
476
X. L]N and J. BALLMANN 7. C O N C L U S I O N S
We can draw the following conclusions. (i) The construction ofa Riemann solver for elastic-plastic solids includes more difficulties than in gasdynamics, because there exist more components of stress and velocity in a sectional plane, which can undergo discontinuous changes. This resulted in complicated loading paths. Three basic elastic-plastic loading paths and a linear tangential approximation method are introduced as useful tools for the numerical analysis. (ii) The second-order Godunov method proposed in this paper can be used to deal with one-dimensional combined tangential and torsional stress waves in thin-walled tubes. The results of numerical tests are in considerable good agreement with the exact solutions. By some proper techniques we will enable this method for applications to more general two-dimensional problems. Therefore, this paper may have a finite significance in the impact dynamics of elastic-plastic solids. Acknowledoements--The authors wish to express their gratitude to the Alexander yon Humboldt Foundation for the research fellowship ofXiao Lin, and to the East China Institute of Technology for providing the opportunity to realize this cooperation.
REFERENCES 1. J. A. AaERSON, J. M. ANDERSON and W. W. KING, Dynamic analysis of cracked structures using singularity finite elements. In Mechanics of Fracture 4, Elastodynamic Crack Problems (edited by G. C. SIH), pp. 249-294. Noordhoff International Publishing, Leyden (1977). 2. H. ANTES, Anwendungen der Methode der Randelemente in der Elastodynamik und der Fluiddynamik. Matheraatische Methoden in der Technik 9 (edited by J. LEHN, H. NEUNZERTand H. WACKER). B. G. Teubner, Stuttgart (1988). 3. R.J. CLIFTON, A difference method for plane problems in dynamic elasticity. Q. appl. Math. 25, 97- l 16 (1967). 4. K.S. KIM, SpannungswellenanGrenzflfJcheninlinearelastischenScheiben. VDIVerlag, Reihe 18, Nr. 91 (1991). 5. S. K. GODUNOV, A finite difference method for the numerical computation of the discontinuous solutions of the equations of fluid dynamics. Mat. Sb. 47, 271-306 (1959). 6. K. NAGAYAMA,Solution of the high-pressure Rieman problem for solids including rigidity effects. J. Phys. Soc. Japan 58, 1631-1638 (1989). 7. J. A. TRANGENSTEINand R. B. PEMRER,The Riemann Problem for longitudinal motion in an elastic-plastic bar. Lawrence Livermore National Laboratory, UCRL-101214 (1989). 8. R. J. CLIftON, An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube. Proc. Fifth U.S. National Congress of Applied Mechanics, University of Minnesota, pp. 465-480 (1966). 9. T. C. T. TING, Elastic-plastic boundaries in the propagation of plane and cylindrical waves of combined stress. Q. appl. Math. 27, 441-449 (1970). 10. T. C. T. TING, Plastic wave propagation in linearly work-hardening materials. J. appl. Mech. 40, 1045-1049 (1973). l 1. T. C. T. Tn~G and NmG NAN, Plane waves due to combined compressive and shear stresses in a half space. J. appl. Mech. 36, 189-197 (1969). 12. R. P. GOEL and L. E. MALVERN,Biaxial plastic simple waves with combined kinematic and isotropic hardening. J. appl. Mech. 37, 1100- 1106 (1970). 13. B. VAN LEER, Towards the ultimate conservative difference scheme, V. A second-order sequel to Godunov's method. J. comp. Phys. 32, 101-136 (1979). 14. G. A. SOD, A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws. J. comp. Phys. 27, 1-31 (1978). 15. M. BEN-ARTZl and J. FALCOVITZ, A second-order Godunov-type scheme for compressible fluid dynamics. J. comp. Phys. 55, 1-32 (1984). 16. A. HARTEN, P. D. LAX and B. VAN LEER, On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Rev. 25, 35-61 (1983).
APPENDIX In order to compare the second-order Godunov method presented in this paper with other methods, we choose a well-known test problem from gasdynamics. The one-dimensional flow of an inviscid and non-heatconducting gas is described by the Euler equations in conservation form Ow
aF
--+--= Ot Ox
O,
(45)
Riemann solver and G o d u n o v method for solids
1.o
477
I.O--
exocL
thLs~;
.... ]-order :)'3 .v ..J .J m C 03 C3
to o
0.5 --
exoc
1 -onden LhLS
0.5.
'~;.. ~--~
work
o.o 0'.2
0'.4
0.0
0~.8 '
0'.6
f
i
0.2
i
0~.6
0.4
0.8
X
X 3.0
k~.%
O t._ D t~ tn m L (3-
.~.r
m
2.5
O3 L
0.5
O C
62
.... I - o r ' d e n • Lh~.s w o r k
•LhLw sork
13.0 0
0.2
2.0
0.6
0.4
,
,
l.S
0.8
0
0.2
X
"
0.6
0'.4
,
0.8
X
FIG. 11. C o m p a r i s o n of the second-order G o d u n o v method with the first-order method and the exact solution for the shock tube problem, C F L = 1, output at N = 40 time step, t = 0.185.
where
w =
,
F =
pu 2 + p
,
\(E + p)u/ with p the density, u the velocity, p the pressure and E the total energy per unit volume. Let e be the internal energy per unit mass, then for a perfect gas with a constant ratio of specific heats 7, P = (7 -
1)pe
= (y - l ) ( E - ½ p u 2 ) .
(46)
7 is a constant. When the initial values wT, F7 in every cell are given, the Riemann problem can be solved at the cell interfaces. The flux at the interfacej + ½ for the second-order G o d u n o v method is taken as an integral average over two half cells (suppose x j + ~ = 0):
lCJ+}
=
Ax
-~/2
F dx,
t = t" +
.
(47)
Since in the Riemann solution, we have determined the distribution of characteristic lines X = - =
(48)
U - - C, U, ;~ 9- C,
t
or shock lines X
= - = u - D, u + D,
(49}
t
in the fan-shaped simple wave region of the (x, 0-plane (where c, D are the sound and shock speed, respectively), we also know the function F(~} for the different ~. Then, Eqn (47) can be rewritten using integration by parts and putting x = ~ A t / 2 At
I F ; '~
Fj+~ = ~(F7 + F~+ i) - 2-~x jF:
~dF.
(50)
Equation (50) would be the same as Eqn (36), if u = 0 in Eqn (48). Having solved the fluxes, they will be used
478
X. LtN and J. BALLMANN
for the updating of w,
w~+! = w~
-~(rj+~pj_;). Ax 2
(5~)
The numerical test problem is taken from Sod [14], where the initialvalues are: u = 0; p ~- p = I for x < 0.5 and g = 0.125, p = 0.I for x > 0.5; ~,= 1.4; ~ = 0.01. Since any second-order method will produce oscillations near shock waves and aecxleration waves of finite amplitudes at the neighbourhood of x = 0.5, a better way is using a hybrid method which suitably combines a first-and a second-order method. O n e stronger point in our paper is that both the first-and the second-order method can be realized by Eqn (51), in which only F is needed to change. In the test example, the fluxes ~i+~ were taken to be first-order accurate in the two interfaces near the head of the acceleration wave at x = 0.5, and one to three interfaces in the shock wave region. Elsewhere, they are all taken as second-order fluxes. With the Riemann solution for the flux, we can use the Courant number CFL = I in the computations. Figure II shows the results for the spatial distribution of p, u, p and e after 40 time steps in comparison to the exact solution and to the results from the first-order G o d u n o v method. W e can see that in regions with continuous solution a significant improvement has been achieved by the second-order G o d u n o v method.