4.21 Nonlinear problems

Back

In the triple {u,E,S} we let now E the Green-Lagrangian strain tensor

and S the second Piola-Kirchhoff stress tensor, and we assume the material

402 4 Plane problems

to be hyperelastic, i.e., there exists a strain-energy function W such that

S = ∂W/∂E. Given volume forces p the elastic state S = {u,E,S} satisfies

at every point x of the undeformed body the system

E(u) − E = 0

1

2

(ui,j +uj ,i +uk,i uk,j ) − εij = 0

W

_(E) − S = 0 ∂W

εij

σij = 0 (4.166)

−div(S + ∇uS) = p −(σij + ui,k σkj),j = pi

and satisfies displacement boundary conditions u = u on a part ΓD of the

boundary and stress boundary conditions t(S,u) =t on the complementary

part ΓN where

t(S,u) := (S + ∇uS)n (4.167)

is the traction vector at a boundary point with outward normal vector n.

With symmetric stress tensors S we have the identity

_

Ω

−div(S + ∇uS) •ˆu dΩ

= −

_

Γ

t(S,u) •ˆu ds +

_

Ω

Euu) •S dΩ (4.168)

where

Euu) :=

1

2

(∇ˆu + ∇ˆuT + ∇uT ∇ˆu + ∇ˆuT ∇u) (4.169)

is the Gateaux derivative of the matrix E(u)

d

dε

[E(u + εˆu)]|ε=0 = Euu) . (4.170)

Collecting terms we can formulate Green’s first identity of the operator A(S),

that is the system (4.166)

G(S, ˆ S) = _A(S), ˆ S_ +

_

Γ

t(S,u) •ˆu ds

_ _ _

δWe

− a(S, ˆ S)

_ _ _

δWi

= 0 (4.171)

where _A(S), ˆ S_ is similar to (4.130) and where

a(S, ˆ S) =

_

Ω

(E(u) − E) •ˆS dΩ

+

_

Ω

(W

_(E) − S) • ˆE dΩ +

_

Ω

Euu) •S dΩ . (4.172)

4.21 Nonlinear problems 403

The identity (4.171) is the basis of many variational principles in nonlinear

mechanics and can be formulated in the same way also for beams and slabs,

[115].

In the case of a pure displacement formulation S = {u,E(u),W

_(E(u))}

and it is ˆu = 0 on ΓD, so that (4.171) reduces to

G(uu) =

_

Ω

p •ˆu dΩ +

_

ΓN

t

•ˆu ds −

_

Ω

Euu) •S dΩ = 0, (4.173)

where S = W

_(E(u)).

Next let uh =

_n

j uj ϕj(x) the FE solution and let ˆu = ϕi a virtual

displacement then

_

Ω

Euh(ϕi) •W

_(E(uh)) dΩ

_ _ _

ki

=

_

Ω

p ϕi dΩ +

_

ΓN

t

ϕi ds

_ _ _

fi

(4.174)

or

k(u) = f (4.175)

where u is the vector of nodal coordinates.

Linearization

For computational purposes a linearization of (4.175) is necessary. Let pΔ and

t

Δ be load increments, and let u+uΔ be the displacement field corresponding

to p + pΔ andt +tΔ then

G(u + uΔu) =

_

Ω

(p + pΔ) •ˆu dΩ +

_

ΓN

(t +tΔ) •ˆu ds

−a(u + uΔu) = 0, (4.176)

where

a(u + uΔu) :=

_

Ω

Eu+uΔu) •W

_(E(u + uΔ)) dΩ . (4.177)

The Gateaux derivative of the strain energy product

a(uu) :=

_

Ω

Euu) •W

_(E(u)) dΩ (4.178)

with respect to a displacement increment uΔ is

aT (u, uΔu) :=

_

d

dε

a(u + εuΔu)

_

ε=0

=

_

Ω

[∇uΔW

_(E(u)) •∇ˆu + Euu) •C[Eu(uΔ)]] dΩ , (4.179)

404 4 Plane problems

Fig. 4.64. Truss element

where the tensor

C = ∂2W

E E

= ∂

E

S (4.180)

is the second derivative of W, which is to be evaluated at u. Note that

aT (u,uΔu) is linear in the second and third argument, uΔ and ˆu.

We then let

a(u + uΔu) _ a(uu) + aT (u,uΔu) , (4.181)

