- Email: [email protected]
A numerical scheme for axisymmetric elastic waves in solids
ELSEVIER
Wave Motion 21 (1995)
115-126
A numerical scheme for axisymmetric elastic waves in solids * Xiao Lin *, Josef Ballmann L.ehr- und Forschungsgebietfiir
Mechanik der RWH Aachen, Templergraben 64. D-52062 Aachen, Germany Received 7 March 1994;
revised 29 August 1994
Abstract The paper considers the axisymmetric stress wave propagation in linear elastic solids which is governed by a system of hyperbolic PDEs with a source term. First, an explicit finite difference scheme is dealt with. Then problems related to the stress wave propagation and focussing in two-dimensional axisymmehic solids are studied using the obtained numerical scheme. Results are presented for the wave propagation in a half space loaded by an impact in a circular surface area, the dynamic stress intensity factor of a penny-shaped crack, the focussing of the von Schmidt wave in an impacted cylinder and the wave focussing in a hemispherically ended cylinder. For reasons of validation, comparisons are made with results gained by other authors. Some new improvements and explanations of the stress increases at cracks and in focussing zones could be attained.
1. Introduction Axis symmetrical stress wave propagation in elastic solids, which may be caused by an impact or another impulsive loading, is governed by a system of hyperbolic partial differential equations with a source term. Solutions of this system are of great practical interest for various initial and boundary conditions with cylindrical symmetry. Some analytical solutions were obtained, e.g. by Laturelle [ 1,2] for a half space using Laplace and Hankel transformations, and by Miklowitz [ 31 for a rod using the approximate Mindlin-Herrmann theory. These show that the source term complicates the analytical solving of the system in a way that does not occur in the plane problem. Therefore a numerical method giving a good approximate solution will be very useful for practical applications. Hirose and Achenbach [ 41 developed a time-domain boundary element method to study the elastic wave interaction in an axisymmetric body. For plane problems without a source term, Lin and Ballmann [5-81 have extended the Godunov-type characteristic-based finite difference methods of gasdynamics for stress waves in elastic-plastic solids and have obtained a great number of results. Problems with cylindrical symmetry can be treated in a similar way. *This paper is dedicated * Corresponding author.
to professor
Alan Jeffrey, University
of Newcastle
0165-2125/95/$09.50 @ 1995 Elsevier Science B.V. All rights reserved SSDIO165-2125(94)00046-8
Upon Tyne, on the occasion
of his 65th birthday
116
X. Lin, J. Balltnann/Wave
Motion 21 (1995) 115-126
Nevertheless, the existing numerical schemes dealing with hyperbolic systems with a source term seem not yet as well-developed as those for the systems without a source term. A widely used method is the time splitting technique which alternately solves a system of conservation laws without any source term and a system of ordinary differential equations modeling the source effect. However, it seems that this technique can produce misleading results, see Westenberger and Ballmann [ 91. For one-dimensional problems, some promising efforts were made by Glimm et al. [ lo], Glaz and Liu [ 111 and Roe [ 121. In this paper, we first propose an explicit finite difference scheme for the numerical integration of hyperbolic PDEs with a source term. Then, this scheme is applied to solve axisymmetric problems of elastic waves, e.g. the wave propagation in a half space due to an impact on a circular area of the boundary, the dynamic stress intensity factor at the tip of a penny-shaped crack and the focussing of stress waves caused by boundary effects. As far as possible the results will be compared with those of other authors in order to show the efficiency of the numerical scheme or to achieve a clearer explanation for the considered problem. 2. The governing equations and the numerical scheme Let r and z be the cylindrical coordinates and t the time. The elastodynamic equations for isotropic, linear elastic solids undergoing torsionless axisymmetric deformations can be written in the following form: au 5
au ar aa, --+x+;, 57ar
aa, 87 =ar+az+~,
1 da, --=PC; at
au ar +a$+&
--=i aa0 PC; at
au au cy$+ff-jy+;’
r
u
7
a24 au -- i aaz =a$+%+“;’ PC; at
u
au -- i a7 = av z+az. PC; at
(1)
where, u and v are the particle velocities in the directions of r and z, gr, a@, gz and r = arz are the stress components, cl and c2 are the longitudinal and transversal wave speeds, respectively. For isotropic materials, LY= I - 2(cz/ct )2 > 0. In order to study the numerical scheme, it is convenient to rewrite ( 1) in matrix form aw T=
af -+$+h. ar
(2)
Introducing Jacobian matrices A, B and C so that
f = Aw,
g = BW and h = Cw, (2) will take the form:
aw at
-=A$+B%+Cw.
(3)
Suppose the solution domain in the (r, z ) plane to be divided into rectangular cells and let w$ = w( nA t, (i i)Ar, (j - 4)A.r) be the values in the cell center (i,j) at the time level P. To obtain a numerical scheme, w is expanded into a Taylor series with respect to time up to the second order, w(t+At,r,z)
=w(t,r,z)
+Atf$+-2-x.
At2 a2w
If the time derivatives in (4) are eliminated by (2) or (3)) one obtains w(t+At,r,z)=w(t,r,z)+AtA;
+AtB&[w+~($+~+h)]
w+~ Z+$+h [ At(af
>I +AtC[w+$($+$+h)].
(5)
117
X. Lin, J. Ballmann/Wave Motion 21 (1995) 115-126
According to (5) and to the fact that quadratic cells in space are adequate bodies, a two-step explicit scheme is obtained as follows: 1 wIf+v2 r+1/2,j+1/2= -4
+ (wi)+ wY+t+l,j +a ( + wEj+l
gY+l.j+l
ti,j+l
in case of linear elastic isotropic
A
+ wY+l,j+l
> ( +
&+1,j - g:j) +
2
fY+l,j+l
-
y(h;+
fEj+l
hy+]qj
+
+
hrj+,
-
fY+l,j
+
ftj)
hy+,,j+l)’
(6) f"+w r+,,2,j_,,2)/2 and analogously for g, and A = At/Ar = At/Az. The scheme (6) represents an extension of ii was’s method [ 131 for the system of hyperbolic PDEs with a source term. In the numerical computation, the domain of the elastic body will be divided into a limited number of cells. The cell centers always lie inside. The boundary conditions are given at the grid points on the surface. Suppose the grid point (i + i, J + i) to be located at an upper boundary z = constant, where the stress components r = T and uz = Q are prescribed. Then the other components, i.e. U, o, ur, ~6, can be calculated using T and Q as well as the function values in the cells (i, J) and (i + 1, J). For example, the schemes for the o and (+r are follows:
wheref:+‘;/lj = (f
:‘:/&+I
+
2
(7) The conditions at other boundaries can be dealt with in a similar way, especially for the boundary r = 0, where u = 0, T = 0 are prescribed, and rl = Ar/2 is used to calculate the other components. In the second step of (6) one has r;_l/2 = 0 in h for the first column cells with the index i = 1. In this case h is treated in such a way, that the singularity l/r is avoided, e.g. approximating the term u/r as At!
r
= At-
U-O
(8)
r-0
The scheme admits everywhere
3. A half-space
the Courant-Friedrichs-Lewy
(CFL)
number cl A = 1.
problem
The first numerical example refers to a half space which is subjected to an impact over a circular area. This problem was considered by Laturelle [ 1,2] who derived the exact solution for the stress cz along the z-axis. We compare our results with those of Laturelle in order to examine the correctness of our numerical scheme.
X. Lin, J. Bahunn/Wave
118
Motion 21 (1995) 115-126
“O 1
-
Laturelle (1991)
t = 1.5
----
This work
r=O
0
1
Fig. 1. Elastic half-space subjected to an impact over a circular area. fig. 2. Comparison of the stress (TVgiven by the present numerical solutionand the analytical solution of Laturelle along the z-axis for the time level t = 1.5.
The material constants of the elastic body are taken as p = 1, ct = 1, c2 = l/d (which corresponds to v = 0.25 for Possion’s ratio). The radius of the loaded arca, see Fig. 1, is set to a = 1, which is divided into 400 cells, i.e., Ar = hr. = At = l/400. The loading is a shock with the amplitude q = 1 and the acting duration of 0.5. The flux components (T, on the boundary z = 0 are set to -q(t) for the grid points located in 0 5 r < 1, while (T, = -q(t) /2 is applied to the grid point on I = 1. The numerical results for (+z along the axis of symmetry for the time levels f = 1.5 and 2 are plotted in Fig. 2 and Fig. 3 in comparison with the exact solution of Laturelle. One can see that the two solutions agree well in the most smooth regions. A strong defect is found near z = 0.5692 at r = 2, because this is a meeting point of all shear wave fronts, in which the stress is undetermined, see Fig. 1 of Ref. [2]. The shock waves appearing in z = 1.108 at t = 1.5, or z = 1.732 and z = 1.115 at t = 2 represent the meeting points of the longitudinal wave fronts. These shocks are modeled over a narrow region, since they move along the axis with a speed c < ct. However, the longitudinal loading and unloading shock fronts with speed ct were calculated correctly. 4. A penny-shaped
crack subjected
to a tension shock
The dynamic stress intensity factor at the tip of a penny-shaped crack was calculated by Chen and Sih [ 141 using integral transformation techniques, and by SMdek and Sl&lek [ 151 solving the boundary integro-differential equations, and by Hirose and Achenbach [4] using a time-domain boundary element method. Zhang and Gross [ 161 have also applied a time-domain boundary element method to this problem and obtained the same results as Ref. [ 41. It is noteworthy that the results from the three methods show remarkable differences. The considered problem is shown in Fig. 4a, where two plane longitudinal tension shock waves travel towards the penny-shaped crack from opposite sides. Because of the symmetry, it is sufficient to solve the equations only in the first quadrant, see Fig. 4b. In our calculation the number of cells along the edge of the crack has been chosen as Zc = 100 and those of the total computational domain are It x Jt = 910 x 810. In normal direction on the line z = 0 the boundary condition is specified by gz = 0 on the crack edge and by v = 0 on the outer part, while in the tangential direction the shear stress is r = 0 along the whole line. The boundary conditions are satisfied by using these values for the flux computation in the boundary grid points. The conditions on the outer boundary of the computational domain can be dealt with in a simple way. Taking the boundary i = Ii as an example, the function values in the cell (It, j) are set to those of the cell (It - 1, j) in every time step. By this approach a plane wave can be kept until the disturbance from the crack tip comes there. Therefore, it will need 2 x 810 time steps for the disturbance wave to take a round trip from the crack tip to boundary and then back to the crack tip, which is sufficient to avoid any influence of the outer boundary on the crack tip during
119
X. Lin, J. Ballmann/Wave Motion 21 (1995) 115-126
I.”
-
Laturelle
t=2
(1991)
,4-z
I_
0,
Cl
-
y
L-l m
J
4-‘-c-‘T-’
Q 7-v
f.
&lo
2
=
‘-‘F
0
3
Fig. 3. Comparison of the stress gL given by the present numerical solution and the analytical solution of Laturelle along the z-axis for t = 2.
the time level
Fig. 4. Sketches for the crack problem. (a) Physical problem, (b) Zoning for the calculation.
the time interval 0 5 czt/a 5 9. The material constants are chosen as cl = 1, c2 = l/fi. The initial conditions were two homogeneous mirror symmetrical states with respect to the line z = 0 and characterized by following values
go=uo= 0 L
0: = a; = -
uo= 70 =
41
v
4”=0.5.
0
l-v4’
0,
(9)
The dynamic stress intensity factor can be obtained in every time step by the following relations (see Fig. 5) (10)
ao z=-zE
4
-312
34 cos -1
2
ao
4 _ -_-E-312
az-
&
sin Z9
2’
(11)
where the values of a@/& and a@/& can be taken in the points C, D and E in order to calculate KI. The results for the normalized dynamic stress intensity factor are shown in Fig. 6 and compared with those of the above-mentioned authors. Over wide time ranges, our results agree well with those of Hirose and Achenbach [4]. But a considerable difference appears for the peak value and its phase position. Any of the other methods shows the peak position clearly at a normalized time less than @/a = 2, while our method gives the peak position at czt/a = 6, which corresponds exactly with the time that a Rayleigh wave needs to propagate along the crack edge from the crack tip to the center point r = 0 and then back to the crack tip. For the discussed material v = 0.25, Q/CR = 1.08766. The result of Chen and Sih [ 141 is based on a material with v = 0.29, Q/CR. = 1.0801. Therefore, there is also 6 > 2. According to the basic analyses for the plane crack by Thau-Lu [ 171, Kim [ 181 and Lin-Ballmann [ 51, the effect of the Rayleigh wave reflected at the crack tip can be assumed as responsible for the steep decrease of the SIF from its first maximum. This explanation supports our present result.
5. The initial wave pattern in an impacted circular md
Impacts on cylindrical rods have been analyzed in many experiments. One typical example is the split Hopkinson bar. In practice, it is mostly assumed that the stress wave propagation inside the rod is a homogeneous plane wave. The wave dispersion due to two-dimensional effects is then neglected. In order to analyse this
X. Lin. J. BallmandWave
120
: -r
Motion 21 (1995) 115-126
lJ.y-z i-2-a
c
--._____._._.-.-.-
KI
K;t
0.5
1 A crack face
’
O.Oi~,,,~,,,,,,,,,,,, 0 1 2
0
This work Hirose and Achenbach (1989) Slddek and Slidek (1986) Chen and Sih (1977) 3
4
5
6
7
8
c2tla Fig. 5. A sketch for the treatment of the crack tip
neighbouring cells.
Fig. 6. Normalized Mode I dynamic stress intensity factor at the tip of a penny-shaped crack under Heaviside tension wave.
phenomenon, the integral transformation technique was applied to the problem by Pochhammer and Chree (see Ref. [ 191) early in last century. They obtained a frequency equation for the approximated calculation of the wave speeds up to the second term of the series. Therefore, the wave dispersion has been taken into account in a homogenized manner. For a practical application with given boundary and initial conditions, Mindlin and Herrmann have proposed an approximate one-dimensional theory, which has been used later by Miklowitz [ 31 to obtain some solutions. The wave propagation in an impacted cylindrical rod can be calculated quite correctly by the numerical method proposed in this paper. In our numerical calculation, the rod has a semi-infinite length and occupies the spatial domain 0 < r 5 1, 0 5 z 5 00. The assumed material constants are ct = 1, c2 = l/d, see Fig. 7. The rod which was originally at rest and stress free is subjected to a tension shock on the limiting surface z = 0 at the time t = 0. The loading condition can be expressed by a Heaviside function gL (t, r, 0) = H(t), which implies that it is uniformly distributed over the surface after the time t > 0. The radius a = 1 of the rod is divided into 50 cells and the number of cells along the z-axis is taken according to the upper time level of the computation with CZ=l. The calculated distributions of the stress components ur and crZ in the rod at time t = 3 are shown in the first part of Fig. 8. The plotted values represent the data in the cell centers. At t = 3 the plane wave front due to the impact has just reached z = 3. The stress u, behind this wave front does not vanish due to the transient plane strain state induced by the impact. Therefore, a von Schmidt wave is produced in the crossing point of the shock wave front z - ct t = 0 and the side boundary r = 1 in order to satisfy the boundary condition g’r = 0 on the surface r = 1. This von Schmidt wave causes also a negative velocity u in the r-direction, which can change the stress state again when the wave propagates towards the axis of the rod, since we have u = 0 along the axis because of symmetry. As a result, a high stress intensity will appear on the axis due to the focussing of the von Schmidt wave. The stress (+Bhas also a focussing peak at the same time and its distribution looks similarly to that of a,; therefore it is not drawn in the picture. In the focussing region, ur and a@ are negative, but uZ is positive. Then the von Mises yield stress K, defined by the second invariant of the deviatoric stress tensor
K2=+.-uz)2+(us -
uz J2+ (ur - ud2] + 72,
(12)
will also have a focussing peak on the axis, which is shown in Fig. 8, too. Theoretically, the focussing point on the axis is a singular point of the elastic solution, where the stresses become infinite. Indeed, the values of the focussing peak become higher if a finer grid is used in the numerical
X. Lin, J. Ballmann/Wave Motion 21 (1995) 115-126
Fig. 7. A circular rod subjected to a tension shock.
max
=
2.06
Fig. 8. The stress distributions in an impacted rod at time t = 3.
121
122
X. Lin, .I. Ballmann/Wave Motion 21 (1995) 115-126
max Hg. 9. The time history of the stress uz (0,
= 2.15
t) along the axis.
2,
2.5 t/a
L 2.0 -
-.-.ii
z/a = 4 (this work)
_____, z/a = 2 (Miklowitz,
1957)
z/a = 4 (Miklowitz,
1957)
1.5 -
I
0.0 0
1 2
I
= 2 (this work)
I 4
I
I 5
I
I 8
I 10
(cd - z)/a Fig. 10. The time history of stress uL at discrete axial points, compared with the approximate results of Miklowitz.
can occur in a region near the axis. The focussing peak appears for the first time when the von Schmidt wave has reached the axis of the rod. Figure 9 shows the time history of the distribution of gZ in the first row of grids along the axis. It can be seen that the maximum of the peak appears at f = 3.54. Its value is cf, = 2.15 at about z = 1.94. Then the peak decreases. At about t = 7, a second focussing peak begins to grow. The wave train formed by the leading tension wave and the von Schmidt wave focussing at the axis of the rod and being reflected at the boundary forming a new tension front maintains always its typical geometry. Figure 10 presents the time histories of the stress CT~ in the fixed axial points z/a = 2 and 4 in order to compare the present numerical results with the analytical solutions found by Miklowitz [ 31 using the approximate Mindlin-Herrmann theory. The input data were chosen as cl = 1 and c2 = 0.5821 in accordance with Possion’s ratio v = 0.24375 used by Miklowitz; the geometric parameters were chosen the same as above. Obviously, there are large differences between the two results, especially with respect to the focussing peak behaviour. The example clarifies the defects of the one-dimensional Mindlin-Hemnann theory in the leading part of the wave over a distance of about 10 radii behind the head of the wave.
X. Lin, J. Balhann/Wave Motion 21 (1995) 135-126
123
Fig. Il. A sketch of a cylindrical bar with a spherical free end.
grid 2 0
r a=1
Fig.
L
12. A sketch of the combination of two grids. The data are exchanged in every time step from grid 1 to grid 2 on those cells near the boundary denoted by X, and from grid 2 to grid 1 on those cells near the boundary denoted by o.
6. The focus of a stress wave in a hemisphere In this last example we consider wave focussing in a cylindrical bar with a spherical free end, see Fig. 11. When a shock wave is induced into the bar at the plane left end, the stress wave will propagate into the cylinder as described in the last section. Since the surface of the hemisphere at the right end is free, the longitudinal stress wave will be first reflected with wave splitting and phase changes and then focused on the axis at a point with the distance of half a radius from the end. A problem similar to that but under plane stress condition has been considered by Ballmann et al. [20] and Niethammer et al. [21]. This calculation is a test for the axisymmetric problem. The hemispherical surface represents a curved boundary with a boundary condition that the rectangular cells used in previous sections fail to satisfy correctly. Of course one could use triangles with edges to fit the boundary geometry, but the discretization of the governing equations on an irregular grid will cause too much numerical dispersion and dissipation. Therefore we apply an overlapping grid technique [ 211 to compute the problem. Suppose the system is still axisymmetric on the z-axis. Then we choose one grid still rectangular in r- and z-direction, while the second grid has to fit the surface, i.e. in the present case we use spherical coordinates (R, p). The arrangement of the two grids is shown in Fig. 12. The governing equations of elastodynamics in spherical coordinates can be found in Ref. [ 221. After some changes of symbols and coordinate directions in accordance to Fig. 12, we may rewrite them for the torsionless case with cylindrical symmetry as follows,
124
X. Lin, J. Balhann/Wave Motion 21 (1995) 115-126
case with cylindrical symmetry as follows, 801 =!%+L!%+ GR asp 1 au1 --=PC; at
2ai -u2-us-ri2tanSp
au1 aR+Zau2+ R a9
au1 aav2 --1 au3 PC; at =“=+w+
&I2 =!Z+L!%+
R
’
2avi - (Yvztan (p R
’
R
3712 - (~2 - (~3) tmqo R
R a5o
I 1av2+u+4vl-~v2t~9 --1 aa2 _avl & at aR R R ap
-v2mp
(~+cu)v~
5-
’
1 aTl2 au2 i avl --= zi+ii&zT PC; at
v2
(13)
where vi and vz are the velocity components in R- and qP-direction, respectively, and ui = (TR,02 = a,+,,us = a0 and 7-12= URq are the stress components. It is evident that this system can be integrated numerically by the scheme (6). The radius of the cylindrical rod is chosen as a = 1. For grid 1 we take 50 cells in r-direction and 100 cells in z-direction (including 50 cells for the hemisphere). The radius of the hemisphere is equal to a. The number of cells of grid 2 in R- and p-direction are set to 11 and 59, respectively, which gives an inner radius of grid 2 of about 0.746. The smallest cell’s length AR of grid 2 is about l/50 in order to obtain approximately the same CFL number as in grid 1. The computations are carried out in the two grids separately, but the data are exchanged in every time step by linear interpolation from grid 1 to grid 2 in the cells near the boundary denoted by x and from grid 2 to grid 1 in the cells near the boundary denoted by o. The material data are The boundary condition is given by uZ = -H(t) at the rod’s plane end chosen such that ct = 1 and cz = l/fi. surface z, = 0. The computation has been carried out to analyze the stress focussing in the hemisphere, which is to be expected by ray arguments at R = a/2 on the axis of the hemisphere. The highest focussing intensity is found at time step N = 122, i.e., tit/a = 2.44. The momentaneous distributions of the stresses ue, u, and u, at this time are presented in Fig. 13. The focussing results shown in Fig. 13 are due to the reflection of a plane shock wave front from the spherical boundary. There exists also a negative peak with the intensity uZ = -1.76 in this picture, which denotes the focus of the von Schmidt wave discussed in the last section. This peak is still moving in z-direction at this time. We have continued the computation until the negative peak was reflected from the hemisphere’s surface in order to see whether another focussing would occur. The result showed indeed another occurrence of focussing, but the stress increase is a lot weaker as presented in Fig. 13 for a, and uL. That is, because the focussed part of the von Schmidt wave is much smaller than the one of the plane shock wave. From Fig. 13 it can be seen that all stress components ue, a, and uZ are positive in the hemisphere’s focussing area. The shear stress T vanishes along the axis because of symmetry. Therefore, according to ( 12), the von Mises yield stress will not show a focussing peak in the discussed case.
7. Conclusions In this paper a finite difference scheme was proposed for the numerical integration of systems of hyperbolic PDEs with a source term, which is typical for cylindrical or spherical geometries. The elastic wave propagation in two-dimensional axisymmetric solids can be handled with the proposed numerical scheme. Several results have been obtained, which show the efficiency of the method compared to analytical methods. For the pennyshaped crack and for cylindrical and spherical geometries some new explanations of stress increases and focussing effects were presented.
X. Lin, J. BallmmdWave
max
=
2.30
max
=
2.37
max
=
2.11
min
= -1.76
Motion 21 (1995) 115-126
Fig. 13. Stress focussing of a plane wave in the hemisphericalend of a cylinder.
125
X. Lin, J. Balbnann/Wave Motion 21 (1995) 115-126
126
Acknowledgements
The research work presented in this paper was supported by the Deutsche Forschungsgemeinschaft under Grant No. Ba 661/ 12-1, which is gratefully acknowledged. References [ 1I PG. Latutelle, “The stresses produced in an elastic half-space by a normal step loading over a circular area. Analytical and numerical results”, Wave Motion 12, 107-127 (1990). 121 EG. Latutelle, ‘The stresses produced in an elastic half-space by a pressure pulse applied uniformly over a circular area: role of the pulse duration”, Wave Motion 14, l-9 (1991). [ 31 J. Miklowitx, “The propagation of compressional waves in a dispersive elastic rod, part 1 - results from the theory”, J. Appl. Mech. 24.231-239 (1957). [4] S. Hirose and J.D. Achenbach, ‘Time-domain boundary element analysis of elastic wave interaction with a crack”, Int. J. Numer. Methods Eng. 28.629-644 ( 1989). [ 51 X. Lin and J. Ballmann, “Reconsideration of Chen’s problem by finite difference method”, Eng. Fracture Mech. 44.735-739 (1993). [ 61 X. Lin and J. Ballmann, “Numerical method for elastic-plastic waves in cracked solids, part 1: anti-plane shear problem”, Arch. Appl. Mech. 63,261-282 (1993). [7] X. Lin and J. Ballmann, “Numerical method for elastic-plastic waves in cracked solids, part 2: plane strain problem”, Arch. Appl. Mech. 63,283-295 ( 1993). [ 81 X. Lin and J. Ballmann, “Improved bicharacteristic schemes for two-dimensional elastodynamic equations”, to appear in Quart. Appl. Math. (1995). [9] H. Westenberger and J. BaJhnann, “The homogeneous homentropic compression or expansion - a test case for analyzing Sod’s operator-splitting “, in J. Ballmann and R. Jeltsch (eds.): Nonlinear Hyperbolic Equations - Theory, Computation Methods, and Applications, 678-687, Vieweg, Braunschweig/Wiesbaden ( 1988). [lo] J. Glimm, G. Marshall and B. Plohr, “A general&d Riemann problem for quasi-one-dimensiotud gas flows”, Adv. Appl. Math. 5,
l-30 (1984).
[ 111 H.M. Glaz and T.-P Liu, “The asymptotic analysis of wave interactions and numerical calculations of transonic nozzle flow”, Adv. Appl. Math. 5, 11 l-146 (1984).
[ 121 P.L. Roe, “Upwind differencing schemes for hyperbolic conservation laws with source terms”, in: C. Carasso, P-A. Raviart and D. Serre (Eds.), Nonlinear Hyperbolic Problems, 41-51, Springer, Berlin (1986). [ 131 B. Eilon, D. Gottlieb and G. Zwas, “Numerical stabilizers and computing time for second-order accurate schemes”, J. Comput. Phys. 9, 387-397 (1972). [ 141 E.F? Chen and G.C. Sih, “Transient response of cracks to impact loads”, in G.C. Sih (cd.): Mechanics of Fracture 4 - Elustodynamic Crack Problems, l-58, Noordhoff, Leyden (1977). [ 151 J. SlPdek and V. Sladek, “Dynamic stress intensity factors studied by boundary integro-differential equations”, ht. J. Numer. Meth. Eng. 23.919-928 (1986). [ 161 Ch. Zhang and D. Gross, “Transient elastodynamic analysis of a penny-shaped crack”, Eng. Fracture Mech. 46.641-654 (1993). [ 171 S.A. Thau and T.H. Lu, “Transient stress intensity factors for a finite crack in an elastic solid caused by a dilatational wave”, ht. J.
Solids Strucmres 7, 731-750 (1971). [ 181 K.S. Kim, “Dynamic propagation of a finite crack”, ht. J. Solids Structure 15,685-699 (1979). [ 191 A.E.H. Love, A Treatise on the Mathematical Theory of Elasticity, 4th ed., 287-292, Dover, New York (1944). [ 201 J. Ballmann, H.J. Raatschen and M. Staat, “High sttess intensities in focussing zones of waves”, in: P Ladeveze (Ed.) : Local Effects in the Analysis ofStructures, 235-252, Elsevier Science Publishers, Amsterdam (1985). [21] R.J. Niethammer, K.-S. Kim and J. Ballmann, “Numerical simulation of shock waves in linear-elastic plates with curvilinear boundaries and material interfaces”, submitted to Int. J. Impact Eng.. [ 221 K.P. Graff, Wave Motion in Elastic Solids, 601-602, Clarendon Press, Oxford ( 1975),
Wave Motion 21 (1995)
115-126
A numerical scheme for axisymmetric elastic waves in solids * Xiao Lin *, Josef Ballmann L.ehr- und Forschungsgebietfiir
Mechanik der RWH Aachen, Templergraben 64. D-52062 Aachen, Germany Received 7 March 1994;
revised 29 August 1994
Abstract The paper considers the axisymmetric stress wave propagation in linear elastic solids which is governed by a system of hyperbolic PDEs with a source term. First, an explicit finite difference scheme is dealt with. Then problems related to the stress wave propagation and focussing in two-dimensional axisymmehic solids are studied using the obtained numerical scheme. Results are presented for the wave propagation in a half space loaded by an impact in a circular surface area, the dynamic stress intensity factor of a penny-shaped crack, the focussing of the von Schmidt wave in an impacted cylinder and the wave focussing in a hemispherically ended cylinder. For reasons of validation, comparisons are made with results gained by other authors. Some new improvements and explanations of the stress increases at cracks and in focussing zones could be attained.
1. Introduction Axis symmetrical stress wave propagation in elastic solids, which may be caused by an impact or another impulsive loading, is governed by a system of hyperbolic partial differential equations with a source term. Solutions of this system are of great practical interest for various initial and boundary conditions with cylindrical symmetry. Some analytical solutions were obtained, e.g. by Laturelle [ 1,2] for a half space using Laplace and Hankel transformations, and by Miklowitz [ 31 for a rod using the approximate Mindlin-Herrmann theory. These show that the source term complicates the analytical solving of the system in a way that does not occur in the plane problem. Therefore a numerical method giving a good approximate solution will be very useful for practical applications. Hirose and Achenbach [ 41 developed a time-domain boundary element method to study the elastic wave interaction in an axisymmetric body. For plane problems without a source term, Lin and Ballmann [5-81 have extended the Godunov-type characteristic-based finite difference methods of gasdynamics for stress waves in elastic-plastic solids and have obtained a great number of results. Problems with cylindrical symmetry can be treated in a similar way. *This paper is dedicated * Corresponding author.
to professor
Alan Jeffrey, University
of Newcastle
0165-2125/95/$09.50 @ 1995 Elsevier Science B.V. All rights reserved SSDIO165-2125(94)00046-8
Upon Tyne, on the occasion
of his 65th birthday
116
X. Lin, J. Balltnann/Wave
Motion 21 (1995) 115-126
Nevertheless, the existing numerical schemes dealing with hyperbolic systems with a source term seem not yet as well-developed as those for the systems without a source term. A widely used method is the time splitting technique which alternately solves a system of conservation laws without any source term and a system of ordinary differential equations modeling the source effect. However, it seems that this technique can produce misleading results, see Westenberger and Ballmann [ 91. For one-dimensional problems, some promising efforts were made by Glimm et al. [ lo], Glaz and Liu [ 111 and Roe [ 121. In this paper, we first propose an explicit finite difference scheme for the numerical integration of hyperbolic PDEs with a source term. Then, this scheme is applied to solve axisymmetric problems of elastic waves, e.g. the wave propagation in a half space due to an impact on a circular area of the boundary, the dynamic stress intensity factor at the tip of a penny-shaped crack and the focussing of stress waves caused by boundary effects. As far as possible the results will be compared with those of other authors in order to show the efficiency of the numerical scheme or to achieve a clearer explanation for the considered problem. 2. The governing equations and the numerical scheme Let r and z be the cylindrical coordinates and t the time. The elastodynamic equations for isotropic, linear elastic solids undergoing torsionless axisymmetric deformations can be written in the following form: au 5
au ar aa, --+x+;, 57ar
aa, 87 =ar+az+~,
1 da, --=PC; at
au ar +a$+&
--=i aa0 PC; at
au au cy$+ff-jy+;’
r
u
7
a24 au -- i aaz =a$+%+“;’ PC; at
u
au -- i a7 = av z+az. PC; at
(1)
where, u and v are the particle velocities in the directions of r and z, gr, a@, gz and r = arz are the stress components, cl and c2 are the longitudinal and transversal wave speeds, respectively. For isotropic materials, LY= I - 2(cz/ct )2 > 0. In order to study the numerical scheme, it is convenient to rewrite ( 1) in matrix form aw T=
af -+$+h. ar
(2)
Introducing Jacobian matrices A, B and C so that
f = Aw,
g = BW and h = Cw, (2) will take the form:
aw at
-=A$+B%+Cw.
(3)
Suppose the solution domain in the (r, z ) plane to be divided into rectangular cells and let w$ = w( nA t, (i i)Ar, (j - 4)A.r) be the values in the cell center (i,j) at the time level P. To obtain a numerical scheme, w is expanded into a Taylor series with respect to time up to the second order, w(t+At,r,z)
=w(t,r,z)
+Atf$+-2-x.
At2 a2w
If the time derivatives in (4) are eliminated by (2) or (3)) one obtains w(t+At,r,z)=w(t,r,z)+AtA;
+AtB&[w+~($+~+h)]
w+~ Z+$+h [ At(af
>I +AtC[w+$($+$+h)].
(5)
117
X. Lin, J. Ballmann/Wave Motion 21 (1995) 115-126
According to (5) and to the fact that quadratic cells in space are adequate bodies, a two-step explicit scheme is obtained as follows: 1 wIf+v2 r+1/2,j+1/2= -4
+ (wi)+ wY+t+l,j +a ( + wEj+l
gY+l.j+l
ti,j+l
in case of linear elastic isotropic
A
+ wY+l,j+l
> ( +
&+1,j - g:j) +
2
fY+l,j+l
-
y(h;+
fEj+l
hy+]qj
+
+
hrj+,
-
fY+l,j
+
ftj)
hy+,,j+l)’
(6) f"+w r+,,2,j_,,2)/2 and analogously for g, and A = At/Ar = At/Az. The scheme (6) represents an extension of ii was’s method [ 131 for the system of hyperbolic PDEs with a source term. In the numerical computation, the domain of the elastic body will be divided into a limited number of cells. The cell centers always lie inside. The boundary conditions are given at the grid points on the surface. Suppose the grid point (i + i, J + i) to be located at an upper boundary z = constant, where the stress components r = T and uz = Q are prescribed. Then the other components, i.e. U, o, ur, ~6, can be calculated using T and Q as well as the function values in the cells (i, J) and (i + 1, J). For example, the schemes for the o and (+r are follows:
wheref:+‘;/lj = (f
:‘:/&+I
+
2
(7) The conditions at other boundaries can be dealt with in a similar way, especially for the boundary r = 0, where u = 0, T = 0 are prescribed, and rl = Ar/2 is used to calculate the other components. In the second step of (6) one has r;_l/2 = 0 in h for the first column cells with the index i = 1. In this case h is treated in such a way, that the singularity l/r is avoided, e.g. approximating the term u/r as At!
r
= At-
U-O
(8)
r-0
The scheme admits everywhere
3. A half-space
the Courant-Friedrichs-Lewy
(CFL)
number cl A = 1.
problem
The first numerical example refers to a half space which is subjected to an impact over a circular area. This problem was considered by Laturelle [ 1,2] who derived the exact solution for the stress cz along the z-axis. We compare our results with those of Laturelle in order to examine the correctness of our numerical scheme.
X. Lin, J. Bahunn/Wave
118
Motion 21 (1995) 115-126
“O 1
-
Laturelle (1991)
t = 1.5
----
This work
r=O
0
1
Fig. 1. Elastic half-space subjected to an impact over a circular area. fig. 2. Comparison of the stress (TVgiven by the present numerical solutionand the analytical solution of Laturelle along the z-axis for the time level t = 1.5.
The material constants of the elastic body are taken as p = 1, ct = 1, c2 = l/d (which corresponds to v = 0.25 for Possion’s ratio). The radius of the loaded arca, see Fig. 1, is set to a = 1, which is divided into 400 cells, i.e., Ar = hr. = At = l/400. The loading is a shock with the amplitude q = 1 and the acting duration of 0.5. The flux components (T, on the boundary z = 0 are set to -q(t) for the grid points located in 0 5 r < 1, while (T, = -q(t) /2 is applied to the grid point on I = 1. The numerical results for (+z along the axis of symmetry for the time levels f = 1.5 and 2 are plotted in Fig. 2 and Fig. 3 in comparison with the exact solution of Laturelle. One can see that the two solutions agree well in the most smooth regions. A strong defect is found near z = 0.5692 at r = 2, because this is a meeting point of all shear wave fronts, in which the stress is undetermined, see Fig. 1 of Ref. [2]. The shock waves appearing in z = 1.108 at t = 1.5, or z = 1.732 and z = 1.115 at t = 2 represent the meeting points of the longitudinal wave fronts. These shocks are modeled over a narrow region, since they move along the axis with a speed c < ct. However, the longitudinal loading and unloading shock fronts with speed ct were calculated correctly. 4. A penny-shaped
crack subjected
to a tension shock
The dynamic stress intensity factor at the tip of a penny-shaped crack was calculated by Chen and Sih [ 141 using integral transformation techniques, and by SMdek and Sl&lek [ 151 solving the boundary integro-differential equations, and by Hirose and Achenbach [4] using a time-domain boundary element method. Zhang and Gross [ 161 have also applied a time-domain boundary element method to this problem and obtained the same results as Ref. [ 41. It is noteworthy that the results from the three methods show remarkable differences. The considered problem is shown in Fig. 4a, where two plane longitudinal tension shock waves travel towards the penny-shaped crack from opposite sides. Because of the symmetry, it is sufficient to solve the equations only in the first quadrant, see Fig. 4b. In our calculation the number of cells along the edge of the crack has been chosen as Zc = 100 and those of the total computational domain are It x Jt = 910 x 810. In normal direction on the line z = 0 the boundary condition is specified by gz = 0 on the crack edge and by v = 0 on the outer part, while in the tangential direction the shear stress is r = 0 along the whole line. The boundary conditions are satisfied by using these values for the flux computation in the boundary grid points. The conditions on the outer boundary of the computational domain can be dealt with in a simple way. Taking the boundary i = Ii as an example, the function values in the cell (It, j) are set to those of the cell (It - 1, j) in every time step. By this approach a plane wave can be kept until the disturbance from the crack tip comes there. Therefore, it will need 2 x 810 time steps for the disturbance wave to take a round trip from the crack tip to boundary and then back to the crack tip, which is sufficient to avoid any influence of the outer boundary on the crack tip during
119
X. Lin, J. Ballmann/Wave Motion 21 (1995) 115-126
I.”
-
Laturelle
t=2
(1991)
,4-z
I_
0,
Cl
-
y
L-l m
J
4-‘-c-‘T-’
Q 7-v
f.
&lo
2
=
‘-‘F
0
3
Fig. 3. Comparison of the stress gL given by the present numerical solution and the analytical solution of Laturelle along the z-axis for t = 2.
the time level
Fig. 4. Sketches for the crack problem. (a) Physical problem, (b) Zoning for the calculation.
the time interval 0 5 czt/a 5 9. The material constants are chosen as cl = 1, c2 = l/fi. The initial conditions were two homogeneous mirror symmetrical states with respect to the line z = 0 and characterized by following values
go=uo= 0 L
0: = a; = -
uo= 70 =
41
v
4”=0.5.
0
l-v4’
0,
(9)
The dynamic stress intensity factor can be obtained in every time step by the following relations (see Fig. 5) (10)
ao z=-zE
4
-312
34 cos -1
2
ao
4 _ -_-E-312
az-
&
sin Z9
2’
(11)
where the values of a@/& and a@/& can be taken in the points C, D and E in order to calculate KI. The results for the normalized dynamic stress intensity factor are shown in Fig. 6 and compared with those of the above-mentioned authors. Over wide time ranges, our results agree well with those of Hirose and Achenbach [4]. But a considerable difference appears for the peak value and its phase position. Any of the other methods shows the peak position clearly at a normalized time less than @/a = 2, while our method gives the peak position at czt/a = 6, which corresponds exactly with the time that a Rayleigh wave needs to propagate along the crack edge from the crack tip to the center point r = 0 and then back to the crack tip. For the discussed material v = 0.25, Q/CR = 1.08766. The result of Chen and Sih [ 141 is based on a material with v = 0.29, Q/CR. = 1.0801. Therefore, there is also 6 > 2. According to the basic analyses for the plane crack by Thau-Lu [ 171, Kim [ 181 and Lin-Ballmann [ 51, the effect of the Rayleigh wave reflected at the crack tip can be assumed as responsible for the steep decrease of the SIF from its first maximum. This explanation supports our present result.
5. The initial wave pattern in an impacted circular md
Impacts on cylindrical rods have been analyzed in many experiments. One typical example is the split Hopkinson bar. In practice, it is mostly assumed that the stress wave propagation inside the rod is a homogeneous plane wave. The wave dispersion due to two-dimensional effects is then neglected. In order to analyse this
X. Lin. J. BallmandWave
120
: -r
Motion 21 (1995) 115-126
lJ.y-z i-2-a
c
--._____._._.-.-.-
KI
K;t
0.5
1 A crack face
’
O.Oi~,,,~,,,,,,,,,,,, 0 1 2
0
This work Hirose and Achenbach (1989) Slddek and Slidek (1986) Chen and Sih (1977) 3
4
5
6
7
8
c2tla Fig. 5. A sketch for the treatment of the crack tip
neighbouring cells.
Fig. 6. Normalized Mode I dynamic stress intensity factor at the tip of a penny-shaped crack under Heaviside tension wave.
phenomenon, the integral transformation technique was applied to the problem by Pochhammer and Chree (see Ref. [ 191) early in last century. They obtained a frequency equation for the approximated calculation of the wave speeds up to the second term of the series. Therefore, the wave dispersion has been taken into account in a homogenized manner. For a practical application with given boundary and initial conditions, Mindlin and Herrmann have proposed an approximate one-dimensional theory, which has been used later by Miklowitz [ 31 to obtain some solutions. The wave propagation in an impacted cylindrical rod can be calculated quite correctly by the numerical method proposed in this paper. In our numerical calculation, the rod has a semi-infinite length and occupies the spatial domain 0 < r 5 1, 0 5 z 5 00. The assumed material constants are ct = 1, c2 = l/d, see Fig. 7. The rod which was originally at rest and stress free is subjected to a tension shock on the limiting surface z = 0 at the time t = 0. The loading condition can be expressed by a Heaviside function gL (t, r, 0) = H(t), which implies that it is uniformly distributed over the surface after the time t > 0. The radius a = 1 of the rod is divided into 50 cells and the number of cells along the z-axis is taken according to the upper time level of the computation with CZ=l. The calculated distributions of the stress components ur and crZ in the rod at time t = 3 are shown in the first part of Fig. 8. The plotted values represent the data in the cell centers. At t = 3 the plane wave front due to the impact has just reached z = 3. The stress u, behind this wave front does not vanish due to the transient plane strain state induced by the impact. Therefore, a von Schmidt wave is produced in the crossing point of the shock wave front z - ct t = 0 and the side boundary r = 1 in order to satisfy the boundary condition g’r = 0 on the surface r = 1. This von Schmidt wave causes also a negative velocity u in the r-direction, which can change the stress state again when the wave propagates towards the axis of the rod, since we have u = 0 along the axis because of symmetry. As a result, a high stress intensity will appear on the axis due to the focussing of the von Schmidt wave. The stress (+Bhas also a focussing peak at the same time and its distribution looks similarly to that of a,; therefore it is not drawn in the picture. In the focussing region, ur and a@ are negative, but uZ is positive. Then the von Mises yield stress K, defined by the second invariant of the deviatoric stress tensor
K2=+.-uz)2+(us -
uz J2+ (ur - ud2] + 72,
(12)
will also have a focussing peak on the axis, which is shown in Fig. 8, too. Theoretically, the focussing point on the axis is a singular point of the elastic solution, where the stresses become infinite. Indeed, the values of the focussing peak become higher if a finer grid is used in the numerical
X. Lin, J. Ballmann/Wave Motion 21 (1995) 115-126
Fig. 7. A circular rod subjected to a tension shock.
max
=
2.06
Fig. 8. The stress distributions in an impacted rod at time t = 3.
121
122
X. Lin, .I. Ballmann/Wave Motion 21 (1995) 115-126
max Hg. 9. The time history of the stress uz (0,
= 2.15
t) along the axis.
2,
2.5 t/a
L 2.0 -
-.-.ii
z/a = 4 (this work)
_____, z/a = 2 (Miklowitz,
1957)
z/a = 4 (Miklowitz,
1957)
1.5 -
I
0.0 0
1 2
I
= 2 (this work)
I 4
I
I 5
I
I 8
I 10
(cd - z)/a Fig. 10. The time history of stress uL at discrete axial points, compared with the approximate results of Miklowitz.
can occur in a region near the axis. The focussing peak appears for the first time when the von Schmidt wave has reached the axis of the rod. Figure 9 shows the time history of the distribution of gZ in the first row of grids along the axis. It can be seen that the maximum of the peak appears at f = 3.54. Its value is cf, = 2.15 at about z = 1.94. Then the peak decreases. At about t = 7, a second focussing peak begins to grow. The wave train formed by the leading tension wave and the von Schmidt wave focussing at the axis of the rod and being reflected at the boundary forming a new tension front maintains always its typical geometry. Figure 10 presents the time histories of the stress CT~ in the fixed axial points z/a = 2 and 4 in order to compare the present numerical results with the analytical solutions found by Miklowitz [ 31 using the approximate Mindlin-Herrmann theory. The input data were chosen as cl = 1 and c2 = 0.5821 in accordance with Possion’s ratio v = 0.24375 used by Miklowitz; the geometric parameters were chosen the same as above. Obviously, there are large differences between the two results, especially with respect to the focussing peak behaviour. The example clarifies the defects of the one-dimensional Mindlin-Hemnann theory in the leading part of the wave over a distance of about 10 radii behind the head of the wave.
X. Lin, J. Balhann/Wave Motion 21 (1995) 135-126
123
Fig. Il. A sketch of a cylindrical bar with a spherical free end.
grid 2 0
r a=1
Fig.
L
12. A sketch of the combination of two grids. The data are exchanged in every time step from grid 1 to grid 2 on those cells near the boundary denoted by X, and from grid 2 to grid 1 on those cells near the boundary denoted by o.
6. The focus of a stress wave in a hemisphere In this last example we consider wave focussing in a cylindrical bar with a spherical free end, see Fig. 11. When a shock wave is induced into the bar at the plane left end, the stress wave will propagate into the cylinder as described in the last section. Since the surface of the hemisphere at the right end is free, the longitudinal stress wave will be first reflected with wave splitting and phase changes and then focused on the axis at a point with the distance of half a radius from the end. A problem similar to that but under plane stress condition has been considered by Ballmann et al. [20] and Niethammer et al. [21]. This calculation is a test for the axisymmetric problem. The hemispherical surface represents a curved boundary with a boundary condition that the rectangular cells used in previous sections fail to satisfy correctly. Of course one could use triangles with edges to fit the boundary geometry, but the discretization of the governing equations on an irregular grid will cause too much numerical dispersion and dissipation. Therefore we apply an overlapping grid technique [ 211 to compute the problem. Suppose the system is still axisymmetric on the z-axis. Then we choose one grid still rectangular in r- and z-direction, while the second grid has to fit the surface, i.e. in the present case we use spherical coordinates (R, p). The arrangement of the two grids is shown in Fig. 12. The governing equations of elastodynamics in spherical coordinates can be found in Ref. [ 221. After some changes of symbols and coordinate directions in accordance to Fig. 12, we may rewrite them for the torsionless case with cylindrical symmetry as follows,
124
X. Lin, J. Balhann/Wave Motion 21 (1995) 115-126
case with cylindrical symmetry as follows, 801 =!%+L!%+ GR asp 1 au1 --=PC; at
2ai -u2-us-ri2tanSp
au1 aR+Zau2+ R a9
au1 aav2 --1 au3 PC; at =“=+w+
&I2 =!Z+L!%+
R
’
2avi - (Yvztan (p R
’
R
3712 - (~2 - (~3) tmqo R
R a5o
I 1av2+u+4vl-~v2t~9 --1 aa2 _avl & at aR R R ap
-v2mp
(~+cu)v~
5-
’
1 aTl2 au2 i avl --= zi+ii&zT PC; at
v2
(13)
where vi and vz are the velocity components in R- and qP-direction, respectively, and ui = (TR,02 = a,+,,us = a0 and 7-12= URq are the stress components. It is evident that this system can be integrated numerically by the scheme (6). The radius of the cylindrical rod is chosen as a = 1. For grid 1 we take 50 cells in r-direction and 100 cells in z-direction (including 50 cells for the hemisphere). The radius of the hemisphere is equal to a. The number of cells of grid 2 in R- and p-direction are set to 11 and 59, respectively, which gives an inner radius of grid 2 of about 0.746. The smallest cell’s length AR of grid 2 is about l/50 in order to obtain approximately the same CFL number as in grid 1. The computations are carried out in the two grids separately, but the data are exchanged in every time step by linear interpolation from grid 1 to grid 2 in the cells near the boundary denoted by x and from grid 2 to grid 1 in the cells near the boundary denoted by o. The material data are The boundary condition is given by uZ = -H(t) at the rod’s plane end chosen such that ct = 1 and cz = l/fi. surface z, = 0. The computation has been carried out to analyze the stress focussing in the hemisphere, which is to be expected by ray arguments at R = a/2 on the axis of the hemisphere. The highest focussing intensity is found at time step N = 122, i.e., tit/a = 2.44. The momentaneous distributions of the stresses ue, u, and u, at this time are presented in Fig. 13. The focussing results shown in Fig. 13 are due to the reflection of a plane shock wave front from the spherical boundary. There exists also a negative peak with the intensity uZ = -1.76 in this picture, which denotes the focus of the von Schmidt wave discussed in the last section. This peak is still moving in z-direction at this time. We have continued the computation until the negative peak was reflected from the hemisphere’s surface in order to see whether another focussing would occur. The result showed indeed another occurrence of focussing, but the stress increase is a lot weaker as presented in Fig. 13 for a, and uL. That is, because the focussed part of the von Schmidt wave is much smaller than the one of the plane shock wave. From Fig. 13 it can be seen that all stress components ue, a, and uZ are positive in the hemisphere’s focussing area. The shear stress T vanishes along the axis because of symmetry. Therefore, according to ( 12), the von Mises yield stress will not show a focussing peak in the discussed case.
7. Conclusions In this paper a finite difference scheme was proposed for the numerical integration of systems of hyperbolic PDEs with a source term, which is typical for cylindrical or spherical geometries. The elastic wave propagation in two-dimensional axisymmetric solids can be handled with the proposed numerical scheme. Several results have been obtained, which show the efficiency of the method compared to analytical methods. For the pennyshaped crack and for cylindrical and spherical geometries some new explanations of stress increases and focussing effects were presented.
X. Lin, J. BallmmdWave
max
=
2.30
max
=
2.37
max
=
2.11
min
= -1.76
Motion 21 (1995) 115-126
Fig. 13. Stress focussing of a plane wave in the hemisphericalend of a cylinder.
125
X. Lin, J. Balbnann/Wave Motion 21 (1995) 115-126
126
Acknowledgements
The research work presented in this paper was supported by the Deutsche Forschungsgemeinschaft under Grant No. Ba 661/ 12-1, which is gratefully acknowledged. References [ 1I PG. Latutelle, “The stresses produced in an elastic half-space by a normal step loading over a circular area. Analytical and numerical results”, Wave Motion 12, 107-127 (1990). 121 EG. Latutelle, ‘The stresses produced in an elastic half-space by a pressure pulse applied uniformly over a circular area: role of the pulse duration”, Wave Motion 14, l-9 (1991). [ 31 J. Miklowitx, “The propagation of compressional waves in a dispersive elastic rod, part 1 - results from the theory”, J. Appl. Mech. 24.231-239 (1957). [4] S. Hirose and J.D. Achenbach, ‘Time-domain boundary element analysis of elastic wave interaction with a crack”, Int. J. Numer. Methods Eng. 28.629-644 ( 1989). [ 51 X. Lin and J. Ballmann, “Reconsideration of Chen’s problem by finite difference method”, Eng. Fracture Mech. 44.735-739 (1993). [ 61 X. Lin and J. Ballmann, “Numerical method for elastic-plastic waves in cracked solids, part 1: anti-plane shear problem”, Arch. Appl. Mech. 63,261-282 (1993). [7] X. Lin and J. Ballmann, “Numerical method for elastic-plastic waves in cracked solids, part 2: plane strain problem”, Arch. Appl. Mech. 63,283-295 ( 1993). [ 81 X. Lin and J. Ballmann, “Improved bicharacteristic schemes for two-dimensional elastodynamic equations”, to appear in Quart. Appl. Math. (1995). [9] H. Westenberger and J. BaJhnann, “The homogeneous homentropic compression or expansion - a test case for analyzing Sod’s operator-splitting “, in J. Ballmann and R. Jeltsch (eds.): Nonlinear Hyperbolic Equations - Theory, Computation Methods, and Applications, 678-687, Vieweg, Braunschweig/Wiesbaden ( 1988). [lo] J. Glimm, G. Marshall and B. Plohr, “A general&d Riemann problem for quasi-one-dimensiotud gas flows”, Adv. Appl. Math. 5,
l-30 (1984).
[ 111 H.M. Glaz and T.-P Liu, “The asymptotic analysis of wave interactions and numerical calculations of transonic nozzle flow”, Adv. Appl. Math. 5, 11 l-146 (1984).
[ 121 P.L. Roe, “Upwind differencing schemes for hyperbolic conservation laws with source terms”, in: C. Carasso, P-A. Raviart and D. Serre (Eds.), Nonlinear Hyperbolic Problems, 41-51, Springer, Berlin (1986). [ 131 B. Eilon, D. Gottlieb and G. Zwas, “Numerical stabilizers and computing time for second-order accurate schemes”, J. Comput. Phys. 9, 387-397 (1972). [ 141 E.F? Chen and G.C. Sih, “Transient response of cracks to impact loads”, in G.C. Sih (cd.): Mechanics of Fracture 4 - Elustodynamic Crack Problems, l-58, Noordhoff, Leyden (1977). [ 151 J. SlPdek and V. Sladek, “Dynamic stress intensity factors studied by boundary integro-differential equations”, ht. J. Numer. Meth. Eng. 23.919-928 (1986). [ 161 Ch. Zhang and D. Gross, “Transient elastodynamic analysis of a penny-shaped crack”, Eng. Fracture Mech. 46.641-654 (1993). [ 171 S.A. Thau and T.H. Lu, “Transient stress intensity factors for a finite crack in an elastic solid caused by a dilatational wave”, ht. J.
Solids Strucmres 7, 731-750 (1971). [ 181 K.S. Kim, “Dynamic propagation of a finite crack”, ht. J. Solids Structure 15,685-699 (1979). [ 191 A.E.H. Love, A Treatise on the Mathematical Theory of Elasticity, 4th ed., 287-292, Dover, New York (1944). [ 201 J. Ballmann, H.J. Raatschen and M. Staat, “High sttess intensities in focussing zones of waves”, in: P Ladeveze (Ed.) : Local Effects in the Analysis ofStructures, 235-252, Elsevier Science Publishers, Amsterdam (1985). [21] R.J. Niethammer, K.-S. Kim and J. Ballmann, “Numerical simulation of shock waves in linear-elastic plates with curvilinear boundaries and material interfaces”, submitted to Int. J. Impact Eng.. [ 221 K.P. Graff, Wave Motion in Elastic Solids, 601-602, Clarendon Press, Oxford ( 1975),