Nonlinear Koiter and Naghdi shells

17. Nonlinear Koiter and Naghdi shells#

The linear shell notebooks introduced Kirchhoff-Love and Reissner-Mindlin shells together with HHJ and TDNNS discretizations. We now keep the nonlinear strain measures. The structure remains familiar:

  • Koiter shells use membrane and bending energy, but no transverse shear.

  • Naghdi shells add a director/shear field and therefore a shear energy.

  • HHJ and TDNNS again avoid globally \(H^2\)-conforming displacement spaces by introducing bending moments and suitable trace variables.

17.1. Nonlinear shell energies#

Let \({\mathcal{S}}\) be the reference midsurface with unit normal \(\hat\nu\) and tangent projection

\[ P_{\mathcal{S}} = I-\hat\nu\otimes\hat\nu. \]

For a displacement \(u:{\mathcal{S}}\to\mathbb{R}^3\) we define the surface deformation gradient and surface Green strain on the reference midsurface by

\[ F_{\mathcal{S}} = P_{\mathcal{S}} + \nabla_{{\mathcal{S}}}u, \qquad E_{\mathcal{S}}(u)=\frac12( F_{\mathcal{S}}^T F_{\mathcal{S}}- P_{\mathcal{S}}). \]

The membrane energy is

\[ \mathcal{W}_{\rm mem}(u) =\frac{t}{2}\int_{{\mathcal{S}}}\| E_{\mathcal{S}}(u)\|_{\mathbb{M}}^2\,d s. \]

Here \(t\) is the thickness and \(\mathbb{M}\) is the plane-stress material tensor with Young modulus \(\bar E\) and Poisson ratio \(\bar\nu\).

membrane

bending

shear

17.2. Bending and shear#

The nonlinear bending energy measures the change of curvature between the reference and deformed midsurfaces. If \(\nu\circ\phi\) is the normal of the deformed shell pulled back to the reference surface, then a compact form is

\[ \mathcal{W}_{\rm bend}(u) =\frac{t^3}{24}\int_{{\mathcal{S}}} \left\| F_{\mathcal{S}}^T\nabla_{{\mathcal{S}}}(\nu\circ\phi)-\nabla_{{\mathcal{S}}}\hat\nu\right\|_{\mathbb{M}}^2\,d s. \]

Koiter shells use membrane plus bending:

\[ \mathcal{W}_{\rm Koiter}(u) =\mathcal{W}_{\rm mem}(u)+\mathcal{W}_{\rm bend}(u)-\int_{{\mathcal{S}}} f\cdot u\,d s. \]

Naghdi shells add a director/shear variable \(\gamma\). With a director \(\tilde\nu\circ\phi\), the shear energy has the form

\[ \mathcal{W}_{\rm shear}(u,\gamma) =\frac{t\kappa G}{2}\int_{{\mathcal{S}}}\| F_{\mathcal{S}}^T\tilde\nu\circ\phi\|^2\,d s, \]

where \(G=\bar E/(2(1+\bar\nu))\) and \(\kappa=5/6\).

17.3. Nonlinear Koiter shell with HHJ#

For affine or piecewise flat surface meshes, curvature information is distributional and lives on edges. We will discuss generalized curvatures tomorrow. The nonlinear HHJ formulation captures this by combining an element curvature term with an edge angle term. The bending moment tensor \(\sigma\) plays the same role as in the linear Kirchhoff-Love shell, but the curvature measure is nonlinear. We start with the following three-field formulation [Neunteufel, Schöberl. The Hellan-Herrmann-Johnson and TDNNS methods for linear and nonlinear shells. Computers & Structures, 305 (2024).] to incorporate the distributional extrinsic curvature difference of the initial and deformed shell configuration:

Find \((u,{\boldsymbol{\kappa}}^{\mathrm{diff}},\boldsymbol{\sigma})\in V_h^k\times M_h^{k-1,\mathrm{dc}}\times M_h^{k-1}\) for the Lagrangian

\[\begin{split} \begin{align*} \mathcal{L}(u,{\boldsymbol{\kappa}}^{\mathrm{diff}},\boldsymbol{\sigma})&=\int_{{\mathcal{S}}}\Big(\frac{t}{2}\|{\boldsymbol{E}}_{\mathcal{S}}(u)\|_{\mathbb{M}}^2+\frac{t^3}{12}\|{\boldsymbol{\kappa}}^{\mathrm{diff}}\|_{\mathbb{M}}^2\Big)\,d s+\sum_{ T\in{\mathcal{T}}}\int_{ T}\big({\boldsymbol{\kappa}}^{\mathrm{diff}}-({\boldsymbol{F}}_{\mathcal{S}}^T\nabla_{{\mathcal{S}}}(\nu\circ\phi)-\nabla_{{\mathcal{S}}}\hat{\nu}) \big):\boldsymbol{\sigma}\,d s \\ &\qquad+\sum_{ E\in{\mathcal{E}}}\int_{ E}(\sphericalangle(\nu_L\circ\phi,\nu_R\circ\phi)-\sphericalangle(\hat{\nu}_L,\hat{\nu}_R))\boldsymbol{\sigma}_{\hat{\mu}\hat{\mu}}\,d l-\int_{{\mathcal{S}}} f\cdot u\,d s. \end{align*} \end{split}\]

17.3.1. Cylindrical shell under volume load#

We start with a cylindrical shell, which is clamped at the left side and free at the right side, under a gradually increased load.

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
import ipywidgets as widgets

order = 2

cyl = Cylinder((0, 0, 0), (1, 0, 0), 0.4, 1).faces[0]
cyl.edges.Min(X).name = "left"
cyl.edges.Min(X).hpref = 1
cyl.edges.Max(X).name = "right"
mesh = Mesh(OCCGeometry(cyl).GenerateMesh(maxh=0.1)).Curve(order)
mesh.RefineHP(2, 0.3)
mesh.Curve(order)

young_modulus, poisson = 2e1, 0.1
thickness = 0.02

Draw(mesh, euler_angles=[0, -50, -10]);

At the clamped boundary we fix the displacement. In the mixed formulation, essential and natural boundary conditions swap for the moment tensor. Thus we impose \(\sigma_{\hat\mu\hat\mu}=0\) at the free boundary, meaning that no bending moment is prescribed there.

fes_sigma = HDivDivSurface(mesh, order=order - 1, dirichlet_bbnd="right")
fes_u = VectorH1(mesh, order=order, dirichlet_bbnd="left")
fes = fes_u * fes_sigma
u, sigma = fes.TrialFunction()
sigma = sigma.Trace()

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

gf_solution = GridFunction(fes, name="solution")

We use \(\hat\nu\), \(\hat t\), and \(\hat\mu=\hat\nu\times\hat t\) for the reference normal, edge tangent, and conormal. Their deformed counterparts depend on the unknown displacement \(u\).

../_images/nv_conv_tang_trig.png 
normal = specialcf.normal(mesh.dim)
tangent = specialcf.tangential(mesh.dim)
conormal = Cross(normal, tangent)

Ptau = Id(mesh.dim) - OuterProduct(normal, normal)
Ftau = Grad(u).Trace() + Ptau
Etautau = 0.5 * (Ftau.trans * Ftau - Ptau)

normal_def = Normalize(Cof(Ftau) * normal)
tangent_def = Normalize(Ftau * tangent)
conormal_def = Cross(normal_def, tangent_def)

normal_hessian_def = (u.Operator("hesseboundary").trans * normal_def).Reshape((3, 3))

For the edge angle term we use an averaged normal vector on facets. This avoids accessing both neighboring element normals at once and gives a more stable formulation using conormals.

../_images/nonsmooth_av_nv_el_nv.png
\[ \begin{align*} \sum_{ E\in{\mathcal{E}}}\int_{ E}(\sphericalangle(\nu_L,\nu_R)\circ\phi-\sphericalangle(\hat{\nu}_L,\hat{\nu}_R))\boldsymbol{\sigma}_{\hat{\mu}\hat{\mu}}\,d l &= \sum_{ T\in{\mathcal{T}}}\int_{\partial T}(\sphericalangle(\mu\circ\phi,P^{\perp}_{\tau\circ\phi}(\{\!\{\nu^n\}\!\}))-\sphericalangle(\hat{\mu},\{\!\{\hat{\nu}\}\!\}))\boldsymbol{\sigma}_{\hat{\mu}\hat{\mu}}\,d l, \end{align*} \]

where

\[ \begin{align*} P^{\perp}_{\tau\circ\phi}(v)= \frac{1}{\|\boldsymbol{P}^{\perp}_{\tau\circ\phi}v\|}\boldsymbol{P}^{\perp}_{\tau\circ\phi}v,\qquad \boldsymbol{P}^{\perp}_{\tau\circ\phi}= \boldsymbol{I}-\tau\circ\phi\otimes\tau\circ\phi \end{align*} \]

denotes the (nonlinear) normalized projection to the plane perpendicular to the deformed edge tangent vector \(\tau\) for measuring the correct angle.

gf_clamped_bnd = GridFunction(FacetSurface(mesh, order=0))
gf_clamped_bnd.Set(1, definedon=mesh.BBoundaries("left"))

cf_normal_cur = Normalize(Cof(Ptau + Grad(gf_solution.components[0])) * normal)

fes_facet_vector = VectorFacetSurface(mesh, order=order)
averaged_normal = GridFunction(fes_facet_vector)
averaged_normal_init = GridFunction(fes_facet_vector)

averaged_normal.Set(normal, dual=True, definedon=mesh.Boundaries(".*"))
averaged_normal_init.Set(normal, dual=True, definedon=mesh.Boundaries(".*"))

cf_averaged_normal_init = Normalize(averaged_normal_init)
cf_projected_averaged_normal = Normalize(
    averaged_normal - (tangent_def * averaged_normal) * tangent_def
)

Define the material and inverse material norms \(\|\cdot\|_{\mathbb{M}}^2\) and \(\|\cdot\|_{\mathbb{M}^{-1}}^2\) with Young modulus \(\bar{E}\) and Poisson’s ratio \(\bar{\nu}\)

\[ \begin{align*} \mathbb{M} {\boldsymbol{E}} = \frac{\bar E}{1-\bar \nu^2}\big((1-\bar \nu){\boldsymbol{E}}+\bar \nu\,\mathrm{tr}({\boldsymbol{E}}) P_{\mathcal{S}}\big),\qquad\mathbb{M}^{-1} \boldsymbol{\sigma} = \frac{1+\bar \nu}{\bar E}\big(\boldsymbol{\sigma}-\frac{\bar \nu}{\bar\nu+1}\,\mathrm{tr}(\boldsymbol{\sigma}) P_{\mathcal{S}}\big). \end{align*} \]
def MaterialNorm(mat, young_modulus, poisson):
    return young_modulus / (1 - poisson**2) * (
        (1 - poisson) * InnerProduct(mat, mat) + poisson * Trace(mat) ** 2
    )


def MaterialNormInv(mat, young_modulus, poisson):
    return (1 + poisson) / young_modulus * (
        InnerProduct(mat, mat) - poisson / (poisson + 1) * Trace(mat) ** 2
    )

Define the bilinear form for the problem including membrane and bending energy

\[\begin{split} \begin{align*} \mathcal{L}(u,\sigma) &=\int_{{\mathcal{S}}}\Big(\frac{t}{2}\|{\boldsymbol{E}}_{\mathcal{S}}(u)\|^2_{\mathbb{M}} -\frac{6}{t^3}\|\boldsymbol{\sigma}\|^2_{\mathbb{M}^{-1}}\Big)\,d s + \sum_{ T\in{\mathcal{T}}}\Big(\int_{ T} \boldsymbol{\sigma}:(\boldsymbol{H}_{\nu\circ\phi}+(1-\hat{\nu}\cdot\nu\circ\phi)\nabla_{{\mathcal{S}}}\hat{\nu})\,d s \\ &\qquad+ \int_{\partial T}(\sphericalangle(\mu\circ\phi,\{\!\{\nu^n\}\!\})-\sphericalangle(\hat{\mu},\{\!\{\hat{\nu}\}\!\}))\boldsymbol{\sigma}_{\hat{\mu}\hat{\mu}}\,d l\Big) - \int_{{\mathcal{S}}} f\cdot u\,d s. \end{align*} \end{split}\]

