# Thin-plate FEM simulation

I want to simulate a vibration of a thin plate (the Kirchhoff-Love model) on triangular meshes. Can you advise me on an introductory-level review of different elements for thin-plate FEM simulation?

I was thinking about Morley elements, but they seem non-conforming (is it really bad? never worked with non-conforming FEM before). Then we have many other elements, and the mixed formulation, and so on - I am astonished by how many approaches there are. Ideally, I would prefer something not very advanced in understanding and doable with open-source FEM libraries.

Any advice is greatly appreciated!

I usually suggest trying Morley first because it's fairly easy to implement and in many cases you just do everything just as when using a conforming finite element. Nonconformity will be a factor in the error analysis but otherwise it is very similar.

If you want to try it out here is a snippet which works, e.g., in Google Colab after running !pip install scikit-fem==5.2.0. It is solving $$\Delta^2 u = \lambda u$$ with $$u = \frac{\partial u}{\partial n}=0$$ on the boundary:

from skfem import *
from skfem.helpers import ddot, dd
from skfem.models import mass
from skfem.visuals.matplotlib import plot

@BilinearForm
def biharmonic(u, v, w):
return ddot(dd(u), dd(v))

m = MeshTri.init_lshaped().refined(4)
basis = Basis(m, ElementTriMorley())

# matrices
A = biharmonic.assemble(basis)
M = mass.assemble(basis)

# solve eigenvalue problem
D = basis.get_dofs()
x = solve(*condense(A, M, D=D))
nth = 4
x0 = x[1][:, nth]

ax = plot(basis, x0, nrefs=2, shading='gouraud', figsize=(4, 4))


Above snippet draws the eigenfunction but I made you an animation of the fourth mode using some extra code and matplotlib.animation:

You really don't want to use conforming elements. Use a penalty approach instead.

The reasoning is described in detail here, along with an implementation based on an open source library (of which, disclaimer!, I am one of the principal authors): https://dealii.org/developer/doxygen/deal.II/step_47.html

An implementation of an alternative approach can be found here, implemented with the same library: https://dealii.org/developer/doxygen/deal.II/step_82.html

This is a brief summary of what was answered.

Roughly speaking, there are 3 options

1. Conforming elements like Argyris triangle, Hsien-Clough-Ticher triangle, and others. People seem to be skeptical about them, perhaps, because the number of DOFs per element is high (21 for AT, 12 for HCT - see https://defelement.com/

2. Non-conforming element with the Morley triangle being the canonical example. The MT has only 6 DOFs. Elements of that kind are difficult to implement, especially in 3D. Another disadvantage is that the MT has only the 1st order of convergence (as a consequence of the non-conformity). However, as soon as such elements are available, they can be quite efficient - see @RobertKirby 's reply (a good paper included!). @knl and @lightxbulb provided interesting examples on the MT. To my knowledge, the MT is implemented in www.firedrakeproject.org, FreeFEM, scikit-fem, FEniCS.

3. Interior-penalty and DG formulations. They are relatively straightforward to implement - see @WolfgangBangerth 's reply with FANTASTIC examples with deal.II package. This strategy can be implemented pretty much in any FE package, provided you know how to efficiently solve the resulting saddle-point system of linear equations.

Firedrake (www.firedrakeproject.org) supports $$C^1$$ triangles (nonconforming Morley plus conforming Bell and Argyris) with high-level syntax for defining bilinear forms. We've found here that doing such elements gives smaller, sparser, better conditioned linear systems than interior penalty-type methods. The major disadvantage to Morley being nonconforming is that it's only first-order accurate, which is also the case for the $$P^2$$ interior penalty method. The punch line is that if you have access to the $$C^1$$ elements, there's not a major reason to prefer interior penalty, but IP methods will do the job if you don't have ready access to $$C^1$$ methods.

The story would be a bit different if you need to to do $$C^1$$ problems in 3D, as we don't have good $$C^1$$ elements implemented (they are far more complicated than in 2D). I imagine that IP methods would work fine.

• As it’s currently written, your answer is unclear. Please edit to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers in the help center.
– Community Bot
Feb 7 at 17:16
• Argyris is not so easy if you want to do nonhomogeneous boundary conditions or the boundary conditions are changing, or even interface conditions (material changes). Suppose you are doing plate on L-shaped domain. In the L-shaped corner, which of the Argyris DOFs should be set to zero? What if the corner is less than 90 degrees? What if the boundary condition changes in the corner from clamped to simply supported?
– knl
Feb 8 at 13:21