{ "cells": [ { "cell_type": "markdown", "id": "bff1b49b", "metadata": {}, "source": [ "Tutorial: lumped mass matrices of the shells with an offset\n", "===\n", "\n", "Date: 24 of October 2026\n", "\n", "Author: Saullo G. P. Castro\n", "\n", "The lumped mass matrices of the shell elements `Quad4`, `Quad4R`, `Tria3R` and `Tria3DSG` (`mtype=1` of `update_M`) are derived here from the kinetic energy of the shell, as a nodal quadrature, and verified with the rigid-body inertia of one element and of a mesh. When the reference surface of the nodes is at an offset from the mid-surface, the first moment of the density through the thickness, `intrhoz`, couples the translations with the rotations. The derivation shows that this coupling stays in the lumped matrix, in the block of each node, and that it is what places the centre of mass of the element at the offset.\n", "\n", "References:\n", "\n", "- Brockman, R. A., Dynamics of the bilinear Mindlin plate element, *International Journal for Numerical Methods in Engineering*, 24(12), 2343-2356, 1987, [doi:10.1002/nme.1620241208](https://doi.org/10.1002/nme.1620241208): the kinetic energy, Eq. (20), the consistent mass matrix, Eq. (21), and the lumped mass matrix by Lobatto integration at the nodes, pp. 2348-2349;\n", "- Cowper, G. R., Gaussian quadrature formulas for triangles, *International Journal for Numerical Methods in Engineering*, 7(3), 405-408, 1973, [doi:10.1002/nme.1620070316](https://doi.org/10.1002/nme.1620070316): the three-point rule of the consistent matrix of `Tria3R`, first row of Table I.\n", "\n", "The helper functions are those of `tests/test_element_mass.py`, imported below." ] }, { "cell_type": "code", "execution_count": 1, "id": "03936d65", "metadata": { "execution": { "iopub.execute_input": "2026-10-25T09:19:31.643847Z", "iopub.status.busy": "2026-10-25T09:19:31.643684Z", "iopub.status.idle": "2026-10-25T09:19:32.452827Z", "shell.execute_reply": "2026-10-25T09:19:32.452160Z" } }, "outputs": [], "source": [ "import os\n", "import sys\n", "sys.path.insert(0, os.path.join('..', '..', 'tests'))\n", "\n", "import numpy as np\n", "import sympy as sp\n", "from scipy.sparse import coo_matrix\n", "\n", "from pyfe3d.shellprop_utils import isotropic_plate\n", "from pyfe3d import DOF, INT, DOUBLE\n", "from test_element_mass import (SHELLS, TRIANGLES, QUADS, area, area_rule,\n", " shell_element, rotation_blocks, dense_M,\n", " rigid_modes, shell_rigid_inertia, quad_N_detJ,\n", " mass_moments)" ] }, { "cell_type": "markdown", "id": "8894a9ab", "metadata": {}, "source": [ "Kinetic energy of the shell\n", "---\n", "\n", "With the displacements $u + z r_y$, $v - z r_x$ and $w$ at the distance $z$ from the reference surface, Eq. (1) of Brockman (1987) and the kinematics of the shells of `pyfe3d`, the kinetic energy per unit area, integrated through the thickness $z \\in [-h/2 + e, h/2 + e]$ of a reference surface at the offset $e$ from the mid-surface, is\n", "\n", "$$\n", "\\mathcal{T} = \\frac{1}{2} \\left[ I_\\rho (\\dot{u}^2 + \\dot{v}^2 + \\dot{w}^2) + 2 I_{\\rho z} (\\dot{u} \\dot{r}_y - \\dot{v} \\dot{r}_x) + I_{\\rho z^2} (\\dot{r}_x^2 + \\dot{r}_y^2) \\right] = \\frac{1}{2} \\dot{\\boldsymbol{d}}^\\top \\boldsymbol{B} \\dot{\\boldsymbol{d}}\n", "$$\n", "\n", "Eq. (20) of Brockman (1987), with $\\boldsymbol{d} = (u, v, w, r_x, r_y, r_z)$ and $I_\\rho$, $I_{\\rho z}$, $I_{\\rho z^2}$ the `intrho`, `intrhoz` and `intrhoz2` of `ShellProp`. The drilling rotation $r_z$ has no inertia. SymPy gives the three integrals and the matrix $\\boldsymbol{B}$, the inertia of the section:" ] }, { "cell_type": "code", "execution_count": 2, "id": "6dc5ae0e", "metadata": { "execution": { "iopub.execute_input": "2026-10-25T09:19:32.455439Z", "iopub.status.busy": "2026-10-25T09:19:32.455074Z", "iopub.status.idle": "2026-10-25T09:19:32.621714Z", "shell.execute_reply": "2026-10-25T09:19:32.620939Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "intrho = h*rho\n", "intrhoz = e*h*rho\n", "intrhoz2 = h*rho*(12*e**2 + h**2)/12\n" ] }, { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}I_{\\rho} & 0 & 0 & 0 & I_{rhoz} & 0\\\\0 & I_{\\rho} & 0 & - I_{rhoz} & 0 & 0\\\\0 & 0 & I_{\\rho} & 0 & 0 & 0\\\\0 & - I_{rhoz} & 0 & I_{rhoz2} & 0 & 0\\\\I_{rhoz} & 0 & 0 & 0 & I_{rhoz2} & 0\\\\0 & 0 & 0 & 0 & 0 & 0\\end{matrix}\\right]$" ], "text/plain": [ "Matrix([\n", "[ I_rho, 0, 0, 0, I_rhoz, 0],\n", "[ 0, I_rho, 0, -I_rhoz, 0, 0],\n", "[ 0, 0, I_rho, 0, 0, 0],\n", "[ 0, -I_rhoz, 0, I_rhoz2, 0, 0],\n", "[I_rhoz, 0, 0, 0, I_rhoz2, 0],\n", "[ 0, 0, 0, 0, 0, 0]])" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "z, h, e, rho = sp.symbols('z h e rho', real=True)\n", "ud, vd, wd, rxd, ryd, rzd = sp.symbols('udot vdot wdot rxdot rydot rzdot')\n", "lim = (z, -h/2 + e, h/2 + e)\n", "I0, I1, I2 = [sp.factor(sp.integrate(rho*z**k, lim)) for k in range(3)]\n", "print('intrho =', I0)\n", "print('intrhoz =', I1)\n", "print('intrhoz2 =', I2)\n", "T = sp.integrate(rho/2*((ud + z*ryd)**2 + (vd - z*rxd)**2 + wd**2), lim)\n", "d = sp.Matrix([ud, vd, wd, rxd, ryd, rzd])\n", "B = sp.hessian(T, d).applyfunc(sp.factor)\n", "Is, Isz, Isz2 = sp.symbols('I_rho I_rhoz I_rhoz2')\n", "B = B.subs({I0: Is, I1: Isz, I2: Isz2})\n", "B" ] }, { "cell_type": "markdown", "id": "1f46271c", "metadata": {}, "source": [ "The offset couples $u$ with $r_y$ and $v$ with $r_x$ through $I_{\\rho z} = \\rho h e$, which is zero only for $e = 0$ (or for a laminate whose density is symmetric about the reference surface).\n", "\n", "The consistent mass matrix\n", "---\n", "\n", "With the nodal velocities $\\dot{\\boldsymbol{q}}$ and $\\dot{\\boldsymbol{d}} = \\boldsymbol{S} \\dot{\\boldsymbol{q}}$, $\\boldsymbol{S} = [S_1 \\boldsymbol{I}_6, S_2 \\boldsymbol{I}_6, \\ldots]$, the kinetic energy of the element is $T = \\frac{1}{2} \\dot{\\boldsymbol{q}}^\\top \\boldsymbol{M} \\dot{\\boldsymbol{q}}$ with\n", "\n", "$$\n", "\\boldsymbol{M} = \\int_A \\boldsymbol{S}^\\top \\boldsymbol{B} \\boldsymbol{S} dA, \\qquad \\boldsymbol{M}_{ij} = h_{ij} \\boldsymbol{B}, \\qquad h_{ij} = \\int_A S_i S_j dA\n", "$$\n", "\n", "the form of Eq. (21) of Brockman (1987): the 6 x 6 block between the nodes $i$ and $j$ is the inertia of the section times $h_{ij}$. For the linear triangle in area coordinates, $h_{ij} = A/6$ for $i = j$ and $A/12$ otherwise, which the three-point rule of Cowper (1973) integrates exactly:" ] }, { "cell_type": "code", "execution_count": 3, "id": "dd3c7926", "metadata": { "execution": { "iopub.execute_input": "2026-10-25T09:19:32.623646Z", "iopub.status.busy": "2026-10-25T09:19:32.623452Z", "iopub.status.idle": "2026-10-25T09:19:32.701369Z", "shell.execute_reply": "2026-10-25T09:19:32.700746Z" } }, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}\\frac{A}{6} & \\frac{A}{12} & \\frac{A}{12}\\\\\\frac{A}{12} & \\frac{A}{6} & \\frac{A}{12}\\\\\\frac{A}{12} & \\frac{A}{12} & \\frac{A}{6}\\end{matrix}\\right]$" ], "text/plain": [ "Matrix([\n", "[ A/6, A/12, A/12],\n", "[A/12, A/6, A/12],\n", "[A/12, A/12, A/6]])" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "N1, N2, A = sp.symbols('N1 N2 A', positive=True)\n", "Ntri = sp.Matrix([N1, N2, 1 - N1 - N2])\n", "# NOTE the reference triangle has the area 1/2, the Jacobian determinant is 2A\n", "H_exact = 2*A*(Ntri*Ntri.T).applyfunc(\n", " lambda f: sp.integrate(f, (N2, 0, 1 - N1), (N1, 0, 1)))\n", "H_exact" ] }, { "cell_type": "markdown", "id": "2fdf096f", "metadata": {}, "source": [ "Lumped mass matrix by nodal quadrature\n", "---\n", "\n", "A quadrature with the points $\\boldsymbol{x}_k$ at the nodes and the weights $w_k$ integrates the whole integrand $\\boldsymbol{S}^\\top \\boldsymbol{B} \\boldsymbol{S}$ as\n", "\n", "$$\n", "\\boldsymbol{M} \\approx \\sum_k w_k \\boldsymbol{S}(\\boldsymbol{x}_k)^\\top \\boldsymbol{B} \\boldsymbol{S}(\\boldsymbol{x}_k)\n", "$$\n", "\n", "Since $S_i(\\boldsymbol{x}_k) = \\delta_{ik}$, each point contributes only to the block of its own node,\n", "\n", "$$\n", "\\boldsymbol{M}_{ij} = \\delta_{ij} w_i \\boldsymbol{B}\n", "$$\n", "\n", "The matrix is block-diagonal per node and each block is the complete inertia of the section, coupling included: lumping acts on the shape functions, never on $\\boldsymbol{B}$. The matrix is diagonal only when $I_{\\rho z} = 0$, which is the condition stated by Brockman (1987), p. 2349, for his lumped matrix of the quadrilateral by Lobatto integration at the nodes. Dropping the coupling with an offset is a further approximation, not a consequence of the lumping.\n", "\n", "The rules of `pyfe3d` are the vertex rule for the triangles, $w_i = A/3$, and the two-point Gauss-Lobatto rule in each direction for the quadrilaterals, $w_i = \\det J(\\xi_i, \\eta_i)$, the Jacobian determinant at node $i$ times the weights $1 \\times 1$. For the triangle, the vertex rule applied symbolically to the integrand gives:" ] }, { "cell_type": "code", "execution_count": 4, "id": "86c99f71", "metadata": { "execution": { "iopub.execute_input": "2026-10-25T09:19:32.703376Z", "iopub.status.busy": "2026-10-25T09:19:32.703217Z", "iopub.status.idle": "2026-10-25T09:19:32.714256Z", "shell.execute_reply": "2026-10-25T09:19:32.713663Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "blocks between different nodes are zero: True\n" ] }, { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}\\frac{A I_{\\rho}}{3} & 0 & 0 & 0 & \\frac{A I_{rhoz}}{3} & 0\\\\0 & \\frac{A I_{\\rho}}{3} & 0 & - \\frac{A I_{rhoz}}{3} & 0 & 0\\\\0 & 0 & \\frac{A I_{\\rho}}{3} & 0 & 0 & 0\\\\0 & - \\frac{A I_{rhoz}}{3} & 0 & \\frac{A I_{rhoz2}}{3} & 0 & 0\\\\\\frac{A I_{rhoz}}{3} & 0 & 0 & 0 & \\frac{A I_{rhoz2}}{3} & 0\\\\0 & 0 & 0 & 0 & 0 & 0\\end{matrix}\\right]$" ], "text/plain": [ "Matrix([\n", "[ A*I_rho/3, 0, 0, 0, A*I_rhoz/3, 0],\n", "[ 0, A*I_rho/3, 0, -A*I_rhoz/3, 0, 0],\n", "[ 0, 0, A*I_rho/3, 0, 0, 0],\n", "[ 0, -A*I_rhoz/3, 0, A*I_rhoz2/3, 0, 0],\n", "[A*I_rhoz/3, 0, 0, 0, A*I_rhoz2/3, 0],\n", "[ 0, 0, 0, 0, 0, 0]])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Bs = sp.Matrix(B)\n", "def integrand(Nvals):\n", " Nmat = sp.Matrix.hstack(*[Ni*sp.eye(6) for Ni in Nvals])\n", " return Nmat.T*Bs*Nmat\n", "\n", "vertices = [(1, 0, 0), (0, 1, 0), (0, 0, 1)]\n", "M_lump = sum((A/3*integrand(p) for p in vertices), sp.zeros(18, 18))\n", "print('blocks between different nodes are zero:',\n", " all(M_lump[6*i:6*i+6, 6*j:6*j+6] == sp.zeros(6, 6)\n", " for i in range(3) for j in range(3) if i != j))\n", "M_lump[0:6, 0:6]" ] }, { "cell_type": "markdown", "id": "543c73ee", "metadata": {}, "source": [ "Row-sum lumping gives the same blocks with $w_i = \\sum_j h_{ij} = \\int_A S_i dA$, since $\\sum_j \\boldsymbol{M}_{ij} = (\\sum_j h_{ij}) \\boldsymbol{B}$ and $\\sum_j S_j = 1$. For the triangle, $\\int_A S_i dA = A/3$ and the two lumpings coincide:" ] }, { "cell_type": "code", "execution_count": 5, "id": "e21f5fb7", "metadata": { "execution": { "iopub.execute_input": "2026-10-25T09:19:32.716806Z", "iopub.status.busy": "2026-10-25T09:19:32.716540Z", "iopub.status.idle": "2026-10-25T09:19:32.728980Z", "shell.execute_reply": "2026-10-25T09:19:32.728299Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "row-sum lumping equal to the vertex rule: True\n" ] } ], "source": [ "M_cons = sp.zeros(18, 18)\n", "for i in range(3):\n", " for j in range(3):\n", " M_cons[6*i:6*i+6, 6*j:6*j+6] = H_exact[i, j]*Bs\n", "row_sum = sp.zeros(18, 18)\n", "for i in range(3):\n", " row_sum[6*i:6*i+6, 6*i:6*i+6] = sum((M_cons[6*i:6*i+6, 6*j:6*j+6]\n", " for j in range(3)), sp.zeros(6, 6))\n", "print('row-sum lumping equal to the vertex rule:',\n", " sp.simplify(row_sum - M_lump) == sp.zeros(18, 18))" ] }, { "cell_type": "markdown", "id": "96c42fc6", "metadata": {}, "source": [ "Rigid-body inertia\n", "---\n", "\n", "The inertia of a rigid body, with the mass $m$, the first moment of the mass $\\boldsymbol{Q} = \\int \\rho \\boldsymbol{p} dV$ and the second moments $\\boldsymbol{I}_2 = \\int \\rho \\boldsymbol{p} \\boldsymbol{p}^\\top dV$, is the 6 x 6 matrix $\\boldsymbol{\\Phi}^\\top \\boldsymbol{M} \\boldsymbol{\\Phi}$, where the columns of $\\boldsymbol{\\Phi}$ are the nodal values of the three translations and of the three rotations about the origin. For a point $\\boldsymbol{p} = \\boldsymbol{x} + z \\boldsymbol{n}$ of the shell, with $\\boldsymbol{x}$ on the reference surface and $\\boldsymbol{n}$ its normal, the exact values are integrals over the area, and the nodal quadrature replaces them with sums:\n", "\n", "$$\n", "m = I_\\rho \\sum_i w_i, \\qquad\n", "\\boldsymbol{Q} = I_\\rho \\sum_i w_i \\boldsymbol{x}_i + I_{\\rho z} \\boldsymbol{n} \\sum_i w_i\n", "$$\n", "\n", "$$\n", "\\boldsymbol{I}_2 = I_\\rho \\sum_i w_i \\boldsymbol{x}_i \\boldsymbol{x}_i^\\top + I_{\\rho z} \\sum_i w_i (\\boldsymbol{x}_i \\boldsymbol{n}^\\top + \\boldsymbol{n} \\boldsymbol{x}_i^\\top) + I_{\\rho z^2} \\boldsymbol{n} \\boldsymbol{n}^\\top \\sum_i w_i\n", "$$\n", "\n", "Hence:\n", "\n", "- the mass is exact when $\\sum_i w_i = A$, true for both rules;\n", "- the first moment through the thickness, $I_{\\rho z} A \\boldsymbol{n}$, is exact under the same condition, and only with the coupling: without it, the centre of mass of the element is on the reference surface instead of at the offset;\n", "- the first moment in the plane, $I_\\rho \\sum_i w_i \\boldsymbol{x}_i$, is exact when $w_i = \\int_A S_i dA$, because $\\int_A \\boldsymbol{x} dA = \\sum_i \\boldsymbol{x}_i \\int_A S_i dA$ for an isoparametric element. This holds for the triangle and for the parallelogram, where $\\det J = A/4$ is constant, but not for a general quadrilateral with the Gauss-Lobatto weights;\n", "- of the second moments, the term $I_{\\rho z^2} A$ is exact, while $I_\\rho \\sum_i w_i \\boldsymbol{x}_i \\boldsymbol{x}_i^\\top$ concentrates the mass at the nodes. About the centroid of the element it is 4 times the exact value for a triangle and 3 times for a parallelogram: the inertia of a single element is overestimated, an error that decreases with the size of the elements, see the last section.\n", "\n", "The weights of the distorted quadrilateral of the tests, with the Gauss-Lobatto rule and with row-sum lumping:" ] }, { "cell_type": "code", "execution_count": 6, "id": "d395f6de", "metadata": { "execution": { "iopub.execute_input": "2026-10-25T09:19:32.731393Z", "iopub.status.busy": "2026-10-25T09:19:32.731223Z", "iopub.status.idle": "2026-10-25T09:19:32.740608Z", "shell.execute_reply": "2026-10-25T09:19:32.739879Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "A = 4.425\n", "Gauss-Lobatto w_i/A = [0.22316 0.19209 0.27684 0.30791] sum = 1.0000000000000002\n", "row-sum w_i/A = [0.24105 0.2307 0.25895 0.2693 ] sum = 1.0000000000000004\n", "centroid of the area [2.12806 1.94652]\n", "centre, Gauss-Lobatto [2.13418 2.08955]\n", "centre, row-sum [2.12806 1.94652]\n" ] } ], "source": [ "xy = QUADS['distorted']\n", "Aq = area(xy)\n", "nodes = [(-1, -1), (1, -1), (1, 1), (-1, 1)]\n", "w_lobatto = np.array([quad_N_detJ(xy, xi, eta)[1] for xi, eta in nodes])\n", "w_rowsum = sum(w*N for N, w in area_rule(xy, 'exact'))\n", "x_exact = sum(w*(N @ xy) for N, w in area_rule(xy, 'exact'))/Aq\n", "print('A =', Aq)\n", "print('Gauss-Lobatto w_i/A =', (w_lobatto/Aq).round(5), ' sum =', w_lobatto.sum()/Aq)\n", "print('row-sum w_i/A =', (w_rowsum/Aq).round(5), ' sum =', w_rowsum.sum()/Aq)\n", "print('centroid of the area ', x_exact.round(5))\n", "print('centre, Gauss-Lobatto ', (w_lobatto @ xy/Aq).round(5))\n", "print('centre, row-sum ', (w_rowsum @ xy/Aq).round(5))\n", "assert np.isclose(w_lobatto.sum(), Aq) and np.allclose(w_rowsum @ xy/Aq, x_exact)" ] }, { "cell_type": "markdown", "id": "01ec98a4", "metadata": {}, "source": [ "`pyfe3d` keeps the Gauss-Lobatto weights of Brockman (1987) for the quadrilaterals: on this element they move the centre of mass, in the plane, by 7 % of the square root of the area, while row-sum lumping keeps it at the centroid. Both give the same matrix for a parallelogram.\n", "\n", "Verification with the elements of pyfe3d\n", "---\n", "\n", "Each element is placed in a generic orientation in space, with the offset $e = 0.02$ of a plate of thickness $h = 0.05$, and $\\boldsymbol{\\Phi}^\\top \\boldsymbol{M} \\boldsymbol{\\Phi}$ of its lumped matrix is compared with the exact inertia of the shell. The comparison is repeated for the same matrix without the coupling terms. The columns are the relative error of the mass, the distance in the plane between the centre of mass and the centroid of the area, over the square root of the area, the distance of the centre of mass from the reference surface, exact $e = 0.02$, and the ratios of the inertia about the centre of mass, about the two axes of the element in its plane and about its normal, to the exact values." ] }, { "cell_type": "code", "execution_count": 7, "id": "55a308bd", "metadata": { "execution": { "iopub.execute_input": "2026-10-25T09:19:32.742982Z", "iopub.status.busy": "2026-10-25T09:19:32.742792Z", "iopub.status.idle": "2026-10-25T09:19:32.765986Z", "shell.execute_reply": "2026-10-25T09:19:32.765128Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "intrho 135, intrhoz 2.7, intrhoz2 0.08212\n", "element geometry coupling mass c plane c offset J x_e J y_e J n\n", "Tria3R right kept 2.2e-16 1.2e-16 0.02000 3.9983 3.9950 4.0000\n", "Tria3R right removed 2.2e-16 3.6e-16 -0.00000 3.9994 3.9982 4.0000\n", "Tria3R distorted kept 4.4e-16 2.4e-16 0.02000 3.9978 3.9981 4.0000\n", "Tria3R distorted removed 4.4e-16 4.2e-16 0.00000 3.9992 3.9993 4.0000\n", "Tria3DSG right kept 4.4e-16 2.8e-17 0.02000 3.9983 3.9950 4.0000\n", "Tria3DSG right removed 6.7e-16 3.2e-16 -0.00000 3.9994 3.9982 4.0000\n", "Tria3DSG distorted kept 4.4e-16 2.4e-16 0.02000 3.9978 3.9981 4.0000\n", "Tria3DSG distorted removed 4.4e-16 4.2e-16 0.00000 3.9992 3.9993 4.0000\n", "Quad4R rectangle kept 4.4e-16 8.6e-16 0.02000 2.9993 2.9978 3.0000\n", "Quad4R rectangle removed 6.7e-16 1.1e-15 -0.00000 3.0000 2.9999 3.0000\n", "Quad4R parallelogram kept 0.0e+00 6.5e-17 0.02000 2.9974 2.9989 3.0000\n", "Quad4R parallelogram removed 0.0e+00 5.6e-16 -0.00000 2.9999 3.0000 3.0000\n", "Quad4R distorted kept 2.2e-16 6.8e-02 0.02000 2.9647 3.0933 3.0139\n", "Quad4R distorted removed 2.2e-16 6.8e-02 -0.00000 2.9655 3.0946 3.0139\n", "Quad4 rectangle kept 4.4e-16 8.6e-16 0.02000 2.9993 2.9978 3.0000\n", "Quad4 rectangle removed 6.7e-16 1.1e-15 -0.00000 3.0000 2.9999 3.0000\n", "Quad4 parallelogram kept 0.0e+00 6.5e-17 0.02000 2.9974 2.9989 3.0000\n", "Quad4 parallelogram removed 0.0e+00 5.6e-16 -0.00000 2.9999 3.0000 3.0000\n", "Quad4 distorted kept 2.2e-16 6.8e-02 0.02000 2.9647 3.0933 3.0139\n", "Quad4 distorted removed 2.2e-16 6.8e-02 -0.00000 2.9655 3.0946 3.0139\n" ] } ], "source": [ "prop = isotropic_plate(thickness=0.05, E=70e9, nu=0.3, rho=2700., offset=0.02)\n", "print('intrho %.4g, intrhoz %.4g, intrhoz2 %.4g' % (prop.intrho, prop.intrhoz, prop.intrhoz2))\n", "\n", "\n", "def without_coupling(M, Rb, nnodes):\n", " Ml = Rb.T @ M @ Rb\n", " for i in range(nnodes):\n", " for a, b in ((0, 4), (1, 3)):\n", " Ml[6*i + a, 6*i + b] = Ml[6*i + b, 6*i + a] = 0.\n", " return Rb @ Ml @ Rb.T\n", "\n", "\n", "print('%-9s %-14s %-8s %9s %9s %9s %8s %8s %8s'\n", " % ('element', 'geometry', 'coupling', 'mass', 'c plane', 'c offset',\n", " 'J x_e', 'J y_e', 'J n'))\n", "for name in SHELLS:\n", " geoms = TRIANGLES if name.startswith('Tria') else QUADS\n", " for gname, xy in geoms.items():\n", " nnodes = xy.shape[0]\n", " e, data, Xg = shell_element(name, xy)\n", " r, Rb = rotation_blocks(e, nnodes)\n", " n = r[:, 2]\n", " Phi = rigid_modes(Xg)\n", " m0, S0, J0 = mass_moments(shell_rigid_inertia(Xg, n, prop, area_rule(xy, 'exact')))\n", " M = dense_M(e, data, prop, 1, nnodes)\n", " for label, Mk in (('kept', M), ('removed', without_coupling(M, Rb, nnodes))):\n", " m, S, Jc = mass_moments(Phi.T @ Mk @ Phi)\n", " dc = S/m - S0/m0\n", " c_plane = np.linalg.norm(dc - (dc @ n)*n)/np.sqrt(area(xy))\n", " c_offset = (S/m - Xg[0]) @ n\n", " J = np.diag(r.T @ Jc @ r)/np.diag(r.T @ J0 @ r)\n", " print('%-9s %-14s %-8s %9.1e %9.1e %9.5f %8.4f %8.4f %8.4f'\n", " % (name, gname, label, abs(m/m0 - 1), c_plane, c_offset, *J))\n", " assert abs(m/m0 - 1) < 1e-13\n", " if name.startswith('Tria') or gname != 'distorted':\n", " assert c_plane < 1e-13\n", " if label == 'kept':\n", " assert np.isclose(c_offset, 0.02, rtol=1e-12)\n", " if gname != 'distorted' or name.startswith('Tria'):\n", " assert np.isclose(J[2], 4. if name.startswith('Tria') else 3.)\n", " else:\n", " assert abs(c_offset) < 1e-12" ] }, { "cell_type": "markdown", "id": "2af24f64", "metadata": {}, "source": [ "With the coupling, the mass is exact and the centre of mass is at the offset $e = 0.02$ for every element; in the plane it is at the centroid of the area for all but the distorted quadrilateral, as derived above. Without the coupling, the centre of mass of every element is on the reference surface. The inertia about the normal through the centre of mass is 4 times the exact value for the triangles and 3 times for the rectangle and the parallelogram. About the axes in the plane the ratio is slightly smaller, because the exact term $I_{\\rho z^2} A$ is part of the inertia, and it differs slightly without the coupling, because the centre of mass, about which it is taken, moves.\n", "\n", "Rotary inertia of a mesh\n", "---\n", "\n", "The error of the rotary inertia belongs to each element, about its own centroid. The inertia of a structure, about its own centre of mass, converges with the mesh. A rectangular plate $a \\times b$, in the plane $z = 0$, with the same property, has the exact inertia about its centre of mass\n", "\n", "$$\n", "J_{xx} = m \\frac{b^2}{12} + \\rho \\frac{h^3}{12} a b, \\qquad\n", "J_{yy} = m \\frac{a^2}{12} + \\rho \\frac{h^3}{12} a b, \\qquad\n", "J_{zz} = m \\frac{a^2 + b^2}{12}, \\qquad m = \\rho h a b\n", "$$\n", "\n", "independent of the offset $e$, since the centre of mass is at $z = e$. The plate is meshed with $n \\times n$ quadrilaterals, or with the triangles of the same mesh, and the lumped matrices are assembled:" ] }, { "cell_type": "code", "execution_count": 8, "id": "326514c6", "metadata": { "execution": { "iopub.execute_input": "2026-10-25T09:19:32.768427Z", "iopub.status.busy": "2026-10-25T09:19:32.768207Z", "iopub.status.idle": "2026-10-25T09:19:33.055410Z", "shell.execute_reply": "2026-10-25T09:19:33.054489Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "element n J_xx error J_yy error J_zz error\n", "Tria3R 1 1.986e+00 1.995e+00 2.000e+00\n", "Tria3R 2 4.966e-01 4.988e-01 5.000e-01\n", "Tria3R 4 1.241e-01 1.247e-01 1.250e-01\n", "Tria3R 8 3.103e-02 3.117e-02 3.125e-02\n", "Tria3R 16 7.759e-03 7.793e-03 7.813e-03\n", "Tria3R 32 1.940e-03 1.948e-03 1.953e-03\n", "Tria3R order of convergence between n = 16 and 32: [2. 2. 2.]\n", "Tria3DSG 1 1.986e+00 1.995e+00 2.000e+00\n", "Tria3DSG 2 4.966e-01 4.988e-01 5.000e-01\n", "Tria3DSG 4 1.241e-01 1.247e-01 1.250e-01\n", "Tria3DSG 8 3.103e-02 3.117e-02 3.125e-02\n", "Tria3DSG 16 7.759e-03 7.793e-03 7.813e-03\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Tria3DSG 32 1.940e-03 1.948e-03 1.953e-03\n", "Tria3DSG order of convergence between n = 16 and 32: [2. 2. 2.]\n", "Quad4R 1 1.986e+00 1.995e+00 2.000e+00\n", "Quad4R 2 4.966e-01 4.988e-01 5.000e-01\n", "Quad4R 4 1.241e-01 1.247e-01 1.250e-01\n", "Quad4R 8 3.103e-02 3.117e-02 3.125e-02\n", "Quad4R 16 7.759e-03 7.793e-03 7.813e-03\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Quad4R 32 1.940e-03 1.948e-03 1.953e-03\n", "Quad4R order of convergence between n = 16 and 32: [2. 2. 2.]\n", "Quad4 1 1.986e+00 1.995e+00 2.000e+00\n", "Quad4 2 4.966e-01 4.988e-01 5.000e-01\n", "Quad4 4 1.241e-01 1.247e-01 1.250e-01\n", "Quad4 8 3.103e-02 3.117e-02 3.125e-02\n", "Quad4 16 7.759e-03 7.793e-03 7.813e-03\n", "Quad4 32 1.940e-03 1.948e-03 1.953e-03\n", "Quad4 order of convergence between n = 16 and 32: [2. 2. 2.]\n" ] } ], "source": [ "a, b = 1.0, 0.6\n", "rho_, h_ = 2700., 0.05\n", "m_plate = rho_*h_*a*b\n", "J_ref = np.array([m_plate*b**2/12 + rho_*h_**3/12*a*b,\n", " m_plate*a**2/12 + rho_*h_**3/12*a*b,\n", " m_plate*(a**2 + b**2)/12])\n", "\n", "\n", "def plate_lumped_inertia(name, n):\n", " Elem, Data, Probe, _ = SHELLS[name]\n", " xs, ys = np.meshgrid(np.linspace(0, a, n + 1), np.linspace(0, b, n + 1))\n", " X = np.column_stack((xs.ravel(), ys.ravel(), np.zeros(xs.size)))\n", " nid = np.arange(xs.size).reshape(n + 1, n + 1)\n", " quads = [(nid[j, i], nid[j, i+1], nid[j+1, i+1], nid[j+1, i])\n", " for j in range(n) for i in range(n)]\n", " conn = quads\n", " if name.startswith('Tria'):\n", " conn = [t for (n1, n2, n3, n4) in quads for t in ((n1, n2, n3), (n1, n3, n4))]\n", " data, probe = Data(), Probe()\n", " size = data.M_SPARSE_SIZE*len(conn)\n", " Mr = np.zeros(size, dtype=INT)\n", " Mc = np.zeros(size, dtype=INT)\n", " Mv = np.zeros(size, dtype=DOUBLE)\n", " x = X.ravel()\n", " for k, nodes in enumerate(conn):\n", " el = Elem(probe)\n", " for i, ni in enumerate(nodes):\n", " setattr(el, 'n%d' % (i + 1), ni)\n", " setattr(el, 'c%d' % (i + 1), DOF*ni)\n", " el.init_k_M = k*data.M_SPARSE_SIZE\n", " el.update_rotation_matrix(x)\n", " el.update_probe_xe(x)\n", " el.update_M(Mr, Mc, Mv, prop, mtype=1)\n", " N = DOF*X.shape[0]\n", " M = coo_matrix((Mv, (Mr, Mc)), shape=(N, N)).tocsr()\n", " Phi = rigid_modes(X)\n", " m, S, Jc = mass_moments(Phi.T @ (M @ Phi))\n", " assert np.isclose(m, m_plate, rtol=1e-12)\n", " # centre of mass at the centre of the plate and at the offset\n", " assert np.allclose(S/m, [a/2, b/2, 0.02], rtol=1e-12)\n", " return np.diag(Jc)\n", "\n", "\n", "print('%-9s %4s %12s %12s %12s' % ('element', 'n', 'J_xx error', 'J_yy error', 'J_zz error'))\n", "for name in SHELLS:\n", " errs = []\n", " for n in (1, 2, 4, 8, 16, 32):\n", " err = plate_lumped_inertia(name, n)/J_ref - 1\n", " errs.append(err)\n", " print('%-9s %4d %12.3e %12.3e %12.3e' % (name, n, *err))\n", " order = np.log(errs[-2]/errs[-1])/np.log(2)\n", " print('%-9s order of convergence between n = 16 and 32: %s' % (name, order.round(3)))\n", " assert np.all(np.abs(order - 2) < 0.01)" ] }, { "cell_type": "markdown", "id": "ff98a9dc", "metadata": {}, "source": [ "The centre of mass of the meshed plate is at the offset, and the error of its rotary inertia decreases with the square of the element size." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.13" } }, "nbformat": 4, "nbformat_minor": 5 }