静动力弹塑性分析的杆系离散单元计算理论研究

许玲玲，叶继红

(1. 河南工业大学土木工程学院，郑州 450001；2. 中国矿业大学，江苏省土木工程环境灾变与结构可靠性重点实验室，徐州 221116；3. 中国矿业大学，徐州市工程结构火灾安全重点实验室，徐州 221116)

以网壳结构为代表的大跨空间结构的静动力弹塑性分析一直是研究热点[1-3]。对于多由梁、杆等长细构件组成的该类结构，有限单元法的弹塑性分析方法有塑性铰法、精细塑性铰法、塑性区法等[4-5]。前两种是以内力为基本未知量建立屈服函数，后者则是以应力作为判断弹塑性状态的依据。塑性铰法[6]作为一种简化算法，无法模拟截面上及沿杆长的塑性分布和发展，计算精度不高。塑性区法[7]将杆件划分为多个区段，在区段间的截面上考虑塑性沿截面和杆长的发展过程，一般被认为是精确算法，然而，因其计算成本高而较少应用于大型复杂结构分析。精细塑性铰法兼顾了上述两种方法的优点，即在计算成本允许的条件下，近似地考虑塑性对单元刚度的削弱。Liew 和Chen 等[8]通过稳定插值函数模拟结构的几何非线性行为，在此基础上建立了精细塑性铰模型考虑残余应力对单元刚度的影响；Kim 等[9]提出了基于位移有限元法的精细塑性铰模型，并对空间钢架的弹塑性行为进行了分析。然而，有限单元法处理弹塑性问题时，结构进入塑性后，单元刚度矩阵会成为非常量，此时平衡条件需在变形后的构形上建立，故需迭代求解非线性方程组[10]。

传统颗粒离散单元[11]已发展40 多年，主要用于解决散体的弹性力学问题，其在结构工程的研究成果集中在地震作用、冲击荷载下筋混凝土或砌体结构的破坏过程[12-13]。基于此，叶继红教授课题组[14-17]提出了适用于空间杆件结构的杆系离散单元法，并率先进行了系统研究。研究成果表明杆系离散单元法在结构的弹塑性行为模拟中具有独特优势，即几何非线性问题和动力响应的求解自动包含在颗粒的运动控制方程中，是一个自然过程；杆系离散单元法中结构的弹塑性行为仅与接触点的内力求解相关，无需改变颗粒的运动控制方程，基本的计算框架不会改变。齐念[14]首先建立了杆系离散单元塑性铰法，该法假设塑性集中在接触点处，一旦接触点的内力状态达到极限屈服面即认为接触截面完全进入塑性，其无法模拟接触截面上及沿接触单元(即等效梁)长的塑性分布和发展。为此，进一步建立了杆系离散单元纤维模型，将颗粒间接触本构模型的单根弹簧等效为分布弹簧，并将分布弹簧看作是接触截面的若干根纤维，单根纤维的受力状态用单轴的应力-应变关系进行描述。叶继红和张梅[16]在其基础上进行了完善，提出了杆系离散单元塑性区法，并推导了接触截面在三维应力状态下的弹塑性本构模型。叶继红和王佳[17]基于CPU-GPU 异构平台，构建了单元级并行、节点级并行的杆系离散单元计算框架，该并行算法加速比最高可达12.7 倍。

然而，分析上述文献发现，对于杆系离散单元法的研究仍存在如下不足：1)现有文献均假定切向弹簧仅用于描述纯剪力引起的纯剪切变形，这与杆件结构常见的受力特征相违背，原因是杆件结构通常长细比较大，可忽略剪切变形的影响，即根据弯曲梁理论可认为切向位移(即挠度)是由剪力产生的弯曲变形引起的，并不是由剪力产生的截面剪切变形引起的。故基于该假定推导出的切向接触刚度系数无法用于杆件结构问题的求解；2)杆系离散单元纤维模型和杆系离散单元塑性区法应用于大型复杂结构的动力弹塑性分析时较为耗时。

