Clamped Circular Plate under Uniform Load (Kirchhoff-Love)

Problem Description

This benchmark tests the Kirchhoff-Love plate model against an analytical solution. The problem is a clamped circular plate of radius \(R\) subject to a uniform out-of-plane load \(F\).

In the Kirchhoff-Love model the rotation field is identified with the gradient of the vertical deformation, \(\boldsymbol{\theta} = \nabla w\), so the fourth-order equilibrium equation is

\begin{align}\nabla\cdot(\nabla\cdot(\mathbf{M}(w))) = F \quad \text{in } \Omega,\end{align}

where the moment tensor for an isotropic material is \(\mathbf{M}(w) = -D(\mathbf{x})\left\{(1-\nu)\nabla\nabla w + \nu\,\Delta w\,\mathbf{I}\right\}\) and \(D(\mathbf{x}) = Et^3/(12(1-\nu^2))\) is the (potentially spatially varying) bending modulus. For uniform \(D\) this simplifies to \(D\,\nabla^4 w = F\). The plate is clamped at the boundary:

\begin{align}w = 0, \quad \nabla w \cdot \boldsymbol{n} = 0 \quad \text{on } \partial\Omega.\end{align}

The fourth-order equation is reformulated as a coupled system of second-order equations using a Helmholtz decomposition of the shear resultant (see the model documentation).

Analytical Solution

For a clamped circular plate of radius \(R\) under uniform load \(F\) the analytical solution is

\begin{align}w(r) = -\frac{F}{64D}\left(R^2 - r^2\right)^2,\end{align}

where \(r = \sqrt{x^2 + y^2}\).

Parameters

Parameter	Symbol	Value	Unit
Young's modulus	\(E\)	\(10^6\)	Pa
Poisson ratio	\(\nu\)	0.3	-
Plate thickness	\(t\)	0.05	m
Bending modulus	\(D = Et^3/(12(1-\nu^2))\)	\(\approx 11.45\)	Pa·m³
Load	\(F\)	1.0	N/m²
Radius	\(R\)	1.0	m

Setup

The domain is a unit disk discretized with an unstructured triangular mesh generated by Gmsh from disk.geo. The mesh size \(h\) controls the resolution. The mixed formulation uses two coupled sub-problems discretized with control-volume finite element (CVFE) methods: the Box scheme (PQ1 CVFE (Box) scheme) for the deformation sub-problem \((\varphi, w, \psi)\) and the PQ1-Bubble scheme (PQ1-B CVFE scheme) for the rotation sub-problem \(\boldsymbol{\theta}\), coupled via the Multidomain framework framework.

Results

To run the convergence test and produce the convergence plot, run:

python3 plot_convergence.py <executable>

The script produces two figures:

convergence.png: relative \(L^2\) error of \(w\) versus mesh size \(h\) on a log-log scale; the expected convergence rate is \(\mathcal{O}(h^2)\).
deformation_3d.png: 3D surface of the deformed plate (requires PyVista).

Convergence plot

Deformation

🛠 Edit above doc on GitLab