Sydney

7. Sydney#

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

Show code cell source

Hide code cell source

# first panel
v1l = Pnt(0, -10, 0.1)
v1t = Pnt(-6, 0, 10.5)
w1 = Wire(ArcOfCircle(v1l, Vec(-0.3, 0.4, 1), v1t))

v12l = Pnt(0.25, -10, 0.1)
v12t = Pnt(-0.5, 0, 11)
w12 = Wire(ArcOfCircle(v12l, Vec(-0.25, 0.3, 1), v12t))

v2l = Pnt(0.5, -10, 0.1)
v2t = Pnt(5 - 0.15, 0, 10)
w2 = Wire(ArcOfCircle(v2l, Vec(0, 0.2, 1), v2t))

# vent
v1lv = Pnt(1, -10.5, 1.5)
v1tv = Pnt(5, 0, 10)
w2v = Wire(ArcOfCircle(v1lv, Vec(0, 0.4, 1), v1tv))

v2lv = Pnt(4.5, -10.5, 0.5)
v2tv = v1tv
w3v = Wire(ArcOfCircle(v2lv, Vec(0, 0.2, 1), v2tv))

v3l = Pnt(5, -10, 1)
v3t = Pnt(5 + 0.15, 0, 10)
w3 = Wire(ArcOfCircle(v3l, Vec(0, 0, 1), v3t))

v4l = Pnt(5.5, -10 - 1, 1)
v4t = Pnt(3.5, 0, 15)
w4 = Wire(ArcOfCircle(v4l, Vec(0, 0, 1), v4t))

# second panel
v212l = Pnt(5.5 + 0.25, -10 - 1, 1)
v212t = Pnt(3.5 + 5.5, 0, 11 + 2)
w212 = Wire(ArcOfCircle(v212l, Vec(0, 0, 1), v212t))

v22l = Pnt(5.5 + 0.5, -10 - 1, 1)
v22t = Pnt(3.5 + 11 - 0.15, 0, 10 + 1)
w22 = Wire(ArcOfCircle(v22l, Vec(0.2, 0, 1), v22t))

# vent
v21lv = Pnt(5.5 + 0.5 + 0.5, -10 - 1.5, 1.5)
v21tv = Pnt(3.5 + 11, 0, 10 + 1)
w22v = Wire(ArcOfCircle(v21lv, Vec(0.3, 0.1, 1), v21tv))

v22lv = Pnt(3.5 + 11 - 1.5, -10 - 1.5, 1)
v22tv = v21tv
w23v = Wire(ArcOfCircle(v22lv, Vec(-0.2, 0, 1), v22tv))


v23l = Pnt(3.5 + 11, -10 - 1, 0)
v23t = Pnt(3.5 + 11 + 0.15, 0, 10 + 1)
w23 = Wire(ArcOfCircle(v23l, Vec(0, 0, 1), v23t))

v24l = Pnt(3.5 + 11, -10 - 1.5, 0)
v24t = Pnt(11, 0, 14 + 8)
w24 = Wire(ArcOfCircle(v24l, Vec(-0.1, 0, 1), v24t))

# third panel
v312l = Pnt(3.5 + 11 + 0.25, -10 - 1.5, 0)
v312t = Pnt(11 + 9, 0, 11 + 6.5)
w312 = Wire(ArcOfCircle(v312l, Vec(0.2, 0, 1), v312t))

v32l = Pnt(3.5 + 11 + 0.5, -10 - 1.5, 0)
v32t = Pnt(11 + 15 - 0.25, 0, 9)
w32 = Wire(ArcOfCircle(v32l, Vec(1, 0, 1), v32t))

# fouth panel

# vent
v32lv = Pnt(3.5 + 11 + 0.5 + 4, -10 - 1.5 - 0.5, 3)
v32tv = Pnt(11 + 15 - 0.25, 0, 9)
w32v = Wire(ArcOfCircle(v32lv, Vec(0.6, 0.3, 1), v32tv))

v33lv = Pnt(3.5 + 11 + 0.5 + 9, -10 - 1.5 - 0.5, 2)
v33tv = Pnt(11 + 15, 0, 9)
w33v = Wire(ArcOfCircle(v33lv, Vec(0, 0.3, 1), v33tv))

v34lv = Pnt(3.5 + 11 + 17 - 2, -10 - 1.5 - 0.5, 3)
v34tv = Pnt(11 + 15, 0, 9)
w34v = Wire(ArcOfCircle(v34lv, Vec(-0.6, 0.3, 1), v34tv))


v41l = Pnt(3.5 + 11 + 17, -10 - 1.5, 0)
v41t = Pnt(11 + 15 + 0.25, 0, 9)
w41 = Wire(ArcOfCircle(v41l, Vec(-0.5, 0, 1), v41t))

v42l = Pnt(3.5 + 11 + 17 + 0.25, -10 - 1.5, 0)
v42t = Pnt(3.5 + 11 + 17 + 0.25 - 0.5, 0, 11.75)
w42 = Wire(ArcOfCircle(v42l, Vec(0, 0, 1), v42t))

v43l = Pnt(3.5 + 11 + 17 + 0.5, -10 - 1.5, 0)
v43t = Pnt(3.5 + 11 + 17 + 0.5 + 5, 0, 12)
w43 = Wire(ArcOfCircle(v43l, Vec(0.25, 0, 1), v43t))

shape = Glue(
    [w1, w12, w2, w2v, w3v, w3, w4]
    + [w212, w22, w22v, w23v, w23, w24]
    + [w312, w32, w32v, w33v, w34v]
    + [w41, w42, w43]
)
Draw(shape);

(model created by Sebastian Hirnschall)

mesh = Mesh(OCCGeometry(shell).GenerateMesh(maxh=2))
mesh.RefineHP(1)
mesh.Curve(3)
Draw(mesh, **ea);

t = 0.1
order = 4

E = 2.85e4
nu = 0.3
mu = E / (2 * (1 + nu))
lam = E * nu / (1 - nu**2)

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 / (nu + 1) * Trace(mat) ** 2)

fesU = VectorH1(mesh, order=order, dirichlet_bbnd="fix")
fesS = HDivDivSurface(mesh, order=order - 1, discontinuous=True)
fesH = NormalFacetSurface(mesh, order=order - 1)  # , dirichlet_bbnd="fix")

fes = fesU * fesS * fesH
(u, sigma, uh) = fes.TrialFunction()
# Need to take traces as we are on the surface
sigma, uh = sigma.Trace(), uh.Trace()

fesRegge = HCurlCurl(mesh, order=order - 1)

# normal, tangential, and co-normal vectors
nv = specialcf.normal(3)
tv = specialcf.tangential(3)
cnv = Cross(nv, tv)

# Projection to tangent space and surface derivatives
Ptau = Id(3) - OuterProduct(nv, nv)
grad_u = Grad(u).Trace()
Ftau = grad_u + Ptau
Etautau = 0.5 * (Ftau.trans * Ftau - Ptau)  # Green strain

# sum_{i=1}^3 hesse(u_i) nv_i
Hn = (u.Operator("hesseboundary").trans * nv).Reshape((3, 3))

bfa = BilinearForm(fes, symmetric=True, condense=True)
# membrane energy
bfa += Variation(t/2 * MaterialNorm(Interpolate(Etautau, fesRegge), E, nu) * ds)

# bending energy
bfa += Variation(-12 / t**3 / 2 * MaterialNormInv(sigma, E, nu) * ds)
bfa += Variation(InnerProduct(sigma, Hn) * ds)
bfa += Variation(sigma[cnv, cnv] * (uh - grad_u[nv, :]) * cnv * ds(element_boundary=True))

force = CF((0, 0, -0.01))
bfa += Variation(-force * u * ds)

gfu = GridFunction(fes)

with TaskManager():
    solvers.Newton(a=bfa, u=gfu, maxerr=1e-8, printing=True, inverse="sparsecholesky")

Newton iteration  0

err =  0.8193186595548007
Newton iteration  1

err =  0.5023462805656246
Newton iteration  2

err =  0.09756148456378141
Newton iteration  3

err =  0.026967332912739853
Newton iteration  4

err =  0.0006861963404850162
Newton iteration  5

err =  7.158281072951311e-06
Newton iteration  6

err =  2.663310371778892e-08
Newton iteration  7

err =  2.572922628578534e-12

Draw(gfu.components[0], scale=20, deformation=True, **ea);

Draw(Norm(gfu.components[1]), mesh, order=3, max=0.1, **ea);