so that (4.175) becomes

KT (u)uΔ = f k(u) , (4.182)

where now uΔ is the vector of nodal displacements of the field uΔ and KT

is the tangential stiffness matrix:

(KT )ij = aT (u,ϕj ,ϕi) . (4.183)

A truss element

We consider a truss element, that extends along the x-axis; see Fig. 4.64. The

deformation of the element is

ϕ(x, y, z) = (x + u(x)) e1 + (y + v(x)) e2 + (z + w(x)) e3 , (4.184)

where u(x) = [u(x), v(x), w(x)]T is the displacement vector. By definition the

deformation gradient is

4.21 Nonlinear problems 405

F = ∇ϕ =

1 + u_ 0 0

v_ 1 0

w_ 0 1

⎦ = I + ∇u, (4.185)

and the Green-Lagrangian strain tensor

E(u) =

1

2

(FT F I) =

⎢⎢⎣

εx

1

2 v_ 1

2 w_

1

2 v_ 0 0

1

2 w_ 0 0

⎥⎥⎦

, (4.186)

where

ε(u) := εx(u) =

1

2

((1 + u

_)2 + (v

_)2 + (w

_)2 − 1)

= u

_ +

1

2

((u

_)2 + (v

_)2 + (w

_)2) ()_ = d

dx

. (4.187)

Green’s first identity reads, if all functions are functions of x only,

G(uu) =

_ l

0

−div (S + ∇uS) •ˆu Adx + [(S + ∇uS) e1 •ˆu A]l

0

_ l

0

Euu) •S Adx = 0. (4.188)

We then let simply σ = σxx = E ε, where E is Young’s modulus, and we let

all other σij = 0, so that with N = A(σ + u_ σ)

G(uu) =

_ l

0

−N

_ ˆudx + [N ˆu ]l

0

_ l

0

εuu)σAdx

_ _ _

a(uu)

= 0, (4.189)

where

εuu) = (1+u

_) ˆu

_ + v

_ ˆv

_ + w

_ ˆ w

_

. (4.190)

The Gateaux derivative of the strain energy product a(uu) is

aT (u,uΔu) :=

_

d

dε

a(u + εuΔu)

_

ε=0

=

_ l

0

[εuΔu) σ + εuu) σ]Adx , (4.191)

where

εuΔu) =

_

d

dε

εu+εuΔu)

_

ε=0

= u

_

Δˆu

_ + v

_

Δ ˆv

_ + w

_

Δ ˆ w

_ (4.192)

406 4 Plane problems

and

σ =

_

d

dε

σ(u + εuΔ)

_

ε=0

= E εu(uΔ) . (4.193)

The elements of the tangential stiffness matrix are therefore

(KT (u))ij = aT (u,ϕi,ϕj) , (4.194)

where the vector-valued functions ϕi(ξ) are the nodal unit displacements

(three at each node).

With linear shape functions on an element Ωe with length le

ue(x) =

_2

i=1

uei

ϕei

(ξ), ve(x) =

_2

i=1

ve

i ϕei

(ξ), we(x) =

_2

i=1

we

i ϕei

(ξ) , (4.195)

we obtain for the nodal vector ue = [u1, v1, w1, u2, v2, w2]T the 6×6 tangential

stiffness matrix [255]

Ke

T (u) =

_

(A1 + A2) −(A1 + A2)

−(A1 + A2) (A1 + A2)

_

, (4.196)

where

A1 = EA

le

(1 + u_

e)2 (1 + u_

e) v_

e (1 + u_

e)w_

e

(1 + u_

e) v_

e (v_

e)2 v_

e w_

e

(1 + u_

e)w_

e v_

e w_

e (w_

e)2

⎦ (4.197)

and

A2 = σA

le

1 0 0

0 1 0

0 0 1

⎦ , (4.198)

so that after the assemblage

KT (u)uΔ = f k(u) . (4.199)

For an extension of these concepts to 3-D beam problems see [104].

Plane problem

Let X denote the initial coordinates of a point and x the coordinates of the

current configuration (see Fig. 4.65)

x = X + u, or xi = Xi + ui , (4.200)

where u is the displacement vector.

The deformation gradient is

4.21 Nonlinear problems 407

Fig. 4.65. Deformations of a

plane element

F = ∂x

X

= I + ∇u =

1 + u1,1 u1,2 0

