FEM analysis of shell structures CMCE Lecture 5/1, Civil Engineering, II cycle, specialty BEC

FEM analysis of shell structures

CMCE Lecture 5/1, Civil Engineering, II cycle, specialty BEC

Adam Wosatko

Institute for Computational Civil Engineering Civil Engineering Department, Cracow University of Technology

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

With thanks to:

Maria Radwańska, Anna Stankiewicz, Jerzy Pamin, Piotr Pluciński

Layout

Membranes (plane stress) Deep cantilever beam

Plates – conforming FE Conforming FE Square plate Corner plate

Shells

Some theoretical aspects Shell FE

Cylinder container

FE Q4 for membranes

Quadrilateral

u3

v₃ v4

u4

b

η = 2 y /b

ξ = 2 x /a X Y

a

u2

v2

v1

u1

4 3

1 2

NDOFN = 2 NNE = 4

NDOFE = NDOFN × NNE = 8

Displacement vectors for node and element:

qn= {un, vn}^T

q^e= {u1, v1|u2, v2|u3, v3|u4, v4}^T for n = 1, ..., NNE , e = 1, ..., NE

For approximation ofboth displacementsu(ξ, η) i v (ξ, η) we use shape functions N_i, i = 1, 2, 3, 4, which are bilinear (liniear in each direction for standard coordinates ξ, η ∈ [−1, +1]):

uⁿ(ξ, η)_(2×1)= {u(ξ, η), v (ξ, η)} = Nⁿ_(2×8)q^{e n}_(8×1) → u(ξ, η) = N₁u₁+ N₂u₂+ N₃u₃+ N₄u₄

v (ξ, η) = N₁v₁+ N₂v₂+ N₃v₃+ N₄v₄

Bilinear shape functions

for quadrilateral Q4

N₁=¹₄(1 − ξ)(1 − η) N₂= ¹₄(1 + ξ)(1 − η)

Stiffness matrix of FE Q4 for membranes

Stiffness matrix:

k^{e n}_(8×8)= Z Z

A^e

B^nT_(8×3)Dⁿ_(3×3)Bⁿ_(3×8)dA

Kinematic relation matrix:

Bⁿ_(3×8)=







∂N₁

∂x 0 ... ^∂N_∂x⁴ 0 0 ^∂N_∂y¹ ... 0 ^∂N_∂y⁴

∂N1

∂y

∂N1

∂x ... ^∂N_∂y^{4 ∂N}_∂x⁴







Constitutive matrix:

Dⁿ_(3×3)= Dⁿ





1 ν 0

ν 1 0

0 0 ^1−ν₂





Membrane stiffness: Dⁿ= _1−ν^{E h}2

Deep cantilever beam (ROBOT)

Membrane example

Mesh 4 x 16 Deformation

Deep cantilever beam

Contour plots

displacement UX [cm] displacement UY [cm]

stress σx [MPa] stress τxy [MPa]

Rectangular FE for plates

Rectangular, conforming

Degrees of freedom of a node:

q_w

[4×1]

= {w , ϕ_x, ϕ_y, χ}_w= {w , ∂w /∂y , −∂w /∂x , ∂²w /∂x ∂y }_w Degrees of freedom of an element:

q^e

[16×1]

= {q₁, q₂, q₃, q₄} = {w1, ϕ_{x 1}, ϕ_{y 1}, χ₁|w2. . . | . . . | . . . χ₄}

Approximation of deflection field

Dimensionless surface coordinates: ξ = 2 ^x_a
− 1, η = 2 ^y_b
− 1 Polynomial approximation:

w (ξ, η) = (α1+ α2ξ + α3ξ²+ α4ξ³)(β1+ β2η + β3η²+ β4η³) = C1+

C2ξ + C3η +

C4ξ²+ C5ξη + C6η²+

C7ξ³+ C8ξ²η + C9ξη²+ C10η³+ C11ξ³η + C12ξ²η²+ C13ξη³+ C₁₄ξ³η²+ C₁₅ξ²η³+

C₁₆ξ³η³ u^e

[1×1]= {w (ξ, η)} = N

[1×16](ξ, η) · q^e

[16×1]

=

= {N₁¹, N₁², N₁³, N₁⁴, |N₂¹. . . | . . . | . . . N₄⁴} · {w₁, ϕ_{x 1}, ϕ_{y 1}, χ₁|w₂. . . | . . . | . . . χ₄}^T

Bicubic shape functions

for first node of rectangular plate FE

N₁¹ corresponds to w₁ N₁²corresponds to ϕ_{x 1} N₁³corresponds to ϕ_{y 1}

N₁⁴corresponds to χ₁

Stiffness matrix of FE Q4 for plates

Stiffness matrix:

k^{e m}_(16×16) = Z Z

A^e

B^mT_(16×3)D^m_(3×3)B^m_(3×16)dA

Kinematic relation matrix:

B^m_(3×16)=









−^∂_∂x^2N2¹¹ −^∂_∂x^2N2¹² −^∂_∂x^2N2¹³ −^∂_∂x^2N2¹⁴ ... −^∂_∂x^2N2⁴⁴

−^∂_∂y^2N₂¹¹ −^∂_∂y^2N₂¹² −^∂_∂y^2N₂¹³ −^∂_∂y^2N₂¹⁴ ... −^∂_∂y^2N₂⁴⁴

−2_{∂x ∂y}^∂2N11 −2_{∂x ∂y}^∂2N21 −2_{∂x ∂y}^∂2N31 −2_{∂x ∂y}^∂2N41 ... −2^∂_{∂x ∂y}^2N44









Constitutive matrix:

D^m_(3×3)= D^m





1 ν 0

ν 1 0

0 0 ^1−ν₂





Bending stiffness: D^m= _12(1−ν^Eh3 2)

Square plate

5 7

4 K

Y

Z

X

Data:

I dimensions LX = LY = 3.0 m

I thickness h = 0.12 m

I load intensity

pz = q = −14.4 kN/m²

I Young modulus E = 32.0 · 10⁶ kPa

I Poisson coefficient ν = 0 Bending stiffness D^m=_12(1−ν^Eh3 ₂₎ = 4608 kNm.

Sign convention – ROBOT

Bending moments and transverse forces

Axis 0Z is oriented upwards, so positive bending moments cause tensile stresses at the top surface.

z

x

y tx

mx my

ty

m_yx mxy

Square plate – results

w7 mx 7 mx 5 my 5 my 4

[mm] _kNm

m

 _kNm

m

 _kNm

m

 _kNm

m

 Eng. comp. -2.303 -11.37 -6.58 -2.19 15.23 Package Mesh

ALGOR 8 × 8 -2.316 -11.27 -6.53 -2.15 13.20 GRAITEC 8 × 8 -2.344 -11.51 -6.64 -2.31 12.27 ROBOT 8 × 8 -2.355 -11.63 -6.73 -2.30 14.96 ALGOR 32 × 32 -2.344 -11.39 -6.59 -2.19 14.75 GRAITEC 32 × 32 -2.392 -11.44 -6.66 -2.18 14.47 ROBOT 32 × 32 -2.394 -11.50 -6.67 -2.21 15.19 ROBOT 128 × 128 -2.410 -11.47 -6.68 -2.16 15.29

Square plate – results

Contour plots mx and my

mx

my

ALGOR GRAITEC ROBOT

Corner plate (ROBOT) – data, applied meshes

I length of segment L = 4.0 m

I thickness h = 0.12 m

I Young modulus E = 10.0 · 10⁶kPa

I Poisson coefficient ν = 0.17 Deformation

Regular mesh Mesh refined at corner

Moments m

Regular mesh Mesh refined at corner

Moments m

Regular mesh Mesh refined at corner

Moments m

Regular mesh Mesh refined at corner

Principal stress directions – top surface

Regular mesh Mesh refined at corner

How many times should mesh be refined?

Element Concave corner mesh mx ,max= my ,maxkNm

m

 regular 38.12

refined once 47.22 refined many times 133.93

Engineering shell structures

Curvilinear roofs

http://www.ketchum.org Drilling vessel supports

Branch shells

Shell classification – slenderness

Shells:

I) thin → apply Kirchhoff-Love theory,

II) moderately thin → apply Mindlin-Reissner theory,

III) thick → compute 3D continuum model, e.g. problem of thick-walled pipe.

Three-parameter shell theory

Kirchhoff-Love – thin shells

Thin shell

Thickness is much smaller than minimum curvature radius:

h R_min < 1

