一种结构稳态非线性动力学响应快速计算方法及系统和存储介质与流程

本发明涉及结构动力分析技术领域，特别是一种结构稳态非线性动力学响应快速计算方法及系统和存储介质。

背景技术：

非线性动力学计算方法是了解和掌握结构非线性动力学特征的基本手段。最为基础的非线性动力学计算方法是以runge-kutta法和中心差分法等为代表的时域算法。可是，时域算法往往仅能计算得出稳定解，却无法获得不稳定解。对于非线性系统而言，非稳定解的确定对于明确结构响应特征具有重要的意义。因此，频域算法便表现出了其优越性。常用的频域算法主要包括多尺度展开法，谐波平衡法等。

多尺度展开法是现阶段非线性动力学领域最为常用也最为经典的算法之一。其主要特点是在对系统进行求解的过程中，可以明确系统的一些非线性动力学基本特征，从而掌握系统发生某些非线性动力响应的原因。但是多尺度展开法由于其推导过程复杂，一般至多适用于2-3个自由度系统的非线性动力学分析。因此对于复杂结构，无法使用此类方法。1981年lau和cheung首次提出了增量谐波平衡法(ihb)，该方法将响应解的形式表达为傅里叶级数展开式，并通过伽辽金积分过程求解傅里叶级数的系数。相较多尺度展开法，该方法计算效率有了大幅的提升。

在这此之后，增量谐波平衡法在不断地发展，有学者在其计算效率方面不断的开展工作。2015年wang和zhu结合快速傅里叶变换(fft)和broyden方法大幅提高了ihb方法的计算效率。然而，尽管如此对于需要上百个自由度表示的复杂系统的非线性动力学问题，现有方法仍然很难开展分析研究。但是，对于很大部分的复杂结构而言，当其发生大幅振动时，往往只有很少部分构件进入非线性工作状态，而其余的大多数构件仍保持线性工作状态。例如，在基础下设置隔震支座的隔震结构，此类结构地震作用下，由于隔震支座的存在，使得上部结构一般仅会发生较小的振动，结构仍处于弹性状态，而隔震支座本身却因承担了大幅的振动变形，进入非线性振动状态。

针对这类情况，尽管已经明确的知道了大部分构件并不需要开展非线性分析，但由于部分构件已经进入非线性状态，现有的计算方法仍然只能通过对整个结构开展非线性动力学分析来掌握结构的动力响应，导致计算效率很低，或无法开展有效的结构计算。

技术实现要素：

有鉴于此，本发明的目的在于提供一种结构稳态非线性动力学响应快速计算方法及系统和存储介质，该方法涉及非线性动力学算法，通过将目标结构分解为线性和非线性子结构并利用增量谐波平衡法实现快速计算方法。

为达到上述目的，本发明提供如下技术方案：

本发明提供的结构稳态非线性动力学响应快速计算方法，包括以下步骤：

获取目标结构及目标结构的初始参数；

将目标结构沿线性构件与非线性构件的交界面分解为线性子结构和非线性子结构；

根据交界面的界面力计算非线性子结构的响应；

计算整体结构响应的解。

进一步，所述线性子结构和非线性子结构的分解是按照以下公式进行的：

fb,l＝-fb,nl

式中，

ql和qnl分别表示为和

进一步，所述非线性子结构的响应是利用弧长法按照以下公式对结构响应进行跟踪求解：

式中，表示等效线性刚度；表示等激励；表示残差；δa表示响应的幅值增量；δω表示外荷载的频率增量。

进一步，所述整体结构响应的解按照以下步骤计算：

计算线性子结构的界面力；

按照以下公式得到界面处个自由度响应的第k个谐波项：

gk＝[a1kcos(kτ),a2kcos(kτ),…,ankcos(kτ),b1ksin(kτ),b2ksin(kτ),…,bnksin(kτ)]^t,

式中：下标n代表界面处所有自由度的数量；

gk表示界面处各自由度第k个谐波项的位移响应时程；ank表示界面处各自由度正弦响应的幅值；bnk表示界面处各自由度余弦响应的幅值；kτ表示相位角；

按照以下公式求出各谐波项qi,l的解hk(k＝1,2,…,m)：

式中：hk表示线性子结构内部自由度第k个谐波项的位移响应幅值；

将解hk带入下式得线性子结构的界面力fb,l的第k个谐波项yk；

式中：yk表示线性子结构的界面力fb,l的第k个谐波项；

循环重复计算求出所有谐波项，得到作用于线性子结构上的界面力fb,l；

通过线性子结构上的界面力fb,l按照以下公式计算残差δf：

通过迭代方法不断减小δf的范数直至小于预设量ξ，则得到整体结构响应的解。

进一步，所述增量δf计算按照以下步骤进行：

将δf整理为以下增量形式：

求解所述增量方程得到δf的伽辽金积分

判断残差值是否小于预设量ξ，如果是则得到整体结构的第i+1组响应的q⁽ⁱ⁺¹⁾,ω⁽ⁱ⁺¹⁾。

本发明提供的结构稳态非线性动力学响应快速计算系统，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现以下步骤：

获取目标结构及目标结构的初始参数；

将目标结构沿线性构件与非线性构件的交界面分解为线性子结构和非线性子结构；

根据交界面的界面力计算非线性子结构的响应；

计算整体结构响应的解。

进一步，所述线性子结构和非线性子结构的分解是按照以下公式进行的：

fb,l＝-fb,nl；

式中，ql和qnl分别表示为和

进一步，所述非线性子结构的响应是利用弧长法按照以下公式对结构响应进行跟踪求解：

式中，表示等效线性刚度；表示等激励；表示残差；δa表示响应的幅值增量；δω表示外荷载的频率增量。

本发明提供的一种存储介质，其上存储有计算机程序，该程序被处理器执行时实现权利要求1-5任一项所述方法的步骤。

本发明的有益效果在于：

本发明提供的结构稳态非线性动力学响应快速计算方法，是针对部分构件进入非线性工作状态的结构体系的非线性动力响应的快速计算方法。通过将目标结构分解为线性和非线性子结构并利用增量谐波平衡法形成一套计算方法，由于仅有非线性子结构需要开展非线性分析，因此计算效率较现有方法有大幅提升。

本实施例所述方法所得结构的稳态响应与runge-kutta法完全吻合，验证了方法的准确性。本实施例所述方法进行一次迭代所需时间为8秒钟。即所述方法计算速度为ihb方法的60倍。说明了本算法较高的计算效率。

本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书来实现和获得。

附图说明

为了使本发明的目的、技术方案和有益效果更加清楚，本发明提供如下附图进行说明：

图1为线性子结构与非线性子结构分解示意图。

图2为带有隔震支座的平面混凝土框架结构示意图。

图3为结构的分解示意图。

图4为结构幅-频响应曲线示意图。

图5为本方法与runge-kutta法对比示意图。

图6为结构稳态非线性动力学响应快速计算方法流程图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步说明，以使本领域的技术人员可以更好的理解本发明并能予以实施，但所举实施例不作为对本发明的限定。

实施例1

如图6所示，本实施例提供的结构稳态非线性动力学响应快速计算方法，是针对部分构件进入非线性工作状态的结构体系的非线性动力响应的快速计算方法。通过将目标结构分解为线性和非线性子结构并利用增量谐波平衡法形成一套计算方法，由于仅有非线性子结构需要开展非线性分析，因此计算效率较现有方法有大幅提升。

本实施例提供的结构稳态非线性动力学响应快速计算方法，包括以下步骤：

步骤1：获取目标结构及目标结构的初始参数如结构的刚度，质量，阻尼，以及隔震支座的动力参数，将目标结构沿线性构件与非线性构件的交界面分解为线性子结构和非线性子结构，如果目标结构可明确哪些构件会进入非线性状态工作，则可将其在两类构件的分界面处分离为线性部分和非线性部分，通过分别在两部分的边界处施加界面力来代替其相互作用。

如图1所示，图1为线性子结构与非线性子结构分解示意图，将目标结构的动力控制方程改写为经过解耦的两个控制方程：

式中：

ml，cl，kl是线性子结构的质量、阻尼和刚度矩阵；

mnl，cnl表示非线性子结构的质量和阻尼矩阵；fr表示非线性系统自身的恢复力；

fl表示作用于线性系统上的外荷载的幅值；ll表示反映界面力位置的映射矩阵；

fnl表示作用于非线性系统上的外荷载的幅值；lnl表示反映界面力位置的映射矩阵；

fb,l表示不同子系统间的界面处作用于线性子系统的界面力；

fb,nl表示不同子系统间的界面处作用于非线性子系统的界面力；

ql，代表线性子结构的位移，速度以及加速度响应；

qnl，代表非线性子结构的位移，速度以及加速度响应；

τ表示相位角；

本实施例中字母中带有“nl”下标的变量的意义与“l”下标的意义区别在于其对应于非线性子结构；其中，nl表示非线性子结构，l表示线性子结构，对线性子结构和非线性子结构的边界处施加相互作用的界面力；由作用力与反作用力定理可知：

fb,l＝-fb,nl。(3)

将式(1-2)中的ql和qnl分别表示为和

其中，qi,l表示线性子结构的内部节点位移向量；qb,l表示线性子结构的界面节点位移向量；

qi,nl表示非线性子结构的内部节点位移向量；qb,nl表示非线性子结构的界面节点位移向量；

下标i和b分别位于子结构内部和界面处的自由度。

步骤2：初步计算非线性子结构的响应

由于本实施例利用弧长法对结构响应进行跟踪求解，因此在寻找新的一组结构响应时，须已知一组已有的结构响应。

设q⁽ⁱ⁾,ω⁽ⁱ⁾为目前已知的结构的第i组响应，现寻找结构的第i+1组响应。根据式(1)及q⁽ⁱ⁾,ω⁽ⁱ⁾可得结构的截面力fb,l；将fb,l带入式(2)，将解的形式表示为增量的形式：

q＝q⁽ⁱ⁾+δq，ω＝ω⁽ⁱ⁾+δω(4)

其中，

q⁽ⁱ⁾表示第i组已知的结构位移响应；ω⁽ⁱ⁾表示第i组已知的结构位移响应所对应的外荷载激励频率；

q表示待求的结构位移响应；δq表示位移响应的增量；

ω表示待求的结构位移响应所对应的外荷载激励频率；δω表示激励频率的增量；

将式(4)带入式(2)可得：

式中，ω0＝ω⁽ⁱ⁾表示第i组已知的结构位移响应所对应的外荷载激励频率；mnl表示非线性子结构的质量矩阵；cnl表示非线性子结构的阻尼矩阵；δq表示激励频率的增量；r为计算残差；f表示作用于非线性子结构的外荷载；fr(q⁽ⁱ⁾)表示非线性子结构的恢复力；

τ＝ωt，变量右上方的撇代表对τ求导；

令响应表示为傅里叶级数的形式，如下所示：

q⁽ⁱ⁾＝sa,δq＝sδa(7)

式中，

其中，ai是非线性子结构第i个自由度所有频率响应成分的幅值所组成的向量；

p表示非线性子结构自由度的数量；

cs＝[1cos(τ)cos(2τ)…cos((m-1)τ)sin(τ)sin(2τ)…sin((m-1)τ)](8)

式中，

m代表响应中所考虑的频率成分的数量；

将式(8)带入式(5)并对所有变量在τ在一个周期内[0,2π]进行积分，得到下式：

式中：

表示质量矩阵的伽辽金积分；表示阻尼矩阵的伽辽金积分；表示刚度矩阵的伽辽金积分；表示残差；s^t表示s矩阵的转置

对式(9)开展弧长法求解，便可初步获得非线性子结构第i+1组响应的“暂定解”。

步骤3：计算整体结构响应的解

尽管在计算时暂时将线性子结构与非线性子结构分离，但实际上在界面处具有完全相同的响应，即

qb,l＝qb,nl.(14)

将式(1)改写为：

式中：

mii，mib表，mbi，mbb分别表示ml的四个子矩阵；

cii，cib，cbi，cbb分别表示cl的四个子矩阵；

kii，kib，kbi，kbb分别表示kl的四个子矩阵；

表示线性子结构内部各自由度的加速度响应；表示线性子结构界面各自由度的加速度响应；表示线性子结构内部各自由度的速度响应；表示线性子结构界面各自由度的速度响应；qi,l表示线性子结构内部各自由度的位移响应；qb,l表示线性子结构界面各自由度的位移响应；fii表示作用于线性子结构内部各自由度的外荷载；fbb表示作用于线性子结构界面各自由度的外荷载。

由式(14)，可知fb,l(qb,l)＝fb,l(qb,nl)。根据式(7)所示响应解的形式，可进一步计算线性子结构的界面力。

将式(15-16)中所有的qb,l由q'b,nl替换，界面处个自由度响应的第k个谐波项可表示为

gk＝[a1kcos(kτ),a2kcos(kτ),…,ankcos(kτ),b1ksin(kτ),b2ksin(kτ),…,bnksin(kτ)]^t,(17)

式中：n代表界面处所有自由度的数量；

式中：下标n代表界面处所有自由度的数量；

gk表示界面处各自由度第k个谐波项的位移响应时程；ank表示界面处各自由度正弦响应的幅值；bnk表示界面处各自由度余弦响应的幅值；kτ表示相位角；

由此，根据式(15)可分别求出qi,l各谐波项的解hk(k＝1,2,…,m)，如下式所示：

式中：hk表示线性子结构内部自由度第k个谐波项的位移响应幅值；

q表示线性子结构内部自由度的等效刚度；p表示线性子结构界面自由度的等效刚度；

将hk带入式(16)，可得线性子结构的界面力fb,l的第k个谐波项yk，如下所示

式中：yk表示线性子结构的界面力fb,l的第k个谐波项；

v表示线性子结构内部自由度的计算刚度；u表示线性子结构界面自由度的计算刚度；

重复式(18-23)将所有谐波项均求出后，便可获得fb,l。

比较分别有线性和非线性子系统所得到的界面力

δfb＝fb,l(qb,l)+fb,nl(qnl)(24)

式中，

δfb表示线性子结构与非线性子结构界面力的残差；

fb,l(qb,l)表示线性子结构的界面力；

fb,nl(qnl)表示非线性子结构的界面力；

由式(1)可知

将式(25)带入式(24)可得

理论上只有当δf＝0时，才意味着两个子系统的响应完全协调，即系统得到了真正的响应。因此，需要通过迭代的方法不断减小δf的范数直至该值小于某一小量ξ。

将式(26)整理为增量形式，并进行伽辽金积分

δaj＝[01×(j-1),δaj,01×(m-j)]^t,(32)

式中：表示非线性子结构的质量矩阵的伽辽金积分；表示非线性子结构的阻尼矩阵的伽辽金积分；表示非线性子结构的刚度矩阵的伽辽金积分；表示残差；表示线性子结构对非线性子结构所提供的等效刚度；表示线性子结构界面处的等效刚度；

表示线性子结构界面处的等效刚度矩阵的第j列；ab表示界面位移向量；

δaj表示假想界面位移的增量；01×(m-j)表示0向量；δaj为一接近于0的正数；

式(27)可再次借助弧长法求解，直至小于ξ，则所得结果即为整体结构的第i+1组响应的q⁽ⁱ⁺¹⁾,ω⁽ⁱ⁺¹⁾。

实施例2

下面通过具体的实例对本发明以及本发明效果进行进一步的说明。

如图2所示，图2为带有各镇支座的平面混凝土框架结构示意图，以某五层三开间的钢筋混凝土结构的一榀框架的非线性动力计算为例，其柱底安装隔震支座；混凝土材料的密度和弹性模量分别按2.5×10³kg/m³and3.0×10⁴mpa计算。框架结构的层高为3m，开间尺寸为6m。柱子和梁的横截面尺寸分别设为0.45m×0.45m和0.3m×0.6m。

本例中每层的柱沿高度被划分为2个单元，同时每跨的梁被划分为3个单元。对上部结构所有构件均采用伯努利梁单元进行模拟。因此，上部结构共由85个单元组成，包括74个节点。对于隔震支座，则选用bouc-wen作为描述其本构关系的数学模型，如式(33-34)所示

fb＝k1x+k2z,(33)

式中：

fb为隔震支座所提供的反力；

x是隔震支座的剪切变形；

z为归一化的滞回力；

k1支座恢复力线性分量的刚度；

k2恢复力非线性分量的模量系数；

β和γ为控制滞回环尺度和形状的参数；

r为反应滞回曲线光滑度的正整数；

kbl＝k1+k2被称为支座的线性化参数。

本实施例中设β＝γ＝50。同时令kbl＝1.0×10³kn/m(其中k1＝100kn/m,k2＝900kn/m)，即支座的线性化刚度为上部柱子刚度的1/45。

由图2可知，当结构承受基底加速度激励后，其上部结构由于隔震支座的存在并不会发生大幅振动，即处于线性工作状态，而各镇支座则会承担大部分结构的位移，导致其非线性工作。由此，以柱底为界面将结构分解为线性子结构和非线性子结构两部分，如图3所示，其中，图3中(3a)为线性子结构；(3b)非线性子结构，图3为结构的分解。

由于每一个隔震支座仅有一个水平自由度，因此非线性子结构共包含4个自由度，而线性子结构包含214个自由度。

设结构承受幅值为0.33m/s²的单频(谐波)基底加速度激励，激励频率由0.25hz逐渐变化至1.5hz(覆盖了结构的第一阶自振频率)。现计算结构在每一个激励频率下的稳态响应。

则结构的动力控制方程可写为

式中：mf表示上部结构的质量矩阵；cf表示上部结构的阻尼矩阵；kf表示上部结构的刚度矩阵；表示隔震支座提供的恢复力；

式中，fn表示；

k1a，k1b，k1c，k1d表示支座a、b、c、d的恢复力线性分量的刚度；

k2a，k2b，k2c，k2d表示支座a、b、c、d的恢复力非线性分量的模量系数；

βa，βb，βc，βd，γa，γb，γc，γd表示支座a、b、c、d的控制滞回环尺度和形状的参数；

表示第1-4自由度的速度响应；

q5，q6，q7，q8表示第5-8自由度的位移响应；

i为单位矩阵；k1x＝100kn/m，k2x＝900kn/m，其中，x＝a,b,c,d；

依据式(13)：

式中：

cs＝[1cos(τ)cos(2τ)…cos((m-1)τ)sin(τ)sin(2τ)…sin((m-1)τ)],

c′s＝[0-sin(τ)-2sin(2τ)…-(m-1)sin((m-1)τ)cos(τ)2cos(2τ)…(m-1)cos((m-1)τ)].

式中，cs表示位移响应的基本谐波成分向量；f′n1表示支座非线性恢复力对求导雅可比矩阵；f′n2表示支座非线性恢复力对q5，q6，q7，q8求导的雅可比矩阵；m表示谐波频率成分的个数；表示非线性子结构的位移响应幅值向量；w表示非线性子结构的恢复力矩阵。

本实施例中将m设为12。依据本实施例中所述方法可得结构屋顶及柱底的幅-频响应曲线如图4所示，图4中的(4a)为第1谐波项，(4b)为第3谐波项。图4表示结构幅频响应曲线。为验证算法的准确性，将本实施例所述算法与经典的时程算法runge-kutta法进行对比。选取激励频率为0.77hz(0.67ω1)和0.32hz，分别通过runge-kutta法计算结构响应。所得屋顶及柱底的时程响应曲线，如图5所示，图5为所提方法与runge-kutta法对比，其中，(5a)为屋顶响应(ω＝4.83rad/s)，(5b)为柱底响应(ω＝4.83rad/s)，(5c)为屋顶响应(ω＝2.01rad/s)，(5d)为柱底响应(ω＝2.01rad/s)，可以看出，本实施例所述方法所得结构的稳态响应与runge-kutta法完全吻合，验证了方法的准确性。

本实施例所述方法与经典的ihb方法对本实施例中的结构开展非线性动力学计算。发现如果令m＝12，即计算结果相对更加准确的情况下，采用ihb方法时由于计算所需时间过长，基本无法看展计算并进行对比。因此，选取m＝6对两种算法进行对比。ihb方法计算一次迭代所需时间为8分钟左右，而本实施例所述方法进行一次迭代所需时间为8秒钟。即所述方法计算速度为ihb方法的60倍。说明了本算法较高的计算效率。

以上所述实施例仅是为充分说明本发明而所举的较佳的实施例，本发明的保护范围不限于此。本技术领域的技术人员在本发明基础上所作的等同替代或变换，均在本发明的保护范围之内。本发明的保护范围以权利要求书为准。