有限元软件ABAQUS在结构弹塑性地震时程中的应用

有限元分析软件ABAQUS在结构弹塑性地震时程分析中的应用摘要：ABAQUS软件是国际上公认的功能强大的通用有限元分析软件，具有丰富的材料、单元模型库，能够处理高度非线性问题，在结构抗震分析领域有着良好的推广前景。本文介绍了ABAQUS在结构弹塑性地震时程分析应用中的两个问题：（1）传统上有限元软件采用的隐式积分算法，求解过程需占用大量的存储空间，随着结构自由度数量和非线性程度的增加，会遭遇到计算瓶颈，而ABAQUS采用的显式积分算法，能否有效地解决这个问题;（2）ABAQUS中的塑性区模型能否通过输入材料的本构关系得到构件各时刻的弹塑性变形和内力，从而有效地进行结构的弹塑性性能评比。关键词：ABAQUS弹塑性时程分析隐式算法显式算法引言ABAQUS是由美国HKS公司研制开发的国际上公认的功能强大的有限元软件之一，可以分析各种固体力学、结构力学系统。特别是能够处理非常复杂的问题和模拟高度非线性问题。近些年来ABAQUS优秀的分析能力和模拟复杂系统的可靠性使得ABAQUS被各国的工业和研究机构广泛采用。需进行弹塑性地震时程分析的结构往往自由度比较多，而且在地震作用下，结构构件反复进入弹塑性受力阶段，具有高度非线性，传统上采用隐式算法求解，求解过程需反复迭代，占用大量存储空间，极易遭遇计算瓶颈。更为严重的是，当结构出现严重的刚度突变（塑性铰模型）时，可能发生即使细分步长也不能收敛的情况[1]。而ABAQUS采用的显式算法能有效地解决存储空间问题，没有计算瓶颈，且ABAQUS可采用塑性区模型通过输入材料的本构关系来进行结构的弹塑性分析，没有刚度突变问题。所以ABAQUS是目前解决大规模结构弹塑性时程分析的最佳方案之一。ABAQUS对时程分析问题的处理众所周知，结构在地震作用下的振动方程式为：（1）式中：[M]、[C]、[K]分别是结构的质量矩阵、阻尼矩阵、刚度矩阵。、、分别是结构的加速度向量、速度向量和位移向量，是地面加速度向量。结构地震反应时程分析的过程即求解振动方程的过程，ABAQUS采用直接积分法来求解振动方程，ABAQUS/Standard模块中采用隐式算法进行积分，ABAQUS/Explicit模块中采用显式算法进行积分。在弹塑性时程分析中，刚度矩阵是随时间变化的函数，选择合适的弹塑性模型是关键问题，ABAQUS可采用塑性区模型等。2.1积分算法2.1.1隐式积分算法ABAQUS/Standard应用Hilber-hughes方法对运动方程进行隐式的时间积分，Hilber-hughes是Newmark-β法的改进法。基本方程如下：（2）（3）（4）式中是数值阻尼系数，如果=0.0就是Newmark-β法。在ABAQUS中，系数γ取1/2，β取1/4，所以积分算法无条件稳定。采用隐式积分进行静力分析和弹性动力时程分析时，拟刚度矩阵只需组集一次便可多次求解，求解效率很高。且在几何非线性和弹塑性非线性分析中，积分步长可以取的比较长，求解速度较快。但拟刚度矩阵求解需消耗很大的存储空间，其大小随结构自由度增加呈几何级数增长；在求解大规模问题中，极易遭遇计算瓶颈。在处理高度非线性问题时，需不断迭代缩小步长求解，每次迭代都需重组刚度矩阵，计算代价很高。2.1.2显式积分算法ABAQUS/Explicit应用中心差分方法对运动方程进行显式时间积分，动力学平衡方程如下：（5）对时间显式积分：（6）（7）式中M是质量矩阵，Pi是所施加外力，Ii 是单元内力。求解过程如下：在当前增量步开始时，用动力学平衡方程(5)计算加速度，对加速度在时间上进行积分时采用中心差分方法，在计算速度的变化时假定加速度为常数。应用这个速度的变化值加上前一个增量步中点的速度来确定当前增量步中点的速度，速度对时间的积分并加上在增量步开始时的位移以确定增量步结束时的位移。采用显式积分进行时程分析时，为满足加速度在时间增量内保持不变的假定，积分步长必须小于模型自振周期中的最小值，积分时间步长一般是隐式步长的1/1000～1/100。以加速度为计算的基本要素，在处理实际不产生加速度的静力加载问题时，需分多步缓慢加载，耗时很长。但其无需进行矩阵求逆，所求解方程均为一元方程，无需迭代求解，每步均可保证收敛。对存储空间的消耗与单元数目成正比，资源占用率低，求解速度取决于CPU浮点计算速度，适合处理大规模问题。隐式算法与显式算法的比较显式算法是有条件稳定，积分步长依赖于模型的最高阶自然频率，而与载荷类型及持续时间无关；隐式算法是无条件稳定，积分步长的大小通常参考精度和收敛情况来决定。显式算法最显著的特点是无需进行隐式方法中所必需的整体切线刚度矩阵求逆。由于是“显式”地前推模型状态，所以不需要迭代和收敛准则。在隐式分析中，每一次迭代都需要求解大型线性方程组，这一过程需要占用相当数量的计算资源、磁盘空间和内存。显式算法的求解速度取决于CPU的浮点计算能力，而隐式算法的求解速度取决于虚拟内存与磁盘之间的I/O运算。2.2弹塑性模型描述单元（梁、柱）的常用弹塑性模型有塑性铰模型和塑性区模型。塑性铰模型将塑性变形限制发生在单元的局部区域，其余部位则永远保持弹性，这是对工程实际情况的简化，该模型的单元本构关系大多采用杆件的内力-变形关系，即与杆件截面相关的本构关系。塑性区模型则是在单元长度方向和截面内逐渐进入塑性，是对工程实际情况的更精确的模拟。该模型的单元本构关系采用材料一点的应力-应变关系，单元刚度阵则由材料一点的应力－应变关系对截面和单元长度积分得到，因而轴力和双向弯矩耦合作用的屈服准则和滞回规律得到较精确的模拟。目前高端的大型通用软件如ABAQUS、ANSYS等均有塑性区模型，虽然塑性区模型计算工作量偏大，但随着现代微机计算能力的快速提高，使得该模型得以进入实用化阶段。ABAQUS有丰富的材料库，可以提供多种材料本构关系及失效准则模型，如塑性模型有金属动力硬化模型、混凝土损伤模型等。同时为了方便用户开发自己需要的材料模型，ABAQUS提供了方便灵活的二次开放平台，用户可以根据需要编写用户子程序（UMAT）。在模拟地震作用下结构的反应时，混凝土材料选择混凝土弹塑性损伤模型，钢材选择动力硬化模型（如图1钢材双线型应力－应变关系曲线图1）。算例分析3.1隐式算法与显式算法的对比一栋10层三跨框架结构[2]，层高均为3.6m，见图2。各楼层质量取为6.21×104kg，梁柱采用Q235钢材，双线性动力强化模型，强化系数0.02，阻尼比0.02。结构构件截面参数见表1。表1构件截面参数构件类别宽度/mm高度/mm翼缘厚度/mm腹板厚度/mm中柱5757252016边柱5006502016所有梁3257002016图2结构分析模型使用ABAQUS计算得出结构的前三阶周期为0.961s、0.389s、0.1705s。考察结构在地震作用在下的响应，选用EI-centro地震记录。分别采用隐式算法和显式算法计算20s内结构在8度罕遇、9度多遇地震下的结构响应。在计算8度罕遇地震结构响应时，模型材料采用金属动力硬化弹塑性模型，在计算9度多遇地震下结构响应时，模型材料采用弹性模型。计算所用机器的CPU为P4-3.1G，内存1G。用隐式算法计算时，使用ABAQUS/Standard模块，分别采用B23单元和B21单元，参数设置：采用ABAQUS固定积分步长为0.02s。用显式算法计算时，使用ABAQUS/Explicit模块，使用B21单元，参数设置：积分步长分别采用2e-6，1e-5两种情况。图3给出了ABAQUS/Standard和ABAQUS/Explicit得到的8度罕遇地震作用下的底部剪力时程曲线的对比图（图中实线表示隐式算法结果，虚线表示显式算法结果）。表2给出的是隐式分析与显式分析的计算费用对比。(a)隐式算法采用B23单元时(b)隐式算法采用B21单元时图38度罕遇地震作用下底部剪力时程曲线表2隐式算法与显式算法的计算费用对比CUP使用率所占内存（MB）耗时(s)磁盘空间（MB）隐式算法8度罕遇地震20％～30％820MB197.5s17.2MB9度多遇地震20％～25％750MB193.6s17.2MB显式算法8度罕遇地震2e-650%~60%670MB5400s1.84MB1e-550%~60%662MB2220s1.77MB9度多遇地震2e-650%~60%665MB5360s1.77MB1e-550%~60%663MB2100s1.80MB从图3、表2可知：（1）显式算法中，积分步长越小，由其得到的结果与由隐式算法得到的结果越接近。一般地，采用系统默认的自动划分积分步长足以保证精度。（2）显式算法采用积分步长2×10e-6时，底部剪力的时程曲线与隐式的在形状和幅值上基本相同（图3-a），差别源于在分析中二者分别采用了B21单元和B23单元。实际工程运算时，基于高精度的目的，通常采用B23单元。如果ABAQUS/Standard采用B21单元，那么由两种算法得到的结果几乎没有差别（图3-b）。（3）从计算成本上看，隐式算法耗费的存储空间比较大，而CPU使用率很低。而显式算法的计算成本较低，有效地使用了CPU资源。随着自由度数量的增加，特别是在弹塑性分析时，隐式算法需要大量迭代，方程求解器对资源的需求迅速增加，因此在实际工程运算上，隐式算法的计算瓶颈在于可用存储空间的大小。（4）由于该模型的单元数较少，没有体现出显式算法在耗时上的优越性；在自由度数目不多，非线性不强烈的情况下，用隐式算法比用显式算法耗时少。为了比较规模庞大的结构在显式算法与隐式算法方面的差异，将上述算例扩展为一个20层15跨的纯钢框架模型，模型其他条件同上例。计算10s内8度罕遇地震下的结构响应。显式算法和隐式算法均采用自动划分积分步长。表3给出的是隐式算法与显式算法的计算费用对比。表3隐式算法与显式算法的计算费用对比CPU使用率所占内存（MB）耗时磁盘空间（MB）隐式算法20%～30％10243小时15分钟716显式算法60％～65％7201小时9分钟59.2从表3可以看出在结构自由度比较多时，显式算法比隐式算法耗时明显减少，存储空间也小得多，所以在处理结构自由度比较多时，宜采用显式算法。3.2塑性区模型与塑性铰模型的对比某5层办公楼采用单跨轻钢框架结构[3]，构件截面尺寸见表4，钢筋混凝土楼面，压型钢板加保温层屋面，墙体为蒸压轻质加气混凝土(ALC)板，（见图4）。第2～第4层每个节点的质量为23000kg，顶层节点质量为3410kg。阻尼系数取0.02，梁柱使用Q235钢，强化系数0.02。选用EI-centro地震记录，时间步长取0.02s，计算抗震设防烈度分别为8、9度罕遇地震下的结构响应。表4构件截面参数构件类别宽度/mm高度/mm翼缘厚度/mm腹板厚度/mm柱3005501812楼面梁300600188屋面梁25035086图4结构分析模型由于模型的自由度较少，使用ABAQUS中的隐式算法进行求解。分别采用B23单元和B21单元进行求解。采用B23单元时每根构件分别划分为1个单元和7个单元两种情况；采用B21单元时每根构件分别划分为1、2、3、7个单元四种情T1=0.976sT2=0.296sT3=0.157s况。材料采用金属动力硬化模型（KinematicPlastic）图5结构的前三个振型，积分步长取0.02s。计算得出结构的前三周期分别为0.976s，0.296s，0.157s。结构的前三个振型如图5所示。图6分别给出了由塑性铰模型[3]和塑性区模型得到的9度罕遇地震下底部剪力时程曲线，图7分别给出了由塑性铰模型[3]和塑性区模型得到9度罕遇地震下顶点位移时程曲线。8、9度罕遇地震下结构的底部剪力与层间位移曲线如图8、9所示。图10给出了由计算结果得到的进入塑性构件的材料本构关系。表5所示的是由塑性区模型与塑性铰模型[3]得到的结果对比。(a)塑性铰模型(b)塑性区模型图69度罕遇地震下的底部剪力时程曲线(a)塑性铰模型(b)塑性区模型图79度罕遇地震下的顶点位移时程曲线.(a)塑性铰模型(b)塑性区模型图88度罕遇地震下的结构底部剪力与层间位移滞回曲线(a)塑性铰模型(b)塑性区模型图99度罕遇地震下的结构底部剪力与层间位移滞回曲线表5塑性区模型与塑性铰模型的对比最大位移和剪力8度罕遇9度罕遇塑性铰模型塑性区模型差值塑性铰模型塑性区模型差值顶点最大位移/ccm19.8815.9219.9%43.7518.557.7%底部最大剪力/kkN656.1377.342.5%727.141043.6%最大层间位移/ccm7.45.328.3%14.636.3556.6%从这些图表中可看出：（1）塑性区模型的计算结果与塑性铰模型的结果对比，顶点位移和底部剪力小很多；并且随着地震最大加速度的提高，位移差距增大，底部剪力差距增加较小，原因可能是因为在塑性区模型中，在单元长度方向和截面内逐渐进入塑性，而在塑性铰模型中，采用集中塑性铰假设，塑性铰的几何长度为0。梁端截面进入塑性后即形成塑性铰，所以在地震反复运动作用下，塑性铰模型的塑性变形大些，这可以从图8和图9的底部剪力与层间位移曲线中看出来。地震运动越剧烈，在地震作用过程中，塑性变形差距越大，两个模型间的顶点位移差距越大。而底部剪力增量等于刚度矩阵乘以位移增量，虽然在9度罕遇地震下，塑性铰模型的位移增量比较大，但由于形成塑性铰时，刚度比较小，所以底部剪力的差距增加较小。（2）从图10可以看出，反演的材料应力-应变关系与金属材料的双线型滞回曲线吻合良好，ABAQUS采用的动力强化塑性模型能很好地模拟钢材等金属材料的滞回性能。所以ABAQUS采用塑性区模型符合工程实际情况。(a)8度罕遇地震(b)9度罕遇地震图10塑性区模型反演的材料应力-应变关系从图11可以看出单元对结果的影响，采用B23单元的模型，网格划分稀疏对结果影响不大；采用B21单元的模型，随着网格划分加密，时程曲