For \(k\geq 2\) so-called membrane locking can occur for small thickness parameters due to the different scaling of the membrane and bending energy. To circumvent this locking problem we can interpolate the membrane strains into Regge finite elements of one order less than the displacement fields [Neunteufel, Schöberl. Avoiding membrane locking with Regge interpolation. Computer Methods in Applied Mechanics and Engineering, (2021).].

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

# Membrane energy. Regge interpolation weakens membrane locking.
bfa += Variation(
    0.5
    * thickness
    * MaterialNorm(Interpolate(Etautau, fes_regge), young_modulus, poisson)
    * ds
).Compile()

# Nonlinear HHJ bending energy.
bfa += Variation(
    (
        -6 / thickness**3 * MaterialNormInv(sigma, young_modulus, poisson)
        + InnerProduct(
            normal_hessian_def + (1 - normal * normal_def) * Grad(normal), sigma
        )
    )
    * ds
).Compile()
bfa += Variation(
    (
        acos(conormal_def * cf_projected_averaged_normal)
        - acos(conormal * cf_averaged_normal_init)
    )
    * sigma[conormal, conormal]
    * ds(element_boundary=True)
).Compile()

par = Parameter(0.0)
bfa += Variation(-thickness * par * 2 * y * u[1] * ds)

gf_solution.vec[:] = 0
gf_history = GridFunction(fes, multidim=0)
gf_u, gf_sigma = gf_solution.components

scene_u = Draw(gf_u, mesh, "displacement", deformation=gf_u, euler_angles=[0, -50, -10])
scene_sigma = Draw(
    Norm(gf_sigma), mesh, "bending_stress", deformation=gf_u, euler_angles=[0, -50, -10]
);

num_steps = 20

tw = widgets.Text(value="step = 0")
display(tw)

with TaskManager():
    for step in range(num_steps):
        par.Set((step + 1) / num_steps)

        averaged_normal.Set(
            (1 - gf_clamped_bnd) * cf_normal_cur + gf_clamped_bnd * normal,
            dual=True,
            definedon=mesh.Boundaries(".*"),
        )

        solvers.Newton(
            bfa,
            gf_solution,
            inverse="",
            printing=False,
            maxerr=1e-5,
        )
        scene_u.Redraw()
        scene_sigma.Redraw()
        gf_history.AddMultiDimComponent(gf_solution.vec)
        tw.value = f"step = {step + 1}/{num_steps}"

Warning: Newton might not converge! Error =  2150.814077266169

Warning: Newton might not converge! Error =  35.60543203465865

Warning: Newton might not converge! Error =  274.34118036908194

Warning: Newton might not converge! Error =  40268.82408729595

Warning: Newton might not converge! Error =  195.43990791716794

Warning: Newton might not converge! Error =  42.619834689204225

Draw(
    gf_history.components[0],
    mesh,
    animate=True,
    min=0,
    max=0.25,
    autoscale=True,
    deformation=True,
    euler_angles=[0, -50, -10],
);

17.4. Hybridized Koiter shell: cantilever with bending moment#

The direct mixed HHJ formulation is a saddle point problem. Hybridization makes \(\sigma\) discontinuous and restores normal-normal continuity by a Lagrange multiplier \(\alpha\). This also makes it easy to prescribe bending moments weakly.

../_images/cant_bend_mom_1d.png

We consider a strip clamped on the left and load it with a bending moment on the right. For Poisson ratio \(\bar\nu=0\), the strip should roll up approximately to a circle.

thickness = 0.1
length = 12
width = 1
young_modulus, poisson = 1.2e6, 0
moment = IfPos(x - length + 1e-6, 1, 0) * 50 * pi / 3

rect = Rectangle(length, width).Face()
rect.edges.Min(X).name = "left"
rect.edges.Max(X).name = "right"
rect.edges.Min(Y).name = "bottom"
rect.edges.Max(Y).name = "top"
mesh = Mesh(OCCGeometry(rect).GenerateMesh(maxh=1))
Draw(mesh);

order = 2

fes_sigma = HDivDivSurface(mesh, order=order - 1, discontinuous=True)
fes_u = VectorH1(
    mesh,
    order=order,
    dirichletx_bbnd="left",
    dirichlety_bbnd="left|bottom|top",
    dirichletz_bbnd="left",
)
fes_hybrid = NormalFacetSurface(mesh, order=order - 1, dirichlet_bbnd="left|bottom|top")
fes = fes_u * fes_sigma * fes_hybrid
u, sigma, alpha = fes.TrialFunction()
sigma, alpha = sigma.Trace(), alpha.Trace()

fes_regge = HCurlCurl(mesh, order=order - 1)
gf_solution = GridFunction(fes, name="solution")

Ptau = Id(mesh.dim) - OuterProduct(normal, normal)
Ftau = Grad(u).Trace() + Ptau
Etautau = 0.5 * (Ftau.trans * Ftau - Ptau)

normal_def = Normalize(Cof(Ftau) * normal)
tangent_def = Normalize(Ftau * tangent)
conormal_def = Cross(normal_def, tangent_def)
normal_hessian_def = (u.Operator("hesseboundary").trans * normal_def).Reshape((3, 3))

gf_clamped_bnd = GridFunction(FacetSurface(mesh, order=0))
gf_clamped_bnd.Set(1, definedon=mesh.BBoundaries("left"))

cf_normal_cur = Normalize(Cof(Ptau + Grad(gf_solution.components[0])) * normal)

fes_facet_vector = VectorFacetSurface(mesh, order=order)
averaged_normal = GridFunction(fes_facet_vector)
averaged_normal_init = GridFunction(fes_facet_vector)

averaged_normal.Set(normal, dual=True, definedon=mesh.Boundaries(".*"))
averaged_normal_init.Set(normal, dual=True, definedon=mesh.Boundaries(".*"))
cf_averaged_normal_init = Normalize(averaged_normal_init)
cf_projected_averaged_normal = Normalize(
    averaged_normal - (tangent_def * averaged_normal) * tangent_def
)

