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

摘要

本发明涉及了一种工程结构在极端荷载下动态响应的复合时域分析方法及系统，所述方法包括如下步骤：采用空间单元离散法获取待求解工程结构的动力学控制模型和结构矩阵；基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；本发明采用复合显式时间积分法，克服了现有的显示时间积分法计算非线性问题时的技术缺陷，相比于现有的显式时间积分方法的计算精度和计算效率都有较大提升。

说明书

技术领域

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

背景技术

随着计算机技术的迅猛发展，时间积分方法已成为求解结构动态响应的重要方法，并大量应用于商业有限元软件中。时间积分法也称逐步积分法，是假定已知时刻t＝0的位移U

在现有的显式方法中，中心差分法(CD法)是一种广泛应用于结构动力学问题的经典的单步显式计算方法，并已应用到ANSYS、ABAQUS和ADINA等著名商用软件中。虽然，目前的CD法通过引入人工阻尼可以有效改进算法的数值耗散特性，但效果仍有待提升。数值耗散的目的是滤除或有效减少离散后的动力系统的虚假高频模态。实际上，显式方法要在高频段获得理想的数值耗散而不在低频段引入较大的耗散是不容易的。Noh和Bathe提出了一种具有两个子步的复合显式方法(称为Noh-Bathe法，Noh G,Bathe K-J.An explicit timeintegration scheme for the analysis of wave propagations.Computers&Structures2013)，这种显式方法采用了子步的策略，可以引入更多的自由参数以控制算法的基本性能，并且能有效控制算法的数值耗散特性，该方法已集成在商用软件ADINA中。尽管如此，Noh-Bathe法对于无阻尼情况可以获得二阶算法精度，对于有阻尼情况的加速度仅获得一阶精度，因而其算法精度仍有待提升，尤其对于极端荷载或非线性问题，低精度的加速度将直接影响结构分析的计算精度与稳定性。当然，相比CD方法，Noh-Bathe方法对于求解波动问题具有更高的计算效率，在计算耗时增加仅15％的情况下，可获得比CD法更高精度的计算结果。

综上所述，现有的显式时间积分方法的计算精度和计算效率有待提高，尤其对于波动、非线性等问题仍有待进一步提升。

发明内容

本发明的目的是提供一种工程结构在极端荷载下动态响应的复合时域分析方法及系统，以提升动力学分析的精度与效率。

为实现上述目的，本发明提供了如下方案：

本发明提供一种工程结构在极端荷载下动态响应的复合时域分析方法，所述方法包括如下步骤：

采用空间单元离散法获取待求解工程结构的动力学控制模型和结构矩阵；所述结构矩阵包括：整体质量矩阵、整体阻尼矩阵和整体刚度矩阵；

基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；所述动态响应物理量包括广义位移矢量、广义速度矢量和广义加速度矢量。

可选的，所述基于所述动力学控制模型和所述结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，具体包括：

利用如下公式确定极端载荷下所述待求解工程结构在第i+1个整体步的前一分步的动态响应物理量：

其中，U

根据待求解工程结构在第i+1个整体步的前一分步的动态响应物理量，利用如下公式确定极端载荷下所述待求解工程结构在第i+1个整体步的动态响应物理量：

其中，

可选的，所述基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，之前还包括：

将初始时刻的广义位移向量和广义速度矢量作为初始广义位移矢量和初始广义速度矢量；

根据所述初始广义位移矢量和初始广义速度矢量，利用所述动力学控制模型，求解极端载荷下所述待求解工程结构的初始广义加速度矢量。

可选的，所述基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，之前还包括：

确定待求解工程结构的自由参数为：

p＝0.5，

可选的，所述基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，之前还包括：

确定待求解工程结构的自由参数为：

可选的，所述基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，之前还包括：

求解方程|K-λM|＝0，获得多个特征根；

根据多个所述特征根中的最大特征值确定整体步的步长的取值范围为

其中，K表示整体刚度矩阵，M表示整体质量矩阵，λ表示特征根，λ

一种工程结构在极端荷载下动态响应的复合时域分析系统，所述系统包括：

动力学控制模型和结构矩阵获取模块，用于采用空间单元离散法获取待求解工程结构的动力学控制模型和结构矩阵；所述结构矩阵包括：整体质量矩阵、整体阻尼矩阵和整体刚度矩阵；

动态响应物理量迭代求解模块，用于基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；所述动态响应物理量包括广义位移矢量、广义速度矢量和广义加速度矢量。

可选的，所述动态响应物理量迭代求解模块，具体包括：

分步长的动态响应物理量迭代求解子模块，用于利用如下公式确定极端载荷下所述待求解工程结构在第i+1个整体步的前一分步的动态响应物理量：

其中，U

整体步的动态响应物理量迭代求解子模块，用于根据待求解工程结构在第i+1个整体步的前一分步的动态响应物理量，利用如下公式确定极端载荷下所述待求解工程结构在第i+1个整体步的动态响应物理量：

其中，

可选的，所述系统还包括：

初始广义位移矢量和初始广义速度矢量获取模块，用于将初始时刻的广义位移向量和广义速度矢量作为初始广义位移矢量和初始广义速度矢量；

初始广义加速度矢量求解模块，用于根据所述初始广义位移矢量和初始广义速度矢量，利用所述动力学控制模型，求解极端载荷下所述待求解工程结构的初始广义加速度矢量。

可选的，所述系统还包括：

自由参数确定模块，用于确定待求解工程结构的自由参数为：

p＝0.5，

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明公开了一种工程结构在极端荷载下动态响应的复合时域分析方法，所述方法包括如下步骤：采用空间单元离散法获取待求解工程结构的动力学控制模型和结构矩阵；基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；本发明采用复合显式时间积分法克服了现有的显示时间积分法计算非线性问题时的技术缺陷，相比于现有的显式时间积分方法的计算精度和计算效率都有较大提升。

本发明在整体步的前一分步的动态响应物理量计算过程中引入了未知参数向量，并引入了一次速度修正，理论分析下的计算误差相比显式Noh-Bathe法有显著降低，进一步提高了计算精度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明提供的一种工程结构在极端荷载下动态响应的复合时域分析方法的流程图；

图2为本发明提供的一种工程结构在极端荷载下动态响应的复合时域分析方法的原理图；

图3为本发明实施例1提供的非线性工程结构的结构图原理图；

图4为本发明实施例1提供的非线性工程结构的不同算法在0-10s内的响应分析结果对比图；

图5为本发明实施例1提供的非线性工程结构的不同算法在795-805s内的响应分析结果对比图；

图6为本发明实施例2提供的八层剪切结构的结构原理图；

图7为本发明实施例2提供的八层剪切结构在南北地震载荷下不同算法获得的响应分析结果对比图；

图8为本发明实施例3提供的空间桁架结构的结构原理图；

图9为本发明实施例3提供的空间桁架结构在竖直方向地震载荷下不同算法获得的响应分析结果对比图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的目的是提供一种工程结构在极端荷载下动态响应的复合时域分析方法及系统，以提高显式时间积分方法的计算精度和计算效率。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

如图1和2所示，本发明提供了本发明提供一种工程结构在极端荷载下动态响应的复合时域分析方法，所述方法包括如下步骤：

步骤101，采用空间单元离散法获取待求解工程结构的动力学控制模型和结构矩阵；所述结构矩阵包括：整体质量矩阵、整体阻尼矩阵和整体刚度矩阵。

以有限元法为例，对工程结构进行有限单元离散，即将工程结构划分为有限个在结点处相连接的小单元(假设)，然后利用在各单元内假设的形函数来分片逼近全求解域上的未知场函数，即单元e内任一点的位移a

a

其中，N

A＝NU

其中

