Shear-locking in thick plate models with quadrilateral elements #
Objectives
This tour implements a thick plate Reissner-Mindlin model and discusses the issue of shear-locking in the thin plate limit. We show how to use selective reduced integration to alleviate this issue.
Download sources
See also
For more details on the Reissner-Mindlin plate model and its implementation, we refer to the Reissner-Mindlin plates introductory tour.
See also
For more advanced Reissner-Mindlin plate elements, we refer to FEniCSx-Shells.
Introduction#
This program solves the Reissner-Mindlin plate equations on the unit square with uniform transverse loading and fully clamped boundary conditions. We will test the use of quadrilateral cells and selective reduced integration (SRI) to remove shear-locking issues in the thin plate limit. Both linear and quadratic interpolation are considered for the transverse deflection
Implementation#
We first create a unit square mesh divided in
import numpy as np
from mpi4py import MPI
import ufl
import basix
from dolfinx import fem, mesh
import dolfinx.fem.petsc
N = 50
domain = mesh.create_unit_square(MPI.COMM_WORLD, N, N, mesh.CellType.quadrilateral)
Material parameters for isotropic linear elastic behavior are first defined:
E = fem.Constant(domain, 210e3)
nu = fem.Constant(domain, 0.3)
Plate bending stiffness
thick = fem.Constant(domain, 1e-3)
D = E * thick**3 / (1 - nu**2) / 12.0
F = E / 2 / (1 + nu) * thick * 5.0 / 6.0
The uniform loading
f = -D / 1.265319087e-3 # with this we have w_Love-Kirchhoff = 1.0
Continuous interpolation using of degree equispaced
or GLL
variants which mostly differ in the dof location up to
name="S"
, it has 8 dofs, spanning polynomials up to order 2, except for
see details here
deg = 2
el_type = "S" # or "Q"
We = basix.ufl.element(el_type, domain.basix_cell(), deg)
Te = basix.ufl.element(el_type, domain.basix_cell(), deg, shape=(2,))
V = fem.functionspace(domain, basix.ufl.mixed_element([We, Te]))
Clamped boundary conditions on the lateral boundary are defined as::
# Boundary of the plate
def border(x):
return np.logical_or(
np.logical_or(np.isclose(x[0], 0), np.isclose(x[0], 1)),
np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1)),
)
facet_dim = 1
clamped_facets = mesh.locate_entities_boundary(domain, facet_dim, border)
clamped_dofs = fem.locate_dofs_topological(V, facet_dim, clamped_facets)
u0 = fem.Function(V)
bcs = [fem.dirichletbc(u0, clamped_dofs)]
Some useful functions for implementing generalized constitutive relations are now defined:
def strain2voigt(eps):
return ufl.as_vector([eps[0, 0], eps[1, 1], 2 * eps[0, 1]])
def voigt2stress(S):
return ufl.as_tensor([[S[0], S[2]], [S[2], S[1]]])
def curv(u):
(w, theta) = ufl.split(u)
return ufl.sym(ufl.grad(theta))
def shear_strain(u):
(w, theta) = ufl.split(u)
return ufl.grad(w) - theta
def bending_moment(u):
DD = ufl.as_tensor([[D, nu * D, 0], [nu * D, D, 0], [0, 0, D * (1 - nu) / 2.0]])
return voigt2stress(ufl.dot(DD, strain2voigt(curv(u))))
def shear_force(u):
return F * shear_strain(u)
The contribution of shear forces to the total energy is under-integrated using
a custom quadrature rule of degree
See also
See the Quadrature schemes tour for more details on the choice of quadrature points.
u = fem.Function(V)
u_ = ufl.TestFunction(V)
du = ufl.TrialFunction(V)
dx = ufl.Measure("dx")
dx_shear = ufl.Measure("dx", metadata={"quadrature_degree": 2 * deg - 2})
L = f * u_[0] * dx
a = (
ufl.inner(bending_moment(u_), curv(du)) * dx
+ ufl.dot(shear_force(u_), shear_strain(du)) * dx_shear
)
We then solve for the solution and print the deflection normalized with respect to the Love-Kirchhoff thin plate analytical solution:
problem = fem.petsc.LinearProblem(
a, L, u=u, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
)
problem.solve()
w = u.sub(0).collapse()
w.name = "Deflection"
print(f"Reissner-Mindlin FE deflection: {max(abs(w.vector.array)):.5f}")
Reissner-Mindlin FE deflection: 1.00002
Results#
We provide here some results for
Type |
||
---|---|---|
|
0.00046 |
0.01116 |
|
0.99261 |
0.99972 |
|
0.72711 |
0.99864 |
|
0.87658 |
1.00002 |
|
0.96450 |
0.99865 |
|
1.00021 |
1.00002 |
The results show that the low-order elements S
exhibits a notable locking behavior which is not necessarily fixed using SRI. On the contrary, for Q
elements, locking is less pronounced and the behavior of the element is improved when using SRI.