61. Fully nonlinear Koiter shell model

[M. Neunteufel and J.S.: The Hellan-Herrmann-Johnson method for nonlinear shells. Computers & Structures, ‘19

[M. Neunteufel and J.S.: Avoiding membrane locking with Regge interpolation, CMAME, ‘21

M. Neunteufel and J.S.: The Hellan-Herrmann-Johnson and TDNNS method for linear and nonlinear shells

Bending energy:

\[ \int | \nabla^2 u_i n_i |^2 \]

with physical normal vector \(n = \operatorname{Cof} (F) n_{ref}\). Discretized by the Hellan-Herrmann-Johnson method. Distributional jump term \([\partial_n w]\) is now angle between deformed faces.

Membrane energy with reduction operator for tangential strain:

\[ \int W( I_{cc} C_{tt}(u) ) \]

from ngsolve import *
from ngsolve.webgui import Draw
from time import sleep

thickness = 0.1   # thickness
L = 12
W = 1
E, nu = 1.2e6, 0

from netgen.occ import *
shape = Rectangle(L,W).Face()
shape.edges.Min(X).name="left"
shape.edges.Max(X).name="right"
shape.edges.Min(Y).name="bottom"
shape.edges.Max(Y).name="top"

mesh = Mesh(OCCGeometry(shape).GenerateMesh(maxh=0.5))
# Draw (mesh)

regge = True   # False -> membrane locking
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 = CF( (u.Operator("hesseboundary").trans*nphys), dims=(3,3) )

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

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

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.AssembleLinearization(averednv.vec)
invF = bfF.mat.Inverse(fesVF.FreeDofs(), inverse="sparsecholesky")

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

UpdateNormal()
averednv_start.vec.data = averednv.vec

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

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

bfA = BilinearForm(fes, symmetric=True, condense=True)

# bending energy
bfA += Variation(-6/thickness**3*MaterialInv(sigma, E, nu)*ds ).Compile()
bfA += Variation(InnerProduct(sigma, Hn + (1-nphys*nsurf)*gradn)*ds ).Compile()
bfA += Variation( -(acos(nel*cfnR)-acos(nelphys*pnaverage)-hyb*nel)*(sigma*nel)*nel*ds(element_boundary=True) ).Compile()

# membrane energy
if regge:
    bfA += Variation( 0.5*thickness*Material(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 )

Apply moment to right edge:

par = Parameter(0.0)
moment = IfPos(x-L+1e-6, 1, 0)*50*pi/3
bfA += Variation( -par*moment*(hyb*nel)*ds(element_boundary=True) )

numsteps=10
scene = Draw(solution.components[0], mesh, deformation=solution.components[0])

with TaskManager():
    for steps in range(0,numsteps):
        par.Set((steps+1)/numsteps)
        print("Loadstep =", steps+1, ", F/Fmax =", (steps+1)/numsteps*100, "%")
        
        UpdateNormal()
        
        solvers.Newton(bfA, solution, inverse="sparsecholesky", printing=False, maxerr=2e-10)
        scene.Redraw()
        sleep(0.5)

Loadstep = 1 , F/Fmax = 10.0 %

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 %