Finite element analysis of incompressible material by residual energy balancing

ht.

.J. Solids

Strucfures,

1974, Vol.

10. pp. 993-1002.

Pergamon

Press. Printed

in Gt. Britain.

FINITE ELEMENT ANALYSIS OF INCOMPRESSIBLE MATERIAL BY RESIDUAL ENERGY BALANCING ISAAC FRIED? Boston University, Department of Mathematics, Boston, Massachusetts 02215, U.S.A. (Received

21 December 1973)

Abstract-Incompressibility is gradually introduced into the finite elements with the mesh refinement in such a way as to balance it with the residual discretization energy in order to ensure fastest convergence to the incompressible solution with best conditioned stiffness matrix to minimize round-off errors in the computations.

INTRODUCTION element analysis of elastic matter, as Poisson’s ratio nears one half or as the material approaches incompressibility two things might happen: the approximation quality of the elements can be entirely lost or the condition of the stiffness matrix may deteriorate indefinitely[l]. Assuming beforehand the material incompressible gives rise to a new independent variable-the pressure and the appearance of an additional equation for zero volume change. Variational principles, from which incompressible finite elements can be formed, including the pressure as a Lagrange multiplier have been brought forth[2] but are no more minimal. Yet, even if the analytical model assumes the material incompressible, computationally it is rarely so owing to the discretization. It is reasonable, therefore, to attempt to introduce the incompressibility only gradually as the mesh is refined in such a manner as to balance it with the discretization error and avoid excessive ill-conditioned stiffness matrix. The purpose of this paper is to demonstrate the practical feasibility of such a scheme. Successful employment of this residual energy balancing technique to the finite element modeling of thin plates and shells is described in[3-51. In the finite

COMPRESSIBILITY

ERROR

In this section we derive the energy error incurred in replacing the incompressible material by a compressible one. The equations of equilibrium of an elastic solid enclosed in D with surface S = S, + .S, are

v2ui + ~~

~ + ~iIG = 0 II

i = 1, 2. 3, in D

with

i?u,

_I

e = 2

+ 1

2x,

&, + -

= Divergence

2x3

t Assistant Professor. 993

U

(2)

994

LSAAC

and with the prescribed

where nl U = (u,. the body terms of

displacements

FRIED

and surface traction

L: = (UT_ UT, pi:)

on

T = CT:: 7-f. r,*>

on

conditions

S,

(3)

S,

(41

= l/v is the inverse of Poisson’s ratio, G the shear modulus. e the volume change. u2, u3) the displacement vector in the Cartesian system (x1, x1, x3), F L=(F;, Fz . FJ I forces, and T = (T,, T, . T,) the surface tractions. The strain tensor aij is given in U by alli Eii

= y--,

Hook’s law which relates solid as

dUi

&I, Eij

zzz -

iixj

“Xi

+

?Xi

the stress gij to the strain

Ti =

ait”,

+

~i2

n, $

energy z(U) of the boundary

[(cl1 :‘,

e2 f -

for the isotropic

i,j=l,2,3

vector n = (n, ,?I~, n3) the traction

where the elastic energy E(U) is explicitly

j&S

cij can be written

cij = Geij

and on any surface with unit normal

The total potential (4) we write as

i,j = 1, 2. 3.

__._A

CIij = 2G

W)=

boundary

is (7)

(Tj3nJ.

vatue problem

in equations

(l)-(3)

given by

E22>’ + (622- e3J2+ (833-

with G set equal 1. The principle of minimum solution to equations (I)-(,4) then

potential

Ei 11~1

energy states that if U is the elastic (10)

E(U) = min n(iij minimization being carried out of the displacement trial field 0 which satisfy equation and for which E(o) < co. In the case of an incompressible material with m = 2 (or v = l/2) the pressure p P = -f(+r -t becomes

independent

and

c22

(111

-t a331

and the direct stress strain relations i=1,22-3

(3)

become with unit G

(13

material by residual energy balancing

Finite element analysis of incompressible

while the equations

of equilibrium

995

become V2u; + g

t

+ Fi = 0

(13)

in which ( )’ refers to the incompressible state. Integration by parts with the use of equations (12) and (13) proves that the same displacement function that minimizes n(U) in equation (8) minimizes also

Ii/(U) =

s,(ml1- e.zz -

+

C&33

+(e3$

which can be concisely

-

811

-

43

&;I

+

+

E>J2

M21

+

+

(E22

m2

-

E33

-

a2

-

+

written

+

(823

E;3y

-

&a2

(14)

-s;1)2]+&Ce-p)2

as

IIe+ /II2

U’11;2++K

$(U>= \I(/-

0

and in which

K=2m+1 Since U minimizes

42

(15)

3m-2’

(16)

ti/(U) 5 $( w

(17)

i&U)

and consequently (18) leading to

which with

/I II e- f

p

2 llello -

0

$llPlo

(20)

becomes

(21) or equivalently

(22) and this is what we were seeking.

996

ISAAC FRIED

To prove that the bound in equation linearly with l/K, we only need find an occurs. An example with an accessible holding an internal unit pressure. In the the pressure are given by

(22) is optimal and that in fact E(U - U’) varies example where equality in equation (22) actually analytical solution is that of a hollow sphere]61 incompressible case the radial displacement U’ and

11’ zz

A

A=

-

r2 in which r is the radial coordinate compressible case

a”h3 4(h3 - n3)’

n=-.&h3 - a3

(23)

and a and b the inner and outer radii, respectively.

A m-2 2f=p+2(m+ l)pr

in the

(24)

from which we find that Efu

assuming

-

u’) =

)p’jK

CT4

unit volume. CONDITION

OF STIFFNESS

According toll] the spectral condition formed from the elastic energy expression

MATRIX

number C,(K) of the global in equation (9) is bounded by

stiffness

matrix

(26) in which N,, is the number of elements per side in the mesh, c1 and t2 numerical constants, i,, the exact lowest eigenvalue of the solid and it, the approximate finite element eigenvalue. Once the approximation is good enough for n1 to be replaced by ;1, the condition number becomes of the form C,(K) = c3

NA

ii,@?8- 2).

127)

Our ability to obtain a good approximation for i_, when I( in equation (16) is large with a reasonable number (say seven) of finite elements may well depend of K because of the term

iK

s

e2du

(28)

D

in the elastic energy expression and we will pay close attention to this crucial approximation question in the following sections. It will be shown that indeed some finite elements are entirely inadequate to cope with K -+ cc while others furnish accurate results regardless of K + m permitting the interchange of i., and pL1in equation (26) for discretization with a moderate number of finite elements. In any event care must be exercised in assigning a value for m = l/v in order to approximate the incompressible state. Lowering v too much away from I/2 might introdu~ a Iarge compressibility error compared with the accuracy which can be obtained with the mesh of finite elements. On the other hand taking v too close to l/2 might either ruin the discretization accuracy or produce excessively ill-conditioned stiffness matrix. Poisson’s ratio v = 1/m ought

Finite element analysis of incompressible

material by residual energy balancing

997

to be increased gradually with the mesh refinement to balance the compressibility error with the best discretization accuracy and in such a way as to keep C,(K) as low as possible. Exploring ways to accomplish this is the subject of the following sections. RESIDUAL

ENERGY

BALANCING

We introduce the residual energy balancing technique by way of carefully examining the problem of a radially pulsating hollow sphere free of surface and body forces. Presently e = u, + 2ulr and to form the element stiffness and mass matrices E and K which we write here as

E= Glab [A

(2% we need the elastic and kinetic energies

(u:+2$)]ridr

e2+

(30)

and K = ;Ibu’r2 (I

dr.

The element stiffness matrix k derived from E in equation k=zk,

(31) (30) we prefer to write in the form

z = G/(m - 2).

+Gk,,

(32)

Assembling the element stiffness matrices k into the global stiffness matrix K and the element mass matrices m into the global mass matrix M produces the algebraic eigenproblem (zK, + GK,)x

= AMx

(33)

and the lowest ,I is given by xTK,x

I=min

xTK2 x

z=+G-

xTMx

x

(34)

The term zxTKlx/xTMx is the compressibility energy and should vanish as z + co, but we have to distinguish here between two distinct cases: the one with a singular Kl and the other with a non singular Kl. Suppose first that discretization is achieved by assuming a polynomial approximation of degree p and in particular let p = 1. For this, the (lumped[7]) element matrix m and the element matrix k, for an element situated between r = rl and r = r2 become k

a2 + 314

ap - 314

up - 314

p’ + 314? p

1

in which ra = (rl + r,)/2, h = r2 - rl, CL= 1 - r,/h and b = 1 + r,/h. The element matrix k, is non singular and since[8] . xTK,x . yTk,y mm, > minx xx Y Y’Y

(35) stiffness

(36)

so will Kl be and we have that min x2

IJSS Vol. 10 No. 9-E

= 0(h2).

(37)

ISAAC FRZED

998

This means that incompressibility is asymptotically achieved as h -+ 0 but for any given h accuracy is ruined by increasing z indefinitely since then zxTKlx/xTMxwill grow without bound making R = 0 or the element perfectly rigid. In boundary value problems the situation is the same. There the energy error E(u - it) in the finite element solution ri (this is also the error in eigenvalues in eigenproblems[9, IO]) i, E(u - 2) = z

lb(e - B)‘r’ Ll

dr +

(.’[(ur- ii,)’ + -$(u -

* il

1

2)’ r2 dr

(38)

and if P can not assume any constant value unless 6 = 6 then as 2 -+ X, c -+ 0, S + 0 but also t; + 0 making the element ever more rigid. This situation is entirely analogous to that in thin inextensional curved beams where if the extensional strain can not assume exactly the rigid body mode (i.e. assume any constant value) then reducing the thickness makes the extensional energy vanish only with vanishing displacements and the element becomes excessively rigid. in view of equation (37) the discretization error can be balanced with the compressibility error by choosing. for first order elements, z = O(h-‘) since then the compressibility error derived before is O(h) and so will be zv*K~s,i.~~M.~. But it is rather a poor approximation since with linear elements we expect quadratic accuracy for the eigenvalues. In boundary value problems the situation is again the same. The second integral in equation (38) can be O(h’) but the first only O(h). Numerica experiments confirm too that no more than O(h) accuracy is obtainable in this case. To avoid the absolute rigidity of the element with z = CC we need singularize Ki and can achieve this by choosing the shape functions to enable 8 to assume any constant value (in curved beams or shells this is analogous to the introduction of exact rigid body modes). Thus we start by assuming a polynomiat assumption for P (as in[S]) and integrating. For

this yields .

bo

rr=~r;?+$i,r

(40)

and with these shape functions both k, and Ki become singular. Now minimization in equation (34) with z + cc will seek the mode corresponding to the zero eigenvalue of K, can become and zxTKlx~xTMx will remain finite for any 2. Since presently :-xTK,/_xTMx arbitrarily small for any z we expect to find a mode that will yield simultaneously 0(/t’) accuracy for the second integral in equation (38) and 0(/z*) also for the first integral in equation (38). Another technically interesting possibility to reduce the rank of K,, or at least hasten the diminishing of its lower spectrum, is by assuming a polynomial expressing for zi and using numerical integration. If ii is assumed a polynomial of degree p then following the rule of[l I, 121 and integrating each energy term by p Gauss points will singularize k, and consequently K,. For instance, with p = I, analytic integration produce the k, (2 x 2) matrix of rank 2, but a one Gauss point integration reduces this to rank I. The beneficial effects of numerical integration in the finite element analysis of thick plates and shells were earlier discovered and described in[13-151.

Finite element analysis of incompressible NUMERICAL

material by residual energy balancing

999

EXPERIMENTS

of these experiments is to verify, at least for the model problem of a freely vibrating hollow sphere made of incompressible material, that by approximating the d~s~~a~~me~tby a polynomial of degree p and integrating (in one d~rne~sio~) the eBnergy terms with p Gauss points the fuh energy rate of convergence O(@) can be obtained by taking The purpose

.z = cNSZp

(41)

where c is a positive constant and ,Vethe number of elements (in mult~d~mensiona1 problems N, should be replaced by N,,, the number of elements per side) A reasonable value for c keeping C,(K) as low as possibte will also emerge from these experiments. The fundamental eigenvalue (square of frequency) /i of the radially vibrating sphere is given by I. = 46

u2 c a6 + 6= &“/,$

142)

which for G = -f (unit elastic modulus), a = _Eand b = 1 becomes I = y = 9.333. First, the sphere is discretised with first order (p = l), two-nodal-point elements and the stiffness matrix created with one Gauss point integration In accordance with equation (41) z is taken in the form z = cA$?and the calculated eigenvalue L is shown in Fig. I as a function of the number of etements in the mesh and for different values of c. The variation of E.with c is also shown in Fig. 2 far N, = 15 and indeed between 2 and 50 the accuracy of the eigenvalue depends only weakly on c A reasonable value for c is seen from Fig. 2 to be 5. Figure 3 shows the reduction of the relative error SA//1.in A vs N, on a logarithmic scale to show that the optimal rate of convergence O(N;‘) is actually achieved. For the choice z = C-N:it is found that C,(K) = 3e.g: (Fig. 4(a) shows C,(K) vs N, for c = 5 on a logarithmic scale) and hence even though overestimating c has Me e%zct on the accuracy, it raises the condition number of K and should be kept down. Actually to obtain three decimals accuracy in 1 with linear (p = 1) elements and with c L=5 about 10 elements are required with C,(K) = 1.5 lo’, which is well within the capabilities of modern computers.

.?

8

$2

fC

Fig. 1. Convergence of computed fundamental eigenvafue .&of a radially pulsating incompressible hollow sphere discretized with first order [I” - 1) elements integrated with one Gauss point as a function of the number of elements in the mesh and for difkent values of c in t = cN,2.

ISAACFRIED

1000

2

5

25

50

Fig. 3. Variation of h from Fig. I with c for N, =--15.

Fig. 3. Convergence of h from Fig. I for c -~ 5 vs N,. Logarithmic scale demonstrates attaj~ment of the optimal rate of convergence o(iV;%

the

Increasing the order of the element must be accompanied by an increase in the rate of growth of z and consequently of C,(K). In general with a polynomial of degree p (and p Gauss integration points) z = cNz* and C,(K) = 0(NzP’2 ) which is unlike the compressible case in which the rate of growth of C,(K) does not depend on p. Nevertheless with higher order elements the error decreases like 0(N;2p) and fewer elements are needed for a given accuracy. Hence in spite of the fact that C,(K) grows faster with higher order elements, less elements are needed for the same accuracy and the value of C,(K) might still be lower than that for the low order elements. The fofiowing numerical example shows that this is actually, the case. It consists of computing A with a quadratic (p = 2), three-nodal-point element, a 2 point Gauss integration and z = CM:. Figure 5 shows the rate of convergence of i, for different values of c and this is seen indeed to be O(NT4). A good choice for c is 10, but again overestimating it changes the discretization accuracy only slightly. The condition number grows now with N;Pas seen from Fig. 4 (6) (drawn for c = IO). Yet per accuracy the

Finite element analysis of incompressible

material by residual energy balancing

1001

Fig. 4. Spectral condition number Cz(K) of global matrix K for hollow sphere discretired with (a) linear elementsintegrated with one Gauss point and with z = 5N,’ and (b) quadratic elements integrated with two Gauss points and with z = lON, 4. Theoretical rates of growth 0(Ne4) and 0(Ne6) are verified.

Fig. 5. Convergence of h with quadratic (p = 2) elements. Optimal rate of convergence 0(Nem4) is obtained with z = cNc4.

higher element leads to a better conditioned global stiffness matrix. To obtain a 3 decimals accuracy in 1 only 3 quadratic elements are needed (they are hence also more economical computationally) and C,(K) = 2 x 104, less than that with linear (p = 1) elements.

1002

ISAAC FRIED REFERENCES

1. I. Fried, Influence of Poisson’s ratio on the condition 2.

3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13.

14. 15.

of the finite element stiffness matrix, Znt. f. Solids Struct. 9, 223-329 (1973). R. L. Taylor, K. S. Pister and L. H. Herrmann, On a variational theorem for incompressible and nearly incompressible orthotropic elasticity, Znt. J. Solids Strurt. 4, 875-883 (1968). 1. Fried, Shear in Co and C’ bending finite elements, Inr. J. Solids Struct. 9, 449460 (1973). 1. Fried and Shok Keng Yang, Triangular, nine-degrees-of-freedom, Co plate bending element of quadratic accuracy, Quart. Appl. Math. 303-312, (October 1973). I. Fried, Shape functions and the accuracy of arch finite elements, AZAA J. 11 (3), 287-291 (1973). A. E. H. Love, The Mathematical Theory of Elasticity, 4th. Edn., Section 98. Dover, New York(1944). I. Fried, Lumping the finite element mass matrix with no accuracy loss. To be published. I. Fried, Bounds on the spectral and maximum norms of the finite element stiffness, flexibility and mass matrices, fnt. 1. Solids Struct. 9, 1013-1034 (1973). I. Fried, Accuracy of finite element eigenproblems, J. Sound V&r. 18 (2), 289-295 11971). I. Fried, Boundary and interior approximation errors in the finite element method, J. App. Mech. 40 (4), 1113-1117 (1973). I. Fried, Accuracy and condition of curved (isoparametric) finite elements. J. Sound V&r. 31(3), 345-355 (1973). I. Fried, Numerical integration in the finite element method, Computers & Structures, to be published. W. P. Doherty, E. L. Wilson and R. L. Taylor, Stress analysis of axi-symmetric solids using higher order quadrilateral finite elements. Structural Engineering Laboratory Report No. SESM 69-3, Univversity of California, Berkeley (1969). 0. C. Zienkiewicz, R. L. Taylor and .I. M. Too, Reduced integration technique in geaeral analysis of plates and shells, Int. J. Num. Meth. Engng 3, 275-290 (1971). S. F. Pawsey and R. W. Clough, Improved numerical integration of thick shell finite elements, ht. J. Num. Meth. Engng 3, 575-586 (1971).