Fourth Order Equation

44. Fourth Order Equation#

We consider the Kirchhoff plate equation: Find \(w \in H^2\), such that

\[ \int \nabla^2 w : \nabla^2 v = \int f v \]

A conforming method requires \(C^1\) continuous finite elements. But there is no good option available, and thus there is no \(H^2\) conforming finite element space in NGSolve.

44.1. Hybridized \(C^0\)-continuous interior penalty method:#

A simple way out is to use continuous elements, and treat the missing \(C^1\)-continuity by a Discontinuous Galerkin method. A DG formulation is

\[ \sum_T \nabla^2 w : \nabla^2 v - \int_{E} \{\nabla^2 w\}_{nn} \, [\partial_n v] - \int_{E} \{\nabla^2 v\}_{nn} \, [\partial_n w] + \alpha \int_E [\partial_n w] [\partial_n v] \]

[Baker 77, Brenner Gudi Sung, 2010]

We consider its hybrid DG version, where the normal derivative is a new, facet-based variable:

\[ \sum_T \nabla^2 w : \nabla^2 v - \int_{\partial T} (\nabla^2 w)_{nn} \, (\partial_n v - \widehat{v_n}) - \int_{\partial T} (\nabla^2 v)_{nn} \, (\partial_n w - \widehat{w_n}) + \alpha \int_E (\partial_n v - \widehat{v_n}) (\partial_n w - \widehat{w_n}) \]

The facet variable is the normal derivative \(n_E \cdot \nabla w\), what is oriented along the arbitrarily chosen edge normal-vector. We cannot use the FacetSpace since it does not have the orientation, but we can use the normal traces of an HDiv space. We don’t need inner basis functions, so we set order inner to 0:

from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh (unit_square.GenerateMesh(maxh=0.1))

order = 5

V1 = H1(mesh, order=order, dirichlet="left|bottom")
V2 = NormalFacetFESpace(mesh, order=order-1) # , dirichlet="left|bottom")

V = FESpace ([V1,V2])

w,what = V.TrialFunction()
v,vhat = V.TestFunction()

Some proxy-functions and gridfunctions provide additional differential operators. We can get them via the Operator function. w.Operator(“hesse”) provides the Hessian, a matrix-valued function. Note that we can use the InnerProduct(.,.) for \(\nabla^2 w : \nabla^2 v\):

n = specialcf.normal(2)
h = specialcf.mesh_size

def jumpdn(v,vhat): 
    return n*(grad(v)-vhat)
def hesse(v):
    return v.Operator("hesse")
def hessenn(v):
    return InnerProduct(n, hesse(v)*n)

a = BilinearForm(V)
a += InnerProduct (hesse(w), hesse(v)) * dx
a += -hessenn(w) * jumpdn(v,vhat) * dx(element_boundary=True)
a += -hessenn(v) * jumpdn(w,what) * dx(element_boundary=True)
a +=  3*order*order/h * jumpdn(w,what) * jumpdn(v,vhat) * dx(element_boundary=True)
a.Assemble()

f = LinearForm(V)
f += 3 * v * dx
f.Assemble();

u = GridFunction(V)
u.vec.data = a.mat.Inverse(V.FreeDofs()) * f.vec

Draw (u.components[0], mesh, "disp_DG")
Draw (grad (u.components[0])[0], mesh, "grad")
Draw (hesse (u.components[0])[0,0], mesh, "hesse");

Exercise: Change the Dirichlet-boundaries for the normal derivative to switch between simply supported and fully clamped plates

Contents