FEM for elastic-plastic problems

FEM for elastic-plastic problems

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

With thanks to:

P. Mika, A. Winnicki, A. Wosatko

TNO DIANA http://www.tnodiana.com FEAP http://www.ce.berkeley.edu/feap

Lecture scope

Physical nonlinearity Plastic flow theory Computational plasticity

Simulation of plastic deformations 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.

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

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

microscopic level

crystal shear dislocation

lattice slip

Plastic flow theory [1,3]

Load-carrying capacity of a material is not infinite, during deformation irreversible strains occur

Notions of plasticity theory

I Yield function f (σ) = 0

- determines the limit of elastic response

I Plastic flow rule ˙^p= ˙λm

- determines the rate of plastic strain

˙λ - plastic multiplier m - direction of plastic flow

(usually associated with the yield function m^T= n^T= _∂σ^∂f)

I Plastic hardening f (σ − α, κ) = 0 kinematic (κ = 0) or isotropic (α = 0)

I Loading/unloading conditions:

f ¬ 0, ˙λ ­ 0, ˙λf = 0 (unloading is elastic)

Plastic flow theory

Response is history-dependent, constitutive relations written in rates Plastic flow when

f = 0 and ˙f = 0

(plastic consistency condition) Additive decomposition

˙ = ˙^e+ ˙^p Bijective mapping

˙

σ = D^e˙^e

Introduce flow rule

˙

σ = D^e( ˙ − ˙λm) Consistency

˙f = _∂σ^∂fσ +˙ _∂κ^∂f ˙κ

Hardening modulus h = −¹_˙

λ

∂f

∂κ˙κ Substitute ˙σ into n^Tσ − h ˙λ = 0˙

Determine plastic multiplier

˙λ =_h+n^nT_T^D_D^e_e^˙_m

Constitutive equation

˙ σ =h

D^e−_h+n^Demn_T^T_D^D_em^ei

˙

Tangent operator D^ep= D^e−_h+n^Demn_TD^TD_em^e Time integration necessary at the point level

Huber-Mises-Hencky plasticity

Most frequently used is the Huber-Mises-Hencky (HMH) plastic flow theory, based on a scalar measure of distortional energy J₂^σ

I Yield function

e.g. with isotropic hardening f (σ, κ) =p3J₂^σ− ¯σ(κ) = 0 κ - plastic strain measure ( ˙κ = ¹_¯σσ^T˙^p= ˙λ)

I Associated flow rule

˙^p= ˙λ_∂σ^∂f

I Hardening rule e.g. linear

¯

σ(κ) = σ_y + hκ h - hardening modulus

Response: force-displacement diagrams

Ideal plasticity Hardening plasticity

Plastic flow theory

Yield functions for metals:

Coulomb-Tresca-Guest i Huber-Mises-Hencky (HMH)

Insensitive to hydrostatic pressure p = ¹₃I₁^σ

Plastic flow theory

Yield functions for soil:

Mohra-Coulomb i Burzyński-Drucker-Prager (BDP)

Sensitive to hydrostatic pressure

‘Yield’ functions for concrete (plane stress)

Kupfer’s experiment

Rankine ‘yield’ function: f (σ, κ) = σ1− ¯σ(κ) = 0 Inelastic strain measure κ = |^p₁|

Computational plasticity

Return mapping algorithm

→ backward Euler algorithm (unconditionally stable) 1) Compute elastic predictor

σtr = σ^t+ D^e∆

2) Check f (σtr, κ^t) > 0 ? If not then elastic compute σ = σtr

If yes then plastic compute plastic corrector σ = σtr− ∆λD^em(σ) f (σ, κ) = 0

(set of 7 nonlinear equations for σ, ∆λ) Determine κ = κ^t+ ∆κ(∆λ)

σ σ

σ

f = 0

Iterative corrections still necessary unless radial return is performed.

Brazilian split test

Elasticity, plane strain

Deformation, vertical stress σyy and stress invariant J₂^σ

Brazilian split test

Elasticity, mesh sensitivity of stresses

Stress σyy for coarse and fine meshes

Stress under the force goes to infinity (results depend on mesh density) - solution at odds with physics

Brazilian split test

Ideal Huber-Mises-Hencky plasticity

Final deformation and stress σyy

Brazilian split test

Ideal Huber-Mises-Hencky plasticity

Final strain _yy and strain invariant J₂^

Brazilian split test

Ideal Huber-Mises-Hencky plasticity

0 0.2 0.4 0.6 0.8 1

Displacement 0

200 400 600 800

Force

0 0.2 0.4 0.6 0.8 1

Displacement 0

200 400 600 800

Force

This is correct!

For four-noded element load-displacement diagram exhibits artificial hardening due to so-called volumetric locking, since HMH flow theory contains kinematic constraint - isochoric plastic behaviour which cannot be reproduced by FEM model.

Eight-noded element does not involve locking.

Brazilian split test

Elasticity, plane strain, eight-noded elements

Deformation, vertical stress σ_yy and stress invariant J₂^σ

Brazilian split test

HMH plasticity

Final deformation and stress σ_yy

Brazilian split test

HMH plasticity

Final strain yy and invariant J₂^

Burzyński-Drucker-Prager plasticity

I Yield function with isotropic hardening f (σ, κ) = q + α p − βcp(κ) = 0 q =√

3J2- deviatoric stress measure p = ¹₃I1- hydrostatic pressure α = _{3−sin ϕ}^{6 sin ϕ} , β = _{3−sin ϕ}^{6 cos ϕ} ϕ - friction angle

cp(κ) - cohesion

I Plastic potential f^p= q + α p α = _{3−sin ψ}^{6 sin ψ}

ψ - dilatancy angle Nonassociated flow rule

˙^p= ˙λm, m = ^∂f_∂σ^p

I Plastic strain measure

˙κ = η ˙λ, η = (1 +²₉ α²)¹²

I Cohesion hardening modulus h(κ) = ηβ ^∂c_∂κ^p

q

p HMH

BDP

ϕ βc_p

Huber-Mises-Hencky yield function is retrieved for sin ϕ = sin ψ = 0

Slope stability simulation

Gradient-enhanced BDP plasticity

Evolution of plastic strain measure

Final remarks

1. In design one usually accepts calculation of stresses (internal forces) based on linear elasticity combined with limit state analysis

considering plasticity or cracking.

2. In nonlinear computations one estimates the load multiplier for which damage/failure/buckling of a structure occurs. The multiplier can be interpreted as a global safety coefficient, hence the computations should be based on medium values of loading and strength.