考虑受压屈曲的圆钢管杆单元等效弹塑性滞回模型

谢道清，沈 金，邓 华，张 瑞

(1.浙江大学 空间结构研究中心，杭州 310058;2.浙江大学 建筑设计研究院，杭州 310027)

网架、网壳等网格结构的弹塑性地震响应分析受到越来越多的重视，主要有两个方面的原因。首先，作为应用最广的一类大跨度结构形式，网格结构被普遍认为有相当好的抗震性能。然而近些年的震害情况表明［1］，网格结构发生各种程度破坏的情况并不少见，而定量分析这些震害发生的原因(包括设计的不合理性)就不能回避弹塑性分析的问题。其次，对于按非地震作用或仅考虑小震作用来设计杆件截面的网格结构，在中震或大震条件下会呈现怎样的动力响应特点以及破坏形态，这也必须进行弹塑性分析。

关于网格结构弹塑性地震响应方面的研究总体上还比较薄弱，这也导致了我国相关抗震设计规范对网格结构在中震和大震下的抗震措施基本上缺乏详细规定。此外在具体进行弹塑性分析时，目前的文献中杆单元普遍采用的是理想弹塑性本构模型(图2)，即拉压极限状态均为轴向屈服［2-4］。然而在中震、大震等强震作用下，网格结构中杆件在受压状态下显然发生的是屈曲破坏，因此理想弹塑性模型并不能客观地反映压杆的实际工作性态。正是出于对此基本问题的考虑，本文将对网格结构中常用规格的圆钢管杆件进行不同长细比条件下的非线性屈曲分析，重点考察这些杆单元的受压极限承载力、平衡路径以及屈曲后卸载路径，并绘制拉压往复作用下杆件轴力P和两端节点相对位移Δ的滞回曲线。进一步统计这些不同规格钢管杆单元的滞回曲线控制点与杆件长细比λ之间的关系，将P-Δ曲线转换成统一描述的应力-应变(σε)曲线，以此提出一个能够同时考虑受拉屈服和受压屈曲的圆钢管杆单元的等效弹塑性滞回模型。文中以一个球面网壳作为算例，应用该等效弹塑性滞回模型进行了网壳在罕遇地震作用下的弹塑性时程响应计算。通过与理想弹塑性模型计算结果的比较，分析了杆件应力应变反应和结构薄弱区域分布上的差异，以此来考察本文提出的等效弹塑性滞回模型的有效性。

1 杆件有限元分析

1.1 考察的钢管规格

根据网格结构常用钢管规格，选择表1所列的46根不同截面和长度的圆钢管杆件作为考察拉压滞回规律的样本杆件。钢管材质为Q235钢，屈服强度fy=235 MPa，弹性模量E=2.06e5 MPa，泊松比为 0.3。表1中样本杆件的长细比范围在40～180之间。

表1 样本杆件的长细比λ及其对应的截面和长度Tab.1 Slenderness ratios of steel circular-tube specimens and their corresponding cross sections and lengths

1.2 有限元模型

钢管杆件的计算简图如图1所示，两端铰接，受轴向荷载P作用。采用LS-DYNA软件建立钢管杆件的有限元模型，管壁采用Shell163壳单元模拟。壳单元分别沿截面周圈和杆长方向按16和80等分划分，单元总数共计1 280个。壳单元采用理想弹塑性本构模型(图2)，符合Von Mises屈服准则。受压屈曲分析时计入杆件初弯曲的影响［5］，初弯曲形状按该壳单元有限元模型的一阶特征值屈曲模态引入，最大挠度取v0=L/1 000［6］。

1.3 杆件加载和卸载

使用LS-DYNA软件进行加卸载计算。加卸载按位移控制，速度为10-3m/s。首先计算杆件受压屈曲前后的平衡路径，并获得受压临界荷载Pcr。受压屈曲后的最大加载位移为8Δcr，其中Δcr为Pcr对应的杆端相对位移。然后再考察杆件在屈曲后平衡路径上2Δcr、4Δcr、8Δcr位移点处的卸载情况，跟踪卸载路径直至受拉屈服。最后再从受拉屈服段反向卸载并重新加压至杆件屈曲，以得到完整的拉压P-Δ滞回曲线。

1.4 拉压滞回特点

在46根样本杆件中选择三根典型长细比杆件来阐述圆钢管杆单元的加载和卸载滞回曲线特点。三根杆件的规格分别为:Ф180×8(L=2.5 m，λ =41)、Ф140×4.5(L=5.0 m ，λ =104)和 Ф75.5 ×3.75(L=4.5 m，λ=177)。图3给出了三根杆件在不同加、卸载情况下的P/Py-Δ/Δcr曲线，其中Py=-Afy，A为杆件截面积。为方便对比，故图3中的纵横坐标均用相对值来表达，如横坐标Δcr/Δ=1对应的纵坐标值为Pcr/Py。

