使用高阶无条件稳定显式算法的实时混合模拟稳定性

摘要

为了提高实时混合模拟（RTHS）计算效率、计算精度和稳定性，基于具有二阶精度的无条件稳定显式Chen‒Ricles（CR）和Chang积分算法，提出具有高阶精度的无条件稳定新型双显式（NDE）和新型半显式（NSE）积分算法。基于离散控制理论，推导出使用两种高阶算法的RTHS系统的闭环离散传递函数，分析了时滞、阻尼比、刚度比例系数、积分算法等对RTHS系统稳定性的影响，并以时滞微分方程推导出的近似解和精确解作为参照。分析结果表明，时滞会降低RTHS系统的稳定性；增大刚度比例系数会降低不同分析结果间的差异；增大阻尼比会提升RTHS系统的稳定性，同时会增大考虑积分算法与不考虑积分算法的稳定性差异，也会增大使用不同积分算法的稳定性差异；使用两种高阶算法的RTHS系统稳定性优于使用两种二阶算法的RTHS系统稳定性。

关键词

实时混合模拟; 积分算法; 稳定性; 离散控制理论

实时混合模拟（RTHS）是一种先进的抗震试验技术^［1-3］。RTHS系统包括试验子结构、数值子结构、积分算法及加载设备。为满足RTHS的实时性，除了对硬件设备 要求高，对积分算法的计算效率也有很高的要求。为提高积分算法的计算效率，近些年来学者们提出一系列“基于模型的积分算法”^［4］，这类算法的位移表达式为显式，且无条件稳定，求解位移相关的非线性动力学问题时无需迭代；当结构自由度数较多时，无需采用很小的时间步长来满足积分算法的稳定。因此，这类算法具有较高的计算效率，已在RTHS^［5-7］中得到了成功的应用。Kolay等^［6］进行了设置非线性黏滞阻尼器的两层混凝土框架的RTHS，非线性黏滞阻尼器和钢筋混凝土框架分别作为试验子结构和数值子结构，采用基于力的纤维单元来模拟框架的梁和柱，使用基于模型的积分算法实现数值子结构的实时分析和运算。

基于模型的积分算法可以进一步分为半显式和双显式算法：半显式算法仅位移为显式，而速度为隐式，而双显式算法的速度和位移均为显式。由Chen和Ricles提出的CR算法^［8］和由Chang提出的Chang算法^［9］分别为半显式和双显式算法的典型例子。近些年，学者们提出了一些新型的半显式^［10］和双显式的基于模型的积分算法^［11-12］，这些算法均具有优异的数值特性。积分算法除了位移表达式的显隐式和稳定性外，数值精度也是一个非常重要的特性。数值精度由阶数来量化，其定义为数值解相对于精确解的收敛速率。与一阶和二阶算法相比，高阶精度算法具有如下优势：能更加准确地求解一些复杂动力学问题，如非线性弹簧摆问题；更适合求解强非线性的长时间动力学问题；使用更大的时间步长可以达到低阶精度算法采用较小时间步长时的精度，从而提升了计算效率。上述提到的基于模型的积分算法^［8-12］的精度均不超过二阶，这无疑限制了算法性能的进一步提升。而传统的高阶精度积分算法（以龙格库塔法为例）的构造相对复杂，且无法同时满足无条件稳定和位移显式，应用于多自由度非线性动力学问题时并不具有显著优势。

为此，本文作者近期提出了一种新型半显式（NSE）^［13］和一种新型双显式（NDE）^［14］的高阶无条件稳定显式算法。这两种新算法的精度可达四阶，且具有极小的周期延长。与已有的基于模型的积分算法相比，具有更高的计算效率和精度。傅博等^［13-14］通过框架有限元模拟地震响应分析、非线性弹簧摆等算例验证了两种新算法在计算效率和精度上的优势，同时也说明两种新算法用于RTHS具有较大的潜力，可以满足RTHS的实时性要求。但是，RTHS除了要满足实时性，同时也要满足稳定性。RTHS的稳定性是保证RTHS成功实施的重要前提，如果RTHS系统发生失稳，将会影响试验结果的准确性，甚至会损坏试验设备。因此，为了探究两种新算法是否能用于RTHS，需要进行使用两种算法的RTHS系统的稳定性分析。已有研究通常采用放大矩阵法^［15-16］和离散控制理论^［17-22］对RTHS系统进行稳定性分析。这两种分析方法本质上是等效的，但是后者在考虑积分算法和时滞影响时更加有效和方便。因此，本文基于离散控制理论，研究使用两种新算法的RTHS系统的稳定性，为推广新算法在RTHS的应用提供理论依据。

