Clamped Kirchhoff-Love plate¶
This demo is implemented in a single Python file
demo_kirchhoff-love-clamped.py.
This demo program solves the out-of-plane Kirchhoff-Love equations on the unit square with uniform transverse loading and fully clamped boundary conditions.
We use the Continuous/Discontinuous Galerkin (CDG) formulation to allow the use of \(H^1\)-conforming elements for this fourth-order, or \(H^2\)-type problem. This demo is very similar to the Biharmonic equation demo in the main DOLFIN repository, and as such we recommend reading that page first. The main differences are:
- we express the stabilisation term in terms of a Lagrangian functional rather than as a bilinear form,
- the Lagrangian function in turn is expressed in terms of the bending moment and rotations rather than the primal field variable,
- and we show how to place Dirichlet boundary conditions on the first-derivative of the solution (the rotations) using a weak approach.
We begin as usual by importing the required modules:
from dolfin import *
from fenics_shells import *
and creating a mesh:
mesh = UnitSquareMesh(64, 64)
We use a second-order scalar-valued element for the transverse displacement field \(w \in CG_2\):
element_W = FiniteElement("Lagrange", triangle, 2)
W = FunctionSpace(mesh, element_W)
and then define Function, TrialFunction and
TestFunction objects to express the variational form:
w_ = Function(W)
w = TrialFunction(W)
w_t = TestFunction(W)
We take constant material properties throughout the domain:
E = Constant(10920.0)
nu = Constant(0.3)
t = Constant(1.0)
The Kirchhoff-Love model, unlike the Reissner-Mindlin model, is a rotation-free model: the rotations \(\theta\) do not appear explicitly as degrees of freedom. Instead, the rotations of the Kirchhoff-Love model are calculated from the transverse displacement field as:
which can be expressed in UFL as:
theta = grad(w_)
The bending tensor can then be calculated from the derived rotation field in exactly the same way as for the Reissner-Mindlin model:
or in UFL:
k = variable(sym(grad(theta)))
The function variable() annotation is important and will allow us to
take differentiate with respect to the bending tensor k to derive the
bending moment tensor, as we will see below.
Again, identically to the Reissner-Mindlin model we can calculate the bending energy density as:
D = (E*t**3)/(24.0*(1.0 - nu**2))
psi_M = D*((1.0 - nu)*tr(k*k) + nu*(tr(k))**2)
For the definition of the CDG stabilisation terms and the (weak) enforcement of the Dirichlet boundary conditions on the rotation field, we need to explicitly derive the moment tensor \(M\). Following standard arguments in elasticity, a stress measure (here, the moment tensor) can be derived from the bending energy density by taking its derivative with respect to the strain measure (here, the bending tensor):
or in UFL:
M = diff(psi_M, k)
We now move onto the CDG stabilisation terms.
Consider a triangulation \(\mathcal{T}\) of the domain \(\omega\) with the set of interior edges is denoted \(\mathcal{E}_h^{\mathrm{int}}\). Normals to the edges of each facet are denoted \(n\). Functions evaluated on opposite sides of a facet are indicated by the subscripts \(+\) and \(-\).
The Lagrangian formulation of the CDG stabilisation term is then:
Furthermore, \(\left< u \right> = \frac{1}{2} (u_{+} + u_{-})\) operator, \([\!\![ u ]\!\!] = u_{+} \cdot n_{+} + u_{-} \cdot n_{-}\), \(\alpha \ge 0\) is a penalty parameter and \(h_E\) is a measure of the cell size. We choose the penalty parameter to be on the order of the norm of the bending stiffness matrix \(\dfrac{Et^3}{12}\).
This can be written in UFL as:
alpha = E*t**3
h = CellDiameter(mesh)
h_avg = (h('+') + h('-'))/2.0
n = FacetNormal(mesh)
M_n = inner(M, outer(n, n))
L_CDG = -inner(jump(theta, n), avg(M_n))*dS + \
(1.0/2.0)*(alpha('+')/h_avg)*inner(jump(theta, n), jump(theta, n))*dS
We now define our Dirichlet boundary conditions on the transverse displacement field:
class AllBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary
# Boundary conditions on displacement
all_boundary = AllBoundary()
bcs_w = [DirichletBC(W, Constant(0.0), all_boundary)]
Because the rotation field \(\theta\) does not enter our weak formulation directly, we must weakly enforce the Dirichlet boundary condition on the derivatives of the transverse displacement \(\nabla w\).
We begin by marking the exterior facets of the mesh where we want to apply boundary conditions on the rotation:
facet_function = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)
facet_function.set_all(0)
all_boundary.mark(facet_function, 1)
and then define an exterior facet Measure object from that
subdomain data:
ds = Measure("ds")(subdomain_data=facet_function)
In this example, we would like \(\theta_d = 0\) everywhere on the boundary:
theta_d = Constant((0.0, 0.0))
The definition of the exterior facets and Dirichlet rotation field were trivial in this demo, but you could extend this code straightforwardly to non-homogeneous Dirichlet conditions.
The weak boundary condition enforcement term can be written:
where \(\theta_e = \theta - \theta_d\) is the effective rotation field, and \(\mathcal{E}_h^{\mathrm{D}}\) is the set of all exterior facets of the triangulation \(\mathcal{T}\) where we would like to apply Dirichlet boundary conditions, or in UFL:
theta_effective = theta - theta_d
L_BC = -inner(inner(theta_effective, n), M_n)*ds(1) + \
(1.0/2.0)*(alpha/h)*inner(inner(theta_effective, n), inner(theta_effective, n))*ds(1)
The remainder of the demo is as usual:
f = Constant(1.0)
W_ext = f*w_*dx
L = psi_M*dx - W_ext + L_CDG + L_BC
F = derivative(L, w_, w_t)
J = derivative(F, w_, w)
A, b = assemble_system(J, -F, bcs=bcs_w)
solver = PETScLUSolver("mumps")
solver.solve(A, w_.vector(), b)
XDMFFile("output/w.xdmf").write(w_)
Unit testing¶
def test_close():
import numpy as np
assert(np.isclose(w_((0.5, 0.5)), 1.265E-6, atol=1E-3, rtol=1E-3))