郑州南站混凝土多跨连续网壳结构非线性稳定分析

周鹏飞 刘传平 吴晓涵，* 宋红召

（1.同济大学结构防灾减灾工程系，上海 200092；2.同济大学建筑设计研究院（集团）有限公司，上海 200092）

0 引言

混凝土网壳结构为一种典型的空间结构，能够充分发挥混凝土材料的受压性能，因此，网壳的稳定问题也随之出现。目前，网壳结构稳定分析方法主要分为两类，即拟壳法和离散法。拟壳法是引用等效刚度条件得出等效壳体的刚度和截面特性，使相应的连续壳单元与离散杆单元表现出相等的变形。国内外许多学者对拟壳法做了大量的研究工作。Wright［1］提出了三角形网格网壳结构的极限承载力计算公式，得出qcr/q0=0.31。胡学仁［2］分别讨论了网壳的上临界荷载和下临界荷载，并提出相应公式，得出特性系数qcr/q0=0.24。离散法是基于计算机技术发展的分析方法，主要有解析半离散、梁柱单元理论和非线性有限元理论。在梁柱单元理论中，Oran［3-4］根据传统梁柱理论推导出用稳定函数表达的框架单元刚度矩阵，其单元内力和位移均为6个，容量需求小，计算效率高，但公式复杂。在非线性有限元理论中，单元为考虑轴力影响的空间梁元［5-6］，每个单元中内力和位移分别为12个，容量需求较大，计算效率相对较低，但方法简洁明了，可以方便考虑初始缺陷和残余应力。本文采用非线性有限元理论对郑州南站雨棚结构进行稳定分析。

大型通用有限元软件ABAQUS在非线性分析方面有着巨大优势，其在混凝土结构分析领域有着广泛的应用。本文采用ABAQUS软件，并借助材料用户子程序UMAT，考虑混凝土网壳结构的几何非线性和材料非线性，对郑州南站雨棚结构进行稳定分析，希望通过本文的工作能为实际工程的建设提供参考。本文主要做了以下工作：

（1）根据规范推荐本构编写的材料用户子程序UMAT和ABAQUS自带混凝土损伤塑性模型进行对比，验证UMAT子程序的合理性

（2）分析荷载分布、初始缺陷和施工方

式等对多跨连续网壳结构的稳定影响。

1 郑州南站工程概况及其有限元模型

郑州南站是郑州“米”字形格局中的主要客站，位于郑州市东南的郑州航空港经济综合实验区境内的临港型商展物流片区，效果图如图1所示，平面布置图如图2所示。

图1 雨棚结构效果图Fig.1 Rendering of canopy

图2 雨棚平面布置图Fig.2 Layout plan of canopy

站台雨棚分布于高架落客车道南北外侧，两侧对称布置。顺轨向长97.9 m，垂轨向宽369.8 m；柱距21.5m×22.8 m，共17跨，采用清水混凝土网壳结构，拱矢高3.9 m。

雨棚柱为内切圆直径为1 000 mm的正八边形，柱内设置Q345十字型钢以增强其承载力，并在边跨采用双柱以提高抗侧移能力，如图3所示。同样，柱间主梁也为型钢混凝土构件，如图4所示。顶棚壳与顶棚梁之间采用非现浇施工方式，因此在计算时可忽略雨棚板对整体刚度的贡献。

图3 中柱与边柱结构图（单位：mm）Fig.3 Structure drawing of middle column and side column（Unit：mm）

图4 主梁结构图Fig.4 Structure drawing of main beam

在设计中，恒载约为5.5 kN/m2，雨棚结构基本雪压为0.45 kN/m2，设计使用荷载为0.5 kN/m2。在使用期间的荷载标准组合值按1.0恒+1.0活+0.7雪计算，为6.32 kN/m2。

本文采用大型通用有限元软件ABAQUS对郑州南站雨棚结构进行稳定分析，在有限元模型中，梁和柱采用B31单元，顶棚壳（本文4.4节所涉及壳体）采用S4R（四边形壳单元）和S3（三角形壳单元）。B31单元采用混凝土材料用户子程序提供材料本构，而S4R和S3单元采用软件自带的损伤塑性本构模型。梁和柱中涉及到的钢筋和型钢采用*rebar关键字赋予相应单元，顶棚壳中的钢筋采用*rebar layer关键字赋予壳单元。有限元模型如图5所示。

图5 雨棚结构有限元模型Fig.5 Finite model of canopy

2 钢筋与混凝土本构模型

2.1 钢筋本构模型[7]