1 高阶无条件稳定显式算法

NSE和NDE算法分别采用Chang和CR算法的速度和位移表达式，对于单自由度（SDOF）体系，Chang和NSE算法的速度位移表达式如下：

$${\dot{x}}_{i+\mathrm{1}}={\dot{x}}_i+\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }t{\ddot{x}}_i+\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }t{\ddot{x}}_{i+\mathrm{1}}$$

（1）

$$x_{i+\mathrm{1}}=x_i+{\alpha }_{\mathrm{1}}\mathrm{\Delta }t{\dot{x}}_i+{\alpha }_{\mathrm{2}}\mathrm{\Delta }{t}^{\mathrm{2}}{\ddot{x}}_i$$

（2）

式中：$\dot{x}$、$\ddot{x}$和$x$分别为速度、加速度和位移；$\mathrm{\Delta }t$为时间步长；下标i+1和i分别代表下一时刻和当前时刻；${\alpha }_{\mathrm{1}}$和${\alpha }_{\mathrm{2}}$为积分参数。

由式（1）、（2）可以看出，Chang和NSE算法的速度表达式为隐式，而其位移表达式为显式，根据定义，这两种算法为半显式算法。

Chang算法的积分参数如下：

$${\alpha }_{\mathrm{1}}=\frac{\mathrm{4}\xi \Omega +\mathrm{4}}{{\Omega }^{\mathrm{2}}+\mathrm{4}\xi \Omega +\mathrm{4}}$$

（3）

$${\alpha }_{\mathrm{2}}=\frac{\mathrm{2}}{{\Omega }^{\mathrm{2}}+\mathrm{4}\xi \Omega +\mathrm{4}}$$

（4）

式中：$\xi$为结构阻尼比；$\Omega ={\omega }_n\mathrm{\Delta }t$，${\omega }_n$为结构圆频率。

NSE算法的积分参数推导过程如下：①将算法的速度位移表达式及SDOF体系的运动方程进行z变换，联立求得算法的离散传递函数（输出和输入分别为位移和外激励的z变换）及特征方程。②采用2阶Pade近似将连续域的极点映射到离散域的极点。③将算法离散传递函数的分母系数和离散域的极点带入特征方程，即可联立求得积分参数。具体推导过程及公式详见文献［13］。NSE算法的积分参数如下：

$${\alpha }_{\mathrm{1}}=\frac{\mathrm{144}\xi \Omega +\mathrm{144}}{{\Omega }^{\mathrm{4}}+\mathrm{12}\xi {\Omega }^{\mathrm{3}}+(\mathrm{48}{\xi }^{\mathrm{2}}+\mathrm{12}){\Omega }^{\mathrm{2}}+\mathrm{144}\xi \Omega +\mathrm{144}}$$

（5）

$${\alpha }_{\mathrm{2}}=\frac{-\mathrm{2}\xi {\Omega }^{\mathrm{3}}+(\mathrm{72}-\mathrm{96}{\xi }^{\mathrm{2}})\xi \Omega +\mathrm{72}}{{\Omega }^{\mathrm{4}}+\mathrm{12}\xi {\Omega }^{\mathrm{3}}+(\mathrm{48}{\xi }^{\mathrm{2}}+\mathrm{12}){\Omega }^{\mathrm{2}}+\mathrm{144}\xi \Omega +\mathrm{144}}$$

（6）

CR和NDE算法的速度位移表达式如下：

$${\dot{x}}_{i+\mathrm{1}}={\dot{x}}_i+{\alpha }_{\mathrm{1}}\mathrm{\Delta }t{\ddot{x}}_i$$

（7）

$$x_{i+\mathrm{1}}=x_i+\mathrm{\Delta }t{\dot{x}}_i+{\alpha }_{\mathrm{2}}\mathrm{\Delta }{t}^{\mathrm{2}}{\ddot{x}}_i$$

（8）

由式（7）、（8）可以看出，CR和NDE算法的速度表达式为隐式，而其位移表达式为显式。根据定义，这两种算法为双显式算法。CR算法的积分参数如下：

$${\alpha }_{\mathrm{1}}={\alpha }_{\mathrm{2}}=\frac{\mathrm{4}}{{\Omega }^{\mathrm{2}}+\mathrm{4}\xi \Omega +\mathrm{4}}$$

（9）

NDE算法的积分参数推导过程与NSE算法相同，具体推导过程及公式详见文献［14］。NDE算法的积分参数如下：

$${\alpha }_{\mathrm{1}}=\frac{\mathrm{144}}{{\Omega }^{\mathrm{4}}+\mathrm{12}\xi {\Omega }^{\mathrm{3}}+(\mathrm{48}{\xi }^{\mathrm{2}}+\mathrm{12}){\Omega }^{\mathrm{2}}+\mathrm{144}\xi \Omega +\mathrm{144}}$$

（10）

$${\alpha }_{\mathrm{2}}=\frac{\mathrm{24}\xi \Omega +\mathrm{144}}{{\Omega }^{\mathrm{4}}+\mathrm{12}\xi {\Omega }^{\mathrm{3}}+(\mathrm{48}{\xi }^{\mathrm{2}}+\mathrm{12}){\Omega }^{\mathrm{2}}+\mathrm{144}\xi \Omega +\mathrm{144}}$$

（11）

综上，NDE和NSE算法的区别主要在于以下两点：①二者速度位移表达式不同，前者速度、位移均为显式，而后者仅位移为显式。应用于位移相关的非线性问题时，两种算法均无需迭代；然而在应用于速度相关的非线性问题时，后者将退化为隐式算法，需要迭代。②二者的积分参数不同，前者的积分参数见式（10）、（11），后者的积分参数见式（5）、（6）。

图1为NDE和NSE算法的谱半径。当阻尼比为0时，谱半径为1；当阻尼比为正值时，谱半径小于1，说明两个算法均为无条件稳定。

图1  积分算法的谱半径

Fig. 1  Spectra radius of integration algorithms

图2为NDE和NSE算法的收敛速率。与CR和Chang算法对比，NDE和NSE算法的斜率为4，说明NDE和NSE算法为四阶精度，而CR和Chang算法只有二阶精度。

图2  积分算法的收敛速率

Fig. 2  Convergence rates of integration algorithms

2 基于离散控制理论的RTHS稳定性分析方法

对于SDOF结构的RTHS系统，其运动方程可表达为

$$m{\ddot{x}}_{i+\mathrm{1}}+c{\dot{x}}_{i+\mathrm{1}}+r_{\mathrm{N},i+\mathrm{1}}+r_{\mathrm{E},i+\mathrm{1}}=F_{i+\mathrm{1}}$$

（12）

式中：m和c分别为结构的质量和阻尼；F为外力；r为恢复力；下标N和E分别代表数值子结构和试验子结构；当结构处于线弹性时，数值子结构和试验子结构的恢复力分别为$k_{\mathrm{N}}{x}_{i+\mathrm{1}}$和$k_{\mathrm{E}}{x}_{i+\mathrm{1}}$，$k_{\mathrm{N}}$和$k_{\mathrm{E}}$分别为数值子结构和试验子结构的刚度，$k_{\mathrm{N}}+k_{\mathrm{E}}$为结构的总刚度k，定义刚度系数$\eta =k_{\mathrm{E}}/k$。

因为只有数值子结构参与计算，式（12）可改写为

$$m{\ddot{x}}_{i+\mathrm{1}}+c{\dot{x}}_{i+\mathrm{1}}+r_{\mathrm{N},i+\mathrm{1}}=F_{i+\mathrm{1}}-r_{\mathrm{E},i+\mathrm{1}}=P_{i+\mathrm{1}}$$

（13）

作动器常用线性斜坡生成器生成命令，较小的时间间隔$\mathrm{\Delta }t$内假定作动器响应是线性的。由于作动器时滞的存在，试验子结构的恢复力需要考虑时滞的影响。作动器需要$\alpha \mathrm{\Delta }t$才能达到位移指令，则时滞$\tau =\left(\alpha -\mathrm{1}\right)\mathrm{\Delta }t$，$\alpha$为大于等于1的值，当$\alpha$=1时，表示无时滞。

