考虑材料非线性的RC 双曲冷却塔风致破坏过程

张军锋，裴 昊，朱 冰，刘庆帅

(郑州大学土木工程学院，郑州 450001)

RC 双曲冷却塔作为一种高耸空间薄壁结构，风荷载是其设计控制荷载[1]。渡桥电厂冷却塔风毁事故[2]和不断攀升的结构高度使结构的抗风安全性始终备受关注，既有研究主要关注结构的风荷载特性、群塔干扰效应、风振效应以及弹性稳定问题[3]。实际上，渡桥电厂事故被一些学者认为是风致失稳破坏，从而引发了近四十年的稳定研究热潮[4 − 8]。然而，Mang 等[9]于1983 年采用有限元方法对Port Gibson 冷却塔首次进行了考虑材料非线性的极限荷载分析，并认为其风致破坏应属于材料破坏而非弹性失稳。之后，冷却塔的静风极限荷载分析受到关注。既有研究基本均对塔筒采用分层壳单元，沿厚度划分为多个层次并用钢板层模拟钢筋。对于荷载主要考虑自重和风荷载[9 − 13]，部分学者亦考虑了温度效应[14 − 15]。但由于问题本身的复杂性以及各学者所采用的软件、单元、材料本构、计算方法以及荷载取值等的差异，即使对Port Gibson 同一冷却塔的计算结果也表现出较大的离散性[16]：所得开裂荷载系数和极限荷载系数分别分布在0.75～1.58 和1.49～3.15，并且其中亦有文献[11]认为是失稳破坏。

近年来，随着国内冷却塔建设规模的快速发展，国内学者亦注意到对冷却塔应考虑材料非线性进行分析。但初期研究往往比较粗略，未能对混凝土的材料非线性进行准确模拟[17 − 18]。文献[19]进行了考虑初始缺陷的极限风荷载分析，但未就关键的材料非线性模型和模拟方法进行介绍。另外，对塔筒破坏的原因也有不同的结论：文献[20]关注了动力极限风荷载，但仅对塔筒的受压损伤进行分析而对更关键的受拉损伤却没有说明，并认为迎风区混凝土的破坏是受压破坏；文献[21]模拟了冷却塔的倒塌过程，其结果显示冷却塔破坏源于混凝土受拉开裂和环向钢筋断裂引起的材料强度失效，这均与大多认为冷却塔源于迎风区子午向受拉破坏的结论[9 − 10, 12]不一致。鉴于塔筒RC 材料非线性模拟的困难，文献[22 − 23]将塔筒视为均质材料并赋予弹塑性材料特性，但这种处理无法准确模拟混凝土开裂后的特性以及钢筋和混凝土之间的协同工作特点，且两者对同一冷却塔所取的弹塑性参数也差异显著。

以一座大型冷却塔为例，采用ABAQUS 软件对其进行极限静风荷载计算。采用分层壳单元模拟塔筒，采用弥撒开裂模型模拟混凝土的非线性受压特性、受拉开裂和拉伸硬化效应，并计入几何非线性。从荷载-位移曲线、变形模式、内力、应变和应力的发展对塔筒的破坏过程进行系统分析，并从中探究其破坏机理。

1 材料非线性模型

1.1 混凝土弥散开裂模型

钢筋混凝土的承载特性表现出高度非线性，合理地反映这些非线性因素是准确获得冷却塔破坏过程的必要前提。本文采用弥散开裂模型模拟混凝土非线性本构，该模型可以准确地模拟混凝土构件受弯时的非线性受力特性[24]，也是既有冷却塔极限风荷载分析的主要方式[9 − 15]。

该模型通过4 个方面的材料行为模拟混凝土的非线性特性：非线性受压特性，破坏面，剪切保留和拉伸硬化，详见文献[16]。

剪切保留模型如式(1)[25]和图1[13]所示。

式中：G0为闭合混凝土的剪切模量；G'为开裂混凝土的剪切模量；β 为裂缝闭合后的剪切刚度系数；εt、εcr和εm分别为混凝土开裂应变、开裂后垂直于裂缝方向的应变和剪切刚度为0 时的应变。参考文献[12 − 13]，下文取β=0.25，εm=0.004。

拉伸硬化效应[10 − 12]用以模拟混凝土开裂后的行为[24]。与文献[10 − 11]一致，下文同样假定混凝土开裂后应力线性卸载(图2)，在应变达到极限应变εu=TSεt时应力降至0，TS 即为拉伸硬化参数。既有研究[10 − 12]表明TS 对结果影响较大，在板壳结构试验结果验证模型的基础上[16]，下文取TS=15。文献[16]还分析了TS 对结果的影响。

图1 剪切保留模型Fig.1 Shear retention model

图2 拉伸硬化模型Fig.2 Tension stiffening model

1.2 钢筋非线性模型

尽管钢筋被模拟为平面应力钢板层，但其实际是单向受力，故采用单轴双线性应力-应变曲线模拟其非线性本构[9,12,15]，如图3 所示。其中：σy和εy分别为屈服应力和屈服应变；Es和Ep分别为弹性模量和塑性模量，下文取Ep=0.05Es。

图3 钢筋双线性模型Fig.3 Bilinear model for reinforcement

2 结构及有限元模型

本文所选冷却塔如图4 所示，塔高177 m。下文采用塔筒相对高度hs/Hs表示所关注的高度位置hs，其中Hs为塔筒总高度(图4)。塔筒为渐变厚度(图5(a))，最小厚度0.271 m，塔筒下缘最大厚度1.40 m，上缘厚度渐增为0.40 m，塔筒通过48 对Φ1.30 m 人字柱与基础连接。塔筒和顶端檐口均采用C40 混凝土和HRB400 级钢筋，塔筒环向和子午向配筋率如图5(b)所示，下支柱采用C45 混凝土。

图4 冷却塔特征尺度示意Fig.4 Characteristic dimensions of cooling tower

图5 塔筒结构参数Fig.5 Structural parameters of tower shell

依规范[26]取基本风速V0=31 m/s，场地类别为B 类，对应风振系数β=1.9。温度和风荷载联合作用下的环向弯矩亦需重视[27]，考虑到冬温和夏温效应接近[28]，故仅考虑冬温作用：依规范[26]塔外设计气温取−9.6 ℃，塔内设计气温取(环梁有挡水设施+单元系统)模式。

由于结构和荷载均对称，为降低计算代价，仅建立半结构模型[9− 12,14]。塔筒采用壳单元S4R模拟，经不同网格密度划分对比将塔筒划分为环向×子午向=96×72 个单元，并沿厚度划分9 个混凝土层加4 个单向受力钢板层(图6)。

图6 分层壳示意图Fig.6 Layered shell element

其中，最外层(也即图6 中的第1 和第9 层)厚30 mm 用以模拟混凝土保护层，其它各层按1∶1∶1.5∶2∶1.5∶1∶1 确定。不同网格和层次划分的对比分析详见文献[16]。钢板层厚度根据配筋面积相等进行确定，其位置均在相应钢筋的中心位置。下支柱和顶端檐口均采用梁单元B31模拟。

考虑两种荷载工况并采用牛顿迭代法进行计算：自重+风荷载(G+λW)和自重+冬温+风荷载(G+T+λW)，其中G、T 和W 分别为规范[26]所列荷载标准值。对于荷载效应，轴力F 以受拉为正，弯矩M 以外表面受拉为正；轴力和弯矩均以X 方向为环向，Y 方向为子午向，如FY和MX分别为子午向轴力和环向弯矩，具体可参考文献[30]。径向位移以内凹为正，外凸为负。环向角度θ 以迎风点为0°，逆时针转动为正。

需要说明，在对本算例进行分析前，还首先根据文献[10, 12]对Port Gibson 冷却塔进行了重复计算[16]。尽管所用软件和单元不同，材料非线性模型和计算方法等也存在差异，但所得结果基本吻合(图7)，从而验证了本文所用软件、单元以及材料非线性模型和计算方法的合理性。

图7 Port Gibson 塔的荷载位移曲线Fig.7 Load displacement curves of Port Gibson tower

3 线性结果分析

作为非线性分析的基础，首先应了解冷却塔在各荷载作用下的内力特征。自重作用下(图8)塔筒以双向受压为主，且子午向压力远大于环向压力，双向弯矩可以忽略。自重产生的轴压力将有助于抵抗风荷载产生的拉力。冬温荷载作用下(图9)塔筒以双向正弯矩为主，双向轴力可以忽略。

图8 自重作用下塔筒内力分布Fig.8 Internal forces of shell under gravity

图9 冬温荷载作用下塔筒内力分布Fig.9 Internal forces of tower shell under gravity

风荷载作用下塔筒内力分布较为复杂(图10)，除MY外，各荷载效应均较显著，并以子午向轴力最为显著。对各内力特征的详细介绍可见文献[27]。

图10 风荷载作用下塔筒内力分布Fig.10 Internal forces of tower shell under wind load

4 非线性结果分析

4.1 变形模式

在两种工况下，施加风荷载前，结构变形可以忽略。施加风荷载后，各截面的变形规律基本一致，并且考虑到对称性仅给出0°≤θ≤180°范围的结果(图11)：即迎风区和侧风区分别产生显著的内凹和外凸位移，且前者比后者增加更为显著，背风区的变形始终较小。在进入非线性阶段后，迎风区内凹位移最大值均始终在θ=0°位置，但随风荷载的增加而有上下移动；而侧风区外凸位移最大点则同时有环向和子午向移动：图11 所示从下到上的4 个断面的外凸位移最大点逐渐由θ=40°区域后移至θ=75°区域，而在子午向则逐渐由hs/Hs=0.65 区域上移至塔筒顶部。

图11 塔筒横向截面变形Fig.11 Deformed shapes of tower shell: sliced views

4.2 荷载-位移曲线

为直观体现结构的非线性效应，选取Ⅰ点(θ=0°、hs/Hs=0.782，迎风点喉部)和Ⅱ点(θ=60°、hs/Hs=0.70)这两个内凹和外凸位移始终较为显著的位置，以及Ⅲ点(θ=120°、hs/Hs=0.70)作为对照位置，给出荷载-位移曲线如图12 所示。可以看出，两种荷载工况下，冷却塔破坏过程中的非线性主要体现在迎风区和侧风区，其他区域的非线性极不明显。

图12 两种荷载工况所得荷载位-移曲线Fig.12 Load-displacement curves of two load cases

对于工况G+λW，以Ⅰ点为例，其荷载位移曲线可划分为3 个阶段。第1 阶段(OA 段，λ<1.384)为线性响应阶段：位移随风荷载线性增加，整个塔筒未出现裂缝。第2 阶段(AB 段，1.384≤λ<1.876)为裂缝扩展阶段：由于风荷载在塔筒迎风区产生的子午向轴拉力逐渐达到并超过自重作用下的子午向轴压力，故λ=1.384 时，迎风子午线上hs/Hs=0.37 位置混凝土应力首先达到峰值并进入卸载阶段，也即该处产生裂缝，且该裂缝为贯穿塔筒壁厚的环向裂缝；随着风荷载的增加，环向裂缝在塔筒迎风区不断形成并扩展，塔筒局部刚度逐渐下降，结构响应进入明显的非线性阶段，迎风区和侧风区位移迅速增大，内力也出现重分布，详见下节。第3 阶段(BC 段，1.876≤λ≤2.007)为硬化阶段：λ=1.876 时，迎风子午线上hs/Hs=0.70位置内侧子午向钢筋首先受拉屈服，结构进入硬化阶段，风荷载又有进一步增加，钢筋屈服区域沿环向逐渐扩展至θ=16°，沿子午向则扩展至hs/Hs=0.52～0.80。并且λ=1.529 时，迎风子午线hs/Hs=0.62高度内侧和θ=75°、hs/Hs=0.95 高度外侧出现子午向裂缝。最终，由于迎风区混凝土的持续开裂和钢筋屈服，λ=2.007 时计算终止(表1)。这一荷载位移曲线形式与文献[9 − 15]所得结果较为类似。

