4. Reissner Mindlin plates

Find: vertical deflection \(w\) and linearized rotations \(\beta\) minimizing total energy

\[ \frac{1}{2} \| \varepsilon(\beta) \|_D^2 + \frac{1}{2t^2} \| \nabla w - \beta \|_0^2 - \int f w \]

First term is bending energy, the second is shear energy.

In general we obtain shear locking for \(t \rightarrow 0\).

Can be avoided if discrete functions

\[ \nabla w_h \quad \text{and} \quad \beta_h \]

match. Choose \(\beta_h\) in the Nedelec - space !

To discretize the bending energy, we use a mixed method. Introduce moments \(\sigma\) such that \(A \sigma = \varepsilon(\beta)\) to obtain mixed formulation:

Find \(\sigma, w, \beta\):

\[\begin{split} \begin{array}{ccccl} (\sigma, \tau)_A & + & \left< \operatorname{div} \sigma, \beta \right> & = & 0 \\ \left< \operatorname{div} \tau, \delta \right> & - & \tfrac{1}{t^2} (\nabla w - \beta, \nabla v - \delta) & = & \int f v \end{array} \end{split}\]

The pairing \(\left< \operatorname{div} \sigma, \beta \right>\) is well defined, continuous and LBB-stable for \(\sigma \in H(\operatorname{div} \operatorname{div})\) and \(\beta \in H(\operatorname{curl})\).

For \(t \rightarrow 0\), the finite element solution converges to the HHJ solution of the Kirchhoff plate.

A. Pechstein+JS: The TDNNS method for Reissner–Mindlin plates, 2017

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

order = 2
Sigma = HDivDiv(mesh, order=order, plus=True)  
W = H1(mesh, order=order+1, dirichlet=".*")  
V = HCurl(mesh, order=order, dirichlet=".*")
X = Sigma * W * V

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

ndof-Sigma: 4085 , ndof-W: 1096 , ndof-X: 6966

\[A(\sigma,u; \tau, v) = a(\sigma,\tau) + b(\sigma,\delta) + b(\tau, \beta) - \frac{1}{t^2} (\nabla w - \beta, \nabla v - \delta) \]

with

\[\begin{eqnarray*} a(\sigma,\tau) & = & \int_\Omega \sigma \tau \, dx \\ b(\sigma, \beta) & = & \sum_T \int_T \text{div} \sigma \, \beta \, dx - \int_{\partial T} \sigma_{nn} \beta_n \, ds \end{eqnarray*}\]

t = 0.01

sigma, w, beta = X.TrialFunction()
tau, v, delta = X.TestFunction()

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

a = BilinearForm(InnerProduct(sigma,tau)*dx + DivDiv(sigma,delta) + DivDiv(tau,beta) \
                 - 1/t**2 * (grad(w)-beta)*(grad(v)-delta)*dx).Assemble()
f = LinearForm(200*v*dx).Assemble()

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

gfsigma, gfw, gfbeta = gfu.components
Draw (gfw, deformation=True)
Draw (gfbeta, vectors= { "grid_size" : 40 } )
Draw (gfu.components[0][0,0],mesh);