r""" Modal participation factors + effective modal mass ================================================== For base-excitation analysis (seismic, shock, hard-mounted random vibration), the **modal participation factor** :math:`\Gamma_{i}` and **effective modal mass** :math:`M_{\mathrm{eff},i}` of each mode tell you which modes will respond when the structure is shaken in a given direction — and which can safely be ignored. Together they answer the practical engineering question: *How many modes do I need to capture 90 % of the response in the y direction?* Definitions (Bathe 2014 §9.3.4; Chopra 2017 §13.1): .. math:: :label: modal-participation \Gamma_{i}^{(d)} \;=\; \frac{\phi_{i}^{T}\, M\, r^{(d)}}{\phi_{i}^{T}\, M\, \phi_{i}}, \qquad M_{\mathrm{eff},i}^{(d)} \;=\; \bigl(\Gamma_{i}^{(d)}\bigr)^{2}\, \phi_{i}^{T}\, M\, \phi_{i}, with :math:`\phi_{i}` the i-th mode shape (``mode_shapes[:, i]``), :math:`M` the consistent mass matrix, and :math:`r^{(d)}` the **influence vector** for direction :math:`d` — a vector of ones on every DOF in direction :math:`d` and zeros elsewhere. When mode shapes are mass-normalised (:math:`\phi^{T} M \phi = I`), the formulas simplify to :math:`\Gamma_{i}^{(d)} = \phi_{i}^{T} M r^{(d)}` and :math:`M_{\mathrm{eff},i}^{(d)} = \Gamma_{i}^{(d) 2}`. Total-mass property: .. math:: \sum_{i=1}^{N_{\mathrm{dof}}} M_{\mathrm{eff},i}^{(d)} \;=\; M_{\mathrm{total}}\quad \text{(every direction; sum over all modes)}. So the running cumulative sum of :math:`M_{\mathrm{eff},i}` gives a clean diagnostic — the first :math:`k` modes capture :math:`\sum_{i \le k} M_{\mathrm{eff},i}^{(d)} / M_{\mathrm{total}}` of the structure's response in direction :math:`d`. Implementation -------------- A slender HEX8-EAS prismatic cantilever (the :class:`~femorph_solver.validation.problems.CantileverHigherModes` geometry). Solve 12 modes, get :math:`M` from :meth:`Model.mass_matrix`, build the X / Y / Z influence vectors, and compute :math:`\Gamma` + :math:`M_{\mathrm{eff}}` for every mode in every direction. Print a participation-factor table sorted by mode number plus the cumulative-mass curve, and render the cumulative bar chart. References ---------- * Bathe, K.-J. (2014) *Finite Element Procedures*, 2nd ed., Prentice Hall, §9.3.4 (mode-superposition + participation factors). * Chopra, A. K. (2017) *Dynamics of Structures*, 5th ed., Pearson, §13.1 — modal effective mass. * Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002) *Concepts and Applications of Finite Element Analysis*, 4th ed., Wiley, §11.6 — base excitation. """ from __future__ import annotations import numpy as np import pyvista as pv import femorph_solver from femorph_solver import ELEMENTS # %% # Build a slender cantilever beam (HEX8 EAS) # ------------------------------------------ E = 2.0e11 NU = 0.30 RHO = 7850.0 L = 4.0 WIDTH = 0.05 HEIGHT = 0.05 NX, NY, NZ = 60, 3, 3 xs = np.linspace(0.0, L, NX + 1) ys = np.linspace(0.0, WIDTH, NY + 1) zs = np.linspace(0.0, HEIGHT, NZ + 1) grid = pv.StructuredGrid(*np.meshgrid(xs, ys, zs, indexing="ij")).cast_to_unstructured_grid() m = femorph_solver.Model.from_grid(grid) m.assign( ELEMENTS.HEX8(integration="enhanced_strain"), material={"EX": E, "PRXY": NU, "DENS": RHO}, ) # Full clamp at x = 0. pts = np.asarray(m.grid.points) clamped = np.where(pts[:, 0] < 1e-9)[0] m.fix(nodes=(clamped + 1).tolist(), dof="ALL") # Total mass for the running-sum normalization. total_mass = RHO * L * WIDTH * HEIGHT print("Cantilever beam — modal participation factors") print(f" L = {L} m, cross = {WIDTH} × {HEIGHT} m, ρ = {RHO} kg/m³") print(f" total mass M_total = ρ·L·A = {total_mass:.3f} kg") # %% # Modal solve + grab the consistent mass matrix # --------------------------------------------- N_MODES = 12 res = m.solve_modal(n_modes=N_MODES) freqs = np.asarray(res.frequency, dtype=np.float64) shapes_dof = np.asarray(res.mode_shapes, dtype=np.float64) # ``mode_shapes`` is shaped ``(n_dof, n_modes)``. Each column is a # mass-normalised eigenvector φ_i with φ_i^T · M · φ_i = 1. M = m.mass_matrix() n_dof = M.shape[0] # %% # Build the X / Y / Z influence vectors # ------------------------------------- # # r^(d) is 1 on every DOF aligned with direction d, 0 elsewhere. # DOFs are interleaved (UX, UY, UZ) per node — read them off the # DOF map. dof_map = m.dof_map() r_x = np.zeros(n_dof) r_y = np.zeros(n_dof) r_z = np.zeros(n_dof) for row, (_node, dof_idx) in enumerate(dof_map.tolist()): if dof_idx == 0: # UX r_x[row] = 1.0 elif dof_idx == 1: # UY r_y[row] = 1.0 elif dof_idx == 2: # UZ r_z[row] = 1.0 # %% # Compute Γ and M_eff for every mode × every direction # ---------------------------------------------------- directions = (("X", r_x), ("Y", r_y), ("Z", r_z)) gamma = np.zeros((N_MODES, 3), dtype=np.float64) m_eff = np.zeros((N_MODES, 3), dtype=np.float64) for i in range(N_MODES): phi = shapes_dof[:, i] # Mass-normalised: φ^T M φ = 1, so Γ = φ^T M r and M_eff = Γ² M_phi = M @ phi modal_mass_i = float(phi @ M_phi) # ≈ 1 for normalised modes for d_idx, (_, r) in enumerate(directions): L_i = float(phi @ M @ r) # = modal-load coefficient gamma_i = L_i / modal_mass_i gamma[i, d_idx] = gamma_i m_eff[i, d_idx] = gamma_i**2 * modal_mass_i cum_eff = np.cumsum(m_eff, axis=0) cum_pct = cum_eff / total_mass * 100.0 print() print( f"{'mode':>4} {'f [Hz]':>10} {'Γ_X':>10} {'Γ_Y':>10} {'Γ_Z':>10} " f"{'M_eff,X [kg]':>12} {'M_eff,Y':>9} {'M_eff,Z':>9} " f"{'cum_X %':>9} {'cum_Y %':>9} {'cum_Z %':>9}" ) print(" " + "-" * 130) for i in range(N_MODES): print( f"{i + 1:>4} {freqs[i]:>10.3f} " f"{gamma[i, 0]:>+10.4f} {gamma[i, 1]:>+10.4f} {gamma[i, 2]:>+10.4f} " f"{m_eff[i, 0]:>11.3f} {m_eff[i, 1]:>8.3f} {m_eff[i, 2]:>8.3f} " f"{cum_pct[i, 0]:>8.2f}% {cum_pct[i, 1]:>8.2f}% {cum_pct[i, 2]:>8.2f}%" ) print() print( f" total M_eff after {N_MODES} modes: " f"X = {cum_pct[-1, 0]:.1f} %, " f"Y = {cum_pct[-1, 1]:.1f} %, " f"Z = {cum_pct[-1, 2]:.1f} % (of M_total = {total_mass:.3f} kg)" ) # %% # Verify the partition-of-unity property # -------------------------------------- # # Sanity check: a sufficient mode count must capture nearly the # full mass in every direction. Twelve modes on a slender beam # usually clear 90 % in the bending directions. assert cum_pct[-1, 1] > 80.0, ( f"Y-direction effective mass {cum_pct[-1, 1]:.1f} % below 80 % — " "increase n_modes or check the mass matrix." ) assert cum_pct[-1, 2] > 80.0, f"Z-direction effective mass {cum_pct[-1, 2]:.1f} % below 80 %." # %% # Render the cumulative effective-mass curve # ------------------------------------------ import matplotlib.pyplot as plt # noqa: E402 fig, ax = plt.subplots(1, 1, figsize=(7.0, 4.0)) mode_idx = np.arange(1, N_MODES + 1) ax.plot(mode_idx, cum_pct[:, 0], "o-", color="#1f77b4", lw=2, label="cum. M_eff,X") ax.plot(mode_idx, cum_pct[:, 1], "s-", color="#d62728", lw=2, label="cum. M_eff,Y") ax.plot(mode_idx, cum_pct[:, 2], "^-", color="#2ca02c", lw=2, label="cum. M_eff,Z") ax.axhline(90.0, color="black", ls="--", lw=1.0, label="90 % design target") ax.set_xlabel("number of modes included") ax.set_ylabel("cumulative effective mass [% of M_total]") ax.set_title("Cantilever — modal effective mass per direction", fontsize=11) ax.set_xticks(mode_idx) ax.set_ylim(0.0, 105.0) ax.legend(loc="lower right", fontsize=9) ax.grid(True, ls=":", alpha=0.5) fig.tight_layout() fig.show() # %% # Take-aways # ---------- print() print("Take-aways:") print( " • Modal participation factor Γ_i = (φ_i^T M r) / (φ_i^T M φ_i) " "answers 'how strongly does mode i respond to base excitation in direction d?'." ) print( " • Effective modal mass M_eff,i = Γ_i² · (φ_i^T M φ_i) sums to the " "total mass of the structure when summed over ALL modes." ) print( " • The first transverse-bending mode dominates one direction (~60-70 %) " "and barely contributes to the others — a single number tells you which " "modes can be skipped per excitation axis." ) print( " • Design codes (ASCE 7, Eurocode 8) usually require ≥ 90 % cumulative " "effective mass per direction — read directly off the table above to " "decide n_modes for a response-spectrum analysis." )