2

’1 ’4 0

7.4.5 TANGENT MATRIX

The linearized virtual work Equation (7.30) can now be discretized for el-

ement (e) linking nodes a and b (see Equation (7.32) and Figure 7.2(a)),

18 DISCRETIZATION AND SOLUTION

δva

δva δva

ub ub

a

a a

b b

(b) (c)

(a)

FIGURE 7.2 Assembly of linearized virtual work

in terms of the total substi¬ness matrix K ab obtained by combining Equa-

tions (7.36), (7.45), and (7.50) to give,

(e)

DδW (e) (φ, Na δv a )[Nb ub ] = δv a · K ab ub ;

(e) (e) (e) (e)

K ab = K c,ab + K σ,ab ’ K p,ab (7.51a,b)

The assembly of the total linearized virtual work can now be accomplished

by establishing (i) the contribution to node a from node b associated with

all elements (e) (1 to ma,b ) containing nodes a and b (see Figure 7.2(b)), (ii)

summing these contributions to node a from all nodes b = 1, na , where na is

the number of nodes connected to node a (see Figure 7.2(c)), (iii) summing

contributions from all nodes a = 1, N . This assembly process is summarized

as,

ma,b

DδW (e) (φ, Na δv a )[Nb ub ] (7.52a)

(i) DδW (φ, Na δv a )[Nb ub ] =

e=1

e a,b

na

(ii) DδW (φ, Na δv a )[u] = DδW (φ, Na δv a )[Nb ub ] (7.52b)

b=1

N

(iii) DδW (φ, δv)[u] = DδW (φ, Na δv a )[u] (7.52c)

a=1

This standard ¬nite element assembly procedure can alternatively be

expressed using the complete virtual velocity vector given in Equation (7.20)

together with the corresponding nodal displacements, uT = [uT , uT , . . . , uT ]

1 2 n

and the assembled tangent sti¬ness matrix K to yield,

DδW (φ, δv)[u] = δvT Ku (7.53a)

where the tangent sti¬ness matrix K is de¬ned by assembling the nodal

19

7.5 MEAN DILATATION METHOD FOR INCOMPRESSIBILITY

components as,

···

K 11 K 12 K 1n

®

···

K K 22 K 2n

21

K= .

. .

..

°. . .»

.

. . .

···

K n1 K n2 K nn

7.5 MEAN DILATATION METHOD FOR

INCOMPRESSIBILITY

The standard discretization presented above is unfortunately not applica-

ble to simulations involving incompressible or nearly incompressible three-

dimensional or plain strain behavior. It is well known that without further

development such a formulation is kinematically overconstrained, resulting

in the oversti¬ phenomenon known as volumetric locking. These de¬cien-

cies in the standard formulation can be overcome using the three-¬eld Hu-

Washizu variational approach together with an appropriate distinction being

made between the discretization of distortional and volumetric components.

¯

The resulting independent volumetric variables p and J can now be interpo-

lated either continuously or discontinuously across element boundaries. In

the former case new nodal unknowns are introduced into the ¬nal solution

process, which leads to a cumbersome formulation that will not be pursued

¯

herein. In the latter case the volumetric variables p and J pertain only to an

element and can be eliminated at the element level. In such a situation the

¯

simplest discontinuous interpolation is to make p and J constant through-

out the element. This is the so-called mean dilatation technique discussed

in Section 6.6.5. Observe, however, that for simple constant stress elements

such as the linear triangle and tetrahedron the mean dilatation method co-

incides with the standard formulation and therefore su¬ers the detrimental

locking phenomenon.

7.5.1 IMPLEMENTATION OF THE MEAN DILATATION

METHOD

Recall from Chapter 6, Section 6.6.5, that the mean dilatation approach for

a given volume v leads to a constant pressure over the volume, as indicated

by Equations (6.50a“b). When this formulation is applied to each element

(e) in a ¬nite element mesh the pressure becomes constant over the element

volume. In particular, assuming for instance that the potential shown in

20 DISCRETIZATION AND SOLUTION

Equation (6.51) is used, the uniform element pressure is given as,

v (e) ’ V (e)

(e)

p =κ (7.54)

V (e)

where V (e) and v (e) are the initial and current element volumes.

The internal equivalent nodal forces for a typical element (e) are given

by Equation (7.15b) where now the Cauchy stress is evaluated from,

σ = σ + p(e) I (7.55)

and the deviatoric stress σ is evaluated using the appropriate constitutive

equation given by (5.51), or (5.103).

Continuing with the discretization, recall the modi¬ed linearized virtual

work Equation (6.60), which for an element (e) is,

(e)

σ : [( u)T

DδWint (φ, δv)[u] = δd : c : µ dv + δv] dv

v (e) v (e)

κv (e)

(e)

+ κv (·u)(·δv);

¯ κ = (e)

¯ (7.56)

V

ˆ ˆ

where the elasticity tensor is c = c +c p and c is the distortional component

that depends upon the material used, and c p is given by Equation (5.55b)

as,

c p = p(I — I ’ 2i ) (7.57)

The average divergences are now rede¬ned for an element (e) as,

n

1 1

·u = ·u dv = Na dv (7.7.58a,b)

ua ·

v (e) v (e)

v (e) v (e) a=1

n

1 1

·δv = ·δv dv = δv a · Na dv

(7.7.58c,d)

v (e) v (e)

(e) (e)

v v a=1

Discretization of the ¬rst two terms in Equation (7.56) is precisely as given

in the previous section, but the ¬nal dilatation term needs further atten-

tion. For element (e) the contribution to the linearized internal virtual work

related to the dilatation and associated, as before, with nodes a and b is,

(e)

DδWκ (φ, Na δv a )[Nb ub ]

κ

= (e) δv a · Na dv Nb dv

ub ·

V v (e) v (e)

κ

= δv a · Na dv — Nb dv ub

V (e) v (e) v (e)

(e)

= δv a · K κ,ab ub (7.59)

21

7.6 NEWTON“RAPHSON ITERATION AND SOLUTION PROCEDURE

where the dilatational tangent sti¬ness component is obtained in terms of

the average Cartesian derivatives of the shape functions* as,

1

(e)

K κ,ab = κv (e)

¯ Na — Nb ; Na = Na dv (7.60)

v (e) v (e)

The complete discretization of Equation (7.56) can now be written in terms

of the total element tangent substi¬ness matrix as,

(e) (e) (e) (e)

K ab = K c,ab + K σ,ab + K κ,ab (7.61)

(e)

into which the surface pressure component K p,ab (see Equation (7.51)), may,

if appropriate, be included. Assembly of the complete linearized virtual

work and hence the tangent matrix follows the procedure given in Equa-

tions (7.52“3).

7.6 NEWTON“RAPHSON ITERATION AND SOLUTION

PROCEDURE

7.6.1 NEWTON“RAPHSON SOLUTION ALGORITHM

In the previous sections it was shown that the equilibrium equation was dis-

cretized as δW (φ, δv) = δvT R, whereas the linearized virtual work term is

expressed in terms of the tangent matrix as DδW (φ, δv)[u] = δvT Ku. Con-

sequently, the Newton“Raphson equation δW (δv, φk )+DδW (φk , δv)[u] = 0

given in Equation (6.2) is expressed in a discretized form as,

δvT Ku = ’δvT R (7.62)

Because the nodal virtual velocities are arbitrary, a discretized Newton“

Raphson scheme is formulated as,

Ku = ’R(xk ); xk+1 = xk + u (7.63)

Although it is theoretically possible to achieve a direct solution for a

given load case, it is however more practical to consider the external load F

as being applied in a series of increments as,

l

F= ∆Fi (7.64)

i=1

where l is the total number of load increments. Clearly, the more increments

taken, the easier it becomes to ¬nd a converged solution for each individual

* The emergence of (v (e) )2 and the average Cartesian derivatives in Equation (7.60) is simply to

conform to expressions occurring in the literature.

22 DISCRETIZATION AND SOLUTION

BOX 7.1: Solution algorithm

r INPUT geometry, material properties and solution parameters

r INITIALIZE F = 0, x = X, R = 0

r LOOP over load increments

r FIND ∆F using (7.15c)

r SET F = F + ∆F

r SET R = R ’ ∆F

r DO WHILE ( R / F > tolerance )

FIND K using (7.51b)

r

SOLVE Ku = ’R

r

UPDATE x = x + u

r

FIND F (e) , b(e) and σ (e) using (7.5), (7.9d) and typically (5.29)

r

FIND T using (7.15b)

r

FIND R = T ’ F

r

r ENDDO

r ENDLOOP

load step. Observe that in the case of a hyperelastic material the ¬nal