20÷ 1 30. Kinematic hypothesis →

independent kinematic fields in shells are three displacements:

u1(ξ1, ξ2), u2(ξ1, ξ2), w (ξ1, ξ2), in local directions (e1, e2, n), introduced at point of middle surface.

Static hypothesis →

stresses σ_n, τ_1n, τ_2n in thin shells are much smaller than others.

Boundary conditions

Boundary conditions (static and/or kinematic) can be imposed:

I at points,

I along lines,

I on surfaces.

At edge points we can apply boundary loads which correspond to local boundary coordinates.

Types of kinematic constraints

Type of kinematic constraints

Free edge → different displacement constraints → fully clamped edge

Membrane state in shells

Membrane state – stress distribution along thickness is homogeneous.

This state is optimal in shell from the viewpoint of material use.

However, we must ensure suitable boundary conditions and loading.

Disturbances of membrane state

induce bending state

I Bending moments and transverse forces increase quickly as their origin is approached

Membrane state and its disturbances

p

boundary conditions boundary loads not smooth changes of thickness

p p

p

Shell FEs – classification

I FEs based on shell structures theory

I FEs 1D – geometrically one-dimensional

I Flat FEs 2D – three- or four-noded

I FEs 2D – curvilinear, based on Kirchhoff-Love theory for thin shells or Mindlin-Reissner theory for moderately thin shells

I so-called degenerate FEs, based on equations for 3D continuum, but with hypotheses as adopted for shells, in accordance with

five-parametric Mindlin-Reissner theory

I FEs for 3D solids – thick shells We consider:

I number of approximated fields in FE,

I number and type of degrees of freedom at node,

I formulas for internal energy and types of integration for stiffness matrix and external load vector.

Shell FE 1D

for axisymmetric problems

Geometrically one-dimensional finite element models meridian of

axisymmetric shell, but all functions are developed in trigonometric series in circumferential direction.

u(ξ, Θ) =



 u1

u2

w



=



 u v w



=

=



















J

P

j =1

u^j(ξ) cos(j Θ)

J

P

j =1

v^j(ξ) sin(j Θ)

J

Pw^j(ξ) cos(j Θ)



















Flat FE 2D for shells

three- and four-noded

Flat finite elements for shells with 3 or 4 nodes are formulated as a result of superposition of membrane and bending states.

+ =

+ =

We obtain the following degrees of freedom in node:

qⁿ_w= {u, v }_w q^m_w= {w , ϕ₁, ϕ₂}w q_w= {qⁿq^m}w= {u, v , w , ϕ₁, ϕ₂}w

Flat FE 2D for shells

implmented in ROBOT package

We should introduce third rotational degree of freedom in shells, especially if we compute shells with sharp edges or branching: ϕ_nis this sixth degree of freedom, hence:

qw= {u, v , w , ϕ1, ϕ2| ϕn}w.

Angleϕnis rotation around the normal to middle surface at node and is useful in tranformation of complete DOF set into other coordinate

Curvilinear shell FE 2D

based on classical three-parameter theory for thin shells

Approximation of displacement field (3 translation variables):

u(ξ₁, ξ₂) = {u₁(ξ₁, ξ₂), u₂(ξ₁, ξ₂), w (ξ₁, ξ₂)} =

= {u(ξ₁, ξ₂), v (ξ₁, ξ₂), w (ξ₁, ξ₂)} = N

[3×NDOFE ]

· q^e

[NDOFE ×1]

DOF vector at node:

qw= {u, v , w , ϕ1, ϕ2}w or qw= {u, v , w , ϕ1, ϕ2| ϕn}w

NDOFN = 6 NNE = 6 NDOFE = 36

NDOFN = 6 NNE = 8 NDOFE = 48

Elastic energy (internal):

U = Uⁿ+ U^m=¹₂RR

A

e^nTsⁿ
dA + ¹₂RR

A

e^mTs^m
dA

Cylinder container

Fourth-order displacement differential equation for axisymmetric membrane-bending state and its solution

w^IV(x ) + 4β⁴w (x ) = 1

D^m[pn+ ν R(

Z

p1dx + C1)]

where: β = r 1

Rh ·p⁴

3(1 − ν²)

w (x ) = w^o(x )

| {z }

general solution

+ w (x )

| {z }

particular solution

=

= e^−βx[B₁cos(βx ) + B₂sin(βx )]

| {z }

quickly increase with x

+ e^+βx[B₃cos(βx ) + B₄sin(βx )]

| {z }

quickly decrease with x

+ w (x )

w (x ) – normal displacement for membrane axisymmetric state w^o(x ) – normal displacement for local bending state

Cylinder container

Geometry, material, loading

I height L = H = 20.0 m

I radius R = 10.0 m

I thickness h = 0.12 m

I Young modulus E = 15.0 · 10⁶kPa

I Poisson coefficient ν = ¹₆≈ 0.1667

I dead load

p₁= −γ_bh = −24 · 0.12 = −2.88 kPa

I hydrostastic pressure

pn= γc(L − x ) = 10(20 − x ) kPa

I bending stiffness

D^m=_12(1−ν^Eh3 2) = 2221.7 kNm

I boundary conditions w (0) = 0, w⁰(0) = 0, m1(L) = 0, t1(L) = 0

Long shell?

β = 1.193 _m¹, λ = ^π_β ≈ 2.63 m λ - length of half-wave 3 · λ < L ?

≈ 7.9 m < 20 m ⇒ YES

Cylinder container – membrane state

Analytical solution

Normal dispalcement w

0 5 10 15 20

0 0.2 0.4 0.6 0.8 1 1.2

x [m]

w(x) [cm]

Meridional force n1 Circumferential force n2

10 15 20

x [m] 10

15 20

x [m]

Normal displacement function:

w (x ) = w (x ) = R

Eh(γ_cR + νγ_bh)(L−x )

Cylinder container – disturbance of membrane state

Analytical solution

Normal displacement w Transverse force t1 Meridional moment m1

0 5 10 15 20

0 0.2 0.4 0.6 0.8 1 1.2

x [m]

w(x) [cm]

0 5 10 15 20

−20 0 20 40 60 80 100 120 140 160

x [m]

t1(x) [kN/m]

0 5 10 15 20

−70−60−50−40−30−20−10 0 10 20

x [m]

m1(x) [kNm/m]

m1(0) = −67.6^kNm_m

Meridional force n1 Circumferential force n2

0 5 10 15 20

−60−50 −40 −30 −20 −10 0

x [m]

n1(x) [kN/m]

0 5 10 15 20

0 500 1000 1500 2000

x [m]

n2(x) [kN/m]

n2(2) = 1841^kN_m

Cylinder container – disturbance of membrane state

Numerical solution – ROBOT

Circumferential force n2 Transverse force t1 Meridional moment m1

Cylinder container – disturbance of membrane state

Numerical solution – ANSYS

Meridional moment m1 Circumferential moment m2

Meridional force n1 Circumferential force n2

Cylinder container – disturbance of membrane state

Summary

Analytical solution

wmax n2,max n2,min m1,max m1,min

[cm] _kN

m

 _kN

m

 _kNm

m

 _kNm

m

 position of x [m] 2.3 2.32 0.0 1.3 0.0

value 1.028 1841.3 -9.60 14.68 -67.64 Numerical solution

ANKA (FE 1D) 1.024 1835.9 -14.25 13.96 -67.67 ROBOT (FE 2D) 1.024 1846.8 -9.97 14.32 -60.52 ANSYS (FE 2.5D) 1.020 1847.8 -7.10 13.72 -58.19

Element type NE Active NDOFS

ANKA SRSK – 1D 20 80

ROBOT flat – 2D 1280 7680

References

M. Radwańska. Ustroje powierzchniowe. Podstawy teoretyczne oraz rozwiązania analityczne i numeryczne. Skrypt PK, Kraków, 2009 (in Polish).

R.D. Cook, D.S. Malkus, M.E. Plesha, R.J Witt. Concepts and Applications of Finite Element Analysis.

University of Wisconsin–Madison, John Wiley&Sons, 2002.

O.C. Zienkiewicz, R.L. Taylor, J.Z. Zhu. The Finite Element Method:

Its Basis and Fundamentals. VI edition, Elsevier Butterworth Heineman, 2005.

M. Radwańska, A. Stankiewicz, A. Wosatko, J. Pamin. Plate and Shell Structures. Selected Analytical and Finite Element Solutions. John Wiley &

Sons, 2017.