本文首先重新定义了杆系离散单元法中接触本构模型的切向弹簧，并严谨推导了面向梁单元的各方向上接触单元刚度系数计算公式。其次，将精细塑性铰法引入杆系离散单元法，借鉴其采用切线模量Et的概念近似考虑由残余应力引起的沿接触单元(即等效梁)长的逐渐屈服效应，接触截面的逐渐屈服则由刚度退化函数μ表示，丰富了杆系离散单元法的弹塑性计算理论，并用一个数值算例验证了所提弹塑性计算理论的优势和精确性。最后，采用杆系离散单元精细塑性铰法对一个1/3.5 缩尺比的单层球面网壳振动台试验进行弹塑性分析，将其与试验结果进行对比，进一步验证本文方法的精度和稳定性。

1 杆系离散单元计算理论的修正

任意空间网壳结构的杆系离散单元模型如图1所示。结构被离散为一组有限的刚性颗粒，颗粒是定义结构的位置坐标、质量、边界条件以及变形等的载体。颗粒间的连接即是接触单元，根据接触单元性质建立的接触力与接触位移的关系即是颗粒间的接触本构模型。因此，杆系离散单元的基本理论主要包括颗粒运动方程的建立与求解、接触本构模型的构建两大部分，可见文献[14 - 16]。所有颗粒的运动遵循牛顿第二定律，则任一颗粒的运动控制方程为：

图1 空间网壳结构离散模型Fig. 1 Discrete element model of a spatial reticulated shell structure

式中：M、I分别为整体坐标系下颗粒的质量和转动惯量；u¨ 、 θ¨分别为整体坐标系下颗粒的平动加速度和转动加速度；F、M分别为整体坐标系下颗粒所受的力和力矩，包括内力和外力。

颗粒所受内力由与之发生接触的颗粒构成的接触单元提供，然后，接触单元处的内力可通过力系平移定理转移至该颗粒形心，而接触点处的内力增量则由接触本构模型计算[14]，其表达式为：

式中： Δf、 Δm为接触点处内力和内力矩增量；K、Kθ分别为接触本构模型的平动和转动刚度矩阵； Δu、 Δθ分别为局部坐标系下接触点处的纯平动位移和转动位移，可通过刚体运动学求得，见文献[14]。

图2 为图1 中任意两相邻颗粒A和B间的接触本构模型，该接触本构模型通过6 个弹簧表征颗粒间的接触力和接触力矩分量，每个弹簧对应一个力或力矩类型，弹簧的变形即是接触点C的位移和转角增量。因此，如何合理而精确地推导各变量的接触刚度系数即式(2)中的K和Kθ是保证杆系离散单元计算精度的关键所在。

图2 面向梁单元的接触本构模型构成Fig. 2 Beam element-oriented construction of contact constitutive model

现有文献中均假定接触本构模型中的切向弹簧仅用于描述纯剪力引起的纯剪切变形，这与杆件结构常见的受力特征相违背，故基于该假定推导出的切向接触刚度系数无法用于杆件结构问题的求解，而现有文献给出的处理方式是直接假定一个数值，然后用大量算例进行验证。简而言之，已有成果中关于切向接触刚度系数有如下困境，详细推导得到的计算公式无法使用，而采用的切向接触刚度系数又缺乏系统的推导过程。