N＝[N

式中，A为单元位移场，N为单元形函数矩阵，U

则单元应变矩阵为

ε＝BU

其中

B＝LN (5)

L为算子矩阵，B为弹性力学几何方程应变与位移关系矩阵，也称应变矩阵。

单元应力矩阵为

σ＝DBU

其中

D为弹性力学物理方程中应力与应变的关系矩阵，也称弹性矩阵，v为泊松比，为弹性模量。

通过多种方法，如最小势能原理，伽辽金法等，工程结构的内力和外力虚功为

其中，K

K

单元结点位移向量可以表示为

U

其中，U＝[U

将式(12)带入式(9)中，可得

其中，K为整体刚度矩阵，F为结构的总体荷载向量，分别表示为

结构的动能为

其中，ρ为工程结构密度，

粘滞力的虚功为

其中，C，C

哈密顿原理公式为

其中，t

将以上各式代入到式(23)中，并进行分部积分，得

位移向量U在t

至此获得了结构的运动方程

对于地震响应分析，外力荷载可表示为

其中，

将式(26)代入到式(25)中，得到地震荷载下结构的运动方程：

步骤102，基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；所述动态响应物理量包括广义位移矢量、广义速度矢量和广义加速度矢量。

由运动方程式(25)可知，要求解有限元系统的位移、速度和加速度矢量，首要步骤是确定有限元系统的整体质量矩阵、整体阻尼矩阵、整体刚度矩阵以及整体载荷向量。如前所述，逐步时间积分方法还需确定初始广义位移矢量U

根据具体工程结构的动力学问题类型和算法稳定区间确定时间步长Δt；确定5个自由参数p,r,α,β和δ。其中，0＜p＜1。

(1)对于求解连续型荷载(如正弦荷载)下的结构响应，需限定关系

(2)对于求解非连续型荷载(如地震荷载)下的结构响应，需限定关系

(3)算法的稳定区间求法:解频率方程|K-λM|＝0，可得n个特征根λ

l为新方法的稳定区间参数，对于情况(1)，l＝4；对于情况(2)，l＝2.497。

基于上述初始化之后，步骤102，所述基于所述动力学控制模型和所述结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，具体包括：

利用如下公式确定极端载荷下所述待求解工程结构在第i+1个整体步的前一分步的动态响应物理量，计算第一分步的广义位移矢量U

其中，U

根据待求解工程结构在第i+1个整体步的前一分步的动态响应物理量，利用如下公式确定极端载荷下所述待求解工程结构在第i+1个整体步的动态响应物理量，即，计算第二分步，得到整体步的广义位移矢量U

其中，

采用上述步骤进行迭代递进计算，将步骤102计算得到的结果进行存储，同时通过图2所示的迭代策略对所有时间步进行递进计算与存储，最后得到所有时间点的计算结果，并将其应用于其他物理量的计算与后处理分析。

相比目前已商用化应用的CD法和经典的Bathe方法，本发明采用了更多的算法参数控制算法特性，能同时控制计算精度、数值阻尼特性、非线性问题的能量保持、计算效率等特性。通过一系列的算例测试验证了本发明对于结构响应计算问题求解的高效性与准确性。三个实施例都将本发明与CD法和Bathe方法进行对比。下面结合算例、附图和实施过程做进一步说明。为实现各种算法的计算时间成本大致相同，采用两个子步的Bathe法和本发明的时间步长是CD法的两倍。所有算例的参考解都采用时间步长为Δt＝1×10

实施例1

本实施例为非线性分析可行性验证的标准算例，主要测试本发明的方法在非线性分析中的可行性和计算精度。如图3所示的两自由度弹性弹簧摆问题。该非线性系统的运动方程写为

式中，θ和r分别是角位移和径向位移。m为摆锤的质量，g为重力加速度，L

为了得到一个高度非线性的情况，令m＝1，g＝9.81，L

实施例2

本实施例为地震响应分析可行性验证的标准基础算例，主要测试本发明的方法在极端荷载下的可行性和计算精度。如图6所示的八层剪切结构受到1940年El Centro南北方向地震载荷作用，作用时长为30秒。该结构为理想简化分析结构，采用层间剪切模型和刚性层假设。采用有限元杆单元离散计算，单元数N＝8，总自由度为8。各楼层的集中质量和楼层抗剪刚度如图6所示。为清楚起见，选取顶层20-22秒时间段内的位移计算结果误差进行分。如图7所示，横坐标为时间，纵坐标为选取不同时间步长的三种方法位移计算结果误差

其中，U

从图7上能明显看出，本发明的方法在同等时间成本下，相比其它两种算法，大部分计算结果点的误差较低。

实施例3

本实施例为验证时间积分方法用于一般结构在极端荷载下的计算效率的典型算例。如图8所示的空间桁架结构受到1940年El Centro竖直方向地震载荷的作用。该结构高32m，长32m，宽20m，每层高为4m，共8层。前四层长、宽分别为32m，20m，每层平均分成4×2个格子；中间两层长、宽分别为24m，20m，每层平均分成3×2个格子；最后两层长、宽分别为16m，20m，每层平均分成2×2个格子。各构件质量密度为7850Kg/m

本发明还提供一种工程结构在极端荷载下动态响应的复合时域分析系统，所述系统包括：

动力学控制模型和结构矩阵获取模块，用于采用空间单元离散法获取待求解工程结构的动力学控制模型和结构矩阵；所述结构矩阵包括：整体质量矩阵、整体阻尼矩阵和整体刚度矩阵。

动态响应物理量迭代求解模块，用于基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；所述动态响应物理量包括广义位移矢量、广义速度矢量和广义加速度矢量。

所述动态响应物理量迭代求解模块，具体包括：

分步长的动态响应物理量迭代求解子模块，用于利用如下公式确定极端载荷下所述待求解工程结构在第i+1个整体步的前一分步的动态响应物理量：

其中，U

整体步的动态响应物理量迭代求解子模块，用于根据待求解工程结构在第i+1个整体步的前一分步的动态响应物理量，利用如下公式确定极端载荷下所述待求解工程结构在第i+1个整体步的动态响应物理量：

其中，

作为一种优选的实施方式，所述系统还包括：初始广义位移矢量和初始广义速度矢量获取模块，用于将初始时刻的广义位移向量和广义速度矢量作为初始广义位移矢量和初始广义速度矢量；初始广义加速度矢量求解模块，用于根据所述初始广义位移矢量和初始广义速度矢量，利用所述动力学控制模型，求解极端载荷下所述待求解工程结构的初始广义加速度矢量。

所述系统还包括：自由参数确定模块，用于确定待求解工程结构的自由参数为：

p＝0.5，

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明一种工程结构在极端荷载下动态响应的复合时域分析方法，所述方法包括如下步骤：采用空间单元离散法获取待求解工程结构的动力学控制模型和结构矩阵；基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；本发明采用复合显式时间积分法克服了现有的显示时间积分法计算非线性问题时的技术缺陷，相比于现有的显式时间积分方法的计算精度和计算效率都有较大提升。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。