Naghdi shell

60. 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);

60.1. unknown field variables:#

  • displacement \(u \in [H^1(S)]^3\)

  • linearized rotation \(\beta \in H(\text{curl},S)\)

  • bending moments \(\sigma \in H(\text{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\):

\[ \int_S W(C_{tt}(u)) \]

Shear energy, \(\beta\) the linearized rotation in tangential plane:

\[ \int_S | n^T \nabla u - \beta |^2 \]

Bending energy (with TDNNS mixed method):

\[ t^2 \int_S | \varepsilon_{tt}(\beta) |^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