一种工程结构在极端荷载下动态响应的复合时域分析方法与流程

1.本发明涉及结构动力学技术领域，特别是涉及一种工程结构在极端荷载下动态响应的复合时域分析方法及系统。


背景技术：

2.随着计算机技术的迅猛发展，时间积分方法已成为求解结构动态响应的重要方法，并大量应用于商业有限元软件中。时间积分法也称逐步积分法，是假定已知时刻t＝0的位移u0、速度并将时间求解时间域平均分成n个时间段0、δt、2δt、
…
、nδt(nδt为总时长)，且已求得0、δt、2δt、
…
和t的解，建立求解t δt时刻解的递推算法。总体来说，时间积分法可以分成两类：显式时间积分法和隐式时间积分法。显式方法通常是条件稳定的，所需的时间步长远小于隐式方法。然而，显式方法比隐式方法花费更少的时间来计算矩阵。显式方法由于其无需进行方程组的迭代求解，在求解大型结构或非线性程度较高的问题时更具优势，因此受到很多外学者的广泛关注和研究。
3.在现有的显式方法中，中心差分法(cd法)是一种广泛应用于结构动力学问题的经典的单步显式计算方法，并已应用到ansys、abaqus和adina等著名商用软件中。虽然，目前的cd法通过引入人工阻尼可以有效改进算法的数值耗散特性，但效果仍有待提升。数值耗散的目的是滤除或有效减少离散后的动力系统的虚假高频模态。实际上，显式方法要在高频段获得理想的数值耗散而不在低频段引入较大的耗散是不容易的。noh和bathe提出了一种具有两个子步的复合显式方法(称为noh
‑
bathe法，noh g,bathe k
‑
j.an explicit time integration scheme for the analysis of wave propagations.computers&structures 2013)，这种显式方法采用了子步的策略，可以引入更多的自由参数以控制算法的基本性能，并且能有效控制算法的数值耗散特性，该方法已集成在商用软件adina中。尽管如此，noh
‑
bathe法对于无阻尼情况可以获得二阶算法精度，对于有阻尼情况的加速度仅获得一阶精度，因而其算法精度仍有待提升，尤其对于极端荷载或非线性问题，低精度的加速度将直接影响结构分析的计算精度与稳定性。当然，相比cd方法，noh
‑
bathe方法对于求解波动问题具有更高的计算效率，在计算耗时增加仅15％的情况下，可获得比cd法更高精度的计算结果。
4.综上所述，现有的显式时间积分方法的计算精度和计算效率有待提高，尤其对于波动、非线性等问题仍有待进一步提升。


技术实现要素：

5.本发明的目的是提供一种工程结构在极端荷载下动态响应的复合时域分析方法及系统，以提升动力学分析的精度与效率。
6.为实现上述目的，本发明提供了如下方案：
7.本发明提供一种工程结构在极端荷载下动态响应的复合时域分析方法，所述方法包括如下步骤：
8.采用空间单元离散法获取待求解工程结构的动力学控制模型和结构矩阵；所述结构矩阵包括：整体质量矩阵、整体阻尼矩阵和整体刚度矩阵；
9.基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；所述动态响应物理量包括广义位移矢量、广义速度矢量和广义加速度矢量。
10.可选的，所述基于所述动力学控制模型和所述结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，具体包括：
11.利用如下公式确定极端载荷下所述待求解工程结构在第i 1个整体步的前一分步的动态响应物理量：
[0012][0013][0014][0015][0016][0017]
其中，u
i p
、和分别表示在第i 1个整体步的前一分步的广义位移矢量、广义速度矢量和广义加速度矢量，u
i
、和分别表示第i个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，p表示整体步的前一分步占整体步的比例，0＜p＜1，δt表示整体步的步长，m表示整体质量矩阵，k表示整体刚度矩阵，c表示整体阻尼矩阵，表示复合显式时间积分法中的未知系数向量，f
i p
表示第i 1个整体步的前一分步的极端载荷；
[0018]
根据待求解工程结构在第i 1个整体步的前一分步的动态响应物理量，利用如下公式确定极端载荷下所述待求解工程结构在第i 1个整体步的动态响应物理量：
[0019][0020][0021][0022][0023]
其中，表示中间过渡向量，用来修正最终的加速度向量，u
i 1
、、分别表示第i 1个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，α表示的计算式中的计算比例系数，β表示u
i 1
的计算式中的比例系数，r和δ分别表示的计算
式中的计算比例系数和的计算比例系数，f
i 1
表示第i 1个整体步的极端载荷。
[0024]
可选的，所述基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，之前还包括：
[0025]
将初始时刻的广义位移向量和广义速度矢量作为初始广义位移矢量和初始广义速度矢量；
[0026]
根据所述初始广义位移矢量和初始广义速度矢量，利用所述动力学控制模型，求解极端载荷下所述待求解工程结构的初始广义加速度矢量。
[0027]
可选的，所述基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，之前还包括：
[0028]
确定待求解工程结构的自由参数为：
[0029]
p＝0.5，β＝p，
[0030]
可选的，所述基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，之前还包括：
[0031]
确定待求解工程结构的自由参数为：
[0032][0033]
可选的，所述基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，之前还包括：
[0034]
求解方程|k
‑
λm|＝0，获得多个特征根；
[0035]
根据多个所述特征根中的最大特征值确定整体步的步长的取值范围为
[0036]
其中，k表示整体刚度矩阵，m表示整体质量矩阵，λ表示特征根，λ
max
表示最大特征值，δt表示整体步的步长。
[0037]
一种工程结构在极端荷载下动态响应的复合时域分析系统，所述系统包括：
[0038]
动力学控制模型和结构矩阵获取模块，用于采用空间单元离散法获取待求解工程结构的动力学控制模型和结构矩阵；所述结构矩阵包括：整体质量矩阵、整体阻尼矩阵和整体刚度矩阵；
[0039]
动态响应物理量迭代求解模块，用于基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；所述动态响应物理量包括广义位移矢量、广义速度矢量和广义加速度矢量。
[0040]
可选的，所述动态响应物理量迭代求解模块，具体包括：
[0041]
分步长的动态响应物理量迭代求解子模块，用于利用如下公式确定极端载荷下所述待求解工程结构在第i 1个整体步的前一分步的动态响应物理量：
[0042][0043][0044][0045][0046]
其中，u
i p
、和分别表示在第i 1个整体步的前一分步的广义位移矢量、广义速度矢量和广义加速度矢量，u
i
、和分别表示第i个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，p表示整体步的前一分步占整体步的比例，0＜p＜1，δt表示整体步的步长，m表示整体质量矩阵，k表示整体刚度矩阵，c表示整体阻尼矩阵，表示复合显式时间积分法中的未知系数向量，f
i p
表示第i 1个整体步的前一分步的极端载荷；
[0047]
整体步的动态响应物理量迭代求解子模块，用于根据待求解工程结构在第i 1个整体步的前一分步的动态响应物理量，利用如下公式确定极端载荷下所述待求解工程结构在第i 1个整体步的动态响应物理量：
[0048][0049][0050][0051][0052]
其中，表示中间过渡向量，用来修正最终的加速度向量，u
i 1
、、分别表示第i 1个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，α表示的计算式中的计算比例系数，β表示u
i 1
的计算式中的比例系数，r和δ分别表示的计算式中的计算比例系数和的计算比例系数，f
i 1
表示第i 1个整体步的极端载荷。
[0053]
可选的，所述系统还包括：
[0054]
初始广义位移矢量和初始广义速度矢量获取模块，用于将初始时刻的广义位移向量和广义速度矢量作为初始广义位移矢量和初始广义速度矢量；
[0055]
初始广义加速度矢量求解模块，用于根据所述初始广义位移矢量和初始广义速度矢量，利用所述动力学控制模型，求解极端载荷下所述待求解工程结构的初始广义加速度矢量。
[0056]
可选的，所述系统还包括：
[0057]
自由参数确定模块，用于确定待求解工程结构的自由参数为：
[0058]
p＝0.5，β＝p，
[0059]
根据本发明提供的具体实施例，本发明公开了以下技术效果：
[0060]
本发明公开了一种工程结构在极端荷载下动态响应的复合时域分析方法，所述方法包括如下步骤：采用空间单元离散法获取待求解工程结构的动力学控制模型和结构矩阵；基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；本发明采用复合显式时间积分法克服了现有的显示时间积分法计算非线性问题时的技术缺陷，相比于现有的显式时间积分方法的计算精度和计算效率都有较大提升。
[0061]
本发明在整体步的前一分步的动态响应物理量计算过程中引入了未知参数向量，并引入了一次速度修正，理论分析下的计算误差相比显式noh
‑
bathe法有显著降低，进一步提高了计算精度。
附图说明
[0062]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
[0063]
图1为本发明提供的一种工程结构在极端荷载下动态响应的复合时域分析方法的流程图；
[0064]
图2为本发明提供的一种工程结构在极端荷载下动态响应的复合时域分析方法的原理图；
[0065]
图3为本发明实施例1提供的非线性工程结构的结构图原理图；
[0066]
图4为本发明实施例1提供的非线性工程结构的不同算法在0
‑
10s内的响应分析结果对比图；
[0067]
图5为本发明实施例1提供的非线性工程结构的不同算法在795
‑
805s内的响应分析结果对比图；
[0068]
图6为本发明实施例2提供的八层剪切结构的结构原理图；
[0069]
图7为本发明实施例2提供的八层剪切结构在南北地震载荷下不同算法获得的响应分析结果对比图；
[0070]
图8为本发明实施例3提供的空间桁架结构的结构原理图；
[0071]
图9为本发明实施例3提供的空间桁架结构在竖直方向地震载荷下不同算法获得的响应分析结果对比图。
具体实施方式
[0072]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于
本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0073]
本发明的目的是提供一种工程结构在极端荷载下动态响应的复合时域分析方法及系统，以提高显式时间积分方法的计算精度和计算效率。
[0074]
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
[0075]
如图1和2所示，本发明提供了本发明提供一种工程结构在极端荷载下动态响应的复合时域分析方法，所述方法包括如下步骤：
[0076]
步骤101，采用空间单元离散法获取待求解工程结构的动力学控制模型和结构矩阵；所述结构矩阵包括：整体质量矩阵、整体阻尼矩阵和整体刚度矩阵。
[0077]
以有限元法为例，对工程结构进行有限单元离散，即将工程结构划分为有限个在结点处相连接的小单元(假设)，然后利用在各单元内假设的形函数来分片逼近全求解域上的未知场函数，即单元e内任一点的位移a
i
可以用该单元的n个结点位移u
ii
(t)表示：
[0078]
a
i
(x,y,z,t)＝n1(x,y,z)u
ii
(t)
ꢀꢀ
(1)
[0079]
其中，n1(x,y,z)为单元形函数，它是取决于单元类型的已知函数，将上式写成矩阵形式
[0080]
a＝nu
e
ꢀꢀ
(2)
[0081]
其中
[0082]
n＝[n
1 n2ꢀ…ꢀ
n
n
]
ꢀꢀ
(3)
[0083]
式中，a为单元位移场，n为单元形函数矩阵，u
e
为单元结点位移。
[0084]
则单元应变矩阵为
[0085]
ε＝bu
e
ꢀꢀ
(4)
[0086]
其中
[0087]
b＝ln
ꢀꢀ
(5)
[0088][0089]
l为算子矩阵，b为弹性力学几何方程应变与位移关系矩阵，也称应变矩阵。
[0090]
单元应力矩阵为
[0091]
σ＝dbu
e
ꢀꢀ
(7)
[0092]
其中
[0093][0094][0095]
d为弹性力学物理方程中应力与应变的关系矩阵，也称弹性矩阵，v为泊松比，为弹性模量。
[0096]
通过多种方法，如最小势能原理，伽辽金法等，工程结构的内力和外力虚功为
[0097][0098]
其中，k
e
，f
e
分别为单元刚度矩阵和单元等效结点荷载向量，表示为
[0099]
k
e
＝∫b
t
dbdv
e
ꢀꢀ
(10)
[0100][0101]
单元结点位移向量可以表示为
[0102]
u
e
＝g
e
u
ꢀꢀ
(12)
[0103]
其中，u＝[u
11 u
21 u
31
ꢀ…ꢀ
u
1n u
2n u
3n
]
t
为由结构各结点位移组成的结构结点位移向量，n为结构的结点总数，g
e
为由1和0组成的选择矩阵。如对于平面三角形单元(单元的3个结点的结点号分别为i、j、m)
[0104]
[0105]
将式(12)带入式(9)中，可得
[0106][0107]
其中，k为整体刚度矩阵，f为结构的总体荷载向量，分别表示为
[0108][0109][0110]
结构的动能为
[0111][0112]
其中，ρ为工程结构密度，为a，u
e
，u对时间求导结果；m，m
e
分别为工程结构的整体质量矩阵和单元质量矩阵，整体质量矩阵m可由单元质量矩阵m
e
组装而成，可表示为
[0113][0114][0115]
粘滞力的虚功为
[0116][0117]
其中，c，c
e
分别为结构的总体阻尼矩阵和单元阻尼矩阵，表示为
[0118]
[0119][0120]
哈密顿原理公式为
[0121][0122]
其中，t1～t2为所求工程结构响应的时间间隔
[0123]
将以上各式代入到式(23)中，并进行分部积分，得
[0124][0125]
位移向量u在t1和t2时刻的值是给定的，即考虑到δu的任意性，则有
[0126][0127]
至此获得了结构的运动方程
[0128]
对于地震响应分析，外力荷载可表示为
[0129][0130]
其中，为地震荷载加速度向量；r为影响向量，它定义了当基础在荷载激励的方向上产生静态的单位位移时，结构中产生的位移。
[0131]
将式(26)代入到式(25)中，得到地震荷载下结构的运动方程：
[0132][0133]
步骤102，基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；所述动态响应物理量包括广义位移矢量、广义速度矢量和广义加速度矢量。
[0134]
由运动方程式(25)可知，要求解有限元系统的位移、速度和加速度矢量，首要步骤是确定有限元系统的整体质量矩阵、整体阻尼矩阵、整体刚度矩阵以及整体载荷向量。如前所述，逐步时间积分方法还需确定初始广义位移矢量u0、广义速度矢量并由式(25)计算广义加速度矢量对于单自由度系统，上述矩阵和矢量退化为标量。具体步骤为：将初始时刻的广义位移向量和广义速度矢量作为初始广义位移矢量和初始广义速度矢量；根据所述初始广义位移矢量和初始广义速度矢量，利用所述动力学控制模型，求解极端载荷下所述待求解工程结构的初始广义加速度矢量。
[0135]
根据具体工程结构的动力学问题类型和算法稳定区间确定时间步长δt；确定5个
自由参数p,r,α,β和δ。其中，0＜p＜1。
[0136]
(1)对于求解连续型荷载(如正弦荷载)下的结构响应，需限定关系
[0137][0138]
(2)对于求解非连续型荷载(如地震荷载)下的结构响应，需限定关系
[0139][0140]
(3)算法的稳定区间求法:解频率方程|k
‑
λm|＝0，可得n个特征根λ
i
(i＝n)，取其中最大值λ
max
，则系统频率最大值稳定区间δt的范围可写为：
[0141][0142]
l为新方法的稳定区间参数，对于情况(1)，l＝4；对于情况(2)，l＝2.497。
[0143]
基于上述初始化之后，步骤102，所述基于所述动力学控制模型和所述结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，具体包括：
[0144]
利用如下公式确定极端载荷下所述待求解工程结构在第i 1个整体步的前一分步的动态响应物理量，计算第一分步的广义位移矢量u
i p
、广义速度矢量广义加速度矢量通过下述公式依次计算：
[0145][0146][0147][0148][0149]
其中，u
i p
、和分别表示在第i 1个整体步的前一分步的广义位移矢量、广义速度矢量和广义加速度矢量，u
i
、和分别表示第i个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，p表示整体步的前一分步占整体步的比例，0＜p＜1，δt表示整体步的步长，m表示整体质量矩阵，k表示整体刚度矩阵，c表示整体阻尼矩阵，表示复合显式时间积分法中的未知系数向量，f
i p
表示第i 1个整体步的前一分步的极端载荷；
[0150]
根据待求解工程结构在第i 1个整体步的前一分步的动态响应物理量，利用如下公式确定极端载荷下所述待求解工程结构在第i 1个整体步的动态响应物理量，即，计算第二分步，得到整体步的广义位移矢量u
i 1
、广义速度矢量广义加速度矢量通过
下述公式依次计算：
[0151][0152][0153][0154][0155]
其中，表示中间过渡向量，用来修正最终的加速度向量，u
i 1
、、分别表示第i 1个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，α表示的计算式中的计算比例系数，β表示u
i 1
的计算式中的比例系数，r和δ分别表示的计算式中的计算比例系数和的计算比例系数，f
i 1
表示第i 1个整体步的极端载荷。
[0156]
采用上述步骤进行迭代递进计算，将步骤102计算得到的结果进行存储，同时通过图2所示的迭代策略对所有时间步进行递进计算与存储，最后得到所有时间点的计算结果，并将其应用于其他物理量的计算与后处理分析。
[0157]
相比目前已商用化应用的cd法和经典的bathe方法，本发明采用了更多的算法参数控制算法特性，能同时控制计算精度、数值阻尼特性、非线性问题的能量保持、计算效率等特性。通过一系列的算例测试验证了本发明对于结构响应计算问题求解的高效性与准确性。三个实施例都将本发明与cd法和bathe方法进行对比。下面结合算例、附图和实施过程做进一步说明。为实现各种算法的计算时间成本大致相同，采用两个子步的bathe法和本发明的时间步长是cd法的两倍。所有算例的参考解都采用时间步长为δt＝1
×
10
‑6s的隐式bathe算法得到。
[0158]
实施例1
[0159]
本实施例为非线性分析可行性验证的标准算例，主要测试本发明的方法在非线性分析中的可行性和计算精度。如图3所示的两自由度弹性弹簧摆问题。该非线性系统的运动方程写为
[0160][0161][0162]
式中，θ和r分别是角位移和径向位移。m为摆锤的质量，g为重力加速度，l0为弹簧未变形时的长度，k为弹簧的弹性常数。系统的初始条件是
[0163][0164]
为了得到一个高度非线性的情况，令m＝1，g＝9.81，l0＝0.5，k
s
＝98.1，r0＝0.25，＝0.25，和摆锤的径向位移r计算结果如图4
‑
5所示。从图4
‑
5上能明显看出，在0
‑
10s，cd法已经出现较大误差，随后出现不收敛的情况，而bathe法和新方法与精确值吻合较好；但在图5所示的795
‑
805s中，bathe法的计算结果已出现明显大于新方法
的误差。总的来说，新方法在非线性问题中的计算精度高于cd法和bathe法。
[0165]
实施例2
[0166]
本实施例为地震响应分析可行性验证的标准基础算例，主要测试本发明的方法在极端荷载下的可行性和计算精度。如图6所示的八层剪切结构受到1940年el centro南北方向地震载荷作用，作用时长为30秒。该结构为理想简化分析结构，采用层间剪切模型和刚性层假设。采用有限元杆单元离散计算，单元数n＝8，总自由度为8。各楼层的集中质量和楼层抗剪刚度如图6所示。为清楚起见，选取顶层20
‑
22秒时间段内的位移计算结果误差进行分。如图7所示，横坐标为时间，纵坐标为选取不同时间步长的三种方法位移计算结果误差0的对数,t时刻的误差
0,t
由下式计算
[0167][0168]
其中，u
exact,t
为t时刻位移精确值，u
t
为t时刻位移计算值。
[0169]
从图7上能明显看出，本发明的方法在同等时间成本下，相比其它两种算法，大部分计算结果点的误差较低。
[0170]
实施例3
[0171]
本实施例为验证时间积分方法用于一般结构在极端荷载下的计算效率的典型算例。如图8所示的空间桁架结构受到1940年el centro竖直方向地震载荷的作用。该结构高32m，长32m，宽20m，每层高为4m，共8层。前四层长、宽分别为32m，20m，每层平均分成4
×
2个格子；中间两层长、宽分别为24m，20m，每层平均分成3
×
2个格子；最后两层长、宽分别为16m，20m，每层平均分成2
×
2个格子。各构件质量密度为7850kg/m3，桁架各杆件轴向抗拉压刚度为2.1
×
105kn。采用有限元梁单元离散计算，单元数n＝314,总自由度为306。选取节点a处的三种算法的位移计算结果误差进行分析。如图9所示，横坐标为时间，纵坐标为选取不同时间步长的三种方法位移计算结果误差0的对数。为清楚起见，选取了20
‑
21s的误差结果进行分析。从图9上能明显看出，新方法在同等时间成本下，相比其它两种算法，误差显著降低。
[0172]
本发明还提供一种工程结构在极端荷载下动态响应的复合时域分析系统，所述系统包括：
[0173]
动力学控制模型和结构矩阵获取模块，用于采用空间单元离散法获取待求解工程结构的动力学控制模型和结构矩阵；所述结构矩阵包括：整体质量矩阵、整体阻尼矩阵和整体刚度矩阵。
[0174]
动态响应物理量迭代求解模块，用于基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；所述动态响应物理量包括广义位移矢量、广义速度矢量和广义加速度矢量。
[0175]
所述动态响应物理量迭代求解模块，具体包括：
[0176]
分步长的动态响应物理量迭代求解子模块，用于利用如下公式确定极端载荷下所述待求解工程结构在第i 1个整体步的前一分步的动态响应物理量：
[0177]
[0178][0179][0180][0181]
其中，u
i p
、和分别表示在第i 1个整体步的前一分步的广义位移矢量、广义速度矢量和广义加速度矢量，u
i
、和分别表示第i个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，p表示整体步的前一分步占整体步的比例，0＜p＜1，δt表示整体步的步长，m表示整体质量矩阵，k表示整体刚度矩阵，c表示整体阻尼矩阵，表示复合显式时间积分法中的未知系数向量，f
i p
表示第i 1个整体步的前一分步的极端载荷。
[0182]
整体步的动态响应物理量迭代求解子模块，用于根据待求解工程结构在第i 1个整体步的前一分步的动态响应物理量，利用如下公式确定极端载荷下所述待求解工程结构在第i 1个整体步的动态响应物理量：
[0183][0184][0185][0186][0187]
其中，表示中间过渡向量，用来修正最终的加速度向量，u
i 1
、、分别表示第i 1个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，α表示的计算式中的计算比例系数，β表示u
i 1
的计算式中的比例系数，r和δ分别表示的计算式中的计算比例系数和的计算比例系数，f
i 1
表示第i 1个整体步的极端载荷。
[0188]
作为一种优选的实施方式，所述系统还包括：初始广义位移矢量和初始广义速度矢量获取模块，用于将初始时刻的广义位移向量和广义速度矢量作为初始广义位移矢量和初始广义速度矢量；初始广义加速度矢量求解模块，用于根据所述初始广义位移矢量和初始广义速度矢量，利用所述动力学控制模型，求解极端载荷下所述待求解工程结构的初始广义加速度矢量。
[0189]
所述系统还包括：自由参数确定模块，用于确定待求解工程结构的自由参数为：
[0190]
p＝0.5，β＝p，
[0191]
根据本发明提供的具体实施例，本发明公开了以下技术效果：
[0192]
本发明一种工程结构在极端荷载下动态响应的复合时域分析方法，所述方法包括
如下步骤：采用空间单元离散法获取待求解工程结构的动力学控制模型和结构矩阵；基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；本发明采用复合显式时间积分法克服了现有的显示时间积分法计算非线性问题时的技术缺陷，相比于现有的显式时间积分方法的计算精度和计算效率都有较大提升。
[0193]
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。
[0194]
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。