本文以任一平面等截面直杆为例给出接触刚度系数的推导过程。该直杆的截面构造和材料特性为：截面积A、惯性矩I、弹性模量E、剪切模量G、截面切应力不均匀系数k。将直杆离散为单排刚性颗粒，取模型中任意两相邻颗粒A和颗粒B，且把颗粒A和B间的接触单元看作一根等效梁，该等效梁的长度L为两颗粒形心的距离，截面构造和材料特性与直杆一致。已有文献[14]通过设置3 根独立的弹簧表征相邻颗粒间的轴力、剪力和弯矩，对应的弹簧类型分别为法向弹簧、切向弹簧和转动弹簧，如图3 所示。其中，法向弹簧和转动弹簧的刚度系数由公式推导获得，可精确地描述颗粒间单向受拉压变形(图3(a))和由相对转动引起的纯弯曲行为(图3(b))。切向弹簧描述的是由剪力引起颗粒间的纯剪切变形(图3(c))，然而，杆件结构的长细比较大，通常可忽略截面剪切变形的影响，即认为切向位移是由剪力产生的弯曲变形引起的，并不是由剪力产生的剪切变形引起的。在上述描述下，文献[14]推导得到的法向接触刚度系数和转动接触刚度系数的计算公式分别为Kθ=EI/L和Kn=EA/L，该部分的假定和推导是合理的，后续本文也会采用。然而，文献[14]推导得到的切向弹簧刚度系数为Kτ=GA/kL，其假定和推导过程不合理，因此对该系数求解的修正即是本文创新点之一。

图3 已有研究中的弹簧等效模型[14]Fig. 3 The spring equivalent model in existing research[14]

本文法向弹簧和转动弹簧的定义和推导过程与已有研究[14]完全一致；对于切向弹簧，其变形包括由剪力产生的弯曲变形引起的切向位移和由剪力产生的剪切变形引起的附加切向位移，如图4 所示。

图4 等效模型中切向弹簧的定义Fig. 4 New definition of tangential spring in the equivalent model

1.1 切向接触刚度系数推导

假设一个时步 Δt内等效梁的相关物理量为：局部坐标系下接触点C的法向位移增量 Δun和切向位移增量 Δuτ；接触点C的转角增量即颗粒间相对转角 Δθ ；接触点C的接触力增量 ΔFn和 ΔFτ，以及接触力矩增量 ΔM；颗粒A和颗粒B的法向位移增量和切向位移增量分别为ui和wi，i=A,B；颗粒A和颗粒B的转角增量为 θi。根据杆系离散单元的基本假定可得接触点C的平动位移和转角增量与两颗粒间的关系为：

如图4 所示，等效模型中切向弹簧考虑了剪切变形影响，即认为切向弹簧的变形包括由剪力产生的弯曲变形引起的切向位移wb和由剪力产生的剪切变形引起的附加切向位移ws，即等效梁上任意点的切向位移w=wb+ws。此时等效梁截面转角、剪切变形和曲率为：

切向弹簧的应变能包括弯曲应变能和剪切应变能两个部分，其表达式为：

式中，k为截面切应力不均匀系数即截面剪切校正因子。

对式(5)的第一项即由剪力产生的弯曲变形在切线方向上引起的等效梁应变能，其的求解依据杆件结构力学中的不考虑剪切变形时弯曲梁单元相关知识[10]，则此时等效梁上任意点的切向位移即挠度函数的插值为：函数。

将式(6)和式(7)代入式(5)可得：

再由平衡条件可得，通过剪切变形得到的剪力与弯矩求导得到的剪力相等，即式(11)等于式(12)，得：

加之，几何关系：

式中，b为剪切变形影响系数，取值如式(14)。

由上述推导过程可见，杆系离散单元法会使用到能量等效原理(即接触单元或等效梁的应变能等于弹簧的变形能)，在计算等效梁的应变能时也会涉及到位移形函数，但这里的位移形函数仅用于推导各方向接触刚度系数，该系数被用于计算相邻颗粒间接触点的接触内力，即位移形函数并不会在杆系离散单元计算程序中。而有限单元法引入位移形函数，是为了利用变分原理得到弱形式的等效积分方程即求解连续体的控制方程。因此，位移形函数的使用亦是杆系离散单元法与有限元法的本质区别之一。此外，式(20)中各方向接触刚度系数的推导过程确保了杆系离散元法求解杆件结构力学行为问题的正确性和精确性。

最后，对式(20)进行扩展，得到面向空间梁单元的杆系离散单元中弹性接触刚度系数矩阵为：

值得说明的是：本文研究对象为单层球面网壳结构，该类结构杆件的长细比较大，通常会忽略剪切变形，即后续理论和算例中均令b=0。

