9. Shells

Shells are the “primadonna of structures” (Ramm, 1986). On the one hand, if properly designed, they can exhibit an optimal ratio of stiffness to weight. On the other hand, especially for thin shells, they can be extremely sensitive to imperfections and develope unstable behavior.

In this chapter we present efficient discretization methods to simulate nonlinear shells.

• Introduction of shell energies and (non-)linear shell models

• Linear shell models

• Nonlinear shell models

• Shell buckling

Example: Sydney opera hall

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw

# first fanel
v1l = Pnt(0, -10, 0.1)
v1t = Pnt(-6, 0, 10.5)
w1 = Wire(ArcOfCircle(v1l, Vec(-0.3, 0.4, 1), v1t))

v12l = Pnt(0.25, -10, 0.1)
v12t = Pnt(-0.5, 0, 11)
w12 = Wire(ArcOfCircle(v12l, Vec(-0.25, 0.3, 1), v12t))

v2l = Pnt(0.5, -10, 0.1)
v2t = Pnt(5 - 0.15, 0, 10)
w2 = Wire(ArcOfCircle(v2l, Vec(0, 0.2, 1), v2t))

# vent
v1lv = Pnt(1, -10.5, 1.5)
v1tv = Pnt(5, 0, 10)
w2v = Wire(ArcOfCircle(v1lv, Vec(0, 0.4, 1), v1tv))

v2lv = Pnt(4.5, -10.5, 0.5)
v2tv = v1tv
w3v = Wire(ArcOfCircle(v2lv, Vec(0, 0.2, 1), v2tv))

v3l = Pnt(5, -10, 1)
v3t = Pnt(5 + 0.15, 0, 10)
w3 = Wire(ArcOfCircle(v3l, Vec(0, 0, 1), v3t))

v4l = Pnt(5.5, -10 - 1, 1)
v4t = Pnt(3.5, 0, 15)
w4 = Wire(ArcOfCircle(v4l, Vec(0, 0, 1), v4t))

# second panel
v212l = Pnt(5.5 + 0.25, -10 - 1, 1)
v212t = Pnt(3.5 + 5.5, 0, 11 + 2)
w212 = Wire(ArcOfCircle(v212l, Vec(0, 0, 1), v212t))

v22l = Pnt(5.5 + 0.5, -10 - 1, 1)
v22t = Pnt(3.5 + 11 - 0.15, 0, 10 + 1)
w22 = Wire(ArcOfCircle(v22l, Vec(0.2, 0, 1), v22t))

# vent
v21lv = Pnt(5.5 + 0.5 + 0.5, -10 - 1.5, 1.5)
v21tv = Pnt(3.5 + 11, 0, 10 + 1)
w22v = Wire(ArcOfCircle(v21lv, Vec(0.3, 0.1, 1), v21tv))

v22lv = Pnt(3.5 + 11 - 1.5, -10 - 1.5, 1)
v22tv = v21tv
w23v = Wire(ArcOfCircle(v22lv, Vec(-0.2, 0, 1), v22tv))


v23l = Pnt(3.5 + 11, -10 - 1, 0)
v23t = Pnt(3.5 + 11 + 0.15, 0, 10 + 1)
w23 = Wire(ArcOfCircle(v23l, Vec(0, 0, 1), v23t))

v24l = Pnt(3.5 + 11, -10 - 1.5, 0)
v24t = Pnt(11, 0, 14 + 8)
w24 = Wire(ArcOfCircle(v24l, Vec(-0.1, 0, 1), v24t))

# third panel
v312l = Pnt(3.5 + 11 + 0.25, -10 - 1.5, 0)
v312t = Pnt(11 + 9, 0, 11 + 6.5)
w312 = Wire(ArcOfCircle(v312l, Vec(0.2, 0, 1), v312t))

v32l = Pnt(3.5 + 11 + 0.5, -10 - 1.5, 0)
v32t = Pnt(11 + 15 - 0.25, 0, 9)
w32 = Wire(ArcOfCircle(v32l, Vec(1, 0, 1), v32t))

# fouth panel

# vent
v32lv = Pnt(3.5 + 11 + 0.5 + 4, -10 - 1.5 - 0.5, 3)
v32tv = Pnt(11 + 15 - 0.25, 0, 9)
w32v = Wire(ArcOfCircle(v32lv, Vec(0.6, 0.3, 1), v32tv))

v33lv = Pnt(3.5 + 11 + 0.5 + 9, -10 - 1.5 - 0.5, 2)
v33tv = Pnt(11 + 15, 0, 9)
w33v = Wire(ArcOfCircle(v33lv, Vec(0, 0.3, 1), v33tv))

v34lv = Pnt(3.5 + 11 + 17 - 2, -10 - 1.5 - 0.5, 3)
v34tv = Pnt(11 + 15, 0, 9)
w34v = Wire(ArcOfCircle(v34lv, Vec(-0.6, 0.3, 1), v34tv))


v41l = Pnt(3.5 + 11 + 17, -10 - 1.5, 0)
v41t = Pnt(11 + 15 + 0.25, 0, 9)
w41 = Wire(ArcOfCircle(v41l, Vec(-0.5, 0, 1), v41t))

v42l = Pnt(3.5 + 11 + 17 + 0.25, -10 - 1.5, 0)
v42t = Pnt(3.5 + 11 + 17 + 0.25 - 0.5, 0, 11.75)
w42 = Wire(ArcOfCircle(v42l, Vec(0, 0, 1), v42t))

v43l = Pnt(3.5 + 11 + 17 + 0.5, -10 - 1.5, 0)
v43t = Pnt(3.5 + 11 + 17 + 0.5 + 5, 0, 12)
w43 = Wire(ArcOfCircle(v43l, Vec(0.25, 0, 1), v43t))

shape = Glue(
    [w1, w12, w2, w2v, w3v, w3, w4]
    + [w212, w22, w22v, w23v, w23, w24]
    + [w312, w32, w32v, w33v, w34v]
    + [w41, w42, w43]
)
Draw(shape);

f1 = ThruSections([w2, w12, w1], solid=False)
f2 = ThruSections([w2v, w2], solid=False)
v1 = ThruSections([w3v, w2v], solid=False)
v2 = ThruSections([w3, w3v], solid=False)
f3 = ThruSections([w4, w3], solid=False)

f4 = ThruSections([w22, w212, w4], solid=False)
f5 = ThruSections([w22v, w22], solid=False)
v3 = ThruSections([w23v, w22v], solid=False)
v4 = ThruSections([w23, w23v], solid=False)
f6 = ThruSections([w24, w23], solid=False)

f7 = ThruSections([w32, w312, w24], solid=False)
v5 = ThruSections([w32v, w32], solid=False)
v6 = ThruSections([w33v, w32v], solid=False)
v7 = ThruSections([w34v, w33v], solid=False)

f8 = ThruSections([w41, w34v], solid=False)
f9 = ThruSections([w43, w42, w41], solid=False)
shell = Glue([f1, f2, f3, f4, f5, f6, f7, f8, f9, v1, v2, v3, v4, v5, v6, v7])

shell = Glue([shell, shell.Mirror(Axes((0, 0, 0), Y))])

shell.edges[Z <= 0.2].name = "fix"

shell.vertices.hpref = 1
shell.edges.hpref = 1
shell.faces.col = (222/255, 224/255, 204/255,1)

Draw(shell)
print ("(model created by Sebastian Hirnschall)")

(model created by Sebastian Hirnschall)

mesh = Mesh(OCCGeometry(shell).GenerateMesh(maxh=2))
mesh.RefineHP(1)
mesh.Curve(3)
Draw(mesh);

thickness = 0.1
order = 4

E = 2.85e4
nu = 0.3
mu = E / (2 * (1 + nu))
lam = E * nu / (1 - nu**2)


def MaterialNorm(mat, E, nu):
    return E / (1 - nu**2) * ((1 - nu) * InnerProduct(mat, mat) + nu * Trace(mat) ** 2)


def MaterialNormInv(mat, E, nu):
    return (1 + nu) / E * (InnerProduct(mat, mat) - nu / (nu + 1) * Trace(mat) ** 2)

fesU = VectorH1(mesh, order=order, dirichlet_bbnd="fix")
fesS = HDivDivSurface(mesh, order=order - 1, discontinuous=True)
fesH = NormalFacetSurface(mesh, order=order - 1)  # , dirichlet_bbnd="fix")

fes = fesU * fesS * fesH
(u, sigma, uh) = fes.TrialFunction()
# Need to take traces as we are on the surface
sigma, uh = sigma.Trace(), uh.Trace()

fesRegge = HCurlCurl(mesh, order=order - 1)

# normal, tangential, and co-normal vectors
nv = specialcf.normal(3)
tv = specialcf.tangential(3)
cnv = Cross(nv, tv)

# Projection to tangent space and surface derivatives
Ptau = Id(3) - OuterProduct(nv, nv)
grad_u = Grad(u).Trace()
Ftau = grad_u + Ptau
Etautau = 0.5 * (Ftau.trans * Ftau - Ptau)

# sum_{i=1}^3 hesse(u_i) nv_i
Hn = (u.Operator("hesseboundary").trans * nv).Reshape((3, 3))

bfa = BilinearForm(fes, symmetric=True, condense=True)
# membrane part
bfa += Variation(
    thickness / 2 * MaterialNorm(Interpolate(Etautau, fesRegge), E, nu) * ds
)
# bending part
bfa += Variation(-12 / thickness**3 / 2 * MaterialNormInv(sigma, E, nu) * ds)
bfa += Variation(InnerProduct(sigma, Hn) * ds)
bfa += Variation(
    sigma[cnv, cnv] * (uh - grad_u[nv, :]) * cnv * ds(element_boundary=True)
)

force = CF((0, 0, -0.01))
bfa += Variation(-force * u * ds)

gfu = GridFunction(fes)

with TaskManager():
    solvers.Newton(a=bfa, u=gfu, maxerr=1e-8, printing=False)

Draw(gfu.components[0], scale=10, deformation=True);

Draw(Norm(gfu.components[1]), mesh, order=3, max=0.2);