2. Governing equations for transient structural dynamics

Transient structural dynamics rests on a single physical statement: at every instant, the internal resisting forces, the inertial forces, and the applied loads must balance. Translating that statement into a solvable numerical problem requires three successive reformulations — the strong form (a pointwise partial differential equation), the weak form (an integral statement), and the semi-discrete equations of motion (a matrix ordinary differential equation in time). Understanding how each level is derived, and what each one assumes, is what allows you to set up an initial-boundary value problem (IBVP) that a finite element or finite difference solver can actually integrate.

Strong form: pointwise equilibrium

Consider a deformable body occupying a domain Ω with boundary Γ. Let u(x, t) denote the displacement field, ρ the mass density, σ the Cauchy stress tensor, and b the body force per unit volume. Newton's second law applied to an infinitesimal material element gives the local momentum balance

ρ ü_i = σ_ij,j + b_i        in Ω × (0, T]

where the double dot denotes the second material time derivative and the comma denotes spatial differentiation. For linear elasticity, the constitutive law closes the system: σ_ij = C_ijkl ε_kl, with the small-strain tensor ε_kl = ½(u_k,l + u_l,k). Substituting yields the Navier–Cauchy equations, a set of coupled second-order hyperbolic PDEs in space and time.

The strong form is completed by boundary and initial conditions. On the Dirichlet (essential) portion Γ_u, displacements are prescribed: u_i = ū_i. On the Neumann (natural) portion Γ_t, tractions are prescribed: σ_ij n_j = t̄_i, where n is the outward unit normal. Because the equation is second order in time, two initial conditions are required: u_i(x, 0) = u_i⁰ and u̇_i(x, 0) = v_i⁰. The strong form demands that u be twice differentiable in space and time — a requirement that real geometries with corners, cracks, or material interfaces routinely violate.

Weak form: integral statement and natural boundary conditions

The weak form relaxes the differentiability requirement by multiplying the strong form with a test function w (also called the virtual displacement δu), integrating over Ω, and applying the divergence theorem. The result is the principle of virtual work:

∫_Ω ρ w_i ü_i dΩ + ∫_Ω w_i,j C_ijkl u_k,l dΩ
    = ∫_Ω w_i b_i dΩ + ∫_Γt w_i t̄_i dΓ

Three properties make this form the foundation of every implicit and explicit FE code. First, it requires only first derivatives of u, so C⁰-continuous shape functions are admissible. Second, the traction boundary condition appears directly in the integral and is therefore satisfied in an average (natural) sense rather than pointwise. Third, the mass and stiffness operators emerge naturally as symmetric bilinear forms, which is what guarantees the symmetry of the assembled matrices.

Semi-discrete equations of motion

Discretizing the weak form with shape functions N_a(x) and nodal degrees of freedom d(t) gives the semi-discrete system

M d̈ + C ḋ + K d = F_ext(t)

Here M is the consistent mass matrix, K the stiffness matrix, C the damping matrix, and F_ext the assembled external load vector. The term "semi-discrete" is precise: space has been discretized, but time remains continuous. The system is a set of N coupled second-order ODEs, where N is the total number of unconstrained degrees of freedom.

The mass matrix deserves particular attention because it determines which time integrator you can use. A consistent mass matrix is obtained by integrating ρ N_a N_b over each element; it is fully populated and couples neighbouring nodes. A lumped mass matrix is diagonal, obtained by row-summing or scaling the consistent matrix. Lumping introduces a small frequency error but makes the inverse M⁻¹ trivial, which is the enabling trick behind explicit central-difference integration. For impact and blast problems, where the stable time step is already tiny, lumped mass is almost always the right choice.

Damping is rarely derived from first principles. Rayleigh damping, C = αM + βK, is the standard engineering model; α controls low-frequency dissipation and β controls high-frequency dissipation. In explicit impact simulations, artificial bulk viscosity is often added instead, to smear shock fronts and suppress spurious oscillations without polluting the physical response.

Setting up the IBVP in practice

A complete, runnable IBVP specification must fix six items. Omitting any one of them produces a solver that either fails to start or returns a plausible but wrong answer.

As a concrete example, consider a 1 m long steel bar (ρ = 7850 kg/m³, E = 210 GPa) struck at one end by a 10 MPa compressive pulse of 100 µs duration. The wave speed is c = √(E/ρ) ≈ 5172 m/s, so the pulse travels the bar in about 193 µs. A mesh with h = 5 mm gives an element transit time of roughly 0.97 µs; with a lumped mass matrix and central difference, the critical step is Δt ≤ h/c ≈ 0.97 µs, and a safety factor of 0.8 gives Δt ≈ 0.77 µs. The bar is fixed at the far end, so a reflected tensile wave is expected at t ≈ 193 µs. If the simulation instead shows the pulse arriving at 100 µs, the most likely cause is a unit error in E or ρ; if it shows no reflection, the fixed boundary was applied to the wrong node set.

From continuous to discrete: the full chain

The diagram below traces how the same physical statement is transformed at each level, and which quantities are carried forward.

The chain is not merely a derivation sequence — it is a design checklist. If your simulation misbehaves, the fault almost always lies at one specific link: an ill-posed strong form (missing BC), an inconsistent weak form (wrong test space), a mis-assembled semi-discrete system (singular K or wrong M), or an unstable time integrator (Δt too large). Isolating which link is broken is far faster than re-running the whole model with adjusted parameters.

Verification before production runs

Before trusting any transient result, run three cheap checks. First, a rigid-body test: apply a uniform initial velocity with no constraints and confirm that the kinetic energy remains constant and no stress develops. Second, a patch test: apply a linear displacement field and verify that the recovered stress is constant to machine precision. Third, a wave-speed check: measure the arrival time of a known pulse and compare it with c = √(E/ρ). These three tests catch the overwhelming majority of IBVP setup errors — wrong units, missing constraints, and incorrect mass scaling — before they contaminate a design decision.