非线性结构动力学系统瞬态响应无条件稳定时间积分方法

阅读：评论：0

1.本发明属于非线性结构动力学系统响应分析领域，尤其涉及一种包含非线性几何项、非线性阻尼项及其耦合项的结构动力学系统的具有二阶精度、可调数值耗散以及无条件稳定性的时间积分方法。
2.背景方法
3.随着新技术、新材料、新结构等前沿成果在航空航天、土木工程、机器人等领域的广泛应用，由此带来的材料非线性、几何非线性和接触非线性因素给结构动力学分析带来了严峻挑战。目前，最为常用的结构动力学分析方案是“空间离散”与“时间离散”相结合，其也被广泛应用在商业软件cae(computer aided engineering)中。
4.灵活有效的有限单元法(finite element method,fem)是目前较为流行的空间离散工具，其作用是将由偏微分方程控制的无穷自由度连续系统转变为由常微分方程控制的有限自由度离散系统。完成空间离散后，下一个环节就是时间离散，基于差分思想建立的时间积分方法(time integration method,tim)是一种强大的时间离散工具。其求解思路是将待求的时间域离散成一系列连续的时间区间，人为假设状态变量在一个时间区间内的变化规律，从上一时刻已知状态变量信息求解当前时刻未知信息。
5.newmark家族方法是时间积分方法中的经典工作，其成员包括著名的梯形法则(trapezoidal rule,tr)和中心差分方法(central difference method,cdm)等。对于线性系统，梯形法则是无条件稳定的，但对于简单的非线性问题如刚性单摆转动问题，梯形法则出现了发散现象。由此引发当前时间积分方法领域的热点和难点问题：如何提高时间积分方法在求解非线性结构动力学系统时的稳定性，以及能否为非线性结构动力学系统设计出无条件稳定的时间积分方法。为了解决上述问题，学者们发展了耗散型方法、保能量方法以及结构依赖型方法。
6.大量的数值结果表明通过引入数值耗散可以有效地推迟时间积分方法在非线性结构动力学系统求解中失稳的时刻。由于newmark方法的耗散格式仅具有一阶精度，因此学者们通过加权newmark方法的平衡方程构造出了同时具有二阶精度和数值耗散的参数类方法，比如广义-alpha方法、wood-bossak-zienkiewicz-alpha方法、hiber-hughes-taylor-alpha方法等。但是，由于平衡方程在时间结点上没有得到满足，导致该类方法的加速度只具有一阶精度。为了设计出具有较高低频精度的耗散型时间积分方法，线性多步法和复合方法被提出，这两类方法均有效地提高了精度。但是，线性多步法和复合方法并不能从根本上改变方法的稳定性。总之，耗散型方法仍然存在失稳问题，并且其低频精度比非耗散型方法低。
7.为了能够给非线性结构动力学系统提供稳定的数值仿真结果，基于能量有界准则的保能量方法被提出。能量有界准则由belytschko提出，其定义是：当前时刻的动能和势能之和小于等于上一时刻的动能、势能和外力功之和。因此，保能量方法是通过强制系统能量守恒来实现无条件稳定性。受构造方法的限制，现有的保能量方法大都只能处理含有非线性几何项的结构动力学问题，比如krenk方法和emm等，而只有少部分如ecm能够处理非线性
几何项和阻尼项独立存在的非线性结构动力学问题。目前，能够处理含有非线性几何项和阻尼项耦合的结构动力学系统的保能量方法尚未被构造出来。综上，保能量方法虽然在求解非线性结构动力学系统时不会发散，但其适用对象有限，且具有额外计算量(如离散能量函数修正运算等)致使计算效率降低。
8.对大型复杂结构进行动力学响应分析时，精度、效率和稳定性三者之间往往难以平衡。为了解决这个问题，2002年chang提出了结构依赖型时间积分方法，这类方法的算法参数与结构初始特性和时间步长大小紧密相关。为了进一步提高chang方法的数值性能，包括精度、稳定性和耗散特性，一系列性能更加优异的结构依赖型时间积分方法被提出，比如kr方法、cr方法、du-yang-zhao方法等。理论分析和数值实验结果表明，该类方法对线性系统是无条件稳定的，同时对非线性软刚度系统具有理想的稳定性。目前，受构造方式限制，这类方法局限于具有对角质量矩阵、正定对称刚度和阻尼矩阵的线性系统和具有幂形式非线性刚度的结构动力学系统。总之，结构依赖型方法虽然有效地平衡了精度、效率和稳定性，但其在求解非线性结构动力学系统时仍有失稳的可能，且适用对象有限。
9.因此，针对包含非线性几何项、非线性阻尼项及其耦合项的结构动力学系统，设计可高精度、高效率、无条件稳定地分析瞬态响应的时间积分方法是现阶段亟待突破的研究内容。其将为非线性结构动力学问题提供性能优越的新解法，且可以促进时间积分方法朝着同时改进稳定性、精度和效率方向发展，服务于国家对复杂结构动力学分析高性能算法的需求，具有重要的理论意义和应用价值。

技术实现要素：