bfa = BilinearForm(fes, symmetric=True, condense=True)
bfa += Variation(
    0.5
    * thickness
    * MaterialNorm(Interpolate(Etautau, fes_regge), young_modulus, poisson)
    * ds
)
bfa += Variation(
    (
        -6 / thickness**3 * MaterialNormInv(sigma, young_modulus, poisson)
        + InnerProduct(
            normal_hessian_def + (1 - normal * normal_def) * Grad(normal), sigma
        )
    )
    * ds
).Compile()
bfa += Variation(
    (
        acos(conormal_def * cf_projected_averaged_normal)
        - acos(conormal * cf_averaged_normal_init)
        + alpha * conormal
    )
    * sigma[conormal, conormal]
    * ds(element_boundary=True)
).Compile()

par = Parameter(0.0)
bfa += Variation(-par * moment * (alpha * conormal) * ds(element_boundary=True))

gf_solution.vec[:] = 0
gf_history = GridFunction(fes, multidim=0)
gf_u, gf_sigma, _ = gf_solution.components
scene_u = Draw(
    gf_u,
    mesh,
    "displacement",
    deformation=gf_u,
    euler_angles=[-20, 0, 0],
);

point_P = (length, width / 2, 0)
result_P = [(0, 0, 0)]
exact_solution = [(0, 0, 0)]


def GetExactSolution(par):
    radius = 100 / (par.Get() * 50 * pi / 3)
    return (radius * sin(length / radius) - length, 0, -radius * cos(length / radius) + radius)


tw = widgets.Text(value="step = 0")
display(tw)
num_steps = 20

with TaskManager():
    for step in range(num_steps):
        par.Set((step + 1) / num_steps)

        averaged_normal.Set(
            (1 - gf_clamped_bnd) * cf_normal_cur + gf_clamped_bnd * normal,
            dual=True,
            definedon=mesh.Boundaries(".*"),
        )

        solvers.Newton(
            bfa,
            gf_solution,
            inverse="sparsecholesky",
            printing=False,
            maxerr=1e-5,
        )
        scene_u.Redraw()

        result_P.append((gf_u(mesh(*point_P, BND))))
        exact_solution.append(GetExactSolution(par))
        gf_history.AddMultiDimComponent(gf_solution.vec)
        tw.value = f"step = {step + 1}/{num_steps}"

Draw(
    gf_history.components[0],
    mesh,
    animate=True,
    min=0,
    max=12,
    autoscale=True,
    deformation=True,
    euler_angles=[-20, 0, 0],
);

import matplotlib.pyplot as plt

u_x, _, u_z = zip(*result_P)
y_axis = [i / num_steps for i in range(len(u_x))]
u_x = [-val for val in u_x]
u_ex_x, _, u_ex_z = zip(*exact_solution)
u_ex_x = [-val for val in u_ex_x]

plt.plot(u_x, y_axis, "-*", label="$-u_x$")
plt.plot(u_z, y_axis, "-x", label="$u_z$")
plt.plot(u_ex_x, y_axis, "--", color="k", label="$-u_{ex,x}$")
plt.plot(u_ex_z, y_axis, "-.", color="k", label="$u_{ex,z}$")
plt.xlabel("displacement")
plt.ylabel(r"$M/M_{\max}$")
plt.legend()
plt.show()
../_images/88e21eae937159bd765f8d98d15a19916bd2d0c82c66eb7306c0d8aced5b6748.png

We can further increase the bending moment to see how many turns the strip performs.

scene_u = Draw(gf_u, mesh, "displacement", deformation=gf_u)
display(tw)

num_steps = 10
with TaskManager():
    for step in range(num_steps):
        par.Set(par.Get() + 1 / num_steps)

        averaged_normal.Set(
            (1 - gf_clamped_bnd) * cf_normal_cur + gf_clamped_bnd * normal,
            dual=True,
            definedon=mesh.Boundaries(".*"),
        )
        solvers.Newton(
            bfa,
            gf_solution,
            inverse="sparsecholesky",
            printing=False,
            maxerr=1e-5,
        )
        scene_u.Redraw()
        tw.value = f"step = {step + 1}/{num_steps}"

17.5. Nonlinear Naghdi shells with TDNNS#

Naghdi shells add shear to the Koiter model. We use a hierarchical shear variable \(\hat\gamma\) on top of the HHJ formulation. The director is

\[ \tilde\nu = \nu + F_{\mathcal{S}}^\dagger\hat\gamma, \]

where \( F_{\mathcal{S}}^\dagger\) is the Moore-Penrose pseudoinverse. The nonlinear TDNNS Lagrangian contains membrane, bending, and shear contributions and reads [Neunteufel, Schöberl. The Hellan-Herrmann-Johnson and TDNNS methods for linear and nonlinear shells. Computers & Structures, 305 (2024).]

\[\begin{align*} \mathcal{L}(u,\sigma,\hat\gamma)&= \int_{{\mathcal{S}}}\Big(\frac{t}{2}\|{\boldsymbol{E}}_{\mathcal{S}}(u)\|_{\mathbb{M}}^2-\frac{6}{t^3}\|\sigma\|^2_{\mathbb{C}^{-1}}+\frac{t\kappa G}{2}\|\hat\gamma\|^2\Big)\,d s+ \sum_{ T\in{\mathcal{T}}}\int_{ T}(\mathcal{H}_{\tilde{\nu}}+(1-\hat{\nu}\cdot\tilde{\nu})\nabla_{{\mathcal{S}}}\hat{\nu}-\nabla_{{\mathcal{S}}}\hat\gamma):\sigma\,d s \\ &\qquad+ \sum_{ E\in{\mathcal{E}}} \int_{ E} (\sphericalangle(\nu_L,\nu_R)\circ\phi-\sphericalangle(\hat{\nu}_L,\hat{\nu}_R)-[\![\hat\gamma_{\hat{\mu}}]\!])\sigma_{\hat{\mu}\hat{\mu}}\,d l. \end{align*}\]

17.6. T-structure#

thickness = 0.1
length = 1
young_modulus, poisson = 6.2e6, 0
G = young_modulus / (2 * (1 + poisson))
kappa = 5 / 6

shear = 3e3 * CF((1, 0, 1))

