A wave approach to finite element analysis of solids

MechanicsResearchCommunications,Vol.26, No. 2, pp. 191-196, 1999 Copyright© 1999 ElsevieTScienceLtd Printed in the USA. All fights reserved 0093-6413/99/S-see front matter

Pergamon PII

S01D3..6413(99)00012-9

A WAVE A P P R O A C H TO F I N I T E E L E M E N T A N A L Y S I S OF SOLIDS

B. F. Shorr Central Institute of Aviation Motors, CIAM, Aviamotornaya St.2, 111250, Moscow, Russia

(Received 10 July 1997; accepted for print 5 January 1999)

Introduction

Numerical simulation of dynamical problems in solids on the basis of the finite element method usually involves the representation Uek(X,y,z,t)= uk(t ) q0ek ( x , y , Z), where uek is the displacement of arbitrary point of the element affected by the displacement of the node k and a function'O depends only on the element type [1]. This leads to excellent solutions for stationary free and forced vibrations but does not describe special wave features of nonstationary processes. A way to deal with this obstacle is to present q~ as a time-dependent function associated with the solution of the wave equation for each element [2]. Unfortunately, this is rather complicated in realization, especially for multi-dimensional problems. An alternative approach is proposed in [3,4] where the finite element representation of displacements is referred only to discrete time instants t i . Strong discontinuities o f stresses and velocities are assumed to propagate during time intervals At i = t i - ti_ I from nodes over elements in accordance with the laws of momentum and energy balance. The procedure allows to describe wave propagation without constructing differential equations of motion and using approximation o f displacements or other functions with respect to time. This approach, called as the method of direct mathematical modelling o f wave propagation, is successfully applied to the multiple of one-dimensional problems of mechanics and engineering [4,5,6]. It may be regarded as a variant of the 'wave' finite element method [7]. In this paper the method is generalized for the first time to multi-dimensional problems.

Theoretical background

Assume that the stress-strain state of a finite element solid system is considered at discrete time instants t i . Let us study its transition from instant

ti_ 1 to t i during finite time interval

At i. Here, i=l,2,...,n i is the

number of time steps. By analysing wave propagation it is necessary to differ infinitely close time instants t~-~ and ti~_~. In doing so, we assume that at the instant t~_I ( including the initial instantt 0 ) : a) the velocity

vector

v j0

o f the element j

and

the displacement vectors Ukj 0 of all its nodes n j, where

k j =1 ..... n j, are known; b) the vector sum o f the forces Fkj 0 applied to the element in its nodes is equal to zero and

these forces are

connected with nodes

displacements by means o f standard rigidity ma-

trix of the finite element method; c) the displacements in common nodes o f neighbouring elements are 191

192

B.F. SHORR

continuous, while its velocities may differ and the vector sum of forces applied to a node may be not equal to zero. As a result, strong discontinuities of velocities and forces arise at nodes. Evidently, such conditions always hold at the initial moment if initial stresses are known and all the elements move with the same initial velocity. Using a) and b), all mechanical parameters of an element ( the kinetic energy

K j0, the

elastic potential energy Pj0, the momentum vector Mj0, the stress and strain tensors ) can be determined. Rotary inertia of elements can be neglected when their number is great.

Strains are assumed to be elastic

and small. +

The discontinuities in a solid must immediately decay, therefore, at a time instant ti_ t the velocity of the element j in its node kj becomes vk~ = Vk :~ Vj0 , where v k is a new value of the node velocity, and the forces applied to the element j in its nodes become Fk, * Fk, 0 . The equilibrium condition for the node k which is applicable to the n k neighbouring elements is ll k

Z F,~ +Fk = 0 , jk=l where Fjk = - F k ,

(1)

Jk = I ..... n k , and F k is the vector of outer surface or volume forces applied to the node

+

at the instant ti_ l . If the outer viscous friction is taken into account then the term - rl kvk, where q k is a viscous coefficient, must be added to the force F k . Disturbances of the velocities

