Naghdi shell

4. Naghdi shell#

Joint work with M. Neunteufel [‘19, ‘21], Phd thesis Neunteufel (2021)

Geometric model and meshing. left surface only needed for boundary conditions.

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

from netgen.occ import *
cyl = Cylinder(Axes((0,0,0), X, Y), r=0.4, h=1).faces.Nearest((0.5,0.4,0))
cyl.edges.Min(X).name="left"

mesh = Mesh(OCCGeometry(cyl).GenerateMesh(maxh=0.1)).Curve(3)
Draw (mesh);

4.1. unknown field variables:#

  • displacement u∈[H1(S)]3

  • linearized rotation β∈H(curl,S)

  • bending moments σ∈H(div div,S)

order = 3

fes1 = HDivDivSurface(mesh, order=order-1)
fes2 = VectorH1(mesh, order=order, dirichlet_bbnd="left")
fes3 = HCurl(mesh, order=order-1, dirichlet_bbnd="left")
fes = fes1 * fes2 * fes3

sigma,u,beta = fes.TrialFunction()
tau,v,delta = fes.TestFunction()

sigma = sigma.Trace()
tau = tau.Trace()
beta = beta.Trace()
delta = delta.Trace()

gradv = Grad(v).Trace()
gradu = Grad(u).Trace()

nsurf = specialcf.normal(3)
t = specialcf.tangential(3)

nel = Cross(nsurf, t)

ngradv = gradv.trans*nsurf
ngradu = gradu.trans*nsurf

sigman = sigma*nel
taun = tau*nel

def tang(u):
    return (u*t)*t

thickness = 0.1

Membrane energy with primal method for displacement vector u:

∫SW(Ctt(u))

Shear energy, β the linearized rotation in tangential plane:

∫S|nT∇u−β|2

Bending energy (with TDNNS mixed method):

t2∫S|εtt(β)|2
a = BilinearForm(fes)

bending = -1/thickness**2*InnerProduct(sigma,tau)*ds + \
    (-div(sigma)*delta - div(tau)*beta)*ds + \
    (sigman * tang(delta) + taun*tang(beta))*ds(element_boundary=True)

shear = (ngradu-beta)*(ngradv-delta)*ds

a += (bending+shear).Compile()

Pt = Id(3) - OuterProduct(nsurf, nsurf)
Ft = gradu + Pt
Ctt = Ft.trans * Ft
Ett = Ctt - Pt
a += Variation((InnerProduct(Ett,Ett)*ds).Compile())

factor = Parameter(0.1)
a += Variation((-factor*y*u[1]*ds).Compile())

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.

sceneu = Draw(Norm(gfu.components[0]), mesh, deformation=gfu.components[1]) 
from ngsolve.solvers import NewtonMinimization

for loadstep in range(4):
    print ("loadstep ", loadstep)
    factor.Set (3*(loadstep+1))
    with TaskManager():
        NewtonMinimization (a, gfu, inverse="sparsecholesky", printing=False)
    sceneu.Redraw()

loadstep  0

loadstep  1

loadstep  2

loadstep  3
Contents