f1 = WorkPlane(Axes((0, 0, 0), n=X, h=Y)).Rectangle(length, length).Face()
f1.edges.Min(Z).name = "bottom"
f1.edges.Max(Z).name = "top"
f1.edges.Min(Y).name = "left"
f1.edges.Max(Y).name = "right"
f2 = WorkPlane(Axes((-length / 2, 0, length), n=Z, h=X)).Rectangle(length, length).Face()
f2.edges.Min(Y).name = "zpbottom"
f2.edges.Max(Y).name = "uptop"
f2.edges.Min(X).name = "upleft"
f2.edges.Max(X).name = "upright"
shape = Glue([f1, f2])
mesh = Mesh(OCCGeometry(shape).GenerateMesh(maxh=0.25))
Draw(mesh, euler_angles=[-60, 0, -15]);

order = 3

fes_u = VectorH1(mesh, order=order, dirichlet_bbnd="bottom")
fes_sigma = HDivDivSurface(mesh, order=order - 1, discontinuous=True)
fes_hybrid = NormalFacetSurface(mesh, order=order - 1, dirichlet_bbnd="bottom")
fes_gamma = HCurl(mesh, order=order - 1, dirichlet_bbnd="bottom")
fes = fes_u * fes_sigma * fes_gamma * fes_hybrid

u, sigma, gamma, alpha = fes.TrialFunction()
sigma, gamma, alpha = sigma.Trace(), gamma.Trace(), alpha.Trace()

fes_regge = HCurlCurl(mesh, order=order - 1, discontinuous=True)
gf_solution = GridFunction(fes, name="solution")

Ptau = Id(mesh.dim) - OuterProduct(normal, normal)
Ftau = Grad(u).Trace() + Ptau
Etautau = 0.5 * (Ftau.trans * Ftau - Ptau)


def PseudoInverse(mat, kernel_vector):
    """Pseudoinverse of a rank-(n-1) matrix.

    The vector kernel_vector must lie in the kernel of mat.
    """
    return Inv(mat.trans * mat + OuterProduct(kernel_vector, kernel_vector)) * mat.trans


normal_def = Normalize(Cof(Ftau) * normal)
tangent_def = Normalize(Ftau * tangent)
conormal_def = Cross(normal_def, tangent_def)
director = normal_def + PseudoInverse(Ftau, normal).trans * gamma

normal_hessian_def = (u.Operator("hesseboundary").trans * director).Reshape((3, 3))

gf_clamped_bnd = GridFunction(FacetSurface(mesh, order=0))
gf_clamped_bnd.Set(1, definedon=mesh.BBoundaries("bottom"))

cf_normal_cur = Normalize(Cof(Ptau + Grad(gf_solution.components[0])) * normal)

fes_facet_vector = VectorFacetSurface(mesh, order=order)
averaged_normal = GridFunction(fes_facet_vector)
averaged_normal_init = GridFunction(fes_facet_vector)

averaged_normal.Set(normal, dual=True, definedon=mesh.Boundaries(".*"))
averaged_normal_init.Set(normal, dual=True, definedon=mesh.Boundaries(".*"))
cf_averaged_normal_init = Normalize(averaged_normal_init)
cf_projected_averaged_normal = Normalize(
    averaged_normal - (tangent_def * averaged_normal) * tangent_def
)

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

# Membrane energy.
bfa += Variation(
    0.5
    * thickness
    * MaterialNorm(Interpolate(Etautau, fes_regge), young_modulus, poisson)
    * ds
).Compile()

# Bending energy.
bfa += Variation(
    (
        -6 / thickness**3 * MaterialNormInv(sigma, young_modulus, poisson)
        + InnerProduct(
            normal_hessian_def
            + (1 - normal * director) * Grad(normal)
            - Grad(gamma),
            sigma,
        )
    )
    * ds
).Compile()
bfa += Variation(
    (
        acos(conormal_def * cf_projected_averaged_normal)
        - acos(conormal * cf_averaged_normal_init)
        + alpha * conormal
        + (PseudoInverse(Ftau, normal).trans * gamma) * conormal_def
    )
    * sigma[conormal, conormal]
    * ds(element_boundary=True)
).Compile()

# Shear energy.
bfa += Variation(0.5 * thickness * kappa * G * gamma * gamma * ds)

par = Parameter(0.0)
bfa += Variation(-par * shear * u * dx(definedon=mesh.BBoundaries("upleft")))

gf_solution.vec[:] = 0
gf_history = GridFunction(fes, multidim=0)

gf_u, gf_sigma, gf_gamma, _ = gf_solution.components

scene = Draw(gf_u, mesh, "displacement", deformation=gf_u, euler_angles=[-60, 0, -10])
scene2 = Draw(
    Norm(gf_sigma),
    mesh,
    "bending_stress",
    deformation=gf_u,
    euler_angles=[-60, 0, -10],
)
scene3 = Draw(
    thickness * kappa * G * gf_gamma,
    mesh,
    "shear_stress",
    deformation=gf_u,
    euler_angles=[-60, 0, -10],
)

num_steps = 20

tw = widgets.Text(value="step = 0")
display(tw)

with TaskManager():
    for step in range(num_steps):
        par.Set((step + 1) / num_steps)

        averaged_normal.Set(
            (1 - gf_clamped_bnd) * cf_normal_cur + gf_clamped_bnd * normal,
            dual=True,
            definedon=mesh.Boundaries(".*"),
        )

        solvers.Newton(
            bfa,
            gf_solution,
            inverse="sparsecholesky",
            printing=False,
            maxerr=1e-5,
        )
        scene.Redraw()
        scene2.Redraw()
        scene3.Redraw()

        gf_history.AddMultiDimComponent(gf_solution.vec)
        tw.value = f"step = {step + 1}/{num_steps}"

Draw(
    gf_history.components[0],
    mesh,
    animate=True,
    min=0,
    max=1.25,
    autoscale=True,
    deformation=True,
    euler_angles=[-60, 0, -10],
);

17.7. Möbius strip EVP with a linear Kirchhoff-Love shell#

We can transfer the HHJ eigenvalue pattern from the flat Kirchhoff-Love plate to a curved surface. The unknown is now the vector-valued shell displacement, the normal facet variable stores the normal derivative trace, and the discontinuous \(H(\mathrm{div}\,\mathrm{div})\) variable stores the bending moment tensor.

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
import ipywidgets as widgets
import numpy as np
import math