表1 两种荷载工况结果Table 1 Results of two load cases

对于工况G+T+λW，由于冬温荷载产生双向正弯矩，使塔筒外表面混凝土在风荷载作用前即处于双向受拉状态，而风荷载作用后侧风区的环向正弯矩效应进一步增加，迎风区子午向也将处于拉弯联合受力状态，所以塔筒提前开裂：λ=1.000 时，塔筒侧风区θ=73°、hs/Hs=0.93 位置外表面即首先出现子午向裂缝；λ=1.080 时，塔筒迎风子午线上hs/Hs=0.08 位置外表面开始产生环向裂缝。因此，该工况的荷载-位移曲线略提前进入非线性阶段，且整体刚度偏小。

对比表1 可知，冬温荷载显著减小了冷却塔的开裂荷载和极限位移，但对极限荷载影响不大。这是由于尽管冬温荷载产生的双向正弯矩导致塔筒提前开裂，但此弯矩效应较小且固定不变，且随着风荷载的持续增加，其与风荷载产生的荷载效应相比也越来越小，故冷却塔的最终破坏依然由风荷载控制。也正因此，两种荷载工况所得荷载-位移曲线几乎一致；二者变形和裂缝开展情况也相似(见4.2 节和4.3 节)，只是在同一λ 下，计入冬温效应的变形和裂缝开展更为严重。因此，对于本例，温度作用只是增加了环向和子午向弯曲损伤，这在文献[14 − 15]中亦有体现。

4.3 内力分布

荷载-位移曲线表现出的明显非线性特征直接源于混凝土开裂导致的结构刚度损伤，这也使结构刚度分布发生变化，进而导致内力重分布并加剧结构开裂。由于两种荷载工况所得荷载-位移曲线变形模式相差不大，为观察塔筒破坏过程中的内力重分布情况，以工况G+λW 为例，给出hs/Hs=0.37 和0.70 这两个断面的内力分布如图13 所示。

由图13 可知，在λ≤1.384 时，整个塔筒所有内力均呈线性变化，并且由于计入了自重效应，FY在迎风区逐渐由轴压力表现为显著的轴拉力，在侧风区始终为显著的轴压力，在背风区亦为轴压力但幅值较小；自重效应的存在也使FX在整个环向均为轴压力；MX在整个环向幅值较小，在迎风区和侧风区分别为负弯矩和正弯矩。

图13 塔筒内力分布Fig.13 Distribution of internal forces of tower shell

λ=1.384 时，由于hs/Hs=0.37 高度首先出现环向裂缝(见下文4.5 节)，θ=0～15°的各内力均开始表现出非线性并逐渐扩展到θ=45°位置，即出现内力重分布。并且环向裂缝的出现不仅影响子午向受力，也影响环向受力：与此同时，hs/Hs=0.70 高度θ=0°的FX和MX也开始出现内力重分布，而该高度FY依然保持线性变化，直至λ=1.48 时开始非线性变化。而此时此处尚未开裂，这同样是内力重分布的结果。此后，随着塔筒开裂的继续发展，塔筒内力重分布显著区域也不断扩展至θ=90°位置并在两个高度表现出相似的规律。

对于子午向的内力重分布，θ=0°位置开裂后，此处的FY略有下降，但其两侧的FY快速增加，用以补偿θ=0°子午线上FY的下降。迎风区大范围受拉开裂后，整个迎风区的FY不再线性增加，还造成侧风区的FY亦不再随荷载线性增加。与环向裂缝主要因子午向轴拉力引起不同，子午向裂缝主要由环向弯矩引起，内力重分布也主要表现为迎风区和侧风区MX和FX幅值和环向波动的急剧增加。

