femlabpy.dynamics#

Time-history solvers, load builders, and response utilities for dynamics.

Workflow role#

This module owns the transient-analysis part of femlabpy. It contains the load callables that turn a spatial force pattern into a function of time, the time integrators that advance the structural state, and the plotting helpers used to inspect histories, energy balance, and frequency-response functions.

Public entry points#

NewmarkParams and TimeHistory hold solver settings and results.
constant_load, ramp_load, harmonic_load, pulse_load, tabulated_load, and seismic_load create time-dependent forcing functions.
solve_newmark, solve_hht, solve_central_diff, and solve_newmark_nl advance the dynamic equilibrium equations.
compute_frf, plot_time_history, plot_energy, and plot_frf support interpretation after the solve.

Functions#

`compute_frf`(M, C, K, input_dof, output_dof, ...)	Compute the Frequency Response Function (FRF), H(omega), for harmonic excitation.
`constant_load`(P)	Return a load function that is constant in time: p(t) = P.
`critical_timestep`(K, M, *[, method, n_iter])	Estimate the critical time step for conditionally stable schemes.
`harmonic_load`(P, omega[, phase])	Sinusoidal: p(t) = P * sin(omega * t + phase).
`plot_energy`(result, *[, ax])	Plot energy components vs time.
`plot_frf`(freq_hz, H, *[, ax, log_scale, ...])	Plot magnitude and phase of a Frequency Response Function.
`plot_time_history`(result, dof_index, *[, ...])	Plot displacement, velocity, or acceleration vs time at specified DOFs.
`pulse_load`(P, t_start, t_duration)	Rectangular pulse: P for t_start <= t <= t_start+t_duration, else 0.
`ramp_load`(P, t_ramp)	Linear ramp: p(t) = P * min(t / t_ramp, 1.0).
`seismic_load`(M, direction, accel_record, ...)	Construct a time-dependent effective force vector for uniform base excitation.
`solve_central_diff`(M_lumped, C, K, p_func, ...)	Explicit central difference time integration.
`solve_hht`(M, C, K, p_func, u0, v0, dt, nsteps, *)	HHT-alpha (Hilber-Hughes-Taylor) time integration.
`solve_newmark`(M, C, K, p_func, u0, v0, dt, ...)	Implicit Newmark-beta time integration for linear structural dynamics.
`solve_newmark_nl`(M, C, tangent_func, ...[, ...])	Nonlinear Newmark-beta with Newton-Raphson iteration at each step.
`tabulated_load`(P, time_table, value_table)	Interpolated load history: p(t) = P * interp(time_table, value_table)(t).

Function Reference#

femlabpy.dynamics.compute_frf(M, C, K, input_dof, output_dof, freq_range, *, n_points=500)[source]#

Compute the Frequency Response Function (FRF), H(omega), for harmonic excitation.

Parameters:

Mndarray or scipy.sparse matrix, shape (ndof, ndof): Global mass matrix.
Cndarray or scipy.sparse matrix, shape (ndof, ndof): Global damping matrix.
Kndarray or scipy.sparse matrix, shape (ndof, ndof): Global stiffness matrix.
input_dofint: The 0-based global degree-of-freedom index where the unit harmonic load P is applied.
output_dofint: The 0-based global degree-of-freedom index where the complex displacement response U is measured.
freq_rangetuple of floats: The frequency spectrum window (f_min, f_max) to evaluate, in Hertz (Hz).
n_pointsint, default 500: Number of discrete frequency sampling points within the freq_range.

Returns:

freq_hzndarray, shape (n_points,)

Array of frequencies evaluated (Hz).

Hndarray, shape (n_points,), dtype=complex

Array of complex scalar frequency response values H(omega). Use numpy.abs(H) for magnitude and numpy.angle(H) for phase.

Frequency values in Hz.

Hndarray, shape (n_points,), complex

Complex FRF values.

Examples

>>> freq, H = compute_frf(M, C, K, input_dof=0, output_dof=0,
...                       freq_range=(0.1, 50.0), n_points=1000)
>>> magnitude = np.abs(H)

femlabpy.dynamics.constant_load(P) → Callable[source]#: Return a load function that is constant in time: p(t) = P.

femlabpy.dynamics.critical_timestep(K, M, *, method: str = 'power', n_iter: int = 50) → float[source]#

Estimate the critical time step for conditionally stable schemes.

Parameters:

K, Mndarray or sparse: Stiffness and mass matrices.
methodstr: 'power' for power iteration to approximate omega_max.
n_iterint: Number of power iterations.

Returns:

dt_crfloat: Critical time step: dt_cr = 2 / omega_max.

femlabpy.dynamics.harmonic_load(P, omega: float, phase: float = 0.0) → Callable[source]#: Sinusoidal: p(t) = P * sin(omega * t + phase).

femlabpy.dynamics.plot_energy(result: TimeHistory, *, ax=None)[source]#

Plot energy components vs time.

Parameters:

resultTimeHistory: Must have compute_energy=True in the solver call.
axmatplotlib Axes, optional

Returns:

matplotlib Axes

femlabpy.dynamics.plot_frf(freq_hz, H, *, ax=None, log_scale=True, mark_peaks=True)[source]#

Plot magnitude and phase of a Frequency Response Function.

Parameters:

freq_hzndarray: Frequency values in Hz.
Hndarray, complex: FRF values.
axmatplotlib Axes or None: If None, creates a new 2-subplot figure.
log_scalebool: If True, use log scale for magnitude.
mark_peaksbool: If True, mark resonance peaks.

Returns:

figmatplotlib Figure

femlabpy.dynamics.plot_time_history(result: TimeHistory, dof_index, *, quantity='displacement', ax=None)[source]#

Plot displacement, velocity, or acceleration vs time at specified DOFs.

Parameters:

resultTimeHistory: Output from any time integrator.
dof_indexint or list of int: DOF index (0-based) or list of DOF indices.
quantitystr: 'displacement', 'velocity', or 'acceleration'.
axmatplotlib Axes, optional: If provided, plot on this axes.

Returns:

matplotlib Axes

femlabpy.dynamics.pulse_load(P, t_start: float, t_duration: float) → Callable[source]#: Rectangular pulse: P for t_start <= t <= t_start+t_duration, else 0.

femlabpy.dynamics.ramp_load(P, t_ramp: float) → Callable[source]#: Linear ramp: p(t) = P * min(t / t_ramp, 1.0).

femlabpy.dynamics.seismic_load(M, direction, accel_record, dt_record: float) → Callable[source]#

Construct a time-dependent effective force vector for uniform base excitation.

Parameters:

Mndarray or scipy.sparse matrix, shape (ndof, ndof): Global mass matrix of the structure.
directionarray_like, shape (ndof,): Unit influence vector r. This vector has a value of 1.0 for DOFs that act parallel to the direction of ground excitation, and 0.0 for orthogonal or unexcited DOFs.
accel_recordarray_like, shape (N,): Discrete time series of ground acceleration. Ensure the units match the mass matrix (e.g., m/s^2, not g’s, unless you pre-scale it).
dt_recordfloat: Time step (seconds) between samples in accel_record.

Returns:

Callable: A function p_eff(t) that returns the effective seismic load vector of shape (ndof, 1) at any time t.

femlabpy.dynamics.solve_central_diff(M_lumped, C, K, p_func: Callable, u0, v0, dt: float, nsteps: int, *, C_bc=None, dof: int = 2, compute_energy: bool = False) → TimeHistory[source]#

Explicit central difference time integration.

Requires a diagonal (lumped) mass matrix.

Parameters:

M_lumpedndarray, shape (ndof,) or (ndof, ndof) diagonal: Lumped (diagonal) mass matrix. Can be a 1D array of diagonal entries or a full diagonal matrix.
Cndarray or sparse, shape (ndof, ndof): Damping matrix. For pure explicit, use C=0.
Kndarray or sparse, shape (ndof, ndof): Stiffness matrix.
p_funccallable: Load function: p_func(t) -> ndarray (ndof, 1).
u0, v0array_like: Initial displacement and velocity.
dtfloat: Time step size (must be < dt_critical).
nstepsint: Number of time steps.
C_bcarray_like, optional: Boundary constraint table.
dofint, default 2: DOFs per node.
compute_energybool: If True, compute energy at each step.

Returns:

TimeHistory: Full time history of u, v, a.

Raises:

ValueError: If the mass matrix appears non-diagonal.

femlabpy.dynamics.solve_hht(M, C, K, p_func: Callable, u0, v0, dt: float, nsteps: int, *, alpha: float = -0.05, C_bc=None, dof: int = 2, compute_energy: bool = False) → TimeHistory[source]#

HHT-alpha (Hilber-Hughes-Taylor) time integration.

A generalization of Newmark-beta that provides controllable numerical dissipation of high-frequency modes while maintaining second-order accuracy.

Parameters:

M, C, Kndarray or sparse: Mass, damping, stiffness matrices.
p_funccallable: Load function.
u0, v0array_like: Initial conditions.
dtfloat: Time step.
nstepsint: Number of time steps.
alphafloat, default -0.05: HHT parameter in [-1/3, 0]. alpha=0 recovers standard Newmark. More negative = more high-frequency dissipation.
C_bc, dof, compute_energysee solve_newmark

Returns:

TimeHistory

femlabpy.dynamics.solve_newmark(M, C, K, p_func: Callable, u0, v0, dt: float, nsteps: int, *, beta: float = 0.25, gamma: float = 0.5, C_bc=None, dof: int = 2, compute_energy: bool = False) → TimeHistory[source]#

Implicit Newmark-beta time integration for linear structural dynamics.

Solves the transient equations of motion:: M * a(t) + C * v(t) + K * u(t) = p(t)

Parameters:

Mndarray or scipy.sparse matrix, shape (ndof, ndof): Global mass matrix.
Cndarray or scipy.sparse matrix, shape (ndof, ndof): Global damping matrix.
Kndarray or scipy.sparse matrix, shape (ndof, ndof): Global stiffness matrix.
p_funccallable: Load function returning the force vector at time t. Signature: p_func(t) -> ndarray of shape (ndof, 1).
u0array_like, shape (ndof, 1): Initial displacement vector.
v0array_like, shape (ndof, 1): Initial velocity vector.
dtfloat: Time step size (seconds).
nstepsint: Number of time steps to compute. Total time = dt * nsteps.
betafloat, default 0.25: Newmark beta parameter. Controls variation of acceleration over the step.
gammafloat, default 0.5: Newmark gamma parameter. Controls artificial numerical damping.
C_bcarray_like, optional: Boundary constraint table. Constrained DOFs are perfectly fixed (u=v=a=0) during the dynamic response.
dofint, default 2: DOFs per node (used to interpret C_bc).
compute_energybool, default False: If True, compute kinetic and strain energy at each step.

Returns:

TimeHistory: Dataclass containing the full time history of displacements u, velocities v, accelerations a, and optionally energy.

Notes

Default parameters (beta=0.25, gamma=0.5) use the constant average acceleration method (trapezoidal rule). This scheme is unconditionally stable for linear systems and produces no artificial numerical dissipation.

femlabpy.dynamics.solve_newmark_nl(M, C, tangent_func: Callable, internal_force_func: Callable, p_func: Callable, u0, v0, dt: float, nsteps: int, *, beta: float = 0.25, gamma: float = 0.5, tol: float = 1e-06, max_iter: int = 20, C_bc=None, dof: int = 2) → TimeHistory[source]#

Nonlinear Newmark-beta with Newton-Raphson iteration at each step.

Parameters:

M, Cndarray or sparse: Mass and damping matrices.
tangent_funccallable: tangent_func(u, state) -> K_t (tangent stiffness).
internal_force_funccallable: internal_force_func(u, state) -> (q, state_new).
p_funccallable: External load function: p_func(t) -> p.
u0, v0array_like: Initial conditions.
dt, nstepsfloat, int: Time step and number of steps.
beta, gammafloat: Newmark parameters.
tolfloat: Newton convergence tolerance.
max_iterint: Max Newton iterations per step.
C_bcarray_like, optional: Boundary constraint table.
dofint: DOFs per node.

Returns:

TimeHistory

femlabpy.dynamics.tabulated_load(P, time_table, value_table) → Callable[source]#

Interpolated load history: p(t) = P * interp(time_table, value_table)(t).

Parameters:

Parray_like: Spatial load pattern (ndof, 1).
time_tablearray_like: Time sample points.
value_tablearray_like: Scalar multiplier values at each time sample.

Classes#

`NewmarkParams`([beta, gamma])	Newmark-beta integration parameters.
`TimeHistory`(t, u, v, a, dt, nsteps[, energy])	Container for time integration results.

Class Reference#

class femlabpy.dynamics.NewmarkParams(beta: float = 0.25, gamma: float = 0.5)[source]#

Bases: object

Newmark-beta integration parameters.

Attributes:

betafloat: Newmark beta parameter (default 0.25, average acceleration).
gammafloat: Newmark gamma parameter (default 0.5).

classmethod average_acceleration() → NewmarkParams[source]#: Unconditionally stable, no numerical dissipation, 2nd order.

classmethod linear_acceleration() → NewmarkParams[source]#: Conditionally stable (dt < 0.551*T_min), 2nd order.

classmethod central_difference() → NewmarkParams[source]#: Explicit, conditionally stable (dt < T_min/pi), 2nd order.

classmethod fox_goodwin() → NewmarkParams[source]#: Conditionally stable, 4th order for undamped SDOF.

class femlabpy.dynamics.TimeHistory(t: ndarray, u: ndarray, v: ndarray, a: ndarray, dt: float, nsteps: int, energy: dict | None = None)[source]#

Bases: object

Container for time integration results.

Attributes:

tndarray, shape (nsteps+1,): Time values at each step.
undarray, shape (nsteps+1, ndof): Displacement history.
vndarray, shape (nsteps+1, ndof): Velocity history.
andarray, shape (nsteps+1, ndof): Acceleration history.
dtfloat: Time step size.
nstepsint: Number of time steps.
energydict or None: Energy history: kinetic, strain, dissipated, external, total.