A numerical method for axisymmetric wave propagation problem of anisotropic solids

A numerical method for axisymmetric wave propagation problem of anisotropic solids

Computer methods in applied mechanics and englnraring Comput. ELSEVIER Methods Appl. Mech. Engrg. 145 (1997) 109-116 A numerical method for ax...

544KB Sizes 0 Downloads 94 Views

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