傅 博， 张付泰， 陈 瑾

（长安大学 建筑工程学院，陕西 西安 710061）

直接逐步积分方法被广泛应用于求解离散化的结构动力学运动方程。根据位移差分方程的表达式不同，可以将积分算法分为显式和隐式积分算法。在求解非线性结构动力学问题时，显式积分算法无需进行迭代，节省了计算时间，因此具有较高的计算效率。根据积分时间步长选取是否有限制，积分算法又可以分为无条件稳定和条件稳定积分算法。当求解具有多自由度的复杂结构动力学问题时，条件稳定积分算法可能需要极小的积分步长以满足稳定，因此极大地降低了计算效率。综上，如果一个算法同时兼具显式和无条件稳定，那么将具有非常高的计算效率。但是，通常显式积分算法是条件稳定，而无条件稳定算法往往是隐式。

为了同时满足显式和无条件稳定这两个特性，近年来学者提出一类新型积分算法，因为其积分参数与模型质量、刚度和阻尼相关，故被称为基于模型的积分算法［1］。基于模型积分算法计算效率高，主要用于混合试验［2］、实时混合模拟［3］、子结构振动台试验［4］及倒塌模拟［5］中。基于模型的积分算法根据其速度差分方程的不同又可以分为双显式和半显式算法，双显式指的是其速度、位移差分方程均为显式，半显式指的是仅位移差分方程是显式，但速度差分方程是隐式。Chang［6］最早提出一种基于模型的Chang 算法，该算法为半显式算法，具有二阶精度。随后，Chen 和Ricles［7］提出了另一种基于模型的CR算法，该算法为双显式，同样具有二阶精度。Gui［8］等基于CR算法的表达式，提出一族双显式的Gui算法，Gui 算法与γ=1/2 的Newmark 族算法具有相同的数值特性，具有二阶精度。Fu［9］等提出一族与Newmark 族算法［10］具有相同数值特性的双显式GCR算法，该族算法的精度不超过二阶。

综上，与传统积分算法相比，基于模型的积分算法虽然在计算效率上有很大优势，但是在精度上和常规方法类似。为此本文提出一种高精度的基于模型的积分算法，该算法既保持了基于模型的积分算法的显式和无条件稳定的数值特性，同时具有更高的数值精度。

1 新算法的建立

1.1 已有积分算法简介

对于线弹性的单自由度（SDOF）结构体系，离散时间系统下结构运动方程可表示为

式中：m，c，k分别为质量，阻尼系数和刚度，c=2mωξ，其中ξ为阻尼比，ω为圆频率，分别为i+1时刻的加速度，速度和位移；Fi+1为i+1时刻的外界激励。

积分算法可用于求解式（1），以经典的Newmark族积分算法为例，其速度、位移差分方程如下：

式中：Δt是积分时间步长，β，γ是积分参数，直接影响积分算法的数值特性。常用的Newmark 算法包括常平均加速度法（CAA，β=1/4，γ=1/2），线性加速度法（LA，β=1/6，γ=1/2），显式算法（NE，β=0，γ=1/2）。由于β和γ与结构模型无关，所以Newmark算法属于模型不相关的积分算法。

与模型不相关的积分算法不同，基于模型的积分算法的积分参数与模型相关，以CR 算法为例，其速度、位移差分方程如下：

式中：α1，α2为与模型相关的积分参数，其值为

类似地，Chang 算法的速度、位移差分方程如下：

式中：α1，α2亦为与模型相关的积分参数，其值为

本文提出一种新型高精度半显式基于模型的积分算法，该算法采用Chang算法相同的速度、位移差分方程，但是具有不同的积分参数。下文将对新算法的积分参数进行推导。

1.2 离散控制理论方法

采用离散控制理论的方法来确定基于模型的积分算法的积分参数，该方法在文献［7-9］已经成功应用。根据离散控制理论，式（7），（8）和式（1）的Z 变换形式如下：

联立式（11）～（13），得

其中，分子分母系数如表1所示。

表1 离散传递函数的分子和分母系数Tab.1 Numerator and denominator coefficients of the discrete transfer function

一种称为“零极点匹配”的离散化方法在控制理论中用来从连续域的极点映射到离散域的极点，这种精确的映射法则规定如下：

式中：Δt为采样周期，可以令其与积分算法的时间步长相等。式（15）的映射法则不太适合实际应用，因此通常采用近似的映射法则。

本文新算法采用2 阶pade 近似其映射规则如下：

采用2阶pade近似与式（15）的逼似程度远优于用于CR算法和Chang算法的1阶pade近似，对比情况如图1 所示，可以看出，2 阶pade 近似的效果要远优于1阶pade近似，尤其当sΔt＜-0.5时。

图1 两种Pade近似与指数函数对比Fig.1 Comparisons of two Pade approximations with exponential function

1.3 新算法积分参数推导（SDOF体系）

采用2阶pade近似的离散传递函数极点如下：

根据式（14）特征方程的零点，特征方程可表达为

将表1中离散传递函数的分母系数d2，d1，d0和式（18），同时代入式（19），可以联立求解得到积分参数α1和α2的值为

1.4 新算法积分参数推导（MDOF体系）

对于具有n个自由度的多自由度（MDOF）线弹性系统，新算法的速度、位移微分方程和运动方程如下：

假定阻尼矩阵是经典阻尼，由于结构振型的正交性，可将式（22）， （23）和（24）写成模态坐标系下的表达式，即

因此，可以进一步得到基于模型的积分参数矩阵α1和α2如下：

2 数值特性分析

2.1 精度分析

首先将式（1）， （7）和（8）写成递推矩阵的形式：

式中：A为放大矩阵，其值为

放大矩阵A的特征方程按照|A—λI|=0进行计算如下：

式中：λ为放大矩阵A的特征值；A1、A2、A3的表达式为

离散时间系统的位移数值解与连续时间系统的位移精确解的差值称为局部截断误差［11］，表示如下：

假设位移任何阶导数均连续可微，则式（36）的有限阶的泰勒级数展开如下所示：

式中：x(l)为位移的第l阶微分，Tl的表达式如下所示：

这里假设L=5，可以得到新算法的局部截断误差表达式如下：

式中：

从系数d1，d2，d3可以看出，新算法具有4 阶精度，而Chang、CR、CAA、LA、NE 等算法均为2 阶精度。为了进一步验证本文提出的新算法具有4 阶精度，用无阻尼线弹性单自由度结构的自由振动来给出新算法的绝对误差收敛速率，如图2所示。取，计算时间tn=1s，计算次数N=10，100，1 000，10 000。

2.2 稳定性分析

本文方法放大矩阵的特征值可由式（35）求得

本文方法的谱半径如图3所示。

图3 新算法的谱半径Fig.3 Spectra radius of new algorithm

由图3可知，当ξ＞0时，谱半径ρ＜1，表明本文方法在有阻尼情况下无条件稳定；当ξ=0 时，谱半径ρ=1，当时，ε=0，意味着有两个重根λ1=λ2=—1，表明在该特例下算法不稳定，而对于其他Ω值，本文方法均稳定。因为真实结构均存在阻尼，因此，本文方法仍可视作无条件稳定的算法。

2.3 周期延长和振幅衰减

算法的周期延长（PE）和振幅衰减（AD）可以反映积分算法的数值准确度。算法的周期延长PE定义如下：

图4 振幅衰减的定义Fig.4 Definition of amplitude decay (AD)

图5 多种积分算法的周期延长PE对比Fig.5 Comparison of period elongation (PE) of various integration algorithms

3 数值算例

3.1 无阻尼线弹性SDOF结构的自由振动

采用本文方法、Chang 和CR 算法对无阻尼线弹性SDOF 结构的自由振动进行求解，采用时间步长为Δt=0.02s 和Δt=0.05s 进行计算。该结构质量m=20kg，k=2 000N·s-1， 因此结构的自振频率ω=10rad·s-1。结构的初始位移x0=0m，初始速度结构模型如图6所示。

图6 无阻尼线弹性SDOF结构的简图Fig.6 Schematic diagram of a linear elastic SDOF system without damping

无阻尼线弹性SDOF 结构的自由振动的位移解析解［12］如下：

图7 为3 种积分算法采用两种不同积分步长计算得到的位移时程曲线。由图7可以看出，本文方法的计算结果与精确解非常接近，而CR、Chang算法的计算结果与精确解存在较为明显的差异，尤其是当Δt=0.05s时，由周期延长引起的相位误差较大。

图7 位移时程曲线Fig.7 Time history curves of displacement

进一步，采用以下两个指标来衡量积分算法的误差：

式中：xN和xR分别代表数值（Numerical）结果和参考（Reference）结果；N是样本数目。NEE和NRMSE分别对幅值和相位误差比较敏感。表2 列出了位移误差指标。

表2 位移误差指标Tab.2 Error indices of displacement%

由表2可以看出，相同时间步长的情况下，Chang算法的NEE小于CR 算法，二者的NRMSE接近，而本文方法的NEE和NRMSE均远小于CR 和Chang 算法；本文方法在较大时间步长（0.05s）时，误差还不到CR 和Chang算法在较小时间步长（0.02s）时的8%。

3.2 两层钢框架有限元模型的地震响应分析

选取某两层钢框架办公楼，该建筑的层高为4m，跨度为9m，侧向力由抗弯框架承担，重力由倚靠柱承担，抗弯框架与倚靠柱之间通过刚性楼板假定进行关联，其结构如图8所示。

图8 两层钢框架结构示意图Fig.8 Schematic diagram of two-story steel frame structure

基于Matlab/Simulink［13］，对该结构进行了有限元建模：材料密度为7.8×103kg·m—3，采用理想弹塑性材料本构关系，弹性模量2×1011Pa，屈服强度为3.45×108Pa。采用基于刚度的纤维梁单元模拟梁柱构件［14］，单元的积分点数为5，截面纤维数为24；采用弹性梁柱单元模拟倚靠柱并考虑P-Δ 效应。该结构模型共包括31 个节点，30 个单元以及83 个自由度。采用集中质量矩阵，阻尼矩阵采用瑞利阻尼，假定前两阶模态阻尼比为2%。

为了验证基于Matlab/Simulink所建立的有限元模型的可靠性，在OpenSEES中对该框架进行建模：材料本构采用uniaxialMaterial ElasticPP，梁、柱均采用dispBeamColumn单元，倚靠柱采用elasticBeamColumn单元，其余模型参数与Matlab/Simulink模型保持一致。Matlab/Simulink模型与OpenSEES模型前5阶周期对比如表3所示。

表3 Matlab/Simulink 模型与OpenSEES 模型前5 阶周期对比Tab.3 Comparison of first 5th order periods of the Matlab/Simulink model and the OpenSEES model

从表3 中可以看出，基于Matlab/Simulink 建立的有限元模型的前5 阶周期与基于OpenSEES 建立的有限元模型的各阶周期非常接近，误差最大仅为1.118 5%，从而说明本文基于Matlab/Simulink 建立的有限元模型的可靠性。

将El Centro NS， 1940作为结构的地震动输入，采用CR、Chang及本文方法进行时程分析，同时采用Δt=0.001s的CAA算法作为参考解。图9给出了不同时间步长（Δt=0.005s， 0.01s 和 0.02s）时各算法的首层位移时程曲线。

图9 首层位移时程曲线Fig.9 Time history curves of displacement at first story

从图9和表4可以看出，相同时间步长下，CR算法和Chang算法的计算结果很接近，误差指标也非常接近，而本文方法的结果与参考解更契合，NEE和NRMSE均小于CR算法和Chang算法，其中本文方法的NRMSE优势非常明显，不及CR 算法和Chang 算法的5%。即便在较大时间步长（Δt=0.02s）时，本文方法的NEE也不到3%，NRMSE不到0.3%，这说明本文方法在较大时间步长时也可以达到较高精度。表5 给出了本算例中各算法在不同时间步长的计算时间，使用的计算机配置为 CPU Intel Core i5-6 300HQ @2.30GHz，内存8G。

表4 位移误差指标Tab.4 Error indices of displacement

表5 各算法的计算时间对比Tab.5 Comparison of computation time of different algorithms

由表5 可以看出，当时间步长相同时，本文方法与CR算法、Chang算法的计算时间接近，但是计算时间明显小于隐式CAA 算法。结合表4、5 可知，Δt=0.02s时本文方法的NEE仍小于Δt=0.01s时CR算法和Chang 算法的NEE；Δt=0.02s 时本文方法的NRMSE仍远小于Δt=0.01s 时CR 算法和Chang 算法的NRMSE，这说明Δt=0.02s时本文方法的计算精度高于Δt=0.01s 时CR 算法和Chang 算法。Δt=0.01s 时CR 算法和Chang 算法分别耗时33.90s 和33.27s，而Δt=0.02s 时本文方法耗时仅为17.19s。这说明本文方法不仅节约了约50%的计算时间，还达到了更高的计算精度。从另一个角度来说，在保证相近的计算精度的前提下，本文方法的计算效率要高于CR算法和Chang算法。

3.3 多自由度非线性质量-弹簧系统的动力学问题

基于模型的积分算法非常重要的特性是无条件稳定性和显式表达式。无条件稳定性意味着无需为了保证算法的稳定性而可以选择较大的时间步长，显式格式则意味着求解非线性问题时无需迭代。因此，基于模型的积分算法计算效率高也对应于上述两个方面：一是较大的时间步长，从而计算步数较少；二是每一积分时间步的计算时间短，因为无需迭代。本节将选取图10中具有N个自由度的非线性质量-弹簧系统来说明本文方法在计算效率和计算精度方面的优势。

图10 质量-弹簧系统示意图Fig.10 Schematic diagram of Mass-Spring system

在图10的质量-弹簧系统中，mi=150kg，ki=5.5×105N·m-1（i=1，2，3···N），ẍg=10sin（10t） m·s—2（t=0···5s）， 阻尼矩阵采用瑞利阻尼，假定前两阶模态阻尼比为2%，采用式（46）的非线性刚度，则

式中：ki，t为t时刻的第i个弹簧的刚度；xi，t为t时刻的第i个质量的位移；xi—1，t为t时刻的第i—1 个质量的位移。考虑N=50，100，400，700，1 000 等5 种不同自由度数的工况，选取Δt=0.001s 的Newmark 显式算法作为参考解，对比Δt=0.05s的CAA算法、振型叠加法、CR 算法、Chang 算法和本文方法的计算结果。表6 给出了各算法的计算时间和x50的误差指标。

表6 各算法的计算时间和位移误差指标对比Tab.6 Comparison of computation time and error indices of displacement of different algorithms

从表6 中可以看出，CAA 算法的计算时间远大于本文方法和其他方法，这是因为CAA算法需要迭代，CAA 算法的NEE和NRMSE均大于本文方法，说明CAA 算法的计算精度和计算效率均逊于本文方法；振型叠加法的计算时间相对于本文方法和其他方法有明显的优势，但是该方法在求解非线性问题时有较大的误差，虽然在100自由度工况时，其NEE略小于本文方法，但是其NRMSE是本文方法的5倍，而在其余4 种自由度工况下，振型叠加法的NEE、NRMSE均大于本文方法；与CR 算法和Chang 算法相比，本文方法计算时间略长，这是因为在计算积分参数矩阵的过程中需要额外的计算时间，但是在计算精度上，本文方法的NEE和NRMSE基本上都要小于CR算法和Chang算法，只有在1 000 自由度工况下，CR 算法的NEE略小于本文方法，但是其NRMSE是本文方法的2.5 倍。综上，本文方法在计算精度上有明显的综合优势。

4 结论

本文提出了一种新型高精度半显式基于模型的积分算法。新算法与半显式Chang 算法具有相同的速度、位移差分方程，采用二阶Pade近似极点映射方法生成新算法的积分参数。对新算法的精度、稳定性和周期延长和振幅衰减等数值特性进行分析，发现本文方法具有四阶精度，周期延长极小，并且没有能量耗散。通过3个具有代表性的数值算例，进一步论证了本文方法的优越性。

作者贡献声明：

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

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

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