FEM for elastic-plastic problems CMCE Lecture 4, Civil Engineering, II cycle, specialty BEC

FEM for elastic-plastic problems

CMCE Lecture 4, Civil Engineering, II cycle, specialty BEC

Jerzy Pamin

Institute for Computational Civil Engineering

Civil Engineering Department, Cracow University of Technology 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

Comp.Meth.Civ.Eng., II cycle

Incremental-iterative analysis

Nonlinear problem:

fext applied in increments

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

n_e

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 =Pn_e

e=1A^{e T}R

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

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

Physical nonlinearity

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

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

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

∂u

^t

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

Discretization: ∆u = N∆u^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

Comp.Meth.Civ.Eng., II cycle

Plastic yielding of material

A B C

displacement force

P

A

+

-

σy

σ_y σy

σ_y σy

σ_y +

- -

+ C B

elastic material plastified

elastic material

microscopic level

crystal shear dislocation

lattice slip

Plastic flow theory [1,2]

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)

Comp.Meth.Civ.Eng., II cycle

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_TD^TD_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

Comp.Meth.Civ.Eng., II cycle

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₁^σ

Comp.Meth.Civ.Eng., II cycle

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 (σ, κ) = σ₁ − ¯σ(κ) = 0 Inelastic strain measure κ = |^p₁|

Comp.Meth.Civ.Eng., II cycle

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₂^σ

Comp.Meth.Civ.Eng., II cycle

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

Comp.Meth.Civ.Eng., II cycle

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.

Comp.Meth.Civ.Eng., II cycle

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

Comp.Meth.Civ.Eng., II cycle

Brazilian split test

HMH plasticity

Final strain _yy and invariant J₂^

Tension of perforated plate

Deformation and invariant J₂^

Ideal plasticity Plasticity with hardening

Comp.Meth.Civ.Eng., II cycle

Tension of perforated plate

Force-displacement diagram

0 200 400 600 800 1000 1200 1400 1600

0 0.5 1 1.5 2

Displacement [mm]

plasticity with hardening

Force[N]

ideal plasticity

Burzyński-Drucker-Prager plasticity

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

q = √

3J₂ - deviatoric stress measure p = ¹₃I₁ - hydrostatic pressure

α = _{3−sin ϕ}^{6 sin ϕ} , β = _{3−sin ϕ}^{6 cos ϕ} ϕ - friction angle

c_p(κ) - 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

Comp.Meth.Civ.Eng., II cycle

Slope stability simulation

Gradient-enhanced BDP plasticity

Interface finite elements in 2D and 3D

Example application: fracture in bending

Comp.Meth.Civ.Eng., II cycle

Interface elements

t = D ∆u

For 2D interface in 3D model:

t = [t_n t_t t_s]^T

Relative displacement ∆u be- tween faces (A) and (B) of inter- face element

∆u = [∆u_n ∆u_t ∆u_s]^T u = [u_n^(A)u_n^(B)u_t^(A)u^(B)_t u^(A)_s u_s^(B)]^T

∆u = Lu , u = N d^e

L =





−1 1 0 0 0 0

0 0 −1 1 0 0

0 0 0 0 −1 1





Interface in 2D model

Shear component

Simulation of cracking in masonry structures [3]

DIANA

Comp.Meth.Civ.Eng., II cycle

Modelling of masonry structures

”Micro”model (discrete interfaces) or ”macro”model (homogenization)

Theory of plasticity

Interface model

Continuum model

Comp.Meth.Civ.Eng., II cycle

Construction of underground in Amsterdam

Collapse of Sleipner A platform, Norway 1991

I RC off-shore platform, founded 82 m below sea level, composed of 24 cells with 12 m diameter (4 support deck)

I Cause of sinking: cracking and leakage of concrete base during

construction due to error in desing related to FEM computations of tricell - element connecting cells (shear force underestimated by o 47%) and insufficient anchoring of reinforcement in critical zone

Pictures from www.ima.umn.edu/∼arnold/disasters/sleipner.html

Comp.Meth.Civ.Eng., II cycle

Collapse of Terminal 2E of Paris Airport, 2004

I Terminal 2E of Paris Charles de Gaulle Airport is a composite

concrete-glass tube-like shell structure, large span freely supported vaulted roof is weakened by many openings.

I Designed by architect Paul Andreu (who also designed Terminal 3 at Dubai International Airport, which collapsed during construction), admitted to service in 2003.

I Combination of causes: too small safety margin in design, probably also mistakes during construction and/or insufficiently good concrete.

World Trade Center diseaster

Comp.Meth.Civ.Eng., II cycle

Collapse of World Trade Center

Built in 1966-77, 110 storeys 3.7m height each, steel frame structure based on „tube in

tubeśtructural concept, 26.5×41.8m core (47 columns connected by short beams, carried 60% of self weight), outre frame (240 box columns 356x356 every 1m along perimeter, carried 40% of self weight), composite floors supported on truss girders connected by hinges with the core and outer frame, hat trusses at levels 107-110 to connect perimeter and core columns. Building weight above ground 3630MN, self weight 2890MN, service load 740MN.

Simplified mechanism of WTC collapse [4,6]

Dynamic effect of high temperature which lowered steel yield strength and caused column buckling under creep con- ditions

1. Structure is weakened at impacted face, fuel fire causes temperature increase to about 600C 2. Stress redistribution takes place, viscoplastic

buckling of columns occurs at critical level 3. Floor trusses sag, buckling of columns grows,

frame joints are destroyed, about half of the columns cease to carry the weight of the storeys above

4. The upper part of the building falls to the floor below with growing kinetic energy, the impact is a dynamic load which cannot be carried by the structure below and progressive collapse begins 5. The upper part falls, its mass and energy grow

gradually

Energetic estimation of over-load coeffi- cient:

P_dyn/mg = 30 − 60.

Comp.Meth.Civ.Eng., II cycle

Response of WTC towers to abnormal load [5]

Simplified dynamic model (110 beam elements) with lumped masses (flood weight 33 MN), bending and shear stiffness determined from the actual space frame model, 3% Rayleigh damping

Results for simplified dynamic model

Displacements and forces did not exceed at impact those due to design wind load, hence the towers sustained the abnormal load at first.

Comp.Meth.Civ.Eng., II cycle

FEM model for detailed damage analysis - LS-DYNA [5]

Model of impact zone and airplane with mass of 140 000 kg, elastic-plastic material.

Critical segment model - results

Normal forces in perimeter columns before and after impact Redistribution of forces after impact, 122/113 destroyed for

WTC1/WTC2 resp., yield strength reached in remaining columns, significant stiffening influence of hat trusses.

Comp.Meth.Civ.Eng., II cycle

References

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

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

Wiley & Sons, Chichester, 2002.

[3] P.B. Lourenco. Computational strategies for masonry structures.

PhD Thesis, Delft University of Technology, 1996.

[4] Z.P. Bazant, Y. Zhou. Why Did the World Trade Center Collapse?

- Simple Analysis. ASCE J. Eng. Mech., 128, 2-6, 2002.

[5] Y. Omika, E. Fukuzawa, N. Koshika, H. Morikawa, R.

Fukuda. Structural Responses of World Trade Center Collapse under Aircraft Attacks. ASCE J. Eng. Mech., 131, 6-15, 2005.

[6] Z.P. Bazant, M. Verdure. Mechanics of Progressive Collapse:

Learning from World Trade Center and Building Demolitions. ASCE J.

Eng. Mech., 133, 308-319, 2007.