QUAD4_PLANE — 2-D 4-node structural solid#

MAPDL-compatible alias: PLANE182.

Kinematics. Bilinear 2-D quadrilateral; two translational DOFs per node (\(u_x, u_y\)); 8 DOFs per element. Lives in the global (x, y) plane.

Stiffness. \(K_e = t \int_A B^\top C\, B\,\mathrm{d}A\) integrated on the reference square with 2 × 2 Gauss.

Mass. Consistent \(M_e = \rho\, t \int_A N^\top N\,\mathrm{d}A\) with 2 × 2 Gauss.

MAPDL compatibility. Reproduces PLANE182 with KEYOPT(3) = 0 (plane stress, default) and KEYOPT(3) = 2 (plane strain). KEYOPT(3) = 1 (axisymmetric) is not yet implemented.

Bilinear shape functions#

On the reference square \((\xi, \eta) \in [-1, 1]^2\), the four corner nodes at \((\xi_i, \eta_i) \in \{\pm 1\}^2\) in CCW order carry the bilinear shape functions [Cook2001] [ZTZ2013]

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

Node order matches VTK_QUAD / MAPDL — the same \((-,-), (+,-), (+,+), (-,+)\) sweep used by SHELL181.

Jacobian, \(B\) matrix, and Voigt strain#

Two-dimensional small strain with engineering shear: \(\varepsilon = [\varepsilon_{xx},\ \varepsilon_{yy},\ \gamma_{xy}]^\top\). The 3 × 8 strain-displacement matrix assembles from two-row blocks per node:

\[\begin{split}B_i = \begin{bmatrix} \partial_x N_i & 0 \\ 0 & \partial_y N_i \\ \partial_y N_i & \partial_x N_i \end{bmatrix}.\end{split}\]

Chain-rule derivatives come from the 2 × 2 Jacobian \(J_{jk} = \partial x_k / \partial \xi_j\).

Plane-stress vs plane-strain constitutive#

KEYOPT(3) selects:

  • Plane stress (KEYOPT(3) = 0, _PLANE_MODE = "plane_stress"): thin in \(\hat{z}\), \(\sigma_{zz} = 0\).

    \[\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}\]

  • Plane strain (KEYOPT(3) = 2, _PLANE_MODE = "plane_strain"): long in \(\hat{z}\), \(\varepsilon_{zz} = 0\).

    \[\begin{split}C_\text{pe} = \frac{E (1 - \nu)}{(1 + \nu)(1 - 2\nu)} \begin{bmatrix} 1 & \tfrac{\nu}{1 - \nu} & 0 \\ \tfrac{\nu}{1 - \nu} & 1 & 0 \\ 0 & 0 & \tfrac{1 - 2\nu}{2 (1 - \nu)} \end{bmatrix}.\end{split}\]

2 × 2 Gauss quadrature#

The 2 × 2 rule with points at \((\pm 1/\sqrt{3}, \pm 1/\sqrt{3})\) and unit weights integrates bilinear \(B^\top C B\) exactly on undistorted quads.

\[K_e = t \sum_{g=1}^{4} B(\xi_g)^\top C\, B(\xi_g)\,\lvert J(\xi_g)\rvert\, w_g, \qquad M_e = \rho\, t \sum_{g=1}^{4} N(\xi_g)^\top N(\xi_g)\,\lvert J(\xi_g)\rvert\, w_g.\]

Thickness \(t\) comes from the element’s real-constant slot real[0] (default 1.0 for plane-strain use).

Numpy walk-through#

import numpy as np

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

g = 1.0 / np.sqrt(3.0)
gp = np.array([(-g,-g), (+g,-g), (+g,+g), (-g,+g)])
w  = np.ones(4)

# Plane stress
E, nu = 2.1e11, 0.30
C = (E / (1 - nu**2)) * np.array([[1, nu, 0], [nu, 1, 0], [0, 0, (1 - nu) / 2]])

# Unit-square nodes in the (x, y) plane.
X = np.array([(0, 0), (1, 0), (1, 1), (0, 1)], dtype=float)

K = np.zeros((8, 8))
thick = 1.0
for (xi, eta), wi in zip(gp, w):
    _, dNdxi = shape_and_derivs(xi, eta)
    J    = dNdxi.T @ X
    detJ = np.linalg.det(J)
    dNdx = np.linalg.solve(J, dNdxi.T).T
    B = np.zeros((3, 8))
    for i in range(4):
        B[0, 2 * i    ] = dNdx[i, 0]
        B[1, 2 * i + 1] = dNdx[i, 1]
        B[2, 2 * i    ] = dNdx[i, 1]
        B[2, 2 * i + 1] = dNdx[i, 0]
    K += thick * (B.T @ C @ B) * detJ * wi

\(K_e\) has 3 rigid-body zeros (two translations + one rotation in the plane) and 5 elastic modes on an undistorted unit square.

Formulation flag#

Name

MAPDL KEYOPT(3)

Description

"plane_stress" (default)

0

Thin in \(\hat{z}\). Default thickness = 1.0.

"plane_strain"

2

Long in \(\hat{z}\). \(\varepsilon_{zz} = 0\).

Axisymmetric (KEYOPT(3) = 1) is not yet in the present kernel.

Validation#

Rigid-body modes. 3 rigid-body zero eigenvalues + 5 elastic modes on an undistorted square.

Patch test. A linear displacement field ramps exactly through every node of any distorted quad mesh — verified in tests.elements.plane182.

Convergence. Mesh refinement drives in-plane cantilever deflection to the plane-stress closed form under both "plane_stress" and "plane_strain" (with the appropriate analytic reference per mode).

API reference#

class femorph_solver.elements.plane182.PLANE182[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#

[Cook2001]

Cook, R. D., Malkus, D. S., Plesha, M. E., and Witt, R. J. (2001). Concepts and Applications of Finite Element Analysis, 4th ed. Wiley. Ch. 3 (plane stress / plane strain). https://www.wiley.com/en-us/Concepts+and+Applications+of+Finite+Element+Analysis%2C+4th+Edition-p-9780471356059

[ZTZ2013]

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

[Bathe2014]

Bathe, K.-J. (2014). Finite Element Procedures, 2nd ed. https://www.klausjurgenbathe.com/fepbook/