HEX8 — 8-node hexahedral solid ============================== **Kinematics.** Trilinear isoparametric 8-node hexahedron; three translational DOFs per node (:math:`u_x, u_y, u_z`); 24 DOFs per element. **Stiffness.** Displacement-based continuum formulation :math:`K_e = \int_{\Omega_e} B^\top C B \, \mathrm{d}V`, integrated on the reference cube :math:`[-1, 1]^3` with 2×2×2 Gauss-Legendre. **Mass.** Consistent :math:`M_e = \int_{\Omega_e} \rho N^\top N \, \mathrm{d}V`, same integration rule. Row-sum lumping available. **Formulations.** Default ``"full"`` is Hughes volumetric B-bar on top of 2×2×2 Gauss; ``"enhanced_strain"`` is Simo–Rifai 9-parameter EAS with Pian–Tong Jacobian scaling; ``"plain_gauss"`` is plain 2×2×2 Gauss for cross-checks. .. contents:: On this page :local: :depth: 2 ---- Kinematics and shape functions ------------------------------ The formulation below is the standard isoparametric continuum brick documented in every FEM textbook [Bathe2014]_ [Hughes1987]_ [ZTZ2013]_ [Cook2001]_; we give it here in a form that maps line-by-line to the implementation in :mod:`femorph_solver.elements.hex8`. On the reference cube with natural coordinates :math:`(\xi, \eta, \zeta) \in [-1, 1]^3`, trilinear shape functions .. math:: N_i(\xi, \eta, \zeta) = \tfrac{1}{8}\,(1 + \xi_i \xi)\,(1 + \eta_i \eta)\,(1 + \zeta_i \zeta) are associated with the eight corner nodes :math:`(\xi_i, \eta_i, \zeta_i) \in \{\pm 1\}^3`. The isoparametric map takes the reference cube to the physical hex: .. math:: x(\xi, \eta, \zeta) = \sum_{i=1}^{8} N_i(\xi, \eta, \zeta)\, x_i, \qquad u(\xi, \eta, \zeta) = \sum_{i=1}^{8} N_i(\xi, \eta, \zeta)\, u_i. The node-ordering convention matches VTK_HEXAHEDRON: corners of the bottom face (:math:`\zeta = -1`) first in CCW order, then the top face. Jacobian and strain-displacement matrix --------------------------------------- Derivatives of :math:`N_i` with respect to physical coordinates are recovered via the Jacobian of the isoparametric map, .. math:: \frac{\partial N_i}{\partial x_k} = \bigl(J^{-1}\bigr)_{kj}\, \frac{\partial N_i}{\partial \xi_j}, \qquad J_{jk} = \frac{\partial x_k}{\partial \xi_j} = \sum_{i=1}^{8} \frac{\partial N_i}{\partial \xi_j}\, x_{i,k}. Small-strain Voigt notation with engineering shear gives :math:`\varepsilon = [\varepsilon_{xx},\ \varepsilon_{yy},\ \varepsilon_{zz},\ \gamma_{xy},\ \gamma_{yz},\ \gamma_{xz}]^\top` and the 6 × 24 strain-displacement matrix .. math:: B_i = \begin{bmatrix} \partial_x N_i & 0 & 0 \\ 0 & \partial_y N_i & 0 \\ 0 & 0 & \partial_z N_i \\ \partial_y N_i & \partial_x N_i & 0 \\ 0 & \partial_z N_i & \partial_y N_i \\ \partial_z N_i & 0 & \partial_x N_i \end{bmatrix}. Stacking :math:`B = [B_1 \,|\, B_2 \,|\, \dots \,|\, B_8]` produces the full 6 × 24 block. Isotropic constitutive matrix ----------------------------- For linear isotropic elasticity with Lamé parameters :math:`\lambda = \frac{E\nu}{(1+\nu)(1-2\nu)}` and shear modulus :math:`\mu = G = \frac{E}{2(1+\nu)}`, .. math:: C = \begin{bmatrix} \lambda + 2\mu & \lambda & \lambda & 0 & 0 & 0 \\ \lambda & \lambda + 2\mu & \lambda & 0 & 0 & 0 \\ \lambda & \lambda & \lambda + 2\mu & 0 & 0 & 0 \\ 0 & 0 & 0 & \mu & 0 & 0 \\ 0 & 0 & 0 & 0 & \mu & 0 \\ 0 & 0 & 0 & 0 & 0 & \mu \end{bmatrix}. Gauss-Legendre quadrature on :math:`[-1, 1]^3` ---------------------------------------------- The 2 × 2 × 2 tensor-product rule uses 8 points at :math:`(\pm 1/\sqrt{3}, \pm 1/\sqrt{3}, \pm 1/\sqrt{3})` with unit weights. The rule is exact for polynomial integrands up to degree 3 in each direction — enough to integrate the :math:`B^\top C B` product exactly on undistorted hexes. Physical volume integration picks up the determinant of :math:`J`: .. math:: K_e = \sum_{g=1}^{8} B(\xi_g)^\top C \, B(\xi_g)\, \lvert J(\xi_g) \rvert\, w_g. Numpy walk-through (plain displacement formulation) --------------------------------------------------- The following pure-numpy listing assembles :math:`K_e` for a unit cube made of generic steel (E = 210 GPa, ν = 0.30). Running it on a 24-DOF element produces a symmetric matrix with six exact rigid-body zero eigenvalues and eighteen non-zero elastic modes. .. code-block:: python import numpy as np # 2x2x2 Gauss-Legendre on [-1, 1]^3. g = 1.0 / np.sqrt(3.0) gp = np.array([(-g,-g,-g), (+g,-g,-g), (+g,+g,-g), (-g,+g,-g), (-g,-g,+g), (+g,-g,+g), (+g,+g,+g), (-g,+g,+g)]) w = np.ones(8) # Natural coordinates of the 8 nodes (VTK_HEXAHEDRON order). xi_n = np.array([(-1,-1,-1), (+1,-1,-1), (+1,+1,-1), (-1,+1,-1), (-1,-1,+1), (+1,-1,+1), (+1,+1,+1), (-1,+1,+1)], dtype=float) def shape_and_derivs(xi, eta, zeta): N = 0.125 * (1 + xi_n[:, 0] * xi) \ * (1 + xi_n[:, 1] * eta) \ * (1 + xi_n[:, 2] * zeta) dN = np.empty((8, 3)) dN[:, 0] = 0.125 * xi_n[:, 0] * (1 + xi_n[:, 1] * eta) * (1 + xi_n[:, 2] * zeta) dN[:, 1] = 0.125 * xi_n[:, 1] * (1 + xi_n[:, 0] * xi) * (1 + xi_n[:, 2] * zeta) dN[:, 2] = 0.125 * xi_n[:, 2] * (1 + xi_n[:, 0] * xi) * (1 + xi_n[:, 1] * eta) return N, dN # Unit cube nodes. X = np.array([(0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1)], dtype=float) E, nu = 2.1e11, 0.30 lam = E*nu / ((1+nu)*(1-2*nu)) mu = E / (2*(1+nu)) C = np.diag([lam+2*mu, lam+2*mu, lam+2*mu, mu, mu, mu]).astype(float) C[:3, :3] += lam * (1 - np.eye(3)) # off-diagonal λ coupling of ε_ii K = np.zeros((24, 24)) for (xi, eta, zeta), wi in zip(gp, w): _, dNdxi = shape_and_derivs(xi, eta, zeta) J = dNdxi.T @ X detJ = np.linalg.det(J) dNdx = np.linalg.solve(J, dNdxi.T).T B = np.zeros((6, 24)) for i in range(8): col = 3 * i B[0, col + 0] = dNdx[i, 0] B[1, col + 1] = dNdx[i, 1] B[2, col + 2] = dNdx[i, 2] B[3, col + 0] = dNdx[i, 1] B[3, col + 1] = dNdx[i, 0] B[4, col + 1] = dNdx[i, 2] B[4, col + 2] = dNdx[i, 1] B[5, col + 0] = dNdx[i, 2] B[5, col + 2] = dNdx[i, 0] K += B.T @ C @ B * detJ * wi eig = np.sort(np.linalg.eigvalsh(K)) # eig[:6] ≈ 0 # six rigid-body modes # eig[6:] > 0 # 18 elastic modes Volumetric locking and Hughes B-bar ----------------------------------- The plain displacement formulation above suffers from **volumetric locking** as :math:`\nu \to 1/2` or when the strain field has a strong volumetric component that the trilinear kinematics cannot accommodate. Hughes's selective reduced integration [Hughes1980]_ projects the volumetric strain component onto the element-constant subspace: .. math:: \bar{B}_i = B_i + \tfrac{1}{3} \begin{bmatrix} \bar{b}_x - b_{i,x} & \bar{b}_y - b_{i,y} & \bar{b}_z - b_{i,z} \\ \bar{b}_x - b_{i,x} & \bar{b}_y - b_{i,y} & \bar{b}_z - b_{i,z} \\ \bar{b}_x - b_{i,x} & \bar{b}_y - b_{i,y} & \bar{b}_z - b_{i,z} \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}, where :math:`b_{i,k} = \partial_{x_k} N_i` and :math:`\bar{b}_k = \frac{1}{V_e} \int_{\Omega_e} b_{i,k} \, \mathrm{d}V` is the volume-averaged derivative. In practice the correction is implemented in closed form as .. math:: K_{\bar B} = K_\text{plain} + K_B \bigl(\tfrac{1}{V_e}\, S S^\top - H\bigr), \qquad K_B = \lambda + \tfrac{2}{3} G, with :math:`S = \sum_g b(\xi_g)\lvert J(\xi_g)\rvert w_g` and :math:`H = \sum_g b(\xi_g)b(\xi_g)^\top \lvert J(\xi_g)\rvert w_g`, where :math:`b(\xi_g)` is the 24-vector of volumetric derivatives at Gauss point :math:`\xi_g`. This is femorph-solver's :ref:`default "full" ` formulation. Enhanced assumed strain (Simo–Rifai) ------------------------------------ An alternative cure for volumetric *and* shear locking replaces the compatible strain :math:`\varepsilon = B u_e` with an enhanced decomposition .. math:: \tilde{\varepsilon} = B u_e + G(\xi)\, \alpha, where :math:`\alpha` is a set of per-element strain parameters and :math:`G(\xi)` is a 9-column enhancement matrix constructed from hierarchic bubble functions [SimoRifai1990]_. Static condensation eliminates :math:`\alpha` at the element level, producing the enhanced stiffness .. math:: \tilde{K}_e = K_e - K_{e\alpha} \bigl(K_{\alpha\alpha}\bigr)^{-1} K_{\alpha e}, with .. math:: K_e &= \int B^\top C B \,\mathrm{d}V, \\ K_{e\alpha} &= \int B^\top C \, \Gamma(\xi) \,\mathrm{d}V, \\ K_{\alpha\alpha} &= \int \Gamma(\xi)^\top C \, \Gamma(\xi) \,\mathrm{d}V, where :math:`\Gamma(\xi) = j_0 \, J_0^{-1}\, G(\xi)\, / \lvert J(\xi)\rvert` is the enhancement in current-configuration coordinates under the Pian–Tong :math:`j_0 / j` scaling. Reach the EAS path via ``ELEMENTS.HEX8(integration="enhanced_strain")`` on :meth:`~femorph_solver.Model.assign`; the kernel keys the choice on ``material["_HEX8_INTEGRATION"]`` internally. Mass matrix ----------- Consistent mass uses the same integration rule: .. math:: M_e = \int_{\Omega_e} \rho \, N^\top N \,\mathrm{d}V = \sum_{g=1}^{8} \rho \, N(\xi_g)^\top N(\xi_g) \, \lvert J(\xi_g)\rvert \, w_g. Row-sum lumping :math:`M_e^{\text{lump}} = \text{diag}\bigl(\sum_j M_{e,ij}\bigr)` is available via :class:`~femorph_solver.Model.mass_matrix` with ``lumped=True``. .. _hex8-formulations: Integration flag ---------------- .. list-table:: :widths: 25 75 :header-rows: 1 * - Name - Description * - ``"full"`` (default) - 2×2×2 Gauss + Hughes volumetric B-bar. Cures volumetric locking at near-incompressible Poisson ratios. * - ``"enhanced_strain"`` - Simo–Rifai EAS with 9 enhancement modes and Pian–Tong :math:`j_0 / j` scaling. Cures both shear and volumetric locking; recovers >95 % of the Euler–Bernoulli tip deflection on a slender cantilever with one element through the thickness (vs <50 % for plain 2×2×2). * - ``"plain_gauss"`` - Vanilla 2×2×2 Gauss with no B-bar / EAS treatment. Locks volumetrically at near-incompressible Poisson and shear- locks on bending-dominated geometries — useful only for cross-solver comparisons with codes that don't ship B-bar. Pick the formulation at :meth:`~femorph_solver.Model.assign` time via the :data:`~femorph_solver.ELEMENTS` namespace: .. code-block:: python import femorph_solver as fs from femorph_solver import ELEMENTS model = fs.Model.from_grid(grid) model.assign( ELEMENTS.HEX8(integration="enhanced_strain"), # or "full" (default, B-bar) {"EX": 2.1e11, "PRXY": 0.30}, ) Validation ---------- **Rigid-body modes.** On any undistorted hex, :math:`K_e` should have exactly six zero eigenvalues (three translations + three rotations) and eighteen non-zero elastic modes. Verify directly: .. code-block:: python import numpy as np from femorph_solver.elements import get HEX8 = get("HEX8") X = np.array([(0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1)], dtype=float) K = HEX8.ke(X, {"EX": 2.1e11, "PRXY": 0.30}, np.empty(0)) eig = np.sort(np.linalg.eigvalsh(K)) assert np.allclose(eig[:6], 0, atol=1e-4 * eig[-1]) # 6 rigid-body modes assert (eig[6:] > 0).all() # 18 elastic modes **Cantilever tip deflection.** For a slender cantilever (:math:`L \gg h`) with tip load :math:`P`, Euler–Bernoulli theory gives :math:`u_\text{tip} = P L^3 / (3 E I)`. A single hex through the thickness reaches this limit to within a few percent under the EAS formulation; the plain-displacement formulation severely under-predicts due to shear locking. The MacNeal–Harder benchmark set [MacNealHarder1985]_ ships the canonical variants of this problem; the :doc:`/gallery/elements/solid185/index` runs the slender-cantilever version end-to-end. **Cross-solver migration safety.** femorph-solver's HEX8 matches established commercial solvers (default B-bar formulation) on the bundled quarter-arc sector to machine precision — see the :doc:`HEX8 cyclic-modal migration check `. API reference ------------- .. autoclass:: femorph_solver.elements.hex8.Hex8 :members: References ---------- .. [Bathe2014] Bathe, K.-J. (2014). *Finite Element Procedures*, 2nd ed. Prentice-Hall / Klaus-Jürgen Bathe. Chapters 4 (isoparametric elements) and 6 (continuum). https://www.klausjurgenbathe.com/fepbook/ .. [Hughes1987] Hughes, T. J. R. (1987). *The Finite Element Method: Linear Static and Dynamic Finite Element Analysis*. Prentice-Hall; Dover reprint 2000 (ISBN 0-486-41181-8). Chapter 3 (isoparametric formulation), §4.5 (selective reduced integration). https://store.doverpublications.com/0486411818.html .. [Hughes1980] Hughes, T. J. R. (1980). "Generalization of selective integration procedures to anisotropic and nonlinear media." *International Journal for Numerical Methods in Engineering* **15**, 1413–1418. https://doi.org/10.1002/nme.1620150914 .. [SimoRifai1990] Simo, J. C. and Rifai, M. S. (1990). "A class of mixed assumed strain methods and the method of incompatible modes." *International Journal for Numerical Methods in Engineering* **29**, 1595–1638. https://doi.org/10.1002/nme.1620290802 .. [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 .. [ZTZ2013] Zienkiewicz, O. C., Taylor, R. L., and Zhu, J. Z. (2013). *The Finite Element Method: Its Basis and Fundamentals*, 7th ed. Butterworth-Heinemann (Elsevier). https://doi.org/10.1016/C2009-0-24909-9 .. [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. ISBN 978-0-471-35605-9. https://www.wiley.com/en-us/Concepts+and+Applications+of+Finite+Element+Analysis%2C+4th+Edition-p-9780471356059