u2,1 1 + u2,2 0

0 0 1

⎦ , (4.201)

so that the Green-Lagrangian strain tensor

E =

1

2

(FTF I) =

_

E11 E12

sym. E22

_

(4.202)

has the components

E11 = u1,1 +

1

2

(u21

,1 + u22

,1)

E22 = u2,2 +

1

2

(u22

,2 + u21

,2) (4.203)

E12 =

1

2

(u1,2 + u2,1) +

1

2

(u1,1 u1,2 + u2,2 u2,1) .

The FE displacement field can be written in two ways

uh =

D_OFS

i=1

ui ϕi(x) ≡

_

i

ui

_

.

.

_

, uh =

N_ODES

i=1

ui ψi(x) ≡

_

i

_

.

.

_

ψi ,

(4.204)

where the sum extends either over the degrees of freedom ui or the nodes

of the structure. The ϕi(x) are the nodal unit displacement fields associated

with the ui, see (4.27) on p. 337, while the ψi in the second formula are the

scalar-valued shape functions of the nodes (ψi(xj) = δij) and the vectors

ui = [u(i)

1 , u(i)

2 ]T are the nodal displacements at node i. Here, we use the

second notation so that, for example, the virtual exterior work of the volume

forces p becomes

408 4 Plane problems

δWe =

_

Ω

p uh dΩ =

_n

i=1

_

Ω

p ui ψi dΩ =

_n

i=1

uTi

_

Ω

p ψi dΩ =

_n

i=1

uTi

fi .

(4.205)

The gradient of the FE displacement field is then

uh =

_n

i=1

ui ψi (4.206)

and the deformation gradient,

F = I + ∇uh = I +

_n

i=1

ui ψi . (4.207)

The gradient of the virtual displacement field ˆuh =

_

i ˆui ψi follows (4.206)

so that Gateaux derivative becomes (we drop the subscript h on uh and ˆuh

for a moment)

Euu) : =

1

2

(∇ˆu + ∇ˆuT + ∇uT ∇ˆu + ∇ˆuT ∇u)

=

1

2

(FT ∇ˆu + ∇ˆuTF)

=

1

2

_n

i=1

[FT (ˆui ψi) + (∇ψi ˆui)F] .

(4.208)

As in the linear theory where E S = ε σ the symmetry of Euu) and S is

motivation to introduce a “Gateaux strain vector”

εuu) :=

(Euu))11

(Euu))22

2 (Euu))12

⎦ =

_n

i=1

BLi ˆui , (4.209)

where

BLi =

F11 ψi,1 F21 ψi,1

F12 ψi,2 F22 ψi,2

F11 ψi,2 + F12 ψi,1 F21 ψi,2 + F22 ψi,1

⎦ . (4.210)

Given a St. Venant type material

S = C[E] = 2μE + λ(trE) I (4.211)

the stress vector of the second Piola-Kirchhoff stress tensor is:

σ =

S11

S22

S12

⎦ =

λ + 2μ λ 0

λ λ+ 2μ 0

0 0 μ

E11

E22

2E12

⎦ , (4.212)

4.21 Nonlinear problems 409

where

μ = E

2(1 + ν), λ= Eν

(1 + ν)(1 − 2ν) . (4.213)

Hence the weak form of the equilibrium conditions is

G(uh,ˆuh) =

_n

i=1

ˆu

T

i [

_

Ω

BT

Li

σ dΩ

_ _ _

ki

_

Ω

p ψi dΩ

_

ΓN

t

ψi ds

_ _ _

fi

] = 0

(4.214)

or

k(u) − f = 0 . (4.215)

As on p. 403 we linearize this equation with the help of a Gateaux derivative

of the strain energy product in the direction of uΔ

aT (u,uΔu) =

_

Ω

[∇uΔ S •∇ˆu + Euu) •C[Eu(uΔ)]] dΩ . (4.216)

The discretization of the first term in (4.216) yields with

uΔh =

_n

j=1

uΔj ψj , ∇ˆuh =

_n

i=1

ˆu

i ψi (4.217)

the so-called initial stress stiffness matrix

_

Ω

uΔ S •∇ˆu dΩ =

_n

i=1

_n

j=1

ˆu

T

i

_

Ω

Gij I dΩ uΔj (4.218)

where

Gij := ∇Tψi S ψj =

_

