一种高频耗散可控的结构动力响应数值计算方法与流程

本发明属于土木工程中的结构振动分析技术领域，具体涉及一种高频耗散可控的结构动力响应数值计算方法。

背景技术：

土木工程结构在地震、风等的时程激励下会产生振动响应，当激励足够大时，可能引起结构的破坏。中国位于环太平洋地震带和欧亚地震带之间，长期受到地震的困扰，给我国带来巨大损失。研究土木工程结构在地震、风等激励作用下的动力响应具有十分重要的意义。

研究结构动力响应的试验方法包括拟静力试验、拟动力试验、振动台试验和实时混合试验。结构动力响应数值计算方法为了解结构的动态性能及进行结构动力学设计提供了强有力的手段。随着拟动力试验和实时混合试验的快速发展，数值算法的精度、稳定性和计算效率成为试验成功开展的关键因素之一。

结构动力响应数值计算方法包括隐式算法和显式算法两类，隐式算法通常为无条件稳定，但需要大量迭代计算，效率不高；显式算法通常为条件稳定，虽不需迭代计算，但时间步长的选择受到限制，对于复杂非线性结构，为满足算法稳定性可能要选择一个较小的时间步长，增加计算量，导致计算效率降低。由于有限元离散不能准确表达结构的高阶振型，会产生虚假的高频响应，影响数值计算结果的准确性。目前，大多数商业有限元使用的仍是newmark(纽马克-β)法，不能过滤高阶虚假振型。

技术实现要素：

本发明提出一种高频耗散可控的结构动力响应数值计算方法，该方法采用离散控制理论推导算法参数和极点，通过选择合适的系数控制算法的精度、非线性稳定界限，也可以控制算法的高频耗散特性，最后对结构的动力响应进行分析，为获得结构动力响应提供了新的、有效的技术手段，有效提高了计算精度和效率，解决了中心差分法和newmark法不能有效过滤高阶虚假振型的问题。

本发明采用以下的技术方案实现：

一种高频耗散可控的结构动力响应数值计算方法，包括以下步骤：

步骤一：输入初始结构特性，建立运动微分方程，选定位移速度递推式的参数，并确定参数表达式；设置时间步长δt和数值计算的总时间tmax；

步骤二：根据初始条件x(0)和进行初始加速度向量计算；

步骤三：数值计算，根据位移速度递推式，计算状态变量在t＝ti+1＝ti+δt时刻的值xi+1、和

步骤四：判断是否进行终止数值积分计算，

若t＜tmax，则令i＝i+1，t＝ti+1＝ti+δt，并转至步骤三，进行下一时刻的数值计算；

若t≥tmax，则转至步骤五；

步骤五：输出结构动力响应计算结果。

优选的，步骤一中，根据初始结构特性，建立的运动微分方程为：

式中，m，c，k为初始结构特性，分别为结构初始时刻的质量、阻尼和刚度矩阵；x，分别为位移向量、速度向量和加速度向量，称为状态变量；f为外部激励，如地震加速度时程；

位移速度递推式为：

式中：δt为时间步长；xi和xi+1为结构位移向量在t＝ti和t＝ti+1时刻的近似值；和为结构速度向量在t＝ti和t＝ti+1时刻的近似值；为结构加速度向量在t＝ti时刻的近似值；α和β为参数；

优选的，结构为单自由度运动体系时，位移速度递推式的参数求解包括以下具体步骤：

(1)建立结构动力系统的传递函数：

式中：x(z)为系统输出结构位移的z变换，f(z)为系统输入外部激励的z变换；n2、n1和n0，及d2、d1和d0为系数，表达式如下：

n2＝0(5)

n1＝αδt²(6)

n0＝(β-α)δt²(7)

d2＝m(8)

d1＝m(αω²+2βξω-2)(9)

d0＝m[(β-α)ω²-2βξω+1](10)

式中：ω＝ω·δt，ω为结构自振频率，ξ为阻尼比；