4.4 应变分布

图14 给出了两种荷载工况所得塔筒最外层和最内层混凝土最大主应变εpr分布。由于冬温和风荷载的弯矩效应有限，两种工况下，塔筒均主要依靠薄膜作用抵抗风荷载，所以塔筒内外层混凝土εpr随风荷载的变化规律基本相同，尤其是工况G+λW。对于工况G+λW，λ=1.50 时，迎风区中下部已有较大面积混凝土产生了环向裂缝，故结构开始表现出明显的非线性特征。λ=1.80 时，迎风区已产生大面积开裂，并存在4 个明显的贯穿壁厚的环向裂缝带(亦可见图15(c))；由于迎风区和侧风区分别存在较大的环向负弯矩和正弯矩，故上述区域内侧和外侧混凝土均产生了子午向裂缝，但其宽度相对较小，直到结构破坏也未贯穿壁厚。λ=2.00 时，迎风区hs/Hs=0.08～0.80 区域全部开裂，且开裂区域与λ=1.90 时近乎一致，仅裂缝宽度急速增大。最终，冷却塔因大面积的混凝土开裂和钢筋屈服而破坏。计入冬温荷载后，冬温荷载产生的双向正弯矩导致塔筒提前开裂，且同一λ 下裂缝开展情况更为严重，主要表现在侧风区外侧。

4.5 应力分布

由于迎风子午线上的子午向受拉开裂和hs/Hs=0.70 断面的环向拉弯开裂最为显著，为进一步了解塔筒破坏过程中的应力分布情况和破坏过程，给出G+λW 工况下迎风子午线上最外层混凝土和外侧钢筋的子午向应力分布和hs/Hs=0.70 断面表层混凝土和钢筋的环向应力分布，同时给出上述两个位置处的应力荷载曲线，如图15 所示。

由图15(a)可以看出，λ=0.50 时，除塔筒顶部和底部边界外，塔筒整个高度范围均处于双向受压状态。λ=1.0 时，风荷载作用下的子午向轴拉力已经超过自重的压力效应，迎风区大部分区域已进入子午向受拉状态。λ=1.50 时，hs/Hs=0.10～0.50范围内混凝土子午向应力已达到抗拉强度，且hs/Hs=0.37 处混凝土已进入卸载阶段。λ=1.80 时，迎风区已有4 个位置混凝土应力完全降至0，即迎风区产生了4 个完全开裂的环向裂缝带，裂缝区荷载完全由钢筋承担导致钢筋应力激增(图15(b))。λ=2.00 时，迎风区子午向开裂区域几乎不再随λ增加，hs/Hs=0.08～0.80 范围混凝土完全退出工作，hs/Hs=0.52～0.80 区域钢筋屈服(图15(b)、图15(c))。

图14 塔筒混凝土最大主应变(仅给出εt=7.354×10−5～5.2×10−3 范围内的拉应变)Fig.14 Maximum principal strain of concrete in tower shell

图15 塔筒外侧混凝土和钢筋应力分布Fig.15 Distribution of stress of exterior concrete and steel in tower shell

另外，从图15(c)更可直观看出，混凝土开裂后其应力急剧减小至0，与此同时裂缝区钢筋的应力则急剧增加。4 个主裂缝带产生的位置顺序依次为hs/Hs=0.37、0.52、0.17 和0.70，对应的λ 依次为1.384、1.451、1.525 和1.631。尽管喉部位置混凝土最后开裂，但由于该处壁厚和配筋均较小，故该处钢筋首先屈服。这也说明，子午向受拉开裂和钢筋屈服是冷却塔在风荷载作用下破坏的直接原因。

