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.
[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.AddSurface(cyl, finitecyl)
geo.NameEdge(cyl,left, "left")
geo.NameEdge(cyl,right, "right")
mesh = Mesh(geo.GenerateMesh(maxh=0.2))
[2]:
mesh.nv
[2]:
84
[3]:
mesh.Curve(order)
Draw(mesh)
[3]:
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\).
[4]:
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)}\]
[5]:
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.
[6]:
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
[7]:
Draw(gfu.components[1], mesh, "rotations", deformation=gfu.components[0])
Draw(gfu.components[0], mesh, "disp")
[7]:
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].
[8]:
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)
[8]:
BaseWebGuiScene
To avoid membrane locking Regge interpolation as in [Neunteufel and Schöberl. Avoiding Membrane Locking with Regge Interpolation] can be used.
[9]:
# 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(fesVF, 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 )
---------------------------------------------------------------------------
NgException Traceback (most recent call last)
Cell In[9], line 44
40 nelphys = Cross(nphys,tphys)
42 Hn = CoefficientFunction( (u.Operator("hesseboundary").trans*nphys), dims=(3,3) )
---> 44 cfnphys = Normalize(Cof(Ptau+grad(solution.components[0]))*nsurf)
46 cfn = Normalize(CoefficientFunction( averednv.components ))
47 cfnR = Normalize(CoefficientFunction( averednv_start.components ))
NgException: Dimensions don't match, op = + dims1 = 0: 3
1: 3
, dims2 = 0: 3
\[\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)
[ ]:
Draw(solution.components[0], mesh, "disp", deformation=solution.components[0])
[ ]:
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)
[ ]:
Draw(solution.components[0], mesh, "disp", deformation=solution.components[0])
[ ]:
ei = ElementId(BND,2)
print ("len dofs", len(fes.GetDofNrs(ei)))
el = fes.GetFE(ei)
print(el.ndof)
[ ]:
a.AssembleLinearization(gfu.vec)
[14]:
fesA = bfA.space
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[14], line 1
----> 1 fesA = bfA.space
NameError: name 'bfA' is not defined
[ ]:
print ("len dofs", len(fesA.GetDofNrs(ei)))
el = fesA.GetFE(ei)
print(el.ndof)
[13]:
print (fesA.components)
for cfes in fesA.components:
print ("len dof", len(cfes.GetDofNrs(ei)))
el = cfes.GetFE(ei)
print (el.ndof)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[13], line 1
----> 1 print (fesA.components)
2 for cfes in fesA.components:
3 print ("len dof", len(cfes.GetDofNrs(ei)))
NameError: name 'fesA' is not defined
[ ]:
print (fesA.components[-1])
[12]:
mesh.GetBoundaries(), mesh.GetBBoundaries()
[12]:
(('default',), ('bottom', 'right', 'top', 'left'))
[ ]: