技术领域
[0001] 本
发明涉及一种超高层结构的动力分析数值模拟方法，具体涉及一种无条件稳定且具有高阶
精度的超高层结构的动力时程分析方法。
背景技术
[0002] 超
高层建筑结构规模庞大，有利于我国城市化建设，近年来，国内超高层建筑结构日趋增多，尽管超高层建筑满足了城市化需求，但却存在
地震风险。经调研发现，目前这些超高层建筑相当一部分坐落在高烈度抗震设防地区，具有经历强烈地震的风险，这些超高层建筑一旦在强
震中发生结构倒塌，必将损失惨重。为保证超高层建筑在强震作用下的安全性，有必要对超高层建筑的损伤特点和安全储备进行研究，但由于缺乏超高层震害资料，超高层建筑结构在强震下的损伤特点并不明确，且超高层结构往往规模巨大，无法对其进行实验分析。因此，采用数值模拟手段，对超高层结构进行大量的动力时程分析，成为了分析超高层结构地震损伤规律的最主要手段。
[0003] 动力时程分析通过建立超高层结构的有限元模型，将地震作用的时间历程细化为一个个时间步，然后利用时域积分进行逐步分析，建立代数方程组并求解，最终获得整个时间历程的数值解。目前
现有技术针对超高层建筑结构进行动力时程分析主要采用的是Newmark法， Wilson法和HHT法等，上述分析方法仅仅具有2阶精度。在使用这些具有2阶精度的分析方法进行动力时程分析时，为保证精度，必须采用较短的时间步长。因为超高层结构构件数量庞大，其精细化的有限元模型具有超多
自由度的特点，时间步长选的过短，地震输入时长不变，总的时间步数量过多，也就造成了分析过程时间耗费巨大，对于一般的超高层结构，一次动力时程分析往往持续十几个小时甚至几天，严重阻碍着超高层结构研究与设计。
[0004] 要想解决效率问题，有必要采用具有高阶精度的时域积分分析方法，可通过大步长分析，保证高阶精度的同时，大幅提升超高层结构动力时程分析的精度和效率。目前现有技术提出了多种高阶精度的动力时程分析方法，但由于这些分析方法在每个分析步所求解方程组的规模至少为2阶精度方法求解方程组规模的2倍甚至数倍，因此尽管提升了每步的求解精度，但付出了单步计算量显著增加的代价，通过增大时间步长所提升的效率大打折扣，因此并未得到实质性推广。
发明内容
[0005] 针对现有超高层结构的动力时程分析方法在求解效率方面存在的不足，本发明的目的在于提供一种无条件稳定且具有高阶精度的超高层结构的动力时程分析方法。
[0006] 本发明的技术方案如下：
[0007] 一种超高层结构的动力时程分析方法，其步骤如下：
[0008] 第一步，对超高层结构进行空间有限元离散，建立超高层结构的有限元模型离散系统，梁柱均采用伯努利欧拉梁单元，并采用Rayleigh阻尼建立单元阻尼矩阵，由单元
刚度矩阵、单元
质量矩阵和单元阻尼矩阵集成整体刚度矩阵、整体质量矩阵和整体阻尼矩阵，并由 Hamilton原理导出离散系统的运动方程组：
[0009]
[0010] 第二步，由
地震波每个时刻
加速度值，计算得到离散系统每个质点各个时刻的等效
惯性力，对每个时刻的数值采用线性插值，进而得到离散系统运动方程组的右端项即等效动力荷载Fg：
[0011] 当在基底输入加速度为ug的地震波时，运动方程式组变为：
[0012]
[0013] 将方程左端的Müg移项，并假设F＝0，则上式变为：
[0014]
[0015] 其中，üg为地面加速度，即空间三个方向的平动加速度ü1、ü2和ü3；
[0016] 第三步，选取时间步长，时间步长取为加速度记录时间间隔的n倍，其中n为整数；
[0017] 第四步，逐时间步计算，计算每个时间步结束时刻的位移、速度和加速度，对于第i个时间步，已知ti-1时刻的位移ui-1和速度vi-1，由下式计算第i个时间步上ti时刻的具有高阶精度的位移 速度 和加速度
[0018] G12ui＝P1+Mvi-1-G11ui-1 (1a)
[0019] Mvi＝P2-G21ui-1-G22ui (1b)
[0020] G12ei＝P1-(H11ui-1+H12ui+H13vi-1+H14vi) (1c)
[0021] Mεi＝P2-(H21ui-1+H22ui+H23vi-1+H24vi)-G22ei (1d)
[0022]
[0023]
[0024] 具体计算顺序为
[0025] a)已知ui-1和vi-1，由式(1a)计算ui；
[0026] b)已知ui-1和ui，由式(1b)计算vi；
[0027] c)已知ui-1、vi-1、ui和vi，由式(1c)计算ei；
[0028] d)已知ui-1、vi-1、ui、vi和ei，由式(1d)计算εi；
[0029] e)已知ui和ei，vi和εi，由式(1e)计算具有高阶精度的位移 速度
[0030] f)已知 和 由(1f)计算高阶精度的加速度
[0031] 其中系数矩阵Gij和Hkl以及Pi的表达式为
[0032]
[0033]
[0034]
[0035]
[0036]
[0037]
[0038]
[0039]
[0040]
[0041]
[0042]
[0043]
[0044]
[0045]
[0046] 优选的，所述的第一步当中，单元刚度矩阵为：
[0047]
[0048] 其中 为空间杆件在局部
坐标系中的单元刚度矩阵，是12×12的对称矩阵，对进行坐标变换，即可得到整体坐标系下的单元刚度矩阵刚度矩阵，即
[0049]
[0050] 优选的，所述的第一步当中，单元质量矩阵为：
[0051]
[0052] 采用HRZ法对 进行对
角化，得到集中质量矩阵，即为
[0053]
[0054] 对 进行坐标变换，即可得到整体坐标系下的单元刚度矩阵刚度矩阵，即[0055]
[0056] 优选的，所述的第一步当中，采用Rayleigh阻尼，建立单元阻尼矩阵[0057] Ce＝a0Me+a1Ke
[0058] 其中
[0059]
[0060] ωi和ωj一般分别取结构的第1阶和第3阶
频率，ζ为阻尼比，一般为0.05。
[0061] 优选的，所述的第二步当中，ü1、ü2和ü3为连续函数，采用线性插值来构造ü1、ü2和ü3。
[0062] 优选的，所述的第三步当中，n取为4-10。
[0063] 优选的，还包括第五步：采用逐单元技术优化分析效率，将式(a)右端项计算如下：
[0064]
[0065] 其中
[0066]
[0067]
[0068]
[0069] 上标e表示单元矩阵和向量。式(b)-(f)中右端项采用同样的技术处理。
[0070] 本发明与现有技术相比，优点在于：
[0071] 1)本发明的动力时程分析方法具有3阶精度，比Newmark法等传统二阶动力时程分析方法相比，精度高一阶。
[0072] 2)本发明的动力时程分析方法中，每个时间步仅需对维数Neq(质点总自由度数目)的矩阵求一次逆，计算中能够保持空间有限元离散矩阵的带状稀疏性质，也可并行计算，计算量和Newmark法等传统二阶动力时程分析方法相当，远低于现有文献中其他高阶动力时程分析方法的计算量，因此，本发明的动力时程分析方法可以采用较长的时间步长，也即较少的时间步，得到与传统动力时程分析方法精度相当的结果，可大幅提升求解效率。
[0073] 3)相比Newmark法，本发明的动力时程分析方法具有一定的数值阻尼，可滤掉由于空间离散引起的虚假的高频模式。
[0074] 4)本发明的动力时程分析方法步骤简单，计算效率高，且容易实施。
附图说明
[0076] 图2空间梁单元示意图；
[0077] 图3El-Centro波示意图；
[0078] 图4不同步长下Newmark法和本发明动力时程分析方法计算的A点位移示意图；
[0079] 图5不同步长下Newmark法和本发明动力时程分析方法计算的A点速度示意图；
[0080] 图6不同步长下Newmark法和本发明动力时程分析方法计算的A点加速度示意图。
具体实施方式
[0081] 下面结合具体
实施例来对本发明进行进一步说明，但并不将本发明局限于这些具体实施方式。本领域技术人员应该认识到，本发明涵盖了
权利要求书范围内所可能包括的所有备选方案、改进方案和等效方案。
[0082] 下面结合附图对本发明的结构原理和工作原理作具体的描述：
[0083] 以一个框架结构为实例，具体阐述本发明的无条件稳定且具有高阶精度的超高层结构的动力时程分析方法，框架结构如图1所示，该框架结构的有限元模型包含10层，
水平方向各两跨，每跨6m，层高4米。梁的截面尺寸为0.2×0.4m，柱截面尺寸为0.4×0.4m，
密度为 2.5×103kg/m3，；图1中框架底部箭头方向为地震波输入方向。所述动力时程分析方法包含如下步骤：
[0084] 第一步，建立超高层结构的有限元模型，梁柱均采用伯努利欧拉梁单元，导出运动方程组；空间梁单元是超高层结构中最典型的构件，如图2所示。单元位移向量和单元杆端位移向量分别为
[0085]
[0086] 单元刚度矩阵为
[0087]
[0088] 其中 为空间杆件在局部坐标系中的单元刚度矩阵，是12×12的对称矩阵。
[0089] 对 进行坐标变换，即可得到整体坐标系下的单元刚度矩阵刚度矩阵，即[0090]
[0091] 单元一致质量矩阵为
[0092]
[0093] 采用HRZ法对 进行对角化，得到集中质量矩阵，即为
[0094]
[0095] 对 进行坐标变换，即可得到整体坐标系下的单元刚度矩阵刚度矩阵，即[0096]
[0097] 采用Rayleigh阻尼，建立单元阻尼矩阵
[0098] Ce＝a0Me+a1Ke (7)
[0099] 其中
[0100]
[0101] 对于一般结构。这里ωi和ωj一般分别取结构的第1阶和第3阶频率。ζ为阻尼比，一般为0.05。
[0102] 以上得到了单元刚度矩阵、质量矩阵和阻尼矩阵之后，集成整体刚度矩阵、质量矩阵和阻尼矩阵，由Hamilton原理，可得到运动方程组如下
[0103]
[0104] 第二步，由地震波每个时刻加速度值，计算得到离散系统每个质点各个时刻的等效惯性力，对每个时刻的数值采用线性插值，进而得到离散系统运动方程组的右端项。
[0105] 当在基底输入加速度为ug的地震波时，运动方程式(9)变为
[0106]
[0107] 将方程左端的 移项，并假设F＝0，则式10变为
[0108]
[0109] 这里-Müg为等效动力荷载Fg。üg为Neq阶向量，其中结构的每个旋转自由度对应的元素为0，而每个平动自由度对应的元素分别为地面x、y、z三个方向的加速度ü1、ü2和ü3。实际工程计算中，地面加速度为空间三个方向的平动加速度ü1、ü2和ü3，且每个方向的地面运动加速度常常为一系列离散时刻的值ü1j、ü2j、ü3j(j＝1…Ng)，本发明的动力时程分析方法要求荷载为随时间变化的连续函数，即要求ü1j、ü2j、ü3j为连续函数，这里只需要对ü1j、ü2j、ü3j分别采用线性插值来构造ü1、ü2和ü3即可。
[0110] 第三步，选取时间步长，由于本发明的动力时程分析方法具有比Newmark法、HHT法等常用2阶动力时程分析方法更高阶的精度，因此可采用较大时间步长以提升求解效率，一般选为加速度记录时间间隔的n倍(n为整数，可取为4-10)。
[0111] 第四步，计算每个时间步结束时刻的位移、速度和加速度。具体的，对于第i个时间步，已知ti-1时刻的位移ui-1和速度vi-1，由式(12a)计算第i个时间步上ti时刻的具有高阶精度的位移 速度 和加速度
[0112] G12ui＝P1+Mvi-1-G11ui-1 (12a)
[0113] Mvi＝P2-G21ui-1-G22ui (12b)
[0114] G12ei＝P1-(H11ui-1+H12ui+H13vi-1+H14vi) (12c)
[0115] Mεi＝P2-(H21ui-1+H22ui+H23vi-1+H24vi)-G22ei (12d)
[0116]
[0117]
[0118] 具体计算顺序为
[0119] g)已知ui-1和vi-1，由式(12a)计算ui；
[0120] h)已知ui-1和ui，由式(12b)计算vi；
[0121] i)已知ui-1、vi-1、ui和vi，由式(12c)计算ei；
[0122] j)已知ui-1、vi-1、ui、vi和ei，由式(12d)计算εi；
[0123] k)已知ui和ei，vi和εi，由式(12e)计算具有高阶精度的位移 速度
[0124] l)已知 和 由(12f)计算高阶精度的加速度
[0125] 其中系数矩阵Gij和Gkl以及Pi的表达式为
[0126]
[0127]
[0128]
[0129]
[0130]
[0131]
[0132]
[0133]
[0134]
[0135]
[0136]
[0137] 对于线弹性问题，若每个时间步长不变，则每个时间步的G12不变，整个计算过程只T T需要对G12做一次LDL分解。对于非线性问题，每个分析步仅需要对G12做一次LDL分解，计算量和Newmark法相当。
[0138] 第五步：采用逐单元以及并行计算优化效率。以下采用一种逐单元(Element by Element，简称EBE)技术优化分析效率，采用EBE技术的好处在于，可以将大规模矩阵分析简化为逐单元分析，即先进行单元级别的矩阵向量乘法运算，后集成整体向量。这样不仅可以减少计算量和存储量，而且单元级别的矩阵向量乘法运算可以并行，可进一步提升计算效率。以下以式(12a)右端项的计算为例进行说明。
[0139]
[0140] 其中
[0141]
[0142]
[0143]
[0144] 上标e表示单元矩阵和向量。式(12b)-(12f)中右端项采用同样的技术处理。
[0145] 在基底输入如图3所示的El Centro地震波，该地震波的加速度记录间隔为0.01s，本发明的动力时程分析方法中，选取分析步长为加速度记录的10倍，即Δt＝0.1s，计算框架顶部A点位移响应。采用集中质量矩阵，阻尼比选取0.05，计算得到的Rayleigh阻尼系数选为a0＝0.3247，a1＝0.005241，各个质点的初始位移和初始速度均为0。
[0146] 为了展示本发明的动力时程分析方法高效性和精确性，首先采用最常用的二阶精度 Newmark法进行分析，Newark法的参数选为α＝1/2，β＝1/4，步长采用Δt＝0.01s，分析结果如图4-6所示，耗费时间为21.4s。仍然采用Newmark法，尝试将步长放大原来10倍，即令Δt＝0.1s，其结果已经背离Δt＝0.01s时的结果，该结果说明，采用Newmark法分析，若放大时间步长，会导致分析结果严重背离。而采用本发明的动力时程分析方法，令步长分别为Δt＝0.01s和Δt＝0.1s，其结果和采用Newmark法，步长为Δt＝0.01s的结果非常接近(3条响应曲线几乎重合)，如图4-6所示。该比较例非常直观的表明，当Δt＝0.1s时，即本发明的动力时程分析方法当步长为Newmark法的10倍时，结果仍然非常精确，耗费时间仅仅为4.4s，约为Newmark法的1/5-1/4。充分验证了本发明超高层结构的动力时程分析方法的精确性和高效性。
[0147] 应当理解的是，本发明描述的方法的步骤仅仅是示例性的描述，对其先后进行的时间顺序没有特殊的要求，除非其本身有必然的先后顺序关系。
[0148] 如上所示，本发明虽然已参照有限的实施例和附图进行了说明，但在本发明所属领域中具备通常知识的人均可以从此记载中进行各种
修改和
变形。由此，其他实施例及权利要求书与等同物均属于权利要求的保护范围。
