Thin-plate FEM simulation

Question

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!

Answer 1

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:

Answer 2

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

Answer 3

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.

Answer 4

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.