网刊加载中。。。

使用Chrome浏览器效果最佳,继续浏览,你可能不会看到最佳的展示效果,

确定继续浏览么?

复制成功,请在其他浏览器进行阅读

一种新型高精度半显式基于模型的积分算法  PDF

  • 傅博
  • 张付泰
  • 陈瑾
长安大学 建筑工程学院,陕西 西安 710061

中图分类号: TU311.3

最近更新:2023-05-09

DOI:10.11908/j.issn.0253-374x.21547

  • 全文
  • 图表
  • 参考文献
  • 作者
  • 出版信息
EN
目录contents

摘要

提出了一种新型高精度基于模型的积分算法,该算法采用了半显式Chang算法的速度、位移表达式,但是具有不同的积分参数。采用离散控制理论中“零极点匹配”和二阶Pade近似推导出新算法的积分参数。对算法的精度、稳定性、周期延长和振幅衰减等数值特性进行分析。与Newmark族积分算法和两种典型的基于模型的积分算法相比,新算法具备四阶精度,并且有非常小的周期误差。最后,通过三个数值算例验证了新算法的高精度特性。

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

为了同时满足显式和无条件稳定这两个特性,近年来学者提出一类新型积分算法,因为其积分参数与模型质量、刚度和阻尼相关,故被称为基于模型的积分算

1。基于模型积分算法计算效率高,主要用于混合试2、实时混合模3、子结构振动台试4及倒塌模5中。基于模型的积分算法根据其速度差分方程的不同又可以分为双显式和半显式算法,双显式指的是其速度、位移差分方程均为显式,半显式指的是仅位移差分方程是显式,但速度差分方程是隐式。Chang6最早提出一种基于模型的Chang算法,该算法为半显式算法,具有二阶精度。随后,Chen和Ricles7提出了另一种基于模型的CR算法,该算法为双显式,同样具有二阶精度。Gui8等基于CR算法的表达式,提出一族双显式的Gui算法,Gui算法与γ=1/2的Newmark族算法具有相同的数值特性,具有二阶精度。Fu9等提出一族与Newmark族算10具有相同数值特性的双显式GCR算法,该族算法的精度不超过二阶。

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

1 新算法的建立

1.1 已有积分算法简介

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

mx¨i+1+cx˙i+1+kxi+1=Fi+1 (1)

式中:m, c, k分别为质量,阻尼系数和刚度,c=2mωξ,其中ξ为阻尼比,ω为圆频率,ω=k/mx¨i+1,x˙i+1,xi+1分别为i+1时刻的加速度,速度和位移;Fi+1i+1时刻的外界激励。

积分算法可用于求解式(1),以经典的Newmark族积分算法为例,其速度、位移差分方程如下:

x˙i+1=x˙i+Δt1-γx¨i+γx¨i+1 (2)
xi+1=xi+Δtx˙i+Δt20.5-βx¨i+βx¨i+1 (3)

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

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

x˙i+1=x˙i+α1Δtx¨i (4)
xi+1=xi+Δtx˙i+α2Δt2x¨i (5)

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

α1=α2=4ω2Δt2+4ξωΔt+4 (6)

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

x˙i+1=x˙i+0.5Δtx¨i+0.5Δtx¨i+1 (7)
xi+1=xi+α1Δtx˙i+α2Δt2x¨i (8)

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

α1=4(1+ξωΔt)ω2Δt2+4ξωΔt+4 (9)
α2=2ω2Δt2+4ξωΔt+4 (10)

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

1.2 离散控制理论方法

采用离散控制理论的方法来确定基于模型的积分算法的积分参数,该方法在文

7-9已经成功应用。

根据离散控制理论,式(7),(8)和式(1)的Z变换形式如下:

X(z)=z-1X(z)+z-1α1ΔtXv(z)+z-1α2Δt2Xa(z) (11)
Xv(z)=z-1Xv(z)+12z-1ΔtXa(z)+12ΔtXa(z) (12)
mXa(z)+cXv(z)+kX(z)=F(z) (13)

式中:X(z)Xv(z)Xa(z)F(z)分别为xi+1,x˙i+1,x¨i+1,fi+1Z变换形式。

联立式(11)~(13),得

G(z)=X(z)F(z)=n2z2+n1z+n0d2z2+d1z+d0 (14)

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

表1  离散传递函数的分子和分母系数
Tab.1  Numerator and denominator coefficients of the discrete transfer function
分子系数分母系数
n2 0 d2 m(2ξωΔt+2)
n1 (α1+2α2)Δt2 d1 m[(α1+2α2)ω2Δt2-4]
n0 (α1-2α2)Δt2 d0 m[(α1-2α2)ω2Δt2-2ξωΔt+2]

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

z=esΔt (15)

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

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

z12+6sΔt+(sΔt)212-6sΔt+(sΔt)2 (16)

采用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近似的离散传递函数极点如下:

z1,212+6sΔt+(sΔt)212-6sΔt+(sΔt)2 (17)

式中:s1,2=-ξω±iω(1-ξ2)。因此离散的z1,2的值如下:

z1,212+6s-ξω±iω(1-ξ2)Δt+-ξω±iω(1-ξ2)Δt212-6s-ξω±iω(1-ξ2)Δt+-ξω±iω(1-ξ2)Δt2 (18)

根据式(14)特征方程的零点,特征方程可表达为

d2z12+d1z1+d0=0d2z22+d1z2+d0=0 (19)

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

α1=144ξΩ+144Ω4+12ξΩ3+(48ξ2+12)Ω2+144ξΩ+144 (20)
α2=72-2ξΩ3+(72ξ-96ξ3)ΩΩ4+12ξΩ3+(48ξ2+12)Ω2+144ξΩ+144 (21)

式中:Ω=ωΔt

1.4 新算法积分参数推导(MDOF体系)

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

X˙i+1=X˙i+12ΔtX¨i+12ΔtX¨i+1 (22)
Xi+1=Xi+Δtα1X˙i+Δt2α2X¨i (23)
MX¨i+1+CX˙i+1+KXi+1=Fi+1 (24)

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

Y˙i+1=Y˙i+12ΔtY¨i+12ΔtY¨i+1 (25)
Yi+1=Yi+Δtα1*Y˙i+Δt2α2*Y¨i (26)
M*Y¨i+1+C*Y˙i+1+K*Yi+1=ΦTFi+1 (27)

式中:Φ=ϕ1ϕ2ϕn是模态矩阵;ϕj(j=1n)是第j阶模态的特征向量;Y是模态坐标中的位移矢量,与X的关系为X=ΦYα1*=Φ-1α1Φα2*=Φ-1α2Φ是对角积分参数矩阵;M*=ΦTMΦ,C*=ΦTCΦK*=ΦTKΦ分别为模态质量,阻尼和刚度矩阵。因此,第j阶模态的积分参数可以表示为SDOF体系的积分参数,其形式如下:

α1 j*=72ΔtCj*Mj*+144Mj*2Δt4Kj*2+6Δt3Kj*Cj*+12Δt2Cj*2+12Δt2Mj*Kj*+72ΔtCj*Mj*+144Mj*2 (28)
α2 j*=72Mj*2-Δt3Kj*Cj*+(36Δt-48ξj2Δt)Cj*Mj*Δt4Kj*2+6Δt3Kj*Cj*+12Δt2Cj*2+12Δt2Mj*Kj*+72ΔtCj*Mj*+144Mj*2 (29)

式中:Mj*,Cj*=2Mj*ξjωj,Kj*=Mj*ωj2分别为第j阶模态的质量,阻尼和刚度矩阵;ξj,ωj分别为第j阶模态的阻尼比和圆频率。因此,对角积分参数矩阵α1*,α2*形式如下:

α1*=α-1[72ΔtC*M*+144M*2] (30)
α2*=α-1[72M*2-Δt3K*C*+(36Δt-48ξ2Δt)C*M*] (31)

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

α1=Φα-1[72ΔtC*M*+144M*2]Φ-1 (32)
α2=Φα-1[72M*2-Δt3K*C*+(36Δt-48ξ2Δt)C*M*]Φ-1 (33)

