A small Python library for the dynamic response of equivalent SDOF systems (beams, two-way slabs) under transient loading. Uses central-difference time integration and the Masing rule with memory for hysteretic restoring-force behaviour.
Required: Python 3, NumPy, Matplotlib. On macOS:
brew install python-matplotlibThat brings in Python 3.14 and NumPy as dependencies. No virtual environment or pip install needed; nothing else to configure. Jupyter is not required.
Verify everything works by running the test suite (see section 4).
sdof/ Library package
solver.py CentralDifferenceSolver
resistance.py LinearElastic, PiecewiseBackbone, Hysteretic
elements.py Element, Beam, TwoWaySlab, natural_period()
loadings.py Loading, TriangularPulse, ConstantLoad, PiecewiseLinearLoad
cases/ Example simulations — copy these to make new ones
beam_blast.py
slab_blast.py
tests/
test_hysteresis.py Constitutive test on sdof.Hysteretic
The library lives in sdof/. Simulations live in cases/. The
constitutive test lives in tests/.
python3 cases/beam_blast.pyFive plots open:
- Load vs time — F(t)
- Restoring force vs time — R(t)
- Acceleration vs time — y''(t)
- Displacement vs time — y(t)
- Restoring force vs displacement — R(y) hysteresis loop (primary diagnostic)
The R(y) plot is the most informative one. A correctly working
elastic-plastic system with reversals shows a hysteresis loop whose
unloading branches are parallel to the initial elastic slope. If the
loop does not close, or unloading slopes are wrong, the constitutive
code is suspect (run test_hysteresis.py to localise it).
To run headless (e.g., on a server, or for a smoke test):
MPLBACKEND=Agg python3 cases/beam_blast.pypython3 tests/test_hysteresis.pyPrints PASS (exit 0) or fails with a diff (exit 1).
test_hysteresis.py drives sdof.Hysteretic with prescribed
displacement histories (no integrator) and verifies the constitutive
rule: elastic unloading at slope k_e, yielding on the negative side
after a 2·yel excursion, monotonic backbone reproduction, partial
reload returns to previous peak, nested inner loops close without
drift, and Masing's scaled-backbone behaviour for large excursions.
If you change the hysteretic code, this test is the one to run first — if it fails, the bug is constitutive, not in the integrator.
Copy an existing case as a template:
cp cases/beam_blast.py cases/my_case.pyEdit the new file. The structure is always the same five steps.
A structural element is reduced to two things — an equivalent mass
M_e and a piecewise-linear (y, R) backbone. The backbone must start
at (0, 0), be strictly monotone in y, and capture the full virgin
loading curve (cracking, yield, hardening, softening, fracture — as
many segments as you like).
Beam from a Biggs-style table:
from sdof import Beam
beam = Beam.from_breakpoints(
M_e=22.5, # equivalent mass (kg)
breakpoints=[
(0.0, 0.0),
(0.0003581, 4375.9), # cracking
(0.0040219, 25615.9), # elastic yield
(0.0137303, 29236.4), # plastic peak
(0.0139645, 28082.9), # ultimate (slight softening)
],
)Bilinear shorthand for an elastic-perfectly-plastic beam:
beam = Beam.from_bilinear(M_e=22.5, yel=0.004, R_y=25_000.0)Two-way slab is the same API with TwoWaySlab.from_breakpoints(...):
from sdof import TwoWaySlab
slab = TwoWaySlab.from_breakpoints(
M_e=180.0,
breakpoints=[(0.0, 0.0), (0.002, 60_000.0), (0.020, 90_000.0)],
)The breakpoints are whatever your geometry / yield-line analysis produces. The library does not derive them — that's deliberately your input.
from sdof import TriangularPulse, ConstantLoad, PiecewiseLinearLoad
load = TriangularPulse(Fm=40_000.0, t0=0.006) # 40 kN peak, 6 ms
# load = ConstantLoad(F=10_000.0) # step
# load = PiecewiseLinearLoad([(0, 0), (0.01, 50_000.0), (0.02, 0)], outside=0.0)For arbitrary load shapes, define your own callable F(t) -> float and
pass it where a Loading is expected.
Pick the total simulation time T_total (long enough to see what you
care about — usually a few natural periods). The time step needs to be
much smaller than the natural period for the central-difference scheme
to be stable and accurate; the legacy convention is 1e-5 · T_n,
which is deeply conservative but matches existing results.
import numpy as np
from sdof import natural_period
T_total = 0.050
T_n = natural_period(beam.M_e, beam.backbone.k_initial)
dt = T_n * 1.0e-5
n = int(T_total / dt)
t = np.arange(0, (n + 1) * dt, dt)from sdof import CentralDifferenceSolver
solver = CentralDifferenceSolver(M=beam.M_e)
R = beam.hysteretic_resistance() # Masing rule + memory
# R = beam.elastic_resistance() # purely linear if you want
resp = solver.solve(t, F_func=load, R_func=R, y0=0.0, ydot0=0.0)resp is a small dataclass with five NumPy arrays of length n+1:
resp.t # time
resp.disp # y(t)
resp.acc # y''(t)
resp.R # restoring force
resp.load # F(t)import matplotlib.pyplot as plt
plt.figure(figsize=(12, 8))
plt.plot(resp.disp * 1000, resp.R / 1000)
plt.xlabel("Displacement (mm)")
plt.ylabel("Restoring force R (kN)")
plt.title("Hysteresis loop")
plt.grid()
plt.show()The R(y) hysteresis plot is the diagnostic for whether plastic behaviour and elastic unloading are correct. The R(t), y(t), and a(t) plots are the response time histories.
Run with:
python3 cases/my_case.pyHysteretic(backbone) implements the Masing rule with memory. On
monotonic loading from rest, R follows the multi-segment backbone
exactly. On a direction reversal at (y_rev, R_rev), the unload-reload
curve is R = R_rev ± 2 · backbone(|y − y_rev| / 2) — the backbone
scaled by 2 in both axes, anchored at the reversal. This guarantees
the unloading slope is the initial elastic stiffness k_e (because the
backbone's slope at zero is k_e), and that yielding on the opposite
side begins after a 2·yel elastic excursion. Each reversal pushes a
pair (y_rev, R_rev) onto a stack; when the active y crosses the
y-coordinate of the second-from-top entry, the inner loop has closed
and that pair is popped, returning the system to the enclosing curve.
This handles arbitrary nested cycling without drift.
For asymmetric materials (different yield strengths in tension and compression), pass a separate negative-side backbone:
from sdof import Hysteretic, PiecewiseBackbone
R = Hysteretic(
backbone=PiecewiseBackbone([(0.0, 0.0), (yel_pos, R_pos)]),
backbone_neg=PiecewiseBackbone([(0.0, 0.0), (yel_neg, R_neg)]),
)If backbone_neg is omitted, the negative envelope mirrors the
positive one.
The R(y) loop does not close. Run test_hysteresis.py. If it
passes, the constitutive code is fine — your time step may be too
large for the integrator to track the trajectory accurately, or your
load function is producing tiny direction flips that fragment the
reversal stack. Try a smaller dt.
The R(y) loop's unloading slope does not match the elastic slope.
This would be a constitutive bug. The hysteresis test asserts this
explicitly; if it passes there but you observe it in your simulation,
plot a few hundred steps near a reversal and compare the slope to
beam.backbone.k_initial.
Solution looks like noise / explodes. The central-difference
scheme is conditionally stable. Reduce reduction_factor (the multiplier
on T_n that gives dt); the legacy convention 1e-5 is well within
the stability limit and is a safe starting point.
Peak R seems too high or too low compared to the backbone. The Masing rule scales the backbone by 2 in both axes from the reversal point, so cyclic peaks can exceed the static backbone's peak by a modest amount under symmetric backbones with hardening. If a peak seems unphysical, plot R(y) and check whether you're on the virgin envelope (stack empty) or on a Masing curve from a deep reversal.