Seismic Isolation Simulator — Theory of Computation

1. Overview of the Method

The Seismic Isolation Simulator models a base-isolated building either as a single-degree-of-freedom (SDOF) system with a rigid superstructure or as a shear-type multi-degree-of-freedom (MDOF) system with an elastic superstructure, and performs a nonlinear time-history response analysis that accounts for the nonlinear restoring force of the isolation layer (bilinear skeleton with the Masing rule). Time integration uses the direct-integration Newmark-β method (β = 1/4, average acceleration) or the Wilson-θ method (θ = 1.4) — choose the model and the method via the "⚙ Settings" button at the top right of the tool (default: Newmark-β). All computation runs in JavaScript inside your browser.

The analysis proceeds as follows.

2. Modeling Assumptions

3. Equation of Motion and Units

The equation of motion of the SDOF system subjected to ground acceleration ÿ_g (relative displacement x) is:

M·ẍ + Q(x, ẋ) = −M·ÿ_g

  M        : mass = ΣW / g
  ΣW       : total story weight [kN] (roof through 1st floor)
  g        : gravitational acceleration = 980.665 cm/s²
  Q(x, ẋ) : restoring force of the isolation layer [kN]
             (bilinear + Masing rule; depends on displacement and history)
  ÿ_g     : ground acceleration [cm/s²]

There is no viscous damping term C·ẋ (C = 0). Energy dissipation in the isolation layer is represented solely by the hysteresis loops of the restoring-force model (hysteretic damping).

4. Restoring Force of the Isolation Layer (Bilinear Skeleton Curve)

The relationship between the lateral force Q and the lateral displacement D of the isolation layer (laminated rubber bearings plus dampers, etc.) is modeled as a bilinear curve made of two straight branches. The three input values are:

The skeleton curve (Q–D relationship under initial loading) is:

if |D| ≤ up1:  Q = K1·D
if |D| > up1:  Q = sign(D) · { K1·up1 + K2·(|D| − up1) }
                 = sign(D)·Qd + K2·D    (second branch)

5. Hysteresis Loops by the Masing Rule

Unloading and reloading paths during the earthquake response follow the Masing rule: measured from each load-reversal point, the response traces the skeleton curve scaled by a factor of two in both axes. This produces spindle-shaped hysteresis loops whose enclosed area corresponds to the dissipated energy (hysteretic damping).

State variables:
  K      : depth of the hysteresis branch (1 = on the skeleton curve)
  D      : scaling factor (skeleton = 1, branch curves = 2)
  U0[·], V0[·] : displacement and force at each reversal point (stack)

Processing at each time step:
1. When the velocity changes sign (ẋ_m · ẋ_m−1 < 0):
     find the reversal point (U_TN, V_TN) by interpolation between steps
     and push it onto the stack (K++); from here on follow the
     twice-scaled curve (D = 2)
2. When the response passes the start of the current branch
   (an inner loop closes):
     return to the outer branch (K −= 2; if K ≤ 3, return to the
     skeleton curve, K = 1)
3. Evaluate the restoring force:
     Q = D · f( (x − U0[K−1]) / D ) + V0[K−1]
     f(·) : bilinear skeleton curve of Section 4

6. Time-History Integration (Newmark-β / Wilson-θ)

The equation of motion is written in incremental form and integrated step by step with the Newmark-β method (β = 1/4) or the Wilson-θ method (θ = 1.4), selectable in "⚙ Settings" (default: Newmark-β). The Wilson-θ method is an unconditionally stable implicit scheme that assumes the acceleration varies linearly over the extended interval t to t+θΔt (unconditionally stable for θ ≥ 1.37). The Newmark-β method with β = 1/4, γ = 1/2 (average-acceleration method) is an unconditionally stable implicit scheme with no numerical damping (Wilson-θ slightly damps high-frequency components numerically).

6.1 Integration constants

θ = 1.4,  Δt = time step of the ground-motion record

a0 = 6 / (θΔt)²      a1 = 3 / (θΔt)
a2 = 6 / (θΔt)       a3 = θΔt / 2
a4 = a0 / θ           a5 = −a2 / θ
a6 = 1 − 3/θ          a7 = Δt / 2         a8 = Δt² / 6

6.2 Computation at each time step

At each step m (time t_m):

1. Effective stiffness:        k̂ = k_m + a0·M
     k_m : instantaneous stiffness of the isolation layer (see 6.3)

2. Effective incremental load: ΔR̂ = M·( −θ·(ÿg_m − ÿg_m−1)
                                          + a2·ẋ_m−1 + 3·ẍ_m−1 )

3. Incremental displacement:   Δx = ΔR̂ / k̂
     (implemented as a Cholesky solve of the linear system;
      equivalent to a scalar division for this single-DOF tool)

4. Update the response:
     ẍ_m = a4·Δx + a5·ẋ_m−1 + a6·ẍ_m−1
     ẋ_m = ẋ_m−1 + a7·(ẍ_m + ẍ_m−1)
     x_m = x_m−1 + Δt·ẋ_m−1 + a8·(ẍ_m + 2·ẍ_m−1)

5. Restoring force:            Q_m from the Masing rule (Section 5)
                               using x_m and ẋ_m

6.3 Updating the instantaneous stiffness

The nonlinearity is tracked by using the secant stiffness computed from the displacement and force increments of the current step as the instantaneous stiffness of the next step (sequential stiffness update).

k_m+1 = |Q_m − Q_m−1| / |x_m − x_m−1|

  If the displacement increment is very small (|Δx| ≤ 10⁻⁶),
  the previous stiffness is kept.

6.4 The Newmark-β method (β = 1/4, average acceleration)

When Newmark-β is selected, the increment is solved directly over Δt without the extended θΔt interval. The effective stiffness, effective incremental load and update formulas are (γ = 1/2):

1. Effective stiffness:        k̂ = k_m + (4/Δt²)·M + (2/Δt)·C

2. Effective incremental load: ΔR̂ = M·( −(ÿg_m − ÿg_m−1)
                                          + (4/Δt)·ẋ_m−1 + 2·ẍ_m−1 )
                                   + C·( 2·ẋ_m−1 )
     (with β = 1/4 the acceleration coefficient of the C term is zero)

3. Incremental displacement:   Δx = ΔR̂ / k̂

4. Update the response:
     ẍ_m = (4/Δt²)·Δx − (4/Δt)·ẋ_m−1 − ẍ_m−1
     ẋ_m = ẋ_m−1 + (Δt/2)·(ẍ_m + ẍ_m−1)
     x_m = x_m−1 + Δx

6.5 Unbalanced-force carry-over (equilibrium correction)

In a plain incremental scheme, the unbalanced force generated in steps that cross a break point of the bilinear curve (the difference between the tangent-stiffness prediction and the actual restoring force) accumulates, and with coarse time steps can noticeably distort the response. This tool eliminates the accumulation by adding the residual of the equation of motion at the previous step to the effective incremental load of the current step.

res_m−1 = −M·(ẍ_m−1 + ÿg_m−1) − C·ẋ_m−1 − Q_m−1   (residual of the equation of motion)

ΔR̂ ← ΔR̂ + res_m−1

6.6 Extension to the MDOF model

For the shear-type MDOF model, M, C and K in the formulas above become matrices (N = stories + 1 DOF). The linear system of each step is solved by Cholesky decomposition. The superstructure story springs are elastic (Q = k·δ); only the isolation-layer spring is bilinear with the Masing rule.

M = diag(m_1, …, m_N)                mass matrix (m_i = W_i/g per floor)
K = shear-type tridiagonal matrix    assembled from story stiffness k_i and
                                     the instantaneous isolation stiffness
C = (2h/ω₁)·K_upper                  stiffness-proportional damping
                                     (isolation-layer rows are zero)

  h   : superstructure damping ratio (from Settings, default 2%)
  ω₁ : 1st circular frequency of the superstructure fixed at the 1st floor
  K_upper : stiffness matrix assembled from the story springs only

This C gives the fixed-base r-th mode a damping ratio ξ_r = h·ω_r/ω₁ (exactly h for the 1st mode). Because energy dissipation of the isolation layer is represented by hysteresis only, the isolation layer contributes nothing to C.

The eigenvalue analysis solves the mass-normalized standard problem A = M^−1/2·K·M^−1/2 with the Jacobi method to obtain natural periods and mode shapes. The displayed eigenvalue results are for the superstructure only, fixed at the 1F floor level; the isolation layer is excluded because its (bilinear) stiffness varies with deformation.

7. Input Ground Motions and Scaling

7.1 Built-in records

The following five acceleration records are built in (units: cm/s² = gal).

7.2 Scaling by peak ground velocity

In seismic isolation design it is customary in Japan to normalize input ground motions by their peak ground velocity (Vmax). The tool stores, for each record, the factor that scales the original wave to Vmax = 25 cm/s, and multiplies the whole acceleration history by a constant according to the selected Vmax (25, 50, 75 or 100 cm/s).

scale factor = (factor normalizing the original record to 25 cm/s)
               × (selected Vmax / 25)

Example: for EL CENTRO-NS the 25 cm/s factor is 0.747
         → selecting Vmax = 50 cm/s multiplies all data by 1.495
A user-specified integer Vmax is handled by the same formula.

7.3 Custom ground motions

8. Output of Results

9. Verification Example

A worked example using the sample "K-Kobe Residence". Selecting the same conditions in the tool reproduces these values.

10. Scope and Limitations

1. Rigid or elastic superstructure — the rigid-SDOF model cannot evaluate story drifts or higher-mode amplification. The shear-type MDOF model evaluates elastic story drifts and higher modes, but yielding of the superstructure is not considered.
2. One horizontal direction only — vertical motion, simultaneous two-directional input and torsion are not considered.
3. No viscous damping in the isolation layer — velocity-dependent devices such as oil dampers cannot be represented; isolation-layer energy dissipation comes from the bilinear hysteresis only. (Stiffness-proportional viscous damping can be applied to the MDOF superstructure.)
4. Constant isolator properties — dependence on axial pressure, velocity or temperature, uplift and hardening are not modeled.
5. No equilibrium iteration within a step — nonlinearity is tracked by sequential secant-stiffness updates with unbalanced-force carry-over (Section 6.5); small numerical errors can occur near the break points of the bilinear curve.
6. Educational use — this tool is intended for studying the behavior of base-isolated structures. It must not be used for actual design.

11. References

• Yorihiko Osaki, Building Vibration Theory (Kenchiku Shindo Riron), Shokokusha.

Related pages: Seismic Isolation Simulator  |  Steel Section Tables  |  Frame Analysis