钢筋采用双折线模型。纵筋（HRB400）弹性模型Es=2.0×105MPa，屈服强度标准值fyk=400 MPa，极限强度标准值fstk=540 MPa，极限拉应变εu=0.01。

2.2 混凝土损伤塑性模型[8]

ABAQUS提供的混凝土损伤塑性模型是基于Lubliner［9］和Lee［10］等提出的混凝土损伤塑性本构模型，对于后续分析中遇到的壳单元（本文4.4节中的壳单元），其可以直接使用软件提供的损伤塑性模型，无须调用混凝土材料用户子程序。壳单元混凝土的非线性阶段应力—应变关系采用中国规范推荐的本构骨架曲线［7］，即混凝土受拉时，应力上升阶段采用线性模型，下降段采用中国规范推荐的应力—应变关系模型，如图6（a）所示；混凝土受压时，在压应力达到σc0之前，混凝土采用线性模型，过σc0后采用中国规范推荐的应力—应变关系模型，如图6（b）所示。本文中εc0取为0.4倍的峰值压应变，即εc0=0.4εc，r。

图6 混凝土应力—应变曲线Fig.6 Stress-strain curve of concrete

ABAQUS中的损伤塑性模型受压时需要输入应力σc、非弹性应变和受压损伤因子Dc（不同于规范给出的dc），受拉时需要输入应力σt、非弹性应变和损伤因子Dt（不同于规范给出的dt）。及可以根据GB 50010—2010附录C计算得到。损伤因子可以按照Sidiroff的能量等价原理［11］求得，也可按照文献［12］采用的方式求得。本文采用文献［12］所述方式。

混凝土受压损伤因子按式（1）确定：

式中：βc为与比值，取0.35～0.7［12］。

混凝土受拉损伤因子按式（2）确定：

式中：βt为与比值，取0.5～0.95［12］。

根据文［13］中的数值模拟与试验结果的对比，取βc=0.5，βt=0.7。

由此，可根据GB50010—2010附录C、式（1）式（2）以及系数βc和βt计算出混凝土的受压、受拉本构参数，后续模型分析中用到的C40标准强度下的计算参数见表1和表2。

除表1和表2中的参数，ABAQUS中混凝土损伤塑性模型还需要确定：膨胀角Ψ、流动势偏移值ϵ、双轴极限抗压强度与单轴极限抗压强度比αf、拉伸子午面与压缩子午面上第二应力不变量之比Kc、黏性系数μ。5个参数用Concrete Damaged Plasticity输入。具体数值［14］见表3。

表1 混凝土受压本构参数Table 1 Constitutive parameters of concrete under compression

表2 混凝土受拉本构参数Table 2 Constitutive parameters of concrete under tension

表3 Concrete Damaged-Plasticity参数Table 3 Parameters of Concrete Damaged-Plasticity

2.3 混凝土材料用户子程序本构

在ABAQUS中，3D梁 单 元（如B31、B32和B33等梁单元）因为会产生扭转剪应力，从而软件直接定义的损伤塑性模型无法应用于3D梁单元［15］。因此，3D梁单元使用规范推荐的混凝土单轴本构模型，需要编写材料用户子程序。本文计算用到的混凝土材料用户子程序是根据《混凝土结构设计规范》（GB 50010—2010）附录C推荐的本构关系编写的，使用时用户只需提供混凝土材料的初始弹性模量、峰值抗压强度和峰值抗拉强度。

不同于3D梁单元，2D梁单元（如B21、B22和B23等梁单元）不会产生扭转剪应力，所以其能够使用通过软件直接定义的损伤塑性模型。因而，建立2层1跨的平面框架结构，分别用混凝土材料子程序和3.2节中所述材料本构，对二者顶层荷载—位移曲线进行对比，以验证本文所用材料用户子程序的合理性。模型中柱为400 mm×400 mm，梁为250 mm×500 mm，层高均为3 000 mm，跨度均为5 000 mm。混凝土强度等级为C40，钢筋为HRB400，分析时均采用标准值。柱顶固定荷载及构件配筋信息见图7。

图7 荷载及配筋Fig.7 Load and reinforcement

在对上述平面框架进行分析时，在ABAQUS中建立两个分析步，分别为static general（step1）和static risks（step2），柱顶荷载从step1创建延续到step2，参考荷载在step2开始创建。梁柱采用B21单元，采用*rebar关键字将钢筋赋予到相应单元中，但到目前为止，ABAQUS的前处理cae模块中不支持rebar关键字，所以需要在inp文件中手动添加rebar关键字到相应单元。

