基于拟力法的增量动力分析方法

2019-02-22郝润霞杨作续余丁浩

振动与冲击 2019年4期

郝润霞，杨作续，李钢，余丁浩，贾硕

(1. 内蒙古科技大学土木工程学院，内蒙古包头 014010；2. 大连理工大学海岸和近海工程国家重点实验室，辽宁大连 116024)

增量动力分析(Incremental Dynamic Analysis,IDA)方法[1-2]是基于性态地震工程(Performance-based Earthquake Engineering,PBEE)中确定结构在不同水准地震动响应的参数化分析方法，可以准确的观察结构从弹性、弹塑性直至动力失稳倒塌的全过程中结构性态的发展变化规律[3]，同时可以克服静力推覆分析(Push-Over Analysis,POA)方法特别是在结构接近倒塌时因近似静力分析而存在的问题[4]。然而其分析计算时运算量大、耗时长等计算效率问题，却严重地影响了其在工程应用中的发展。

近年来，国内外学者围绕IDA的计算效率问题开展了一系列的研究。在简化分析方法方面，如基于静力推覆分析的简化近似IDA(SPO2IDA)方法[5]，该方法主要适用于基本振型控制的结构体系；基于模态推覆分析(Modal Pushover Analysis ,MPA)的简化快速模态增量分析(Modal Incremental Dynamic Analysis, MIDA)方法[6-7]，该方法考虑了模态选取、结构延性和滞回模型对其计算精度的影响，当考虑的高阶振型贡献越多时[8]，其计算精度越高，且保证精度的前提下，能够有效提高计算效率，但对于不规则结构体系适用性较差[9]；基于等效二自由度模型(Equivalent Dual Degree of Freedom ,EDDOF)的简化IDA方法[10]，该方法以单向偏心结构的等效二自由度模型应用于IDA，对单向偏心结构进行抗震性能评估，在保证较高的精度前提下，计算效率得到了大幅提高。对比分析SPO2IDA方法与MIDA方法，SPO2IDA方法计算效率更为突出，但其计算精度较差[11]。在优化非线性求解方面，以快速非线性时程分析(Fast Nonlinear Analysis,FNA)方法[12]作为计算工具提出的快速增量动力分析(Fast Incremental Dynamic Analysis ,FIDA)方法[13]，与传统IDA方法、MIDA方法相比，FIDA方法在计算效率与精度上优势明显。此外，如基于Nataf变换点估计法与单地震动记录IDA相结合的单地震动记录随机IDA方法[14]，该方法可以使IDA分析的计算精度和计算效率达到平衡。以IDA曲线混合统计方法进行统计分析，并直接以结构离散IDA曲线进行地震易损性曲线及地震失效概率计算，可有效避免多余的统计与拟合[15]。IDA方法在结构非线性动力时程分析时常采用传统变刚度方法，该方法在任一时刻非线性迭代步求解时，结构整体刚度矩阵需要实时合成与分解，计算量大，效率低。

拟力法(Force Analogy Method,FAM)[16-17]作为一种高效的结构非线性分析新方法，自提出以来，被广泛应用到结构振动控制[18]、消能减震[19]、抗震性能评估[20]等领域。之后，李钢等[21]将其形变分解思想扩展到材料层面，与有限元相结合，提出了基于FAM的精细化有限元非线性分析方法。在求解过程中结构整体刚度矩阵始终保持不变，非线性信息通过局部塑性矩阵予以体现，迭代计算时避免了传统变刚度方法中结构整体刚度矩阵的实时合成与分解，且局部塑性矩阵相比传统变刚度法的结构整体刚度矩阵规模小的多，可有效改善传统变刚度法的计算效率问题[22]。

本文提出了基于拟力法的IDA计算方法，建立了基于FAM的IDA运动方程，考虑运用FAM非线性求解时整体刚度矩阵始终保持不变的特性以及迭代求解平衡方程时步步更新位移增量的特点，实现特定刚度矩阵存储与调用，并给出本文方法与传统IDA方法的时间复杂度函数，通过数值算例验证本文方法的准确性和高效性。

1 基本理论

1.1 增量动力分析方法

IDA法是一种基于非线性动力时程分析的参数化分析法，其实质是将单一强度的动力时程分析扩展成为多重强度的动力时程分析，以获取结构在不同强度地震动作用下的动力响应。

图1为IDA方法的基本流程图，即按照一定的调幅法则对待施加的一条(或多条)地震动记录的地震动强度(Intensity Measure,IM)进行调幅，使之成为涵盖结构从弹性、弹塑性直至动力失稳倒塌所需的多重强度；然后对给定的结构模型在多重强度地震动作用下分别进行非线性动力时程分析，得到每一强度下的地震最大响应结果即结构的工程需求参数(Engineering Demand Parameters,EDP)；将每一条地震动记录下一系列的离散点(EDP,IM)点绘制于二维坐标系中，并选用合适的插值方法拟合EDP-IM关系曲线，生成一条(或多条)IDA曲线；通过在IDA曲线定义结构不同性态水准的极限状态如：正常使用(Normal Occupancy,NO)极限状态，立即使用(Immediate Occupancy,IO)极限状态，生命安全(Life Safe,LS)极限状态，防止倒塌(Collapse Prevention,CP)极限状态，并对其趋势及离散状况进行统计分析，可以较为准确的观察结构在弹塑性响应历程中强度、刚度及变形的发展变化规律，能够反映出结构在不同水准地震动作用下地震需求能力和整体抗倒塌能力。

图1 IDA方法流程图Fig.1 Flow chart of IDA method

1.2 拟力法

FAM作为结构非线性分析新方法，其核心思想是在结构非线性分析过程中引入局部塑形机制，将结构的变形进行弹塑性分解。通过在结构中设置塑性插值点来获取塑性应变分布形式，使结构的非线性状态可由塑性插值点处塑性变形予以体现；结构控制方程则由虚功原理及塑性插值点处的内力平衡条件给出。同时在给出结构控制方程时仅考虑塑性插值点处进入非线性状态时塑性插值点(对于处于弹性状态的塑性插值点处的塑性自由度为零，可进行消元处理)。对于含有n个结点位移自由度、d个塑性自由度(塑性插值点处产生的塑性自由度数)的结构体系，其非线性求解时每一迭代步中，FAM的结构控制方程可表示为

(1)

结合Newton-Raphson迭代方法，对于式(1)进一步整理，得到每一非线性迭代步的求解方程可表示为

(2)

其中，

(3)

式中：Kf为局部塑性刚度矩阵，反映结构中进入塑性状态的局部非线性信息，其规模与结构中参与计算的塑性自由度数量有关。

对于传统变刚度法每一非线性迭代步的求解方程可表示为

(4)

对比式(2)和式(4)可知，式(2)中系数矩阵项实为式(4)中结构切线刚度矩阵的逆矩阵的等效表达。由此可知，FAM通过引入Kf，不仅避免了传统变刚度法对切线刚度矩阵实时合成与分解时计算效率低的问题，而且保留了Newton-Raphson迭代收敛速度快的特性。式(2)即为Woodbury求逆公式，求解过程中仅需对Kf进行实时的合成与分解即可。

2 基于拟力法的增量动力分析

2.1 运动方程

结构在地震动作用下的运动方程可表示为

(5)

对于式(1)进一步展开可写成

(6)

(7)

对于式(6)以全量形式可进一步表示为

(8)

将式(8)代入式(5)，并考虑引入地震动强度调幅系数λi(i=1,2,...,n)，用以表示不同的调幅强度，即可得到多重强度下的基于FAM的IDA运动方程为

(9)

2.2 方程求解

这里以Newmark-β方法[23]为例，以增量形式对式(9)运动方程进行逐步积分求解，得到λi调幅强度下时间步长为Δtj(j=1,2,...,n表示不同的时间步长)的位移增量向量的等效平衡方程可表示为

(10)

其中，

(11)

(12)

(13)

式(13)为Woodbury求逆公式，对于式(13)可进一步表示为

ΔX=ΔX1+ΔX2

(14)

其中，

(15)

(16)

(17)

2.3 刚度矩阵存储与调用

本文方法的特定刚度矩阵存储与调用流程示意如图2所示。

图2 本文方法的特定刚度矩阵存储与调用流程图Fig.2 Flow chart of storage and call of a particular stiffness matrix of this method

3 计算效率评价

算法时间复杂度分析[24]作为一种行之有效的评价算法效率的分析方法，其将运行的算法从计算机中抽象出来，使其量化结果不受计算条件的影响，依据算法时间复杂度与算法计算效率之间的关系(算法时间复杂度越高，算法计算效率越低)，可直接用于定量的评价算法的计算效率。文献[25]中采用算法时间复杂度理论对比分析了FAM与传统变刚度方法的计算效率。基于文献[25]已有理论，本文将算法时间复杂度理论应用于IDA分析的效率评价中，考虑特定矩阵的调用特性，忽略矩阵调用等次要运算的时间复杂度，给出本文方法与传统IDA方法的算法时间复杂度函数，从数学理论的角度用以揭示本文方法较传统IDA方法在IDA分析时所具有的计算高效性。本节从IDA的实质出发，考虑非线性动力时程分析中，结点位移的求解最为耗时，且本文方法与传统IDA方法主要区别也在于此，故本节仅对结构结点位移求解部分进行算法时间复杂度分析讨论。

本节以单一地震动记录的IDA分析为例，分别对本文方法与传统IDA方法的算法时间复杂度进行分析比较。对于本文1.2节给出的结点位移自由度数为n，塑性自由度数为d的结构体系，考虑因算法的迭代次数对复杂度分析的影响，假定两种方法在计算结点位移时均采用Newton-Raphson迭代法，即保证相同的收敛速度，记IDA分析每一强度下的非线性动力时程分析所需迭代次数为Ni(i对应于IDA分析时调幅系数的下标，用以表征不同调幅的强度)。

3.1 算法时间复杂度

(1)本文方法

考虑本文方法在求解结点位移时，各系数矩阵项的特殊性质(如稀疏性、带状性、对称性、可存储性等)以及计算次序对算法时间复杂度的影响，忽略矩阵调用等次要运算的复杂度，以最合适的计算方式进行算法时间复杂度分析。则本文方法的一次非线性动力时程分析的时间复杂度，如表1所示。

表1 本文方法的时间复杂度

本文方法在单一地震动记录下总时间复杂度可表示为

(18)

(2)传统IDA方法

传统IDA方法在单一地震动记录下总时间复杂度可表示为

(19)

3.2 效率对比

(20)

从上述分析可知，di值的大小是影响两种方法的时间复杂度差异的重要参数。取ηmax=dimax/n为塑性自由度数理论最大值dimax与结点位移自由度数比值，以结点位移自由度n=500，n=1 000，n=2 000为例，分别给出如图3所示两种方法的时间复杂度随塑性自由度比值ηmax的变化的关系曲线。

从图3中可以看出，本文方法的时间复杂度随ηmax的增大而增大；时间复杂度增长速率随着结点位移自由度数的增大而逐渐增强；达到dimax时的塑性自由度比例随着结点位移自由度数的增大而减小。对于不同结点位移自由度的结构来说，当ηmax较小时，本文方法的效率优势明显；随着ηmax的逐渐增大，优势逐渐减弱；当d≥dimax，本文方法的效率将劣于传统方法。

图3 塑性自由度比值与时间复杂度关系曲线Fig.3 The relationship between the ratio of plastic degree of freedom and time complexity

事实上，随着λi不断增大，结构实现了从线弹性、弹塑性直至动力失稳倒塌的全过程，在这一过程中结构往往只有少部分截面或区域(薄弱部位)产生了非线性变形即此部位塑性插值点处的塑性自由度被激活，而这部位相对于整个结构所占比例较小，如框架结构的非线性变形一般集中于梁或柱的端部，构件的中间大部分区域并不会进入非线性状态。因此，通常情况下，结构出现di≥dimax情况的概率较小，不会影响整个IDA分析过程的计算效率评价。

4 数值算例

4.1 工程概况

本文所选结构模型为8层钢筋混凝土框架结构，抗震设防烈度Ⅷ度(0.30g),设计地震分组为第一组，场地类别为Ⅱ类，结构阻尼比为5%。结构平面布置如图4所示，楼面恒荷载标准值为5.0 kN/m2，活荷载标准值为2.0 kN/m2，屋面恒荷载为6.0 kN/m2，屋面活荷载为0.5 kN/m2。楼层层重按1.0×恒载+0.5×活载折算。首层层高3.9 m，其余各层层高3.3 m，结构总高度27 m。结构梁柱尺寸及配筋见表2。混凝土强度等级，柱、梁、楼板均为C30；梁柱纵向钢筋为HRB400，箍筋为HPB300。由于结构平面对称，选取一榀横向框架(3轴线)进行建模与分析，结构立面模型如图5所示。

图4 结构平面布置简图(单位：mm)Fig.4 Sketch of structural plane layout(unit:mm)

图5 结构立面布置简图(单位：mm)Fig.5 Sketch of structural vertical layout(unit:mm)

楼层梁截面/mm配筋/mm2中柱截面/mm配筋/mm2边柱截面/mm配筋/mm21300×6004 688(BL)6 541(ZL)700×7008 758700×70011 7842300×6004 688(BL)6 541(ZL)700×7005 024700×7005 0243300×6004 688(BL)6 541(ZL)700×7005 024700×7005 0244300×6004 688(BL)6 541(ZL)700×7005 024700×7005 0245300×6004 112600×6003 786600×6003 7866300×6004 112600×6003 786600×6003 7867300×6004 112600×6003 786600×6003 7868300×6002 518600×6003 786600×6003 786注:BL为边梁;ZL为中梁

4.2 有限元模型

分别采用有限元分析软件ABAQUS和FAM程序建立分析模型进行结构的动力非线性分析，采用迭代法，不考虑P-delta效应，两种方法的梁柱均选用纤维梁单元予以模拟，单元类型为B32。其中每个构件划分3个单元，选取的一榀横向框架结构共计168个单元，864个结点位移自由度，1 008个塑性自由度，梁柱截面纤维划分均为10×10，每根钢筋相当于1根纤维；钢筋选用弹塑性随动硬化单轴本构模型；混凝土本构的骨架曲线选用Kent-Park[28]模型，采用Yassin[29]提出的滞回规则，其中ABAQUS中材料本构模型调用清华大学开发的子程序PQ-Fiber进行定义。结合两种方法对此框架结构进行单一地震下的IDA分析。

4.3 IDA分析及对比

为便于量化地震动强度等级，本文选用PGA作为地震动强度指标，地震动如图5所示方向输入，考虑场地类型、震级、峰值加速度等因素，以FEMA P695[30]推荐的22组远场地震动记录(共44条水平分量)中随机选取一条地震动记录Northridge波(NORTHR-MUL279)并归一化处理后作为IDA分析的地震动输入，持时取20 s。根据Vamvatsikos等的建议，依据hunt & fill调幅法则进行PGA调幅，给定PGA初值为0.035g，调幅步长取0.2g，步长增量取0.05g。考虑到层间位移角与楼层层间的变形能力以及节点的转动能力有关，能够反映结构的层间位移延性和整体位移延性，选取最大层间位移角θmax作为工程需求参数以及倒塌点的判定指标。结合文献[31]给出钢筋混凝土框架结构的不同性态水准下的最大层间位移角限值，如表3所示，当给定最小Δtj下数值算法仍不收敛或结构的最大层间位移角θmax>4%时，可认为结构倒塌并结束计算。

表3 RC框架结构不同性态水准下层间位移角限值

通过IDA分析，得到该结构的工程需求参数最大层间位移角θmax与地震动强度指标PGA的关系曲线，即IDA曲线，如图6所示。

由图6可以看出，对于结构的变形而言，本文方法与传统方法的计算结果基本一致，误差0.8%～5.5%，两者的IDA曲线基本吻合；在PGA处于0～0.2g内时，此时调幅强度较小，结构处于线弹性状态，结构的动力响应与地震动强度基本成线性关系，计算结果误差在1%左右，两种方法的IDA曲线基本一致；随着调幅强度的逐渐增大，结构达到屈服，非线性状态逐步发展，在PGA处于0.2～0.6g内时，计算结果误差在3.5%左右；当继续调幅地震动强度(PGA≈0.625g后)，曲线开始表现出了较小的差异，计算结果误差在5.5%左右，究其原因是两种方法在计算过程中引入的计算累积误差的差异性以及结构动力响应与地震动强度之间显著的非线性关系等造成的；在PGA处于0.8～0.9g内时，结构的性能变化特征表现出强化现象，使得曲线的斜率出现上升；当PGA>0.9g后，结构在微小的地震动强度激励下即可引起较大的地震动力响应，此时结构的最大层间位移角θmax趋于无穷大，结构体系已达到整体动力失稳状态。结合IDA方法的最终目的是通过IDA曲线簇进行统计分析，从概率的意义上用以评估结构的抗震性能。综上所述，本文方法在计算精度方面是可予以保证的，对结构的IDA分析是合理准确的。

图6 结构的IDA曲线Fig.6 IDA curve of the structure

4.4 计算效率分析

通过4.3节IDA分析，可获取该结构在每一强度下每一非线性迭代步产生的塑性自由度数的曲线。由于篇幅有限，以调幅强度为0.635g为例，现给出该强度下每一时间增量步的塑性自由度数曲线，如图7所示。

图7 塑性自由度数曲线Fig.7 Number of plastic degrees of freedom curve

由图8可以看出，在整个时程分析的过程中，传统方法的时间复杂度均高于本文方法；本文方法产生的时间复杂度的最大值约为1.94×107，但仅仅发生在最大迭代次数所在的时间增量步内，其每一时间增量步的时间复杂度的平均值约为2.16×106，将每一时间增量步下的时间复杂度汇总可得到本文方法的总时间复杂度约为4.32×109；对于传统IDA方法产生的时间复杂度的最大值约为6.65×107，其每一时间增量步的时间复杂度的平均值约为9.29×106。两种方法时间复杂度对比可知，本文方法的时间复杂度的平均值(约为2.16×106)仅为传统IDA方法时间复杂度的平均值(约为9.29×106)的23.25%，即该强度下本文方法的总时间复杂度仅为传统IDA方法的总时间复杂度的23.25%。

图8 时间复杂度曲线Fig.8 Time complexity curve

现给出地震动记录Northridge波下的IDA的不同调幅强度的时间复杂度以及复杂度比，如表4所示。

表4 本文方法与传统IDA方法的时间复杂度对比

由表4中可知，两种方法的时间复杂度随着调幅强度的不断增大，呈现增大的趋势。当调幅强度较小时，如0.035g，0.135g，0.235g，0.36g，此时结构处于线弹性状态或仅有极少部分进入塑性状态，结构中被激活的塑性自由度数远小于结构的结点位移自由度数，此时本文方法的时间复杂度相对于传统方法的优势较为明显，仅为传统方法时间复杂度的10%以下。随着调幅强度的不断增大，结构进入塑性状态的程度加大，此时塑性自由度数不断增加，使得本文方法的时间复杂度相对于传统方法的优势逐渐缩减。当调幅强度为0.035g时，本文方法的时间复杂度量级仅为108，复杂度比约为6%；而当调幅系数为0.916 25g时，此时结构处于倒塌极限状态，时间复杂度量级为1011，复杂度比约为56%。两强度之间的时间复杂度量级相差近103之多，复杂度比相差近50%。将各个调幅强度下的时间复杂度累积求和，本文方法的总的时间复杂度(1.597 8×1011)仅约为传统IDA方法(3.854 8×1011)的41.45%。由上述分析可知，本文方法较传统IDA方法具有较明显的计算效率优势且随着调幅强度的逐渐增大而减少。