ψi,1 ψi,2

_ _

S11 S12

S21 S22

__

ψj,1

ψj,2

_

. (4.219)

With

Euh(uΔ) =

1

2

_n

j=1

[FT (uΔj ψj) + (∇ψj uΔj)F] =

_n

j=1

BLjuΔj

(4.220)

the second term in (4.216) becomes:

_

Ω

Euu) •C[Eu(uΔ)] dΩ =

_n

i=1

_n

j=1

ˆu

T

i

_

Ω

BT

LiCBLj dΩ uΔj . (4.221)

410 4 Plane problems

Fig. 4.66. Nonlinear analysis of a cantilever plate that carries an equivalent nodal

force; large displacements, small strains

Adding all terms we obtain

aT (u,uΔu) =

_n

i=1

_n

j=1

ˆu

T

i KTij uΔj (4.222)

where

KTij :=

_

Ω

_

Gij I + BT

LiCBLj

_

dΩ . (4.223)

In the case of a bilinear element the (2×2) sub matrices KTij of the tangential

stiffness matrix of a single element are arranged as follows:

Ke

T =

⎢⎢⎣

KT11 KT12 KT13 KT14

KT21 KT22 KT23 KT24

KT31 KT32 KT33 KT34

KT41 KT42 KT43 KT44

⎥⎥⎦

. (4.224)

After the assemblage, the resulting system of equations

KT (u)uΔ = f k(u) , (4.225)

is solved with the Newton-Raphson method.

A cantilever plate

A cantilever plate of length l = 10 m, width h = 1 m and thickness t = 1 m,

is loaded at its end with a vertical point force P = 50, 000 kN; see Fig. 4.66

and 4.68. The parameter of the St. Venant type material are E = 107 kN/m2

and ν = 0. The analysis was done with plane bilinear elements. The vertical

4.21 Nonlinear problems 411

Table 4.14. Deflection w (m) at the foot of the point load

element linear nonlinear

mixed shell elements 19.950 7.533

bilinear plate elements 13.413 7.307

0 2 4 6 8 10 12 14

1

2

3

4

5

6

Fig. 4.67. Deflection

point load

2 4 6 8 10

1

2

3

4

5

x 10 kN

4

deflection(m)

Fig. 4.68. Maxwells

true

deflection w at the base of the point load is plotted in Fig. 4.67. For small

values of P the linear and the nonlinear results coincide. But if the load P

increases the increase in the deflection slows down, i.e., in linear analysis the

deflections are overestimated! This was confirmed in independent tests with

ADINA. The reference solution in Table 4.66 is based on a mixed four-node

shell element [45]. The nonlinear results agree quite well, while in the linear

case—one is tempted to say: as usual—the deviations are larger. Obviously the

bilinear element fares better in nonlinear problems than in linear problems!

at the foot of the

theorem is no longer

412 4 Plane problems

Goal-oriented refinement

In nonlinear problems the dual problem for the generalized Green’s function

z is formulated at the equilibrium point u, see Sect. 7.5, p. 526,

z ∈V a

T (u; z, v) = J

_(u; v) ∀ v ∈ V . (4.226)

and a

T corresponds to the tangential form aT of the original problem at the

equilibrium point, see (7.181) p. 528, so that zh is the solution of

KT (uh) z = j (4.227)

where KT (uh) is the tangential stiffness matrix at uh and the components

jk = J

_(uh;ϕk) . (4.228)

of j are the equivalent nodal forces. While in linear problems the equivalent

nodal forces jk are simply the values J(ϕk) of the different shape functions we

now must evaluate the Gateaux derivative (with respect to u) of the functional

J(.) for each single ϕk.

Let us assume we calculate the stress σij(x) at a point x, that is

J(u) = σij(u)(x) . (4.229)

The Gateaux derivative of this functional is

J

_(u; v) :=

_

d

dε

J(u + ε v)

_

ε=0

=

_

d

dε

σij(u + ε v)

_

ε=0

(4.230)

or with S = C[E(u)],

_

d

dε

S(u + ε v)

_

ε=0

= C

_

d

dε

E(u + ε v)

_

ε=0

= C[Eu(v)] =: ˆ S

(4.231)

where Eu(v) is the Gateaux derivative of the Green-Lagrangian strain tensor,

see (4.208). Hence the equivalent nodal forces

jk = J

