摘要
提出了一种新型高精度基于模型的积分算法,该算法采用了半显式Chang算法的速度、位移表达式,但是具有不同的积分参数。采用离散控制理论中“零极点匹配”和二阶Pade近似推导出新算法的积分参数。对算法的精度、稳定性、周期延长和振幅衰减等数值特性进行分析。与Newmark族积分算法和两种典型的基于模型的积分算法相比,新算法具备四阶精度,并且有非常小的周期误差。最后,通过三个数值算例验证了新算法的高精度特性。
直接逐步积分方法被广泛应用于求解离散化的结构动力学运动方程。根据位移差分方程的表达式不同,可以将积分算法分为显式和隐式积分算法。在求解非线性结构动力学问题时,显式积分算法无需进行迭代,节省了计算时间,因此具有较高的计算效率。根据积分时间步长选取是否有限制,积分算法又可以分为无条件稳定和条件稳定积分算法。当求解具有多自由度的复杂结构动力学问题时,条件稳定积分算法可能需要极小的积分步长以满足稳定,因此极大地降低了计算效率。综上,如果一个算法同时兼具显式和无条件稳定,那么将具有非常高的计算效率。但是,通常显式积分算法是条件稳定,而无条件稳定算法往往是隐式。
为了同时满足显式和无条件稳定这两个特性,近年来学者提出一类新型积分算法,因为其积分参数与模型质量、刚度和阻尼相关,故被称为基于模型的积分算
综上,与传统积分算法相比,基于模型的积分算法虽然在计算效率上有很大优势,但是在精度上和常规方法类似。为此本文提出一种高精度的基于模型的积分算法,该算法既保持了基于模型的积分算法的显式和无条件稳定的数值特性,同时具有更高的数值精度。
对于线弹性的单自由度(SDOF)结构体系,离散时间系统下结构运动方程可表示为
(1) |
式中:m, c, k分别为质量,阻尼系数和刚度,,其中为阻尼比,为圆频率,;分别为i+1时刻的加速度,速度和位移;为i+1时刻的外界激励。
积分算法可用于求解
(2) |
(3) |
式中:是积分时间步长,, 是积分参数,直接影响积分算法的数值特性。常用的Newmark算法包括常平均加速度法(CAA, =1/4,=1/2),线性加速度法(LA, =1/6,=1/2),显式算法(NE, =0,=1/2)。由于和与结构模型无关,所以Newmark算法属于模型不相关的积分算法。
与模型不相关的积分算法不同,基于模型的积分算法的积分参数与模型相关,以CR算法为例,其速度、位移差分方程如下:
(4) |
(5) |
式中:为与模型相关的积分参数,其值为
(6) |
类似地,Chang算法的速度、位移差分方程如下:
(7) |
(8) |
式中:亦为与模型相关的积分参数,其值为
(9) |
(10) |
本文提出一种新型高精度半显式基于模型的积分算法,该算法采用Chang算法相同的速度、位移差分方程,但是具有不同的积分参数。下文将对新算法的积分参数进行推导。
采用离散控制理论的方法来确定基于模型的积分算法的积分参数,该方法在文
根据离散控制理论,式(
(11) |
(12) |
(13) |
式中:,,,分别为, 的Z变换形式。
(14) |
其中,分子分母系数如
分子系数 | 值 | 分母系数 | 值 |
---|---|---|---|
0 | |||
一种称为“零极点匹配”的离散化方法在控制理论中用来从连续域的极点映射到离散域的极点,这种精确的映射法则规定如下:
(15) |
式中:为采样周期,可以令其与积分算法的时间步长相等。
本文新算法采用2阶pade近似其映射规则如下:
(16) |
采用2阶pade近似与
图1 两种Pade近似与指数函数对比
Fig.1 Comparisons of two Pade approximations with exponential function
采用2阶pade近似的离散传递函数极点如下:
(17) |
式中:。因此离散的的值如下:
(18) |
根据
(19) |
将
(20) |
(21) |
式中:。
对于具有n个自由度的多自由度(MDOF)线弹性系统,新算法的速度、位移微分方程和运动方程如下:
(22) |
(23) |
(24) |
假定阻尼矩阵是经典阻尼,由于结构振型的正交性,可将式(
(25) |
(26) |
(27) |
式中:是模态矩阵;是第阶模态的特征向量;是模态坐标中的位移矢量,与的关系为;和是对角积分参数矩阵;分别为模态质量,阻尼和刚度矩阵。因此,第阶模态的积分参数可以表示为SDOF体系的积分参数,其形式如下:
(28) |
(29) |
式中:分别为第阶模态的质量,阻尼和刚度矩阵;分别为第阶模态的阻尼比和圆频率。因此,对角积分参数矩阵形式如下:
(30) |
(31) |
因此,可以进一步得到基于模型的积分参数矩阵和如下:
(32) |
(33) |
式中:
(34) |
式中:为放大矩阵,其值为
放大矩阵的特征方程按照进行计算如下:
(35) |
式中:为放大矩阵的特征值;A1、A2、A3的表达式为
离散时间系统的位移数值解与连续时间系统的位移精确解的差值称为局部截断误
(36) |
假设位移任何阶导数均连续可微,则
(37) |
式中:为位移的第阶微分,的表达式如下所示:
(38) |
这里假设L=5,可以得到新算法的局部截断误差表达式如下:
(39) |
式中:
从系数可以看出,新算法具有4阶精度,而Chang、CR、CAA、LA、NE等算法均为2阶精度。为了进一步验证本文提出的新算法具有4阶精度,用无阻尼线弹性单自由度结构的自由振动来给出新算法的绝对误差收敛速率,如
图2 CR, Chang和本文方法对位移、速度和加速度的收敛速率对比
Fig.2 Comparison of the convergence rates of CR, Chang and this paper's methods for displacement, velocity and acceleration
如果积分算法放大矩阵的谱半径,那么该算法就是稳定的,但是如果三个特征根中有重根时,那么仅当时,算法才是稳定的。
本文方法放大矩阵的特征值可由
(40) |
本文方法的谱半径如
图3 新算法的谱半径
Fig.3 Spectra radius of new algorithm
由
算法的周期延长(PE)和振幅衰减(AD)可以反映积分算法的数值准确度。算法的周期延长PE定义如下:
(41) |
式中:,为SDOF结构体系无阻尼自振周期。振幅衰减AD如
图4 振幅衰减的定义
Fig.4 Definition of amplitude decay (AD)
振幅衰减AD还可以通过算法的数值阻尼来量化,其表达式如下:
(42) |
正是由于数值阻尼的存在才造成了振幅衰减现象。当时,新算法的数值阻尼为,故不存在振幅衰减,该数值特性与Chang、CR、CAA、LA、NE等算法一致。
图5 多种积分算法的周期延长PE对比
Fig.5 Comparison of period elongation (PE) of various integration algorithms
采用本文方法、Chang和CR算法对无阻尼线弹性SDOF结构的自由振动进行求解,采用时间步长为s和s进行计算。该结构质量m=20kg, k=2 000N·
图6 无阻尼线弹性SDOF结构的简图
Fig.6 Schematic diagram of a linear elastic SDOF system without damping
无阻尼线弹性SDOF结构的自由振动的位移解析
(43) |
图7 位移时程曲线
Fig.7 Time history curves of displacement
进一步,采用以下两个指标来衡量积分算法的误差:
(44) |
(45) |
式中:和分别代表数值(Numerical)结果和参考(Reference)结果;N是样本数目。NEE和NRMSE分别对幅值和相位误差比较敏感。
Δt | 算法 | NRMSE | |
---|---|---|---|
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 |
由
选取某两层钢框架办公楼,该建筑的层高为4m,跨度为9m,侧向力由抗弯框架承担,重力由倚靠柱承担,抗弯框架与倚靠柱之间通过刚性楼板假定进行关联,其结构如
图8 两层钢框架结构示意图
Fig.8 Schematic diagram of two-story steel frame structure
基于Matlab/Simulin
为了验证基于Matlab/Simulink所建立的有限元模型的可靠性,在OpenSEES中对该框架进行建模:材料本构采用uniaxialMaterial ElasticPP,梁、柱均采用dispBeamColumn单元,倚靠柱采用elasticBeamColumn单元,其余模型参数与Matlab/Simulink模型保持一致。Matlab/Simulink模型与OpenSEES模型前5阶周期对比如
阶数 | 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 |
从
将El Centro NS, 1940作为结构的地震动输入,采用CR、Chang及本文方法进行时程分析,同时采用Δt=0.001s的CAA算法作为参考解。
图9 首层位移时程曲线
Fig.9 Time history curves of displacement at first story
从
时间步长Δt | 算法 | NEE | NRMSE |
---|---|---|---|
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 |
算法 | 时间步长Δt | ||
---|---|---|---|
0.005 | 0.01 | 0.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。
由
基于模型的积分算法非常重要的特性是无条件稳定性和显式表达式。无条件稳定性意味着无需为了保证算法的稳定性而可以选择较大的时间步长,显式格式则意味着求解非线性问题时无需迭代。因此,基于模型的积分算法计算效率高也对应于上述两个方面:一是较大的时间步长,从而计算步数较少;二是每一积分时间步的计算时间短,因为无需迭代。本节将选取
图10 质量-弹簧系统示意图
Fig.10 Schematic diagram of Mass-Spring system
在
(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算法和本文方法的计算结果。
自由度数 | 方法 | 计算时间 | x50计算误差/ % | |
---|---|---|---|---|
NEE | NRMSE | |||
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 |
从
本文提出了一种新型高精度半显式基于模型的积分算法。新算法与半显式Chang算法具有相同的速度、位移差分方程,采用二阶Pade近似极点映射方法生成新算法的积分参数。对新算法的精度、稳定性和周期延长和振幅衰减等数值特性进行分析,发现本文方法具有四阶精度,周期延长极小,并且没有能量耗散。通过3个具有代表性的数值算例,进一步论证了本文方法的优越性。
作者贡献声明
傅博:论文修改,数据分析;
张付泰:论文撰写,编程计算;
陈瑾:论文修改,数据校核。
参考文献
KOLAY C, RICLES 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 Engineering, 2016, 107(1):49. [百度学术]
彭天波, 谢馨, 曾忠, 等. 采用Chang方法的混合试验的稳定性和精度[J]. 同济大学学报(自然科学版), 2014, 42(12): 1790. [百度学术]
PENG Tianbo, XIE Xin, ZENG Zhong, et al. Stability and accuracy of sharking table-actuator hybrid test with Chang method [J]. Journal of Tongji University (Natural Science), 2014, 42(12): 1790. [百度学术]
KOLAY C, RICLES J M. Force-based frame element implementation for real-time hybrid simulation using explicit direct integration algorithms [J]. Journal of Structural Engineering, 2018, 144(2): 04017191. [百度学术]
傅博, 蒋欢军. 具有高稳定性的剪切型子结构振动台试验[J]. 同济大学学报(自然科学版), 2017, 45(11): 1602. [百度学术]
FU Bo, JIANG Huanjun. Shear-type substructure sharking table testing method with high stability [J]. Journal of Tongji University (Natural Science), 2017, 45(11): 1602. [百度学术]
FENG D C, KOLAY C, RICLES J M, et al. Collapse simulation of reinforced concrete frame structures [J]. The Structural Design of Tall and Special Buildings, 2016, 25(12): 578. [百度学术]
CHANG S Y. Explicit pseudodynamic algorithm with unconditional stability [J]. Journal of Engineering Mechanics, 2002, 128(9): 935. [百度学术]
CHEN C, RICLES J M. Development of direct integration algorithms for structural dynamics using discrete control theory [J]. Journal of Engineering Mechanics, 2008, 134(8): 676. [百度学术]
GUI Y, WANG J T, JIN F, et al. Development of a family of explicit algorithms for structural dynamics with unconditional stability [J]. Nonlinear Dynamics, 2014, 77(4): 1157. [百度学术]
FU B, FENG D C, JIANG H. A new family of explicit model-based integration algorithms for structural dynamic analysis [J], International Journal of Structural Stability and Dynamics, 19(6): 1950053. [百度学术]
NEWMARK N M. A method of computation for structural dynamics [J]. Journal of the Engineering Mechanics 1959, 85(3): 67―94. [百度学术]
HILBER H M. Analysis and design of numerical integration methods in structural dynamics [R]. Berkeley: Earthquake Engineering Research Center, 1976. [百度学术]
CHOPRA A K. Dynamics of structures: theory and applications to earthquake engineering [M]. 4th ed. Upper Saddle River: Prentice Hall, 2011. [百度学术]
The MathWorks Inc. Matlab R2014b [EB/OL].[2021-08-29] http://www.mathworks.com. [百度学术]