The viscoelastic space-time element

CHl45-794996 Pergamon

THE VISCOELASTIC

SPACE-TIME

53 00 - .oO Journals Ltd.

ELEMENT

ADAM PODHORECKI Bydgoszcz, Poland

(Received 4 June 1985) Abstract-The way of modeling of viscoelastic medium in the method of finite space-time elements has been presented in this work. This work contains the general algorithm of stiffness matrix building of space-time elements for any state of stress. The presented algorithm is characterised by the great simplicity preserving the values of classical method of finite elements.

where

I. INTRODUCTION

The methods of space-time elements (MECZ) have been presented in the papers of Kaczkowski[l-31. The above-mentioned papers have illustrated the practical utilization of that method (bar and string vibrations). Besides, the influence of the external damping of viscous type with longitudinal vibration of the bar[4] and the internal damping according to the Kelvin-Voigt model have been analyzed[3]. The internal damping (material) in the finite element method and other numerical methods are frequently described by means of viscoelastic KelvinVoigt medium. This model does not thoroughly describe the characteristics of viscoelastic real bodies and this is connected with the necessity of forming more complex models. Several models can be distinguished in the linear viscoelasticity[5]: Hook’s model [Fig. I(a)]. Kelvin-Voigt’s [Fig. I(b)], IMaxwell’s model [Fig. l(c)], Zener’s model of the first and second kind [Figs. l(d) and I(e)], and Burgers model [Fig. l(f)]. The general formula describing the stiffness matrix of space-time element for the arbitrary linear viscoelastic medium has been introduced in this work. In particular, such well-known models as Hook’s, Kelvin-Voigt’s. Maxwell’s, Zener’s. and Burgers have been considered.

Sij

=

Ui,

-

(+kkq

(2.2) yjj

=

E;j

-

f&/jr!&.

The tensors Sij and y;j in succession are stress and strain deviator stresses Sij cause of normal stresses (a” = l/3aIk_), which causes the change of volume. There are two material constants of nondilational strain: modulus TVand compressibility modulus K in eqn (2.1). K=

E 3(1 - 2”)

(2.3)

(E-Young’s modulus, u-Poisson’s ratio). Analogically to formulas (2. I), a general relationship between the state of stress and strain for viscoelastic medium can be presented[5. 61. PI(D)Si,(X,

I) = ~Pz(D)yjj(~. f), D=;,

P,(Dbu(x. 2. THE EQUATION

f 6ij

(2.4)

1) = P.~(Dku(x. t).

OF STATE

where P(D) arc linear differential operators in relation to time I. Equation (2.4) can be written in the developed form utilizing the indicator calculus[j.

The generalized Hook’s law can be written down in the deviation form[6].

61. al

&--i--&:,& Q2

e*

Q3

II

‘11

E,

E2

4

‘l2 E3

Fig. I. The diagram of linear viscoelasticity of models: a) Hooke’s model. b) Kelvin-Voigt’s model. model. d) Zener’s model of the first kind, e) Zener’s model of the second king. fl Burger’s model.

c) Maxwell’s

535

536

ADAYPODHORECKI

P,(D)P,(Dbi, +

= PL(D)P,(Dk,,

!JBij[PI(D)P4(D)

-

P?(D)P3(D)lf!&.

(2.5)

or the matrix calculus P,(D)Pj(D)a,

3. THE VIRTUAL

(2.6)

= E*e,.

where E* = fP,(D)P,(D)L,

We can also accept other assumptions[5. 61. for example, Poisson’s ratio is constant and does not depend on time or viscoelastic body and the characteristic of incompressibility (then PI = 0).

+ #‘?(D)P2(D)L2

(2.7)

is the viscoelastic matrix. The descriptions of matrices LI and LZ and a, and l, depend on the analyzed state of stress and strain; for example, for the spatial state of stress we have

FOL’R-WORK

EQUATIOS

Let us consider the viscoelastic body, which undergoes the influence of external forces pi(x, I) which change in time. Under the influence of these causes the displacement field 11,(x, r), and strain field eij(x, t), and stress field oi,(x, f) are formed. The equation of motion does not depend on mechanical characteristics of material[j, 61: Uij,j+

Xi -

The boundary character:

0,

pii; =

conditions

x E V,

I 2 0.

(3.1)

have static and kinematic

L, =

s,,.

ui,(x, f)nj(X) = pi(X, f).

x E

I(;(X. t) = ir;(x. f),

x E s,,.

(3.2)

The foregoing equation can be written in the matrix form

(2.8)

a,;u.,+ x na,

a, = co1 E, =

{uII, U22r U33, UI2t U137 U23 19

Assuming that during the extensive distension or compressing the viscoelastic medium behaves in the same way as the elastic body, we have the equation which has not changed eqn (2. I): =

1,

P&I)

=

3K.

(2.9)

x E S,,.

aI

I

I

a1

al

1: I I-I-I ax,! ’ I ax> , dX3 , 4___~__t_-T_7_--

a,'=

a, I a l1-I I ax, 1 ax2 ; T-ay----1-1-1 , ax, , --I

model

Pl(D1

P2(Dl

1

2/J,

Hooke’s

Ketwn

-Voqt’s

1

Maxwett’s

‘*&x

Zener’s Of the

first

Zener’s of the

second

kind kmd

1*X,%+ l.l~~-a,l-&

2N 3

b&-)

AZ= * A,=*

2XZ*rl&2e41*a,$$

zGJ,($

) *a&

A,:% A,+.

Id

I 1 ah ys-t-aL I---I--I a.r, I a.rl

Table I. The description of the operators P, for different viscoelastic models of

(3.3)

x E s,,,

= p,

The dimensions and descriptions of particular matrices which occur in compounds depend on analyzed states of stress. In the case of spatial state of stress we have in succession:

Differential operators P,(D) and P?(D) are described in this case in the same way as the uniaxial state of stress (compare Fig. 1) which is included in Table 1.

Nome

x E v.

u = u,

co1 {El,* Et29 E33, El?, El39 f231.

P3(D)

pii = 0,

.

(3.3)

The viscoelastic space-time element

531

and accepting the arbitrary virtual displacement ii. in that way, that the kinematic conditions on the surface S, are not changed we receive the equation

ET[f ,P,X - P,P.,(pii)l

x =

-I-

COl{X,‘ x2. X,],

Let us multiply by P,(D)P$D).

p =

the first two equations

P,P&,o,

+ P,PJX

P,P,na,

i

COliPI.P2. PSI.

- PIP&

bilaterally

= 0.

(3 5)

= P,P2p.

In order to derive the virtual four-work equation we have to set the virtual displacement 611(x, t) = U(x. I) on the body of volume V which observed in time from f; to ta. Let us multiply eqns (3.5)$ by ii(x, f) and let us integrate this product on the body field V and according to time. Let us also multiply the boundary static condition (3.5) by li(x, t) and let us integrate on the body field S,, according to time. From these integrations we create the expression

+

s,, U’(P,Psp)d.s

z:(P,PJa,) I $’

” IfL

- P,f3p]

iiTar(P, P,a,)

= @W(P,

[iifP,P>(pbk)

=

PJU,)l G = a.~~~

(3.7)

eqn (3.6) can be brought to form:

J I‘ {J 1 +

iir[P,P3X

-

-

iiTPiP3(n&i)]dV ir[P,P,(pi)]

Finally,

P,Pdpii)j

dV

@]iiT(P,P~a,)l

- j- ZI(P,Psa.,) $’

(3.11)

we receive the equation - ii:P,P,(pid] cT(P,P,X)dV [Z:(P,&u.r)

dV +

- iirP,P,(pti)]dVdt.

(3.12) The left-hand side of the foregoing equation can be interpreted as a “work of external impulses” performed on the virtual displacement while the right-hand side expresses the “work of internal impulses” performed on virtual strain. it is called “the virtual four-work equation” by Kaczkowski. 4. THE ELEMENTAL

STIFFNESS

VLSCOELASTIC

I$

dV dt.

dS

- z.:iffP,P3u,),

I‘

(3. IO)

‘iiT[PrPJ(pii)] dV dt

=

dr = 0. (3.6) 1 After emploving . - _ the formula on the calculation of the derivation of the product of function

dr = 0. 1

v

+

js,iiTIP,Ps(no,)

dV

The second expression of the first integration (3.10) can be partially integrated according to time:

[$P,Ps(pi;)

- P,P,(pii)l dV -

dV

MATRIX OF THE

ELEMENT

The area of displacements, strains, and displacement velocity is described in the way which is characteristic for the method of finite elementsf71:

dV

dV

u(x, 1) = N(x, r)&. -i

ii7[P,P&ia,)] $8,

dS

E, = &u = a,NS,. = Bg,., (3.8)

After making a Gauss-Ostrogradski transrormation[6] on the second relation derivative

fv

~~~~~~P~P~u~~~dV =

~~~P~P,(nu~~~ dS,

(3.9)

i = (&It)u = (&r)N&

(4.1)

= I%, ,

when N is the shape matrix containing time, space functions, and 6, contains nodal displacement of the finite space-time element (SKECZ). We employ the equation of the state of viscoelastic medium (2.6): P,P3a,y

= E*e., = E*Bg,..

(4.2)

AD;\.w PODHORECKI

538

Then. we give the virtual displacement 3,’ to the finite viscoelastic element (SKECZ). We can get the virtual displacements and strains from formulas (4.1) substituting their virtual increments. The virtual four-work equation takes new form for finite viscoelastic element of R = V.t area or R, = S,;t (compare Ref. [2]). iiT(P,P3X) dR + =

Gr(P,P;p)

J-Jo

[g:(P,P,u.,)

dR,

- I’P,P,(pti)]

do.

but Table 1 contains the factors CY. Assuming that the medium has the same viscoelastic characteristics with dilatational strain and nondilatational strain. the stiffness matrix has formed (compare Ref. [3])

(4.3)

After taking advantage of eqn (4. I ), we receive the already known relation F, = K,6,.

-

F:,

(4.4)

K,. =

II0

[a,B’EB

+ a2BFEh

+ a?B’EBj dR

in which the expression K, =

describes F,P =

III1

R II

[B’E’B

-

rirT(P,P,pr;r,] dfl

(4.5)

the stiffness matrix, and NT(P,P,X) +

where E is the elastic matrix (compare Ref. [7]). The matrix F,P is described by formulas (4.6) or (4.7).

da

IIfl,

NrU’,P,p) d%

s.

(4.6)

THE SHAPE FUNCTION SKECZ

The shape function for finite viscoelastic element is selected in the characteristic way for FEM[7]. One can conclude from formulas (4.7) or (4.8) that the variation of shape function according to time must differ in degree depending on the accepted model of viscoelasticity. So in case of Hook’s and Kelvin-Voigt’s models the linear variation according to time is sufficient. But if the vis-

-the nodal impulses SKECZ from the volume loading X and the surface loading p. The viscoelastic matrix E* is characterized by the relationship (2.7). If the assumption about the dilatational strain (9.9) is obligatory based on Table I containing operators P,(D), we can express these formulas in another way:

Table 2. The description of coefficient factors a, for different viscoelastic models

I

I

I Name of model

a1

‘s

1

Ketvln - Votgt ‘5

1

Hooke

Maxwell’s ZWlW’S of the first

d:

Al

A2

kind

81

81X

x3 ^

_

The viscoelastic

Fig. 2. Examples of finite space-time elements: a). b) space-time elements of bar (beams. columns): cl. d). e) space-time elements of space girders (plates. disks, shells).

space-time element

indicate one after another the static stiffness matrix and inertia matrix. (2) The case of elastic quahties with nondilatational strain (4.7): tN:N, + a,N:l;1’, + acN;ti,)dt@So

K, = +

coelastic

medium

is described

according

539

to Max-

f I

(cu,NfN,

+ a2N:&,

f

a,N:&fdr9S,,

well’s or Zener’s models at least square functions (ti:li, + a,&‘:& -t cqfi:N,) dr B hl,.. must be accepted, and in the case of four-parameter models (Burgers’) they must be cubical. (5.4) The shape and number of nodes can be arbitrary (Fig. 1). where In the case of uniform digitization along the time $0 = K~L/,%J%,dV. coordinate (Fig. 31, the shape matrix N(x, t) can be (5.5) divided, using the tensor product of the matrix Q9 B.:L2B, d V (compare Ref. [3]), into two matrices, one of them S, = P. depends only on the time variable N, and the other one depends on space variables N,. indicate one after another dilatational and nondilatational static stiffness matrices. For the viscoelastic element showed in Fig. 3(b) N(x, I) = N, C3N.,. (5.1) the matrix of time is described with the help of square functions (Fig. 4) Formula (5. I) gives us the possibility of partial differentiation according to the time made on the shape matrix. So stiffness matrices can be written in a different form: -1STSl. , (5.6) ( I) The case of identical viscoelastic characteristics with di~atationai and nondilatational strain:

FE

Fig. 3. The idealization

of space-time:

a) the example element.

tn

.1*1'

moment

of time digitization:

b)

the finite space-time

540

ADAM

Fig. 4. Graphic illustration Now we can calculate

derivatives

PODHORECKI

of time function

of shape and their derivatives.

N,: N;N,dt

=

(5.7)

N:N, dt =

and we can integrate products depending on shape matrices N,, which are included in formulas (5.2) and (5.4). N,T&dt =

fi:ti,h

d-i = 3

NfN,dt = (5.8) The adoption of square shape functions allows describe the viscoelastic medium using the models of Hooke. Kelvin-Voigt, Maxwell. and Zener. Utilizing the expressions (5.81, the stiffness matrices SKECZ (5.2) and (5.4) adopt the following forms successively:

us to

N:N,dt

=

%I

The viscoelastic space-time element

I

-

__--_

-

&a,h

?a:+ ia3

&a,h

ja, - -ji;a, ------________

-----__f--__‘-__

- fa2 +

’

4

1 _____I

z

a3

_-----s-e-_

- halt1

+

I

_-- --_--

K, =

I&a,h

+ Aa, + $

8 fJa,h - 3h

-_

I a3 1 &alh

___

I-

-

aa,f

-a,

&a,h

+

%a2+

Aa3 8 S,

I I

a3

____

- ~+a?- i

1

I

&a,h

3h

311 I__--___-_--__ I I &a,h + ia? + - a3 a3 / 3h i

-

K, =

(5.9) ----_----_--__

-

7 -_a I ; --_+4 6h h’ ’ ; 3h _------_-+--_--I 8 -- 4 3h

_-------+-----_

!

z

Matrixes S, and M, do not depend on the viscoelastic model. They are set in the typical way for the classical method of finite elements, that is for finite space element (FE). The structure of a stiffness matrix for adopted SKECZ as in Fig. 3(b) is as follows:

(5. IO)

These factors are calculated directly from formulas (5.9). The dimensions S, and M, and also K:’ depend on the number of nodes FE on the number of degrees of freedom of the node. The time dimension SKECZ should be selected in such a way that the stability condition is fulfilled. The stability condition depends on self-prevailing frequency. (5.11)

h<% Wmax

6. THE EQUATIOSS

where Ky(r, s = n, p, q) is the submatrix of matrix K, containing the stiffness factors at the moment “T” which are made from unitary displacements the moment “s”.

OF MOTION

For the digitized time space as in Fig. 3, the system of equations has the form

at K6 = F,

(6.1)

ADAM PODHORECKI

where

t6.21

K=

is the global

stiffness

matrix,

the stiffness

matrix

of the

r from

moment ment

the factors whole

the unitary

K”

indicate

structure

displacements

where

at the p”

at the mo-

s,

= (I(W) -

1F” _ (K’,‘,) - , F”. 1.

* = _(E;‘,‘,)-lK”” 6 = co1{6”,

&“.I,

6’.

6’,‘.

.}.

6’,

+ ,Ii”4)-l~‘..i~

B = (KW) - 1K”” _ (h;‘,‘,) - 1Kti’,P.

(6.3)

C = (K”‘,)- 1K,‘L’. the

matrix-the

ment

column

vectors

containing

of all the structure

the

displace-

in successive

,, = (K”“) - 1K,‘P

time

F’-’

moments.

F = col{F’,

F”.‘, F’, F’,‘, F’, :.

= (KC!“)-IF’-,

_ (Kc,‘,)-IF’-1,

_ (KU,) - I ~“‘,6’- 2.1- I,

.}

(6.4)

60.1 = B;I(Kw-IFO contains

nodal

impulses

For the initial

conditions

system

of equations

formula

(compare i

from

loading. known

(6.1) changes Ref.

beforehand

+ B;‘A,6”,

the

into the recurrent

6’-

1.;

_ _

(6.8)

B;‘(K”‘J)-IF’-’

[3])

= &-‘p’

_ ~;l(~qp~-I~O.l

_ B;‘(K‘,J-)-IF’-,_’

_ B;‘Cl&‘-2

i = I, 2, 3,

,

_ B;‘(K”,‘)-‘K”“6’-‘.‘-’

(6.5)

+ B;‘(A,

where

+ D,)S’-‘,

Al = (K’,“) - 1K‘,” _ ( K”” ) - 1K”“. A’ = [_;;!rj,

p!

=

r=

K.

F’ _ rKPYfj’_--_-__--------~~---~-

I

for

i=O

for

i2

change

(6.6).

can

if.

for

of medium

ary conditions

BI = (K”“)-‘K”‘,

,

,.KPYfj-

‘.i(KfJ(J ---_

c,

+ ,-KPP)fj’ --_------

way,

that

displacements P, . from

(KC,,,)- , K,“’ , ( K””

) - 1I<,‘,‘,

7. ILLUSTRATIVE

submatrices,

which

be

or

constant

example,

or inertia

manner

_ (K‘,,‘)-‘K’,“,

(6.6)

the

occur

they

can

The

bound-

straight

section

The

tudinal

is, we must

in intermediate the sequence

first

eliminate

moments:

6”,‘.

in the

form

space

equal

to A, and density

operating

vibration

The can also be established

EXAMPLE

bar is of length

area equal

5(a)].

viscoelastic

or kinematic

change.

The recurrent

=

D, = _

1.

Ku and other

in time

equalities

[_;;-/_I;]

F i.i+ I - KVU&

I

formulas

another

1 _

0

The matrix in

=

loading

p(.r,

6’=

= R-‘p B-IF’-’

-

t) evokes

longi-

u(x.. t). matrix

of Hermite’s

of shape

is accepted

multinominal

(compare

in the Refs.

[Z

and 41).

6’.‘,

of equations -I

6’

to 1, of cross

equal to p [Fig.


I.

(7.1)

R-‘A#‘, _ R-‘(A

+ D)#-!

_ B-‘CSj-‘, It has been

i = 2, 3. 4,

.

(6.7)

viscoelastic

decided qualities

to consider with

the case of equal

dilatational

and

non-

The viscoelastic al

p(r.tl

543

space-time element

scribed on the of formula (5.9):

X.”

7---.--

I

[EA.?

1

:

:

I

b) 0 !

U

12 :

20

" :

:

(

20

Fig. 5. Example: Longitudinal bar vibration: a) the sector of bar of length 1: b) space digitization of bar: c) finite element.

dilatational

The static matrix of stiffness and are derivated from eqns (5.3).

strain.

inertia matrix

EA

B;EB.,ll d( = c

(7.2)

Formula (4.7) defines the matrix of node impulses operating at the moment ‘y’ from external loading: KS? =

&a,h

- ia,

(

F& =

z - G

a3

Sik >

4

x

(pj + a,$

+ a5p’)h dr 9 d{.

(

(7.3)

If in the area SKECZY p = const. then impulses for the accepted stiffness matrix has the form

K?! =

(

-3h

1

+ $a2 - -a3

&u,h

- ($-

)

2 - 7 a4 Mik3 h> S;k

+ ;:)7 4 Mi!. 1 i,k

= 1.2 (7.6)

and formulas

s,, = s::

(7.2):

= E

,

M,, = Mz2 = SpAa.

s 2, -__E

s,: =

The structure of matrix K’” is built in the same way as the global matrix for the static problem

K’”

(7.7)

.

=

20 ’

M,? = Mz, = 4pAa.

(7.5)

nxn

1, s = a, 9, p n -number The particular terms

FE [compare Fig. 5(b)]. of matrix (7.5) are de-

The coefficient factors a depend on viscoelastic model and are contained in Table 2.

54-l

ADAM

PODHORECKI

I “1

6

0

Ii

Fig. 6. Longitudinal and Kelvin-Voigt’s

18

Z’L

vibrations of the bar: a) The course of vibration of the braced ending for Hook’s models; b) the course of braced vibrations for Maxwell‘s and Zener’s models; c) the manner of acting of the power on the braced.

Now we can start counting displacements in the successive time moments according to recurrent formula (6.5) or (6.7). The numerical example concerns longitudinal vibrations of bracket of length equal to 1 = 16, 0 [ml and of the cross-sectional area equal to A = 0, I2 [m’], density equal to p = 7500 [kg/m’), and modulus of elasticity equal to E = 2 x IO” [N/m’]. This bracket has been divided into finite elements of length equal to 8, 0 [ml each (a = 4, 0 [ml). The time dimension SKECZ has been established according to the formula (5. I I):

h < YL

= $

= 1, 7

t 10~3Esl

3b

X

10-3[Sec],

8. CONCLUSION

The way of modeling the viscoelastic medium in the method of finite time-space elements has been presented in this work. The work contains the general algorithm of the stiffness matrix building SKECZ for any state of tension. The presented manner of calculating is characterized by great simplicity, preserving the values of the classical method of finite elements. After performing more detailed quantitative analysis we can build computer programs for more complex engineering problems. REFERENCES

wmax

(7.8) For the calculations h = I. IO-’ [set] has been accepted. The loading consists of suddenly applied force P = 30 [kN], at moment r0 = 0 and detached at the moment t3” = 30 x IO-’ = 0, 03 [set] [Fig. 6(c)]. The vibrations of the end of the brace for different viscoelastic models are shown in Figs. 6(a) and 6(b). The received results of the calculation in the case of Hook’s and Kelvin-Voigt’s models are highly exact referring to exact solutions, both according to the quantity the error of displacements amplitude and vibration period of about 2%~~and according to quality.

The method of finite space-time eleI. Z. Kaczkowski, ments of structures. J. Tech. Phys. 16, 69-84 (1975). The method of time-dependent finite 2. Z. Kgczkowski, elements. Arch. Civil Engng 22, 365-378 (1976). 3. Z. Kaczkowski, General formulation of the stiffness matrix for the space-time finite elements. Arch. Civil Engng 25, 351-357 (1979). 4. Z. Kaczkowski and M. Witkowski. Accounting for external damping in the method of finite space-time elements. Arch. Civil Engng 23, 247-254 (1977). W. Nowacki, Theory of Creep. Arkady, Warsaw (1963). W. Nowacki, Dynamics of Strucrures. Arkady. Warsaw (1972). 0. C. Zienkiewicz. The Finife Nement Merhod. Arkady. Warsaw (1972). J. Langer, Spurious damping in computer-aided solutions of the equations of motion. Arch. Civil Engng 25, 359-369 ( 1979). 9. W. Derski and S. Ziemba, Analysis of Rheological Models. PWN, Warsaw (1968).