阻尼为系统的固有特征，对结构的动力反应有显著影响。根据材料的耗能特性，可建立相应的阻尼模型^[1]。对于土木工程而言，粘滞阻尼和滞后阻尼是两个广泛使用的模型。粘滞阻尼力与速度成正比，具有计算简便的特点而成为广泛使用的模型^[2]，如Rayleigh阻尼模型和速度相关型的阻尼器^[3-5]。但是，粘滞阻尼的耗能和频率相关。已有实验表明，很多材料的耗能与频率无关，如型钢混凝土梁^[6]、土^[7-8]和粘弹性夹层梁和板^[9-10]等，此时采用阻尼力与位移成正比相位差π/2的滞后阻尼模型更符合实验结果。滞后阻尼将导致复刚度运动方程，而常采用频域方法进行求解。频域计算方法以Fourier变换为基础，理论上仅适用于线弹性体系，对于非线性系统，常采用等效线性化进行计算^[11-13]。为进行真非线性计算，需要在时域中进行直接积分法进行计算。

为在时域中进行复本构运动方程的求解，朱镜清^[14-15]、何钟怡^[16]等根据对偶原则，建立与实荷载对应的对偶虚荷载，完善了滞后阻尼体系的输入理论。但是，采用时域直接积分法计算滞后阻尼体系动力反应，即使是无条件稳定的直接积分方法也易于出现发散的结果^[17-19]。微分方程数值解不稳定的原因有两种：一种是数值方法的不稳定；另一种是由于方程本身具有发散解。复阻尼运动方程的逐步积分法的发散解是由后一种原因造成的^{[17, 20]}。为使复阻尼运动方程的逐步积分法计算结果稳定，一种常用的方法是将滞后阻尼等效为近似的粘滞阻尼^[21-22]，但计算结果的误差较大。

在直接积分法得到滞后阻尼体系稳定解方面，孙攀旭等分别基于激励插值方法^[23]和滞变阻尼时域理论^[24]，建立了稳定的直接积分计算方法。事实上，滞后阻尼体系的特征值为互为相反数的复数特征对^[25]，必然有一个特征值的实部是正的，由此导致强迫荷载下补解中的一项是没有物理意义的发散项，在解析解时人工删除发散项而使计算结果稳定。而直接积分法的计算结果是包含补解的，且滞后阻尼体系直接积分法是收敛到包含发散项补解^[17]而导致计算结果不稳定。Pan等^[26]首先提出虚初始条件的概念，然后基于无条件稳定的Newmark法建立使滞后阻尼直接积分解不出现发散项的计算方法。本文进一步建立了恒载下的虚初始条件，以及基于有条件稳定的中心差分法建立滞后阻尼的中心差分虚初始条件法。在此基础上，通过算例验证验证所提方法的稳定性、计算精度和计算效率。

1   滞后阻尼体系的虚初始条件

1.1   自由振动

单自由度滞后阻尼体系的运动方程为^[25]：

方程的初始条件为：

式中：m和k分别为系统的质量和刚度；η为复阻尼系数；u₀和v₀分别为初始位移和初始速度，是实数。式(1)的不稳定的通解为：

稳定的通解为^[25]：

式中：$s \!=\! \omega (\alpha - {\rm{i}}\mu )$，$\omega\! =\!\! \sqrt {k/m}$，$\alpha \!=\! \sqrt {\dfrac{{ \!-\! 1\! +\!\! \sqrt {1\! +\!\! {\eta ^2}} }}{2}}$，$\mu \!= \!\sqrt {\dfrac{{1 + \sqrt {1 + {\eta ^2}} }}{2}}$，${C_1} \!=\! \dfrac{{s{u_0} + {v_0}}}{{2s}}$，${C_2}\! =\! \dfrac{{s{u_0} - {v_0}}}{{2s}}$，$C =$${u_0} - {\rm{i}}\dfrac{{{v_0} + \omega \alpha {u_0}}}{{\omega \mu }}$。由于s的实部大于零，式(3)中右边第一项${C_1}{{\rm{e}}^{st}}$将产生无物理意义的发散解。在解析解的求解过程中，可以人为去除发散项而得到稳定解，但直接积分法将收敛到式(3)而产生发散解。因此，为使直接积分法稳定的关键是使式(3)和式(4)相同。为此，引入虚初始条件$u_0^{v1}$和$v_0^{v{\rm{1}}}$达到这个目的^[26]。在引入虚初始条件后，总初始条件为：

以式(5)为初始条件，并使式(3)等于式(4)，可得：

求解式(6)和式(7)的代数方程组可得：

式(8)和式(9)表明，虚初始条件$u_0^{v1}$和$v_0^{v{\rm{1}}}$为纯虚数。因此，对于总初始条件的实部可解释为初始条件中可观察和测量的部分，与此同时还存在与实初始条件相对应的虚部。在这样的复初始条件下，式(1)的通解为式(4)的稳定解而无需人为干涉。

1.2   简谐荷载作用

在简谐荷载$A{{\rm{e}}^{{\rm{i}}\theta t}}$作用下，单自由度体系的运动方程为：

式中：A为简谐荷载的幅值；θ为迫振频率。方程(10)不稳定的通解为：

稳定的通解为：

式中：$X = \dfrac{A}{{ - {\theta ^2}m + (1 + {\rm{i}}\eta )k}}$为稳态解的幅值；D₁、D₂和C_f为常数，与初始条件有关。当$u(0) = {\rm{0}}$和$\dot u(0)\! \!=\!\! {\rm{0}}$时，${D_1}\!\! = \!\!\dfrac{{ \!-\! s\! -\! {\rm{i}}\theta }}{{2s}}X$，${D_2} \!\!=\! \!\dfrac{{{\rm{i}}\theta X \!\!-\!\! sX}}{{2s}}$，${C_{\rm{f}}}\! =\!- {\rm{Re}} (X)-$${\rm{i}}\dfrac{{\theta {\rm{Im}} (X) - \alpha \omega {\rm{Re}} (X)}}{{\mu \omega }}$，Re表示取实部，Im表示取虚步。直接积分法解收敛到式(11)而导致计算结果发散。为使式(11)等于式(12)，引入虚初始位移$u_0^{v2}$和虚初始速度$v_0^{v2}$。并令与虚初始条件$u_0^{v2}$和$v_0^{v2}$对应的D₁=0和D₂=C_f，即：

由此可得：

$u_0^{v2}$和$v_0^{v2}$也是纯虚数。式(12)中的${C_{\rm{f}}}{{\rm{e}}^{ - st}}$是由荷载引起的伴生自由振动，与初始条件无关。虚初始条件$u_0^{v2}$和$v_0^{v2}$也是与荷载相关的纯虚数，不改变实部可测量部分，但可使收敛到式(11)的直接积分法稳定。

式(15)和式(16)是关于零初始条件下的虚初始条件。当$u(0) = {u_0}$, $\dot u(0) = {v_0}$时，由叠加原理可得简谐荷载$A{{\rm{e}}^{{\rm{i}}\theta t}}$下稳定的通解为：

而不稳定的通解为：

同理引入虚初始位移$u_0^v$和虚初始速度$v_0^v$。由式(8)、式(9)、式(15)和式(16)可知，当$u_0^v = u_0^{v1} + u_0^{v{\rm{2}}}$，$v_0^v = v_0^{v1} + v_0^{v{\rm{2}}}$时，式(17)和式(18)相同，从而使收敛到式(18)的直接积分法解也稳定。即非零初始条件下，简谐荷载作用的总初始条件为：

1.3   恒载作用

在恒载${A_0}(1 + {\rm{i}}\eta )/2$作用下，单自由度体系的运动方程为：

式(20)不稳定的通解为：

当$u(0) = {\rm{0}}$，$\dot u(0) = {\rm{0}}$时，可得E₁和E₂为：

零初始条件下稳定的通解为：

式中，${C_{\rm{f}{0}}} = - {A_0}{\rm{/}}2k + {\rm{i}}\alpha {A_0}{\rm{/}}2\mu k$。为使式(21)等于式(23)，引入虚初始位移$u_0^{v3}$和虚初始速度$v_0^{v3}$。在虚初始条件$u_0^{v3}$和$v_0^{v3}$下，可得E₁和E₂的解为：

令式(24)中的${E_1} = 0$，${E_2} = {C_{{\rm{f}}0}}$可得：

如果将恒载看成激振频率为零的简谐荷载，且稳态解的幅值$X = {A_0}/2k$，则式(25)和式(26)可以直接由式(15)和式(16)计算得到。因此，非零初始条件$u(0) = {u_0}$, $\dot u(0) = {v_0}$下，恒载作用的总初始条件的计算方法与式(19)相同。

1.4   地震作用

若已知地面运动加速度为a(t)，将a(t)采用Fourier级数展开：

式中，${\theta _j} = 2\pi j{\rm{/}}N\Delta {t_{\rm{s}}}$(j=0,1,2,···, N/2), $\Delta {t_{\rm{s}}}$是加速度时程的采样时间间隔，${t_l} = l\Delta {t_{\rm{s}}}$(l=0, 1 , 2,···, N−1)。根据对偶荷载原则，可构造a(t)的对偶函数b(t)^[15]：

则地震作用下单自由度体系的运动方程为：

式中：$f(t) = - m[a(t) + {\rm{i}}b(t)]$。对于$u(0) = {\rm{0}}$和$\dot u(0) = {\rm{0}}$，对f(t)中每一项简谐荷载对应的虚初始条件求和，可得地震作用的虚初始条件为：

式中，${C_{{\rm{f}}j}} = - {\rm{Re}} ({X_j}) -{\rm{i}}\dfrac{{{\theta _j}{\rm{Im}} ({X_j}) - \alpha \omega {\rm{Re}} ({X_j})}}{{\mu \omega }}$(j=0, 1, 2,···, N/2)，${X_0} = \dfrac{{m{A_0}}}{{2k}}$，${X_j} = \dfrac{{m({A_j} - {\rm{i}}{B_j})}}{{ - \theta _j^2m + (1 + {\rm{i}}\eta )k}}$(j=1, 2, ···, N/2−1)，${X_{N/2}} = \dfrac{{m{A_{N/2}}}}{{2( - \theta _{N/2}^2m + (1 + {\rm{i}}\eta )k)}}$。

因此，非零初始条件$u(0) = {u_0}$, $\dot u(0) = {v_0}$下，由叠加原理可得体系的总初始条件为：

2   虚初始条件的中心差分计算方法

若已知t_n时刻及以前时刻的反应已知的情况下，根据中心差分法，求解t_n+1时刻的位移方程为：

式中：$\Delta t$为计算时间步长；$\hat k = m{\rm{/}}\Delta {t^2}$。虚初始条件计算方法，仅需根据实初始条件计算相应的虚部作为初始条件，其他计算过程与粘滞阻尼体系的中心差分法^[3]相同。同时，为避免数值积分过程中舍入误差的影响，对每个荷载步计算完成后，由平衡条件进行速度和加速度的计算，可提高计算精度。即由式(35)得到t_n+1时刻的位移后，根据t_n+1时刻的平衡条件可得${\ddot u_{n + 1}}$为：

为计算${\dot u_{n + 1}}$，可根据t_n+1时刻的差分条件：

将式(38)、式(39)消去${u_{n + {\rm{2}}}}$，可得：

由此可得t_n+1时刻的位移虚部修正解：

则t_n+1时刻的解为：

为此，中心差分法具体的计算步骤如表1所示。

表  1  任意荷载下中心差分虚初始条件法计算步骤

Table  1.  Calculation procedure of central differential virtual initial condition method under arbitrary load

1) 准备工作
1.1 对加速度时程a(t)根据式(28)和式(29)计算系数
1.2 计算虚初始条件，形成总初始条件$u_0^t = u_0^{} + u_0^{v1} + u_0^{v{\rm{4}}}$, $v_0^t = v_0^{} + v_0^{v1} + v_0^{v{\rm{4}}}$
1.3 $\ddot u_0^t = \left[ { {f_0} - k(1 + {\rm{i}}\eta )u_0^t} \right]/m$，$u_{ - 1}^t = u_0^t - \Delta t\dot u_0^t - {(\Delta t)^2}\ddot u_0^t{\rm{/2}}$
1.4 选择Δt
1.5 $\hat k = m{\rm{/}}\Delta {t^2}$
2) 每个荷载步计算
2.1 ${\hat f_{n + 1} } = {f_n} - \left[ {k(1 + {\rm{i}}\eta ) - \dfrac{ {2m} }{ {\Delta {t^2} } } } \right]u_n^t - \dfrac{m}{ {\Delta {t^2} } }u_{n - 1}^t$
2.2 ${u_{n + 1}} = {\hat f_{n + 1}}/\hat k$
2.3 ${\ddot u_{n + 1} } = \dfrac{1}{m}[ { {f_{n + 1} } - k(1 + {\rm{i}}\eta ){u_{n + 1} } } ]$，${\dot u_{n{\rm{ + 1} } } } = \dfrac{ {\Delta t} }{2}{\ddot u_{n + {\rm{1} } } }{\rm{ + } }\dfrac{ { {u_{n + {\rm{1} } } } - {u_n} } }{ {\Delta t} }$
2.4 $u_{n{\rm{ + } }1}^{v1} = - {\rm{i}}\dfrac{ { {\rm{Re} } ({ {\dot u}_{n + 1} }) + \omega \alpha {\rm{Re} } ({u_{n + 1} })} }{ {\omega \mu } }$
2.5 $u_{n + 1}^t = {\rm{Re}} ({u_{n{\rm{ + 1}}}}) + u_{n + 1}^{v1} + u_0^{v{\rm{4}}}$
3) 用n+1代替n，重复2.1~2.5进行下一荷载步的计算。

下载: 导出CSV 

| 显示表格

3   算例

为验证本文算法的有效性，下面以文献[17]分析过的算例进行本文算法的验证。已知体系的滞后阻尼系数η=0.1，对无阻尼自振频率f=0.1 Hz、1 Hz、10 Hz的三个体系进行地震反应计算。以表2中的3条地震波分别作为体系的地震输入。输入地震波的加速度时程如图1所示。在进行中心差分法计算时，根据中心差分法的稳定性条件和输入地震波的采样时间间隔，计算时间步长$\Delta t$取为：

式中：T_n为体系的自振周期；Δt_s为地震波加速度时程的采样时间间隔。

表  2  地震波主要参数

Table  2.  The main parameters of ground motions

地震波	时间	地震	峰值加速度/(m/s2)
El Centro	1940.5.18	California	3.417
迁安波	1976.8.31	唐山余震	0.974
天津波	1976.11.25	唐山余震	1.458

下载: 导出CSV 

| 显示表格

图  1  输入地震波的加速度时程

Figure  1.  Acceleration time histories of input seismic waves

下载: 全尺寸图片 幻灯片

作为对比，对地震作用形成式(31)的方程后，计算各项简谐荷载的解析解并求和，所得结果作为解析解。当$u(0) = {u_0}$, $\dot u(0) = {v_0}$时，位移、速度和加速度的解析解为：

对式(31)两边进行Fourier变换^[27]，可得：

式中，$F({\theta _0}) = m{A_0}/2$, $F({\theta _j}) = m({A_j} - {\rm{i}}{B_j})$(j=1, 2, ···, N/2−1)，$F({\theta _{N/2}}) = m{A_{N/2}}/2$，则：

利用式(46)进行Fourier逆变换，可得位移、速度和加速度的频域解为：

对比式(44)和式(47)可知，解析解包括由初始条件引起的自由振动，由荷载引起的伴生自由振动和荷载所引起的纯强迫振动。而频域解仅包含荷载所引起的纯强迫振动。

当u₀=v₀=0，f=10 Hz时，三条地震波作用下体系的位移反应如图2所示，f=10 Hz时天津波作用下的速度和加速度反应如图3所示，f=1 Hz和0.1 Hz时天津波作用下的位移反应如图4所示。

图  2  三条地震作用下位移反应(f=10 Hz, u₀=v₀=0)

Figure  2.  Displacement responses under three seismic excitations (f=10 Hz, u₀=v₀=0)

下载: 全尺寸图片 幻灯片

图  3  天津波作用下的速度和加速度反应(f=10 Hz, u₀=v₀=0)

Figure  3.  Velocity and acceleration response under Tianjin wave excitation

下载: 全尺寸图片 幻灯片

图  4  天津波作用下不同自振频率体系的位移反应(u₀=v₀=0)

Figure  4.  Displacement response for various natural vibration frequencies under Tianjin wave excitation

下载: 全尺寸图片 幻灯片

由图2~图4可知，对于不同的地震输入和自振频率体系，中心差分虚初始条件法计算的位移、速度和加速度都和解析解基本一致，计算结果稳定而不存在临界自振周期的问题^[17]。而对于频域解，当f=10 Hz时与解析解几乎重合；当f=0.1 Hz时，频域解与解析解存在明显差别。这是因为频域解得到的是体系的稳态反应，而解析解是包含稳态和伴生自由振动的瞬态反应，自振频率越高，瞬态反应衰减越快，因此，对于自振频率较高的体系，瞬态反应对体系的总反应影响很小而可直接采用稳态解进行描述，但是，对于自振频率低的体系，瞬态反应衰减需要较长的时间，当自振周期和地震波持时为同一量级时，则瞬态反应在整个地震反应过程中都有显著影响，此时采用仅包含稳态解的频域分析方法将产生显著误差。本文所提的中心差分虚初始条件法，计算结果包含瞬态反应和稳态反应，对于不同自振周期体系的反应都和解析解吻合的很好。

当体系初始条件非零时，在天津波作用下不同自振频率的位移反应如图5所示，f=10 Hz体系的速度与加速度反应如图6所示。初始条件非零情况与零初始条件下的计算结果对比可知，由于初始位移和初始速度瞬态振动的进一步影响，导致频域解在初始阶段误差增大。对于自振频率较高的10 Hz体系，由于瞬态振动很快消失，此时忽略初始条件对体系总反应的影响较小。但是对于自振周期与输入地震波持时为同一量级的0.1 Hz体系，初始位移和初始速度引起的瞬态振动将对总反应产生显著影响，此时忽略初始条件将引起显著误差。

