This page was generated from unit-6.1.4-shells/shell.ipynb.

6.1.4 Shell model

Simple Naghdi shell model

Geometric model and meshing. Clamped on left boundary.

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.AddSurface(cyl, finitecyl)
geo.NameEdge(cyl,left, "left")
geo.NameEdge(cyl,right, "right")

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

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\).

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)}\]

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

eps_beta = Sym(Ptau*grad(beta).Trace())
gradu = grad(u).Trace()
ngradu = gradu.trans*nsurf
#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
a += Variation( 0.5*thickness**3*InnerProduct(eps_beta-Sym(gradu.trans*grad(gfn)),eps_beta-Sym(gradu.trans*grad(gfn)))*ds )
#shearing energy
a += Variation( thickness*(ngradu-beta)*(ngradu-beta)*ds )

# 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.

with TaskManager():
    for loadstep in range(6):
        print("loadstep ", loadstep)
        factor.Set (1.5*(loadstep+1))
        solvers.NewtonMinimization(a, gfu, printing=False)

loadstep  0
loadstep  1
loadstep  2
loadstep  3
loadstep  4
loadstep  5

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

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].

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)

BaseWebGuiScene

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

# 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)
Ftau    = grad(u).Trace() + Ptau
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) )

cfnphys = Normalize(Cof(Ptau+grad(solution.components[0]))*nsurf)

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\]

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

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)

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) )

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

Newton iteration  0
err =  1.8137993642278916
Newton iteration  1
err =  106.02124724338891
Newton iteration  2
err =  6.430451091512385
Newton iteration  3
err =  0.07240855382933263
Newton iteration  4
err =  4.902465485748893e-06
Newton iteration  5
err =  1.416202596960721e-13

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

BaseWebGuiScene

numsteps=10
with TaskManager():
    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 %

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

BaseWebGuiScene