图8为平面框架的荷载—位移图，从图中可以看出，使用材料用户子程和损伤塑性模型的荷载位移曲线在初始阶段基本重合，到达非线性阶段后，二者的荷载位移曲线出现一定程度的偏差，但在可接受的范围，从而证实用户子程序编写的本构是合理的，并且相较损伤塑性模型应用于工程分析中更保守。

图8 平面框架荷载—位移曲线Fig.8 Load-displacement curves of plane frame

3 雨棚稳定分析

3.1 稳定分析方法

特征值屈曲分析用于预测理想弹性体的理论屈曲强度，但实际工程中，结构初始缺陷、几何非线性和材料非线性等使得结构实际屈曲荷载远低于理论屈曲荷载。尽管如此，在实际分析中通常首先采用特征值屈曲分析初步判断结构的失稳模式，并且特征值屈曲分析得到的屈曲形状可以作为非线性分析的初始缺陷。在ABAQUS中，可直接定义buckle分析步对结构进行特征值屈曲分析。

非线性屈曲分析通过逐渐增加结构的非线性静力来求得结构失稳时的临界荷载。在非线性屈曲分析中可以考虑结构的材料非线性和几何非线性，因此非线性屈曲分析得到的极限承载力会比特征值屈曲分析得到的结果低很多。

在ABAQUS中，非线性屈曲分析采取以下步骤：

（1）定义特征值屈曲分析工况，施加荷载，进行屈曲分析；

（2）得到结构的一阶特征屈曲模态，找到结果中最大位移点；

（3）将最大位移点的初始缺陷取合理值［16］，然后将此初始缺陷与一阶特征屈曲模态最大位移点相比得到比例系数a，用a乘以所有节点屈曲位移，得到模型中所有节点的初始缺陷，用具有初始缺陷的节点替换模型原节点，至此整个模型都被施加上初始缺陷；

（4）建 立static general（step1）和static risks（step2）分析步，得到结构荷载位移曲线，最终求出结构极限承载力。

影响网壳结构极限承载力的因素有很多，了解清楚这些影响因素，有助于避免网壳失稳破坏。本文主要从荷载分布、初始缺陷、顶棚壳和顶棚梁施工方式三个方面讨论雨棚结构的极限承载力。

3.2 荷载分布对网壳整体稳定性的影响

荷载分布对网壳整体稳定性的影响是不可忽略的因素，但是在实际工程，对于17跨的雨棚，其荷载分布类型较多，本文考虑两种常见工况：工况1，恒载+全跨活载，如图9（a）所示；工况2，恒载+半跨活载，如图9（b）所示。对以上两种工况进行特征值屈曲分析，分析时建立static general（step1）和buckle（step2）分析步，在step1中施加恒载，step2中施加活荷载，大小为0.972 N/mm。得到两种工况下的一阶屈曲模态及其对应特征值λ1，如图10所示。由图10结果可知，工况1和工况2的屈曲特征值分别为206.08和410.52，将线荷载转换为面荷载，可求得两种工况下的承载力特征值（Pe）分别为195.86 kN/m2和390.16 kN/m2。两种工况下结构的一阶屈曲模态很类似，都为沿跨度方向侧倾失稳模式，这与网壳跨度方向刚度比垂直跨度方向刚度小很多是相符合的。通过对结构进行特征值屈曲分析，可以初步得到结构在荷载作用下的失稳模式，但是为了考虑结构在实际情况下的稳定性，需要对结构进行非线性分析。

图9 荷载分布Fig.9 Load distribution

图10 各工况屈曲模态及特征值Fig.10 Buckling modes and eigenvalues under various loading conditions

采取4.1节中非线性屈曲分析的方法对上述两种工况下结构的非线性屈曲进行分析，在分析的时候未引入初始缺陷，其荷载-位移曲线如图11所示。

图11 工况1及工况2荷载-位移曲线Fig.11 Load-displacement curves of load-1 and load-2

通过图11，一方面可以看出，工况1的极限承载力低于工况2的极限承载力，因此，同样也可得到工况1先发生失稳；另一方面可以看出，两种工况下结构的极限承载力（Pu）较特征值屈曲分析得到的特征荷载值小很多，具体见表4。

表4 Pe及Pu值Table 4 Value of Pe and Pu

导致上述两方面出现的原因是由于特征值屈曲分析是线弹性分析，未考虑结构非线性行为，而在非线性屈曲分析中，考虑了几何非线性和材料非线性。

不论是特征值屈曲分析还是非线性屈曲分析，都可看出荷载的分布对结构的稳定性影响是较大的。

3.3 初始缺陷对网壳整体稳定性的影响

在实际网壳结构中，不可避免有各种初始缺陷，如网壳结构与理论几何模型的偏差以及构件对节点的偏心等，这些初始缺陷会极大降低网壳结构的极限承载力。因此网壳结构在进行设计时必须考虑初始缺陷的影响。本文根据一致缺陷法来考虑郑州南站网壳结构的初始缺陷，即假定结构的初始缺陷是按照最低阶屈曲模态分布的。本文取满跨活载分布时结构的最低阶屈曲模态作为结构的初始缺陷分布。在ABAQUS中用，用*imperfection关键字输入初始缺陷，分别取最大位移点初始缺陷为L/20、L/25、L/35、L/40和L/50对结构进行非线性稳定分析，其中L=21.5 m，得到结果荷载-位移曲线如图12所示。

图12 不同初始缺陷下荷载-位移曲线Fig.12 Load-displacement curves under different initial defects

从图12可以看出，在曲线上升阶段，初始缺陷越大，相同荷载下，结构的位移会更大，这是由于初始缺陷使得结构在受到外部竖向荷载时，结构内部会主要产生压应力刚度，初始缺陷越大压应力刚度就会越大（此刚度起到削弱结构刚度的作用），所以相同荷载下，初始缺陷越大位移会越大。当荷载达到一定程度，结构没有了净刚度，就会发生失稳。

为方便观察极限承载力与初始缺陷关系，将图12中每条曲线峰值与对应初始缺陷绘制成曲线，如图13所示。图13显示，当初始缺陷小于L/50时，雨棚结构的极限能力随缺陷增大降低速率较慢，当超过L/50后，极限承载力下降得较快，因此本文建议郑州南站雨棚结构施工时初始缺陷控制在L/50以内。根据文献［17］的第4.3.4节，取安全系数K=2，则在L/50初始缺陷下，雨棚容许承载力（标准值）为

图13 不同初始缺陷下极限承载力Fig.13 Ultimate bearing capacity under variable initial defects

3.4 施工方式对网壳整体稳定性的影响

顶棚壳的施工方式也会对雨棚结构整体刚度产生影响，当顶棚壳采用装配方式，其对整体刚度的贡献可以忽略，但当采用现浇时，其刚度贡献不可忽略。本节采用顶棚壳与顶棚梁共节点建模分析方法，对比有无顶棚壳以及不同厚度顶棚壳对结构整体稳定性的影响，分析时未引入初始缺陷，图14为雨棚带壳有限元模型。图15为雨棚壳厚度不同时结构的荷载-位移曲线。将极限承载力与壳厚关系绘制成曲线，如图16所示。

图14 带顶棚壳有限元模型Fig.14 Finite element model with shell elements

图15 不同壳厚下荷载-位移曲线Fig.15 Load-displacement curves of variable thickness of shell elements

图16 不同壳厚下极限承载力Fig.16 Ultimate bearing capacity of variable thickness of shell elements

根据图16，从整体来看，壳体的存在对多跨连续网壳结构的极限承载力并无明显增强作用，180 mm厚壳体的极限承载力为45.22 kN/m2，无壳体的极限承载力为42.04 kN/m2，相较之下，极限承载力仅仅增加7.6%。因此，通过图16中的极限承载力变化可知，顶棚壳采用现浇或装配方式（装配方式在模型分析中体现为壳体厚度为0）对本文的网壳整体稳定性影响相对较小。

4 结论

本文通过ABAQUS对郑州南站17跨连续网壳结构稳定性能进行了细致分析，结合相关文献，主要得出以下结论：

（1）特征值屈曲分析结果与非线性稳定分析得到的结果往往相差很大，但是可以通过特征值屈曲分析感性了解结构可能的失稳模式，要得到结构体的极限承载力需要通过非线性屈曲分析求得。对于多跨连续网壳结构，其失稳形式较为常见的是沿跨度方向的侧倾失稳，要增加其稳定性，需增加结构沿跨度方向的侧向刚度。

（2）荷载分布对多跨连续网壳结构的稳定性影响较为显著，因此在进行此类结构的设计时应充分考虑结构的最不利荷载分布。

（3）初始缺陷对多跨连续网壳结构的极限承载力的影响是显著的，因此施工时应将初始缺陷控制在一定范围内，具体到郑州南站雨棚结构，本文建议将初始缺陷控制在L/50以内。

（4）对于多跨连续网壳结构，顶棚梁和板是否现浇对结构整体稳定性影响不大，原因在于现浇板并不能为连续网壳结构提供有力的侧向刚度。

（5）初始缺陷为L/50时的容许承载力大于雨棚结构设计使用荷载的标准组合，郑州南站雨棚结构整体稳定性满足规范要求。