Computer methods in applied mechanics and englnraring Comput.
ELSEVIER
Methods
Appl.
Mech.
Engrg.
145 (1997)
109-116
A numerical method for axisymmetric wave propagation problem of anisotropic solids Kaishin Liu*, Xin Li, Xiannian Research Instifute of Engineering
Mechanics,
Received
3 April
Sun
Dalian University of Technology, 1996; revised
17 September
Dalian 116024, P.R. China
1996
Abstract A numerical algorithm based on the method of characteristics is presented for calculating axisymmetric wave propagation problems of anisotropic solids. The algorithm is applied to study the dynamic behaviour of a transversely isotropic tube due to impact loading. The stability of the present algorithm is investigated. The limitation in the treatment of boundaries and the advantages of present method are discussed.
1. Introduction
The general axisymmetric problem of stress wave propagation in anisotropic elastic media is a mixed initial and boundary value problem. From the axisymmetric properties, the material property of the media can be described by a transversely isotropic solid. The governing equations form a system of hyperbolic partial differential equation in three independent variables. Because of the well-known difficulties associated with obtaining exact solution to such a problem, it is important to develop an accurate and efficient numerical algorithm for achieving its approximate solutions. Among the numerical procedures, the method of characteristics [1,2] seems to be the most suitable for analysing the stress wave propagation problems. This method has been successfully applied to solve twodimensional stress wave propagation problems in isotropic solids, such as the plane problems and the axisymmetric problems in elastodynamics [3-61, and dynamic elastoplasticity [7]. This method also can be applied to simulate elastic wave propagation in the anisotropic media [8]. As one-dimensional problems, Greif and Chou [9] investigated the propagation of radially symmetric stress waves in anisotropic nonhomogeneous elastic media, and Mengi and McNiven [lo] calculated the transient response of a semi-infinite transversely isotropic rod due to an impact loading by means of an approximate theory, by using the method of characteristics. In this investigation, a numerical algorithm for obtaining the numerical solution for typical symmetric stress wave propagation problems in transversely isotropic solids is presented. The numerical technique is basically equivalent to the method of finite difference along bicharacteristics. The algorithm is an explicit second-order accurate procedure. The stability of the method is evaluated by both results on the von Neumann necessary condition and examining the practical numerical calculation. Moreover, the dynamic response of a transversely isotropic tube subjected to an impact loading is calculated.
* Corresponding
author.
0045.7825/97/$17.00 f$J 1997 Elsevier PII SOO45-7825(96)01204-2
Science
S.A.
All rights
reserved
2. Analysis For an axisymmetric cylindrical co-ordinates written as
wave propagation problem of anisotropic (1.0, z). Because of the symmetric property,
solid, it is convenient to take the equations of motion can be
(1)
where p denotes the velocity components, is homogeneous and material are given as
V, and ZJ~the particle mass density, Us. w#, (T_ and Us; the stress components. and t the time. It is assumed that the deformations are very small and the material elastic. As a result, the constitutive equations for the transversely isotropic follows:
(2)
where c,,, c,:, c,~, ci3 and cl4 are elastic coefficients learned from [ 1 I]. The governing equations A’lJ, f AWr
whose
(1) and (2) can be expressed
+ A’U,,
relationship
in the following
with elastic symmetric
constants matrix
+ BU = 0
can be form: (3)
where
in which U,,, U,,, U,; are the partial derivatives of U with respect to t, r and Z, respectively. and B are symmetric matrices with the following non-zero components:
A’,,=Ai2=p,
2c3,
A$= CJ&,,
A;, =
BI=$
CII + Cl? C&I,
+ c,z) - 24,
B2h=B3,=$.
’
+ c,z> - 2& A& =
A:, = ’
A’, A’, A’,
2 CII -cIz
-2c,, c&,,
+ c,T!f - 24,
y
B,, =+.
Eq. (3) represents hyperbolic system of first-order partial differential equations with constant coefficients. The solution of Eq. (3) can be obtained by means of the method of characteristics. By analysing the characteristic properties of the governing equation (3), we can obtain the desired differential relations along bicharacteristics which are the generators of the characteristic cone.
K. Liu et al. I Comput.
Methods Appl. Mech. Engrg.
111
145 (1997) 109-116
It is convenient to use the bicharacteristics on the surfaces r = 0 and z = 0. So the differential relations along the eight selected bicharacteristics may be obtained as
da,
dvr
dt-w,dt=
da;
dv,
-dt-Pcldt=
(ur ( ___
-cl
utl
-
--Cl
au,,, a2 >
+
r
+
u,-u~ r
_TI
au
av
(
3+,,,-g
clzr
+
az >
(
>
av
c
**
V’fc,,$ > r
da,, dt
do;* dt
dun
avz
0; - “6 r
dt
dcr,,
&
dvz
do-2 7-Pc3dt=-C3
duz
-dt-Pc3dt=-C3
ru -u
dv ---r-=-
--PC2
dt
dvz
+ c44 7
-
r
(aa,, a,, (& @f-z> ( ar
+r
ar
)
+r
+
(
-
Cl3
Iav
c44 ar ayl+c,3: ar
av
c13*+c*37
)
v, )
(4)
wnere c, = vc,, Ip, c2 = vc441p and c3 = Vc33 lp are normal wave speeds corresponding to coordinate axis directions, respectively. Integrating Eqs. (4) along the corresponding bicharacteristics and Eq. (3) along the time axis, and performing proper linear combinations of those equations, we can successively obtain the following expressions for six unknown increments 6v,, au,,, 6v,, Zb,, 6u8 and 6~~: tiv, =
2r2k 4pr2 + k2(c,, - c12) +
2u,,,, + 4 c;k2
+ k”..,‘.’
“”
2u,,, + ; c;k2 (v)
,,, + c,zk(:)
., + cnkV,.,z + c,lkv,,,,
u -u z ,z + v,,,,) + 2+
,2;,,+;$;j,,+u,z.,,+;cllk($),,,+u,,,,)
~,,=~~2yil+~(~,.,~+u,‘z7~~‘z+~,~,~~)+2~~.,+~(u~.,~+(~)~,+u,~~,,)}+O(k3) k
6% = 2p
1 a,, ZI ,,, + c&V,,,, + c44kv~:,,,+ 2u,., + 7 c:k2 1.
-2T,l+O(k’)
where time.
- k, x,,, y,, ). k is the time step.
Sf =f(f,,, s,,. J’,,) -f(/,,
(t,,, s,,, _v,,) is an arbitrary
point
in space-
The quantities on the right-hand side of Eqs. (5) are the values at the point (t,, - k, ra, z,,), which arc all known. and the partial differences f:,-, .f:_. j’,-,. f’_._. ./‘,-,. at point (r, - k. r,,. z,,) which can be approximated by central differences [3]. When the point (t,, - k. Y,,. z,, ) is located on the boundary. which includes some surfaces normal to the plane z (or r) = constant, some bicharacteristics on the characteristic surface intersect the plane 2 (or r) = constant at points outside the region of interest. Therefore, the relations along those bicharacteristics cannot be used. Those unused relations, however, are substituted with admissible boundary conditions for construction of the numerical solution at the point located on the boundary.
3. Numerical
examples
By using the numerical technique developed in the present paper, we calculate several numerical examples simulating stress wave propagation for a semi-infinite tube made of E-type glass-epoxy composite. The material constants of E-type glass are: E = 72.5 GPa, G = 30.4 GPa, Y = 0.20, p = 2.55 X 10” kg/m’. The material constants of epoxy are: E = 3.11 GPa, G = 1.17 GPa, v = 0.35. p = 1.214 x lo3 kg/m’. The relationships between the engineering elastic constants and the properties of fiber and matrix constituents are shown in [I 11. Four cases with different ratio c,/c, = 1.0, 1.5, 1.X and 2.0, which are, respectively. corresponding to “; = 0, 0.1010. 0.2185 and 0.4350. have been investigated, where V, denotes the fiber volume fraction [ 111. The initial and boundary conditions are given by a, = (T#= cr: = CT;= -P(t),
‘J,:
=
0.
v,
=
u,
=
0
*
fort = 0 . for z = 0 .
‘7,; = 0 .
for r = R, or r = R, _
1 g = a,, = 0 , where
P(f) =
12 i 1 -co,,+ { PC,
,)
u,,
for j S t, for t > t,
in which I?, and R, are the inner and outer radius of the tube, respectively, the rise time tr = 10 p.s. the pressure P,, = 150 MPa. The spatial mesh size is taken equal to h = 0.00162 m, the thickness of the tube is taken as d = 2Oh, and the time step size is taken as k = 0.5 h/c,. Fig. 1 shows the distributions of the longitudinal stress CT:for different materials in the middle surface of the tube r = (R, + I??)/2 at time I = 100 ~J-s.From Fig. 1. it can be seen that the difference of material properties can only change the field of the longitudinal stress wave propagation (because the wave the tendency of the distributions of the longitudinal stress is not speeds are different). However, changed very much. Fig. 2 indicates the distributions of the longitudinal stress CT:.in the middle surfaces of the tube at t = 100 t_~sfor different inner radii. It is clear that, for the distributtons of the longitudinal stress, the difference between the thick tube and the middle-thick one is very large, but the difference between the middle-thick tube (d/R, = 0.5) and the thin one (d/R? = 0.0s) is very small. Figs. 3 and 4 show the distributions of the longitudinal stress ~7~near the impact end of the tubes. As we expected,
K. Liu et al.
113
i Comput. Methods Appl. Mech. Engrg. 145 (1997) 109-116
250
E
d/F&=0.6
-_cs/cr --823/c*
t=lCQlrs
--c+sq
= 1 .o = 1.5
= 3.8
50
50
o-
0
0
0.1
0.2
0.3
0.4
0.5
0
0.1
0.2
2 (m) Fig. 1. Longitudinal Fig. 2. Longitudinal
stresses stresses
0.3
0.4
0.5
z (m) at the middle at the middle
surface surface
at time t = 100 ys for the different at time I = 100 p.s for the different
material. inner
d/R2=0.20 c3/c,=1.5
Fig. 3. Distribution
of longitudinal
stress
of the thick tube for c,ic,
= 1.5.
Fig. 4. Distribution
of longitudinal
stress of the thin tube for c,ic,
= 1.5.
radii
t=ioow
the distribution of crz along the r direction varies violently in the thick tube (d/R, = 0.95), which seems a two-dimensional dynamic behaviour; while gently in the thin one (d/R, = 0.20), which is almost like one-dimensional dynamic response. When taking E(1 - u) Cl1 =c33 = (1 + V)(l - 2V) ’
C,? = c,3 = (1 + V;(: - 2V) ’
E c41 = 2(1 + V) ’
the problem is reduced to axisymmetric wave propagation problem of an isotropic tube. Fig. 5 shows the comparison between the present numerical result and the Reeker’s numerical result [5]. The material and shape of the tube and the impact condition are all the same as those in reference [5]. As shown in Fig. 5, almost the same results are obtained.
K. Liu et (11. I Compur. Methods Appl.
II4
Fig. 5. Comparison present
solution;
4. Stability
between
points.
present
Reeker’s
numerical
solution
of the difference
solution
Mech. Engrg. 145 (1997) 109-116
and Reeker’s
solution
for isotropic
tube with
IJ = 0.333.
Solid curve.
[5].
equations
As the inherent complexities to evaluate the stability and convergence of the numerical solution for mixed initial and boundary value problems of hyperbolic systems of partial differential equations, we will not make a theoretical investigation of the stability and convergence of the proposed numerical procedure. However, we can get some useful results by use of the von Neumann necessary condition. With the help of these results, we can obtain the approximate necessary convergent conditions for the numerical analysis. It is easily verified that Eqs. (5) are equivalent to the Lax-Wendroff difference equations [12] for Eq. (3). The stability conditions for the method of characteristics can be obtained by analysing the stability of the corresponding Lax-Wendroff difference equations [2,3,5]. The Lax-Wendroff difference equations for Eq. (3) are based on a Taylor’s expansion of the solution about the point (nk, ro, za). The Taylor’s expansion about the point (nk. ro, z,) is 7
A'U""
= A’U”
+ kA’U;‘,
+;
A/U”,, + 0(k3)
i6)
where I/” = U(nk, ro, z,,) ,
U”’
’ = U(nk + k, rn. z,,) .
By applying Fourier analysis, the amplification
matrix can be derived as
S(e,,&)=R+iJ
(7)
where R=I-k~+~~~+(~)-[(~‘)‘(cos~~-l)+(A’)’(cos~~-1)~ - $($
‘(A ‘A ’ + A ‘A ‘) sin 0, sin O2 ,
J-(X),-r A A’ = (A’)
3
sin 0, + A ’ sin 02] + i
‘A’ ,
‘q: = (A’)-‘A’
‘[ (A’l? (
+ &‘)h
i? = (A’)-?3
sin 8, + (A ‘B + BA ‘)/I sin O,] ,
K. Liu et al. I Comput.
.
m-
. ’
0
Methods Appl. Mech. Engrg.
d/&=0.5 ca/c1=1.5
-csk/h --c&/h --cak/h
2 /2
r=(Rl+FW
z-(RrRl)
100
50
145 (1997)
109-116
115
-0.8 - 0.92 = 0.98
150
tw4
Fig. 6. Effect of ratio c,k/h
on stability of numerical solutions for longitudinal stress.
where 8,, 8, are arbitrary parameters that correspond to arbitrary wavelengths in the r, z directions and I is the unit matrix. The strongest restriction on k/h that ensures satisfaction of the von Neumann necessary condition is 8, = 13,= 7~, for the case that r is very large (i.e. L?= 0, such as a thin tube), the condition for spectral radius of S less than or equal to unity is
It should be noted that the stability analysis described above does not concern the effect of boundary condition and other factors. Consequently, it is just a rough approximation and the results can only give us a guiding rule to select the appropriate value of c,klh. In fact, the value of c,klh for taking a stable numerical result in the present calculation is less than the value obtained above. For example, for the material of c,/c, = 1.5 (i.e. Vf = O.lOl), the stable condition is that c,k/h would be slightly smaller than 0.946 according to the stability analysis described above. However, the calculation shown in Fig. 6 indicates that, for c,klh = 0.8, the numerical solution is convergent and reasonable; for c,klh = 0.92, the numerical solution is divergent.
5. Conclusions
A numerical algorithm, equivalent to the method of finite difference along bicharacteristics, has been developed which is capable of analysing a variety of axisymmetric stress wave propagation problems in the transversely isotropic media. However, in view of the fact that the treatment for general curved boundaries is difficult by the present technique, the numerical algorithm is weak in such axisymmetric problems. So, the computer program is made for analysing the impact problems in the tube, the cylinder, and its coaxial jointers. On the other hand, the present algorithm can be expressed explicitly and no iterative procedure is needed, therefore, it is suitable for a personal computer to calculate. All calculations for the numerical examples presented in the preceding section were carried out by using a personal computer with 486DX CPU (33 MHz). For example, computation time on this machine was 35 min for the numerical result shown in Fig. 4 (400 time steps).
Acknowledgment The authors of China.
gratefully
acknowledge
the financial
support
of the National
Natural
Science
Foundation
References [II
R. C’uurant and I>. Hilbert.
[2] D.S.
Butler.
Proc.
[3]
R.J.
R.
The
Sec. Lend
Clifton.
Structures (51 W.W.
H.J.
ol Mathcmatlcal of hyperbolic
Phy”cs,
II
(Intcrxience
Publishera.
system\ of partial differential
Inc..
equations
NY,
lY62).
in three independent
method for plant
Elsevitx
A numerical
variahlcs,
(1960) 232-X problem\ 111dynamic elasticity.
0.
Appl.
Math..
25
Y7- I 16.
( lY67)
Raatzchen and M. Staat. High stress tntenGtie\ in focusing zones of waves. Local Effects
(P.. ED.
Reeker.
solution
A 25.5
A differcncc
[1] J. Ballmann.
Methods
numerical
Science quhlishers. solution
Amsterdam.
of arisymmetric
problems
1Y7l
)
in the Analysi\
cl!
235-252.
in Elastodynamics.
Int. J. Numcr.
Method5
Engrg.
3
( lY71)
x-377. [h] K.
Liu.
JSME
S. Tanimura. Int.
171 .I. Bejda. Mcch..
181T
CT.
[ 1I]
[ 121P.D. Appl.
of an elasllc arcular
\trcbs wave\ 111Lin clastlc, vlscopl;istic
1.ing. Characteristic
tube due to longitudinal
~mpxr.
form\ of differcntlal
cyuations
matcl-ial.
for wave pi-opagation II) nonlinear
Proc.
12th Int.
media. ASME
(‘ongr.
AppI
J. Appl.
Mech
743%73X. (‘hou. The propagation ot radiall!
and S.C. J. Appl.
Mech. 3X (1971)
and H.D.
McNiven.
\ymmctrlc \trcs\ wa\c\ 111anlsotropic
nonhomogeneous
elastic med~,~.
51-V
Analysis
ol‘ the transient
excitation
of a transverxly
Isotropic
rod. .I. Acoust.
Sot.
Am. SO( I )
25X-265.
Y. Shindo Struct.
of two-dimensional
(1968) 121-131.
[IO] Y. Magi (lY71)
The dynamic hcha\loul
S3.V53Y
LJniverslry
R. Greit ASME
H. Igaki and K. Ka~u.
I 32 (lY8Y)
Propagation
Stanford
48 ( IYXI) ISI
J. Series
and H.
23(l)
(lY87)
Nozaki.
Lax and B. Wcndroff. Math
Impact response
of
;I
transverselp
Isotropic
cylinder
with a penny-shaped
crack.
Int.
J. Solid
l87-IYY.
17 (1064)
Difference
3X1-380.
scheme< wtth hlph-order
ot accuracy for halving hyperbolic
equations.
&mm.
Pure