式中:α=Δt4K*2+6Δt3K*C*+12Δt2C*2+12Δt2M*Κ*+72ΔtC*M*+144M*2

2 数值特性分析

2.1 精度分析

首先将式(1), (7)和(8)写成递推矩阵的形式:

xi+1Δtx˙i+1Δt2x¨i+1=AxiΔtx˙iΔt2x¨i+00Δt2/mfi+1 (34)

式中:A为放大矩阵,其值为

Α=1α1α2-ω2Δt22(ξωΔt+1)-α1ω2Δt2-22(ξωΔt+1)-α2ω2Δt2-12(ξωΔt+1)-ω2Δt2ξωΔt+1-α1ω2Δt2+2ξωΔtξωΔt+1-α2ω2Δt2+ξωΔtξωΔt+1

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

λ3-2A1λ2+A2λ-A3=0 (35)

式中:λ为放大矩阵A的特征值;A1A2A3的表达式为

A1=(4-α1Ω2-2α2Ω2)/4(ξΩ+1)A2=(α1Ω2-2α2Ω2-2ξΩ+2)/2(ξΩ+1)A3=0

离散时间系统的位移数值解与连续时间系统的位移精确解的差值称为局部截断误

11,表示如下:

τ=[x(t+Δt)-2A1x(t)+A2x(t-Δt)-A3x(t-2Δt)]/Δt2 (36)

假设位移任何阶导数均连续可微,则式(36)的有限阶的泰勒级数展开如下所示:

τ=l=0LTlΔtl-2x(l)+OΔtL-1) (37)

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

T0=1-2A1+A2-A3Tl=1l![1+(-1)lA2-(-2)lA3],l1 (38)

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

τ=d1d3x(t)+d2d3x˙(t)+O(Δt4) (39)

式中:

d1=96ω2[(ξ2-58)(ξ2-112)Ω2+2ξ4-ξ2]Ω4d2=192ωξ[(ξ4-2324ξ2+16)Ω2+2ξ4-32ξ2+18]Ω4d3=60Ω4+720ξΩ3+(2 880ξ2+720)Ω2+8 640ξΩ+8 640

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

图2  CR, Chang和本文方法对位移、速度和加速度的收敛速率对比

Fig.2  Comparison of the convergence rates of CR, Chang and this paper's methods for displacement, velocity and acceleration

2.2 稳定性分析

如果积分算法放大矩阵的谱半径ρ=max{λ1,λ2,λ3}1,那么该算法就是稳定的,但是如果三个特征根中有重根时,那么仅当ρ<1时,算法才是稳定的。

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

λ1,2=σ±εi,λ3=0ρ=σ2+ε2σ=144+Ω4+(48ξ2-60)Ω2144+Ω4+12ξΩ3+(48ξ2+12)Ω2+144ξΩε=12Ω(Ω2-12)1-ξ2144+Ω4+12ξΩ3+(48ξ2+12)Ω2+144ξΩ (40)

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

图3  新算法的谱半径

Fig.3  Spectra radius of new algorithm

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

2.3 周期延长和振幅衰减

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

PE=(T¯n-Tn)/Tn=Ω/Ω¯-1 (41)

式中:T¯n=2π/ω¯, ω¯=Ω¯/ΔtTn为SDOF结构体系无阻尼自振周期。振幅衰减AD如图4定义。

图4  振幅衰减的定义

Fig.4  Definition of amplitude decay (AD)

振幅衰减AD还可以通过算法的数值阻尼ξ¯来量化,其表达式如下:

ξ¯=-ln(σ2+ε2)2Ω¯Ω¯=tan-1(ε/σ)1-ξ¯2 (42)

正是由于数值阻尼ξ¯的存在才造成了振幅衰减现象。当ξ=0时,新算法的数值阻尼为ξ¯=0,故不存在振幅衰减,该数值特性与Chang、CR、CAA、LA、NE等算法一致。图5为本文方法与几种已有算法的周期延长对比,可以看出,本文方法相较于几种已有算法,具有极低的周期延长,这也从侧面反映出本文方法的精度高。

图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,初始速度x˙0=1s-1。结构模型如图6所示。

图6  无阻尼线弹性SDOF结构的简图

Fig.6  Schematic diagram of a linear elastic SDOF system without damping

无阻尼线弹性SDOF结构的自由振动的位移解析

12如下:

x(t)=x(0)cosωt+x˙(0)ωsinωt (43)

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

图7  位移时程曲线

Fig.7  Time history curves of displacement

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

NEE=i=1nxR,i2-i=1nxN,i2i=1nxN,i2 (44)
NRMSE=i=1n(xR,i-xN,i)2nmax(xN)-min(xN) (45)

式中:xNxR分别代表数值(Numerical)结果和参考(Reference)结果;N是样本数目。NEENRMSE分别对幅值和相位误差比较敏感。表2列出了位移误差指标。

表2  位移误差指标
Tab.2  Error indices of displacement ( % )
Δt算法NEENRMSE
0.02s CR 1.459 3 10.170 0
Chang 1.522 3 10.217 9
本文方法 0.002 8 0.053 8
0.05s CR 7.588 7 44.888 4
Chang 4.323 7 46.255 2
本文方法 0.041 3 0.727 6

表2可以看出,相同时间步长的情况下,Chang算法的NEE小于CR算法,二者的NRMSE接近,而本文方法的NEENRMSE均远小于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
阶数Matlab/Simulink模型/s

OpenSEES

模型/ s

误差/%
1 0.306 945 0.306 991 0.015 0
2 0.076 473 0.076 019 0.597 2
3 0.019 860 0.020 079 1.090 7
4 0.017 239 0.017 434 1.118 5
5 0.008 600 0.008 620 0.232 0

表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算法的计算结果很接近,误差指标也非常接近,而本文方法的结果与参考解更契合,NEENRMSE均小于CR算法和Chang算法,其中本文方法的NRMSE优势非常明显,不及CR算法和Chang算法的5%。即便在较大时间步长(Δt=0.02s)时,本文方法的NEE也不到3%,NRMSE不到0.3%,这说明本文方法在较大时间步长时也可以达到较高精度。表5给出了本算例中各算法在不同时间步长的计算时间,使用的计算机配置为CPU Intel Core i5-

表4  位移误差指标
Tab.4  Error indices of displacement
时间步长Δt算法NEENRMSE
0.005s CR 1.055 4 0.441 2
Chang 1.046 1 0.442 4
本文方法 0.328 7 0.032 9
0.01s CR 3.216 1 1.765 4
Chang 3.191 1 1.769 9
本文方法 0.823 2 0.067 6
0.02s CR 4.818 7 6.272 1
Chang 4.723 0 6.289 1
本文方法 2.833 8 0.245 9
表5  各算法的计算时间对比
Tab.5  Comparison of computation time of different algorithms
算法时间步长Δt
0.0050.010.02
CAA 212.18 107.50 79.07
CR 67.718 33.90 17.32
Chang 66.19 33.27 17.26
本文方法 65.37 34.45 17.19

6 300HQ @2.30GHz,内存8G。

表5可以看出,当时间步长相同时,本文方法与CR算法、Chang算法的计算时间接近,但是计算时间明显小于隐式CAA算法。结合表45可知,Δ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×105m-1i=1,2,3···N), x¨g=10sin(10t) m·s-2 t=0···5s), 阻尼矩阵采用瑞利阻尼,假定前两阶模态阻尼比为2%,采用式(46)的非线性刚度,则

ki,t=ki,01-xi,t-xi-1,t (46)