_(uh;ϕk) = ˆσij(ϕk)(x) (4.232)

are the component ˆσij of the “tangential” stress tensor C[Euh(ϕk)]. These

are the stress increments at the actual equilibrium point uh resulting from

the displacement field ϕk.

For a second example we consider the integral of the stresses along a cross

section A-A

J(u) =

_

A−A

σij(u) ds . (4.233)

4.21 Nonlinear problems 413

The Gateaux derivative of this functional is

J

_(u; v) =

_

A−A

ˆσij(v) ds (4.234)

so that the equivalent nodal forces are the integrals of the stresses ˆσij resulting

from the unit displacement fields ϕk

jk = J

_(uh;ϕk) =

_

A−A

ˆσij(ϕk) ds (4.235)

The aim of goal oriented methods is to improve the accuracy of the output

functional J(u) by minimizing the energy error in the associated generalized

Green’s functional and in the solution u itself.

Now how do we proceed? We apply the first load increment p1 (let p =

p1+p2+. . .) and we find the equilibrium point u1 and the generalized Green’s

function z1. Then the mesh is adaptively refined so that both errors (in u1

as well as in z1) are below a certain threshold value. This completes the cycle

and we apply the next increment p2 and repeat the process, etc. At the end

we have a mesh which is optimal for the output value J(u).

As in the linear case the final values of the stresses σij(x) are computed

directly by differentiating the FE solution. Actually, we have no other choice.

Unlike linear problems where

J(u) =

_

Ω

p z dΩ (4.236)

no such formula exists in nonlinear problems. And also such formulations as

J(u) _

_

Ω

p1

z1 dΩ +

_

Ω

p2

z2 dΩ +

_

Ω

p3

z3 dΩ + . . . (4.237)

or

J(u) _

_

Ω

p zn dΩ zn = final Green’s function (4.238)

lead to nowhere.

The results in Fig. 4.69 illustrate nicely how the actions of the Green’s

functions automatically (!) follow the movement of the structure—a cantilever

plate to which a nodal force is applied. Depicted is the orientation of the nodal

forces of the FE Green’s functions for the two functionals

J(u) = σij(u)(x) J(u) =

_

A−A

σij(u) ds (4.239)

at the final stage. For not to complicate the drawings the meshes were not

refined. Just a normal nonlinear iterative analysis was performed for each of

the two solutions u and z.

414 4 Plane problems

Fig. 4.69. Nonlinear analysis of a cantilever plate a) original load b) deformed

structure c) equivalent nodal forces of the dual problem for σyy(x) at the equilibrium

point d) FE approximation of the dual solution for σyy(x) e) equivalent nodal forces

for

2

2 A σxy ds in cross-section A-A f) FE approximation of the dual problem for

A σxy ds, [162]

5. Slabs

The bending theory of plates can be viewed as an extension of beam theory.

In an Euler–Bernoulli-beam, EIwIV = p, shear deformations are neglected,

i.e., a straight line that is initially normal to the neutral axis remains so after

the load is applied. In a Timoshenko beam, the line instead rotates by an

angle γ; see Fig. 5.1 b.

The extension of the Euler–Bernoulli-beam to plate theory is the Kirchhoff

plate, KΔΔw = p, and the extension of the Timoshenko beam is the

Reissner–Mindlin plate. The first finite elements developed for plate bending

problems were based on the Kirchhoff plate theory. But the problem that the

shape functions must be C1 and must be easy to be implement at the same

time soon favored Reissner–Mindlin plate elements, where C0 suffices for the

shape functions, see Fig. 5.2.

Normally slabs are relatively thin with negligible shear deformations, so

that good Reissner–Mindlin plate elements tend to produce the same results

Fig. 5.1. Cantilever beams: a) EulerBernoulli beam b) Timoshenko beam

Fig. 5.2. A Kirchhoff plate cannot be

folded like sheet metal

416 5 Slabs

Fig. 5.3. Resultant stresses in a slab

as Kirchhoff plate elements, with the exception possible of zones close to the

boundary. The switch to the Reissner–Mindlin plate theory in FE programming

therefore probably remained unnoticed in the engineering community,

even though the most popular plate elements, like the DKT element and the

Bathe–Dvorkin element are not genuine Reissner–Mindlin plate elements but

rather an ingenious mixture of Kirchhoff and Reissner–Mindlin plate theories.