def mobius_point(u, v, R=10.0):
    x = (R + v * math.cos(u / 2.0)) * math.cos(u)
    y = (R + v * math.cos(u / 2.0)) * math.sin(u)
    z = v * math.sin(u / 2.0)
    return gp_Pnt(x, y, z)


num_u, num_v = 10, 10
pts = np.zeros((num_u + 1, num_v + 1, 3), dtype=np.float64)

for i in range(num_u + 1):
    upar = math.pi * i / num_u
    for j in range(num_v + 1):
        vpar = j / num_v - 0.5
        p = mobius_point(upar, vpar, 3)
        pts[i, j, 0] = p[0]
        pts[i, j, 1] = p[1]
        pts[i, j, 2] = p[2]

face1 = SplineSurfaceApproximation(pts)

for i in range(num_u + 1):
    upar = math.pi + math.pi * i / num_u
    for j in range(num_v + 1):
        vpar = j / num_v - 0.5
        p = mobius_point(upar, vpar, 3)
        pts[i, j, 0] = p[0]
        pts[i, j, 1] = p[1]
        pts[i, j, 2] = p[2]

face2 = SplineSurfaceApproximation(pts)
face = Glue([face1, face2])
face.vertices.Min(X).name = "pnt"

mesh_mobius = face.GenerateMesh(maxh=0.25)
mesh_mobius.Curve(3)
Draw(mesh_mobius, euler_angles=[-70, 0, -25]);

import ngsolve.solvers


def SolveMobiusShellEVP(mesh, order=2, num=20, maxit=40):
    fes_u = VectorH1(mesh, order=order, dirichlet_bbbnd="")
    fes_normal_trace = NormalFacetSurface(mesh, order=order - 1)
    fes_sigma = PrivateSpace(
        HDivDivSurface(mesh, order=order - 1, discontinuous=True)
    )
    fes = fes_u * fes_normal_trace * fes_sigma

    (u, normal_trace, sigma), (du, dnormal_trace, dsigma) = fes.TnT()
    sigma, dsigma = sigma.Trace(), dsigma.Trace()
    normal_trace = normal_trace.Trace()
    dnormal_trace = dnormal_trace.Trace()

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

    normal = specialcf.normal(3)
    tangent = specialcf.tangential(3)
    conormal = Cross(normal, tangent)
    Ptau = Id(3) - OuterProduct(normal, normal)

    grad_u = Grad(u).Trace()
    grad_du = Grad(du).Trace()
    eps_u = 0.5 * (Ptau * grad_u + grad_u.trans * Ptau)
    eps_du = 0.5 * (Ptau * grad_du + grad_du.trans * Ptau)

    hesse_u = (u.Operator("hesseboundary").trans * normal).Reshape((3, 3))
    hesse_du = (du.Operator("hesseboundary").trans * normal).Reshape((3, 3))

    thickness = 1
    young_modulus = 1.0
    poisson = 0.3
    shift = 1e-3

    def material_stress(mat):
        return young_modulus / (1 - poisson**2) * (
            (1 - poisson) * mat + poisson * Trace(mat) * Ptau
        )

    def material_stress_inv(mat):
        return (1 + poisson) / young_modulus * (
            mat - poisson / (poisson + 1) * Trace(mat) * Ptau
        )

    stiffness = BilinearForm(fes, symmetric=True, condense=True)
    stiffness += (
        thickness
        * InnerProduct(Interpolate(material_stress(eps_u), fes_regge), Interpolate(eps_du, fes_regge))
        * ds
    )
    stiffness += (
        -12.0
        / thickness**3
        * InnerProduct(material_stress_inv(sigma), dsigma)
        * ds
    )
    stiffness += (InnerProduct(sigma, hesse_du) + InnerProduct(dsigma, hesse_u)) * ds
    stiffness += (
        sigma[conormal, conormal]
        * (dnormal_trace - grad_du[normal, :])
        * conormal
        + dsigma[conormal, conormal]
        * (normal_trace - grad_u[normal, :])
        * conormal
    ) * ds(element_boundary=True)
    stiffness += shift * InnerProduct(u, du) * ds
    stiffness.Assemble()

    mass = BilinearForm(InnerProduct(u, du) * ds, symmetric=True).Assemble()
    preconditioner = stiffness.mat.Inverse(
        freedofs=fes.FreeDofs(True), inverse="sparsecholesky"
    )

    evals, evecs = ngsolve.solvers.LOBPCG(
        stiffness.mat,
        mass.mat,
        pre=preconditioner,
        num=num,
        maxit=maxit,
        printrates=False,
    )

    modes = GridFunction(fes, multidim=0)
    for evec in evecs[6:]:
        modes.AddMultiDimComponent(evec)

    return fes, evals, modes


fes_mobius, evals_mobius, modes_mobius = SolveMobiusShellEVP(
    mesh_mobius, order=2, num=30, maxit=30
)
print("first computed eigenvalues:")
for i, lam in enumerate(evals_mobius[:10]):
    print(f"  {i:2d}: {lam:.6e}")

Draw(
    modes_mobius.components[0],
    mesh_mobius,
    name="Mobius shell eigenmodes",
    min=-0.1,
    max=0.35,
    deformation=True,
    scale=1,
    animate=True,
    euler_angles=[-70, 0, -25],
);

first computed eigenvalues:
   0: 1.000000e-03
   1: 1.000000e-03
   2: 1.000000e-03
   3: 1.000000e-03
   4: 1.000000e-03
   5: 1.000000e-03
   6: 7.839178e-03
   7: 7.839894e-03
   8: 8.330478e-03
   9: 8.331817e-03

17.8. Origami - pure plate bending#

To bend plate structures at a given folding line, we need to enforce that the membrane energy is (close to) zero. Using the Regge interpolation to avoid membrane locking is crucial.

Reference: [Bartels, Bonito, Hornung, Neunteufel, Babuska’s paradox in a nonlinear bending-folding model. Interfaces and Free Boundaries (accepted).]

../_images/wrinkling_experiment.png
../_images/results_h0.05_p3_hp1_vis2.png 
from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw
from ngsolve.solvers import Newton

thickness = sqrt(1e-6)
E, nu = 1e2, 0

order = 3

force_len = 0.2
force_ind = IfPos(force_len - x, 1, 0)
uD = IfPos(1 - x, 1 - x, 0) * 1 * CF((0, -y * 0.1, 0))

y_coo = -0.23
x_coo = 1
p1 = (sqrt(x_coo**2 - 0.5**2), -0.5, 0)
p2 = (sqrt(x_coo**2 - y_coo**2), y_coo, 0)
p3 = (sqrt(x_coo**2 - y_coo**2), 0, 0)
L = sqrt((p2[0] - p1[0]) ** 2 + (p2[1] - p1[1]) ** 2 + (p2[2] - p1[2]) ** 2)
angle1 = acos((p2[0] - p1[0]) / L) / pi * 180

V1 = Vertex(((p1[0] + p2[0]) / 2, (p1[1] + p2[1]) / 2, 0))
V1.name = "P"
V2 = Vertex(((p2[0] + p3[0]) / 2, (p2[1] + p3[1]) / 2, 0))
V2.name = "P2"
V3 = Vertex(p3)
V3.name = "P3"

rec = MoveTo(0, -0.5).Rectangle(2, 0.5).Face()

domain_left = (
    WorkPlane()
    .MoveTo(0, -0.5)
    .Line(p1[0])
    .Rotate(angle1)
    .Line(L)
    .Rotate(90 - angle1)
    .Line(-y_coo)
    .Rotate(90)
    .Line(p2[0])
    .Close()
    .Face()
)

rec2 = MoveTo(0, -0.5).Rectangle(force_len, 0.5).Face()
rec2.faces.name = "left_domain"
rec2.edges.Min(Y).name = "force"
rec2.edges.Max(Y).name = "sym"
rec2.edges.Max(X).name = "interface"
rec2.edges.Min(X).name = "left"

domain_left.faces.name = "left_domain"
domain_left.edges.name = "kink"
domain_left.edges.Min(X).name = "left"
domain_left.edges.Max(Y).name = "sym"
domain_left.edges.Min(Y).name = "bot"
v = domain_left.vertices.Nearest(p2)

domain_right = Glue([rec - domain_left, V1, V2, V3])
domain_right.faces.name = "right_domain"
domain_right.edges.name = "kink"
domain_right.edges.Max(X).name = "right"
domain_right.edges.Max(Y).name = "sym"
domain_right.edges.Min(Y).name = "bot"

domain_left1 = domain_left * rec2
domain_left1.edges.Max(X).name = "interface"
domain_left1.edges.Max(Y).name = "sym"
domain_left1.edges.Min(Y).name = "force"
domain_left2 = domain_left - rec2
domain_left2.edges.Min(X).name = "interface"
domain_left2.edges.Max(Y).name = "sym"

domain = Compound(
    [Glue([domain_left1, domain_left2]), domain_right], separate_layers=False
)
mesh = Mesh(OCCGeometry(domain).GenerateMesh(maxh=0.2))

fes = H1(mesh, order=1)
gf = GridFunction(fes)
gf.Set(1, definedon=mesh.Boundaries("left_domain"))
Draw(gf);

fesU = VectorH1(mesh, order=order, dirichlety_bbnd="sym")
fesM = HDivDivSurface(mesh, order=order - 1, dirichlet_bbnd="left|right|bot|force|kink")
fes = fesU * fesM

u, m = fes.TrialFunction()
m = m.Trace()

Regge = HCurlCurl(mesh, order=order - 1, discontinuous=True)

gf_solution = GridFunction(fes, name="solution")

# Surface unit normal, edge-tangential, and co-normal vectors on initial configuration
nv = specialcf.normal(mesh.dim)
tv = specialcf.tangential(mesh.dim)
cnv = Cross(nv, tv)

# Projection to the surface tangent space
Ptau = Id(mesh.dim) - OuterProduct(nv, nv)
# Surface deformation gradient
Ftau = Grad(u).Trace() + Ptau
# Surface Cauchy-Green tensor
Ctau = Ftau.trans * Ftau
# Surface Green-Lagrange strain tensor
Etautau = 0.5 * (Ctau - Ptau)
# surface determinant
J = Norm(Cof(Ftau) * nv)

# Surface unit normal, edge-tangential, and co-normal vectors on deformed configuration
nv_def = Normalize(Cof(Ftau) * nv)
tv_def = Normalize(Ftau * tv)
cnv_def = Cross(nv_def, tv_def)

# Surface Hessian weighted with unit normal vector on deformed configuration
H_nv_def = (u.Operator("hesseboundary").trans * nv_def).Reshape((3, 3))

# Save clamped boundary for updating averaged unit normal vector during load-steps
gf_clamped_bnd = GridFunction(FacetSurface(mesh, order=0))
gf_clamped_bnd.Set(1, definedon=mesh.BBoundaries(""))
gf_sym_bnd = GridFunction(FacetSurface(mesh, order=0))
gf_sym_bnd.Set(1, definedon=mesh.BBoundaries("sym"))

# Unit normal vector on current configuration
# Used to update averaged unit normal vector during load-steps
cf_nv_cur = Normalize(Cof(Ptau + Grad(gf_solution.components[0])) * nv)

# FESpace for averaged normal vectors living only on the edges of the mesh
fesVF = VectorFacetSurface(mesh, order=order - 1)
# GridFunctions to save averaged normal vectors on deformed and initial configurations
averaged_nv = GridFunction(fesVF)
averaged_nv_init = GridFunction(fesVF)

# Initialize by averaging unit normal vectors on initial configuration
# definedon=mesh.Boundaries(".*") is needed as we interpolate on the surface mesh
averaged_nv.Set(nv, dual=True, definedon=mesh.Boundaries(".*"))
averaged_nv_init.Set(nv, dual=True, definedon=mesh.Boundaries(".*"))
# Normalize averaged normal vector on initial configuration
cf_averaged_nv_init_norm = Normalize(averaged_nv_init)
# Project averaged unit normal vector being perpendicular to deformed edge-tangent vector
# to measure correct angle on deformed configuration
cf_proj_averaged_nv = Normalize(averaged_nv - (tv_def * averaged_nv) * tv_def)

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


# Material stress
def MaterialStress(mat, E, nu):
    return E / (1 - nu**2) * ((1 - nu) * mat + nu * Trace(mat) * Ptau)


# Inverse of the material norm
def MaterialNormInv(mat, E, nu):
    return (1 + nu) / E * (InnerProduct(mat, mat) - nu / (nu + 1) * Trace(mat) ** 2)

dS = ds(element_boundary=True)

# Bilinear form for problem
# We define the Lagrangian of the HHJ formulation. Therefore,
# we use Variation() such that Newton's method knows to build the first variation
# for the residual and the second variation for the stiffness matrix. The
# stiffness matrix  will be symmetric.
# We use .Compile() to simplify (linearize) the coefficient expression tree.

par_iso = Parameter(0)
bfA = BilinearForm(fes, symmetric=True)
# Membrane energy = Isometry constraint
# Interpolate the membrane strains Etautau into the Regge space
# to avoid membrane locking
bfA += Variation(
    par_iso * 0.5 / thickness**2 * MaterialNorm(Interpolate(Etautau, Regge), E, nu) * ds
).Compile()
# Bending energy
# Element terms of bending energy
bfA += Variation(
    (-6 * MaterialNormInv(m, E, nu) + InnerProduct(H_nv_def, 1 / J * m)) * ds
).Compile()
# Boundary terms of bending energy
bfA += Variation(
    (
        acos(cnv_def * cf_proj_averaged_nv) - pi / 2
    )  # acos(cnv * cf_averaged_nv_init_norm)=pi/2 for plates
    * m[cnv, cnv]
    * dS
).Compile()

# Parameter par for load-stepping below.
par = Parameter(0.0)
bfA += Variation(
    -(1 - par**2 + 1e-10) * 4e-1 * CF((0, 0, 1)) * u * ds
    + par_iso * 1e6 * (u - par * uD) ** 2 * dx(definedon=mesh.BBoundaries("force"))
)

# Reset the solution to zero
gf_solution.vec[:] = 0.0
gf_Dirichlet = GridFunction(fes)


gf_history = GridFunction(fes, multidim=0)
# Extract components of solution
gfu, gfm = gf_solution.components
gfC = (Grad(gfu) + Ptau).trans * (Grad(gfu) + Ptau)

# Draw displacement, bending moment tensor, and displacement gradient
scene = Draw(gfu, mesh, "displacement", deformation=gfu)
scene2 = Draw(Norm(gfm), mesh, "bending_stress", deformation=gfu)
scene3 = Draw(Grad(gfu)[0], mesh, "gradient_u", deformation=gfu)

# Use num_steps uniform load-steps in [0,1]
num_steps = 20
num_iso_step = 5

with TaskManager():
    for steps in range(num_steps):
        par.Set((steps + 1) / num_steps)
        print("Loadstep =", steps + 1, ", F/Fmax =", (steps + 1) / num_steps * 100, "%")
        # Update averaged normal vector
        # On clamped boundary it remains the unit normal vector of the initial configuration
        # and on the rest of the boundary it is updated by the current unit normal vector
        averaged_nv.Set(
            (1 - gf_clamped_bnd - gf_sym_bnd) * cf_nv_cur
            + gf_sym_bnd * (cf_nv_cur - (cf_nv_cur * cnv) * cnv)
            + gf_clamped_bnd * nv,
            dual=True,
            definedon=mesh.Boundaries(".*"),
        )
        for step_iso in range(num_iso_step):
            if steps == 0:
                par_iso.Set(sqrt((step_iso + 1) / num_iso_step))
            else:
                par_iso.Set(1)
            print("iso_constraint =", par_iso.Get() * 0.5 / thickness**2)

            # Use Newton solver with residual tolerance 1e-8 and maximal 100 iterations
            res, it = solvers.Newton(
                bfA,
                gf_solution,
                inverse="umfpack",
                printing=False,
                dampfactor=0.05 * par_iso.Get(),
                maxerr=1e-6,
                maxit=100,
            )
            if res != 0:
                raise Exception("did not converge")
            # Redraw solutions
            scene.Redraw()
            scene2.Redraw()
            scene3.Redraw()
            if steps > 0:
                break

        gf_history.AddMultiDimComponent(gf_solution.vec)

Loadstep = 1 , F/Fmax = 5.0 %
iso_constraint = 223606.79774997898

iso_constraint = 316227.76601683797

iso_constraint = 387298.3346207417

iso_constraint = 447213.59549995797

iso_constraint = 500000.0

Loadstep = 2 , F/Fmax = 10.0 %
iso_constraint = 500000.0

Loadstep = 3 , F/Fmax = 15.0 %
iso_constraint = 500000.0

Loadstep = 4 , F/Fmax = 20.0 %
iso_constraint = 500000.0

Loadstep = 5 , F/Fmax = 25.0 %
iso_constraint = 500000.0

Loadstep = 6 , F/Fmax = 30.0 %
iso_constraint = 500000.0

Loadstep = 7 , F/Fmax = 35.0 %
iso_constraint = 500000.0

Loadstep = 8 , F/Fmax = 40.0 %
iso_constraint = 500000.0

Loadstep = 9 , F/Fmax = 45.0 %
iso_constraint = 500000.0

Loadstep = 10 , F/Fmax = 50.0 %
iso_constraint = 500000.0

Loadstep = 11 , F/Fmax = 55.00000000000001 %
iso_constraint = 500000.0

Loadstep = 12 , F/Fmax = 60.0 %
iso_constraint = 500000.0

Loadstep = 13 , F/Fmax = 65.0 %
iso_constraint = 500000.0

Loadstep = 14 , F/Fmax = 70.0 %
iso_constraint = 500000.0

Loadstep = 15 , F/Fmax = 75.0 %
iso_constraint = 500000.0

Loadstep = 16 , F/Fmax = 80.0 %
iso_constraint = 500000.0

Loadstep = 17 , F/Fmax = 85.0 %
iso_constraint = 500000.0

Loadstep = 18 , F/Fmax = 90.0 %
iso_constraint = 500000.0

Loadstep = 19 , F/Fmax = 95.0 %
iso_constraint = 500000.0

Loadstep = 20 , F/Fmax = 100.0 %
iso_constraint = 500000.0

Draw(
    gf_history.components[0],
    mesh,
    max=0.9,
    animate=True,
    deformation=gf_history.components[0],
    euler_angles=[-60, 0, -10],
);
Contents