(2)令z1,2为特征方程d2z²+d1z+d0＝0的根，由式(8)～(10)则得到：

(3)算法耗散特性由谱半径ρ决定，若结构具有高频耗散能力，则需满足：

ω→0,ρ＝1(12)

ω→∞,ρ＝ρ∞(13)

式中：ρ∞为ω取值无穷大时对应的谱半径值，为区间[0,1]内的常数；

(4)令ξ＝0，将式(11)代入式(13)，则得到：

式中：a为ω的一次函数，为引入结构特性和ω＝0时参数分母有意义，a＝a(ξω+1)，a为系数；

联立式(11)和式(14)即得参数表达式：

由以上参数表达式可得到算法谱半径示意图，如附图2所示，当ω→0时，谱半径ρ＝1，当ω→∞时，谱半径ρ＝ρ∞，满足式(12)和式(13)的要求。

优选的，结构为多自由度运动体系时，参数表达式为：

β＝a·(δt²·k+a/2δt·c+a·m)^-1·m(18)

参数确定后，根据式(2)和式(3)即可进行后续步骤的数值计算。

优选的，步骤二中，所述的x(0)是指结构初始时刻的位移向量，所述的是指结构初始时刻的速度向量，根据初始条件x(0)和进行初始数值计算是指计算结构初始加速度向量：土木工程结构的位移x(0)和速度通常为0，根据式(1)得到结构初始加速度向量计算公式为：

优选的，步骤五中，所述的结构动力响应计算结果依据式(1)～(3)计算而得，包括位移时程曲线、速度时程曲线和加速度时程曲线。

本发明一种高频耗散可控的结构动力响应数值计算方法，具有以下性质：

1)显式算法；

2)对线性结构和非线性刚度软化结构无条件稳定，对非线性刚度硬化结构为条件稳定，稳定界限为

3)通过系数a的取值控制算法精度，最佳精度取值为a＝11；

4)通过系数ρ∞的取值控制算法数值耗散能力，ρ∞＝1时算法无数值耗散，ρ∞＝0时算法数值耗散能力最强。

与现有技术相比，本发明的有益效果为：

1)本发明利用离散控制理论，以算法特性为条件推导参数，最终获得算法极点，是一种新的结构动力响应数值计算方法，该新的数值计算方法为分析结构在地震、风等作用下的动力响应提供了新的技术手段；

2)在计算结构动力响应的时候，可通过参数取值控制算法的精度、稳定性，以及控制数值耗散能力，计算过程无需迭代计算，计算效率高，能过滤高阶虚假振型对计算结果的影响；与传统的中心差分法和newmark法相比，不仅能有效过滤高阶虚假振型，而且具有更高的精度和更高的稳定性。

附图说明

图1为本发明参数推导流程图；

图2为本发明谱半径示意图；

图3为本发明实施例双层剪切框架结构；

图4为中心差分法得到的底层位移时程曲线；

图5为中心差分法得到的底层速度时程曲线；

图6为newmark线加速度法得到的底层位移时程曲线；

图7为newmark线加速度法得到的底层速度时程曲线；

图8为本发明实施例方法得到的底层位移时程曲线；

图9为本发明实施例方法得到的底层速度时程曲线；

图10为本发明的算法步骤流程框图。

具体实施方式

为更加清楚地理解本发明所提出的一种高频耗散可控的结构动力响应数值计算方法，下面结合附图及实施例对本发明做进一步说明。

本实施例给出一种高频耗散可控的结构动力响应数值计算方法，采用离散控制理论获得算法传递函数，以高频耗散特性为条件推导算法参数和极点，得到算法最终位移和速度递推式，为获得结构动力响应提供了新的、有效的技术手段，有效提高了计算精度和效率。通过位移速度递推式中的系数取值，可控制算法的精度和刚度硬化系统稳定界限，以及高频耗散能力，并且算法为显式算法，不需迭代计算，计算效率高，具体实施分析过程如下：