式中:kitt时刻的第i个弹簧的刚度;xitt时刻的第i个质量的位移;xi-1,tt时刻的第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
自由度数方法计算时间x50计算误差/ %
NEENRMSE
50 CAA 0.059 2 35.469 1 7.469 6
振型叠加法 0.009 3 28.464 8 5.013 8
CR 0.011 1 22.297 2 7.166 4
Chang 0.016 2 22.697 2 7.314 5
本文方法 0.023 0 3.410 9 0.618 0
100 CAA 0.196 1 5.334 3 6.170 1
振型叠加法 0.008 9 1.237 4 4.341 4
CR 0.027 1 16.991 1 6.851 9
Chang 0.041 5 16.548 4 6.958 3
本文方法 0.047 0 1.448 9 0.817 3
400 CAA 2.336 4 16.058 5 4.411 9
振型叠加法 0.018 9 18.450 3 3.542 3
CR 0.216 0 4.738 6 4.101 9
Chang 0.474 3 7.760 8 4.379 1
本文方法 0.536 6 3.688 4 0.779 3
700 CAA 6.384 6 12.099 5 3.307 9
振型叠加法 0.021 7 16.462 1 2.483 7
CR 0.727 4 4.466 1 3.000 4
Chang 1.523 7 7.660 1 3.314 4
本文方法 1.737 2 3.700 3 1.196 0
1 000 CAA 15.463 9 8.951 9 2.275 6
振型叠加法 0.032 2 13.633 1 1.961 1
CR 1.410 1 3.704 6 2.012 7
Chang 3.090 7 6.541 5 2.011 2
本文方法 3.685 1 3.856 8 0.789 7

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

4 结论

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

作者贡献声明

傅博:论文修改,数据分析;

张付泰:论文撰写,编程计算;

陈瑾:论文修改,数据校核。

参考文献

1

KOLAY CRICLES J M. Assessment of explicit and semi-explicit classes of model-based algorithms for direct integration in structural dynamics [J]. International Journal for Numerical Methods in Engineering20161071):49. [百度学术] 

2

彭天波谢馨曾忠. 采用Chang方法的混合试验的稳定性和精度[J]. 同济大学学报(自然科学版)20144212): 1790. [百度学术] 

PENG TianboXIE XinZENG Zhonget al. Stability and accuracy of sharking table-actuator hybrid test with Chang method [J]. Journal of Tongji University (Natural Science)20144212): 1790. [百度学术] 

3

KOLAY CRICLES J M. Force-based frame element implementation for real-time hybrid simulation using explicit direct integration algorithms [J]. Journal of Structural Engineering20181442): 04017191. [百度学术] 

4

傅博蒋欢军. 具有高稳定性的剪切型子结构振动台试验[J]. 同济大学学报(自然科学版)20174511): 1602. [百度学术] 

FU BoJIANG Huanjun. Shear-type substructure sharking table testing method with high stability [J]. Journal of Tongji University (Natural Science)20174511): 1602. [百度学术] 

5

FENG D CKOLAY CRICLES J Met al. Collapse simulation of reinforced concrete frame structures [J]. The Structural Design of Tall and Special Buildings20162512): 578. [百度学术] 

6

CHANG S Y. Explicit pseudodynamic algorithm with unconditional stability [J]. Journal of Engineering Mechanics20021289): 935. [百度学术] 

7

CHEN CRICLES J M. Development of direct integration algorithms for structural dynamics using discrete control theory [J]. Journal of Engineering Mechanics20081348): 676. [百度学术] 

8

GUI YWANG J TJIN Fet al. Development of a family of explicit algorithms for structural dynamics with unconditional stability [J]. Nonlinear Dynamics2014774): 1157. [百度学术] 

9

FU BFENG D CJIANG H. A new family of explicit model-based integration algorithms for structural dynamic analysis [J], International Journal of Structural Stability and Dynamics196): 1950053. [百度学术] 

10

NEWMARK N M. A method of computation for structural dynamics [J]. Journal of the Engineering Mechanics 1959853): 6794. [百度学术] 

11

HILBER H M. Analysis and design of numerical integration methods in structural dynamics [R]. BerkeleyEarthquake Engineering Research Center1976. [百度学术] 

12

CHOPRA A K. Dynamics of structures: theory and applications to earthquake engineering [M]. 4th ed. Upper Saddle RiverPrentice Hall2011. [百度学术] 

13

The MathWorks Inc. Matlab R2014b [EB/OL].[2021-08-29http://www.mathworks.com. [百度学术] 

14

傅博蒋欢军. 使用基于模型的积分算法的钢框架非线性时程分析[J]. 土木工程学报201851S1): 63. [百度学术] 

FU BoJIANG Huanjun. Nonlinear time history analysis of steel frame using model-based integration algorithm [J]. China Civil Engineering Journal201851S1): 63. [百度学术]