2 杆系离散单元精细塑性铰模型

本文将精细塑性铰模型引入到三维杆系离散元中，即借鉴CRC 切线模量Et的概念近似考虑由残余应力引起的沿杆长的逐渐屈服效应，而截面的逐渐屈服则由抛物线函数表示。在现有数值方法中，最常采用精细塑性铰模型考虑材料非线性的是非线性有限元法[8]，具体的使用步骤为：1) 通过稳定函数考虑单元二阶效应，由此得到考虑几何非线性的空间梁单元弹性增量力-变形关系，可见式(22)；2) 在弹性增量力-变形关系的基础上引入精细化塑性铰模型以得到空间梁单元弹塑性增量力-变形关系，可见式(23)。

借用上述非线性有限元法的思路，文中将精细塑性铰模型引入三维杆系离散元法以考虑材料非线性。与非线性有限元法相比，杆系离散元法无需特殊处理即可考虑几何非线性，故直接在弹性接触本构模型即式(21)的基础上引入精细塑性铰模型，进而建立适用于空间杆系结构的弹塑性接触本构模型。对比式(21)与式(22)知，两种方法单元刚度矩阵的表达式相似，但有以下两点区别：1)三维杆系离散单元即式(21)无需特别考虑轴力与弯矩的耦合效应，故在杆系离散元中认为s1y=1.0,s1z=1.0,s2y=0,s2z=0；2) 式(22)对应的是节点变形，而式(21)对应的是接触点(图2中C点)变形，故在杆系离散元中需将式(23)中的单元两端截面刚度退化函数修正为接触点处截面刚度退化函数。最后将上述参数取值代入式(23)并修正刚度退化函数，得到适用于空间杆系结构的杆系离散元弹塑性接触本构模型的表达式为：

3 大型网壳结构算例验证与分析

3.1 K6 型单层网壳结构

K6 型单层网壳结构的几何构造如图5 所示。跨度和矢高分别为30 m 和2 m，杆件均为φ180 mm×5 mm 的圆钢管，弹性模型E=210 GPa，屈服应力σy=240 MPa ，剪切模量G=850 GPa，材料为理想弹塑性。底部支座为铰接，结构仅受顶部节点A处的竖向集中荷载作用。文献[20]结合了塑性铰和塑性区法的基本概念，即认为单元弹塑性仅产生在杆端截面，杆件其他部位仍为弹性，同时将杆端截面划分成若干小面积模拟截面的塑性发展。文献[16]采用杆系离散元塑性区法对该结构进行了弹塑性分析，每个接触截面被划分为16 个小面积。通用有限元软件ANSYS 中，每根杆件被划分为4 个beam189 单元，采用弧长法进行位移追踪，且考虑残余应力的影响。

图5 K6 型单层网壳结构Fig. 5 K6 single-layer reticulated shell structure

为便于文献对比，文中杆系离散单元将结构离散为529 个颗粒，每个杆件有4 个接触单元，计算时步为10-4s，阻尼系数为10。采用杆系离散单元位移控制法进行加载，位移控制速度为1 mm/s。材料非线性选取Orbison[18]提出的截面极限屈服函数作为接触截面的屈服面方程，切线模量为CRC切线模量即式(25)，截面的退化刚度系数的取值按式(26)。

杆系离散单元得到的节点A和节点B的荷载-竖向位移曲线如图6 和图7 所示。由图可见，文中杆系离散单元得到的荷载-位移曲线与已有文献、ANSYS 计算结果在弹性阶段均吻合较好；相比于已有文献[20]，杆系离散单元得到的3 条荷载-位移曲线在塑性阶段与ANSYS 计算结果更吻合；当杆件进入塑性后，塑性铰法的计算结果明显大于塑性区法[16]和ANSYS，与两者的最大误差分别为9.09%和14.78%，这是因为塑性铰法无法考虑塑性在接触截面上和沿接触单元(即等效梁)长的分布和发展，使得结构的整体刚度偏大；相比于塑性铰法，杆系离散单元精细塑性铰法与ANSYS 计算结果和已有文献[16]在塑性阶段吻合更好，且与文献[16]的最大误差仅为2.2%，表明精细塑性铰法的计算精度较塑性法有显著提高。