当考虑时滞时，作动器在i+1时刻能达到的实际位移$x_{i+\mathrm{1}}^{\mathrm{\text{'}}}$为

$$x_{i+\mathrm{1}}^{\mathrm{\text{'}}}=x_i^{\mathrm{\text{'}}}+\left(x_{i+\mathrm{1}}-x_i\right)/\alpha$$

（14）

对应地，试验子结构的恢复力$r_{\mathrm{E},i+\mathrm{1}}^{\mathrm{\text{'}}}$为

$$r_{\mathrm{E},i+\mathrm{1}}^{\mathrm{\text{'}}}=k_{\mathrm{E}}{x}_{i+\mathrm{1}}^{\mathrm{\text{'}}}$$

（15）

综上，考虑时滞的RTHS系统的方框图如图3所示。

图3  RTHS系统的方框图

Fig. 3  Block diagram of RTHS system

图3中，F（z）、P（z）、X（z）和$R_{\mathrm{E}}^{\mathrm{\text{'}}}\left(z\right)$分别由F_i+1、P_i+1、x_i+1和$r_{\mathrm{E},i+\mathrm{1}}^{\mathrm{\text{'}}}$经过z变换得到；G（z）=X（z）/P（z）是反映积分算法的传递函数，可以通过对式（13）和积分算法的位移、速度表达式进行z变换，然后联立求得；H（z）可通过对式（14）、（15）进行z变换，然后联立求得，其表达式为

$$H(z)=\frac{R_{\mathrm{E}}^{\mathrm{\text{'}}}\left(z\right)}{X\left(z\right)}=\frac{k_{\mathrm{E}}z}{\alpha z-\alpha +\mathrm{1}}$$

（16）

最终，RTHS系统的闭环离散传递函数可以表达为

$$G_{\mathrm{c}}(z)=\frac{X\left(z\right)}{F\left(z\right)}=\frac{G\left(z\right)}{\mathrm{1}+G\left(z\right)H\left(z\right)}=\frac{n_{\mathrm{2}}{z}^{\mathrm{2}}+n_{\mathrm{1}}z+n_{\mathrm{0}}}{d_{\mathrm{3}}{z}^{\mathrm{3}}+d_{\mathrm{2}}{z}^{\mathrm{2}}+d_{\mathrm{1}}z+d_{\mathrm{0}}}$$

（17）

式中：$n_i$和$d_i$分别为分子和分母的多项式系数。

CR和NDE算法的系数见表1；Chang和NSE算法的系数见表2。

表1  CR和NDE算法的闭环离散传递函数系数

Tab. 1  Coefficients of closed-loop discrete transfer function of CR and NDE algorithms

系数	计算式
$n_2$	$\alpha {\alpha }_2\mathrm{\Delta }{t}^2$
$n_1$	$\left[\left(1-\alpha \right){\alpha }_2+\alpha \left({\alpha }_1-{\alpha }_2\right)\right]\mathrm{\Delta }{t}^2$
$n_0$	$\left(1-\alpha \right)\left({\alpha }_1-{\alpha }_2\right)\mathrm{\Delta }{t}^2$
$d_3$	$\alpha m$
$d_2$	${\alpha }_2\left(\alpha -\alpha \eta +\eta \right)k\mathrm{\Delta }{t}^2+\alpha {\alpha }_{1}c\mathrm{\Delta }t+\left(1-3\alpha \right)m$
$d_1$	$\left[\left(\alpha {\alpha }_1-2\alpha {\alpha }_2+{\alpha }_2\right)\left(1-\eta \right)+\left({\alpha }_1-{\alpha }_2\right)\eta \right]k\mathrm{\Delta }{t}^2+\left(1-2\alpha \right){\alpha }_{1}c\mathrm{\Delta }t+\left(3\alpha -2\right)m$
$d_0$	$\left(1-\alpha \right)\left[\left({\alpha }_1-{\alpha }_2\right)\left(1-\eta \right)k\mathrm{\Delta }{t}^2-{\alpha }_{1}c\mathrm{\Delta }t+m\right]$

表2  Chang和NSE算法的闭环离散传递函数系数

Tab. 2  Coefficients of closed-loop discrete transfer function of Chang and NSE algorithms

系数	计算式
$n_2$	$\alpha \left({\alpha }_1+2{\alpha }_2\right)\mathrm{\Delta }{t}^2$
$n_1$	$\left[\left(1-\alpha \right)\left({\alpha }_1+2{\alpha }_2\right)+\alpha \left({\alpha }_1-2{\alpha }_2\right)\right]\mathrm{\Delta }{t}^2$
$n_0$	$\left(1-\alpha \right)\left({\alpha }_1-2{\alpha }_2\right)\mathrm{\Delta }{t}^2$
$d_3$	$\alpha \left(2m+c\mathrm{\Delta }t\right)$
$d_2$	$\left({\alpha }_1+2{\alpha }_2\right)\left(\alpha -\alpha \eta +\eta \right)k\mathrm{\Delta }{t}^2+\left(1-\alpha \right)c\mathrm{\Delta }t+2\left(1-3\alpha \right)m$
$d_1$	$\left[\left(-4\alpha {\alpha }_2+{\alpha }_1+2{\alpha }_2\right)\left(1-\eta \right)+\left({\alpha }_1-2{\alpha }_2\right)\eta \right]k\mathrm{\Delta }{t}^2-\alpha c\mathrm{\Delta }t+2\left(3\alpha -2\right)m$
$d_0$	$\left(1-\alpha \right)\left[\left({\alpha }_1-2{\alpha }_2\right)\left(1-\eta \right)k\mathrm{\Delta }{t}^2-c\mathrm{\Delta }t+2m\right]$

根据闭环离散传递函数的表达式，可以得到不同$\alpha$对应的稳定限值${\Omega }_{\mathrm{m}\mathrm{a}\mathrm{x}}$变化规律。

3 稳定性分析结果

Wallace等^［23］采用时滞微分方程推导出具有常时滞RTHS系统的稳定条件，包括近似解和精确解。需要注意的是，这两个解均没有考虑积分算法对RTHS系统稳定性的影响。

近似解的稳定限值如下：

$${\Omega }_{\mathrm{m}\mathrm{a}\mathrm{x}}=\frac{\mathrm{2}\xi }{\left(\alpha -\mathrm{1}\right)\eta }$$

（18）

精确解的稳定限值如下：

$$\widehat{\tau }=\frac{\mathrm{1}}{\widehat{\omega }}\left[\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{{\widehat{\omega }}^{\mathrm{2}}-\mathrm{1}}{\eta }+\mathrm{1}\right)+\mathrm{2}n\mathrm{\pi }\right]$$

（19）

$$\xi =\frac{\sqrt[]{{\eta }^{\mathrm{2}}-{\left({\widehat{\omega }}^{\mathrm{2}}-\mathrm{1}+\eta \right)}^{\mathrm{2}}}}{\mathrm{2}\widehat{\omega }}$$

（20）

$${\Omega }_{\mathrm{m}\mathrm{a}\mathrm{x}}=\frac{{\widehat{\tau }}_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\alpha -\mathrm{1}}$$

（21）

式（19）~（21）中：$\widehat{\omega }=\omega /{\omega }_n$；$\widehat{\tau }={\omega }_n\tau$；n为非负整数。

考虑两种阻尼比ξ=2%和5%及4种刚度比例系数η=0.25、0.50、0.75和1.00，基于离散控制理论，对使用CR、Chang、NDE和NSE 4种算法的RTHS系统进行稳定性分析，以Wallace的近似解、精确解作为对比。

图4为不同阻尼比和刚度比例系数时的稳定限值的对比。通过图4的对比分析，可以得出以下结论：

图4  稳定限值对比

Fig. 4  Comparison of stability limit

（1）时滞相当于负阻尼，对RTHS的稳定性有不利的影响。随着时滞的增大，图4所有工况下的稳定限值均会下降。

（2）Wallace的近似解和精确解没有考虑积分算法的影响，与考虑积分算法的结果可能存在较大差异（图4e）。因此，在进行RTHS稳定性分析时，需要考虑积分算法的影响。

（3）当阻尼比ξ相同，刚度比例系数η不同时，如图4a~4d所示，可以发现：随着刚度比例系数η的增大，试验子结构的刚度占比增大，意味着数值积分算法的影响变小，因此考虑积分算法与不考虑积分算法的稳定限值差异会减小。

（4）当阻尼比ξ不同，刚度比例系数η相同时，如图4a和4e所示，可以发现：结构阻尼对提高RTHS的稳定性有利。随着阻尼比ξ的增大，稳定限值增大，考虑积分算法与不考虑积分算法的稳定限值差异增大，不同积分算法的稳定限值差异增大。

（5）图4的所有工况均表明：CR和Chang算法稳定性接近，Chang算法的稳定限值略高于CR算法；NDE和NSE算法的稳定限值接近，NSE算法稳定性略高于NDE算法。

（6）图4的所有工况均表明：高阶无条件稳定显式NDE和NSE算法的稳定限值高于二阶无条件稳定显式CR和Chang算法。

4 虚拟RTHS数值仿真验证

4.1　线性RTHS数值仿真验证

为了进一步验证稳定性分析结果，进行虚拟RTHS数值仿真验证。首先进行线性RTHS数值仿真验证，结构对象为某SDOF结构，结构质量m为 1 000 kg，圆频率ω为20π rad·s^-1，阻尼比ξ为5%，刚度比例系数η为0.25，即对应图4e的工况。El Centro地震波作为该结构的输入。首先采用时间步长为0.01 s的CR算法对整体结构进行分析，分析结果作为参考解。接着采用CR、Chang、NDE和NSE 4个算法进行虚拟RTHS数值仿真，使用两种时间步长0.01 s和0.02 s，对应的Ω分别为0.628和1.257，通过调整α的大小来反映时滞的影响。

图5为结构的位移时程曲线，可以看出：由于时滞的存在，相当于系统中引入了负阻尼，因此虚拟RTHS的结果均大于参考解。图5a表明4个算法得到的位移响应均未失稳；图5b表明CR算法出现明显的失稳现象；图5c表明CR和Chang算法均失稳，CR的失稳现象更加严重，而NDE和NSE均未失稳。

图5  线性RTHS的位移时程曲线

Fig. 5  Displacement time history curves of linear RTHS

结合图4e的结果可知：对于图5a，当α=2，Ω=0.628时，CR算法的稳定限值最小，为0.8，但是仍高于0.628，因此没有出现失稳；对于图5b，当α=2，Ω=1.257时，CR算法的稳定限值虽然没有发生变化，但是该值已经小于1.257，于是发生失稳，而其他3个算法的稳定限值仍高于1.257，因此没有失稳；对于图5c，当α=2.25，Ω=0.628时，CR和Chang算法的稳定限值分别降低到0.430和0.582，均小于0.628，因此这两个算法出现失稳，而NDE和NSE的稳定限值在此时仍为无穷大，因此没有发生失稳。

4.2　非线性RTHS数值仿真验证

当试验子结构具有非线性时，应用RTHS具有一定的优势。因此，本文进行非线性RTHS的数值仿真验证。非线性结构的初始状态与4.1节的线性结构相同，区别在于试验子结构具有非线性刚度，其瞬时刚度表达式为$k_{\mathrm{E},i+\mathrm{1}}=k_{\mathrm{E},\mathrm{0}}\left[\mathrm{1}+\theta \sqrt[]{\left|x_{i+\mathrm{1}}^{\mathrm{\text{'}}}\right|}\right]$，其中$k_{\mathrm{E},\mathrm{0}}$为试验子结构的初始刚度，与4.1节的试验子结构刚度相同；θ为反映非线性程度的系数：当θ为正值时，试验子结构为刚度硬化，当θ为负值时，试验子结构为刚度软化，本文取θ为5与-5分别表示刚度硬化和刚度软化。刚度硬化意味着Ω的增大，对RTHS系统的稳定性不利；刚度软化则意味着Ω的减小，对RTHS系统的稳定性有利。

图6为两种不同非线性情形时RTHS的数值仿真结果。图6a和6b分别为当α=2.5， Δt=0.02 s时刚度软化和当α=1.5， Δt=0.02 s时刚度硬化的结果。由4.1节的结果可知，当α=2.25， Δt=0.01 s时，CR和Chang算法已经发生失稳，对于图6a中更大的时滞（α=2.5）和时间步长（Δt=0.02 s），这两个算法如果应用于线性RTHS时会更加失稳。然而，图6a显示刚度软化时，CR和Chang算法与NDE和NSE算法一样均未失稳，这证明了刚度软化在一定程度上提升了RTHS的稳定性。图6b显示在较小时滞（α=1.5）时，CR和Chang算法应用于刚度硬化的RTHS发生失稳，证明了刚度硬化会降低RTHS的稳定性。此外，NDE和NSE算法应用于该工况下，RTHS系统依然保持稳定，说明对于非线性RTHS系统，NDE和NSE算法比CR和Chang算法具有更高的稳定性。

图6  非线性RTHS的位移时程曲线

Fig. 6  Displacement time history curves of nonlinear RTHS

5 结论

（1）NDE和NSE算法分别与二阶无条件稳定显式CR和Chang算法具有相同的速度、位移表达式，而NDE和NSE算法具有更高阶精度。

（2）基于离散控制理论可以分析RTHS系统的稳定性，与时滞微分方程推导的解相比，该方法可以考虑积分算法对RTHS系统稳定性的影响。

（3）与使用CR和Chang算法的RTHS系统相比，使用NDE和NSE算法的RTHS系统具有更高的稳定限值，结合二者的高精度，意味着两种新算法在应用于RTHS时有着较大的优势和潜力。

作者贡献声明

傅 博：论文修改，数据分析。

张付泰：论文撰写，编程计算。

陈 瑾：论文修改，数据校核。