10.本发明公开的非线性结构动力学系统瞬态响应无条件稳定时间积分方法主要目的是：利用有限单元方法建立包含非线性几何项、非线性阻尼项以及它们的耦合项的结构动力学方程；结合多分步技术将一个时间单元离散成三个分步，并设计各分步状态变量(包括位移、速度和加速度)的更新公式；基于bn稳定性理论、传递因子分析方法、局部截断误差分析方法，以效率最大化、二阶精度、可控耗散和无条件稳定性为目标设计分步公式中的参数；在此基础上，建立三分步时间积分方法在包含非线性几何项、非线性阻尼项以及它们的耦合项在结构动力学系统中的求解流程，实现对包含非线性几何项、非线性阻尼项以及它们的耦合项的结构动力学系统瞬态响应的高效、精确且稳定地时域分析。
11.本发明的目的是通过如下技术方法实现。
12.本发明公开的非线性结构动力学系统瞬态响应无条件稳定时间积分方法，通过利用有限单元方法建立包含非线性几何项、非线性阻尼项以及它们耦合项的结构动力学方程；将一个时间单元划分成三个分步，基于差分格式建立状态变量在一个时间单元内的演化规律，依次使用n点后向插值格式(n＝1,2,3)，建立便于数值性能调控的多参数优化时间步进方程，由于调控参数增多，从而提高参数调节的可控性和精度。所述数值性能包括效率、精度、耗散、稳定性。基于构建的多参数优化时间步进方程，推导便于数值性能调控的三个分步的有效刚度矩阵，以效率最大化为目标建立具有统一系数的参数约束关系，利用建立参数约束关系将所述三个分步的有效刚度矩阵变换为具有统一形式的有效刚度矩阵，提高非线性系统的结构动力学响应时域分析效率。基于构建的多参数优化时间步进方程，针对包含非线性几何项、非线性阻尼项以及它们的耦合项的结构动力学系统，以二阶精度为
目标建立局部截断误差为o(δt3)的参数约束关系，以数值耗散可控为目标建立高频段传递因子可调的参数约束关系，以无条件稳定性为目标建立满足bn稳定性条件的参数约束关系，从而提高非线性系统的结构动力学响应时域分析精度、耗散性和稳定性。将得到的参数约束关系代入建立的多参数优化时间步进方程，再将所述多参数优化时间步进方程代入构建的包含非线性几何项、非线性阻尼项及其耦合项的结构动力学控制方程，实现无条件稳定地分析包含非线性几何项、非线性阻尼项及其耦合项的结构动力学系统的瞬态响应。
13.本发明公开的非线性结构系统瞬态响应无条件稳定时间积分方法，包括以下步骤：
14.步骤1：基于有限单元方法fem建立空间离散的包含非线性几何项、非线性阻尼项及其耦合项的结构动力学控制方程。
15.基于有限单元方法fem(finite element method)建立空间离散的包含非线性几何项、非线性阻尼项及其耦合项的结构动力学控制方程，如下所示：
[0016][0017]
式中，m为质量矩阵；为包含非线性内力和外部激励项的向量；t为时间。对于线性系统，向量退化为其中k为刚度矩阵；c为阻尼矩阵；q(t)为外部激励向量。
[0018]
步骤2：基于差分格式建立状态变量在一个时间单元内的演化规律，将一个时间单元划分成三个分步，依次使用n点后向插值格式(n＝1,2,3)，建立便于数值性能调控的多参数优化时间步进方程，由于调控参数增多，提高参数调节的可控性和精度。所述数值性能包括效率、精度、耗散、稳定性。
[0019]
基于差分格式，建立状态变量，包括位移、速度和加速度，在一个时间单元[t,t+δt]内随时间t的演化规律。将一个时间单元[t,t+δt]划分成三个分步，即[t,t+c1δt]、[t+c1δt,t+c2δt]以及[t+c2δt,t+c3δt]。第一个分步[t,t+c1δt]的步进公式为
[0020][0021]
其中，x
t+c1δt
、和代表t+c1δt时刻的位移、速度和加速度；x
t
和代表t时刻的位移和速度；δt代表时间步长大小；c1为第一个分步引入的算法参数。第二个分步[t+c1δt,t+c2δt]的步进公式为
[0022][0023]
其中，x
t+c2δt
、和代表t+c2δt时刻的位移、速度和加速度；c2和α为第二个分步引入的算法参数。第三个分步[t+c2δt,t+c3δt]的步进公式为
[0024][0025]
其中，x
t+c3δt
、和代表t+c3δt时刻的位移、速度和加速度；c3、β和η为第三个分步引入的算法参数。时间单元终点时刻t+δt信息通过对三个已知时刻信息t+c1δt、t+c2δt以及t+c3δt进行加权得到，对应的步进方程为
[0026][0027]
其中，b1和b2为待定加权参数。
[0028]
步骤3：基于步骤2构建的多参数优化时间步进方程，推导便于数值性能调控的三个分步的有效刚度矩阵，以效率最大化为目标建立具有统一系数的参数约束关系，利用建立参数约束关系将所述三个分步的有效刚度矩阵变换为具有统一形式的有效刚度矩阵，提高非线性系统的结构动力学响应时域分析效率。根据具有统一形式的有效刚度矩阵得到统一系数的参数约束关系。
[0029]
当时，整理三个分步的步进方程得
[0030][0031]
其中，和分别代表有效刚度矩阵和有效载荷向量。三个分步的有效刚度矩阵的表达式为
[0032][0033][0034][0035]
从公式(7)至(9)中发现，当参数满足公式(10)中的关系时，三个分步的有效刚度矩阵完全相同。此时，三个分步使用同一个有效刚度矩阵，进而降低计算量。
[0036]
c1＝c2α＝c3η
ꢀꢀꢀ
(10)
[0037]
另外，要求时间点t+c1δt和t+c3δt关于时间单元中点t+0.5δt对称，且对应的加权参数相同。
[0038]
由此可得，
[0039]
c1+c3＝1
ꢀꢀꢀ
(11)
[0040]
b1＝1-b
1-b2ꢀꢀꢀ
(12)
[0041]
步骤4：基于步骤2构建的多参数优化时间步进方程，针对包含非线性几何项、非线性阻尼项以及它们耦合项的结构动力学系统，以二阶精度为目标建立局部截断误差为o(δt3)的参数约束关系，以数值耗散可控为目标建立高频段传递因子可调的参数约束关系，以无条件稳定性为目标建立满足bn稳定性条件的参数约束关系，从而提高非线性系统的结构动力学响应时域分析的精度、耗散性和稳定性。
[0042]
步骤4.1：针对包含非线性几何项、非线性阻尼项以及它们耦合项的结构动力学系统，推导便于数值性能调控的一个时间单元内的butcher表。首先，将公式(2)至(5)等价地写为基于一阶方程的形式，如下所示
[0043][0044][0045][0046][0047]
式中，基于butcher表，参数为
[0048][0049]
步骤4.2：在步骤4.1的基础上，建立的基于一阶方程的多参数优化时间步进方程推导其在一个时间单元内的传递方程，如下所示
[0050]zt+δt
＝a(τ)z
t
,τ＝λδt
ꢀꢀꢀ
(18)
[0051]
其中，λ代表特征根。式中，传递因子a的表达式为
[0052][0053]
步骤4.3：在步骤4.1的基础上，基于局部截断误差确定多参数优化时间步进方程的精度阶次，局部截断误差的定义为
[0054]
σ＝a(τ)-a
exact
(τ)，a
exact
(τ)＝exp(τ)
ꢀꢀꢀ
(20)
[0055]
根据二阶精度，要求
[0056]
a(0)＝a
(1)
(0)＝a
(2)
(0)＝1
ꢀꢀꢀ
(21)
[0057]
根据式(21)求解得到如下满足二阶精度要求的时间步进方程的参数关系，即
[0058][0059]
步骤4.4：在步骤4.1的基础上，以数值耗散可控为目标建立高频段传递因子可调的参数约束关系，如下所示，
[0060][0061]
对于非耗散格式ρ
∞
＝1，根据式(23)建立如下参数关系，
[0062][0063]
对于耗散格式0≤ρ
∞
《1，根据式(23)建立如下参数关系
[0064][0065]
步骤4.5：在步骤4.1的基础上，以无条件稳定性为目标建立满足bn稳定性条件的参数约束关系。根据bn稳定性定义，建立如下不等式约束条件，即
[0066]
s1＝b1(2c
1-b1)≥0
ꢀꢀꢀ
(26)
[0067]
s2＝b2(2c2α-b2)≥0
ꢀꢀꢀ
(27)
[0068]
s3＝(1-b
1-b2)[2c3η-(1-b
1-b2)]≥0
ꢀꢀꢀ
(28)
[0069]
s5＝s1s
3-(1-b
1-b2)2[c3(1-β-η)-b1]2≥0
ꢀꢀꢀ
(29)
[0070]
s6＝s2s
3-(1-b
1-b2)2(c3β-b2)2≥0
ꢀꢀꢀ
(30)
[0071][0072]
上述不等式约束条件，即公式(26)至公式(31)，能够保证在任意时间步长下稳定地计算包含非线性几何项、非线性阻尼项以及它们耦合项的结构动力学系统。根据公式(26)至(31)，可建立β和c3之间的关系。对于非耗散格式ρ
∞
＝1，根据公式(26)至公式(31)，建立如下参数关系，
[0073][0074]
对于耗散格式0≤ρ
∞
《1，根据公式(26)—公式(31)，建立如下参数关系，
[0075][0076]
式中，c3＝f(ρ
∞
)为关于高频段耗散量的已知函数。
[0077]
步骤5：根据步骤3得到具有统一系数的参数约束关系，根据步骤4得到满足二阶精度、可控数值耗散和无条件稳定性的参数约束关系，将步骤3和步骤4得到的参数约束关系代入步骤2建立的多参数优化时间步进方程，将所述多参数优化时间步进方程结合newton迭代技术，代入步骤1构建的包含非线性几何项、非线性阻尼项及其耦合项的结构动力学控制方程，实现高效、精确、且无条件稳定地计算包含非线性几何项、非线性阻尼项及其耦合项的结构动力学系统瞬态响应。
[0078]
还包括步骤6：根据步骤1至步骤5得到的适用于包含非线性几何项、非线性阻尼项及其耦合项的结构动力学系统的无条件稳定时间积分方法，能够在结构动力学工程应用中广泛应用于非线性结构动力学系统的瞬态响应分析，改善并预测非线性结构动力学结构和
性能，解决相关工程技术问题。
[0079]
有益效果：
[0080]
1、本发明公开的一种非线性结构动力学系统瞬态响应的无条件稳定时间积分方法，基于差分格式建立状态变量在一个时间单元内的演化规律，将一个时间单元划分成三个分步，依次使用n点后向插值格式(n＝1,2,3)，建立便于数值性能(包括效率、精度、耗散性和稳定性)调控的多参数优化时间步进方程，由于调控参数增多，提高参数调节的可控性和精度。
[0081]
2、本发明公开的一种非线性结构动力学系统瞬态响应的无条件稳定时间积分方法，基于步骤2构建的多参数优化时间步进方程，推导便于数值性能调控的三个分步的有效刚度矩阵，以效率最大化为目标建立具有统一系数的参数约束关系，利用建立参数约束关系将所述三个分步的有效刚度矩阵变换为具有统一形式的有效刚度矩阵，提高非线性结构动力学系统瞬态响应的时域分析效率。
[0082]
3、本发明公开的一种非线性结构动力学系统瞬态响应的无条件稳定时间积分方法，基于步骤2构建的多参数优化时间步进方程，针对包含非线性几何项、非线性阻尼项以及它们耦合项的非线性结构动力学系统，以二阶精度为目标建立满足局部截断误差为o(δt3)的参数约束条件；以耗散可控为目标建立高频段传递因子可调的参数约束关系；以无条件稳定性为目标建立满足bn稳定性条件的参数约束关系，提高包含非线性几何项、非线性阻尼项以及它们耦合项的非线性结构动力学系统瞬态响应的时域分析的精度、耗散性和稳定性。
[0083]
4、本发明公开的一种非线性结构动力学系统瞬态响应的无条件稳定时间积分方法，相较于已有的耗散型方法、保能量方法以及结构依赖型方法，计算量小、精度高、稳定性强，方便使用。即使在缺乏专业知识背景的情况下也能进行操作，能够在结构动力学工程应用中广泛应用于非线性结构动力学系统的瞬态响应分析，改善并预测非线性结构动力学结构和性能，解决相关工程技术问题。
附图说明：
[0084]
图1为本发明公开的非线性结构动力学系统瞬态响应无条件稳定时间积分方法的计算流程图；
[0085]
图2为考虑非线性格林应变的悬臂梁冲击模型；
[0086]
图3为广义-alpha方法、ρ
∞-bathe方法和本发明方法在悬臂梁自由端沿z方向位移和速度随时间演化规律(初始速度为20m/s)；
[0087]
图4为广义-alpha方法、ρ
∞-bathe方法和本发明方法在悬臂梁自由端沿z方向位移和速度随时间演化规律(初始速度为40m/s)；
[0088]
图5为弹塑性悬臂板冲击模型；
[0089]
图6为广义-alpha方法、ρ
∞-bathe方法和本发明方法在悬臂板自由端沿x方向位移随时间演化规律；
[0090]
图7为广义-alpha方法、ρ
∞-bathe方法和本发明方法在悬臂板自由端沿z方向位移随时间演化规律。
具体实施方式
[0091]
为了更好地说明本发明的目的和优点，下面结合附图和实例对发明内容做进一步说明。其中，自始至终相同或类似的符号表示相同或类似功能。
[0092]
实施例1：
[0093]
如图1所示，本实施例公开的非线性结构动力学系统瞬态响应无条件稳定时间积分方法，应用于预测考虑非线性格林应变的悬臂梁冲击模型的动力学响应，具体实现步骤如下：
[0094]
步骤1：以考虑非线性格林应变的悬臂梁冲击模型为例，如图2所示，来说明本发明方法的稳定性和效率优势。基于有限单元法，利用10个非线性双节点梁单元对考虑非线性格林应变悬臂梁冲击模型进行空间离散，进而可得如下非线性结构动力学方程
[0095][0096]
其中，v0为在自由端施加的初始速度。
[0097]
步骤2：在步骤1的基础上，将待求时间域[0,0.5]离散成n个连续的时间单元，即[t,t+δt]、[t+δt,t+2δt]、
…
、[t+(n-1)δt,t+nδt]，其中δt代表一个时间单元的尺寸。基于差分格式，建立状态变量，包括位移、速度和加速度，在一个是时间单元[t,t+δt]内随时间t的演化规律。将一个时间单元[t,t+δt]划分成三个分步，即[t,t+c1δt]、[t+c1δt,t+c2δt]以及[t+c2δt,t+c3δt]。
[0098]
第一个分步[t,t+c1δt]的步进公式为
[0099][0100]
第二个分步[t+c1δt,t+c2δt]的步进公式为
[0101][0102]
第三个分步[t+c2δt,t+c3δt]的步进公式为
[0103][0104]
时间单元终点时刻t+δt信息通过对三个已知时刻信息t+c1δt、t+c2δt以及t+c3δt进行加权得到，对应的步进方程为
[0105][0106]
步骤3：根据具有统一形式的有效刚度矩阵得到统一系数的参数约束关系，如下所示
[0107]
c1＝c2α＝c3η
ꢀꢀꢀ
(39)
[0108]
另外，要求时间点t+c1δt和t+c3δt关于时间单元中点t+0.5δt对称，且对应的加权参数相同。由此可得，
[0109]
c1+c3＝1
ꢀꢀꢀ
(40)
[0110]
b1＝1-b
1-b2ꢀꢀꢀ
(41)
[0111]
步骤4：基于步骤2构建的多参数优化时间步进方程，推导便于数值性能调控的一个时间单元内的butcher表，如下所示，
[0112][0113]
推导便于数值性能调控的一个时间单元内的传递因子公式，如下所示
[0114][0115]
基于局部截断误差确定多参数优化时间步进方程的精度阶次，局部截断误差的定义为
[0116]
σ＝a(τ)-a
exact
(τ)，a
exact
(τ)＝exp(τ)
ꢀꢀꢀ
(44)
[0117]
根据二阶精度，要求
[0118]
a(0)＝a
(1)
(0)＝a
(2)
(0)＝1
ꢀꢀꢀ
(45)
[0119]
根据式(21)求解得到如下满足二阶精度要求的时间步进方程的参数关系，即
[0120][0121]
以耗散可控为目标建立高频段传递因子可调的参数约束关系，如下所示
[0122][0123]
对于非耗散格式ρ
∞
＝1，根据式(47)可建立如下参数关系，
[0124][0125]
对于耗散格式0≤ρ
∞
《1，根据式(47)可建立如下参数关系
[0126][0127]
以无条件稳定性为目标建立满足bn稳定性条件的参数约束关系，如下所示
[0128]
s1＝b1(2c
1-b1)≥0
ꢀꢀꢀ
(50)
[0129]
s2＝b2(2c2α-b2)≥0
ꢀꢀꢀ
(51)
[0130]
s3＝(1-b
1-b2)[2c3η-(1-b
1-b2)]≥0
ꢀꢀꢀ
(52)
[0131]
s5＝s1s
3-(1-b
1-b2)2[c3(1-β-η)-b1]2≥0
ꢀꢀꢀ
(53)
[0132]
s6＝s2s
3-(1-b
1-b2)2(c3β-b2)2≥0
ꢀꢀꢀ
(54)
[0133][0134]
根据公式(51)—(55)，可建立β和c3之间的关系。对于非耗散格式ρ
∞
＝1，整理可得
[0135][0136]
对于耗散格式0≤ρ
∞
《1，整理可得
[0137][0138]
式中，c3＝f(ρ
∞
)为关于高频段耗散量的已知函数。
[0139]
考虑非线性格林应变的悬臂梁冲击模型是保守系统。因此，根据动力学特性，非耗散格式ρ
∞
＝1被采用，由公式(39)、(40)、(41)、(46)、(48)和(56)可确定其余参数为
[0140][0141]
步骤5：将步骤3和步骤4得到的参数约束关系代入步骤2建立的多参数优化时间步进方程，结合newton迭代技术和所述多参数优化时间步进方程由初始时刻响应依次递推求解步骤1建立的考虑非线性格林应变的悬臂梁冲击模型在每个时间离散点处的状态变量信息。
[0142]
在第一个分步[t,t+c1δt]中，首先假定
[0143][0144]
其中，上标“0”表示迭代次数。将公式(59)代入第一个分步的步进方程(35)，可得
[0145][0146]
若小于等于给定的容许误差，则进入下一个分步，否则进行迭代修正直到满足容许误差要求。根据公式(60)中的第一个动力学平衡方程可建立第一个分步的修正增量即
[0147][0148]
进而可得修正的加速度为
[0149][0150]
将公式(62)代入公式(60)即可得到修正的位移和速度满足容许误差要求的第一个分步的响应作为已知量用于第二个分步的计算。
[0151]
在第二个分步[t+c1δt,t+c2δt]中，首先假定
[0152][0153]
其中，上标“0”表示迭代次数。将公式(63)代入第二个分步的步进方程(36)，可得
[0154][0155]
若小于等于给定的容许误差，则进入下一个分步，否则进行迭代修正直到满足容许误差要求。根据公式(64)中的第一个动力学平衡方程可建立第二个分步的修正增量即
[0156][0157]
进而可得修正的加速度为
[0158][0159]
将公式(66)代入公式(64)即可得到修正的位移和速度满足容许误差要求的第二个分步的响应作为已知量用于第三个分步的计算。
[0160]
在第三个分步[t+c2δt,t+c3δt]中，首先假定
[0161][0162]
其中，上标“0”表示迭代次数。将公式(67)代入第三个分步的步进方程(37)，可得
[0163][0164]
若小于等于给定的容许误差，则进入下一个分步，否则进行迭代修正直到满足容许误差要求。根据公式(68)中的第一个动力学平衡方程可建立第三个分步的修正增量即
[0165][0166]
进而可得修正的加速度为
[0167][0168]
将公式(70)代入公式(68)即可得到修正的位移和速度将三个分步满足容许误差的状态变量值代入公式(38)可得时间单元终点的响应。
[0169]
图3和图4分别表示v0＝20m/s和40m/s时自由端沿y方向的位移和速度，随着v0的增大，该动力学模型能够激发出更多的高频模式从而导致时间积分方法更容易失稳。由数值模拟结果可以得出：本发明公开的非线性结构动力学系统瞬态响应无条件稳定时间积分方法在整个仿真过程中都给出了稳定且精确的结果，而常用的ρ
∞-bathe方法在早期仿真阶段得到了发散的结果。此外，表1将已有方法和本发明方法的计算效率进行了比较，可以发现本发明方法的计算成本最低。由实施例1可以说明，本发明公开的非线性结构动力学系统瞬态响应无条件稳定时间积分方法比传统方法具有更强的稳定性和更高的计算效率。
[0170]
表1广义-alpha方法、ρ
∞-bathe方法和本发明方法的cpu和迭代次数(容许误差ε＝10e-12)
[0171][0172][0173]
实施例2：
[0174]
如图1所示，本实施例公开的非线性结构动力学系统瞬态响应无条件稳定时间积分方法，应用于预测弹塑性悬臂板在集中冲击力作用下的动力学响应，具体实现步骤如下：
[0175]
步骤1：以弹塑性悬臂板为例，如图5所示，来说明本发明方法的精度优势。基于有限单元法，利用1296个考虑j2双线性次弹塑性材料板单元对悬臂板模型进行空间离散，进而可得如下非线性结构动力学方程
[0176][0177]
步骤2：在步骤1的基础上，将待求时间域[0,1]离散成n个连续的时间单元，即[t,t+δt]、[t+δt,t+2δt]、
…
、[t+(n-1)δt,t+nδt]，其中δt代表一个时间单元的尺寸。基于差分格式，建立状态变量，包括位移、速度和加速度，在一个是时间单元[t,t+δt]内随时间t的演化规律。将一个时间单元[t,t+δt]划分成三个分步，即[t,t+c1δt]、[t+c1δt,t+c2δt]以及[t+c2δt,t+c3δt]。
[0178]
第一个分步[t,t+c1δt]的步进公式为
[0179][0180]
第二个分步[t+c1δt,t+c2δt]的步进公式为
[0181][0182]
第三个分步[t+c2δt,t+c3δt]的步进公式为
[0183][0184]
时间单元终点时刻t+δt信息通过对三个已知时刻信息t+c1δt、t+c2δt以及t+c3δt进行加权得到，对应的步进方程为
[0185][0186]
步骤3：根据具有统一形式的有效刚度矩阵得到统一系数的参数约束关系，如下所示
[0187]
c1＝c2α＝c3η
ꢀꢀꢀ
(76)
[0188]
另外，要求时间点t+c1δt和t+c3δt关于时间单元中点t+0.5δt对称，且对应的加权参数相同。由此可得，
[0189]
c1+c3＝1
ꢀꢀꢀ
(77)
[0190]
b1＝1-b
1-b2ꢀꢀꢀ
(78)
[0191]
步骤4：基于步骤2构建的多参数优化时间步进方程，推导便于数值性能调控的一个时间单元内的butcher表，如下所示，
[0192][0193]
推导便于数值性能调控的一个时间单元内的传递因子公式，如下所示
[0194][0195]
基于局部截断误差确定多参数优化时间步进方程的精度阶次，局部截断误差的定义为
[0196]
σ＝a(τ)-a
exact
(τ)，a
exact
(τ)＝exp(τ)
ꢀꢀꢀ
(81)
[0197]
根据二阶精度，要求
[0198]
a(0)＝a
(1)
(0)＝a
(2)
(0)＝1
ꢀꢀꢀ
(82)
[0199]
根据式(82)求解得到如下满足二阶精度要求的时间步进方程的参数关系，即
[0200][0201]
以耗散可控为优化目标建立高频段传递因子可调的参数约束关系，如下所示
[0202][0203]
对于非耗散格式ρ
∞
＝1，根据式(84)可建立如下参数关系，
[0204][0205]
对于耗散格式0≤ρ
∞
《1，根据式(84)可建立如下参数关系
[0206][0207]
以无条件稳定性为优化目标建立满足bn稳定性条件的参数约束关系，如下所示
[0208]
s1＝b1(2c
1-b1)≥0
ꢀꢀꢀ
(87)
[0209]
s2＝b2(2c2α-b2)≥0
ꢀꢀꢀ
(88)
[0210]
s3＝(1-b
1-b2)[2c3η-(1-b
1-b2)]≥0
ꢀꢀꢀ
(89)
[0211]
s5＝s1s
3-(1-b
1-b2)2[c3(1-β-η)-b1]2≥0
ꢀꢀꢀ
(90)
[0212]
s6＝s2s
3-(1-b
1-b2)2(c3β-b2)2≥0
ꢀꢀꢀ
(91)
[0213][0214]
根据公式(87)—(92)，可建立β和c3之间的关系。对于非耗散格式ρ
∞
＝1，整理可得
[0215][0216]
对于耗散格式0≤ρ
∞
《1，整理可得
[0217][0218]
式中，c3＝f(ρ
∞
)为关于高频段耗散量的已知函数。
[0219]
考虑非线性格林应变的悬臂梁冲击模型是保守系统。因此，根据动力学特性，非耗散格式ρ
∞
＝1被采用，由公式(76)、(77)、(78)、(83)、(85)和(93)可确定其余参数为
[0220][0221]
步骤5：将步骤3和步骤4得到的参数约束关系代入步骤2建立的多参数优化时间步进方程，结合newton迭代技术和所述多参数优化时间步进方程由初始时刻响应依次递推求解步骤1建立的考虑非线性格林应变的悬臂梁冲击模型在每个时间离散点处的状态变量信息。
[0222]
在第一个分步[t,t+c1δt]中，首先假定
[0223][0224]
其中，上标“0”表示迭代次数。将公式(96)代入第一个分步的步进方程(72)，可得
[0225][0226]
若小于等于给定的容许误差，则进入下一个分步，否则进行迭代修正直到满足容许误差要求。根据公式(97)中的第一个动力学平衡方程可建立第一个分步的修正增量即
[0227][0228]
进而可得修正的加速度为
[0229][0230]
将公式(99)代入公式(97)即可得到修正的位移和速度满足容许误差要求的第一个分步的响应作为已知量用于第二个分步的计算。
[0231]
在第二个分步[t+c1δt,t+c2δt]中，首先假定
[0232][0233]
其中，上标“0”表示迭代次数。将公式(100)代入第二个分步的步进方程(73)，可得
[0234][0235]
若小于等于给定的容许误差，则进入下一个分步，否则进行迭代修正直到满足容许误差要求。根据公式(101)中的第一个动力学平衡方程可建立第二个分步的修正增量即
[0236][0237]
进而可得修正的加速度为
[0238][0239]
将公式(103)代入公式(101)即可得到修正的位移和速度满足容许误差要求的第二个分步的响应作为已知量用于第三个分步的计算。
[0240]
在第三个分步[t+c2δt,t+c3δt]中，首先假定
[0241][0242]
其中，上标“0”表示迭代次数。将公式(104)代入第三个分步的步进方程(74)，可得
[0243][0244]
若小于等于给定的容许误差，则进入下一个分步，否则进行迭代修正直到满足容许误差要求。根据公式(105)中的第一个动力学平衡方程可建立第三个分步的修正增量即
[0245][0246]
进而可得修正的加速度为
[0247][0248]
将公式(107)代入公式(105)即可得到修正的位移和速度将三个分步满足容许误差的状态变量值代入公式(75)可得时间单元终点的响应。
[0249]
图6和图7分别绘制了a点(自由端中点)在x和z方向的位移，由数值模拟结果可以得出：本发明公开的非线性结构动力学系统瞬态响应无条件稳定时间积分方法在整个仿真过程中都给出了最高的计算精度。由实施例2表明，本发明公开的非线性结构动力学系统瞬态响应无条件稳定时间积分方法比传统方法具有更高的计算精度。
[0250]
以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

技术特征：

1.非线性结构动力学系统瞬态响应无条件稳定时间积分方法，其特征在于：包括以下步骤，步骤1：基于有限单元方法fem建立空间离散的包含非线性几何项、非线性阻尼项及其耦合项的结构动力学控制方程；步骤2：基于差分格式建立状态变量在一个时间单元内的演化规律，将一个时间单元划分成三个分步，依次使用n点后向插值格式(n＝1,2,3)，建立便于数值性能调控的多参数优化时间步进方程，由于调控参数增多，提高参数调节的可控性和精度；所述数值性能包括效率、精度、耗散、稳定性；步骤3：基于步骤2构建的多参数优化时间步进方程，推导便于数值性能调控的三个分步的有效刚度矩阵，以效率最大化为目标建立具有统一系数的参数约束关系，利用建立参数约束关系将所述三个分步的有效刚度矩阵变换为具有统一形式的有效刚度矩阵，提高非线性系统的结构动力学响应时域分析效率；根据具有统一形式的有效刚度矩阵得到统一系数的参数约束关系；步骤4：基于步骤2构建的多参数优化时间步进方程，针对包含非线性几何项、非线性阻尼项以及它们耦合项的结构动力学系统，以二阶精度为目标建立局部截断误差为o(δt3)的参数约束关系，以数值耗散可控为目标建立高频段传递因子可调的参数约束关系，以无条件稳定性为目标建立满足bn稳定性条件的参数约束关系，从而提高非线性系统的结构动力学响应时域分析的精度、耗散性和稳定性；步骤5：根据步骤3得到具有统一系数的参数约束关系，根据步骤4得到满足二阶精度、可控数值耗散和无条件稳定性的参数约束关系，将步骤3和步骤4得到的参数约束关系代入步骤2建立的多参数优化时间步进方程，将所述多参数优化时间步进方程结合newton迭代技术，代入步骤1构建的包含非线性几何项、非线性阻尼项及其耦合项的结构动力学控制方程，实现高效、精确、且无条件稳定地计算包含非线性几何项、非线性阻尼项及其耦合项的结构动力学系统瞬态响应。2.如权利要求1所述的非线性结构动力学系统瞬态响应无条件稳定时间积分方法，其特征在于：步骤1实现方法为，基于有限单元方法fem建立空间离散的包含非线性几何项、非线性阻尼项及其耦合项的结构动力学控制方程，如下所示：式中，m为质量矩阵；为包含非线性内力和外部激励项的向量；t为时间；对于线性系统，向量退化为其中k为刚度矩阵；c为阻尼矩阵；q(t)为外部激励向量。3.如权利要求2所述的非线性结构动力学系统瞬态响应无条件稳定时间积分方法，其特征在于：步骤2实现方法为，基于差分格式，建立状态变量，包括位移、速度和加速度，在一个是时间单元[t,t+δt]内随时间t的演化规律；将一个时间单元[t,t+δt]划分成三个分步，即[t,t+c1δt]、[t+c1δt,t+c2δt]以及[t+c2δt,t+c3δt]；第一个分步[t,t+c1δt]的步进公式为
其中，x
t+c1δt
、和代表t+c1δt时刻的位移、速度和加速度；x
t
和代表t时刻的位移和速度；δt代表时间步长大小；c1为第一个分步引入的算法参数；第二个分步[t+c1δt,t+c2δt]的步进公式为其中，x
t+c2δt
、和代表t+c2δt时刻的位移、速度和加速度；c2和α为第二个分步引入的算法参数；第三个分步[t+c2δt,t+c3δt]的步进公式为其中，x
t+c3δt
、和代表t+c3δt时刻的位移、速度和加速度；c3、β和η为第三个分步引入的算法参数；时间单元终点时刻t+δt信息通过对三个已知时刻信息t+c1δt、t+c2δt以及t+c3δt进行加权得到，对应的步进方程为其中，b1和b2为待定加权参数。4.如权利要求3所述的非线性结构动力学系统瞬态响应无条件稳定时间积分方法，其特征在于：步骤3实现方法为，当时，整理三个分步的步进方程得其中，和分别代表有效刚度矩阵和有效载荷向量；三个分步的有效刚度矩阵的表达式为
从公式(7)至(9)中发现，当参数满足公式(10)中的关系时，三个分步的有效刚度矩阵完全相同；此时，三个分步使用同一个有效刚度矩阵，进而降低计算量；c1＝c2α＝c3η
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(10)另外，要求时间点t+c1δt和t+c3δt关于时间单元中点t+0.5δt对称，且对应的加权参数相同；由此可得，c1+c3＝1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(11)b1＝1-b
1-b2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)5.如权利要求4所述的非线性结构动力学系统瞬态响应无条件稳定时间积分方法，其特征在于：步骤4实现方法为，步骤4.1：针对包含非线性几何项、非线性阻尼项以及它们耦合项的结构动力学系统，推导便于数值性能调控的一个时间单元内的butcher表；首先，将公式(2)至(5)等价地写为基于一阶方程的形式，如下所示基于一阶方程的形式，如下所示基于一阶方程的形式，如下所示基于一阶方程的形式，如下所示式中，基于butcher表，参数为步骤4.2：在步骤4.1的基础上，建立的基于一阶方程的多参数优化时间步进方程推导其在一个时间单元内的传递方程，如下所示z
t+δt
＝a(τ)z
t
,τ＝λδt
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(18)其中，λ代表特征根；式中，传递因子a的表达式为步骤4.3：在步骤4.1的基础上，基于局部截断误差确定多参数优化时间步进方程的精度阶次，局部截断误差的定义为σ＝a(τ)-a
exact
(τ)，a
exact
(τ)＝exp(τ)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(20)
根据二阶精度，要求a(0)＝a
(1)
(0)＝a
(2)
(0)＝1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(21)根据式(21)求解得到如下满足二阶精度要求的时间步进方程的参数关系，即步骤4.4：在步骤4.1的基础上，以数值耗散可控为目标建立高频段传递因子可调的参数约束关系，如下所示，对于非耗散格式ρ
∞
＝1，根据式(23)建立如下参数关系，对于耗散格式0≤ρ
∞
<1，根据式(23)建立如下参数关系步骤4.5：在步骤4.1的基础上，以无条件稳定性为目标建立满足bn稳定性条件的参数约束关系；根据bn稳定性定义，建立如下不等式约束条件，即s1＝b1(2c
1-b1)≥0
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(26)s2＝b2(2c2α-b2)≥0
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(27)s3＝(1-b
1-b2)[2c3η-(1-b
1-b2)]≥0
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(28)s5＝s1s
3-(1-b
1-b2)2[c3(1-β-η)-b1]2≥0
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(29)s6＝s2s
3-(1-b
1-b2)2(c3β-b2)2≥0
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(30)上述不等式约束条件，即公式(26)至公式(31)，能够保证在任意时间步长下稳定地计算包含非线性几何项、非线性阻尼项以及它们耦合项的结构动力学系统；根据公式(26)至(31)，可建立β和c3之间的关系；对于非耗散格式ρ
∞
＝1，根据公式(26)至公式(31)，建立如下参数关系，对于耗散格式0≤ρ
∞
<1，根据公式(26)—公式(31)，建立如下参数关系，式中，c3＝f(ρ
∞
)为关于高频段耗散量的已知函数。

技术总结

本发明公开的非线性结构动力学系统瞬态响应无条件稳定时间积分方法，属于非线性结构动力学系统响应分析领域。本发明实现方法如下：建立包含非线性几何项、非线性阻尼项以及它们的耦合项的结构动力学方程；结合多分步技术将时间单元离散成三个分步，并设计各分步状态变量更新公式；基于BN稳定性理论、传递因子分析方法、局部截断误差分析方法，以效率最大化、二阶精度、可控耗散和无条件稳定性为目标设计分步公式中的参数；建立三分步时间积分方法在包含非线性几何项、非线性阻尼项以及它们的耦合项在结构动力学系统中的求解流程，实现对包含非线性几何项、非线性阻尼项以及它们的耦合项的结构动力学系统瞬态响应的高效、精确且稳定地时域分析。且稳定地时域分析。且稳定地时域分析。