图1 钢管杆件的计算模型Fig.1 The computational model of circular-tube members

图2 理想弹塑性模型Fig.2 The ideal elasto-plastic model

(1)单向受压加载

图3(a)为三根杆件仅在压力加载下的P/Py-Δ/Δcr曲线。可以看出，小长细比(λ=41)杆件的屈曲临界荷载与极限屈服荷载之比(Pcr/Py)高达0.988 7，到达临界点前基本呈线弹性反应，但屈服后承载力较快下降。对于中等长细比(λ=104)杆件，Pcr/Py下降至0.641 6，屈曲前还基本呈线弹性反应，屈曲后承载力也较快下降，但后段的下降速度相对趋缓。对于λ高达177的杆件，Pcr/Py下降到0.278 6，屈曲前后曲线过渡平缓，临界点并不明显，且屈曲后路径下降平缓。

(2)2Δcr处卸载

图3 三根典型杆件在加卸载作用下的P/Py-Δ/Δcr曲线Fig.3 The P/Py-Δ/Δcrcurves of three circular-tube specimens under loading and unloading

在三根杆件受压屈曲后变形至2Δcr时，进行卸载并拉至屈服，然后再反向卸载和加压以形成一个拉压循环，其P/Py-Δ/Δcr滞回曲线如图3(b)所示。可以看出，屈曲后路径上的卸载曲线并不像理想弹塑性模型那样呈直线变化。对于小长细比(λ=41)杆件，卸载初期基本还按直线卸载，曲线斜率接近弹性加载斜率，但是在受拉状态后期至屈服前，曲线斜率快速减小。而大长细比(λ=177)杆件的卸载特点却相反，卸载曲线前期斜率小，在进入受拉区后斜率增大并总体呈线弹性变化。对于中等长细比(λ=104)杆件，卸载曲线在到达受拉屈服前斜率虽有所变化，但幅度很小。

当卸载至受拉屈服后，再反向卸载，此时受拉区的P/Py-Δ/Δcr曲线基本上和理想弹塑性模型一样呈线性变化。继续反向卸载进入受压状态时，屈曲前后曲线形状也基本和初始的加压曲线相同，无非是由于受拉屈服使得曲线沿横坐标向左平移了一个相对伸长值。

(3)4Δcr和 8Δcr处卸载

图3(c)和(d)分别为三根杆件在受压屈曲后路径上的4Δcr和8Δcr位移点进行卸载，并反向加载的P/Py-Δ/Δcr滞回曲线图。总体上看，滞回曲线的变化规律和2Δcr处卸载的滞回曲线相差不大。小长细比杆件卸载曲线在受压区斜率较大，在受拉区后期斜率变小。而大长细比杆件在卸载初期曲线斜率较小，进入受拉区后斜率增大。中等长细比杆件的卸载曲线斜率总体变化依然相对较小。此外，8Δcr位移点的卸载曲线斜率相对4Δcr位移点的卸载曲线变化平缓些。当杆件由受拉屈服段反向卸载至受压屈曲时，P/Py-Δ/Δcr曲线变化规律和2Δcr处卸载的滞回曲线基本一致。

2 等效弹塑性滞回模型

2.1 基本思路

从图3中可以看出，考虑受压屈曲的圆管杆单元P-Δ滞回曲线与理想弹塑性模型的主要不同之处是受压区的临界点、屈曲后路径以及屈曲后的卸载路径，而且以上三方面内容与杆件长细比λ密切相关。而在受拉区，P-Δ曲线是基本符合理想弹塑性模型的。因此为便于计算机处理，可采用如图4所示(以拉为正、压为负)的分段线性化办法来描述受压屈曲后路径和卸载路径。首先，将杆件屈曲后路径用 Δcr、2Δcr、4Δcr、8Δcr对应的坐标点来分段模拟，即图4中的3～6点。其次，考虑到不同长细比杆件屈曲后卸载曲线的斜率变化情况，将 2Δcr、4Δcr和 8Δcr处的卸载曲线按两个直线段来模拟，即通过增加如图4所示的7～9三个过渡点来定义。如果屈曲后的卸载点在3～6点之间，可相应在2、7～9点间按等比例线性插值来确定该卸载曲线的过渡点，并最终定义其卸载路径。至此，考虑受压屈曲的圆钢管杆单元滞回模型实际上就可通过图4中的1～9号点来描述。

图4 等效弹塑性滞回模型Fig.4 The equivalent elasto-plastic hysteretic model

进一步将表1所列46根样本杆件的P-Δ曲线按图4形状进行无量纲归一化处理，即纵坐标表示杆件的名义应力σ=(P/A)和fy的比值，横坐标表示杆件的名义应变ε=(Δ/L)和εy的比值，其中εy=fy/E。通过归一化，易知九个控制点中，点1和点2的坐标分别为(1，1)、(0，0)。设临界点 3 的坐标为(γ1，α1)，则点4 ～6的坐标可表示为(2γ1，α2)、(4γ1，α3)和(8γ1，α4)。根据对46根杆件的归一化滞回曲线分析，确定卸载路径上7～9过渡点的相对应变值取(1-2α1)/2、(1-8α1)/3、(1-24α1)/4比较恰当，而此三点的相对应力值分别定义为 β2、β3、β4，如图 4。

2.2 控制参数的拟合

在图4定义的圆钢管杆单元等效弹塑性滞回模型中，有八个控制参数需要确定，即 γ1、α1、α2、α3、α4、β2、β3、β4。于是在46根样本杆件的归一化滞回曲线上，确定该八个控制参数并统计其与长细比λ的关系，见图5(a)～(h)所示。可以看出，有些控制参数和长细比λ的关系为单调递减，而有些为先递减再递增，因此统一采用简单的多项式形式进行拟合，并发现三次曲线就能达到较好的拟合精度。于是，进一步可得出以上八个参数和长细比λ∈［40，180］的近似关系为:

图5 八个控制参数和长细比λ的关系Fig.5 The relationships between 8 governing parameters and slenderness ratios λ

2.3 等效弹塑性模型的校验

用表1所列样本杆件以外的多根圆钢管杆件来考察等效弹塑性滞回模型的有效性。首先采用壳单元有限元法计算出这些杆件的拉压滞回曲线，然后将计算结果和等效弹塑性滞回模型的拟合结果进行比较。结果表明，等效弹塑性模型的计算结果和有限元法的计算结果非常吻合，特别在控制点处。考虑到篇幅限制，本文仅给出了一根校验杆件Ф127×4.5(L=3.25 m，λ=75)的计算结果，见图6。

其次，临界点控制参数α1实际上与我国《钢结构设计规范》(GB50017-2003)中的压杆稳定系数φ有非常相近的含义，尽管后者是根据边缘屈服准则确定的。由于本文杆件的初始缺陷仅考虑了初弯曲，接近GB50017-2003规定的A类截面［7］，故将 α1和A类截面的稳定系数φ进行比较，如图7。可以发现，α1和φ在长细比λ∈［40，180］区域的误差在5%之内，吻合很好。

3 罕遇地震作用下的球面网壳计算

3.1 网壳模型

图8所示为一个Kiewitt型双层球面网壳屋盖。网壳跨度为60 m，矢高6 m，厚度2 m。屋盖均匀支承在周边12根10 m高的C30混凝土柱上，柱截面为1.2 m×0.8 m。网壳承受的荷载标准值为:静荷载 0.5 kN/m2(不包括结构自重);活荷载为1.0 kN/m2;温度差为±20℃;风荷载0.5 kN/m2。考虑以上非地震荷载的工况组合，采用Mstcad软件对该网壳进行满应力设计以确定杆件的截面，其中满应力设计时所选择的钢管截面规格如表1。

3.2 地震波

参照我国《建筑抗震设计规范》(GB50011-2001)，对结构进行8度罕遇地震下的弹塑性时程分析。地震加速度时程采用Big-bear波［8］，按三向输入(X为主方向，如图9)。截取了该波的前38 s记录(共1 900个点)，X向加速度最大值按8度罕遇地震取400 cm/s2，三向加速度峰值的比例调整系数为1∶0.85∶0.65，结构阻尼比取0.035。在输入Big-bear波之前的2 s内，先以渐增竖向加速度的方式施加重力荷载代表值(1.0静+0.5活)，然后再叠加上Big-bear波的加速度时程，以此同时考虑结构的重力荷载代表值和地震动的作用。

3.3 结构响应计算

采用Ansys软件进行结构在Big-bear波作用下的弹塑性时程计算。网壳杆件考虑两种材料模型:模型一为图2所示的理想弹塑性模型;模型二为图4所示的考虑压杆屈曲的等效弹塑性模型。在采用模型二计算时，利用Ansys软件的APDL语言控制计算流程，即不断记录每一个荷载步所求得的杆件变形量，然后通过重启动(restart)命令并利用前一荷载步的杆件变形量来修改材料本构模型，以此实现等效弹塑性模型加、卸载迭代求解。

(1)失效杆件分布和薄弱区域分析

从理想弹塑性模型的计算结果发现，在罕遇地震作用下网壳杆件始终处于线弹性阶段，即最大拉压应力绝对值均没有超过屈服应力fy。而等效弹塑性模型的计算结果中，有110根杆件出现了受压屈曲，并最终存在残余塑性应变(指名义应变Δ/L)。发生受压屈曲的杆件主要分布在由外至内的2～4环的上下弦层内，尤其是下弦层以及部分支座附近(如图9)，其中最大名义塑性应变(总名义应变扣除弹性应变)为-1.015×103με。进一步考察这些杆件的非地震作用控制内力，发现其中大多数杆件位于壳体边缘弯矩效应和薄膜内力效应相互抵消的区域，杆件(特别是弦杆)轴力较小，因此在满应力设计时其截面配置主要由长细比控制，截面普遍较小。

(2)杆件的应力和应变

图10和11分别给出了名义塑性应变最大的两根杆件1 191、1 201(如图9)的应力、应变时程曲线。可以看出，采用等效弹塑性模型所求得的杆件名义压应变要大于理想弹塑性模型的计算结果，两者最大值相差2倍以上。而采用等效弹塑性模型所求得的名义压应力要小于理想弹塑性模型的结果。究其原因，主要是因为等效弹塑性模型中当杆件发生受压屈曲后，其承载力和轴向刚度将大大低于不考虑屈曲的理想弹塑性模型。此外，理想弹塑性模型中构件的时程响应始终处于弹性状态，而等效弹塑性模型的结果是构件最终存在明显的残余应变和应力。

4 结论

要研究网格结构对应于中震、大震设防水准的抗震措施，就不能回避对结构进行动力弹塑性性能分析。而现代数值计算技术的快速发展，为网格结构动力弹塑性响应的精细化分析提供了有力支持。本文的目的就是试图将强震作用下网格结构中杆件受压屈曲这一基本问题尽可能细致考虑，以克服理想弹塑性模型高估了杆件受压承载力的缺点。较为有利的是，无论是受拉屈服还是受压屈曲，杆单元实际上仅需通过杆件轴力和伸长量这两个简单参数来定义其受力状态，因此形式上就可以将屈服和屈曲问题进行统一描述，而且是在构件层面上。此外，计算分析也发现尽管考察的杆件规格不同，但它们的受压屈曲路径和屈曲后卸载路径从大的趋势上看是相同的，P-Δ曲线的控制参数也主要跟长细比相关，这也为建立能同时考虑受拉屈服和受压屈曲的圆钢管杆单元等效弹塑性滞回模型提供了条件。

本文提出的等效弹塑性滞回模型与理想弹塑性模型的差别主要反映在杆件受压时的加、卸载曲线上，特别是屈曲后的卸载路径随杆件长细比不同存在一些差异，即小长细比杆件的初期卸载刚度(曲线斜率)大而后期小，而大长细比杆件正好相反，因此在等效弹塑性模型中引入了八个控制参数来加以描述，形式上略显复杂。但是，由于这八个参数仅与长细比相关，在实际计算(包括插值计算)时并不复杂，用较短的程序段就可以实现。更为重要的是，只要在常规的结构弹塑性分析中对材料的本构模型进行修改就能将杆件屈曲问题加以考虑，总体上应该是方便的。此外，与理想弹塑性模型相比，本文算例结果表明了采用等效弹塑性滞回模型确实能够有效地反映网格结构在地震作用下的受压屈曲杆件分布情况，并有助于分析结构的薄弱区域。

［1］中国建筑科学研究院编.2008年汶川地震建筑震害图片集［M］.北京:中国建筑工业出版社，2008.

［2］范 峰.空间网壳结构弹塑性地震响应及抗震性能分析［J］.哈尔滨建筑大学学报，1999，32(1):132-37.

［3］薛素铎，王建宁，曹 资，等.钢网壳弹塑性地震反应分析［J］.北京工业大学学报，2001，27(1):50-53.

［4］沈世钊，支旭东.球面网壳结构在强震下的失效机理［J］.土木工程学报，2005，38(1):11-20.

［5］ Kala Z. Sensitivity assessmentofsteelmembersunder compression［J］.Engineering Structures，2009，31:1344-1348.

［6］陈 骥.钢结构稳定-理论与设计［M］.北京:科学出版社，2001.

［7］夏志斌，姚 谏.钢结构——原理与设计［M］.北京:中国建筑工业出版社，2010.

［8］ PEER Strongmotion database［EB OL］. http://peer.berkeley.edu/products/strong_ground_motion_db.html.