以图3所示的双层剪切结构的动力响应计算为例，本发明的具体实施步骤归纳如下：

根据结构初始特性建立运动微分方程：

式中，m，c，k分别为结构的质量、阻尼和刚度矩阵；x，分别为位移向量、速度向量和加速度向量，称为状态变量；f为外部激励，如地震加速度时程；

在所述的数值计算方法中，位移速度递推式为：

式中：δt为时间步长；xi和xi+1为结构位移向量在t＝ti和t＝ti+1时刻的近似值；和为结构速度向量在t＝ti和t＝ti+1时刻的近似值；为结构加速度向量在t＝ti时刻的近似值；α和β为参数；

设置时间步长δt和数值计算的总时间tmax，选定系数；

进一步，数值计算方法中位移速度递推式的参数求解包括以下步骤：

建立结构动力系统的传递函数：

式中：x(z)为系统输出结构位移的z变换，f(z)为系统输入外部激励的z变换；n2、n1和n0，及d2、d1和d0为系数：

n2＝0(5)

n1＝αδt²(6)

n0＝(β-α)δt²(7)

d2＝m(8)

d1＝m(αω²+2βξω-2)(9)

d0＝m[(β-α)ω²-2βξω+1](10)

令z1,2为特征方程d2z²+d1z+d0＝0的根，由式(8)～(10)则得到：

算法耗散特性由谱半径ρ决定，若结构具有高频耗散能力，则需满足：

ω→0,ρ＝1(12)

ω→∞,ρ＝ρ∞(13)

式中：ρ∞为区间[0,1]内的常数；.

令ξ＝0，将式(11)代入式(13)，则有：

式中：a为ω的一次函数，为引入结构特性和ω＝0时参数分母有意义，a＝a(ξω+1)，a为系数；

联立式(11)和式(14)即得参数表达式：

对应多自由度结构为：

β＝a·(δt²·k+a/2δt·c+a·m)^-1·m(18)

参数确定后，根据式(2)和式(3)即可进行后续步骤的数值计算。

根据初始条件计算结构初始加速度向量：

初始条件为结构初始时刻的位移x(0)和速度代入式(1)则得到初始时刻的加速度向量：

实施例双层剪切框架结构各物理量均为国际单位，质量为：第一层m1＝1，第二层m2＝1；刚度为：第一层k1＝10⁴，第二层k2＝1；结构阻尼比为0；外荷载作用于第一层，f1＝10000sin(1.2t)，f2＝0；初始位移和速度均为0；

建立结构运动微分方程：

选择时间步长δt＝0.02s，数值计算总时间中心差分法为tmax＝5.0s，newmark线加速度法和本发明实施例算法为tmax＝20.0s，系数a＝11，ρ∞＝0；

利用本发明的式(17)和式(18)，获得算法参数：

α＝12·(δt²·k+11/2δt·c+11·m)^-1·m

β＝11·(δt²·k+11/2δt·c+11·m)^-1·m

根据式(19)得到初始加速度为：

使用本发明式(1)～(3)即可获得实施例双层剪切框架结构各层的位移、速度的时程响应，如图8和图9所示。

图4～图7分别是中心差分法和newmark线加速度法获得的实施例双层剪切框架结构底层位移时程曲线和速度时程曲线。从图4和图5中可以看出，中心差分法不能得到位移和速度时程的正确解，计算结果随着时间增长逐渐放大，存在误差逐步累积的现象，有失稳的趋势；图6和图7的结果表明，newmark线加速度法不能过滤高阶振型对结构动力响应的影响，高阶振型对结构位移、速度时程曲线的计算结果影响较大，与图8和图9相比，本发明提出的高频耗散可控的数值计算方法可以很好地过滤高阶振型，得到正确结果，且具有很好的稳定性。

本发明的重点在于利用离散控制理论获得传递函数和特征方程，以算法高频耗散能力为条件，推导算法参数和极点，最终算法既能有效过滤结构高频响应，又具有较高的精度和稳定性。