FEA Formulations & Mass Matrices

5.1 Spatial Discretization and Variational Foundations

In structural dynamics, noise, vibration, and harshness (NVH) engineering, we are routinely confronted with systems whose geometric complexity, material anisotropy, and boundary conditions render analytical closed-form solutions of their partial differential equations (PDEs) impossible. Chapters 2 and 3 established the methods to solve discrete Multi-Degree-of-Freedom (MDOF) systems, and Chapter 4 formulated the exact solutions for one-dimensional continuous systems like strings, bars, and Euler-Bernoulli beams. However, analyzing a three-dimensional automotive body-in-white (BIW), an engine block, or a passenger cabin cavity requires a robust, general, and mathematically rigorous framework to discretize continuous systems into solvable discrete models. This is the role of the Finite Element Method (FEM).

The Finite Element Method is not merely a tool for dividing a structure into a mesh of small blocks; it is a variational approximation technique that translates the continuous infinite-dimensional boundary value problem into a finite-dimensional algebraic problem. To understand the dynamics of a discretized structure, we must begin with the fundamental variational principles of mechanics, specifically Hamilton's principle.

5.1.1 Continuous System State Variables and Kinetic Energy

Consider a three-dimensional elastic body occupying a volume $V$ bounded by a surface $S$. Let the vector displacement field of a point within the body be denoted by:

where $\mathbf{x} = [x, y, z]^T$ represents the spatial coordinates in the reference configuration, and $t$ represents time.

The kinetic energy $T$ of the continuous body is the integral of the kinetic energy of its infinitesimal mass elements $dm = \rho(\mathbf{x}) dV$:

where $\rho(\mathbf{x})$ is the mass density of the material.

5.1.2 Potential Energy, Stress-Strain Relations, and Virtual Work

The potential energy of the body is composed of the strain energy $U$ stored in the deformed elastic medium. Assuming a linear elastic, isotropic material, the strain energy is:

where $\boldsymbol{\epsilon}$ is the vector of engineering strains and $\boldsymbol{\sigma}$ is the vector of corresponding stresses:

The strain-displacement relationship is given by the linear kinematic relations:

which can be compactly written using a differential operator matrix $\mathbf{L}$:

The stress-strain constitutive relation is governed by Hooke's law:

where $\mathbf{D}$ is the symmetric elasticity (material stiffness) matrix containing the elastic moduli. For an isotropic material, $\mathbf{D}$ is determined solely by Young's modulus $E$ and Poisson's ratio $\nu$.

The virtual work $\delta W$ done by external forces consists of body forces $\mathbf{f}_b$ (such as gravity or electromagnetic fields) acting over the volume $V$ and surface tractions $\mathbf{f}_s$ (such as pressure or friction) acting on the boundary surface portion $S_{\sigma}$:

where $\delta \mathbf{u}$ represents a kinematically admissible virtual displacement field.

5.1.3 Hamilton's Principle and the Equations of Motion

Hamilton's principle states that the actual trajectory of an elastic body between two fixed times $t_1$ and $t_2$ extremizes the action integral of the Lagrangian $L = T - U$. Mathematically:

subject to the condition that the virtual displacements vanish at the temporal boundaries:

We can derive the equations of motion by taking the variation of each term. First, let us evaluate the variation of the kinetic energy:

Using integration by parts over time for the action integral:

Since the virtual displacement $\delta \mathbf{u}$ is zero at $t_1$ and $t_2$, the boundary term vanishes:

Next, we evaluate the variation of the strain energy:

Substituting these variations into Hamilton's principle yields:

Since this integral must hold for any arbitrary, kinematically admissible virtual displacement field $\delta \mathbf{u}$, the terms inside the brackets must satisfy the weak variational form of the dynamic equilibrium equation:

This is the principle of virtual displacements for dynamic systems. It serves as the foundation for the finite element formulation, as it does not require the displacement field to possess the higher-order derivatives needed for the structural Navier-Cauchy equations.

5.1.4 The Ritz Method: Trial Functions and Minimization

To appreciate how spatial discretization works, it is instructive to examine the classical Ritz Method, which predates the modern finite element method. In the Ritz method, we approximate the continuous displacement field $u(x)$ using a linear combination of globally defined, linearly independent trial functions $\psi_i(x)$ multiplied by unknown generalized coordinates $q_i$:

These trial functions must satisfy the essential (geometric) boundary conditions of the system. Let us consider a static one-dimensional elastic bar of length $L$, fixed at $x=0$, under a uniform distributed axial load $q(x) = q_0$. The total potential energy functional of the bar is:

The essential boundary condition is $u(0) = 0$. We select a two-term polynomial approximation:

which automatically satisfies the boundary condition at $x=0$. Here, the trial functions are $\psi_1(x) = x$ and $\psi_2(x) = x^2$, and the generalized coordinates are $c_1$ and $c_2$. The derivative of the displacement is:

Substituting this into the energy functional:

To minimize the potential energy with respect to the coordinates $c_1$ and $c_2$, we take the partial derivatives and set them to zero:

We solve this system of equations for $c_1$ and $c_2$: Subtracting the first equation from the second:

Substituting $c_2$ back into the first equation:

Thus, the Ritz approximation for the displacement field is:

In this specific case, the two-term Ritz approximation matches the exact closed-form solution derived from the differential equation. However, if the load $q(x)$ were non-uniform or the cross-sectional area $A(x)$ were variable, the Ritz approximation would yield an approximate solution whose accuracy depends on the number and choice of trial functions.

The Finite Element Method is a generalization of the Ritz method. Instead of using global trial functions that span the entire structure (which are difficult to construct for complex geometries), FEM divides the domain into small elements and defines simple, localized shape functions within each element.

---

5.2 The Concept of Mass in Dynamic FEA

To solve the weak variational form using a computer, we must discretize the continuous domain. We divide the volume $V$ into a set of $M$ non-overlapping subdomains, called finite elements:

Within each element $V_e$, we approximate the continuous displacement field $\mathbf{u}(\mathbf{x}, t)$ by interpolating from the displacements at discrete points on the element boundary, known as nodes.

5.2.1 Element-Level Shape Functions and Interpolation

Let the vector of nodal displacements (degrees of freedom) for a single element $e$ be:

where $n$ is the number of degrees of freedom of the element. The continuous displacement field within the element is interpolated using the element shape function matrix $\mathbf{N}(\mathbf{x})$:

The shape functions must satisfy several critical mathematical properties:

$N_i(\mathbf{x}_j)$ must satisfy:

This ensures that the displacement field at node $j$ matches the nodal degree of freedom $d_j$.

element must equal one:

This is a requirement for representing rigid-body translation. If all nodal displacements are equal (e.g., $d_i = c$), the displacement at any point within the element must be $u(\mathbf{x}) = c \sum N_i(\mathbf{x}) = c$.

ensure that the element stiffness and mass matrices are non-singular.

across element boundaries. For elements involving only first-order derivatives in the strain energy (e.g., axial bars, 2D planes, and 3D solids), the displacements themselves must be continuous across boundaries. This is known as **$C^0$ continuity**. For elements involving second-order derivatives in the strain energy (e.g., Euler-Bernoulli beams and thin plates), both the displacements and their first derivatives (slopes) must be continuous across boundaries. This is known as **$C^1$ continuity**. If compatibility is violated, fictitious energy sinks or sources are created at the boundaries, and the solution will fail to converge as the mesh is refined.

5.2.2 Mathematical Derivation of consistent and Lumped Matrices

By substituting the shape function approximation into the weak variational form, we can transform the integral equation into a set of ordinary differential equations. Let us express the virtual displacement field within the element using the same shape functions (the Ritz-Galerkin approach):

Substituting these expressions into the left-hand side of the weak form for a single element $e$:

1. Kinetic Energy Term

We define the element Consistent Mass Matrix $\mathbf{m}^e_c$ as:

The term "consistent" is used because this matrix is derived using the same shape functions $\mathbf{N}$ that are used to formulate the stiffness matrix. It represents the mathematically exact distribution of kinetic energy under the assumed interpolation field.

2. Strain Energy Term

We define the strain-displacement matrix $\mathbf{B} = \mathbf{L} \mathbf{N}$. Substituting this in:

We define the element Stiffness Matrix $\mathbf{k}^e$ as:

3. External Force Term

We define the element Consistent Nodal Force Vector $\mathbf{f}^e$ as:

Substituting these definitions back into the variational statement, we obtain the element-level equations of motion:

These element-level matrices are then assembled into the global equations of motion using the direct stiffness method:

where $\mathbf{M}_c$, $\mathbf{K}$, and $\mathbf{F}$ represent the global consistent mass matrix, global stiffness matrix, and global force vector, respectively, and $\mathbf{U}$ is the vector of global nodal coordinates.

---

5.3 Lumped Mass Matrix Formulations

While the consistent mass matrix is mathematically rigorous, it is not always the most computationally efficient choice. An alternative approach is to construct a Lumped Mass Matrix (LMM), denoted by $\mathbf{m}^e_l$. The physical concept behind lumping is to concentrate the distributed mass of the element at its nodal points, creating a set of discrete point masses. As a result, the lumped mass matrix is strictly diagonal:

This diagonal structure has profound implications for computational efficiency, particularly in transient dynamics.

5.3.1 Row-Sum Technique

The simplest method to construct a lumped mass matrix is the Row-Sum Technique. In this method, the diagonal term for each degree of freedom is obtained by summing all the terms in the corresponding row of the consistent mass matrix:

The physical justification is that it preserves the total mass of the element. For example, for a linear bar element with consistent mass matrix:

Applying the row-sum technique:

which yields the lumped mass matrix:

This represents exactly half of the bar's mass concentrated at each end.

However, the row-sum technique is not universally applicable. For higher-order elements or elements with rotational degrees of freedom (like beams and plates), row-summing can result in physical anomalies. For instance, summing the rows of a beam element's consistent mass matrix can yield negative or zero values for the rotational degrees of freedom. This would mean that the element has zero or negative rotational inertia, which violates physical laws and leads to a singular mass matrix.

5.3.2 Hinton-Rock-Zienkiewicz (HRZ) Lumping Scheme

To overcome the limitations of the row-sum technique, Hinton, Rock, and Zienkiewicz proposed the HRZ Lumping Scheme in 1976. This scheme ensures that the lumped mass matrix is always positive-definite, preserves the total mass of the element, and is applicable to any element type, including those with rotational degrees of freedom.

The HRZ lumping algorithm proceeds as follows:

where $\mathcal{T}$ represents the set of indices corresponding to translational DOFs.

this sum:

This scaling is applied to both translational and rotational diagonal terms.

The HRZ scheme preserves the relative inertia distribution predicted by the consistent mass matrix while guaranteeing a diagonal, positive-definite matrix. Because the diagonal terms of the consistent mass matrix are always positive, the resulting lumped terms are also guaranteed to be positive.

5.3.3 Computational and Physical Comparison of CMM and LMM

The choice between a consistent mass matrix and a lumped mass matrix involves a trade-off between mathematical accuracy, frequency range of interest, and the numerical solver used for the simulation.

1. Eigenvalue Boundedness and Frequency Convergence

The consistent mass matrix is derived from the Ritz-Galerkin approach, which imposes an implicit constraint on the displacement field (limiting it to the shapes represented by the interpolation functions $\mathbf{N}$). This variational constraint artificially stiffens the system. Consequently, the natural frequencies calculated using CMM are mathematically proven to be upper bounds (or equal) to the exact continuous eigenvalues (assuming the mesh is sufficiently refined to capture the modes).

Conversely, the lumped mass matrix represents a physical system with concentrated masses and massless structural members. This arrangement does not satisfy the variational stiffness constraints in the same manner. As a result, the lumped mass model is typically more flexible (less stiff) than the consistent mass model, and the natural frequencies calculated using LMM are generally lower bounds to the exact continuous values.

This bounding behavior is highly useful. By running a modal analysis with CMM and then with LMM, an engineer can establish a mathematical bracket that contains the exact natural frequencies of the continuous system.

2. Wave Propagation and Grid Dispersion

In high-frequency dynamics and structural acoustics, we analyze the propagation of elastic waves through a mesh. The choice of mass matrix determines the wave speed and numerical dispersion properties.

Due to the artificial stiffening effect of CMM, high-frequency waves propagate through a CMM mesh at speeds that are faster than the physical speed of sound in the material. This is known as a positive phase velocity error (leading phase error). Conversely, the LMM mesh underestimates the high-frequency wave speeds, causing waves to propagate slower than the physical speed. This is a **negative phase velocity error** (lagging phase error).

We can derive the numerical dispersion relations for axial wave propagation in a uniform 1D bar discretized with elements of length $h$. The exact wave relation is $\omega(k) = c k$, where $c = \sqrt{E/\rho}$ is the physical sound speed.

node $j$ is:

Assuming a harmonic wave solution $u_j(t) = U e^{i(k x_j - \omega t)}$ where $x_j = j h$:

Dividing by $U e^{i(jkh - \omega t)}$:

Using the trigonometric identity $1 - \cos(kh) = 2 \sin^2(kh/2)$:

We perform a Taylor series expansion for small $kh$ (long wavelengths relative to mesh size):

This shows that the numerical wave speed $c_{\text{num}} = \omega_{\text{LMM}}/k$ is:

The term $-\frac{(kh)^2}{24}$ is a negative phase velocity error (lagging phase error).

at node $j$ is:

Substituting the harmonic wave solution:

Expanding this in Taylor series for small $kh$:

This yields a positive phase velocity error:

The corresponding numerical wave speed is:

By setting $\theta = 0.5$ (or $\theta = 1/2$), the second-order error term cancels out:

This cancels the leading dispersion error, providing $O(h^4)$ accuracy and making the blended mass matrix highly effective for structural acoustics simulations.

3. Explicit Solver Time Step and Numerical Stability

For transient dynamics solved with explicit schemes (such as the Central Difference Method), the simulation must satisfy the Courant-Friedrichs-Lewy (CFL) stability criterion. The time step $\Delta t$ must be smaller than a critical value:

where $\omega_{\text{max}}$ is the highest natural frequency of the mesh.

For a 1D linear bar element, the highest element-level natural frequency is:

Since CMM makes the mesh artificially stiffer, it increases the highest natural frequency by a factor of $\sqrt{3}$. Consequently, using CMM reduces the critical time step by a factor of $\sqrt{3}$:

This means that using CMM requires roughly 73% more time steps than LMM to solve the same duration of transient response, in addition to the overhead of solving a non-diagonal system of equations at each step. This makes LMM the clear choice for explicit dynamic solvers in commercial codes.

5.4 One-Dimensional Bar and Beam Elements

To illustrate the mathematical steps of finite element discretization, we formulate the stiffness and mass matrices for two of the most widely used elements in structural dynamics: the one-dimensional bar (axial vibration) and the Euler-Bernoulli beam (flexural vibration). We also examine the Timoshenko beam element to see how shear deformation affects the formulation.

5.4.1 Linear Bar Element (Axial Vibration)

Consider a uniform thin bar of length $L$, cross-sectional area $A$, Young's modulus $E$, and mass density $\rho$, aligned with the local coordinate $x$ from $x=0$ to $x=L$. The bar undergoes only axial displacement $u(x, t)$.

1. Shape Functions

The bar is discretized using a two-node element. At Node 1 ($x=0$), the displacement is $u_1(t)$. At Node 2 ($x=L$), the displacement is $u_2(t)$. The vector of nodal degrees of freedom is:

We interpolate the displacement $u(x, t)$ linearly between the nodes:

The linear shape functions are:

These functions satisfy the Kronecker delta property ($N_1(0)=1$, $N_1(L)=0$, etc.) and sum to unity at any point, which ensures that they can represent rigid-body translation.

2. Element Stiffness Matrix Derivation

The strain in the bar is:

The strain-displacement matrix $\mathbf{B}(x)$ contains the spatial derivatives of the shape functions:

Substituting this into the stiffness matrix definition:

This is the standard axial stiffness matrix.

3. Consistent Mass Matrix Derivation

Using the shape functions, we compute the consistent mass matrix:

Evaluating each integral term-by-term:

\rho A \left[ x - \frac{x^2}{L} + \frac{x^3}{3L^2} \right]_0^L = \frac{\rho A L}{3}$

\rho A \left[ \frac{x^2}{2L} - \frac{x^3}{3L^2} \right]_0^L = \frac{\rho A L}{6}$

\right]_0^L = \frac{\rho A L}{3}$ Assembling these terms yields:

5.4.2 Euler-Bernoulli Beam Element (Flexural Vibration)

The Euler-Bernoulli beam theory governs the flexural vibration of slender beams. The theory assumes that plane sections of the beam remain plane, rigid, and perpendicular to the longitudinal neutral axis after deformation. This implies that shear deformation is neglected, which is a highly accurate assumption for beams whose length-to-height ratio is large ($L/h > 10$).

1. Beam Kinematics and Degrees of Freedom

Let the transverse displacement of the beam's neutral axis be $w(x, t)$. Under the Euler-Bernoulli assumptions, the rotation of the cross-section is equal to the slope of the neutral axis:

The strain energy stored in a beam element depends on the curvature (second spatial derivative of $w$):

Because the variational formulation involves the second derivative of displacement, the shape functions must maintain continuity of both displacement and slope at the nodes to ensure a finite strain energy. This requires $C^1$ continuity.

A 2-node beam element must therefore have at least two degrees of freedom at each node: transverse displacement $w$ and slope $\theta$. The vector of element degrees of freedom is:

2. Cubic Hermite Shape Functions

To interpolate four nodal degrees of freedom, we use a cubic polynomial displacement field:

By solving for the coefficients $a_i$ in terms of the nodal coordinates $w_1, \theta_1, w_2, \theta_2$, we obtain the four Cubic Hermite Shape Functions:

displacement $w_1$)

slope $\theta_1$)

displacement $w_2$)

$\theta_2$)

These shape functions are plotted below inside the element coordinate system:

3. Consistent Stiffness Matrix

The curvature of the beam is:

where $\mathbf{B}_b(x)$ contains the second spatial derivatives of the Hermite shape functions:

Evaluating the derivatives:

Substituting these into the stiffness matrix formulation $\mathbf{k}^e = \int_0^L \mathbf{B}_b^T EI \mathbf{B}_b \, dx$ and integrating yields the consistent stiffness matrix:

4. Consistent Translational Mass Matrix

The consistent translational mass matrix is calculated using the Hermite shape functions:

Performing this integration (which is derived step-by-step in the Derivations section) yields:

5. Consistent Rotatory Inertia Mass Matrix

In Euler-Bernoulli beam theory, the rotational kinetic energy of the beam elements is often neglected because it is small compared to the translational kinetic energy. However, for higher-frequency modes, the slope of the beam changes rapidly, and the rotatory inertia of the cross-sections begins to contribute.

The kinetic energy due to cross-sectional rotation is:

where $I$ is the second moment of area. This leads to the consistent rotatory inertia mass matrix:

Performing this integration yields:

To include both translation and rotation, we add the two matrices:

5.4.3 Timoshenko Beam Element and Shear Locking

For beams that are thick or short ($L/h < 10$), the Euler-Bernoulli assumptions are no longer valid because shear deformation becomes significant. Timoshenko beam theory incorporates shear deformation by treating the transverse displacement $w(x, t)$ and the cross-sectional rotation $\phi(x, t)$ as independent variables.

1. Timoshenko Kinematics

The shear strain $\gamma_s$ is the difference between the slope of the neutral axis and the rotation of the cross-section:

The strain energy now includes a bending term and a shear term:

where $G$ is the shear modulus and $\kappa$ is the shear correction factor (which accounts for the non-uniform distribution of shear stress across the cross-section, typically $5/6$ for rectangular sections).

2. Shear Locking and the Shear Parameter $\Phi_s$

To model Timoshenko beams without encountering shear locking, we define a dimensionless shear parameter $\Phi_s$, which represents the ratio of bending stiffness to shear stiffness:

For very slender beams ($L/h \to \infty$), the shear stiffness becomes infinite, which forces $\Phi_s \to 0$. In this limit, the Timoshenko element must converge to the Euler-Bernoulli element.

By using interdependent interpolation fields for $w$ and $\phi$ to satisfy the shear strain relations, we can derive the exact, locking-free element matrices.

3. Timoshenko Element Stiffness Matrix ($\mathbf{k}^e_{\text{Timo}}$)

The $4 \times 4$ element stiffness matrix for a 2-node Timoshenko beam element in terms of the shear parameter $\Phi_s$ is:

Note that as $\Phi_s \to 0$, this matrix converges exactly to the Euler-Bernoulli stiffness matrix.

4. Timoshenko Element Consistent Mass Matrix ($\mathbf{m}^e_{c, \text{Timo}}$)

The consistent translational mass matrix is also modified by the shear parameter $\Phi_s$:

where the individual components are:

This formulation allows for high-fidelity modeling of short brackets and thick support plates common in automotive powertrain mounts.

---

5.5 Eigensolvers and Modal Analysis in FEA

Discretization yields the global equations of motion:

To analyze the natural frequencies and mode shapes of the structure, we set $\mathbf{F}(t) = \mathbf{0}$ and assume a harmonic free response:

Substituting this into the equations of motion yields the **Generalized Eigenvalue Problem**:

where $\omega$ represents the natural frequencies (eigenvalues) and $\boldsymbol{\Phi}$ represents the mode shapes (eigenvectors).

5.5.1 Mathematical Proofs of Orthogonality

The mode shapes possess a critical mathematical property: they are orthogonal with respect to the mass and stiffness matrices. This property allows us to decouple the multi-degree-of-freedom equations of motion into a set of independent single-degree-of-freedom equations.

Let $\omega_i^2$ and $\boldsymbol{\phi}_i$ be the $i$-th eigenvalue and eigenvector, and let $\omega_j^2$ and $\boldsymbol{\phi}_j$ be the $j$-th eigenvalue and eigenvector. They satisfy:

Premultiplying Eq. 1 by $\boldsymbol{\phi}_j^T$:

Premultiplying Eq. 2 by $\boldsymbol{\phi}_i^T$ and taking the transpose (noting that $\mathbf{K}$ and $\mathbf{M}$ are symmetric):

Subtracting Eq. 4 from Eq. 3:

If the eigenvalues are distinct ($\omega_i^2 \neq \omega_j^2$), we must have:

Substituting this back into Eq. 3 yields:

This completes the proof of orthogonality.

5.5.2 Mass-Normalization

The amplitude of eigenvectors is arbitrary. In finite element analysis, it is standard practice to normalize the eigenvectors such that the modal mass is equal to unity:

Under this mass-normalization scheme, the orthogonality relations simplify to:

where $\boldsymbol{\Phi}$ is the modal matrix containing the eigenvectors as columns, $\mathbf{I}$ is the identity matrix, and $\boldsymbol{\Omega}^2 = \text{diag}(\omega_1^2, \omega_2^2, \dots, \omega_N^2)$.

5.5.3 Eigensolution Algorithms: Subspace Iteration and Lanczos Methods

For large finite element models containing millions of degrees of freedom, calculating the eigenvalues by solving the characteristic equation $\det(\mathbf{K} - \omega^2 \mathbf{M}) = 0$ is impossible. Specialized numerical algorithms are required to calculate a subset of the lowest eigenvalues:

the stiffness and mass matrices. It is highly robust but computationally expensive, limiting its use to small systems.

projects the large eigenvalue problem onto a smaller subspace spanned by a set of iteration vectors, solving the reduced problem using the Jacobi method.

orthogonal basis of the Krylov subspace generated by the operator $\mathbf{K}^{-1}\mathbf{M}$. This tridiagonalizes the system matrices, allowing for the rapid extraction of the lowest eigenvalues.

---

5.6 Modal Participation and Effective Mass

In vehicle NVH and aerospace design, structures are frequently subjected to base excitation (e.g., a car driving over a rough road or a rocket launch). To evaluate how strongly the structure responds to this base motion, we use the concepts of modal participation factors and modal effective mass.

5.6.1 Base Excitation Formulation

Let the global displacement vector $\mathbf{U}(t)$ be partitioned into the displacement relative to the base, $\mathbf{U}_{\text{rel}}(t)$, and the rigid-body displacement due to base motion:

where $u_g(t)$ is the base displacement and $\mathbf{r}$ is the influence vector (or excitation vector). The influence vector represents the static displacement of the degrees of freedom due to a unit displacement of the base. For example, for base excitation in the $x$-direction, $\mathbf{r}$ has values of $1.0$ for all $x$-translational DOFs and $0.0$ for all other DOFs.

Substituting this partition into the undamped equations of motion:

Since a rigid-body displacement does not store strain energy, $\mathbf{K} \mathbf{r} = \mathbf{0}$. The equations simplify to:

5.6.2 Modal Transformation

We transform the equations into modal coordinates using the modal matrix $\boldsymbol{\Phi}$:

Substituting this and premultiplying by $\boldsymbol{\Phi}^T$:

Using the mass-normalization relations:

where the vector of Modal Participation Factors $\boldsymbol{\Gamma}$ is defined as:

For the $i$-th mode, the participation factor is:

This factor measures how strongly the $i$-th mode is excited by the base acceleration.

5.6.3 Modal Effective Mass

The Modal Effective Mass of the $i$-th mode, $m_{\text{eff}, i}$, is defined as:

The physical meaning of modal effective mass is the portion of the system's total mass that participates in the dynamic response of that mode.

A fundamental theorem in structural dynamics states that the sum of the modal effective masses over all modes equals the total mass of the system in the direction of the base excitation:

(The proof of this identity is derived in the Derivations section).

5.6.4 Modal Truncation and the 90% Rule

For large models with millions of degrees of freedom, calculating all modes is computationally impossible. We must truncate the modal expansion, retaining only the lowest $P$ modes.

To determine if we have retained enough modes to accurately capture the dynamic response under base excitation, we monitor the cumulative modal effective mass:

According to industry standards (such as NASA-STD-5002 and automotive NVH specifications), the modal expansion must include a sufficient number of modes to capture at least 90% of the total structural mass in the direction of excitation. If the cumulative mass is less than 90%, it indicates that key modes have been omitted, which would lead to significant errors in the predicted dynamic stresses and cabin accelerations.

---

5.7 Dynamic Reduction Schemes

Even with high-performance computing, conducting transient or harmonic analyses on detailed automotive body-in-white models (which contain millions of DOFs) is computationally intensive. To make these simulations practical, we use dynamic reduction schemes to condense the large system of equations into a much smaller, equivalent system.

5.7.1 Guyan Reduction (Static Condensation)

Guyan reduction partitions the structural degrees of freedom into two sets:

nodes, connection points, or nodes where forces are applied).

The partitioned undamped equations of motion are:

Assuming that no external forces act on the slave DOFs.

In Guyan reduction, we neglect the inertia terms in the second row of the equations:

Solving for the slave displacements:

We define the static transformation matrix $\mathbf{T}_G$:

Using this static coordinate transformation, the slave degrees of freedom are eliminated. The reduced system matrices are:

Guyan reduction is mathematically exact for static problems. However, because it completely neglects the slave inertia terms $\mathbf{M}_{ss} \mathbf{\ddot{U}}_s$ during the construction of the transformation matrix $\mathbf{T}_G$, it is only accurate for low frequencies where the structural response is dominated by stiffness. At higher frequencies, the neglected inertia terms lead to significant errors in the natural frequencies, typically overestimating the system stiffness.

The relative error in the calculated eigenvalues $\lambda_i = \omega_i^2$ can be estimated using:

where $\omega_{s, 1}$ is the lowest natural frequency of the slave structure when the master DOFs are clamped ($\mathbf{U}_a = \mathbf{0}$). This demonstrates that Guyan reduction is highly accurate when the excitation frequency is well below the lowest resonant frequency of the clamped slave substructure.

5.7.2 Craig-Bampton Method (Component Mode Synthesis)

The Craig-Bampton method is the industry standard for dynamic reduction. It resolves the limitations of Guyan reduction by incorporating the internal dynamic modes of the sub-assembly (fixed-interface normal modes) to supplement the static constraint modes.

The degrees of freedom are partitioned into:

rest of the structure.

The Craig-Bampton transformation matrix $\mathbf{T}_{CB}$ is composed of two sets of modes:

the internal DOFs due to a unit displacement at each boundary DOF while all other boundary DOFs are fixed:

internal structure obtained by holding the boundary DOFs fixed ($\mathbf{U}_b = \mathbf{0}$):

We retain only the lowest $P$ modes.

The displacement coordinates are transformed as:

where $\mathbf{q}_n$ is the vector of modal coordinates for the fixed-interface modes.

By projecting the global matrices using $\mathbf{T}_{CB}$, we obtain the Craig-Bampton reduced matrices:

where $\boldsymbol{\Omega}_n^2$ is the diagonal matrix of the fixed-interface natural frequencies.

This formulation yields a decoupled stiffness matrix while retaining the dynamic contributions of the internal modes through the mass coupling terms $\mathbf{M}_{bn}^*$. The Craig-Bampton method allows for the assembly of dynamic sub-assemblies (components) into a complete vehicle model with high accuracy across a wide frequency range.

---

5.8 Damping Formulations in FEA

Damping represents the dissipation of mechanical energy in a vibrating structure. In finite element models, incorporating damping is essential to predict the amplitude of vibrations at resonance.

5.8.1 Rayleigh (Proportional) Damping

The most common damping model in linear structural dynamics is Rayleigh damping. In this model, the damping matrix $\mathbf{C}$ is assumed to be a linear combination of the mass and stiffness matrices:

where $\alpha$ and $\beta$ are the Rayleigh coefficients.

1. Modal Damping Ratios

By applying the modal transformation, the damping matrix is also diagonalized:

For the $i$-th mode, the modal damping ratio $\zeta_i$ is:

2. Physical Behavior of Rayleigh Coefficients

dominates at low frequencies. It behaves like a drag force acting against the absolute velocity of the nodes, which can be interpreted as external fluid damping.

dominates at high frequencies. It is proportional to the strain rate, representing internal material hysteresis.

3. Calculating Rayleigh Coefficients

If we specify target damping ratios $\zeta_1$ and $\zeta_2$ at two critical frequencies $\omega_1$ and $\omega_2$, we can solve for $\alpha$ and $\beta$ by setting up a system of two equations:

Solving this system yields:

Rayleigh damping is highly popular because it allows us to solve the damped equations of motion in modal coordinates without introducing coupling between the modes. However, it has a significant limitation: modes at frequencies outside the interval $[\omega_1, \omega_2]$ are heavily over-damped, which can lead to artificial attenuation of high-frequency responses.

5.9 Advanced Numerical Eigensolvers in Finite Element Analysis

For structural systems discretized with millions of degrees of freedom, solving the characteristic polynomial $\det(\mathbf{K} - \omega^2 \mathbf{M}) = 0$ is computationally impossible. Instead, modern commercial FEA codes utilize advanced numerical methods to extract the lowest natural frequencies and mode shapes.

5.9.1 The Subspace Iteration Method

The Subspace Iteration Method is a classical, highly robust algorithm designed to compute the lowest $p$ eigenvalues and eigenvectors of the generalized eigenvalue problem:

The method is a multi-vector generalization of the inverse power iteration. Rather than iterating a single vector, it operates on a subspace of dimension $q$ (where $q > p$, typically $q = \min(2p, p + 8)$ to ensure fast convergence of the higher modes of interest).

1. The Algorithm Steps

Let $\mathbf{X}_0$ be an starting iteration matrix of size $N \times q$, where $N$ is the number of degrees of freedom in the FEA model. The columns of $\mathbf{X}_0$ represent starting vectors. For each iteration step $k = 1, 2, \dots$:

$\mathbf{\bar{X}}_k$:

This step represents an inverse power projection. Since $\mathbf{K}$ is symmetric and positive-definite, it is factorized once using Cholesky decomposition ($\mathbf{K} = \mathbf{L}\mathbf{L}^T$) or LDL decomposition, and the solve step is performed via forward and backward substitutions, which is highly efficient.

the columns of $\mathbf{\bar{X}}_k$:

The projected matrices $\mathbf{K}_k^*$ and $\mathbf{M}_k^*$ are symmetric and of small dimension $q \times q$.

Because $q$ is small (typically less than 100), this reduced problem is solved quickly using direct methods like the Jacobi method or QZ algorithm. The diagonal matrix $\boldsymbol{\Omega}_k^2$ contains the eigenvalues at iteration step $k$.

to the original space:

This step rotates the iteration vectors to align them with the principal axes of the projected subspace, preventing them from collapsing onto the same dominant eigenvector.

where $\lambda_i^{(k)}$ is the $i$-th diagonal element of $\boldsymbol{\Omega}_k^2$.

The subspace iteration method is highly stable and guaranteed to converge to the lowest $p$ eigenvalues, even in the presence of multiple or closely spaced eigenvalues. Its convergence rate is proportional to:

where $\lambda_i$ is the $i$-th eigenvalue and $\lambda_{q+1}$ is the first eigenvalue outside the iteration subspace. This explains why we choose the subspace dimension $q$ to be larger than the number of requested eigenvalues $p$, as it accelerates the convergence of the $p$-th mode.

5.9.2 The Lanczos Method

The Lanczos Method is the industry standard for large-scale dynamic FEA models due to its superior computational speed compared to subspace iteration. The method constructs an orthonormal basis of the Krylov subspace:

where $\mathbf{A} = \mathbf{K}^{-1}\mathbf{M}$ and $\mathbf{v}_1$ is a starting vector.

1. Lanczos Recursion

The Lanczos algorithm generates a sequence of vectors $\mathbf{v}_j$ (called Lanczos vectors) that are $\mathbf{M}$-orthonormal:

The recursion relation is defined as:

where the coefficients $\alpha_j$ and $\beta_j$ are calculated to enforce orthonormality:

In matrix form, after $m$ steps, this recursion relation can be written as:

where $\mathbf{V}_m = [\mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_m]$ is the $N \times m$ matrix of Lanczos vectors, $\mathbf{e}_m = [0, 0, \dots, 1]^T$, and $\mathbf{T}_m$ is a symmetric, tridiagonal matrix of size $m \times m$:

2. Solving the Tridiagonal System

The eigenvalues $\theta_i$ (called Ritz values) and eigenvectors $\mathbf{y}_i$ of the tridiagonal matrix $\mathbf{T}_m$:

serve as approximations to the eigenvalues and eigenvectors of the original system. The approximate eigenvalues are $\lambda_i \approx 1/\theta_i$, and the approximate eigenvectors are:

Because $\mathbf{T}_m$ is tridiagonal, its eigenvalues can be computed using QR algorithms in a fraction of a second, even for $m = 1000$.

3. Loss of Orthogonality and Reorthogonalization

In finite precision arithmetic, the generated Lanczos vectors quickly lose their orthogonality due to round-off errors. Once orthogonality is lost, the algorithm starts generating redundant copies of already converged eigenvalues (called duplicate or spurious Ritz values).

To resolve this issue, modern codes implement reorthogonalization strategies:

$\mathbf{v}_{j+1}$ against all previously computed vectors. This is mathematically robust but computationally expensive, as it requires storing and processing all Lanczos vectors.

Ritz value converges, the corresponding Ritz vector is used to selectively orthogonalize the Lanczos vectors. This provides a balance between mathematical robustness and computational speed, and is the method of choice in commercial FEA solvers.

---

5.10 Advanced Component Mode Synthesis (CMS)

Component Mode Synthesis is a form of substructuring where a large finite element model is partitioned into several independent sub-assemblies (components). Each component is condensed using dynamic reduction, and the reduced components are then assembled at their boundaries to solve the global dynamic system.

5.10.1 Assembly and Coupling of Craig-Bampton Components

Let us consider a global structure consisting of two substructures, Component $A$ and Component $B$, which share a common boundary. The Craig-Bampton reduced equations of motion for each component are:

where the generalized coordinates are partitioned into boundary displacements and internal modal coordinates:

The uncoupled system of equations is:

To assemble the components, we apply the boundary compatibility and equilibrium conditions:

external forces:

We define a coordinate transformation matrix $\mathbf{T}_{\text{assembly}}$ that maps the uncoupled coordinates to the assembled coordinates:

The mapping is:

Using this transformation, the assembled global mass and stiffness matrices are:

Evaluating the products:

This assembled system represents the dynamics of the coupled structure. The Craig-Bampton approach is highly efficient for NVH engineering: the internal modes of components $A$ and $B$ can be computed once and stored, allowing for fast design iterations on the boundary connections (e.g., optimizing engine mount stiffnesses or suspension bushings).

5.10.2 MacNeal-Rubin Component Mode Synthesis (Free-Interface CMS)

While the Craig-Bampton method is the industry standard for numerical FEA models, it is difficult to implement experimentally. Clamping the boundaries of a heavy component (such as an engine block or suspension subframe) to enforce fixed-interface normal modes ($\mathbf{U}_b = \mathbf{0}$) is technically challenging and introduces test fixture errors.

To address this, MacNeal and Rubin developed the Free-Interface CMS method. In this method, the normal modes of the component are computed with free boundaries at the connection interfaces.

1. Residual Flexibility

Eliminating high-frequency modes (modal truncation) introduces errors in the static boundary compliance. Free-Interface CMS compensates for this by adding a **residual flexibility matrix** $\mathbf{G}_{\text{res}}$, which represents the static flexibility contribution of the truncated high-frequency modes.

Let the total static flexibility matrix be $\mathbf{G} = \mathbf{K}^{-1}$. We partition it using the modal expansion:

The second term on the right-hand side represents the residual flexibility:

By incorporating $\mathbf{G}_{\text{res}}$ at the boundary degrees of freedom, the free-interface CMS model achieves the same level of accuracy as the Craig-Bampton method while using modes that can be measured experimentally in a test lab.

5.11 Advanced Damping Formulations and Damping Treatments

In NVH engineering, predicting the amplitude of vibrations at resonance is a primary design task. Because the resonant response is limited only by energy dissipation, accurate damping models are essential.

5.11.1 Hysteretic Damping and Complex Stiffness

For steady-state harmonic vibrations, material damping is often modeled using Hysteretic Damping (also known as structural damping). Hysteretic damping assumes that the energy dissipated per cycle is independent of frequency and proportional to the square of the displacement amplitude.

This is implemented mathematically in the frequency domain by introducing a complex stiffness matrix $\mathbf{K}^*$:

where $\eta$ is the structural damping coefficient (loss factor) and $i = \sqrt{-1}$.

Substituting this into the equations of motion for harmonic excitation $\mathbf{F}(t) = \mathbf{F}_0 e^{i \omega t}$:

The imaginary term $i \eta \mathbf{K}$ represents the damping force, which is in phase with the velocity but proportional to the displacement.

The Causality Violation in the Time Domain

While the complex stiffness model is highly useful for frequency-domain harmonic analysis, it cannot be easily used in time-domain transient analysis. If we attempt to transform the complex stiffness term back to the time domain, we obtain a convolution integral:

This integral requires integrating over future times ($\tau > t$) to compute the damping force at the current time $t$. This violates the physical principle of causality (the cause must precede the effect). Thus, time-domain simulations must use viscoelastic or viscous damping models (such as Rayleigh damping or modal damping) instead.

5.11.2 Viscoelastic Damping Models

Viscoelastic materials exhibit both viscous (fluid-like) and elastic (solid-like) behaviors. We model this behavior by combining springs and dampers:

relation is:

In FEA, this leads to a damping matrix that is directly proportional to the stiffness matrix:

This model incorporates a relaxation time constant $\tau_{\text{rel}} = \eta_v / E$. In FEA, the Maxwell model requires introducing auxiliary state variables to track the internal stresses of the dampers, resulting in a system of differential-algebraic equations (DAEs).

5.11.3 Constrained Layer Damping (CLD) and Rao-Nakra Formulations

To control the vibration of thin sheet metal panels (such as vehicle doors, floors, and roofs), NVH engineers apply Constrained Layer Damping (CLD) treatments. A CLD treatment is a sandwich structure consisting of:

rubber-based polymer) with a high shear loss factor.

When the base structure bends, the stiff constraining layer forces the soft viscoelastic core to undergo shear deformation. This shear strain dissipates the vibrational energy as heat.

The analysis of CLD treatments is governed by the Rao-Nakra and Mead-Markus formulations:

longitudinal deformation, and the core undergoes only shear deformation. It leads to a sixth-order differential equation for the transverse displacement.

of all three layers, making it highly accurate for thick panels and high frequencies.

In finite element analysis, CLD treatments are modeled using a layer of solid elements (for the shearing viscoelastic core) sandwiched between shell elements (for the base and constraining layers), or by using specialized laminated composite shell elements that incorporate shear damping theories directly.

---

5.12 Convergence, Locking, and Mesh Density

The accuracy of dynamic finite element simulations depends on the element formulation, mesh density, and the frequency range of interest.

5.12.1 Element Locking Phenomena

Locking occurs when the numerical integration scheme or shape function interpolation introduces artificial constraints that make the element excessively stiff, preventing it from converging.

beams or Reissner-Mindlin plates) when the shape functions cannot represent zero shear strain under pure bending. It is resolved using Reduced Integration (using fewer Gauss integration points for the shear terms) or **Assumed Natural Strain (ANS)** methods.

bending and membrane deformations is uncoupled, causing bending to trigger spurious membrane stresses. It is resolved using higher-order shape functions or covariant formulations.

or metals undergoing plastic deformation, where Poisson's ratio $\nu \to 0.5$). The element stiffness matrix becomes ill-conditioned because the volume preservation constraint ($u_{i, i} = 0$) cannot be satisfied locally. It is resolved using selective-reduced integration (integrating the volumetric strain at the element center and the shear/deviatoric strain at the normal Gauss points) or **B-bar formulations**.

5.12.2 Mathematical Rates of Convergence of Eigenvalues

According to variational theory (the Rayleigh-Ritz theorem), the convergence rate of the calculated eigenvalues $\lambda_i^{(h)}$ as the element size $h \to 0$ is governed by the polynomial degree $p$ of the shape functions.

For a consistent mass matrix model, the eigenvalue error converges as:

where $p$ is the degree of the polynomial interpolation, and $m$ is the order of the highest derivative in the strain energy functional (e.g., $m=1$ for bars, $m=2$ for Euler-Bernoulli beams).

This explains why the 2-element CMM beam model (cubic Hermite shape functions) achieved an error of only +0.05% for the fundamental frequency: it converges at a rate of $O(h^4)$, meaning that halving the element length reduces the error by a factor of 16!

For a lumped mass matrix model, the convergence rate of eigenvalues is limited to:

regardless of the shape function degree. This occurs because the lumping process introduces a second-order approximation error ($O(h^2)$) to the kinetic energy representation, which dominates the overall convergence rate.

5.12.3 Mesh Refinement Strategies

To improve the accuracy of finite element simulations, engineers use three refinement strategies:

degree $p$ constant. This is the most common refinement strategy.

(e.g., changing from linear elements to quadratic or cubic elements) while keeping the mesh geometry constant. It offers faster convergence rates for smooth fields but is computationally expensive.

gradient/curvature (such as stress concentrations or wave fronts) without changing the number of elements or nodes.

5.13 Global coordinate transformations and boundary conditions

To analyze structures with complex geometries, individual elements are formulated in their own local coordinate systems. To assemble these elements into a complete global model, we must transform their stiffness and mass matrices into a single, unified global coordinate system.

5.13.1 Local vs. Global Coordinates

element (e.g., the local $x$-axis lies along the longitudinal axis of a bar or beam). This makes the integration of shape functions straightforward.

structure (e.g., the coordinate system of the vehicle chassis).

5.13.2 Coordinate Transformation Matrix in 2D

Consider a 2-node bar element in a 2D plane. In its local coordinate system, the element has two degrees of freedom: local axial displacements $u_{l1}$ and $u_{l2}$. In the global coordinate system, each node has two degrees of freedom: horizontal displacement $U$ and vertical displacement $V$. The vector of global nodal coordinates is:

Let the local longitudinal axis $x$ of the element be oriented at an angle $\theta$ relative to the global $X$-axis. The local axial displacement at node $i$ is the projection of the global displacements $U_i$ and $V_i$:

We can write this in matrix form:

where $\mathbf{T}_e$ is the coordinate transformation matrix of size $2 \times 4$.

1. Transformation of Element Matrices

The strain energy of the element must be invariant (equal) regardless of the coordinate system:

Thus, the element stiffness matrix in global coordinates $\mathbf{k}^g$ is:

Following the same energy invariance principle for kinetic energy:

2. Frame Rotation for Beams in 2D

For a beam element in 2D, each node has three degrees of freedom: axial displacement $U_i$, transverse displacement $V_i$, and bending rotation $\theta_i$ (which is invariant under 2D rotation). The transformation matrix for a 2-node beam element (6 DOFs total) is:

This allows for the assembly of frame structures (such as vehicle roll cages or space frames).

5.13.3 Boundary Condition Implementation Methods

Once the global stiffness matrix $\mathbf{K}$ and mass matrix $\mathbf{M}$ are assembled, they are singular (meaning their determinants are zero). This indicates that the structure has rigid-body modes and cannot support static loads. To solve the equations, we must apply boundary conditions to eliminate these rigid-body modes.

1. Direct Method (Row and Column Deletion)

If a degree of freedom $i$ is fixed (e.g., $u_i = 0$), we delete the $i$-th row and $i$-th column of the stiffness and mass matrices. This reduces the size of the system of equations. While mathematically exact, this method requires re-indexing the degrees of freedom, which can be computationally expensive to implement in parallel codes.

2. The Penalty Method

The penalty method enforces boundary conditions without changing the size of the global matrices. If the degree of freedom $u_i$ is constrained to a value $g_i$ (where $g_i = 0$ for a fixed boundary), we add a very large artificial stiffness value $C_p$ (called the penalty parameter) to the diagonal term $K_{ii}$ of the global stiffness matrix:

We also modify the global force vector by adding $C_p g_i$ to the $i$-th term:

The parameter $C_p$ is chosen to be much larger than the structural stiffness terms, typically:

This effectively models a very stiff spring connected to DOF $i$. The displacement converges to $g_i$ with an error of $O(1/C_p)$.

conditioning errors (the matrix becomes ill-conditioned, leading to round-off errors). If $C_p$ is too small, the boundary condition is not enforced accurately.

3. Lagrange Multiplier Method

The Lagrange multiplier method enforces boundary conditions exactly by introducing auxiliary variables (Lagrange multipliers $\lambda_j$, which physically represent the boundary reaction forces).

If we have $m$ linear constraints of the form $\mathbf{C} \mathbf{U} = \mathbf{g}$, we augment the system of equations:

This method enforces the constraints exactly and yields the reaction forces $\boldsymbol{\lambda}$. However, it increases the size of the system by $m$ equations and introduces zero diagonal blocks, which prevents the use of standard Cholesky solvers.

5.14 Modal Truncation and Dynamic Correction Methods

When solving the dynamic equations of motion of a large finite element model using modal coordinates, we must truncate the modal expansion to include only the lowest $P$ modes. This truncation introduces errors, particularly in the local stress and strain fields. To mitigate these errors, we use dynamic correction methods.

5.14.1 The Mode Displacement Method

In the Mode Displacement Method, the physical displacement vector $\mathbf{U}(t)$ is approximated using the truncated modal matrix $\boldsymbol{\Phi}_p = [\boldsymbol{\phi}_1, \boldsymbol{\phi}_2, \dots, \boldsymbol{\phi}_p]$ of size $N \times P$:

where the modal coordinates $q_i(t)$ are solved from the decoupled equations:

Limitations of the Mode Displacement Method

While the lowest modes accurately capture the global displacements, the high-frequency modes (which are truncated) act as localized stiffness springs under low-frequency excitation. By omitting these modes, the Mode Displacement Method fails to capture local deformations and stress concentrations near force application points. To achieve convergence for local stresses, an excessively large number of modes must be calculated, which is computationally expensive.

5.14.2 The Mode Acceleration Method

The Mode Acceleration Method (MAM) addresses the limitations of the Mode Displacement Method by incorporating the static flexibility contribution of the truncated high-frequency modes.

Step-by-Step Derivation of MAM

We begin with the original equation of motion:

For simplicity, let us assume an undamped structure ($\mathbf{C} = \mathbf{0}$). We can write the displacement vector by inverting the stiffness matrix:

Now, we approximate the acceleration vector $\mathbf{\ddot{U}}(t)$ using the truncated modal expansion:

Substituting this approximation into the equation for $\mathbf{U}(t)$:

We use the eigenvalue relationship for the mass-normalized eigenvectors:

Substituting this identity back into our expression:

This is the Mode-Acceleration Equation.

Physical and Computational Significance of MAM

The Mode-Acceleration equation contains two distinct terms:

static response of the structure to the external force $\mathbf{F}(t)$. This term is computed using the full, unreduced stiffness matrix, which means it implicitly incorporates the static contribution of all modes, including the truncated high-frequency ones.

\ddot{q}_i(t)$)**: Adjusts the response by subtracting the static contribution of the retained modes and adding back their dynamic inertial effects.

Because the static response is captured exactly, the Mode Acceleration Method converges much faster than the Mode Displacement Method. It is particularly effective for passenger cabin NVH models where structural panels behave quasi-statically in response to low-frequency engine and road excitations.

5.15 The Inertia Relief Method for Free-Free Structures

In automotive and aerospace NVH engineering, we frequently need to analyze the dynamic response of structures that are completely unsupported in space. Examples include an aircraft in flight subjected to aerodynamic gusts, a satellite in orbit under thruster loads, or a vehicle chassis suspended on soft bungee cords in an experimental shaker test rig.

Such structures are called free-free structures. Because they lack boundary constraints, their global stiffness matrix $\mathbf{K}$ is singular and cannot be inverted. Specifically, the nullspace of $\mathbf{K}$ has a dimension of $6$ (corresponding to the six rigid-body modes of three-dimensional translation and rotation):

where $\boldsymbol{\Phi}_r$ is the $N \times 6$ matrix of rigid-body modes.

To calculate the elastic deformations of a free-free structure under external dynamic forces $\mathbf{F}(t)$, we use the Inertia Relief Method.

5.15.1 Displacement Partitioning

We decompose the total displacement vector $\mathbf{U}(t)$ into a rigid-body motion component $\mathbf{U}_r(t)$ and an elastic deformation component $\mathbf{U}_e(t)$:

The rigid-body displacement is represented as a linear combination of the rigid-body modes:

where $\mathbf{q}_r(t)$ is the vector of six rigid-body coordinates.

The elastic deformation $\mathbf{U}_e(t)$ represents the deformation relative to the moving rigid-body reference frame. To make this decomposition unique, we impose the mathematical constraint that the elastic displacements must be mass-orthogonal to the rigid-body modes:

This constraint ensures that the elastic deformations do not introduce any linear or angular momentum relative to the rigid-body frame.

5.15.2 Derivation of the Inertia Relief Projection Matrix

We write the undamped equations of motion:

Substituting the partitioned displacements and noting that $\mathbf{K} \mathbf{U}_r = \mathbf{K} \boldsymbol{\Phi}_r \mathbf{q}_r = \mathbf{0}$:

To determine the rigid-body accelerations $\mathbf{\ddot{q}}_r$, we premultiply the equation of motion by the transpose of the rigid-body modes $\boldsymbol{\Phi}_r^T$:

We apply the orthogonality and boundary conditions:

energy).

derivative of the mass-orthogonality constraint).

The equation simplifies to:

The term $\mathbf{M}_{rr} = \boldsymbol{\Phi}_r^T \mathbf{M} \boldsymbol{\Phi}_r$ represents the $6 \times 6$ rigid-body mass matrix of the structure. Solving for the rigid-body accelerations:

The physical inertial force vector $\mathbf{F}_i$ opposing the rigid-body acceleration is:

We define the net force vector $\mathbf{F}_{\text{el}}$ causing elastic deformation as the sum of the external forces and the opposing inertial forces:

where the Inertia Relief Projection Matrix $\mathbf{P}$ is:

The projection matrix $\mathbf{P}$ satisfies the projection property $\mathbf{P}^2 = \mathbf{P}$. It filters out any components of the external forces that would cause rigid-body acceleration, leaving a self-equilibrated force system $\mathbf{F}_{\text{el}}$ that satisfies static equilibrium:

5.15.3 Solving for Elastic Deformations

The governing equation for the elastic deformations is:

subject to the constraint $\boldsymbol{\Phi}_r^T \mathbf{M} \mathbf{U}_e = \mathbf{0}$.

Because $\mathbf{K}$ is singular, we cannot solve this system directly. To obtain a solution, we temporarily constrain a set of six degrees of freedom (called support DOFs, denoted by $\mathbf{U}_s$) to eliminate rigid-body motion. We partition the coordinates into active master DOFs ($\mathbf{U}_a$) and support DOFs ($\mathbf{U}_s$):

Since $\mathbf{U}_s = \mathbf{0}$, we solve the first row of equations:

This yields a temporary displacement field $\mathbf{U}^* = [\mathbf{U}_a^*, \mathbf{0}]^T$.

Finally, we project $\mathbf{U}^*$ to satisfy the mass-orthogonality constraint, yielding the true elastic displacements:

The inertia relief method is a powerful tool in NVH, allowing engineers to calculate the static stress and deformation states of flying aircraft or cruising vehicles without introducing fictitious boundary constraints that would distort the results.

5.16 The B-bar Method for Volumetric Locking

In vehicle NVH engineering, elastomeric components—such as suspension bushings, engine mounts, and exhaust hangers—are widely used to isolate high-frequency vibrations. These components are made of rubber, which is nearly incompressible.

For an incompressible material, Poisson's ratio $\nu$ approaches $0.5$, which causes the bulk modulus $K$ to approach infinity:

Under these conditions, the strain energy associated with volumetric deformation becomes extremely large. The variational formulation forces the volumetric strain (dilatation) to vanish:

In standard finite element formulations, this incompressibility constraint is enforced at every integration point. For low-order elements, the number of volumetric constraints exceeds the number of nodal degrees of freedom, causing the element to lock (volumetric locking). The element exhibits artificially high stiffness, preventing it from deforming.

To resolve volumetric locking, we use the B-bar ($\bar{\mathbf{B}}$) Method, which decouples the deviatoric (shear) and volumetric (dilatational) deformations.

5.16.1 Strain and B-Matrix Splitting

We decompose the strain vector $\boldsymbol{\epsilon}$ into a deviatoric component $\boldsymbol{\epsilon}_{\text{dev}}$ and a volumetric component $\boldsymbol{\epsilon}_{\text{vol}}$:

where:

and the vector $\mathbf{m}$ is defined for 3D stress states as:

The dilatation $\theta$ is the sum of the normal strains:

The strain-displacement relationship is $\boldsymbol{\epsilon} = \mathbf{B} \mathbf{U}$. We split the strain-displacement matrix $\mathbf{B}$ into deviatoric and volumetric parts:

where:

5.16.2 Formulation of the Modified B-bar Matrix

To prevent locking, the B-bar method replaces the volumetric strain-displacement matrix $\mathbf{B}_{\text{vol}}$ with a modified volumetric matrix $\bar{\mathbf{B}}_{\text{vol}}$. This modified matrix is obtained by projecting the volumetric strain onto a lower-order approximation space, typically by evaluating the volumetric strains at the element center (centroid) or by using a lower-order integration scheme.

The modified strain-displacement matrix $\bar{\mathbf{B}}$ is defined as:

where:

and $\bar{\mathbf{m}}^T \mathbf{B}$ represents the spatially averaged volumetric strain. For a 4-node quadrilateral element, $\bar{\mathbf{B}}_{\text{vol}}$ is constant over the element, evaluated at the centroid:

The stiffness matrix is then computed by separating the deviatoric and volumetric strain energies:

where $\mathbf{D}_{\text{dev}}$ contains the shear moduli and $\mathbf{D}_{\text{vol}}$ contains the bulk modulus $K$.

By averaging the volumetric strain over the element, we reduce the number of independent incompressibility constraints, softening the volumetric stiffness and resolving volumetric locking. This allows for accurate modeling of elastomeric bushings and engine mounts in vehicle dynamic simulations.

---

5.17 Craig-Bampton Substructuring Assembly with Elastic Connectors

In vehicle system dynamics, dynamic components (such as the engine block and the vehicle body-in-white) are represented by Craig-Bampton reduced-order models. These components are connected not by rigid links, but by elastic and dissipative elements (such as rubber engine mounts).

5.17.1 Component Formulation

Let Component $A$ represent the engine block and Component $B$ represent the vehicle body. Each component is reduced using the Craig-Bampton method:

where $\mathbf{R}_c^A$ and $\mathbf{R}_c^B$ are the reaction forces at the connection boundaries, and the coordinates are:

5.17.2 Modeling the Elastic Connector

Let the two components be connected between boundary degree of freedom $\mathbf{u}_b^A$ and boundary degree of freedom $\mathbf{u}_b^B$ by an elastic mount of stiffness $k_m$ and damping coefficient $c_m$.

The interaction forces acting on the boundaries are proportional to the relative displacement and velocity:

We can write this interaction force system in matrix form relative to the uncoupled coordinates:

where the mount stiffness matrix $\mathbf{K}_m$ acts only on the boundary degrees of freedom:

acting on the partition $[\mathbf{u}_b^A, \mathbf{u}_b^B]^T$.

5.17.3 Assembled System Equations

By moving the interaction forces to the left-hand side, the assembled equations of motion for the coupled engine-mount-body system are:

where the coupling stiffness matrix $\mathbf{K}_{\text{mount}}$ is:

with the rows and columns mapped to the coordinates $[\mathbf{u}_b^A, \mathbf{q}_n^A, \mathbf{u}_b^B, \mathbf{q}_n^B]^T$.

This formulation allows NVH engineers to analyze the vibration isolation efficiency of the engine mounts, evaluate cabin noise levels, and optimize mount properties using a highly condensed, computationally efficient model.

5.18 Modal Effective Mass in Multi-Directional 3D Space

In practical NVH design, structures are subjected to multi-directional dynamic loads. A vehicle body panel is excited by longitudinal, lateral, and vertical engine forces, and a satellite is subjected to three-axis accelerations during launch. To evaluate the dynamic participation of modes under multi-directional base excitation, we extend the modal effective mass concept to three-dimensional space.

5.18.1 Multi-Directional Influence Matrix

Instead of a single influence vector $\mathbf{r}$, we define a multi-directional influence matrix $\mathbf{R}$ of size $N \times 6$ (or $N \times 3$ for translation only), where each column represents the static displacement of the degrees of freedom due to a unit displacement or rotation of the base in one of the six rigid-body directions:

where:

rotational influence vectors.

5.18.2 Modal Participation Vector

For the $i$-th mass-normalized mode, we define the Modal Participation Vector $\boldsymbol{\Gamma}_i$ of size $6 \times 1$ as:

where:

5.18.3 The $6 \times 6$ Modal Effective Mass Matrix

For the $i$-th mass-normalized mode, the Modal Effective Mass Matrix $\mathbf{M}_{\text{eff}, i}$ of size $6 \times 6$ is defined as:

Physical Interpretation of Matrix Terms

effective mass participating in mode $i$ due to excitation in the $k$-th coordinate direction.

Represent the inertial coupling between the $k$-th and $l$-th directions in mode $i$. If the center of gravity of the structure is offset from the coordinate origin, a translation of the base in the $X$-direction excites bending and rotational modes, generating off-diagonal coupling terms.

5.18.4 Bounding of the Rigid-Body Mass Properties

By applying the completeness relation of the eigenvectors, we can prove that the sum of the modal effective mass matrices over all modes equals the global **Rigid-Body Mass Matrix** $\mathbf{M}_{\text{rb}}$ of the structure:

The matrix $\mathbf{M}_{\text{rb}}$ is a $6 \times 6$ symmetric matrix containing the fundamental rigid-body mass properties of the structure:

where:

coordinates of the center of gravity ($x_g, y_g, z_g$).

I_{yy}, I_{zz}$) and products of inertia ($I_{xy}, I_{yz}, I_{zx}$) relative to the origin.

This relationship connects modal dynamics to fundamental rigid-body physics. In launch vehicle design, this identity is monitored to ensure that the truncated modal models retain at least 90% of the rigid-body mass and moments of inertia in all six coordinate axes.

5.19 Powertrain Mounting Systems and Torque Roll Axis (TRA) Decoupling

The vehicle powertrain (engine and transmission assembly) is a major source of dynamic forces. These forces arise from reciprocating piston unbalances, combustion pressure fluctuations, and gear whine. To isolate the passenger cabin from these vibrations, the powertrain is supported by elastic engine mounts.

From a structural dynamics perspective, the engine block is modeled as a rigid body suspended on $M$ discrete rubber mounts. A rigid body in space has six degrees of freedom. The displacement vector of the engine's center of gravity (CG) is:

where:

5.19.1 Engine Mass Matrix Formulation

The $6 \times 6$ mass matrix $\mathbf{M}_{\text{eng}}$ relative to the CG coordinate system is:

where:

align with the principal inertial axes of the engine, the products of inertia vanish, and the mass matrix becomes diagonal.

5.19.2 Mount Stiffness Transformation

Consider the $k$-th mount located at position $\mathbf{r}_k = [x_k, y_k, z_k]^T$ relative to the CG. In the local coordinate system of the mount, the stiffness is diagonal:

where $u, v, w$ represent the local axes of the mount.

Let $\mathbf{R}_k$ be the $3 \times 3$ rotation matrix that transforms coordinates from the global CG frame to the local mount frame. The local translation vector at the mount $\mathbf{u}_k$ is related to the global CG displacements through the kinematic relation:

Writing the cross-product using a skew-symmetric matrix:

where $\mathbf{B}_k$ is the $3 \times 6$ kinematic transformation matrix:

The contribution of the $k$-th mount to the global stiffness matrix at the CG is:

The total stiffness matrix of the mounted powertrain is the sum of the contributions of all $M$ mounts:

5.19.3 Torque Roll Axis (TRA) Decoupling

When the engine fires, the fluctuating pressure on the pistons exerts a dynamic torque about a specific axis, known as the Torque Roll Axis (TRA). If the rigid-body roll mode of the engine is coupled with its translation modes, the dynamic torque will excite lateral and vertical translations, which are transmitted through the mounts as forces. This generates cabin hum and vibration.

To minimize this vibration transmission, NVH engineers optimize the mounting system to achieve Torque Roll Axis Decoupling.

Decoupling is achieved when the stiffness coupling terms in the total stiffness matrix $\mathbf{K}_{\text{total}}$ between the roll coordinate (rotation about the TRA) and the translation coordinates are zero:

By aligning one of the system's principal stiffness axes with the Torque Roll Axis, the dynamic engine torque excites only the roll mode, preventing lateral or vertical accelerations. This is achieved by tuning:

hardnesses.

optimal coordinates.

relative to the global coordinates.

Decoupling the powertrain modes represents a critical NVH design milestone, significantly improving vehicle refinement and cabin comfort.

5.20 The Eigensystem Realization Algorithm (ERA) for Experimental Modal Identification

In NVH engineering, we must validate our finite element models using experimental test data. While FEA provides numerical natural frequencies and mode shapes, dynamic testing measures physical acceleration, velocity, or displacement responses. To bridge the gap, we use modal identification algorithms to extract experimental modal parameters (frequencies, damping ratios, and mode shapes) from measured time-domain impulse response functions or frequency response functions.

The Eigensystem Realization Algorithm (ERA), developed by Juang and Pappa in 1985, is a powerful, system-identification technique widely used to extract structural modal parameters.

5.20.1 Hankel Matrix Construction

ERA begins with the discrete-time state-space representation of the structural dynamic system:

where $\mathbf{x}(k)$ is the state vector of size $2N \times 1$ at time step $k$, $\mathbf{u}(k)$ is the input vector (excitation force), and $\mathbf{y}(k)$ is the output vector (sensor responses).

For an impulse excitation at $k=0$ (e.g., an impact hammer hit), the output response matrix at subsequent time steps (Markov parameters) is:

where $\mathbf{Y}(k)$ is a matrix of size $N_y \times N_u$ ($N_y$ output sensors, $N_u$ excitation inputs).

We construct two large Hankel Matrices, $\mathbf{H}_0$ and $\mathbf{H}_1$, of size $r N_y \times s N_u$:

where $\mathbf{H}_1$ is shifted by one time step relative to $\mathbf{H}_0$.

5.20.2 Singular Value Decomposition and Model Realization

We perform a Singular Value Decomposition (SVD) on the baseline Hankel matrix $\mathbf{H}_0$:

where $\mathbf{U}$ and $\mathbf{V}$ are orthonormal matrices, and $\mathbf{\Sigma}$ is the diagonal matrix of singular values sorted in descending order:

If the structure has $M$ active modes within the measured frequency band, the state-space system has order $n = 2M$. In the absence of noise, only the first $n$ singular values are non-zero. We truncate the SVD matrices to this active subspace:

where $\mathbf{\Sigma}_n$ is of size $n \times n$.

The realized state-space matrices $\mathbf{A}, \mathbf{B}, \mathbf{C}$ are computed as:

where $\mathbf{E}_u$ and $\mathbf{E}_y$ are selection matrices that extract the input and output blocks.

5.20.3 Extracting Modal Parameters

To find the natural frequencies and damping ratios, we solve the eigenvalue problem of the realized state matrix $\mathbf{A}$:

where $z_i$ are the discrete-time complex eigenvalues (poles).

We map these discrete-time poles back to continuous-time poles $s_i$:

where $\Delta t$ is the sampling time step. The modal parameters are extracted from the continuous pole $s_i$:

The columns of the matrix product $\mathbf{C} \boldsymbol{\Psi}$ represent the experimental mode shapes (where $\boldsymbol{\Psi}$ contains the eigenvectors of $\mathbf{A}$).

---

5.21 Modal Assurance Criterion (MAC) for FEA-Test Model Correlation

Once experimental mode shapes are extracted using ERA, we must correlate them with the mode shapes predicted by the FEA model. Since experimental sensors are only placed at a subset of the structural nodes, we define a set of correlation degrees of freedom.

The primary tool to correlate mode shapes is the Modal Assurance Criterion (MAC).

5.21.1 Mathematical Formulation of MAC

The MAC value between the $i$-th FEA mode shape $\boldsymbol{\phi}_{\text{FEA}, i}$ and the $j$-th test mode shape $\boldsymbol{\phi}_{\text{Test}, j}$ is defined as:

The MAC is a scalar value between $0.0$ and $1.0$:

linearly dependent, differing only by a scaling factor).

orthogonal (uncorrelated).

In practice, a MAC value above 0.9 for diagonal terms ($i=j$) indicates acceptable correlation between the FEA model and the test data. The off-diagonal terms ($i \neq j$) should be below 0.1, which confirms that the modes are orthogonal and decouples the dynamic responses.

5.21.2 Spatial Aliasing and Sensor Placement

If too few sensor points are used in the test, the discrete mode shapes may lose their orthogonality. This numerical anomaly is called Spatial Aliasing. High-frequency modes may appear identical to low-frequency modes because the sensors cannot resolve the spatial variation between them. This results in spurious off-diagonal MAC values close to 1.0.

To prevent spatial aliasing, engineers optimize the sensor layout using the Coordinate Modal Assurance Criterion (CoMAC), which identifies which degrees of freedom contribute most to the correlation errors, guiding the placement of additional accelerometers.

5.22 State-Space Formulation and Complex Modes for Non-Proportional Damping

In classical structural dynamics, we frequently assume that the damping matrix $\mathbf{C}$ is proportional to a linear combination of the mass and stiffness matrices (Rayleigh damping). Under this assumption, the real normal mode shapes $\boldsymbol{\phi}_i$ of the undamped system diagonalize all three global matrices: mass $\mathbf{M}$, stiffness $\mathbf{K}$, and damping $\mathbf{C}$. This results in a set of decoupled second-order differential equations that can be solved independently.

However, in many practical engineering designs, damping is highly non-proportional. Examples include:

engine mounts, or localized viscoelastic dampers.

capacity (viscoelastic layers) and others have low damping (carbon fiber or steel).

acoustic damping at the boundaries.

friction generates unsymmetric joint stiffness and damping.

When the damping is non-proportional, the modal damping matrix $\mathbf{C}_q = \boldsymbol{\Phi}^T \mathbf{C} \boldsymbol{\Phi}$ contains non-zero off-diagonal terms. These terms couple the classical normal modes, meaning that a force exciting one mode will transfer energy to other modes through the damping mechanism. To analyze these systems, we must use a first-order state-space formulation.

---

5.22.1 The Duncan State-Space Formulation

We begin with the second-order equation of motion for an $N$-degree-of-freedom system:

To convert this into a set of $2N$ first-order differential equations, we define the state vector $\mathbf{y}(t)$ of size $2N \times 1$:

We also write a trivial identity matching the velocity vector:

Combining the original equation of motion and this identity, we obtain the state-space system:

Let us define the state-space matrices $\mathbf{A}$ and $\mathbf{B}$ and the state excitation vector $\mathbf{F}(t)$:

Both $\mathbf{A}$ and $\mathbf{B}$ are symmetric matrices of size $2N \times 2N$, assuming that $\mathbf{M}$, $\mathbf{C}$, and $\mathbf{K}$ are symmetric. The state-space equation of motion simplifies to:

---

5.22.2 The Complex Eigenvalue Problem

For free vibrations ($\mathbf{F}(t) = \mathbf{0}$), we assume a solution of the form:

Substituting this into the homogeneous state-space equation yields the generalized eigenvalue problem:

This system has $2N$ eigenvalues $\lambda_r$ and corresponding eigenvectors $\boldsymbol{\psi}_r$. Since the matrices $\mathbf{A}$ and $\mathbf{B}$ are real, any complex eigenvalues and eigenvectors must occur in complex conjugate pairs.

For a stable, underdamped system, the eigenvalues are:

where:

to the modal damping ratio.

The undamped natural frequency $\omega_r$ and damping ratio $\zeta_r$ are extracted via:

---

5.22.3 Physical Interpretation of Complex Mode Shapes

The state-space eigenvector $\boldsymbol{\psi}_r$ has the structure:

where $\boldsymbol{\phi}_r$ is the complex displacement mode shape of size $N \times 1$.

Unlike real mode shapes where all degrees of freedom oscillate in phase or $180^\circ$ out of phase (producing standing waves with fixed nodes), complex mode shapes introduce phase angles between different degrees of freedom:

where $| \phi_{r, j} |$ is the amplitude of the $j$-th degree of freedom in the $r$-th mode, and $\theta_{r, j}$ is its phase angle.

This phase difference has profound physical implications:

stationary; they travel along the structure.

through the structure from regions of low damping to regions of high damping, where it is dissipated.

displacement at any instant. Instead, it exhibits a continuous waving motion.

---

5.22.4 Orthogonality and Modal Decoupling

Because the state-space matrices $\mathbf{A}$ and $\mathbf{B}$ are symmetric, the complex eigenvectors satisfy the following orthogonality relations:

We normalize the eigenvectors with respect to the matrix $\mathbf{A}$ such that:

Using the modal matrix $\boldsymbol{\Psi} = [\boldsymbol{\psi}_1, \boldsymbol{\psi}_2, \dots, \boldsymbol{\psi}_{2N}]$ and modal coordinates $\mathbf{q}(t)$, we transform the state vector:

Substituting this transformation into the state-space equation of motion and pre-multiplying by $\boldsymbol{\Psi}^T$ yields $2N$ uncoupled first-order differential equations:

where the generalized modal force is:

Integrating these equations yields the dynamic response of the system in terms of its complex modes:

---

5.22.5 Applications in NVH: Brake Squeal Analysis

Complex mode analysis is the industry standard for studying Brake Squeal (high-frequency friction-induced vibration). In a disc brake assembly, the friction forces between the rotor and pad introduce unsymmetric terms in the stiffness matrix ($\mathbf{K} \neq \mathbf{K}^T$).

This asymmetry modifies the eigenvalue problem. The real parts of some eigenvalues can become positive ($\sigma_r > 0$). In the time domain, a positive real part indicates an exponentially growing oscillation:

This dynamic instability is the root cause of squeal. NVH engineers use complex eigenvalue analysis to: