一种大规模建筑结构的完全显式的动力时程分析方法与流程

本发明涉及一种大规模建筑结构的动力分析数值模拟方法，具体涉及一种大规模建筑结构的无需阻尼矩阵对角化、具有一定方法阻尼、具有大稳定步长、具有3阶精度、自起步的完全显式的动力时程分析方法。

背景技术：

近年来，国内大规模建筑结构日趋增多，尽管大规模建筑结构满足了城市化需求，但却存在地震风险。为规避大规模建筑结构的地震风险，在设计阶段，常常采用动力分析数值模拟手段，对大规模建筑结构进行大量的动力时程分析，以保证大规模建筑结构在突发强震作用下的安全。为便于计算，人们提出了一些例如底部剪力法、反应谱法等方法。然而随着土木工程结构形式的不断发展和对时程分析可靠度越来越高的要求，我们需要对结构在动力荷载下的行为进行更为精确的预测。加上有限元方法的发展与计算机性能的迅猛提高，人们对时程分析法的关注越来越多。动力时程分析法要求我们计算出结构位移、速度等量随时间的定量变化，是描述结构在动力作用下力学行为的最精准的分析方法。对于普遍的非线性结构，直接积分方法是进行动力时程分析的最强有力工具。

直接积分方法可以分为隐式方法与显式方法。常用的隐式方法有平均加速度法、wilson-θ法、hht-α法等等，这些隐式方法均为无条件稳定方法，计算中可以采用的时间步长不受模型网格尺寸限制，效率较高，且有一定的数值阻尼，能够过滤掉结构空间有限元划分造成的虚假高阶模态，因此常常被应用于结构动力时程分析这种动力响应以低阶模态为主的分析。但由于隐式方法在每一增量步内都需要对静态平衡方程进行迭代求解，并且每次迭代都需要求解大型的线性方程组，这一过程需要占用相当数量的计算资源、磁盘空间和内存，反而表现出较低的计算效率，且对于非线性比较强的复杂结构模型，非线性迭代不容易收敛。因此，近些年以中心差分法为代表的显式方法逐渐被应用于结构动力时程分析，如国际知名软件abaqus、ls-dyna等等。相比隐式方法，显式方法的每步计算中无需求解大规模方程组，仅需要进行单元级别的矩阵和向量乘法运算以及整体向量加法运算，在计算效率上有明显的优势，且显式方法不存在无法收敛的问题，因此在计算强震作用下建筑结构的动力响应、倒塌模拟这类复杂模型的强非线性响应方面具有先天的优势，近年来国内也逐渐的涌现出一批采用显式方法进行结构动力分析的软件，如pkpm-sausage等。

尽管中心差分方法具有效率上的优势，但仍然存在以下几个非常棘手的问题：(1)为避免求解大规模的整体平衡方程组，中心差分法要求阻尼矩阵为对角矩阵，而结构动力时程分析中，阻尼矩阵的对角化往往通过采用质量阻尼或者振型阻尼来实现，这种处理往往会引入较大的误差，限制了该方法的应用范围。若采用对速度的偏心差分代替中心差分，可避免阻尼矩阵对角化，但偏心差分方法只有一阶精度，计算结果可靠性收到影响；(2)当结构中存在非线性粘滞阻尼器，使得阻尼矩阵为非线性矩阵时，若采用中心差分法，每个时间步需要像隐式方法一样进行非线性迭代并求解大规模方程组，显式方法的效率优势和收敛优势则无从体现；(3)结构的空间有限元离散，不可避免的引入了虚假的高阶模态，由于中心差分法无数值耗散特性，无法抑制这些虚假高频响应，会造成很大的误差；(4)显式方法普遍存在稳定性问题，分析步长不允许超过临界步长，临界步长正比于结构最短周期，中心差分法的临界步长为结构最短周期的1/π倍，普遍被认为是显式方法中最长的，但临界步长仍有拓展的空间。以上几大问题制约了中心差分方法在结构动力时程分析中的应用。除此之外，由于中心差分法为多步法，还存在方法启动问题，造成了编程上的麻烦。另外，中心差分方法具有二阶精度，尽管已经满足大部分工程计算的需求，但是对于超长持时分析，误差的累积会造成计算结果偏离理论解。

技术实现要素：

针对中心差分法在大规模复杂建筑结构非线性动力时程分析方面的缺陷，本发明的目的在于提供一种大规模建筑结构的无需阻尼矩阵对角化、具有一定方法阻尼、具有大稳定步长、具有3阶精度、自起步的完全显式的动力时程分析方法。

本发明的技术方案如下：

一种大规模建筑结构的显式动力分析方法，步骤如下：

第一步，对大规模建筑结构进行空间有限元离散，建立大规模建筑结构的有限元模型离散系统，梁柱均采用伯努利欧拉梁单元，并采用rayleigh阻尼建立单元阻尼矩阵，由单元刚度矩阵、单元质量矩阵和单元阻尼矩阵集成整体刚度矩阵k、整体质量矩阵m和整体阻尼矩阵c，并由hamilton原理导出离散系统的运动方程组，其中f为右端项：

第二步，根据体系最短周期，选取分析步长；计算建筑结构有限元模型的最短周期，记为tmin；选取时间步长δt＜0.45tmin，令δt＝0.4tmin；

第三步，逐时间步计算，对于第i个时间步，已知ti-1时刻的位移ui-1和速度vi-1，由下式计算第i个时间步上ti时刻的位移ui、速度vi；

其中，pi的表达式为

第四步，对ti时刻的位移ui和速度vi进行修正，获得具有高阶精度的位移速度并计算加速度

ei＝m^-1(p1-h11ui-1-h12ui-h13vi-1-h14vi)δt

其中系数矩阵hij的表达式为

优选地，所述的第一步当中，单元刚度矩阵为：

其中为空间杆件在局部坐标系中的单元刚度矩阵，是12×12的对称矩阵，对进行坐标变换，即可得到整体坐标系下的单元刚度矩阵刚度矩阵，即

优选地，所述的第一步当中，单元质量矩阵为：

采用hrz法对进行对角化，得到集中质量矩阵，即为

对进行坐标变换，即可得到整体坐标系下的单元刚度矩阵刚度矩阵，即

优选地，所述的第一步当中，采用rayleigh阻尼，建立单元阻尼矩阵

c^e＝a0m^e+a1k^e

其中

ωi和ωj一般分别取结构的第1阶和第3阶频率，ζ为阻尼比，一般为0.05。

优选地，所述的第一步当中，对于单元右端项f^e的计算，需对每个时刻的地震波加速度记录采用线性插值得到由下式计算

其中m^e为单元质量矩阵。

优选地，首先计算并存贮各个单元的单元刚度矩阵k^e、质量矩阵m^e和阻尼矩阵c^e以及右端项f^e；其中m^e为对角矩阵，仅需要对质量矩阵m^e进行集成得到整体质量矩阵m，而无需计算整体刚度矩阵c和整体阻尼矩阵m以及右端项f。

优选地，所述的第二步当中，建筑结构有限元模型的最短周期tmin可采用估计的方法，无需实际计算。

优选地，所述的第三步当中，采用单元层级并行方法，ui的计算方法为：记首先计算各个单元的此处可采用各个单元并行方法，然后由集成得到r1，的计算采用同样方式处理。

优选地，所述的第四步当中，同样采用单元层级并行方法，计算ei、εi和

本发明与现有技术相比，优点在于：

(1)本发明的动力时程分析方法只需要对质量矩阵对角化，而无需对阻尼矩阵对角化，因此可以采用常用的rayleigh阻尼等非对角阻尼矩阵；

(2)本发明的动力时程分析方法同样适用于存在非线性粘滞阻尼器的结构体系的显式动力时程分析，每个时间步上无需迭代求解；

(3)本发明的动力时程分析方法对于高阶模态具有较强的数值耗散特性，可以屏蔽掉由于空间离散造成的虚假的高频响应；

(4)本发明的动力时程分析方法具有较长的稳定区间，对于无阻尼系统，方法临界步长为结构最短周期的0.4502倍，比中心差分法的临界步长增大近50％，谱半径曲线，如图1所示；

(5)本发明的动力时程分析方法具有高阶精度，对于无阻尼体系，方法具有4阶精度，对于阻尼体系，方法具有3阶精度，均高于中心差分法；其中加速度收敛阶与中心差分法的对比如图2所示；

(6)本发明的动力时程分析方法为单步法，不存在起步问题；

(7)本发明的动力时程分析方法实施简便，效率高，只需在单元层级进行矩阵向量乘法，而无需集成整体刚度、质量、阻尼矩阵进行计算，可并行。

附图说明

图1本发明的动力时程分析方法与中心差分法的谱半径曲线对比，其中虚线为中心差分法，点划线为本发明方法。

图2本发明的动力时程分析方法与中心差分法的加速度收敛性对比，其中虚线为中心差分法，点划线为本发明方法。

图3某2层框架结构。

图4el-centro波。

图5本发明的动力时程分析方法与中心差分法的加速度精度对比，其中点划线为中心差分法，实线为本发明方法。

具体实施方式

下面结合具体实施例来对本发明进行进一步说明，但并不将本发明局限于这些具体实施方式。本领域技术人员应该认识到，本发明涵盖了权利要求书范围内所可能包括的所有备选方案、改进方案和等效方案。

下面结合附图对本发明的结构原理和工作原理作具体的描述：

以一个2层框架结构为实例，具体阐述本发明所涉及的完全显式动力时程分析方法，框架结构如图3所示，下部箭头方向为地震波，该框架结构的有限元模型包含2层，水平两个方向跨度各为6m，层高4米。梁的截面尺寸为0.2×0.4m，柱截面尺寸为0.4×0.4m，密度为2.5×10³kg/m³。所述的完全显式的动力时程分析方法包含如下步骤：

第一步，对上述框架结构进行空间有限元离散，建立上述框架结构的有限元模型离散系统，梁柱均采用伯努利欧拉梁单元，并采用rayleigh阻尼建立单元阻尼矩阵，由单元刚度矩阵、单元质量矩阵和单元阻尼矩阵集成整体刚度矩阵k、整体质量矩阵m和整体阻尼矩阵c，并由hamilton原理导出离散系统的运动方程组，其中f为右端项：

第二步，选取分析步长。为与中心差分法对比，选取满足中心差分法稳定性的步长，令δt＝0.002s；

第三步，逐时间步计算，对于第i个时间步，已知ti-1时刻的位移ui-1和速度vi-1，由下式计算第i个时间步上ti时刻的位移ui、速度vi；

其中，pi的表达式为

第四步，对ti时刻的位移ui和速度vi进行修正，获得具有高阶精度的位移速度并计算加速度

其中系数矩阵hij的表达式为

在基底输入如图4所示的elcentro地震波，该地震波的加速度记录间隔为0.01s，计算框架一层顶部a点位移响应。采用集中质量矩阵，阻尼比选取0，各个质点的初始位移和初始速度均为0。

为了展示本发明的完全显式动力时程分析方法高效性和精确性，首先采用本发明方法进行一次动力分析，分析步长选取第二步中所提到的δt＝0.002s。同时采用中心差分法对该框架进行动力分析，同样选取步长同样选取δt＝0.002s，为方便对比，选取更小的步长δt＝0.0002s，采用中心差分法进行计算，作为精确解的参考解。由本发明方法和中心差分法计算得到的a点加速度响应如图5所示，图中只给出39秒-40秒的响应。经对比发现，采用本发明的动力时程分析方法的计算结果与参考解几乎重合，而中心差分法的计算结果严重偏离参考解，该例充分验证了本发明显式动力时程分析方法高效性和精确性。

应当理解的是，本发明描述的方法的步骤仅仅是示例性的描述，对其先后进行的时间顺序没有特殊的要求，除非其本身有必然的先后顺序关系。

如上所示，本发明虽然已参照有限的实施例和附图进行了说明，但在本发明所属领域中具备通常知识的人均可以从此记载中进行各种修改和变形。由此，其他实施例及权利要求书与等同物均属于权利要求的保护范围。