图6 节点A 的荷载-位移曲线Fig. 6 Load-displacement curve of node A

图7 节点B 的荷载-位移曲线Fig. 7 Load-displacement curve of node B

杆系离散单元塑性铰法、本文精细塑性铰法以及塑性区法各自的计算耗时如表1 所示。由表可见，杆系离散单元精细塑性铰法和塑性区法的计算时间分别是塑性铰法的1.21 倍和2.88 倍，可见精细塑性铰法未显著增加杆系离散单元的计算量；与杆系离散单元塑性区法[16]相比，在误差可接受的情况下计算耗时被大大减小，可以推测在大型结构弹塑性分析时，计算耗时的对比会更加明显。因此，沿杆长塑性分布不明显时，采用杆系离散单元精细塑性铰法求解最为划算。

表1 杆系离散元中三种弹塑性分析方法的计算效率对比Table 1 Comparison of the calculation efficiency of three elastoplastic analysis methods using MDEM

3.2 大型单层球面网壳振动台试验

本次振动台试验由同济大学多点振动台试验系统完成，该系统由4 个振动台组成，试验模型底圈共设40 个支座，平均落在4 个振动台上，其余节点处于悬空状态，如图8 所示。试验模型的缩尺比为1/3.5，结构类型为K6 型单层球面网壳，试验模型严格按照原型的拓扑关系，不作任何简化。试验模型的跨度为23.4 m，矢跨比为0.5，焊接空心球数为1261，杆件数为3660，模型节点的配重为30 kg。杆件截面共有4 种规格，且满足满应力设计准则，其中，截面尺寸为φ23 mm×1 mm 、φ38 mm×2 mm 的杆件为采用经过真空退火处理的20#钢，以模拟原型材料Q235 的应力-应变关系，其余杆件采用Q235 热轧钢管，每种规格杆件的各项力学指标详见文献[21 - 22]。

图8 振动台布置与试验模型Fig. 8 Shaking table layout and test model

试验采用EI-Centro 地震波(1940)，根据模型相似比将该波的持时压缩为28 s，时间间隔为0.010 69 s。每个振动台实际输入的地震波加速度时程曲线和速度时程曲线、以及加载工况详见文献[21 - 22]。

试验模型被离散为14 581 个颗粒和16 980 个接触单元，得到的杆系离散模型如图9 所示。采用位移法施加多点激励的地震作用，具体加载的步骤和公式见文献[22]。结构的阻尼比和基频根据白噪声扫频结果计算得到，数值如图10 所示。材料本构为理想弹塑性模型，仍选取Orbison[18]提出的截面极限屈服函数作为接触截面的屈服面方程。

图9 试验模型的杆系离散单元数值模型Fig. 9 Discrete element model of test model

图10 试验模型的阻尼比和基频Fig. 10 Damping ratio and fundamental frequency of test model

图11 给出了试验模型的PGA-Δmax曲线对比(荷载-位移选取点见图9)，由图可见，数值仿真结果和试验结果的变化趋势相同，即随着峰值加速度的提高，最大位移显著增加，相应倒塌的前一级加载工况，模型节点最大位移达到70.71 mm(试验为78.76 mm)，为跨度的1 /330，倒塌征兆显著，PGA-Δmax曲线斜率逐渐减小，未有突变现象，结构表现出良好的延性，塑性发展充分，说明结构的破坏属于承载力破坏。此外，从图11 的对比结果亦可得出与算例3.1 一致的结论，即当结构进入动力弹塑性行为后，杆系离散元精细塑性铰法的计算精度明显高于塑性铰法，进一步验证了前者的可靠性。

图11 PGA-Δmax 曲线对比Fig. 11 Comparison of PGA-Δmax curves

将试验模型在不同工况下位移时程曲线与试验结果进行对比如图12～图13。由图可见，杆系离散单元精细塑性铰法计算得到的位移时程反应曲线与试验结果几乎重合，仅在峰值上有些许误差，可能原因是：1)试验中各工况下使用模型的初始状态为前一级工况下变形后的试验模型，而数值仿真中各工况下采用的计算结构均为未变形的原始模型；2)数值仿真中未考虑地震荷载下钢材应变率效应的影响，使得部分杆件的屈服强度低于实际值，使得数值计算位移大于试验值，这一现象在变形大的工况较为明显，如图12(b)、图12(c)所示。

图12 不同工况X 向最大位移测点的时程曲线位移对比Fig. 12 Comparison of displacement time history curves of the maximum displacement of measuring point in X-direction in different cases

图13 不同工况Y 向最大位移测点的时程曲线位移对比Fig. 13 Comparison of displacement time history curves of the maximum displacement of measuring point in Y-direction in different cases

图14 给出了400 gal 工况下测点的位移时程曲线。由图可知：数值仿真得到的模型极限荷载为400 gal，与试验结果一致，但二者完全倒塌时间点略有差别，试验中完全倒塌时刻为6.62 s，数值模拟中为5.69 s，误差较小。

图14 400 gal 工况下测点的位移时程曲线Fig. 14 Displacement time history curve of a measuring point in the case of (PGA=400 gal)

图15 为400 gal 极限荷载下试验模型破坏过程的数值仿真结果，该仿真结果与试验结果[21-22]相一致。由图可知：随着输入峰值加速度的提高，各台面上方第2 圈、第3 圈斜杆首先发生弯曲；之后，结构整体变形程度逐渐加大、破坏区域逐渐扩展，同时与地震波输入方向垂直的台面之间第3 圈～第5 圈杆件变形明显，但上部仍处于弹性状态，无明显变形；最终，模型支座上方及与地震波输入方向垂直的台面之间第1 圈～第6 圈杆件均变形显著，整个结构向下坍塌。可见，试验模型的杆件屈曲过程是由下而上、由各台面上方破坏区域向两侧扩展；对于正常设计的网壳结构，支座上方杆件和与地震波输入方向垂直的主肋区下部杆件及其附近区域为结构关键部位。

图15 试验模型的破坏过程数值仿真结果Fig. 15 Numerical simulation results of the failure process of the test model

4 结论

本文对杆系离散单元法的基本计算理论进行了修正，并将精细塑性铰法引入其中，最后通过两个大型网壳结构的弹塑性行为验证了所提方法的精确性，得到如下结论：

(1)重新定义了杆系离散单元中接触本构模型的切向弹簧，并根据能量等效原理推导了不考虑剪切变形和考虑剪切变形时切向接触刚度系数的统一计算公式，实现了杆系离散单元基本理论系统化，为结构力学行为模拟提供更为严谨的理论支撑。

(2)提出了杆系离散单元精细塑性铰法，并推导了精细塑性铰法的弹塑性接触本构模型。其通过切线模量和截面刚度退化系数近似考虑残余应力对接触单元刚度的削弱，前者考虑了沿接触单元长的逐渐屈服效应；后者表示接触截面的逐渐屈服。静力弹塑性分析算例表明，杆系离散单元精细塑性铰法可近似考虑构件的塑性发展，其计算精度明显高于塑性铰法，且不会显著增加杆系离散单元的计算量；当截面分布塑性不明显时，相比于塑性区法，采用杆系离散单元精细塑性铰法“性价比”更高。

(3)将修正后的杆系离散单元法应用于多点激励单层球面网壳弹塑性行为模拟，并与试验结果对比，分析结果表明，二者极限荷载、倒塌模式一致，且数值仿真得到的位移时程曲线与试验值几乎重合，仅在幅值上有微小差异，进一步验证了本文方法处理大跨空间结构动力弹塑性问题的精确性。