FEM - solution of nonlinear problems

FEM - solution of nonlinear problems

Jerzy Pamin e-mail: JPamin@L5.pk.edu.pl

With thanks to:

A. Wosatko, A. Winnicki

ADINA R&D, Inc.http://www.adina.com ANSYS, Inc. http://www.ansys.com TNO DIANA http://www.tnodiana.com

Lecture scope

Nonlinear problems Geometrical nonlinearity Physical nonlinearity Cracking

Final remarks References

[1] R. de Borst and L.J. Sluys. Computational Methods in Nonlinear Solid Mechanics. Lecture notes, Delft University of Technology, 1999.

[2] G. Rakowski, Z. Kacprzyk. Metoda elementow skończonych w mechanice kostrukcji. Oficyna Wyd. PW, Warszawa, 2005.

[3] M. Jir´asek and Z.P. Baˇzant. Inelastic Analysis of Structures. J. Wiley &

Sons, Chichester, 2002.

[4] M. Kwasek Advanced static analysis and design of reinforced concrete deep beams. Diploma work, Politechnika Krakowska, 2004.

Nonlinearity sources

Caused by change of geometry of (deformable) body

I large strains (e.g. 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 nonliear constitutive relations

I plasticity (irreversible strains)

I damage (degradation of elastic properties)

I fracture (continuous representation of cracks)

I . . . Remarks:

I Superposition principle does not hold.

I It is possible to describe discontinuum in which components are connected by interfaces (e.g. composite structures) or (discrete) cracks occur. The interfaces usually have nonlinear features which represent for instance friction, adhesion, cracking.

Nonlinear continuum [1,2,3]

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 δWint = δWext

Z

V

(Lδu)^Tσ dV = Z

V

δu^Tb dV + Z

S

δu^Tˆt dS

Galerkin method

Displacement-based finite elements

u ≈ u^h=

nw

X

i =1

N_i(ξ, η, ζ)u_i= Nd^e

where: N - shape function matrix, d^e - element vector of degrees of freedom (dofs),

n_w - number of nodes

Transformation of nodal degrees of freedom d^e = A^ed where: d - global vector of dofs

Weak form of equilibrium for discretized system

n_e

X

e=1

A^{e T} Z

V^e

B^Tσ dV = f_ext, B = LN

Isoparametric approach, numerical integration

Linear elasticity

Hooke’s law

Tensor notation: σ = D^e : , σij = D_ijkl^e kl

Matrix notation:

σ = D^e, σ =









σx

σy

σz

τxy

τyz

τzx







 ,  =









x

y

z

γxy

γyz

γzx







 Material isotropy: D^e= D^e(E , ν)

E 1

σ



loading unloading

Linear kinematic equations

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

Hence stress tensor: σ = D^e = D^eLu = D^eLNd^e = D^eBA^ed Equilibrium equations for discretized system

ne

X

e=1

A^{e T} Z

V^e

B^TD^eB dV A^ed = f_ext, Kd = f_ext

Incremental-iterative analysis

Nonlinear problem:

fext applied in increments

t → t + ∆t → σ^t+∆t= σ^t+ ∆σ Equilibrium at time t + ∆t:

ne

X

e=1

A^{e T} Z

V^e

B^Tσ^t+∆tdV = f_ext^t+∆t

n_e

X

e=1

A^{e T} Z

V^e

B^T∆σ dV = f_ext^t+∆t− f_int^t where: f_int^t =Pne

e=1A^{e T}R

V^eB^Tσ^tdV Linearization of the left-hand side at time t

∆σ = ∆σ(∆(∆u)) Equation set for an increment:

K ∆d = f_ext^t+∆t− f_int^t

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: Rj = f_ext^t+∆t− f_int,j^t+∆t→ 0 f

R1

f_int,1

u u^t+∆t

u^t f_ext^t f_ext^t+∆t

∆u₁ du₂

∆f_ext

K ∆d = f_ext^t+∆t− fint^t

K- tangent operator First iteration:

∆d1=K⁻¹₀ (f_ext^t+∆t− fint^t ) σ1 → f_int,1^t+∆t6= f_ext^t+∆t Corrections:

ddj +1= K⁻¹_j (f_ext^t+∆t−f_int,j^t+∆t) σj +1 → f_{int,j +1}^t+∆t

Convergence criterion:

kf^t+∆t_ext −f^t+∆t int,j k k∆f_extk ¬ δ Modified algorithm:

Kj= K0

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: Rj = f_ext^t+∆t− f_int,j^t+∆t→ 0 f

R1

f_int,1

u u^t+∆t

u^t f_ext^t f_ext^t+∆t

∆u₁ du₂

∆f_ext

K ∆d = f_ext^t+∆t− fint^t

K - tangent operator First iteration:

∆d1= K⁻¹₀ (f_ext^t+∆t− f_int,0^t+∆t) σ1 → f_int,1^t+∆t6= f_ext^t+∆t Corrections:

ddj +1= K⁻¹_j (f_ext^t+∆t−f_int,j^t+∆t) σj +1 → f_{int,j +1}^t+∆t

Convergence criterion:

kf^t+∆t_ext −f^t+∆t int,j k k∆f_extk ¬ δ Modified algorithm:

Kj= K0

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: Rj = f_ext^t+∆t− f_int,j^t+∆t→ 0 f

R1

f_int,1

u u^t+∆t

u^t f_ext^t f_ext^t+∆t

∆u₁ du₂

∆f_ext

Kjdd = f_ext^t+∆t− f_int,j^t+∆t K - tangent operator First iteration:

∆d1= K⁻¹₀ (f_ext^t+∆t− f_int,0^t+∆t) σ1 → f_int,1^t+∆t6= f_ext^t+∆t Corrections:

ddj +1= K⁻¹_j (f_ext^t+∆t−f_int,j^t+∆t) σj +1 → f_{int,j +1}^t+∆t

Convergence criterion:

kf^t+∆t_ext −f^t+∆t int,j k k∆f_extk ¬ δ Modified algorithm:

Kj= K0

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: Rj = 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

Kjdd = f_ext^t+∆t− f_int,j^t+∆t K - tangent operator First iteration:

∆d1= K⁻¹₀ (f_ext^t+∆t− f_int,0^t+∆t) σ1 → f_int,1^t+∆t6= f_ext^t+∆t Corrections:

ddj +1= K⁻¹_j (f_ext^t+∆t−f_int,j^t+∆t) σj +1 → f_{int,j +1}^t+∆t

Convergence criterion:

kf^t+∆t_ext −f^t+∆t int,j k k∆f_extk ¬ δ Modified algorithm:

Kj= K0

Options of incremental loading control

Force or displacement control

Arc length control

Geometrical nonlinearity

X

u

x φ(X, t)

V V₀

S S0

x1, X1

x2, X2

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

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

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

E = 1

2(F^TF − I) = 1

2[∇_Xu + (∇_Xu)^T+ (∇_Xu)^T∇_Xu]

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 [4]:

K ∆d = f_ext^t+∆t− f_int^t

where tangent stiffness matrix:

K = K0+Ku+Kσ

K0- linear stiffness matrix Ku - initial displacement matrix

(discrete kinematic relations matrix B dependent on displacements) K_σ - initial stress matrix (dependent on generalized stresses)

Physical nonlinearity

K ∆d = f_ext^t+∆t− f_int^t Linearization of LHS at time t:

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

∆σ = ^∂σ_∂
^t _∂

∂u

^t

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

Discretization: ∆u = N∆d^e

Linear geometrical relations → Matrix of discrete kinematic relations B = LN independent of displacements

Tangent stiffness matrix

K =

n_e

X

e=1

A^{e T} Z

V^e

B^TDB dV A^e

Plastic yielding of material

A B C

displacement force

P

A

+

-

σy

σy

σy

σy

σy

σy

+

- -

+ C B

elastic material

equivalent plastic strain distribution plastified

elastic material

Discrete and smeared cracks

Fracture energy Gf (used to create unit surface area of a crack)

Simulation of cracking in RC panel

with ATENA [4]

Final remarks

1. Computer simulations offer priceless possibilities, but should be performed only by conscious FEM users.

2. Three-dimensional (3D) modelling starts to dominate in prediction of nonlinear behaviour of structures.

3. Consistent linearization of equations guarantees quadratic convergence of Newton-Raphson algorithm.

4. In order to improve the quality of FEM approximation adaptive mesh refinenment based on evaluation of discretization error is advisable.