NONLINEAR ANALYSIS IN FEM

NONLINEAR ANALYSIS IN FEM

INTRODUCTION TO COMPUTATIONAL MECHANICS OF MATERIALS Civil Engineering, 1st cycle studies, 7th semester

elective subject academic year 2014/2015

Institute L-5, Faculty of Civil Engineering, Cracow University of Technology

Adam Wosatko

Jerzy Pamin

Lecture contents

FEM discretization

Nonlinear problems

Geometrical nonlinearity Physical nonlinearity Incremental-iterative analysis References:

R. de Borst and L.J. Sluys.

Computational methods in nonlinear solid mechanics.

Lecture Notes, Delft University of Technology, Delft, 1999.

DIANA Finite Element Analysis - User’s manual, release 7.2.

TNO Building and Construction Research, Delft, 1999.

Equilibrium equations for deformable continuum

Equilibrium equations + static boundary conditions:

L ^T σ + ˆ b = 0 w V , σν = ˆ t na S where:

L – differential operator matrix

σ – tensor/vector of generalized stresses b – body force vector

ν

V S ˆ t

Weak form of equilibrium equations:

Z

V

δu ^T (L ^T σ + ˆ b) dV = 0 ∀δu Virtual work principle δW int = δW ext :

Z

V

(Lδu) ^T σ dV = Z

V

δu ^T ˆ b dV + Z

S

δu ^T ˆ t dS

where

FEM discretization

Displacement-based finite elements:

u ^h =

n _w

X

i =1

N i (ξ, η, ζ)q i = Nq ^e

where: N - shape function matrix

q ^e - element vector of degrees of freedom (dofs) n _w - number of element nodes

Transformation of nodal degrees of freedom:

q ^e = I T ^e Q where: Q - global vector of dofs.

Weak form of equilibrium for discretized system

n e

X

e=1

I T ^{e T}

Z

V ^e

B ^T σ dV = F _ext , B = LN

Isoparametric approach and numerical integration are applied.

When is FE isoparametric?

Geometry approximation:

∀P ∈ Ω ^e : x P (ξ, η) = N(ξ, η) x ^e , where:

x P =

 x P

y P



, x ^e =























 x i

y i

x j

y j

x k

y _k x _l y _l























 FE is isoparametric if we use the same nodes and the same shape functions for geometry and displacement field approximation.

x(ξ, η) = N(ξ, η) x ^e u(ξ, η) = N(ξ, η) q ^e

Linear elasticity

Hooke’s law

Tensor notation: σ = D : , σ _ij = D _ijkl  _kl Matrix notation:

σ = D, σ =









σ _x σ _y σ _z τ _xy τ _yz τ _zx







 ,  =









 _x

 _y

 _z γ _xy γ _yz γ _zx









E 1 σ



Linear kinematic equations

Tensor notation:  = ¹ ₂ [∇u + (∇u) ^T ],  ij = ¹ ₂ (u i ,j + u j ,i ) Matrix notation:  = Lu

Hence stress tensor:

σ = D = DLu = DLNu = DB I T ^e Q Equilibrium equations for discretized system:

n _e

X

e=1

I T ^{e T}

Z

V ^e

B ^T DB dV I T ^e Q = F ext , K Q = F ext

Numerical integration for FE

Gauss quadrature (2D)

k ^e = Z Z

A ^e

B ^T DB h dA =

1

Z

−1 1

Z

−1

B(ξ, η) ^T D B(ξ, η) h det J dξdη

≈

n

X

i =1 m

X

j =1

w _i w _j B ^T _{(i ,j )} D B _{(i ,j )} h det J _{(i ,j )}

Integration Q4 Q8

full (FI)

reduced (RI)

Nonlinear problem

F ext applied in increments:

t → t + ∆t σ ^t+∆t = σ ^t + ∆σ

Equilibrium at time t + ∆t:

n e

X

e=1

I T ^{e T}

Z

V ^e

B ^T σ ^t+∆t dV = F ^t+∆t _ext

n _e

X

e=1

I T ^{e T}

Z

V ^e

B ^T ∆σ dV = F ^t+∆t _ext − F ^t _int where: F ^t _int = P n e

e=1 T I ^{e T} R

V ^e B ^T σ ^t dV Linearization of the left-hand side at time t:

∆σ = ∆σ(∆(∆u))

So far no assumptions about kinematics/constitutive properties!

Nonlinearity sources

Caused by change of geometry of (deformable) body:

I large strains (rubber, metal forming),

I large displacements (e.g. slender, thin-walled structures),

I contact (interaction of bodies in contact),

I follower load (varying with deformation) Caused by nonlinear constitutive relations:

I plasticity (irreversible strains)

I damage (degradation of elastic properties)

I fracture (continuous representation of cracks)

I . . .

Superposition principle does not hold.

Geometrical nonlinearity

X

u

x φ(X, t)

V V 0

S S ₀

x 1 , X 1

x 2 , X 2

Initial and current configuration Motion function: x = φ(X, t)

Displacement vector: u(X, t) = x − X

Deformation gradient (main deformation measure): F = ^∂φ _∂X = ∇ X x Strain tensor (one of possible strain measure):

E = 1

2 (F ^T F − I) = 1

2 [∇ _X u + (∇ _X u) ^T + (∇ _X u) ^T ∇ X u]

Geometrical nonlinearity

Nonlinear kinematic equations, e.g. ε _x = ε ^L _x + ε ^N _x = ^∂u _∂x + ¹ ₂ ^∂w _∂x
²

∆σ = ∆σ(∆(∆u))

I Balance equations describe the equilibrium of deformed body. The virtual work principle can be written for the initial or current configuration.

I Different stress measures are associated with different strain measures.

I Small strains: E ≈  = ¹ ₂ [∇u + (∇u) ^T ] < 2%.

I Small displacements (and rotations): V ≈ V ₀ (one description,

equilibrium equations for undeformed configuration).

Geometrical nonlinearity

Equilibrium of discretized system:

K ∆Q = F ^t+∆t _ext − F ^t _int where tangent stiffness matrix:

K = K 0 + K u + K σ

K 0 - linear stiffness matrix K u - initial displacement matrix

(discrete kinematic relations matrix B dependent on displacements)

K _σ - initial stress matrix (dependent on generalized stresses)

Physical nonlinearity

K ∆Q = F ^t+∆t _ext − F ^t _int Linearization of LHS at time t:

∆σ = ∆σ(∆(∆u))

∆σ = ^∂σ _∂
^t ∂

∂u

^t

∆u D = ^∂σ _∂ , L = ^∂ _∂u

Discretization: ∆u = N∆q ^e

Linear geometrical relations → Matrix of discrete kinematic relations

B = LN independent of displacements

Physical nonlinearity

K ∆Q = F ^t+∆t _ext − F ^t _int Tangent stiffness matrix

K =

n e

X

e=1

I T ^{e T}

Z

V ^e

B ^T D B dV A ^e

We limit interest to physical nonlinearity!

Iterative corrections necessary in order to reach equilibrium

at time t + ∆t → Newton-Raphson algorithm

Plastic state in a beam

A B C

displacement P

A

+

-

σ y

σ y

σ y

σ y

σ y

σ y

+

- -

+ C B

force

elastic material

plastified

elastic material

Incremental analysis

Explicit algorithm

u f

Real equilibrium path Numerical solution

Load increment ∆f

How can one improve the solution in subsequent increments?

Algorithm of incremental-iterative method

Iterative corrections necessary in order to reach equilibrium at time t + ∆t → Newton-Raphson algorithm

Out-of-balance (residual) forces: R j = f _ext ^t+∆t − f _int,j ^t+∆t → 0 f

R 1

f _int,1 f _ext ^t

f _ext ^t+∆t

∆f _ext

K ∆Q = f _ext ^t+∆t − f _int ^t K - tangent operator First iteration:

∆Q 1 = K ⁻¹ ₀ (f _ext ^t+∆t − f int,0 ) σ 1 → f int,1 6= f _ext ^t+∆t Corrections:

dQ j = K ⁻¹ _{j −1} (f _ext ^t+∆t − f _{int,j −1} ^t+∆t ) σ j → f int,j

Convergence criterion:

kf ^t+∆t _ext −f ^t+∆t

int,j k

k∆f _ext k ≤ δ

Algorithm of incremental-iterative method

Iterative corrections necessary in order to reach equilibrium at time t + ∆t → Newton-Raphson algorithm

Out-of-balance (residual) forces: R j = f _ext ^t+∆t − f _int,j ^t+∆t → 0 f

R 1

f _int,1

u u ^t+∆t

u ^t f _ext ^t f _ext ^t+∆t

∆u ₁ du ₂

∆f _ext

K ∆Q = f _ext ^t+∆t − f _int ^t K - tangent operator First iteration:

∆Q 1 = K ⁻¹ ₀ (f _ext ^t+∆t − f int,0 ) σ 1 → f int,1 6= f _ext ^t+∆t Corrections:

dQ j = K ⁻¹ _{j −1} (f _ext ^t+∆t − f _{int,j −1} ^t+∆t ) σ j → f int,j

Convergence criterion:

kf ^t+∆t _ext −f ^t+∆t

int,j k

k∆f _ext k ≤ δ

Algorithm of incremental-iterative method

Iterative corrections necessary in order to reach equilibrium at time t + ∆t → Newton-Raphson algorithm

Out-of-balance (residual) forces: R j = f _ext ^t+∆t − f _int,j ^t+∆t → 0 f

R 1

f _int,1 f _ext ^t

f _ext ^t+∆t

∆f _ext

K ∆Q = f _ext ^t+∆t − f _int ^t K - tangent operator First iteration:

∆Q 1 = K ⁻¹ ₀ (f _ext ^t+∆t − f int,0 ) σ 1 → f int,1 6= f _ext ^t+∆t Corrections:

dQ j = K ⁻¹ _{j −1} (f _ext ^t+∆t − f _{int,j −1} ^t+∆t ) σ j → f int,j

Convergence criterion:

kf ^t+∆t _ext −f ^t+∆t

int,j k

k∆f _ext k ≤ δ

Algorithm of incremental-iterative method

Iterative corrections necessary in order to reach equilibrium at time t + ∆t → Newton-Raphson algorithm

Out-of-balance (residual) forces: R j = f _ext ^t+∆t − f _int,j ^t+∆t → 0 f

u u ^t+∆t

u ^t f _ext ^t f _ext ^t+∆t

∆u ₁ du ₂

∆f _ext

R ₂

f int,2

K ∆Q = f _ext ^t+∆t − f _int ^t K - tangent operator First iteration:

∆Q 1 = K ⁻¹ ₀ (f _ext ^t+∆t − f int,0 ) σ 1 → f int,1 6= f _ext ^t+∆t Corrections:

dQ j = K ⁻¹ _{j −1} (f _ext ^t+∆t − f _{int,j −1} ^t+∆t ) σ j → f int,j

Convergence criterion:

kf ^t+∆t _ext −f ^t+∆t int,j k k∆f _ext k ≤ δ Modified algorithm:

K j = K 0

Newton-Raphson algorithm

Possible iteration schemes

u f

u f

Standard algorithm Modified algorithm

K j K j = K 0

Update in each iteration Update at the start of increment

Options of incremental control

Force or displacement control

Arc length control