Vkj -- Vj0 and of the forces

Fkj - Fkj 0 propagate into element j from its

nodes in all directions interacting with each other in a complicated manner. The following assumptions are introduced in the approach proposed: a) disturbances entirely cover all the elements within the time interval At i = t~- - t ~ 1 ; b) disturbances from a given node arrive other nodes only at the end of a step at the instant t~- ; c) outer forces and b o u n d a r y conditions refer only to nodes and do not vary during the time interval Att;

d) equal

parts of the element mass Amkj = Amj/nj are related to each node k i.

Values of At, vary over time and are determined through calculation. There are two phases of the wave process. During the first half of the interval At~ we do not assume interaction between the disturbances travelling from different nodes. Thus, wave fronts spread f r o m nodes into the u n p e r t u r b e d domain of an element affecting the velocities Vk, and the inner forces Fk, of the parts of the mass Amk,. It follows from the law of m o m e n t u m balance of these parts that 2Amj Vkj = Fkp + ~ ( v k j

-Vjo) .

(2)

Inserting Eq.(2) into Eq.(l) for all the elements surrounding the node k and taking into account the equalities vj~ o = vk, o , Fj~ o = -Fkj 0 , we calculate the new value of the node velocity. It becomes 11k

II k

Vjk0Amjk /n~k +0.5Ati(Fk + S~ Fjk0) Vk

Jk = t

Jk = 1

(3)

IIk

Y~ Amjk / n Jk = l

+0.5qkAt i Jk

Hence we define all the velocities vk, = v k and, using Eq.(2), all the forces Fkj . At the instant

ti+.~ + 0.5At: all the disturbances converge to a conditional 'point of contact' which may be

regarded as an inner element node in which F k = 0. The

new value of velocity at this point becomes

vj and the new values of the forces which arise at this point and are reflected from it become Fkj r . The fronts of reflected waves propagate

into

previously

disturbed

domains characterized by the

parameters

WAVE

APPROACH

TO FINITE ELEMENTS

193

Vk ~ and F k . New disturbances reach the outer (real) nodes during the second half of the interval At i. By using the above procedure inserting velocity vk~ instead of vj0 and forces Fk, instead of Fkj 0 we obtain Fkj ~ = Fk,

2Amj + njAt-------~(vj- Vk, )

(4)

and (after a simple transformation) the new velocity of the element ns

vj = 2 ( ~

Vkj)/n j-vj0

.

(5)

kj=l

It can be shown that the vector sum of the forces Fkj r is equal to zero. The node velocity v k does not vary during the whole interval At i . Therefore, the new values of nodes displacements are Ukj = Uk~0 + VkjAt i .

(6)

Thus, at the time instant t~- = t~_I + At i we arrive again at the mechanical state of the element, analogous to the previous one, in which its velocity and node displacements take new values. This allows formulating a simple recurrence procedure for the numerical simulation of the problem under consideration. The correlation between the work of the forces applied to an element during the interval At i and the change of the total mechanical energy of an element depends on the value of At L. This relates also to its sum over all the elements n e .The law of the energy balance holds for a single value of At,. Therefore, it is necessary to determine the time intervals At i for each calculation step starting from the equation of the energy balance AEi~ - AWiz : 0 ,

(7) _

ne

where W~ is the work of all outer forces applied to nodes and AEiz - ~j=t(AKi + APi) is the change of the total mechanical energy of the system during the time interval At i . The value of At i , which is a minimal positive root of non-linear Eq. (7), can be easily calculated using successive approximations.

Peculiarities of modelling of one- and multi-dimensional phenomena

Consider propagation of one-dimensional (ID) longitudinal waves in a rod of a variable cross-sectional area Aj. We obtain from Eq.(7) ne

Z A j A x j ( v 2 j - vlj)2(19j - E j A t ~ / A x e ) = 0 ,

(8)

j=l

where Ej is Young's m o d u l u s , pj is density and

