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

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

Info

Publication number
CN113792461B
CN113792461B CN202111072634.4A CN202111072634A CN113792461B CN 113792461 B CN113792461 B CN 113792461B CN 202111072634 A CN202111072634 A CN 202111072634A CN 113792461 B CN113792461 B CN 113792461B
Authority
CN
China
Prior art keywords
engineering structure
generalized
integral
dynamic response
matrix
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202111072634.4A
Other languages
English (en)
Other versions
CN113792461A (zh
Inventor
温伟斌
何旭辉
段晟昱
刘添豪
黄方林
方岱宁
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN202111072634.4A priority Critical patent/CN113792461B/zh
Publication of CN113792461A publication Critical patent/CN113792461A/zh
Application granted granted Critical
Publication of CN113792461B publication Critical patent/CN113792461B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Geometry (AREA)
  • Computational Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computer Hardware Design (AREA)
  • Pure & Applied Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

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

Description

一种工程结构在极端荷载下动态响应的复合时域分析方法
技术领域
本发明涉及结构动力学技术领域，特别是涉及一种工程结构在极端荷载下动态响应的复合时域分析方法及系统。
背景技术
随着计算机技术的迅猛发展，时间积分方法已成为求解结构动态响应的重要方法，并大量应用于商业有限元软件中。时间积分法也称逐步积分法，是假定已知时刻t＝0的位移U0、速度
Figure BDA0003260989530000011
并将时间求解时间域平均分成n个时间段0、Δt、2Δt、…、nΔt(nΔt为总时长)，且已求得0、Δt、2Δt、…和t的解，建立求解t+Δt时刻解的递推算法。总体来说，时间积分法可以分成两类：显式时间积分法和隐式时间积分法。显式方法通常是条件稳定的，所需的时间步长远小于隐式方法。然而，显式方法比隐式方法花费更少的时间来计算矩阵。显式方法由于其无需进行方程组的迭代求解，在求解大型结构或非线性程度较高的问题时更具优势，因此受到很多外学者的广泛关注和研究。
在现有的显式方法中，中心差分法(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个整体步的前一分步的动态响应物理量：
Figure BDA0003260989530000021
Figure BDA0003260989530000022
Figure BDA0003260989530000023
Figure BDA0003260989530000031
Figure BDA0003260989530000032
其中，Ui+p、
Figure BDA0003260989530000033
和
Figure BDA0003260989530000034
分别表示在第i+1个整体步的前一分步的广义位移矢量、广义速度矢量和广义加速度矢量，Ui、
Figure BDA0003260989530000035
和
Figure BDA0003260989530000036
分别表示第i个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，p表示整体步的前一分步占整体步的比例，0＜p＜1，Δt表示整体步的步长，M表示整体质量矩阵，K表示整体刚度矩阵，C表示整体阻尼矩阵，
Figure BDA0003260989530000037
表示复合显式时间积分法中的未知系数向量，Fi+p表示第i+1个整体步的前一分步的极端载荷；
根据待求解工程结构在第i+1个整体步的前一分步的动态响应物理量，利用如下公式确定极端载荷下所述待求解工程结构在第i+1个整体步的动态响应物理量：
Figure BDA0003260989530000038
Figure BDA0003260989530000039
Figure BDA00032609895300000310
Figure BDA00032609895300000311
其中，
Figure BDA00032609895300000312
表示中间过渡向量，用来修正最终的加速度向量，Ui+1、
Figure BDA00032609895300000313
Figure BDA00032609895300000314
分别表示第i+1个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，α表示
Figure BDA00032609895300000315
的计算式中
Figure BDA00032609895300000316
的计算比例系数，β表示Ui+1的计算式中
Figure BDA00032609895300000317
的比例系数，r和δ分别表示
Figure BDA00032609895300000318
的计算式中
Figure BDA00032609895300000319
的计算比例系数和
Figure BDA00032609895300000320
的计算比例系数，Fi+1表示第i+1个整体步的极端载荷。
可选的，所述基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，之前还包括：
将初始时刻的广义位移向量和广义速度矢量作为初始广义位移矢量和初始广义速度矢量；
根据所述初始广义位移矢量和初始广义速度矢量，利用所述动力学控制模型，求解极端载荷下所述待求解工程结构的初始广义加速度矢量。
可选的，所述基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，之前还包括：
确定待求解工程结构的自由参数为：
p＝0.5，
Figure BDA0003260989530000041
β＝p，
Figure BDA0003260989530000042
可选的，所述基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，之前还包括：
确定待求解工程结构的自由参数为：
Figure BDA0003260989530000043
可选的，所述基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，之前还包括：
求解方程|K-λM|＝0，获得多个特征根；
根据多个所述特征根中的最大特征值确定整体步的步长的取值范围为
Figure BDA0003260989530000044
其中，K表示整体刚度矩阵，M表示整体质量矩阵，λ表示特征根，λmax表示最大特征值，Δt表示整体步的步长。
一种工程结构在极端荷载下动态响应的复合时域分析系统，所述系统包括：
动力学控制模型和结构矩阵获取模块，用于采用空间单元离散法获取待求解工程结构的动力学控制模型和结构矩阵；所述结构矩阵包括：整体质量矩阵、整体阻尼矩阵和整体刚度矩阵；
动态响应物理量迭代求解模块，用于基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；所述动态响应物理量包括广义位移矢量、广义速度矢量和广义加速度矢量。
可选的，所述动态响应物理量迭代求解模块，具体包括：
分步长的动态响应物理量迭代求解子模块，用于利用如下公式确定极端载荷下所述待求解工程结构在第i+1个整体步的前一分步的动态响应物理量：
Figure BDA0003260989530000051
Figure BDA0003260989530000052
Figure BDA0003260989530000053
Figure BDA0003260989530000054
其中，Ui+p、
Figure BDA0003260989530000055
和
Figure BDA0003260989530000056
分别表示在第i+1个整体步的前一分步的广义位移矢量、广义速度矢量和广义加速度矢量，Ui、
Figure BDA0003260989530000057
和
Figure BDA0003260989530000058
分别表示第i个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，p表示整体步的前一分步占整体步的比例，0＜p＜1，Δt表示整体步的步长，M表示整体质量矩阵，K表示整体刚度矩阵，C表示整体阻尼矩阵，
Figure BDA0003260989530000059
表示复合显式时间积分法中的未知系数向量，Fi+p表示第i+1个整体步的前一分步的极端载荷；
整体步的动态响应物理量迭代求解子模块，用于根据待求解工程结构在第i+1个整体步的前一分步的动态响应物理量，利用如下公式确定极端载荷下所述待求解工程结构在第i+1个整体步的动态响应物理量：
Figure BDA00032609895300000510
Figure BDA00032609895300000511
Figure BDA00032609895300000512
Figure BDA00032609895300000513
其中，
Figure BDA00032609895300000514
表示中间过渡向量，用来修正最终的加速度向量，Ui+1、
Figure BDA00032609895300000515
Figure BDA00032609895300000516
分别表示第i+1个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，α表示
Figure BDA00032609895300000517
的计算式中
Figure BDA00032609895300000518
的计算比例系数，β表示Ui+1的计算式中
Figure BDA00032609895300000519
的比例系数，r和δ分别表示
Figure BDA0003260989530000061
的计算式中
Figure BDA0003260989530000062
的计算比例系数和
Figure BDA0003260989530000063
的计算比例系数，Fi+1表示第i+1个整体步的极端载荷。
可选的，所述系统还包括：
初始广义位移矢量和初始广义速度矢量获取模块，用于将初始时刻的广义位移向量和广义速度矢量作为初始广义位移矢量和初始广义速度矢量；
初始广义加速度矢量求解模块，用于根据所述初始广义位移矢量和初始广义速度矢量，利用所述动力学控制模型，求解极端载荷下所述待求解工程结构的初始广义加速度矢量。
可选的，所述系统还包括：
自由参数确定模块，用于确定待求解工程结构的自由参数为：
p＝0.5，
Figure BDA0003260989530000064
β＝p，
Figure BDA0003260989530000065
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明公开了一种工程结构在极端荷载下动态响应的复合时域分析方法，所述方法包括如下步骤：采用空间单元离散法获取待求解工程结构的动力学控制模型和结构矩阵；基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；本发明采用复合显式时间积分法克服了现有的显示时间积分法计算非线性问题时的技术缺陷，相比于现有的显式时间积分方法的计算精度和计算效率都有较大提升。
本发明在整体步的前一分步的动态响应物理量计算过程中引入了未知参数向量，并引入了一次速度修正，理论分析下的计算误差相比显式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内任一点的位移ai可以用该单元的n个结点位移uiI(t)表示：
ai(x,y,z,t)＝N1(x,y,z)uiI(t) (1)
其中，N1(x,y,z)为单元形函数，它是取决于单元类型的已知函数，将上式写成矩阵形式
A＝NUe (2)
其中
N＝[N1 N2 … Nn] (3)
式中，A为单元位移场，N为单元形函数矩阵，Ue为单元结点位移。
则单元应变矩阵为
ε＝BUe (4)
其中
B＝LN (5)
Figure BDA0003260989530000081
L为算子矩阵，B为弹性力学几何方程应变与位移关系矩阵，也称应变矩阵。
单元应力矩阵为
σ＝DBUe (7)
其中
Figure BDA0003260989530000091
Figure BDA0003260989530000092
D为弹性力学物理方程中应力与应变的关系矩阵，也称弹性矩阵，v为泊松比，为弹性模量。
通过多种方法，如最小势能原理，伽辽金法等，工程结构的内力和外力虚功为
Figure BDA0003260989530000093
其中，Ke，Fe分别为单元刚度矩阵和单元等效结点荷载向量，表示为
Ke＝∫BTDBdVe (10)
Figure BDA0003260989530000094
单元结点位移向量可以表示为
Ue＝GeU (12)
其中，U＝[U11 U21 U31 … U1N U2N U3N]T为由结构各结点位移组成的结构结点位移向量，N为结构的结点总数，Ge为由1和0组成的选择矩阵。如对于平面三角形单元(单元的3个结点的结点号分别为I、J、M)
Figure BDA0003260989530000101
将式(12)带入式(9)中，可得
Figure BDA0003260989530000102
其中，K为整体刚度矩阵，F为结构的总体荷载向量，分别表示为
Figure BDA0003260989530000103
Figure BDA0003260989530000104
结构的动能为
Figure BDA0003260989530000105
其中，ρ为工程结构密度，
Figure BDA0003260989530000106
为A，Ue，U对时间求导结果；M，Me分别为工程结构的整体质量矩阵和单元质量矩阵，整体质量矩阵M可由单元质量矩阵Me组装而成，可表示为
Figure BDA0003260989530000107
Figure BDA0003260989530000108
粘滞力的虚功为
Figure BDA0003260989530000111
其中，C，Ce分别为结构的总体阻尼矩阵和单元阻尼矩阵，表示为
Figure BDA0003260989530000112
Figure BDA0003260989530000113
哈密顿原理公式为
Figure BDA0003260989530000114
其中，t1～t2为所求工程结构响应的时间间隔
将以上各式代入到式(23)中，并进行分部积分，得
Figure BDA0003260989530000115
位移向量U在t1和t2时刻的值是给定的，即
Figure BDA0003260989530000116
考虑到δU的任意性，则有
Figure BDA0003260989530000117
至此获得了结构的运动方程
对于地震响应分析，外力荷载可表示为
Figure BDA0003260989530000118
其中，
Figure BDA0003260989530000119
为地震荷载加速度向量；R为影响向量，它定义了当基础在荷载激励的方向上产生静态的单位位移时，结构中产生的位移。
将式(26)代入到式(25)中，得到地震荷载下结构的运动方程：
Figure BDA0003260989530000121
步骤102，基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；所述动态响应物理量包括广义位移矢量、广义速度矢量和广义加速度矢量。
由运动方程式(25)可知，要求解有限元系统的位移、速度和加速度矢量，首要步骤是确定有限元系统的整体质量矩阵、整体阻尼矩阵、整体刚度矩阵以及整体载荷向量。如前所述，逐步时间积分方法还需确定初始广义位移矢量U0、广义速度矢量
Figure BDA0003260989530000122
并由式(25)计算广义加速度矢量
Figure BDA0003260989530000123
对于单自由度系统，上述矩阵和矢量退化为标量。具体步骤为：将初始时刻的广义位移向量和广义速度矢量作为初始广义位移矢量和初始广义速度矢量；根据所述初始广义位移矢量和初始广义速度矢量，利用所述动力学控制模型，求解极端载荷下所述待求解工程结构的初始广义加速度矢量。
根据具体工程结构的动力学问题类型和算法稳定区间确定时间步长Δt；确定5个自由参数p,r,α,β和δ。其中，0＜p＜1。
(1)对于求解连续型荷载(如正弦荷载)下的结构响应，需限定关系
Figure BDA0003260989530000124
(2)对于求解非连续型荷载(如地震荷载)下的结构响应，需限定关系
Figure BDA0003260989530000125
(3)算法的稳定区间求法:解频率方程|K-λM|＝0，可得n个特征根λi(i＝n)，取其中最大值λmax，则系统频率最大值
Figure BDA0003260989530000126
稳定区间Δt的范围可写为：
Figure BDA0003260989530000127
l为新方法的稳定区间参数，对于情况(1)，l＝4；对于情况(2)，l＝2.497。
基于上述初始化之后，步骤102，所述基于所述动力学控制模型和所述结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，具体包括：
利用如下公式确定极端载荷下所述待求解工程结构在第i+1个整体步的前一分步的动态响应物理量，计算第一分步的广义位移矢量Ui+p、广义速度矢量
Figure BDA0003260989530000131
广义加速度矢量
Figure BDA0003260989530000132
通过下述公式依次计算：
Figure BDA0003260989530000133
Figure BDA0003260989530000134
Figure BDA0003260989530000135
Figure BDA0003260989530000136
其中，Ui+p、
Figure BDA0003260989530000137
和
Figure BDA0003260989530000138
分别表示在第i+1个整体步的前一分步的广义位移矢量、广义速度矢量和广义加速度矢量，Ui、
Figure BDA0003260989530000139
和
Figure BDA00032609895300001310
分别表示第i个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，p表示整体步的前一分步占整体步的比例，0＜p＜1，Δt表示整体步的步长，M表示整体质量矩阵，K表示整体刚度矩阵，C表示整体阻尼矩阵，
Figure BDA00032609895300001311
表示复合显式时间积分法中的未知系数向量，Fi+p表示第i+1个整体步的前一分步的极端载荷；
根据待求解工程结构在第i+1个整体步的前一分步的动态响应物理量，利用如下公式确定极端载荷下所述待求解工程结构在第i+1个整体步的动态响应物理量，即，计算第二分步，得到整体步的广义位移矢量Ui+1、广义速度矢量
Figure BDA00032609895300001312
广义加速度矢量
Figure BDA00032609895300001313
通过下述公式依次计算：
Figure BDA00032609895300001314
Figure BDA00032609895300001315
Figure BDA00032609895300001316
Figure BDA00032609895300001317
其中，
Figure BDA00032609895300001318
表示中间过渡向量，用来修正最终的加速度向量，Ui+1、
Figure BDA00032609895300001319
Figure BDA0003260989530000141
分别表示第i+1个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，α表示
Figure BDA0003260989530000142
的计算式中
Figure BDA0003260989530000143
的计算比例系数，β表示Ui+1的计算式中
Figure BDA0003260989530000144
的比例系数，r和δ分别表示
Figure BDA0003260989530000145
的计算式中
Figure BDA0003260989530000146
的计算比例系数和
Figure BDA0003260989530000147
的计算比例系数，Fi+1表示第i+1个整体步的极端载荷。
采用上述步骤进行迭代递进计算，将步骤102计算得到的结果进行存储，同时通过图2所示的迭代策略对所有时间步进行递进计算与存储，最后得到所有时间点的计算结果，并将其应用于其他物理量的计算与后处理分析。
相比目前已商用化应用的CD法和经典的Bathe方法，本发明采用了更多的算法参数控制算法特性，能同时控制计算精度、数值阻尼特性、非线性问题的能量保持、计算效率等特性。通过一系列的算例测试验证了本发明对于结构响应计算问题求解的高效性与准确性。三个实施例都将本发明与CD法和Bathe方法进行对比。下面结合算例、附图和实施过程做进一步说明。为实现各种算法的计算时间成本大致相同，采用两个子步的Bathe法和本发明的时间步长是CD法的两倍。所有算例的参考解都采用时间步长为Δt＝1×10-6s的隐式Bathe算法得到。
实施例1
本实施例为非线性分析可行性验证的标准算例，主要测试本发明的方法在非线性分析中的可行性和计算精度。如图3所示的两自由度弹性弹簧摆问题。该非线性系统的运动方程写为
Figure BDA0003260989530000148
Figure BDA0003260989530000149
式中，θ和r分别是角位移和径向位移。m为摆锤的质量，g为重力加速度，L0为弹簧未变形时的长度，k为弹簧的弹性常数。系统的初始条件是
Figure BDA00032609895300001410
为了得到一个高度非线性的情况，令m＝1，g＝9.81，L0＝0.5，ks＝98.1，r0＝0.25，
Figure BDA00032609895300001411
Figure BDA00032609895300001412
和
Figure BDA00032609895300001413
摆锤的径向位移r计算结果如图4-5所示。从图4-5上能明显看出，在0-10s，CD法已经出现较大误差，随后出现不收敛的情况，而Bathe法和新方法与精确值吻合较好；但在图5所示的795-805s中，Bathe法的计算结果已出现明显大于新方法的误差。总的来说，新方法在非线性问题中的计算精度高于CD法和Bathe法。
实施例2
本实施例为地震响应分析可行性验证的标准基础算例，主要测试本发明的方法在极端荷载下的可行性和计算精度。如图6所示的八层剪切结构受到1940年El Centro南北方向地震载荷作用，作用时长为30秒。该结构为理想简化分析结构，采用层间剪切模型和刚性层假设。采用有限元杆单元离散计算，单元数N＝8，总自由度为8。各楼层的集中质量和楼层抗剪刚度如图6所示。为清楚起见，选取顶层20-22秒时间段内的位移计算结果误差进行分。如图7所示，横坐标为时间，纵坐标为选取不同时间步长的三种方法位移计算结果误差0的对数,t时刻的误差0,t由下式计算
Figure BDA0003260989530000151
其中，Uexact,t为t时刻位移精确值，Ut为t时刻位移计算值。
从图7上能明显看出，本发明的方法在同等时间成本下，相比其它两种算法，大部分计算结果点的误差较低。
实施例3
本实施例为验证时间积分方法用于一般结构在极端荷载下的计算效率的典型算例。如图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上能明显看出，新方法在同等时间成本下，相比其它两种算法，误差显著降低。
本发明还提供一种工程结构在极端荷载下动态响应的复合时域分析系统，所述系统包括：
动力学控制模型和结构矩阵获取模块，用于采用空间单元离散法获取待求解工程结构的动力学控制模型和结构矩阵；所述结构矩阵包括：整体质量矩阵、整体阻尼矩阵和整体刚度矩阵。
动态响应物理量迭代求解模块，用于基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；所述动态响应物理量包括广义位移矢量、广义速度矢量和广义加速度矢量。
所述动态响应物理量迭代求解模块，具体包括：
分步长的动态响应物理量迭代求解子模块，用于利用如下公式确定极端载荷下所述待求解工程结构在第i+1个整体步的前一分步的动态响应物理量：
Figure BDA0003260989530000161
Figure BDA0003260989530000162
Figure BDA0003260989530000163
Figure BDA0003260989530000164
其中，Ui+p、
Figure BDA0003260989530000165
和
Figure BDA0003260989530000166
分别表示在第i+1个整体步的前一分步的广义位移矢量、广义速度矢量和广义加速度矢量，Ui、
Figure BDA0003260989530000167
和
Figure BDA0003260989530000168
分别表示第i个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，p表示整体步的前一分步占整体步的比例，0＜p＜1，Δt表示整体步的步长，M表示整体质量矩阵，K表示整体刚度矩阵，C表示整体阻尼矩阵，
Figure BDA0003260989530000169
表示复合显式时间积分法中的未知系数向量，Fi+p表示第i+1个整体步的前一分步的极端载荷。
整体步的动态响应物理量迭代求解子模块，用于根据待求解工程结构在第i+1个整体步的前一分步的动态响应物理量，利用如下公式确定极端载荷下所述待求解工程结构在第i+1个整体步的动态响应物理量：
Figure BDA0003260989530000171
Figure BDA0003260989530000172
Figure BDA0003260989530000173
Figure BDA0003260989530000174
其中，
Figure BDA0003260989530000175
表示中间过渡向量，用来修正最终的加速度向量，Ui+1、
Figure BDA0003260989530000176
Figure BDA0003260989530000177
分别表示第i+1个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，α表示
Figure BDA0003260989530000178
的计算式中
Figure BDA0003260989530000179
的计算比例系数，β表示Ui+1的计算式中
Figure BDA00032609895300001710
的比例系数，r和δ分别表示
Figure BDA00032609895300001711
的计算式中
Figure BDA00032609895300001712
的计算比例系数和
Figure BDA00032609895300001713
的计算比例系数，Fi+1表示第i+1个整体步的极端载荷。
作为一种优选的实施方式，所述系统还包括：初始广义位移矢量和初始广义速度矢量获取模块，用于将初始时刻的广义位移向量和广义速度矢量作为初始广义位移矢量和初始广义速度矢量；初始广义加速度矢量求解模块，用于根据所述初始广义位移矢量和初始广义速度矢量，利用所述动力学控制模型，求解极端载荷下所述待求解工程结构的初始广义加速度矢量。
所述系统还包括：自由参数确定模块，用于确定待求解工程结构的自由参数为：
p＝0.5，
Figure BDA00032609895300001714
β＝p，
Figure BDA00032609895300001715
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明一种工程结构在极端荷载下动态响应的复合时域分析方法，所述方法包括如下步骤：采用空间单元离散法获取待求解工程结构的动力学控制模型和结构矩阵；基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；本发明采用复合显式时间积分法克服了现有的显示时间积分法计算非线性问题时的技术缺陷，相比于现有的显式时间积分方法的计算精度和计算效率都有较大提升。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (8)

1.一种工程结构在极端荷载下动态响应的复合时域分析方法，其特征在于，所述方法包括如下步骤：
采用空间单元离散法获取待求解工程结构的动力学控制模型和结构矩阵；所述结构矩阵包括：整体质量矩阵、整体阻尼矩阵和整体刚度矩阵；
基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；所述动态响应物理量包括广义位移矢量、广义速度矢量和广义加速度矢量；
所述基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，具体包括：
利用如下公式确定极端载荷下所述待求解工程结构在第i+1个整体步的前一分步的动态响应物理量：
Figure FDA0003771436530000011
Figure FDA0003771436530000012
Figure FDA0003771436530000013
Figure FDA0003771436530000014
其中，Ui+p、
Figure FDA0003771436530000015
和
Figure FDA0003771436530000016
分别表示在第i+1个整体步的前一分步的广义位移矢量、广义速度矢量和广义加速度矢量，Ui、
Figure FDA0003771436530000017
和
Figure FDA0003771436530000018
分别表示第i个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，p表示整体步的前一分步占整体步的比例，0＜p＜1，Δt表示整体步的步长，M表示整体质量矩阵，K表示整体刚度矩阵，C表示整体阻尼矩阵，
Figure FDA0003771436530000019
表示复合显式时间积分法中的未知系数向量，Fi+p表示第i+1个整体步的前一分步的极端载荷；
根据待求解工程结构在第i+1个整体步的前一分步的动态响应物理量，利用如下公式确定极端载荷下所述待求解工程结构在第i+1个整体步的动态响应物理量：
Figure FDA0003771436530000021
Figure FDA0003771436530000022
Figure FDA0003771436530000023
Figure FDA0003771436530000024
其中，
Figure FDA0003771436530000025
表示中间过渡向量，Ui+1、
Figure FDA0003771436530000026
分别表示第i+1个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，α表示
Figure FDA0003771436530000027
的计算式中
Figure FDA0003771436530000028
的计算比例系数，β表示Ui+1的计算式中
Figure FDA0003771436530000029
的比例系数，r和δ分别表示
Figure FDA00037714365300000210
的计算式中
Figure FDA00037714365300000211
的计算比例系数和
Figure FDA00037714365300000212
的计算比例系数，Fi+1表示第i+1个整体步的极端载荷。
2.根据权利要求1所述的工程结构在极端荷载下动态响应的复合时域分析方法，其特征在于，所述基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，之前还包括：
将初始时刻的广义位移向量和广义速度矢量作为初始广义位移矢量和初始广义速度矢量；
根据所述初始广义位移矢量和初始广义速度矢量，利用所述动力学控制模型，求解极端载荷下所述待求解工程结构的初始广义加速度矢量。
3.根据权利要求1所述的工程结构在极端荷载下动态响应的复合时域分析方法，其特征在于，所述基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，之前还包括：
确定待求解工程结构的自由参数为：
p＝0.5，
Figure FDA00037714365300000213
β＝p，
Figure FDA00037714365300000214
4.根据权利要求1所述的工程结构在极端荷载下动态响应的复合时域分析方法，其特征在于，所述基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，之前还包括：
确定待求解工程结构的自由参数为：
Figure FDA0003771436530000031
5.根据权利要求1所述的工程结构在极端荷载下动态响应的复合时域分析方法，其特征在于，所述基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量，之前还包括：
求解方程
Figure FDA0003771436530000032
获得多个特征根；
根据多个所述特征根中的最大特征值确定整体步的步长的取值范围为
Figure FDA0003771436530000033
其中，K表示整体刚度矩阵，M表示整体质量矩阵，λ表示特征根，λmax表示最大特征值，Δt表示整体步的步长，l表示稳定区间参数。
6.一种工程结构在极端荷载下动态响应的复合时域分析系统，其特征在于，所述系统包括：
动力学控制模型和结构矩阵获取模块，用于采用空间单元离散法获取待求解工程结构的动力学控制模型和结构矩阵；所述结构矩阵包括：整体质量矩阵、整体阻尼矩阵和整体刚度矩阵；
动态响应物理量迭代求解模块，用于基于所述动力学控制模型和结构矩阵，采用复合显式时间积分法对所述待求解工程结构进行分步长迭代求解，确定极端载荷下所述待求解工程结构在每个整体步的动态响应物理量；所述动态响应物理量包括广义位移矢量、广义速度矢量和广义加速度矢量；
所述动态响应物理量迭代求解模块，具体包括：
分步长的动态响应物理量迭代求解子模块，用于利用如下公式确定极端载荷下所述待求解工程结构在第i+1个整体步的前一分步的动态响应物理量：
Figure FDA0003771436530000034
Figure FDA0003771436530000041
Figure FDA0003771436530000042
Figure FDA0003771436530000043
其中，Ui+p、
Figure FDA0003771436530000044
和
Figure FDA0003771436530000045
分别表示在第i+1个整体步的前一分步的广义位移矢量、广义速度矢量和广义加速度矢量，Ui、
Figure FDA0003771436530000046
和
Figure FDA0003771436530000047
分别表示第i个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，p表示整体步的前一分步占整体步的比例，0＜p＜1，Δt表示整体步的步长，M表示整体质量矩阵，K表示整体刚度矩阵，C表示整体阻尼矩阵，
Figure FDA0003771436530000048
表示复合显式时间积分法中的未知系数向量，Fi+p表示第i+1个整体步的前一分步的极端载荷；
整体步的动态响应物理量迭代求解子模块，用于根据待求解工程结构在第i+1个整体步的前一分步的动态响应物理量，利用如下公式确定极端载荷下所述待求解工程结构在第i+1个整体步的动态响应物理量：
Figure FDA0003771436530000049
Figure FDA00037714365300000410
Figure FDA00037714365300000411
Figure FDA00037714365300000412
其中，
Figure FDA00037714365300000413
表示中间过渡向量，Ui+1、
Figure FDA00037714365300000414
分别表示第i+1个整体步的广义位移矢量、广义速度矢量和广义加速度矢量，α表示
Figure FDA00037714365300000415
的计算式中
Figure FDA00037714365300000416
的计算比例系数，β表示Ui+1的计算式中
Figure FDA00037714365300000417
的比例系数，r和δ分别表示
Figure FDA00037714365300000418
的计算式中
Figure FDA00037714365300000419
的计算比例系数和
Figure FDA00037714365300000420
的计算比例系数，Fi+1表示第i+1个整体步的极端载荷。
7.根据权利要求6所述的工程结构在极端荷载下动态响应的复合时域分析系统，其特征在于，所述系统还包括：
初始广义位移矢量和初始广义速度矢量获取模块，用于将初始时刻的广义位移向量和广义速度矢量作为初始广义位移矢量和初始广义速度矢量；
初始广义加速度矢量求解模块，用于根据所述初始广义位移矢量和初始广义速度矢量，利用所述动力学控制模型，求解极端载荷下所述待求解工程结构的初始广义加速度矢量。
8.根据权利要求6所述的工程结构在极端荷载下动态响应的复合时域分析系统，其特征在于，所述系统还包括：
自由参数确定模块，用于确定待求解工程结构的自由参数为：
p＝0.5，
Figure FDA0003771436530000051
β＝p，
Figure FDA0003771436530000052
CN202111072634.4A 2021-09-14 2021-09-14 一种工程结构在极端荷载下动态响应的复合时域分析方法 Active CN113792461B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111072634.4A CN113792461B (zh) 2021-09-14 2021-09-14 一种工程结构在极端荷载下动态响应的复合时域分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111072634.4A CN113792461B (zh) 2021-09-14 2021-09-14 一种工程结构在极端荷载下动态响应的复合时域分析方法

Publications (2)

Publication Number Publication Date
CN113792461A CN113792461A (zh) 2021-12-14
CN113792461B true CN113792461B (zh) 2022-11-01

Family

ID=79183220

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111072634.4A Active CN113792461B (zh) 2021-09-14 2021-09-14 一种工程结构在极端荷载下动态响应的复合时域分析方法

Country Status (1)

Country Link
CN (1) CN113792461B (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104615888A (zh) * 2015-02-06 2015-05-13 华北水利水电大学 一种基于广义最小残差方法的桥梁移动车辆荷载识别方法
CN106096257A (zh) * 2016-06-06 2016-11-09 武汉理工大学 一种非线性索单元分析方法及系统

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR100989190B1 (ko) * 2008-08-29 2010-10-20 한양대학교 산학협력단 등가정하중을 이용한 위상최적설계방법
CN101510230A (zh) * 2009-03-11 2009-08-19 同济大学 车辆道路载荷的仿真方法
CN105136592A (zh) * 2015-05-14 2015-12-09 华北水利水电大学 一种判断桥墩抗震性能的方法
CN112199799B (zh) * 2020-10-29 2022-06-28 中南大学 一种工程结构振动响应的时域分析方法和系统

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104615888A (zh) * 2015-02-06 2015-05-13 华北水利水电大学 一种基于广义最小残差方法的桥梁移动车辆荷载识别方法
CN106096257A (zh) * 2016-06-06 2016-11-09 武汉理工大学 一种非线性索单元分析方法及系统

Also Published As

Publication number Publication date
CN113792461A (zh) 2021-12-14

Similar Documents

Publication Publication Date Title
CN109902404B (zh) 不同阻尼形式的结构时程响应积分的统一递推计算方法
CN109460622B (zh) 一种大规模建筑结构的完全显式的动力时程分析方法
CN113111547A (zh) 基于缩减基的频域有限元模型修正方法
Matthies et al. Nonlinear Galerkin methods for the model reduction of nonlinear dynamical systems
CN106354954A (zh) 一种基于叠层基函数的三维力学模态仿真模拟方法
CN109902357B (zh) 一种复变差分的光滑非线性结构动响应灵敏度分析方法
CN111640296A (zh) 交通流预测方法、系统、存储介质及终端
CN111090942A (zh) 基于拓扑优化的高灵敏度压阻式单轴力传感器设计方法
CN113792461B (zh) 一种工程结构在极端荷载下动态响应的复合时域分析方法
Chwalowski et al. Flutter Prediction Report in Support of the High Angle Working Group at the Third Aeroelastic Prediction Workshop
CN112199799B (zh) 一种工程结构振动响应的时域分析方法和系统
CN110147571B (zh) 一种组件结构的拓扑优化方法及装置
CN111241728A (zh) 一种欧拉方程的间断伽辽金有限元数值求解方法
Yang et al. Direct‐Iterative Hybrid Solution in Nonlinear Dynamic Structural Analysis
CN115455779A (zh) 一种工程结构振动响应的复合显式时域分析方法
CN115906333A (zh) 一种桁架结构的几何非线性等效板动力学建模与响应分析方法
CN110348158B (zh) 一种基于分区异步长解的地震波动分析方法
CN112487689B (zh) 基于统计ckf模型更新混合试验方法
Chang et al. A one-parameter controlled dissipative unconditionally stable explicit algorithm for time history analysis
CN114004117A (zh) 考虑土体参数空间变异性的边坡地震滑移概率分析方法
Rezayibana The effect of soil type on seismic response of tall telecommunication towers with random vibration analysis
CN113722859B (zh) 一种基于凸多面体模型的不确定性结构静力响应确定方法
JP2903098B2 (ja) 構造設計方法
CN115995277B (zh) 一种材料动力学特性评估方法、装置、设备及介质
CN110569611B (zh) 一种基于多复变量法的结构频响函数灵敏度分析方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant