Buckling of a heated von-Kármán plate

This demo is implemented in the single Python file demo_von-karman-mansfield.py.

This demo program solves the von-Kármán equations on a circular plate with a lenticular cross section free on the boundary. The plate is heated, causing it to bifurcate. Bifurcation occurs with a striking shape transition: before the critical threshold the plate assumes a cup-shaped configurations (left); above it tends to a cylindrical shape (right). An analytical solution has been found by Mansfield, see [1].


Pre-critical and post-critical plate configuration.

The standard von-Kármán theory gives rise to a fourth-order PDE which requires the transverse displacement field \(w\) to be sought in the space \(H^2(\Omega)\). We relax this requirement in the same manner as the Kirchoff-Love plate theory can be relaxed to the Reissner-Mindlin theory, resulting in seeking a transverse displacement field \(w\) in \(H^1(\Omega)\) and a rotation field \(\theta\) in \([H^2(\Omega)]^2\). To alleviate the resulting shear- and membrane-locking issues we use Partial Selective Reduced Integration (PSRI), see [2].

To follow this demo you should know how to:

  • Define a MixedElement and a FunctionSpace from it.
  • Write variational forms using the Unified Form Language.
  • Automatically derive Jabobian and residuals using derivative().
  • Apply Dirichlet boundary conditions using DirichletBC and apply().
  • Assemble forms using assemble().

This demo then illustrates how to:

  • Define the Reissner-Mindlin-von-Kármán plate equations using UFL.
  • Use the PSRI approach to simultaneously cure shear- and membrane-locking issues.

We start with importing the required modules, setting matplolib as plotting backend, and generically set the integration order to 4 to avoid the automatic setting of FEniCS which would lead to unreasonably high integration orders for complex forms.

from dolfin import *
from fenics_shells import *
import mshr
import numpy as np
import matplotlib.pyplot as plt

parameters["form_compiler"]["quadrature_degree"] = 4

The mid-plane of the plate is a circular domain with radius \(a = 1\). We generate a mesh of the domain using the package mshr:

a = 1.
ndiv = 8
domain_area = np.pi*a**2

centre = Point(0.,0.)

geom = mshr.Circle(centre, a)
mesh = mshr.generate_mesh(geom, ndiv)

The lenticular thinning of the plate can be modelled directly through the thickness parameter in the plate model:

t = 1E-2
ts = interpolate(Expression('t*(1.0 - (x[0]*x[0] + x[1]*x[1])/(a*a))', t=t, a=a, degree=2), FunctionSpace(mesh, 'CG', 2))

We assume the plate is isotropic with constant material parameters, Young’s modulus \(E\), Poisson’s ratio \(\nu\); shear correction factor \(\kappa\):

E = Constant(1.0)
nu = Constant(0.3)
kappa = Constant(5.0/6.0)

Then, we compute the (scaled) membrane, \(A/t^3\), bending, \(D/t^3\), and shear, \(S/t^3\), plate elastic stiffnesses:

A = (E*ts/t**3/(1. - nu**2))*as_tensor([[1., nu, 0.],[nu, 1., 0.],[0., 0., (1. - nu)/2]])
D = (E*ts**3/t**3/(12.*(1. - nu**2)))*as_tensor([[1., nu, 0.],[nu, 1., 0.],[0., 0., (1. - nu)/2]])
S = E*kappa*ts/t**3/(2*(1. + nu))

We use a \(CG_2\) element for the in-plane and transverse displacements \(u\) and \(w\), and the enriched element \([CG_1 + B_3]\) for the rotations \(\theta\). We collect the variables in the state vector \(q =(u,w,\theta)\):

P1 = FiniteElement("Lagrange", triangle, degree = 1)
P2 = FiniteElement("Lagrange", triangle, degree = 2)
bubble = FiniteElement("B", triangle, degree = 3)
enriched = P1 + bubble

element = MixedElement([VectorElement(P2, dim=2), P2, VectorElement(enriched, dim=2)])

Q = FunctionSpace(mesh, element)

Then, we define Function, TrialFunction and TestFunction objects to express the variational forms and we split them into each individual component function:

q_, q, q_t = Function(Q), TrialFunction(Q), TestFunction(Q)
v_, w_, theta_ = split(q_)

The membrane strain tensor \(e\) for the von-Kármán plate takes into account the nonlinear contribution of the transverse displacement in the approximate form:

\[e(v, w) = \mathrm{sym}\nabla v + \frac{\nabla w \otimes \nabla w}{2}\]

which can be expressed in UFL as:

e = sym(grad(v_)) + 0.5*outer(grad(w_), grad(w_))

The membrane energy density \(\psi_m\) is a quadratic function of the membrane strain tensor \(e\). For convenience, we use our function strain_to_voigt() to express \(e\) in Voigt notation \(e_V = \{e_1, e_2, 2 e_{12}\}\):

ev = strain_to_voigt(e)
psi_m = 0.5*dot(A*ev, ev)

The bending strain tensor \(k\) and shear strain vector \(\gamma\) are identical to the standard Reissner-Mindlin model. The shear energy density \(\psi_s\) is a quadratic function of the shear strain vector:

gamma = grad(w_) - theta_
psi_s = 0.5*dot(S*gamma, gamma)

The bending energy density \(\psi_b\) is a quadratic function of the bending strain tensor. Here, the temperature profile on the plate is not modelled directly. Instead, it gives rise to an inelastic (initial) bending strain tensor \(k_T\) which can be incoporated directly in the Lagrangian:

k_T = as_tensor(Expression((("c/imp","0.0"),("0.0","c*imp")), c=1., imp=.999,  degree=0))
k = sym(grad(theta_)) - k_T
kv = strain_to_voigt(k)
psi_b = 0.5*dot(D*kv, kv)


The standard von-Kármán model can be recovered by substituting in the Kirchoff constraint \(\theta = \nabla w\).

Shear- and membrane-locking are treated using the partial reduced selective integration proposed by Arnold and Brezzi, see [2]. In this approach shear and membrane energy are splitted as a sum of two contributions weighted by a factor \(\alpha\). One of the two contributions is integrated with a reduced integration. We use \(2\times 2\)-points Gauss integration for a portion \(1-\alpha\) of the energy, whilst the rest is integrated with a \(4\times 4\) scheme. We adopt an optimized weighting factor \(\alpha=(t/h)^2\), where \(h\) is the mesh size.:

dx_h = dx(metadata={'quadrature_degree': 2})
h = CellDiameter(mesh)
alpha = project(t**2/h**2, FunctionSpace(mesh,'DG',0))

Pi_PSRI = psi_b*dx + alpha*psi_m*dx + alpha*psi_s*dx + (1.0 - alpha)*psi_s*dx_h + (1.0 - alpha)*psi_m*dx_h

Then, we compute the total elastic energy and its first and second derivatives:

Pi = Pi_PSRI
dPi = derivative(Pi, q_, q_t)
J = derivative(dPi, q_, q)

This problem is a pure Neumann problem. This leads to a nullspace in the solution. To remove this nullspace, we fix all the variables in the central point of the plate and the displacements in the \(x_0\) and \(x_1\) direction at \((0, a)\) and \((a, 0)\), respectively:

zero_v1 = project(Constant((0.)), Q.sub(0).sub(0).collapse())
zero_v2 = project(Constant((0.)), Q.sub(0).sub(1).collapse())
zero = project(Constant((0.,0.,0.,0.,0.)), Q)

bc = DirichletBC(Q, zero, "near(x[0], 0.) and near(x[1], 0.)", method="pointwise")
bc_v1 = DirichletBC(Q.sub(0).sub(0), zero_v1, "near(x[0], 0.) and near(x[1], 1.)", method="pointwise")
bc_v2 = DirichletBC(Q.sub(0).sub(1), zero_v2, "near(x[0], 1.) and near(x[1], 0.)", method="pointwise")
bcs = [bc, bc_v1, bc_v2]

Then, we define the nonlinear variational problem and the solver settings:

init = Function(Q)
problem = NonlinearVariationalProblem(dPi, q_, bcs, J = J)
solver = NonlinearVariationalSolver(problem)
solver.parameters["newton_solver"]["absolute_tolerance"] = 1E-8

Finally, we choose the continuation steps (the critical loading c_cr is taken from the Mansfield analytical solution [1]):

c_cr = 0.0516
cs = np.linspace(0.0, 1.5*c_cr, 30)

and we solve as usual:

defplots_dir = "output/3dplots-psri/"
file = File(defplots_dir + "sol.pvd")

ls_kx = []
ls_ky = []
ls_kxy = []
ls_kT = []

for count, i in enumerate(cs):
    k_T.c = i
    v_h, w_h, theta_h = q_.split(deepcopy=True)

To visualise the solution we assemble the bending strain tensor:

K_h = project(sym(grad(theta_h)), TensorFunctionSpace(mesh, 'DG', 0))

we compute the average bending strain:

Kxx = assemble(K_h[0,0]*dx)/domain_area
Kyy = assemble(K_h[1,1]*dx)/domain_area
Kxy = assemble(K_h[0,1]*dx)/domain_area


and output the results at each continuation step:

v1_h, v2_h = v_h.split()
u_h = as_vector([v1_h, v2_h, w_h])
u_h_pro = project(u_h, VectorFunctionSpace(mesh, 'CG', 1, dim=3))
file << u_h_pro

Finally, we plot the average curvatures as a function of the inelastic curvature:

fig = plt.figure(figsize=(5.0, 5.0/1.648))
plt.plot(ls_kT, ls_kx, "o", color='orange', label=r"$k_{1h}$")
plt.plot(ls_kT, ls_ky, "x", color='red', label=r"$k_{2h}$")
plt.xlabel(r"inelastic curvature $\eta$")
plt.ylabel(r"curvature $k_{1,2}$")

Comparison with the analytical solution.

Unit testing

import pytest
def test_close():


[1] E. H. Mansfield, “Bending, Buckling and Curling of a Heated Thin Plate. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. Vol. 268. No. 1334. The Royal Society, 1962.

[2] D. Arnold and F.Brezzi, Mathematics of Computation, 66(217): 1-14, 1997. https://www.ima.umn.edu/~arnold//papers/shellelt.pdf