对于环向应力，由图15(d)可以看出，hs/Hs=0.70高度外侧混凝土环向应力的分布与图13 所示MX的分布规律一致，这也表明冷却塔的环向应力状态受弯矩控制。λ≤1.50 时，除迎风子午线局部区域，环向应力基本处于线性状态；进入非线性状态后，受环向正弯矩影响，θ=25°～65°区域外侧混凝土应力达到抗拉强度，但此后其应力并未快速卸载，而是维持在一个较高的水平，这是由于该区域的裂缝宽度在塔筒破坏过程中始终较小(图14)，这也体现在该区域钢筋应力在混凝土开裂后增长不大(图15(e))。而由图15(f)可知，λ=1.529和1.696 时，hs/Hs=0.70 高度θ=0°内侧和θ=50°外侧混凝土环向应力先后达到抗拉强度并进入卸载阶段，这两处混凝土产生子午向裂缝，且钢筋的拉应力迅速增加。尽管这两处环向钢筋均未屈服，但λ=1.98 时迎风子午线上hs/Hs=0.57 处内侧环向钢筋屈服，并迅速扩展至hs/Hs=0.53～0.61 范围。这也表明，尽管结构破坏源于子午向受拉破坏，但环向受力亦需关注。

5 结论

以一座大型RC 冷却塔为例，通过数值模拟研究了冷却塔在两种荷载工况下的非线性静风响应，从荷载位移曲线、变形模式、内力、应变和应力分布对其破坏过程进行了系统阐述，并从中探究了其破坏机理。主要结论如下：

(1)通过对Port Gibson 冷却塔的计算和与既有结果对比，表明ABAQUS 的分层壳单元配合弥散开裂模型可以有效进行冷却塔极限静风荷载分析。

(2)在自重+风荷载作用下，在λ=1.384 时，迎风点hs/Hs=0.37 位置在子午向轴拉力作用下首先开裂并产生贯穿厚度的环向裂缝，之后裂缝随风荷载增加而不断形成并沿环向扩展至整个迎风区，沿子午向几乎扩展至整个高度，开裂位置的钢筋应力迅速增加；λ=1.529 时，迎风子午线hs/Hs=0.62高度内侧和θ=75°、hs/Hs=0.95 高度外侧混凝土因环向弯矩而开裂；λ=1.876 时，hs/Hs=0.70 处子午向钢筋首先受拉屈服，之后钢筋屈服区域沿环向逐渐扩展至θ=16°，沿子午向则扩展至hs/Hs=0.52～0.80；λ=1.98 时迎风子午线上hs/Hs=0.57 处内侧环向钢筋屈服，并迅速扩展至hs/Hs=0.53～0.61 范围。λ=2.007 时，因混凝土持续开裂和钢筋屈服，冷却塔达到抗拉极限状态而破坏。

(3)在自重+风荷载作用下，塔筒开裂之前各响应均表现为线性特征，塔筒开裂之后，塔筒局部刚度下降，荷载-位移曲线也由线性进入非线性并且迎风区和侧风区位移迅速增加，同时迎风区和侧风区内力也表现出明显的重分布特征。结构破坏时，迎风区子午向受拉导致大面积的受拉贯通环向裂缝和钢筋屈服，而环向仅有塔筒中上部受弯开裂且裂缝宽度有限，钢筋屈服范围也较小，因此可认为迎风子午向受拉导致的塔筒开裂和钢筋屈服是结构破坏的首要原因，但环向受力亦需关注。

(4)在自重+冬温荷载+风荷载作用下，由于冬温荷载产生的双向正弯矩效应增加了塔筒的弯曲损伤，导致塔筒侧风区在λ=1.0 时提前产生子午向裂缝。但由于此弯矩效应固定不变且远小于风荷载作用下的子午向受拉效应，故结构的破坏依然由风荷载控制，极限荷载系数为λ=1.842。

(5)两种荷载工况下结构均在λ=1.50 时明显进入非线性段，且之后塔筒变形急速增加，此荷载系数这也与结构设计时的风荷载分项系数1.4 相当，说明在设计风荷载下结构的冗余安全度较为有限。