Solution of nonlinear problems

Solution of nonlinear problems

CMCE Lecture 3, 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:

A. Wosatko, A. Winnicki

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

Comp.Meth.Civ.Eng., II cycle

Lecture scope

Nonlinear problems Geometrical nonlinearity

Theory of moderately large deflections Physical nonlinearity

Cracking

Nonlinear analysis of RC shells Final remarks

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.

Comp.Meth.Civ.Eng., II cycle

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 δW_int = δW_ext

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 =

n_w

X

i =1

N_i(ξ, η, ζ)d_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

ne

X

e=1

A^{e T} Z

V^e

B^Tσ dV = f_ext, B = LN

Isoparametric approach, numerical integration

Comp.Meth.Civ.Eng., II cycle

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 = ¹₂(u_{i ,j} + u_{j ,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^e d = f_ext, Kd = f_ext

Incremental-iterative analysis

Nonlinear problem:

f_ext 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

ne

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σ^t dV Linearization of the left-hand side at time t

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

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

Comp.Meth.Civ.Eng., II cycle

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₁

f_int,1

u u^t+∆t

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

∆u₁ du2

∆f_ext

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

∆d1=K⁻¹₀ (f_ext^t+∆t− f_int,0^t ) σ1 → f_int,1^t+∆t 6= 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_ext^t+∆t−f_int,j^t+∆tk 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: R_j = f_ext^t+∆t − f_int,j^t+∆t → 0

f

R₁

f_int,1

u u^t+∆t

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

∆u₁ du₂

∆fext

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

∆d1 = K⁻¹₀ (f_ext^t+∆t− f_int,0^t ) σ1 → f_int,1^t+∆t 6= 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_ext^t+∆t−f^t+∆t int,j k k∆f_extk ¬ δ Modified algorithm:

Kj = K0

Comp.Meth.Civ.Eng., II cycle

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₁

f_int,1

u u^t+∆t

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

∆u₁ du2

∆f_ext

K ∆d = f_ext^t+∆t − f_int^t K - tangent operator First iteration:

∆d1 = K⁻¹₀ (f_ext^t+∆t− fint,0) σ1 → f_int,1^t+∆t 6= 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_ext^t+∆t−f_int,j^t+∆tk 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: 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₂

∆fext

R2

fint,2

K ∆d = f_ext^t+∆t − f_int^t K - tangent operator First iteration:

∆d1 = K⁻¹₀ (f_ext^t+∆t− fint,0) σ1 → f_int,1^t+∆t 6= 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_ext^t+∆t−f^t+∆t int,j k k∆f_extk ¬ δ Modified algorithm:

Kj = K0

Comp.Meth.Civ.Eng., II cycle

Options of incremental loading control

Force or displacement control

Arc length control

Geometrical nonlinearity

X

u

x φ(X, t)

V

V₀

S S₀

x₁, X₁ x₂, X₂

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]

Comp.Meth.Civ.Eng., II cycle

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

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

where tangent stiffness matrix:

K = K₀+K_u +K_σ K₀ - linear stiffness matrix

K_u - initial displacement matrix

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

Comp.Meth.Civ.Eng., II cycle

Karman theory of moderately large deflections [2]

Deflection of the order of thickness admitted Medium plane deflections

ε_x = ε^L_x + ε^N_x = ^∂u_∂x + ¹₂ ^∂w_∂x
² ε_y = ε^L_y + ε^N_y = ^∂v_∂y + ¹₂ 

∂w

∂y

² γ_xy = γ_xy^L + γ_xy^N = ^∂u_∂y + ^∂v_∂x + ^∂w_∂x ^∂w_∂y Curvatures and twist as in linear theory

κ_x = κ^L_x = −^∂_∂x^2w2 , κ_y = κ^L_y = −^∂_∂y^2w2, χ_xy = χ^L_xy = −2_{∂x ∂y}^∂2w

Two equations of Karman theory for moderately large plate deflections

∇²∇²F (x , y ) + ^Eh₂ L(w , w ) = 0 D^m∇²∇²w (x , y ) − L(w , F ) − ˆp_z = 0 where:

F (x , y ) - stress function (n_x = F_,yy, n_y = F_,xx, n_xy = −F_,xy) L(a, b) = a_,xxb_,yy − 2a_,xyb_,xy + a_,yyb_,xx

Theory of moderately large deflections

FEM approximation [2]

Total potential energy (additional membrane state energy) Π˜^m = U^m + ˜Uⁿ − W^m

Discretization uⁿ =

 u(x , y ) v (x , y )



=

 N_u 0 N_v 0

  dⁿ d^m



w^m = [w (x , y )] =

0 N_w 

 dⁿ d^m



Nonlinear kinematic equations B_{(6×LSSE )} = B^L_{(6×LSSE )}+ B^N_{(6×LSSE )}

B^L_{(6×LSSE )} =

 Bⁿ 0 0 B^m



, B^N_{(6×LSSE )} =

 0 Bⁿ_w 0 0



Comp.Meth.Civ.Eng., II cycle

Theory of moderately large deflections

FEM approximation [2]

Matrices of discrete kinematic relations

Bⁿ =





N_u,x N_{v ,y} N_u,y + N_{v ,x}



, B^m =





−N_{w ,xx}

−N_{w ,yy}

−2N_{w ,xy}





Bⁿ_w =





w_,xN_{w ,x} w_,yN_{w ,y} w_,xN_{w ,y} + w_,yN_{w ,x}



 Deflection gradients

g =

 w_,x w,y



=

 N_{w ,x} Nw ,y



d^m = G^md^m

Theory of moderately large deflections

FEM approximation [2]

Element tangent stiffness k^e_T = k^e₀+ k^e_u + k^e_σ

k^e₀ = Z Z

A^e

 B^nTDⁿBⁿ 0 0 B^mTD^mB^m

 dA

k^e_u = Z Z

A^e

 0 B^nTDⁿBⁿ_w B^nT_w DⁿBⁿ 0

 dA

k^e_σ = Z Z

A^e

 0 0

0 G^mTSⁿG^m



dA, Sⁿ =

 n_x n_xy nxy ny



Vector of nodal forces representing stresses f_int^e =

Z Z

A^e

 B^nT 0 B^nT_w B^mT

  Sⁿ S^m

 dA

Comp.Meth.Civ.Eng., II cycle

Square plate [2]

Moderately large deflections

ˆ p_z

x

z, w y

a

a C

Dane:

L = Lx = Ly = 2a = 1.0 m h = 0.002 m

E = 200000 MPa, ν = 0.25

Relation deflection-loading:

0.001 0.003

0.30

0.15

0.002

w_C [m]

ˆ p_z [kPa]

FEM

(geometrical nonlinearity)

wC = 0.00406^pˆz_D^L4 linear solution:

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

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

equivalent plastic strain distribution plastified

elastic material

Discrete and smeared cracks (DIANA)

Fracture energy G_f (used for a unit surface area of a crack to be formed)

Comp.Meth.Civ.Eng., II cycle

Simulation of cracking in RC panel

with ATENA [3]

Load-carrying capacity of multi-span panels (DIANA) [4]

Comp.Meth.Civ.Eng., II cycle

Geometrically and physically nonlinear analysis [5]

Scheme of computation strategy at the levels of:

I structure I finite element I layer

I point

Effects considered:

I stress evolution in cross-section I elastic-cracking concrete I elastic-plastic reinforcement I large displacements and their

gradients Aims:

I computation of displacement evolution

I determination of damage mechanism

I estimation of load-carrying capacity

RC shell model

I Degenerated 8-noded shell element (Mindlin-Reissner theory)

I Layered RC shell model (5 concrete layers, 4 steel layers representing two reinforcement grids)

I Continuum elastic-cracking model for concrete layer (concrete softening, reduction of shear stiffness)

I Elastic-plastic model for steel layer

Comp.Meth.Civ.Eng., II cycle

Numerical analysis of cooling tower shell [5]

Numerical analysis of cooling tower shell

Diagrams λ − w_K obtained using two FEM packages using force or displacement control for loading g + λ(w + s)

Cooling tower loads:

I self-weight g

I wind w

I internal suction s

I temperature variations

I subsidence

Comp.Meth.Civ.Eng., II cycle

Numerical analysis of RC shell

Results for shell with technological opening

Deformation Membrane forces along longitude lines

Numerical analysis of RC shell

Directions of principal stresses in external layer Smeared cracks visualization

Comp.Meth.Civ.Eng., II cycle

Damaged cooling tower shell [6,7]

Analized cases for load combination g + λ(w + s)

I designed shell

I built shell with zones of weak concrete (f_cm=11 MPa)

I shell with two circumferential openings (25m and 14m in length)

I repaired and strengthened shell (5cm reinforced shotcrete in height zone 18-40m)

Comp.Meth.Civ.Eng., II cycle

Nonlinear analysis results

The construction error did not have a significant influence on the

short-term load carrying capacity of the cooling tower, but it affected its durability due to local concrete overload and reinforcement corrosion.

Cracking zone prediction (DIANA)

Cracking zones in inner and outer concrete layer for λ = 3.2 (DIANA)

Comp.Meth.Civ.Eng., II cycle

Model of shell with holes (DIANA)

For smeared cracking computations diverged for λ ≈ 1.0 - it is necessary to perform mesh refinement and use a more stable material model.

Final remarks

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

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

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

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

considering plasticity or cracking.

5. 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 (or characteristic) values of loading and strength.

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. Radwańska. Ustroje powierzchniowe, podstawy teoretyczne oraz rozwiązania analityczne i numeryczne. Wydawnictwo PK, Kraków, 2009.

[3] M. Kwasek Advanced static analysis and design of reinforced concrete deep beams.

Diploma work, Politechnika Krakowska, 2004.

[4] M. Asin The behaviour of reinforced concrete continuous deep beams. PhD Thesis, Delft University of Technology, 2000.

[5] Z. Waszczyszyn, E. Pabisek, J. Pamin, M. Radwańska. Nonlinear analysis of a RC cooling tower with geometrical imperfections and a technological cut-out. Engineering Structures, 2, 480-489, 2000.

[6] A. Moroński. Analiza zarysowania i utraty stateczności uszkodzonej powłoki żelbetowej chłodni kominowej. Praca dyplomowa, Politechnika Krakowska, Kraków, 1996.

[7] A. Moroński, J. Pamin, M. Płachecki, Z. Waszczyszyn. Fracture and loss of stability of a partly-damaged cooling tower shell. Proc. 2nd Int. DIANA Conf. on Finite Elements in Engineering and Science, Eds M.A.N. Hendriks et al, 107-110, Kluwer Academic Publishers, Dordrecht, 1997.