今天整理电脑里的资料，翻到了二十多年前写的一个 C++ 程序（我二十多年前居然还能写这个），是用 Nigam-Jennings 方法求解精确的地震波反应谱。这个算法我最初是在李宏男的《结构多维抗震理论与设计方法》^[1] 中看到的。当时在写这个程序的时候还专门作了个帮助文件来记录这个算法，可惜这个帮助文件是 hlp 格式，Windows 10 以上已经不支持了，也就没法打开查看内容，于是我重新搜集了一下资料，复习了一下这个算法的相关知识，并把它记录下来。

基本理论

Nigam 和 Jennings 提出的数值积分算法直接从单自由度体系运动微分方程出发：

\ddot{y}(t) + 2 \xi \omega \dot{y}(t) + \omega^2 y(t) =-\ddot{X}(t)

这里， $\ddot{X}(t)$ 是以数字形式给出的，而且数字化密度一般很高，有理由假定相邻两个数字加速度间加速度是按线性变化，于是在 $t_i\leqslant t \leqslant t_{i+1}$ 区间内， $\ddot{X}(t)$ 可写成：

\ddot{X}(t) = \ddot{X}(t_i) + \frac{\ddot{X}(t_{i+1}) - \ddot{X}(t_i) }{t_{i+1}-t_i}(t-t_i) = \ddot{X}(t_i) + \frac{\Delta \ddot{X}(t_i)}{\Delta t_i} (t-t_i)

单自由度体系的运动微分方程可以写成为：

\ddot{y}(t) + 2 \xi \omega \dot{y}(t) + \omega^2 y(t) = -\ddot{X}(t_i) - \frac{\Delta \ddot{X}(t_i)}{\Delta t_i} (t-t_i)

在 $[t_i,t_{i+1}]$ 区间内，上式的解为：

\begin{aligned} y(t) = {} & \mathrm{e}^{-\xi\omega(t-t_i)} \Bigl[ c_1 \cos \omega'(t-t_i) + c_2 \sin\omega'(t-t_i)\Bigr]\\ & {}- \frac{\ddot{X}(t_i)}{\omega^2}+\frac{2\xi}{\omega^3}\cdot\frac{\Delta\ddot{X}(t_i)}{\Delta t_i}-\frac{1}{\omega^2}\frac{\Delta\ddot{X}(t_i)}{\Delta t_i}(t-t_i) \end{aligned}

式中：

\left\{\begin{aligned} \omega' & = \sqrt{1-\xi^2} \omega \\ c_1 & = y(t_i) - \frac{2\xi}{\omega^3}\cdot\frac{\Delta \ddot{X}(t_i)}{\Delta t_i} + \frac{\ddot{X}(t_i)}{\omega^2}\\ c_2 & = \frac{1}{\omega'}\left(\xi \omega y(t_i) + \dot{y}(t_i) - \frac{2\xi^2-1}{\omega^2} \cdot \frac{\Delta\ddot{X}(t_i)}{\Delta t_i} + \frac{\xi}{\omega} \ddot{X}(t_i)\right) \end{aligned}\right.

根据微分方程的解可以得到体系在 $t_n$ 和 $t_{n+1}$ 时刻的运动状态之间的递推公式为：

\left\{\begin{aligned} y(t_{i+1}) & = a_{11} y(t_i) + a_{12} \dot{y}(t_i) +b_{11} \ddot{X}(t_i) + b_{12} \ddot{X}(t_{i+1}) \\ \dot{y}(t_{i+1}) & = a_{21} y(t_i) + a_{22} \dot{y}(t_i) +b_{21} \ddot{X}(t_i) + b_{22} \ddot{X}(t_{i+1}) \end{aligned}\right.

当数字记录 $\ddot{X}(t_i)$ 的时间间隔一定时，式中参数为：

\begin{aligned} a_{11}= {} & \mathrm{e}^{-\xi\omega\Delta t} \left(\frac{\xi}{\sqrt{1-\xi^2}}\sin \omega'\Delta t +\cos \omega' \Delta t \right)\\ a_{12}= {} & \mathrm{e}^{-\xi\omega\Delta t} \cdot \frac{\sin \omega'\Delta t}{\omega'} \\ a_{21}= {} & - \frac{\omega}{\sqrt{1-\xi^2}} \cdot \mathrm{e}^{-\xi\omega\Delta t} \cdot \sin \omega'\Delta t \\ a_{22}= {} & \mathrm{e}^{-\xi\omega\Delta t} \left(\cos \omega'\Delta t - \frac{\xi}{\sqrt{1-\xi^2}} \sin\omega'\Delta t\right)\\ b_{11}= {} & \mathrm{e}^{-\xi\omega\Delta t} \left[\left( \frac{2\xi^2-1}{\omega^2\Delta t} + \frac{\xi}{\omega} \right)\frac{\sin \omega'\Delta t}{\omega'}+\left( \frac{2\xi}{\omega^3\Delta t} + \frac{1}{\omega^2} \right) \cos \omega'\Delta t\right]-\frac{2\xi}{\omega^3\Delta t}\\ b_{12}= {} & -\mathrm{e}^{-\xi\omega\Delta t} \left( \frac{2\xi^2-1}{\omega^2\Delta t} \cdot \frac{\sin \omega' \Delta t}{\omega'} +\frac{2\xi}{\omega^3\Delta t} \cos \omega'\Delta t \right) - \frac{1}{\omega^2} + \frac{2\xi}{\omega^3\Delta t} \\ b_{21}= {} & \mathrm{e}^{-\xi\omega\Delta t} \Biggl[ \left( \frac{2\xi^2-1}{\omega^2 \Delta t}+\frac{\xi}{\omega} \right) \left( \cos \omega' \Delta t -\frac{\xi}{\sqrt{1-\xi^2}} \sin\omega'\Delta t \right) \\ & - \left( \frac{2\xi}{\omega^3\Delta t} + \frac{1}{\omega^2} \right) \left( \omega'\sin\omega'\Delta t+ \xi \omega \cos \omega' \Delta t\right) \Biggr]+ \frac{1}{\omega^2\Delta t}\\ b_{22}= {} & \mathrm{e}^{-\xi\omega\Delta t} \left[\frac{1}{\omega^2\Delta t} \cos \omega'\Delta t + \frac{\xi}{\omega^2 \sqrt{1-\xi^2}\Delta t} \sin \omega' \Delta t\right] - \frac{1}{\omega^2\Delta t} \end{aligned}

而体系在 $t=t_{i+1}$ 时刻的绝对加速度可直接由运动微分方程得到：

\ddot{X}(t_{i+1})+\ddot{y}(t_{i+1}) = -2\xi\omega y(t_{i+1})-\omega^2y(t_{i+1})

根据递推公式就可以得到整个时间历程的位移、速度、和加速度时程反应，对不同周期的体系分别进行计算并取最大值即为地震波的位移、速度和加速度反应谱。

本方法以运动方程的精确解为基础，为单自由度体系在时间步长的起终点时刻的运动状态建立矩阵形式的递推关系。如果没有舍入误差，该方法给出的是精确解，且不会因为截断误差而发散。

代码实现

二十多年前我是用 C++ 实现的这个算法，现在来看，那个代码写得就是一堆 shit mountain，不堪入目，于是我用现在更熟悉一些的 AutoLisp/VisualLisp 来实现，代码如下：

;;; 单自由度体系时程计算
(defun CAL:SdofTh (AccLst Dt Period xi / a11 a12 a21 a22 b11 b12 b21 b22 w wd wdt w2 w3 ee xi2 xi3 xi4 tha thv thd aa vv dd)
  (setq xi2 (* xi xi)
        xi3 (sqrt (- 1.0 xi2))
        xi4 (/ xi xi3)
        w   (/ pi 0.5 Period)
        wd  (* xi3 w)
        wdt (* wd Dt)
        w2  (* w w)
        w3  (* w w w)
        ee  (exp (* -1.0 xi w dt))
        thd (list 0.0)
        thv (list 0.0)
        tha (list 0.0)
        Lst Acclst)
  (setq a11 (* ee (+ (* xi4 (sin wdt)) (cos wdt)))
        a12 (* ee (/ (sin wdt) wd))
        a21 (* ee -1.0 (/ w xi3) (sin wdt))
        a22 (* ee (- (cos wdt) (* xi4 (sin wdt))))
        b11 (- (* ee (+ (* (+ (/ (- (* 2.0 xi2) 1.0) w2 dt) (/ xi w)) (/ (sin wdt) wd)) (* (+ (/ xi 0.5 w3 Dt) (/ 1.0 w2)) (cos wdt)) )) (/ xi 0.5 w3 Dt))
        b12 (- (/ xi 0.5 w3 Dt) (* ee (+ (/ (* (- (* 2.0 xi2) 1.0) (sin wdt)) w2 wd Dt) (* (/ xi 0.5 w3 Dt) (cos wdt)))) (/ 1.0 w2))
        b21 (+ (* ee (- (* (+ (/ (- (* 2.0 xi2) 1.0) w2 Dt) (/ xi w)) (- (cos wdt) (* xi4 (sin wdt)))) (* (+ (/ xi 0.5 w3 Dt) (/ 1.0 w2)) (+ (* wd (sin wdt)) (* xi w (cos wdt)))))) (/ 1.0 w2 Dt))
        b22 (- (* ee (+ (/ (cos wdt) w2 Dt) (* (/ xi4 w2 Dt) (sin wdt)))) (/ 1.0 w2 Dt)))
  (while (> (length Lst) 1)
    (setq dd  (+ (* (car thd) a11) (* (car thv) a12) (* (car lst) b11) (* (cadr lst) b12))
          vv  (+ (* (car thd) a21) (* (car thv) a22) (* (car lst) b21) (* (cadr lst) b22))
          aa  (+ (* -2.0 xi w vv) (* -1.0 w2 dd))
          tha (cons aa tha)
          thv (cons vv thv)
          thd (cons dd thd)
          Lst (cdr Lst)))
  (mapcar 'reverse (list tha thv thd)))
;;; 地震动反应谱计算
(defun CAL:Spectra (Acclst Dt Tlst xi / RLst)
  (foreach Period Tlst
    (setq Rlst 
      (cons 
				(cons Period 
					(mapcar 
						'(lambda(x)
							(apply 'max x)) 
						(CAL:SdofTh Acclst Dt Period xi)))
				RLst)))
  (reverse RLst))

参考文献

[1] 李宏男. 结构多维抗震理论与设计方法 [M]. 科学出版社, 北京, 1998.

极客斋

求解地震波反应谱的 Nigam-Jennings 法

基本理论

代码实现

参考文献