Transient time integrators ========================== Modal superposition (:doc:`eigenproblem`) handles the harmonic case efficiently, but transient analysis — a load history with arbitrary time content, an initial-condition sweep, a shock- response step — requires direct integration of the second-order ODE in time. This page derives the Newmark-:math:`\beta` family femorph-solver ships, motivates the default parameter choice, and lists the stability / accuracy regions a user needs to pick the right :math:`\beta` and :math:`\gamma` for the problem at hand. The driving API is :meth:`~femorph_solver.Model.solve_transient`; see :doc:`/user-guide/solving/transient` for the user-facing walk-through. Setting ------- The semi-discrete equation of motion after FE discretisation is .. math:: \mathbf{M}\, \ddot{\mathbf{u}}(t) + \mathbf{C}\, \dot{\mathbf{u}}(t) + \mathbf{K}\, \mathbf{u}(t) \;=\; \mathbf{F}(t), with :math:`\mathbf{M}, \mathbf{K} \in \mathbb{R}^{N \times N}` the assembled mass and stiffness, :math:`\mathbf{C}` an optional damping matrix (Rayleigh, modal, or directly assembled), and :math:`\mathbf{F}(t)` the time-varying load. We seek the discrete history :math:`\mathbf{u}^n \approx \mathbf{u}(t_n)` on the time grid :math:`t_n = n \Delta t`. Newmark-:math:`\beta` family ---------------------------- Newmark 1959 (*ASCE J. Eng. Mech. Div.* 85, 67–94) proposed the implicit single-step update .. math:: \mathbf{u}^{n+1} \;=\; \mathbf{u}^n + \Delta t\, \dot{\mathbf{u}}^n + \tfrac{\Delta t^2}{2}\, \bigl[ (1 - 2\beta)\, \ddot{\mathbf{u}}^n + 2\beta\, \ddot{\mathbf{u}}^{n+1} \bigr], .. math:: \dot{\mathbf{u}}^{n+1} \;=\; \dot{\mathbf{u}}^n + \Delta t\, \bigl[ (1 - \gamma)\, \ddot{\mathbf{u}}^n + \gamma\, \ddot{\mathbf{u}}^{n+1} \bigr], with two free parameters :math:`(\beta, \gamma)` that pick the specific scheme. Substituting back into the equation of motion at :math:`t_{n+1}` and solving for :math:`\ddot{\mathbf{u}}^{n+1}` gives a single sparse linear system per step: .. math:: \bigl( \mathbf{M} + \gamma \Delta t\, \mathbf{C} + \beta \Delta t^2\, \mathbf{K} \bigr)\, \ddot{\mathbf{u}}^{n+1} \;=\; \mathbf{F}^{n+1} - \mathbf{C}\, \tilde{\dot{\mathbf{u}}}^{n+1} - \mathbf{K}\, \tilde{\mathbf{u}}^{n+1}, where :math:`\tilde{\mathbf{u}}^{n+1}` and :math:`\tilde{\dot{\mathbf{u}}}^{n+1}` are the predictor terms that depend only on :math:`(\mathbf{u}^n, \dot{\mathbf{u}}^n, \ddot{\mathbf{u}}^n)`. The "effective" :math:`(\mathbf{M} + \gamma \Delta t\, \mathbf{C} + \beta \Delta t^2\, \mathbf{K})` matrix is the same at every step (when :math:`\Delta t` is constant), so it factorises *once* and back-solves for the remaining :math:`n_\text{steps} - 1` updates — the dominant cost is the single up-front factorisation. That single-factorisation property is what makes implicit Newmark cheap on long histories, and is why femorph-solver's :meth:`Model.solve_transient` factorises once and reuses the factor; passing a non-uniform time step would force a re-factorisation per step, which the API doesn't currently expose. Common parameter choices ------------------------ Two named members of the family cover almost every engineering case: * **Average acceleration** — :math:`\beta = 1/4, \gamma = 1/2`. Unconditionally stable for *any* :math:`\Delta t`, second-order accurate, no numerical damping. Trapezoidal rule from another angle. This is the **femorph-solver default** because the unconditional stability matches the user expectation that ``transient_solve(dt=...)`` doesn't silently blow up on a too-large step. * **Linear acceleration** — :math:`\beta = 1/6, \gamma = 1/2`. Conditionally stable (:math:`\Delta t < 2\sqrt{3} / \omega_\text{max}`), second-order accurate, slightly less numerical dispersion than average-acceleration on smooth problems, but the stability cliff makes it dangerous as a default. Two more from the same family that show up in special-purpose contexts: * **Central difference** — :math:`\beta = 0, \gamma = 1/2`. Explicit (:math:`\mathbf{M}` only on the LHS), conditionally stable. Used in explicit-dynamics codes (LS-DYNA, Abaqus/Explicit) but not in femorph-solver — the implicit path is the right tool for the structural-dynamics regime the project targets. * **Backward Euler** — :math:`\beta = 1, \gamma = 1`. Unconditionally stable, first-order accurate, heavy numerical damping. Used to "freeze" a transient at a quasi- static state for diagnostic purposes. The :meth:`Model.solve_transient` API exposes ``beta=0.25, gamma=0.5`` as the defaults; pass other values through the keyword arguments. Stability and accuracy regions ------------------------------ For the undamped scalar case (:math:`\mathbf{C} = 0`), the update operator's amplification factor at frequency :math:`\omega` is .. math:: A(\omega \Delta t) \;=\; \frac{1 - (\tfrac{1}{2} - \beta + \gamma) \Omega^2} {1 + \beta \Omega^2}, \qquad \Omega \equiv \omega \Delta t. Stability (no exponential growth) requires :math:`|A(\Omega)| \le 1` for every :math:`\Omega = \omega_i \Delta t` in the spectrum. The classical result (Hughes 2000, *The Finite Element Method: Linear Static and Dynamic FE Analysis*, §9): * Unconditionally stable iff :math:`2 \beta \ge \gamma \ge 1/2`. * Second-order accurate iff :math:`\gamma = 1/2` (otherwise first-order, with numerical damping at amplitude :math:`O(\Omega^2 \Delta t / 2)`). The "Newmark unconditional + 2nd order" sweet spot is therefore the locus :math:`\gamma = 1/2, 1/4 \le \beta \le 1/2`. Average-acceleration sits at one corner of that locus; linear-acceleration just outside it (so conditionally stable). Numerical damping at :math:`\gamma > 1/2` is the *per-frequency* damping that high-frequency content sees: .. math:: \zeta_\text{num}(\omega) \;\approx\; \tfrac{1}{2}\, \bigl(\gamma - \tfrac{1}{2}\bigr)\, \omega \Delta t. Choosing :math:`\gamma = 0.6` and a relatively coarse :math:`\Delta t` puts ~5 % numerical damping on every mode above :math:`\omega = 1/\Delta t` — handy for shock-response problems where the high-frequency content is parasitic, but disastrous when the high-frequency content is the answer. Time-step selection ------------------- Even an unconditionally stable scheme has an *accuracy* limit. Rule of thumb (Hughes 2000 §9.5): .. math:: \Delta t \;\lesssim\; \frac{T_\text{shortest}}{20}, where :math:`T_\text{shortest} = 2\pi/\omega_\text{max}` is the period of the highest mode the user wants to resolve. Sampling 20 points per period suppresses the 2nd-order phase error below ~1 %. The "highest mode the user wants to resolve" is *not* the mesh's highest natural frequency — that goes to MHz on a typical structural mesh and would force absurd time steps. It is the highest mode in the band of engineering interest, set by the load's frequency content and the analyst's judgment. Damping models -------------- femorph-solver supports two damping forms in :meth:`~femorph_solver.Model.solve_transient`: * **Pre-assembled damping matrix** — pass a sparse :math:`\mathbf{C}` directly. Lets the user assemble whatever damping model they trust (Rayleigh, structural, frequency- dependent at a fixed band centre, …). * **Rayleigh form** — pass ``rayleigh=(α, β)`` and the integrator builds :math:`\mathbf{C} = \alpha \mathbf{M} + \beta \mathbf{K}` inline using the same K and M it already needs to assemble for the implicit step. Avoids the second K + M assembly that pre-built Rayleigh would otherwise pay. The two forms are mutually exclusive — pass at most one. Modal damping (per-mode :math:`\zeta_i`) is *not* a transient- integrator option because it requires the eigendecomposition that the modal-superposition path uses; for problems that work with modal damping, prefer :meth:`~femorph_solver.Model.solve_harmonic` if the loading is pure-harmonic, or assemble a Rayleigh fit to the measured per- mode damping ratios for the transient. Roadmap: generalised-:math:`\alpha` ------------------------------------ Newmark-:math:`\beta` doesn't damp the highest-frequency end of the spectrum without also damping the engineering-relevant modes. Generalised-:math:`\alpha` (Chung & Hulbert 1993, *ASME J. Appl. Mech.* 60, 371–375) extends the family with two extra parameters that decouple high-frequency damping from low-frequency accuracy. It's the right tool for problems where the parasitic high-frequency content is significant; it isn't shipped today (see :doc:`/getting-started/roadmap`). References ---------- * **N. M. Newmark**, "A method of computation for structural dynamics," *ASCE Journal of the Engineering Mechanics Division* 85 (1959), 67–94 — the original Newmark family. * **T. J. R. Hughes**, *The Finite Element Method: Linear Static and Dynamic Finite Element Analysis*, Dover 2000, §9 — comprehensive treatment of stability / accuracy / numerical damping for the family. * **K.-J. Bathe**, *Finite Element Procedures*, 2nd ed., 2014, §9 — companion treatment, including the predictor / corrector form. * **R. D. Cook, M. E. Plesha, M. S. Witt**, *Concepts and Applications of Finite Element Analysis*, 4th ed., Wiley 2002, §11 — engineering-oriented treatment with worked examples. * **J. Chung & G. M. Hulbert**, "A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-α method," *ASME J. Appl. Mech.* 60 (1993), 371–375 — the follow-on family that's on the femorph-solver roadmap. See also -------- * :doc:`eigenproblem` — modal-superposition alternative for problems where the loading is harmonic or modal-damping is required. * :doc:`/user-guide/solving/transient` — driving API. * :doc:`/getting-started/known-limitations` — explicit dynamics, generalised-:math:`\alpha`, and non-linear transient are out of scope today.