Transient analysis#

femorph_solver.solvers.transient.solve_transient() integrates the structural equation of motion

\[M \, \ddot{u} + C \, \dot{u} + K \, u = F(t)\]

over a uniform timestep grid using the Newmark-β scheme. With the default parameters \(\beta = 1/4,\ \gamma = 1/2\) the integrator reduces to the unconditionally-stable average-acceleration method; any consistent pair (β, γ) can be used.

Algorithm#

Newmark-β updates:

\[\begin{split}u_{n+1} &= u_n + \Delta t\, \dot{u}_n + \Delta t^{2}\bigl[\tfrac{1}{2}(1-2\beta)\ddot{u}_n + \beta\,\ddot{u}_{n+1}\bigr] \\ \dot{u}_{n+1} &= \dot{u}_n + \Delta t\bigl[(1-\gamma)\ddot{u}_n + \gamma\,\ddot{u}_{n+1}\bigr]\end{split}\]

Rearranging the equation of motion at \(t_{n+1}\) gives an effective linear system

\[A \, u_{n+1} = b_{n+1}\]

with \(A = a_0 M + a_1 C + K\) and \(a_0 = 1/(\beta \Delta t^2)\), \(a_1 = \gamma / (\beta \Delta t)\). \(A\) is factored once outside the loop; every subsequent time step is a cheap back-solve.

API#

from femorph_solver.solvers.transient import solve_transient

K = m.stiffness_matrix()
M = m.mass_matrix()

# Constant step load (no damping):
result = solve_transient(
    K, M, F=load_vector,
    dt=1e-4,
    n_steps=1000,
    prescribed={d: 0.0 for d in fixed_dofs},
)

result.time             # (n_steps + 1,) time grid
result.displacement     # (n_steps + 1, N)
result.velocity         # (n_steps + 1, N)
result.acceleration     # (n_steps + 1, N)

Time-varying loads pass a callable instead of an array:

def load(t):
    return force_template * np.sin(2 * np.pi * 50 * t)

result = solve_transient(K, M, F=load, dt=1e-4, n_steps=2000,
                         prescribed=prescribed)

Damping#

Pass a sparse damping matrix C if you have one:

C = 0.01 * M + 0.001 * K       # Rayleigh damping
result = solve_transient(K, M, F=load, dt=1e-4, n_steps=1000, C=C,
                         prescribed=prescribed)

Without a C, the integrator runs the undamped form. femorph-solver does not (yet) synthesise damping from material records; build the matrix yourself with whatever model is appropriate.

Initial conditions#

u0 and v0 set the starting displacement / velocity. Omitting them zero-initialises both. The initial acceleration is derived from the equation of motion at \(t=0\) so M ü₀ = F(0) - K u₀ - C v₀ is satisfied exactly.

Explicit central-difference#

Setting beta=0, gamma=0.5 would give the explicit central-difference method, but the current implementation factors \(A\) with scipy.sparse.linalg.splu and thus requires \(\beta > 0\). Explicit paths need a separate mass-factoring code path not yet implemented.

Thread control#

thread_limit wraps both the one-off \(A\) factorisation and every back-solve. See Thread control.

Stability notes#

Unconditional stability only applies to \(\gamma \geq 1/2\) and \(\beta \geq (1/4)(1/2 + \gamma)^2\). The default \((1/4, 1/2)\) sits exactly on that bound.
Numerical damping appears if \(\gamma > 1/2\). The trapezoid variant \((1/4, 1/2)\) is energy-conservative on linear problems.