Axj is the length of the rod element j. Generally, the ve-

locities of the nodes vlj and v2j depend on the interval At i corresponding to Eq. (3). However, choose the length Axj of elements so that the expression Eq.(8) is satisfied at

we can

pjAx~ / Ej = At 2 becomes constant. In this case

At i = At = const for any values of vlj and v2j and, consequently, for all the types of

loading and boundary conditions. The necessary time intervals based on the law of energy balance for each element are equal to the value At i for the whole system and coincide with the analytical value of the wave speed

Cj =

A X j / A t = ( E j / p j ) " ~ . Analogous situation is typical for other

1D problems with a single

wave speed (torsion of a rod, transverse motion of a string and others ). It is shown that in these cases the method described leads to the exact solution for a given discrete system [4]. By decreasing the dimension of elements Axj and the appropriate time interval A t

=

Axj/cj

we

approximate

the solution for the corre-

sponding elastic solid. Stress-strain relations are satisfied not only at the discrete time instances t i but also

194

B.F.

SHORR

at intermediate moments that allows analysing wave propagation of the strong discontinuities of stresses and velocities over an element. The time interval At i for 2D or 3D wave processes can be determined only for the whole system. By calculating the intervals At i at each step on the basis of the energy balance, we get implicit unconditionally stable solution. Making use of the governing equations of elasticity only at discrete time instants does not cause considerable inaccuracy in the case of short time steps. Since the time intervals determined from the energy balance for individual elements may differ from the interval At i , local oscillations of the period 2At i can appear. They do not affects the stability of computations. However, a smoothing out procedure is utilized to avoid them from final results.. The influence of elements shapes, mesh nonregularity and other

local peculiarities on numerical results

requires a separate consideration.. Using the approach proposed in [4,5], the method can be applied to the analysis of elastic-viscoplastic waves.

Numerical example

Consider generalized plane stress arising in a rectangular plate of the thickness h, the height H and the width L ( Figure 1). Let a short-term point force F0 of the duration

t,mp be applied to the middle of the

fiee side y=0. Three other sides are assumed to be fixed. The quadratic two-dimensional mesh A x = Ay = H / n y , where ny is the number of elements within the height H, is used. We introduce the non-dimensional parameters ~ y OFO H 4.

1

~=x/H,

y=y/H,

t=t/t,

~=u/u,,

- v 2) t. = H

k

,

E

where F0(I- v 2)

u. -

(9)

Eh

It is evident that the minimal numerical value of the impulse duration is equal tO tim p = Ati=I, o r timp ~ 1/ n y . This value is taken below. In calculations, FIG.1 Calculated scheme to the numerical example

damping is neglected. Figures 2 and 3 relate to the plate of the side ratio L/H=2; at the considered instant the waves have not approach yet the fixed sides ( theoretically, at t = 1 ). This corresponds to 2D wave propagation

in a semi-infinite plate. A similar problem in application to plane or axi-symmetric strain is known as Lamb's problem [8]. Its various aspects referred to impact or harmonic loading has

been analyzed

analytically by many authors ( e.g. see [9,10]). Figure 2 shows the displacement fly (x,y) of the plate at the instant t = I. It is seen that a calculated dilatation wave-front takes the form of semicircle in spite of the chosen quadratic mesh. The Rayleigh wave is also clearly indicated on the free surface. It follows from the theoretical analysis of Lamb's problem for an infinitesimal impulse duration timp --+

0 that its solution written as fi = f(~,~) becomes automodelling, e.g. it does not depend on

time [9]. Here,

x = 2 / t , ~ = y / t , "fi = ~ . The displacement ~,(0,~) along the axis y, calculated for

ny=const at three time instants, is shown in Figure 3a. Analogously, the displacement ~ ( 2 , 0 ) free surface is presented in Figure 3b. The analytical

along the

values of the positions of the Rayleigh wave xR

and of the distortion wave-front Yc are also labelled in these Figures. A correlation between numerical

WAVE APPROACH TO FINITE ELEMENTS

195

and analytical results is rather good, excluding a narrow domain near the dilatation wave-front where, theoretically, c~uv / c~y ~ ~ if t,, o ~ 0.

~y 0.10 0.05 0.0

FIG.2 The displacement ~y (x,y) without reflection from boundaries; % = 36, i = ~.~6, -t

1

Changing the number of elements ny, while the impulse duration F0/n~ stays constant, we obtain the same points on the curves as they are shown in Fig.3a,b. Thus, when decreasing element dimensions, numerical results tend to the limit curve which corresponds to the analytical solution in this 2D wave problem.

fi'y 2 3~

0.06

b) 1.03 }.02

0.04

mlm|. 1.01 i'

0.02

).00 0.00 O.

0.2

0.4

YG

O.

1.01

0,8

).2

I

FIG.3 The relative displacement "Uy('~,'y) without reflection from boundaries: a) ~y (0,~) - along the axis y, b) "~y(~,0) - along the free surface; 1) - t ~ 1, 2) - ~ ~ 075, 3) - ~ ~ 05: n> = 48

In Figures 4a,b

the displacement

Uy at the impulse application point versus time i is shown. These Fig-

ures illustrate the effect of repeated reflections from the boundaries o f the plate for different value of the ratio L/H. The first peak

value Oy = 0.942 as well as the two next ones (more than 0.3) are not given. At

the instants t ~ 2, 4,... waves reflected from

the

bottom

side return to the above point. The dilata-

tion and the Rayleigh waves reflected from

the lateral sides also return back to this point. The corre-

sponding time intervals depend on the side ratio L/H. The first return instants coincide with the

196

B.F. SHORR

theoretical values marked by t× and

tR, respectively.

The

displacement peaks

gradually decrease as a

consequence of 2D waves dispersion caused by repeated reflections, however oscillations do not die out as

'

a

0.15

oo io

L

-0,1 5

!

-0.30

fly 015 0.0 -015

0~xR~" t

~J't4

~.-~

~Vl~

"w

',V;

"W'-

"v.

16

.

.

.

" - -

r 7-

-030

FIG.4 The displacement fit' in the impulse application point vs timet- : a - L/H = 4.5, b - L / H = 0.7; n> = 20:

the total mechanical energy of the given undamped system remains constant.

Acknowledgement

The author is grateful to Prof. J.D.Kaplunov for useful discussion of numerical analysis of Lamb's problem.

References

l .O.C.Zienkiewicz. The Finite Element Method,4th ed.,McGraw-Hill, London ( 1991 ) 2.J.T.Oden, L.C.Wellford and C.T.Reddy. A Study of Convergence and Stability of Finite Element Approximation of Shock and Acceleration Waves in Non-linear Materials. US Army Research Office (1976) 3.B.F.Shorr, Mekhanika Tverdogo Te[a 4,144 (1984),(in Russian). Trans.Allerton Press,141(1984) 4.B.F.Shorr and G.V.Mel'nikova. Calculation of Mechanical Systems Using Method of Direct Mathematical Modelling (in Russian),Mashinostroenie, Moscow (1988) 5.B.F Shorr, Archive of Applied Mechanics 65,537 (1995) 6.B.F.Shorr, Mathematical Modelling of Wave Propagation in a Rod of Time-Dependent Material, to appear in 'Mechanics of Time-Dependent Materials' (1998) 7.B.F.Shorr,3rd Euromech Solid Mechanics Conference, Book of Abstracts,329,Stockholm (1997) 8.H.Lamb,Philosophical Transactions of the Royal Society, A 203,1 (1904) 9.G.I.Petrashen,G.I.Marchuk,K.I.Ogurtzov. Trans.of Leningrad State University,21,71 ( 1950),(in Russian) 10.M.J.Kuhn, Geophysical Prospectings, 33, 1103 (1985)