图  5  天津波作用下不同自振频率体系的位移反应

Figure  5.  Displacement response of various natural vibration frequencies under Tianjin wave excitation

下载: 全尺寸图片 幻灯片

图  6  天津波作用下体系的速度和加速度反应(f=10 Hz)

Figure  6.  Velocity and acceleration response under Tianjin wave excitation

下载: 全尺寸图片 幻灯片

为定量研究不同方法的精确性，不同方法相对解析解的峰值相对误差为：

式中：${| {r_{}^*(t)} |_{\max }}$表示精确解的峰值；${\left| {r(t)} \right|_{\max }}$表示近似解的峰值。

零初始条件下本文方法及频域解的峰值相对误差如表3所示。由表3可知，当体系的自振频率较高时，频域解的计算精度高；而当体系的自振频率较低时，频域解将产生显著误差；而中心差分虚初始条件法对于不同自振周期体系反应的峰值相对误差都小于5%，显示了良好的精度。

表4为u₀=v₀=0，f=1 Hz时，天津波作用下三种方法位移反应的计算时间。解析解为输入地震波的Fourier变换时间和式(44a)各项直接求和的计算时间，中心差分虚初始条件法的计算时间包括输入地震波的Fourier变换时间、计算虚初始条件的时间和每个时刻位移、速度和加速度的递推计算的时间，频域法则为Fourier变换和Fourier逆变换的总时间。在进行解析解计算时，由于自由振动和伴生自由振动的影响，无法采用快速Fourier变换进行计算而采用逐项求和，因此，计算时间最长。频域采用快速Fourier变换，计算时间最短。中心差分虚初始条件法通过前一步的计算结果递推后一步的计算结果，计算效率仅次于频域方法，且计算时间显著小于解析解。结合表3和表4的数据分析结果表明，中心差分虚初始条件法兼顾了计算精度和计算效率。

表  3  本文方法与频域解的峰值相对误差

Table  3.  Peak relative errors of the proposed method and frequency domain solution

地震波	自振频率/Hz	${e_u}$/(%)			${e_{\dot u}}$/(%)			${e_{\ddot u}}$(%)
		本文方法	频域解		本文方法	频域解		本文方法	频域解
El Centro波	10.0	1.46	0.0002		0.210	0.000		1.37	0.0000
	1.0	0.64	0.0029		0.540	0.020		0.08	0.0218
	0.1	0.69	17.1300		0.360	37.960		0.39	1.7400
迁安波	10.0	0.08	0.0000		0.004	0.000		0.27	0.0000
	1.0	1.23	0.0104		1.540	0.006		0.11	0.0033
	0.1	3.36	18.0800		0.930	2.570		0.12	0.1000
天津波	10.0	0.76	0.0002		1.930	0.000		2.21	0.0000
	1.0	0.19	0.4200		0.770	0.420		0.06	0.2500
	0.1	0.53	6.0800		0.320	0.170		0.18	0.0700

下载: 导出CSV 

| 显示表格

表  4  不同算法计算时间比较

Table  4.  Comparison of calculation time of various methods /s

地震波	解析解	本文方法	频域解
天津波	1.40	0.003	0.0015

下载: 导出CSV 

| 显示表格

4   结论

针对滞后阻尼体系直接积分法收敛到包含发散项而不稳定的问题，提出了中心差分虚初始条件计算方法。由理论分析和数值计算结果可知：

(1)对于不同峰值加速度、卓越频率的地震输入，及不同自振频率体系，本文所提的中心差分虚初始条件法对位移、速度和加速度均能得到稳定的计算结果，不存在临界自振周期的问题。实初始条件是可观测和测量部分，虚初始条件是伴随实初始条件而存在的，可无需人为干涉而使滞后阻尼体系计算结果稳定。由此进行直接积分法计算，即使有条件稳定的中心差分法，依然可以得到稳定的计算结果。

(2)对于自振频率较高的体系，瞬态反应对体系的总反应影响很小而可直接采用稳态解进行描述，而对于自振频率低的体系，如自振周期和地震波持时为同一量级时，瞬态反应在整个地震反应过程中都有显著影响，此时采用仅包含稳态解的频域分析方法将产生显著误差。

(3)中心差分虚初始条件法的计算结果包括瞬态反应和稳态反应，计算误差和体系的自振频率无关，可适用于不同的自振频率体系。

(4)中心差分虚初始条件法计算的峰值相对误差小于5%。同时，计算时间显著小于解析解，因此，这种方法兼顾了计算精度和计算效率。