Kirchhoff plate equation

58. Kirchhoff plate equation#

\(\DeclareMathOperator{\operatorname{div}}{div}\) Find vertical deflection \(w\) such that

\[ \operatorname{div} \operatorname{div} \nabla^2 w = f \qquad \text{ + boundary conditions} \]

Reduce to second order system with bending moments \(\sigma\):

\[\begin{eqnarray*} \sigma - \nabla^2 w & = & 0 \\ \operatorname{div} \operatorname{div} \sigma & = & f \end{eqnarray*}\]

Variational formulation with \(\sigma \in H(\operatorname{div})^{SYM}\) and \(w \in H^1\):

\[\begin{split} \begin{array}{ccccll} \int \sigma \tau & + & \int \operatorname{div} \tau \, \nabla w & = & 0 & \forall \, \tau \\[0.5em] \int \operatorname{div} \sigma \nabla v & & & = & \int f v & \forall \, v \end{array} \end{split}\]

\(H(\operatorname{div})^{SYM}\) requires finite elements which are symmetric, and the normal vector \(\sigma n\) is continuous. This is non-trivial, but possible [Arnold-Winther, …]

58.1. Hellan Herrmann Johnson Method#

Hellan 67, Herrmann 67, Johnson 73, Brezzi-Raviart 77, Arnold+Brezzi 85, Comodi 89, Krendl+Rafetseder+Zulehner 16, Chen+Hu+Huang 16, Braess+Pechstein+JS 20

Is an arbitrary order (\(k \geq 0\)) mixed discretization method on possibly curved, mapped triangular elements.

Allows hybridization to obtain a positive definite system matrix, lowest order case corresponds to the Morley element.

Function space $\( H(\operatorname{div} \operatorname{div}) = \{ \sigma \in L_2^{2\times 2, SYM} : \operatorname{div} \operatorname{div} \sigma \in H^{-1} \} \)$

HHJ finite element space: symmetric matrices, piecewise polynomial, continuous nn-component

Find \(\sigma_h \in \Sigma_h \subset H(\operatorname{div} \operatorname{div})\) and \(w_h \in W_h \subset H^1\) such that

\[\begin{split} \begin{array}{ccccll} \int \sigma_h \tau & + & \sum_T \int_T \operatorname{div} \tau \, \nabla w_h + \int_{\partial T} \tau_{nt} \nabla_t w_h & = & 0 & \forall \, \tau \\[0.5em] \sum_T \int_T \operatorname{div} \sigma_h \nabla v + \int_{\partial T} \sigma_{nt} \nabla_t v & & & = & \int f v & \forall \, v \end{array} \end{split}\]

This mixed method satisfies the magic discrete kernel inclusion

\[ V_{h,0} \subset V_0 \]

leading to the best-approximation property of the bending moments \(\sigma\):

\[ \| \sigma - \sigma_h \|_{L_2} \leq \inf_{\tau_h \in \Sigma_h} \| \sigma - \tau_h \|_{L_2} + \| f - I_h f \| \]

58.1.1. Distributional derivatives:#

for continuous, piece-wise smooth \(w\), and smooth \(\tau\):

\[ \left< \nabla^2 w , \tau \right> = \sum_T \int_T \nabla^2 w : \tau + \sum_E [\partial_n w] \tau_{nn} \]

for \(nn\)-continuous, piece-wise smooth \(\sigma\), and smooth \(v\):

\[\begin{eqnarray*} \left< \operatorname{div} \operatorname{div} \sigma, v \right> & = & \sum_T \int_T \operatorname{div} \operatorname{div} \sigma v + \sum_E ([\operatorname{div} \sigma]_n - \operatorname{div}_t [\sigma n] ) \, v \\ & + & \sum_V \sum_{T, V\in T} (\sigma_{n_1,t_1} - \sigma_{n_2,t_2}) \, v \end{eqnarray*}\]
from ngsolve import *
from ngsolve.webgui import Draw

mesh = Mesh (unit_square.GenerateMesh(maxh=0.1))

Define product finite element space for \(H(\text{div div}) \times H^1\)

order = 2
Sigma = HDivDiv(mesh, order=order)  # dirichlet="left|right|top")
W = H1(mesh, order=order+1, dirichlet=".*")   # "bottom|left|right")
X = Sigma * W

print ("ndof-Sigma:", Sigma.ndof, ", ndof-W:", W.ndof, ", ndof-X:", X.ndof)

ndof-Sigma: 3192 , ndof-W: 1105 , ndof-X: 4297
\[A(\sigma,w; \tau, v) = a(\sigma,\tau) + b(\sigma,v) + b(\tau, w)\]

with

\[\begin{eqnarray*} a(\sigma,\tau) & = & \int_\Omega \sigma \tau \, dx \\ b(\sigma, v) & = & \sum_T \int_T \text{div} \sigma \, \nabla v \, dx - \int_{\partial T} \sigma_n \nabla_t v \, ds \end{eqnarray*}\]
sigma, w = X.TrialFunction()
tau, v = X.TestFunction()

n = specialcf.normal(2)
def tang(u): return u-(u*n)*n
    
def DivDiv(sigma,v): 
    return div(sigma)*grad(v)*dx - (sigma*n)*tang(grad(v))*dx(element_boundary=True)

a = BilinearForm(InnerProduct(sigma,tau)*dx \
                 + DivDiv(sigma,v) + DivDiv(tau,w) \
                 - 1e-10*w*v*dx).Assemble()
f = LinearForm(200*v*dx).Assemble()

gfu = GridFunction(X)
gfu.vec.data = a.mat.Inverse(X.FreeDofs(), inverse="sparsecholesky") * f.vec

the vertical deflection:

gfsigma, gfw = gfu.components
Draw (gfw, mesh, deformation=True);

the bending moments \(\sigma_{xx}\):

Draw (gfsigma[0,0], mesh);

58.2. Hybridization:#

Lagrange-parameter $\( \widehat{w_n} \approx \partial_n w \)\( enforces continuity of \)\sigma_{nn}$.

order = 2
Sigma = Discontinuous(HDivDiv(mesh, order=order)) 
W = H1(mesh, order=order+1, dirichlet=".*")   # "bottom|left|right")
What = NormalFacetFESpace (mesh, order=order, dirichlet=".*")
X = Sigma * W * What
print (X.ndof)

6385

sigma, w, what = X.TrialFunction()
tau, v, vhat = X.TestFunction()

n = specialcf.normal(2)
def tang(u): return u-(u*n)*n
def normal(u): return (u*n)*n

def DivDiv(sigma,v,vhat): 
    return div(sigma)*grad(v)*dx - (sigma*n)*(tang(grad(v))-normal(vhat))*dx(element_boundary=True)

a = BilinearForm(InnerProduct(sigma,tau)*dx + DivDiv(sigma,v,vhat) + DivDiv(tau,w,what), \
                 condense=True).Assemble()
f = LinearForm(200*v*dx).Assemble()

# inverse from static conensation
invSchur = a.mat.Inverse(X.FreeDofs(True), inverse="sparsecholesky")
ext = IdentityMatrix()+a.harmonic_extension
inv =  ext @ invSchur @ ext.T + a.inner_solve

gfu = GridFunction(X)
gfu.vec.data = inv * f.vec

gfsigma, gfw, gfwhat = gfu.components
Draw (gfu.components[1], mesh, deformation=True)
Draw (gfu.components[0][0,0], mesh);
Contents