# 6.1.4 Shell model¶

## Simple Naghdi shell model¶

Geometric model and meshing. Clamped on left boundary.

[1]:

from netgen.csg import *
from ngsolve import *
from ngsolve.internal import visoptions
from ngsolve.webgui import Draw

order = 3

geo = CSGeometry()
cyl   = Cylinder(Pnt(0,0,0),Pnt(1,0,0),0.4).bc("cyl")
left  = Plane(Pnt(0,0,0), Vec(-1,0,0))
right = Plane(Pnt(1,0,0), Vec(1,0,0))
finitecyl = cyl * left * right
geo.NameEdge(cyl,left, "left")
geo.NameEdge(cyl,right, "right")

mesh = Mesh(geo.GenerateMesh(maxh=0.2))
mesh.Curve(order)
Draw(mesh)

[1]:

BaseWebGuiScene


Use Lagrangian elements for displacement $$u \in [H^1(S)]^3$$ and the rotation $$\beta \in [H^1(S)]^3$$. It might lock for small thickness $$t$$.

[2]:

fes1 = VectorH1(mesh, order=order, dirichlet_bbnd="left")
fes = fes1*fes1
u,beta = fes.TrialFunction()

nsurf = specialcf.normal(3)

thickness = 0.1


Membrane energy

$t\|E_{tt}(u)\|^2_{L^2(S)}$

Shear energy

$t\int_S | \nabla u^\top n - \beta |^2$

Bending energy

$\frac{t^3}{2}\|\boldsymbol{\varepsilon}(\beta)-\text{Sym}(\nabla u^\top\nabla\nu)\|^2_{L^2(S)}$
[3]:

Ptau = Id(3) - OuterProduct(nsurf,nsurf)
Ctautau = Ftau.trans * Ftau
Etautau = 0.5*(Ctautau - Ptau)

#Average normal vector for affine geometry
if order == 1:
gfn = GridFunction(fes1)
gfn.Set(nsurf,definedon=mesh.Boundaries(".*"))
else:
gfn = nsurf

a = BilinearForm(fes, symmetric=True)
#membrane energy
a += Variation( thickness*InnerProduct(Etautau, Etautau)*ds )
#bending energy
#shearing energy

# external force
factor = Parameter(0.0)
a += Variation( -thickness*factor*y*u[1]*ds )

gfu = GridFunction(fes)


Increase the load step-wise, solve the non-linear problem by Newton’s method. First and second order derivatives are computed by automatic differentiation.

[4]:

with TaskManager():
solvers.NewtonMinimization(a, gfu, printing=False)

loadstep  0

[5]:

Draw(gfu.components[1], mesh, "rotations", deformation=gfu.components[0])
Draw(gfu.components[0], mesh, "disp")

[5]:

BaseWebGuiScene


## Nonlinear Koiter shell model¶

We present the method described in [Neunteufel and Schöberl. The Hellan-Herrmann-Johnson method for nonlinear shells. Computers & Structures , 225 (2019), 106109].

[6]:

from math import pi
from ngsolve.meshes import MakeStructuredSurfaceMesh
thickness = 0.1
L = 12
W = 1
E, nu = 1.2e6, 0
moment = IfPos(x-L+1e-6, 1, 0)*50*pi/3

mapping = lambda x,y,z : (L*x, W*y,0)
mesh = MakeStructuredSurfaceMesh(quads=False, nx=10, ny=1, mapping=mapping)
Draw(mesh)

[6]:

BaseWebGuiScene


To avoid membrane locking Regge interpolation as in [Neunteufel and Schöberl. Avoiding Membrane Locking with Regge Interpolation] can be used.

[7]:

# False -> membrane locking
regge = True
order = 2

fes1 = HDivDivSurface(mesh, order=order-1, discontinuous=True)
fes2 = VectorH1(mesh, order=order, dirichletx_bbnd="left", dirichlety_bbnd="left|bottom", dirichletz_bbnd="left")
fes3 = HDivSurface(mesh, order=order-1, orderinner=0, dirichlet_bbnd="left")
if regge:
fes4 = HCurlCurl(mesh, order=order-1, discontinuous=True)
fes  = fes2*fes1*fes3*fes4*fes4
u,sigma,hyb,C,R = fes.TrialFunction()
sigma, hyb, C, R = sigma.Trace(), hyb.Trace(), C.Trace(), R.Operator("dualbnd")
else:
fes  = fes2*fes1*fes3
u,sigma,hyb = fes.TrialFunction()
sigma, hyb = sigma.Trace(), hyb.Trace()

fesVF = VectorFacetSurface(mesh, order=order)
b = fesVF.TrialFunction()

gfclamped = GridFunction(FacetSurface(mesh,order=0))
gfclamped.Set(1,definedon=mesh.BBoundaries("left"))

solution = GridFunction(fes, name="solution")
averednv = GridFunction(fesVF)
averednv_start = GridFunction(fesVF)

nsurf = specialcf.normal(mesh.dim)
t     = specialcf.tangential(mesh.dim)
nel   = Cross(nsurf, t)

Ptau    = Id(mesh.dim) - OuterProduct(nsurf,nsurf)
Ctau    = Ftau.trans*Ftau
Etautau = 0.5*(Ctau - Ptau)

nphys   = Normalize(Cof(Ftau)*nsurf)
tphys   = Normalize(Ftau*t)
nelphys = Cross(nphys,tphys)

Hn = CoefficientFunction( (u.Operator("hesseboundary").trans*nphys), dims=(3,3) )

cfn  = Normalize(CoefficientFunction( averednv.components ))
cfnR = Normalize(CoefficientFunction( averednv_start.components ))
pnaverage = Normalize( cfn - (tphys*cfn)*tphys )

$\sum_{T\in \mathcal{T}_h}\int_{\partial T} b\cdot\delta b\,ds = \sum_{T\in \mathcal{T}_h}\int_{\partial T} \nu^n\cdot\delta b\,ds,\qquad \forall \delta b$
[8]:

bfF = BilinearForm(fesVF, symmetric=True)
bfF += Variation( (0.5*b*b - ((1-gfclamped)*cfnphys+gfclamped*nsurf)*b)*ds(element_boundary=True))
rf = averednv.vec.CreateVector()
bfF.Apply(averednv.vec, rf)
bfF.AssembleLinearization(averednv.vec)
invF = bfF.mat.Inverse(fesVF.FreeDofs(), inverse="sparsecholesky")
averednv.vec.data -= invF*rf
averednv_start.vec.data = averednv.vec

[9]:

gradn = specialcf.Weingarten(3) #grad(nsurf)

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/(2*nu+1)*Trace(mat)**2)

[10]:

bfA = BilinearForm(fes, symmetric=True, condense=True)
bfA += Variation( (-6/thickness**3*MaterialNormInv(sigma, E, nu) \
+ InnerProduct(sigma, Hn + (1-nphys*nsurf)*gradn))*ds ).Compile()
if regge:
bfA += Variation( 0.5*thickness*MaterialNorm(C, E, nu)*ds )
bfA += Variation( InnerProduct(C-Etautau, R)*ds(element_vb=BND) )
bfA += Variation( InnerProduct(C-Etautau, R)*ds(element_vb=VOL) )
else:
bfA += Variation( 0.5*thickness*MaterialNorm(Etautau, E, nu)*ds )
bfA += Variation( -(acos(nel*cfnR)-acos(nelphys*pnaverage)-hyb*nel)*(sigma*nel)*nel*ds(element_boundary=True) ).Compile()
par = Parameter(0.0)
bfA += Variation( -par*moment*(hyb*nel)*ds(element_boundary=True) )

[11]:

par.Set(0.1)
bfF.Apply(averednv.vec, rf)
bfF.AssembleLinearization(averednv.vec)
invF.Update()
averednv.vec.data -= invF*rf
solvers.Newton(bfA, solution, inverse="sparsecholesky", maxerr=1e-10, maxit=20)

Newton iteration  0
err =  1.81379936422795
Newton iteration  1
err =  106.02124724341046
Newton iteration  2
err =  6.430451092532192
Newton iteration  3
err =  0.07240855358895455
Newton iteration  4
err =  4.902465394446003e-06
Newton iteration  5
err =  9.450264246561326e-14

[12]:

Draw(solution.components[0], mesh, "disp", deformation=solution.components[0])

[12]:

BaseWebGuiScene

[13]:

numsteps=10
for steps in range(1,numsteps):
par.Set((steps+1)/numsteps)
print("Loadstep =", steps+1, ", F/Fmax =", (steps+1)/numsteps*100, "%")

bfF.Apply(averednv.vec, rf)
bfF.AssembleLinearization(averednv.vec)
invF.Update()
averednv.vec.data -= invF*rf

(res,numit) = solvers.Newton(bfA, solution, inverse="sparsecholesky", printing=False, maxerr=2e-10)

Loadstep = 2 , F/Fmax = 20.0 %
Loadstep = 3 , F/Fmax = 30.0 %
Loadstep = 4 , F/Fmax = 40.0 %
Loadstep = 5 , F/Fmax = 50.0 %
Loadstep = 6 , F/Fmax = 60.0 %
Loadstep = 7 , F/Fmax = 70.0 %
Loadstep = 8 , F/Fmax = 80.0 %
Loadstep = 9 , F/Fmax = 90.0 %
Loadstep = 10 , F/Fmax = 100.0 %

[14]:

Draw(solution.components[0], mesh, "disp", deformation=solution.components[0])

[14]:

BaseWebGuiScene

[ ]: