2. Governing equations of structural dynamics
Every vibration qualification problem — whether it is a sine sweep on an avionics bracket or a random vibration test on a satellite panel — ultimately reduces to one matrix equation. Understanding where that equation comes from, what each matrix physically represents, and how it is solved numerically is the foundation for interpreting any modal or harmonic FEA result. This section develops the semi-discrete equation of motion, examines the mass, damping, and stiffness matrices individually, and then shows how modal superposition transforms a large coupled system into a set of independent single-DOF problems that can be evaluated analytically.
From Continuous Structures to Discrete Matrices
For a linear elastic body, the strong form of dynamic equilibrium is written in terms of a continuous displacement field u(x, t):
ρ ü(x,t) + c u̇(x,t) + ∇·σ(u) = f(x,t)
where ρ is density, c is a distributed damping coefficient, σ is the stress tensor, and f is the applied body force. The finite element method discretizes this field into nodal degrees of freedom u(t) ∈ ℝn using shape functions N(x), so that u(x,t) ≈ N(x)u(t). Applying the Galerkin weighted-residual procedure and integrating by parts yields the semi-discrete equation of motion:
M ü(t) + C u̇(t) + K u(t) = F(t)
This is the governing equation for all structural dynamic analysis in FEA. The term "semi-discrete" means space is discretized but time remains continuous — time integration or frequency-domain transformation is applied afterward. The three matrices M, C, and K are each n×n and symmetric for standard formulations.
Physical Meaning of the Three Matrices
Mass matrix M. The consistent mass matrix is assembled from element contributions Me = ∫Ve ρ NTN dV. It is positive definite and fully populated within each element, coupling all nodal DOFs through the shape functions. In production solvers, a lumped (diagonal) mass matrix is often used instead because it reduces storage and enables explicit time integration. The trade-off is slightly lower accuracy in higher modes. For qualification work, the consistent formulation is preferred when accurate modal frequencies above roughly the fifth mode are needed.
Stiffness matrix K. Assembled as Ke = ∫Ve BTDB dV, where B is the strain-displacement matrix and D is the constitutive matrix. It is symmetric and positive semi-definite; rigid-body modes produce zero eigenvalues. Boundary conditions remove these singularities. In qualification models, stiffness dominates the low-frequency response, and any error in joint stiffness or material modulus directly shifts the predicted natural frequencies.
Damping matrix C. This is the most problematic term. Physical damping mechanisms — material hysteresis, joint friction, acoustic radiation — are not easily expressed in a finite element form. The exact damping matrix would require knowledge of energy dissipation at every material point, which is rarely available. Instead, engineering practice relies on phenomenological models, most commonly Rayleigh damping:
C = α M + β K
where α and β are calibrated from two target modal damping ratios ζi and ζj at frequencies ωi and ωj:
α = 2ω_i ω_j (ζ_j ω_i − ζ_i ω_j) / (ω_i² − ω_j²)
β = 2(ζ_i ω_i − ζ_j ω_j) / (ω_i² − ω_j²)
For a typical aluminum spacecraft structure, ζ = 0.02 to 0.05 is common; for bolted joints with friction, ζ can reach 0.05–0.10. The choice of the two anchor frequencies matters: they should bracket the frequency range of interest. A common mistake is to anchor at the first two modes, which produces excessive damping at higher frequencies and under-predicts resonant response in the qualification band.
Modal Superposition
The undamped eigenvalue problem is obtained by setting C = 0 and F = 0 in the equation of motion, assuming harmonic motion u = φejωt:
(K − ω² M) φ = 0
This generalized eigenvalue problem yields n eigenvalues ωr² and n eigenvectors φr. The eigenvectors are M-orthogonal and K-orthogonal:
φ_r^T M φ_s = 0, φ_r^T K φ_s = 0 for r ≠ s
φ_r^T M φ_r = m_r (modal mass)
φ_r^T K φ_r = k_r = ω_r² m_r (modal stiffness)
Modal superposition exploits this orthogonality by expressing the physical response as a linear combination of mode shapes:
u(t) = Σ_r φ_r q_r(t) = Φ q(t)
Substituting into the equation of motion and pre-multiplying by ΦT decouples the system into n independent SDOF equations:
m_r q̈_r + c_r q̇_r + k_r q_r = f_r(t), r = 1, 2, ..., n
where fr(t) = φrTF(t) is the modal force. In practice, only the first m modes (m ≪ n) are retained because higher modes contribute negligibly to the response in the frequency band of interest. The truncation error depends on the modal participation of the discarded modes and the frequency range of the excitation.
Harmonic Response via Modal Superposition
For harmonic excitation F(t) = F0ejωt, the steady-state response is u(t) = U(ω)ejωt. Substituting into the equation of motion gives the frequency-domain formulation:
(−ω² M + jω C + K) U(ω) = F₀
The term in parentheses is the dynamic stiffness matrix Z(ω). Direct solution of this complex linear system at each frequency is possible but expensive for large models. Modal superposition offers a far more efficient alternative: project the force onto the modal basis, solve each SDOF equation analytically, and sum the contributions.
For mode r, the modal coordinate response is:
q_r(ω) = f_r(ω) / (k_r − ω² m_r + jω c_r)
Dividing numerator and denominator by mr and introducing the modal damping ratio ζr = cr/(2mrωr) gives the standard form:
q_r(ω) = (f_r / m_r) / (ω_r² − ω² + 2jζ_r ω_r ω)
The physical response is then recovered as U(ω) = Σr φr qr(ω). This is the modal frequency response function (FRF). Each mode contributes a resonant peak at ω ≈ ωr with half-power bandwidth Δω = 2ζrωr. The quality factor Qr = 1/(2ζr) determines the sharpness of the peak — a critical parameter when planning a sine sweep test, because the sweep rate must be slow enough to allow the structure to reach steady state at each resonance.
Practical Implications for Qualification Testing
Three consequences of this formulation directly affect how a qualification test is specified and interpreted:
- Mode truncation. If the excitation spectrum has significant energy above the highest retained mode, the missing static contribution must be recovered. The residual flexibility method adds a static correction term K−1F minus the static contribution of the retained modes. Without this correction, the response at low frequencies can be under-predicted by 5–15% in typical aerospace structures.
- Damping calibration. The Rayleigh coefficients α and β must be chosen so that the modal damping ratios match measured or specified values across the qualification band. A common approach is to anchor ζ at the first mode and at a frequency near the upper bound of the test spectrum, then verify that intermediate modes have physically reasonable damping.
- Resonance tracking. During a sine sweep, the response amplitude at each resonance depends on ζr. If the FEA model uses ζ = 0.02 but the actual structure exhibits ζ = 0.01 due to a bolted joint that is tighter than modeled, the predicted peak response will be half the actual value — a non-conservative error that can lead to under-testing.
Summary of the Solution Chain
The complete workflow from physical structure to qualification-ready response prediction follows this chain: assemble M, C, and K from element contributions; solve the undamped eigenvalue problem for ωr and φr; calibrate Rayleigh damping coefficients from target modal damping ratios; project the harmonic force onto the modal basis; solve each SDOF equation analytically; and superpose the modal responses to recover the physical displacement field. Each step introduces assumptions — consistent vs. lumped mass, proportional vs. non-proportional damping, mode truncation — that must be documented and justified in the qualification report. The next section applies this framework to a concrete bracket model and compares the modal superposition result against direct frequency response.