QUAD4_SHELL — 4-node Mindlin–Reissner shell#

MAPDL-compatible alias: SHELL181.

Kinematics. Flat four-node first-order shear-deformation shell (Mindlin–Reissner theory). Six DOFs per node (\(u_x, u_y, u_z, \theta_x, \theta_y, \theta_z\)); 24 DOFs per element in the local frame.

Stiffness. \(K_e = K_m + K_b + K_s + K_\text{drill}\) — membrane + bending + transverse shear + small drilling stabilisation.

Mass. Consistent \(M_e\) includes both translational and (through-thickness) rotational inertia.

MAPDL compatibility. Matches SHELL181 at default settings (flat shell, selective-reduced transverse shear, drilling-DOF stabilisation).

Mindlin–Reissner kinematics#

First-order shear-deformation theory [Hughes1987] [Bathe2014] describes the displacement of a point at through-thickness coordinate \(z\) in terms of the mid-surface displacement \((u, v, w)\) and two rotation DOFs \((\theta_x, \theta_y)\):

\[\begin{split}u_z(x, y, z) = u(x, y) + z\,\theta_y(x, y), \\ v_z(x, y, z) = v(x, y) - z\,\theta_x(x, y), \\ w_z(x, y, z) = w(x, y).\end{split}\]

The transverse normal remains straight but not necessarily perpendicular to the deformed mid-surface — the angle between the normal and the mid-surface is the transverse shear strain, captured by

\[\gamma_{xz} = \partial_x w + \theta_y, \qquad \gamma_{yz} = \partial_y w - \theta_x.\]

The in-plane strains split into membrane (\(z^0\)) and bending (\(z^1\)) contributions after through-thickness integration:

\[\begin{split}\varepsilon_m = \begin{bmatrix} \partial_x u \\ \partial_y v \\ \partial_x v + \partial_y u \end{bmatrix}, \qquad \varepsilon_b = \begin{bmatrix} \partial_x \theta_y \\ -\partial_y \theta_x \\ \partial_x \theta_x - \partial_y \theta_y \end{bmatrix}.\end{split}\]

Local frame and bilinear interpolation#

The element is assembled in a local 2-D frame with \(\hat{x}\) along edge 1→2, \(\hat{z}\) the best-fit normal, and \(\hat{y} = \hat{z} \times \hat{x}\). Nodes project to this plane and the element uses standard bilinear shape functions on \((\xi, \eta) \in [-1, 1]^2\):

\[N_i(\xi, \eta) = \tfrac{1}{4}\,(1 + \xi_i \xi)\,(1 + \eta_i \eta),\]

with corner nodes at \(\{\pm 1\}^2\) in CCW order — same as PLANE182 / VTK_QUAD.

Selective reduced integration (membrane, bending, shear)#

Locking management splits the stiffness into three integrals with different Gauss rules [Hughes1987] [ZTZ2013]:

  • Membrane: 2 × 2 Gauss (full).

  • Bending: 2 × 2 Gauss (full).

  • Transverse shear: 1 × 1 Gauss (reduced).

The reduced transverse-shear integration is the classical cure for shear locking on thin shells: full integration over-stiffens the element dramatically as the thickness-to-span ratio \(t/L \to 0\). MITC4 [BatheDvorkin1985] is the alternative cure (via tying of the shear strain at specific points); femorph-solver uses selective reduced integration, which matches MAPDL’s default SHELL181 behaviour.

Through-thickness integration uses the Reissner–Mindlin shear correction factor \(\kappa = 5/6\) on the transverse-shear stiffness.

Plane-stress constitutive#

SHELL181 runs in the plane-stress regime (\(\sigma_{zz} = 0\) in the local frame). Membrane and bending share the 3 × 3 plane-stress elastic matrix:

\[\begin{split}C_\text{ps} = \frac{E}{1 - \nu^2} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & (1 - \nu) / 2 \end{bmatrix}.\end{split}\]

Through-thickness constants:

\[D_m = C_\text{ps}\, t, \qquad D_b = \frac{t^3}{12}\, C_\text{ps}, \qquad D_s = \kappa\, G\, t\, I_2,\]

with \(G = E / (2 (1 + \nu))\).

Drilling-DOF stabilisation#

A flat plate has no physical stiffness for the local \(\theta_z\) rotation — the drilling DOF. Without stabilisation, the 24 × 24 local \(K\) has four extra zero-energy modes (one per corner) and global assembly can become singular. femorph-solver adds a small penalty term

\[K_\text{drill} = \alpha\, G\, t\, A_e \, \text{diag}_i\!\bigl(\delta_{i, \theta_z^\text{local}}\bigr),\]

with \(\alpha \approx 10^{-3}\) — enough to stabilise, small enough not to bleed into physical bending / membrane stiffness.

A block-diagonal 3 × 3 direction-cosine matrix then rotates every node’s 6-DOF block from the local frame back to global.

Numpy walk-through (skeleton)#

import numpy as np

# Bilinear shape functions on the reference quad.
xi_n = np.array([(-1,-1), (+1,-1), (+1,+1), (-1,+1)], dtype=float)

def shape_and_derivs(xi, eta):
    N  = 0.25 * (1 + xi_n[:, 0] * xi) * (1 + xi_n[:, 1] * eta)
    dN = np.empty((4, 2))
    dN[:, 0] = 0.25 * xi_n[:, 0] * (1 + xi_n[:, 1] * eta)
    dN[:, 1] = 0.25 * xi_n[:, 1] * (1 + xi_n[:, 0] * xi)
    return N, dN

# Full 2x2 for membrane + bending, single centre point for shear.
g = 1.0 / np.sqrt(3.0)
gp_full  = np.array([(-g,-g), (+g,-g), (+g,+g), (-g,+g)])
w_full   = np.ones(4)
gp_shear = np.array([(0.0, 0.0)])
w_shear  = np.array([4.0])             # full square area in natural coords

# Assemble K_m and K_b at gp_full with (D_m, D_b);
# assemble K_s at gp_shear with D_s; sum the three;
# add K_drill; rotate the 24x24 matrix to global.

Real constants and KEYOPT#

Slot

Meaning

real[0]

Through-thickness dimension \(t\) (mandatory).

The broader MAPDL SHELL181 KEYOPT table (composite lay-ups, shell offsets, incompatible-modes extension) is not yet exposed; the present implementation targets KEYOPT(3) = 0 (flat, uniform thickness).

Validation#

Rigid-body modes. 6 rigid-body zeros + 18 elastic eigenvalues on an undistorted square shell.

Bending patch. A slender cantilever strip (length \(L\), width \(b\), thickness \(t\)) with tip moment \(M\) should deflect to \(u_\text{tip} = M L^2 / (2 E I)\) with \(I = b t^3 / 12\). One-element-wide mesh reaches this limit to within a few percent via selective-reduced shear; plain full integration would be off by orders of magnitude on the thin-shell limit.

MacNeal–Harder benchmarks. The twisted-beam and pinched-cylinder problems from [MacNealHarder1985] are the community-standard shell stress tests; they’re on the roadmap for the verification gallery.

API reference#

class femorph_solver.elements.shell181.SHELL181[source]#

Bases: ElementBase

static ke(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray[source]#

Return the element stiffness matrix in global DOF ordering.

Parameters:

  • coords ((n_nodes, 3) float64)

  • material (mapping with MAPDL property names as keys (EX, PRXY, DENS, ...))

  • real (1-D array of real constants)

static me(coords: ndarray, material: dict[str, float], real: ndarray, lumped: bool = False) → ndarray[source]#

Return the element mass matrix in global DOF ordering.

References#

[Hughes1987] (1,2)

Hughes, T. J. R. (1987 / Dover 2000). The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Ch. 5 (Mindlin plates / shells), §5.2 (selective integration). https://store.doverpublications.com/0486411818.html

[Bathe2014]

Bathe, K.-J. (2014). Finite Element Procedures, 2nd ed. Ch. 5 & 6 on plates and shells. https://www.klausjurgenbathe.com/fepbook/

[ZTZ2013]

Zienkiewicz, O. C., Taylor, R. L., and Zhu, J. Z. (2013). The Finite Element Method: Its Basis and Fundamentals, 7th ed. Butterworth-Heinemann. https://doi.org/10.1016/C2009-0-24909-9

[BatheDvorkin1985]

Bathe, K.-J. and Dvorkin, E. N. (1985). “A four-node plate bending element based on Mindlin/Reissner plate theory and a mixed interpolation.” International Journal for Numerical Methods in Engineering 21 (2), 367–383. https://doi.org/10.1002/nme.1620210213

[MacNealHarder1985]

MacNeal, R. H. and Harder, R. L. (1985). “A proposed standard set of problems to test finite element accuracy.” Finite Elements in Analysis and Design 1 (1), 3–20. https://doi.org/10.1016/0168-874X(85)90003-4