超大型冷卻塔風致倒塌全過程數(shù)值仿真與受力性能分析

吳鴻鑫，柯世堂, ，王飛天，葛耀君

(1.南京航空航天大學土木工程系，江蘇，南京，210016；2.同濟大學土木工程防災國家重點實驗室，上海 200092)

大型冷卻塔屬于鋼筋混凝土旋轉(zhuǎn)薄殼結(jié)構(gòu)，風荷載是其結(jié)構(gòu)設計的控制荷載。歷史上國內(nèi)外大型冷卻塔的風毀事件多次發(fā)生，事后調(diào)查[1－4]表明風致倒塌均是由局部風荷載過大使得塔筒局部受拉損傷，進而引發(fā)整個塔筒結(jié)構(gòu)強度破壞和屈曲失穩(wěn)倒塌。事實上，此類由局部損壞引發(fā)的整體風毀破壞均屬于典型的結(jié)構(gòu)連續(xù)倒塌破壞形式。結(jié)構(gòu)的連續(xù)倒塌[5]是指由于偶然荷載造成結(jié)構(gòu)的局部破壞，并導致結(jié)構(gòu)產(chǎn)生非穩(wěn)定的連續(xù)性破壞。國外相關(guān)設計指南[6－7]及國內(nèi)《建筑結(jié)構(gòu)抗倒塌設計規(guī)范》[8]提供了具體的抗連續(xù)倒塌設計方法，其主要目標是針對RC框架結(jié)構(gòu)進行抗震設計，均缺少此類特種高聳薄殼結(jié)構(gòu)的風致倒塌相關(guān)條款。

已有研究[9－12]表明大型冷卻塔風致倒塌的主要評價指標有屈曲失穩(wěn)、強度破壞、疲勞損傷等，目前高度突破規(guī)范限值的200 m級超大型冷卻塔風致倒塌關(guān)注的重要內(nèi)容是結(jié)構(gòu)穩(wěn)定性能?，F(xiàn)有針對大型冷卻塔的抗風穩(wěn)定性研究，主要集中在結(jié)構(gòu)靜風失穩(wěn)[13]、整體與局部屈曲失穩(wěn)[14－15]、施工全過程穩(wěn)定性[16－17]和幾何缺陷所致失穩(wěn)[18]等方面。相關(guān)成果很好地指導了此類超大型冷卻塔抗風穩(wěn)定性設計，并對大型冷卻塔倒塌過程認識具有一定的借鑒意義。然而，均忽略了由局部失穩(wěn)引發(fā)的整體連續(xù)倒塌破壞的后續(xù)現(xiàn)象，難以揭示此類超大型冷卻塔風致倒塌破壞過程及內(nèi)在受力機制。

鑒于此，以山西潞安電廠世界最高220 m超大型冷卻塔結(jié)構(gòu)為例，基于CFD與LS-DYNA技術(shù)提出了大型冷卻塔風致倒塌全過程數(shù)值仿真模擬方法，數(shù)值再現(xiàn)了超大型冷卻塔風致倒塌全過程。對比研究了塔筒應力分布變化規(guī)律與倒塌全過程的扭曲變形姿態(tài)等特征，提煉出冷卻塔風致倒塌的受力特點與作用機制，并討論了單元失效參數(shù)的影響。研究結(jié)論為此類超大型冷卻塔結(jié)構(gòu)抗倒塌設計提供了科學依據(jù)和參考價值。

1 工程概況

本文以山西潞安電廠已建成世界最高220 m超大型冷卻塔為工程背景，主要結(jié)構(gòu)參數(shù)見圖1。該鋼筋混凝土雙曲線自然通風間接空冷塔高為220 m，塔頂直徑為128.1 m，喉部離地高度為165 m，喉部直徑為123.0 m，通風口離地高度為30.75 m，塔底直徑為 185.0 m，冷卻塔風筒塔體采用指數(shù)變厚，最小壁厚為0.38 m，最大壁厚為1.85 m。由64榀X型支柱與環(huán)板基礎連接，混凝土X型柱采用矩形截面(1.7 m×1.0 m)。環(huán)板基礎為現(xiàn)澆鋼筋混凝土結(jié)構(gòu)，寬為10.5 m，高為2.2 m。

圖1 超大型冷卻塔主要結(jié)構(gòu)尺寸示意圖Fig.1 Dimensions of super large cooling tower

2 數(shù)值仿真方法與有效性驗證

2.1 CFD數(shù)值模擬

為了保證冷卻塔尾流能充分發(fā)展，計算域尺寸設置為：x=16D，y=20D，z=4H(x為橫風向，y為順風向，z為高度方向)，其中，D為塔底直徑，H為塔高，模型的最大堵塞率不超過 5%，滿足規(guī)范要求。為兼顧計算效率與精度，并且考慮冷卻塔形狀復雜，接近結(jié)構(gòu)的內(nèi)部區(qū)域流場較復雜，整個計算域劃分采用混合網(wǎng)格的離散形式，將其分割為外部區(qū)域與內(nèi)部加密區(qū)域，其中外部區(qū)域采用高質(zhì)量的六面體結(jié)構(gòu)化網(wǎng)格進行劃分，內(nèi)部加密區(qū)域采用四面體非結(jié)構(gòu)化網(wǎng)格進行劃分。內(nèi)部加密區(qū)域中，結(jié)構(gòu)模型的最小網(wǎng)格尺寸為1 m，除結(jié)構(gòu)之外的計算域網(wǎng)格尺寸為 15 m，整體計算模型總網(wǎng)格數(shù)量為1200萬，網(wǎng)格的最大歪斜率為0.7(不大于0.9)，模型壁面y+值為35.2，可保證底層網(wǎng)格位于對數(shù)律層中，網(wǎng)格的質(zhì)量與數(shù)量滿足計算要求。計算域邊界條件采用速度入口、壓力出口、自由滑移壁面(側(cè)面及頂面)與無滑移壁面(地面及塔壁)。圖2給出CFD網(wǎng)格劃分、測點布置及邊界條件示意圖。

圖2 CFD網(wǎng)格劃分及邊界條件示意圖Fig.2 CFD meshing and boundary conditions

所在地為B類地貌，大氣邊界層采用對應分布形式，如式(1)和式(2)所示：

式中：Z為距地面高度；α=0.15為B類地貌的地面粗糙度指數(shù)；U0=23.7 m/s為該地區(qū)Z0=10 m高度處50年10 min最大平均風速；I10=0.14為該地區(qū)Z0=10 m高度處名義湍流強度。

通過用戶自定義函數(shù)(UDF)定義上述脈動風場速度入口，在模型前10 m處設置10個測點，以此測量UDF所定義速度入口是否達到穩(wěn)定，圖3給出平均風速、湍流度剖面模擬結(jié)果與理論值的對比曲線，表明入流風場與理論值吻合良好，風場模擬標準滿足工程要求。

圖3 速度及湍流度剖面示意圖Fig.3 Velocity and turbulence profile

數(shù)值計算采用3D單精度、Pressuredbased求解器，空氣流場采用不可壓縮空氣流場，湍流模型采用RNG的k-ε湍流模型、壓力速度耦合方程組求解采用SIMPLEC格式。梯度離散選用Least Squares Cell Based格式，壓力離散采用Standard格式，動力離散采用二階迎風格式，控制方程的計算殘差設置為 1×10-6。

圖4給出了冷卻塔風場三維速度流線圖及湍動能分布，由圖4可知：來流流經(jīng)塔筒在迎風面產(chǎn)生分流，進而沿塔筒兩側(cè)外壁加速流動至背風面形成尺寸不同的渦旋，部分氣流透過百葉窗進入塔筒內(nèi)部，氣流產(chǎn)生回流并在喉部形成較大尺度的渦旋；塔筒外部湍動能強度顯著大于內(nèi)部，峰值主要位于塔筒迎風區(qū)頂部、背風區(qū)喉部、百葉窗迎風區(qū)內(nèi)部及背風區(qū)外部。

圖5給出了冷卻塔外表面壓力系數(shù)云圖，由圖5可知：冷卻塔每層的最大風壓系數(shù)出現(xiàn)在0°子午向上，其中塔頂?shù)妮^大，塔底較??；而冷卻塔每層的風壓系數(shù)最小值出現(xiàn)在子午向72°與288°上，最大負壓由塔底向上逐增到喉部位置達到最大負壓，越過喉部后逐減至塔頂。

圖4 冷卻塔風場模擬結(jié)果示意圖Fig.4 Schematic diagram of cooling tower wind field simulation results

圖5 冷卻塔風壓系數(shù)云圖Fig.5 Cloud diagram of cooling tower wind pressure coefficient

圖6給出典型喉部薄弱層(第8層~第10層)平均風壓數(shù)值分析結(jié)果，作為后續(xù)風荷載輸入?yún)?shù)，并將其與規(guī)范曲線[19]、風洞試驗結(jié)果[20]及實測西熱曲線[21]進行對比。分析可知，迎風和背風區(qū)域風壓系數(shù)吻合較好，側(cè)風區(qū)略大于規(guī)范，但與西熱曲線較為一致。

圖6 風壓分布系數(shù)曲線模擬結(jié)果示意圖Fig.6 Schematic diagram of simulation results of wind pressure distribution coefficient curve

2.2 LS-DYNA顯示動力分析

基于宏模型建模方法建立該超大型冷卻塔三維有限元模型。其中塔筒采用殼單元Shell163，塔筒沿高度方向劃分為128層，沿環(huán)向劃分為256個單元；X型支柱采用梁單元beam161，上部塔體殼單元與下部柱梁單元的連接，采用殼單元與梁單元節(jié)點剛域耦合，以達到位移協(xié)調(diào)的目的；地面采用剛性墻平面，圖7給出冷卻塔有限元模型示意圖。

圖7 冷卻塔有限元模型示意圖Fig.7 Schematic diagram of element model of cooling tower finite

超大型冷卻塔筒壁較薄，其縱向和環(huán)向配筋率較大，鋼筋的強大拉接力使得塔壁形成了一個整體受力良好的聯(lián)合體，具有明顯的塑性隨動特點。相關(guān)研究[22－23]發(fā)現(xiàn)可將混凝土和鋼筋進行綜合考慮，以單元失效準則來模擬倒塌過程中的構(gòu)件破壞行為。本文采用的具有失效機制的塑性隨動模型(PLAW)整體考慮塔筒的材料特性，表1給出材料模型各參數(shù)。

表1 材料模型參數(shù)列表Table 1 Material model parameter list

采用LS-DYNA提供的接觸計算模型單面接觸(*Contact_Automatic_Single_Surface)，材料的摩擦系數(shù)統(tǒng)一設定為0.25。在LS-DYNA中采用單點積分算法，提高了計算效率，為了削弱不利的零能模式(沙漏模式)，通過粘性阻尼法與小剛度法控制沙漏模式消耗的能量。

實際上，大型冷卻塔在風荷載作用下，其結(jié)構(gòu)的破壞形式是類似屈曲模態(tài)下的材料破壞，可以視風荷載為一種擬靜力荷載。為了減低加載過程中的動力效應，防止出現(xiàn)類似沖擊荷載的效應，將平均風壓在20 s內(nèi)線性遞增形成時程分析輸入所需的風荷載，再以點荷載加載到塔筒相應的加載點。

為了避免施加重力荷載引起的塔筒附加振動，荷載分為兩個階段施加：首先，施加重力荷載，同時施加給模型一個較大的阻尼，計算直至冷卻塔在重力作用下達到穩(wěn)定為止；然后，刪除大阻尼設置，利用“完全重啟動”技術(shù)重新定義冷卻塔等效綜合阻尼比為 2%[24]，并調(diào)用第一階段結(jié)束時冷卻塔在重力荷載作用下的受力狀態(tài)，之后施加風荷載，進行迭代計算。

2.3 動力特性分析及有效性驗證

基于分塊Lanczos方法進行風致冷卻塔倒塌數(shù)值模型的動力特性分析，圖8給出LS-DYNA下冷卻塔前200階自振頻率隨振型階數(shù)變化曲線圖。由圖8分析：冷卻塔的奇數(shù)階和偶數(shù)階頻率相同，具有對稱性；自振頻率分布非常密集，基頻僅為0.68 Hz，前50階振型頻率主要集中分布在0.6 Hz~ 1.9 Hz。

圖8 冷卻塔自振頻率分布示意圖Fig.8 Schematic diagram of cooling tower natural frequency distribution

表2給出冷卻塔結(jié)構(gòu)典型模態(tài)振型圖，由表2分析發(fā)現(xiàn)：塔筒振型在環(huán)向和子午向差異較大，隨著階數(shù)的增大，環(huán)向和子午向諧波數(shù)顯著增加；而在低階振型中，子午向表現(xiàn)為1~3個諧波，環(huán)向表現(xiàn)為 6~12個諧波；其中，在前幾階振型中，多以塔頂振動為主，以塔底振動為主振型中，環(huán)向諧波數(shù)較多且塔頂振幅較小。

為了驗證基于CFD與LS-DYNA技術(shù)的大型冷卻塔風致倒塌全過程數(shù)值仿真方法的有效性，以文獻[25]開展的大型冷卻塔整體穩(wěn)定性風洞試驗為驗證對象進行數(shù)值仿真對比。試驗冷卻塔為塔高212 m的間接空冷冷卻塔，出口半徑為111.38 m，喉部高度為159.11 m，喉部直徑為105.4 m，進風口高度為32.2 m，進風口直徑為142.97 m，最小壁厚為 0.36 m，關(guān)于風洞試驗的具體介紹，詳見文獻[25]。

表2 冷卻塔典型模態(tài)固有頻率及振型列表Table 2 Cooling tower typical modal natural frequency and vibration mode list

圖9給出了該大型冷卻塔風致失穩(wěn)破壞形態(tài)的對比示意圖。由圖9對比發(fā)現(xiàn)兩者具有相同的倒塌趨勢，均是從喉部變橢圓到出現(xiàn)褶痕，最后，由剛性環(huán)撕扯破壞，兩者僅具體畸變形狀略有差異，本文數(shù)值仿真獲得的破壞形態(tài)和作用過程均能有效反映風洞試驗整體倒塌。

基于流動相似理論，由試驗[25]屈曲風速換算得實際倒塌風速為88.48 m/s，與本文數(shù)值模擬結(jié)果相差 9.6%；圖10給出數(shù)值仿真與文獻[25]的穩(wěn)定性分析結(jié)果對比示意圖，數(shù)值仿真的無量綱數(shù)log(qcr/E)與非線性穩(wěn)定性結(jié)果相差僅為 0.15%，與本文顯式動力分析的基本假設吻合，說明本文數(shù)值仿真方法具有一定有效性。

圖10 冷卻塔穩(wěn)定性分析結(jié)果對比示意圖Fig.10 Comparison of cooling tower stability analysis results

3 倒塌全過程分析

3.1 倒塌形態(tài)與塑性應變

基于增量動力分析(IDA)方法對結(jié)構(gòu)進行非線性分析，以25 m/s為起始基本風速進行逐級加載，加載風速步長為 5 m/s。當風速增大至混凝土受拉破壞(C40，ftk≥1.71 MPa)時，局部區(qū)域混凝土開裂而鋼筋受拉，隨著風速進一步增大，塔筒受壓區(qū)接近極限受力狀態(tài)(C40，fck≥19.1 MPa)，冷卻塔風致響應顯著增大。因此，本文以冷卻塔最大位移變化率急劇增大作為大型冷卻塔倒塌的判斷依據(jù)，圖11給出冷卻塔塔頂與喉部最大位移隨風速變化示意圖。

從圖11中可以發(fā)現(xiàn)：該位移-風速曲線存在一個明顯臨界風速(80 m/s)，在達到該臨界風速之前，位移與風速基本呈線性關(guān)系，冷卻塔塔筒發(fā)生小變形，處于彈性階段；超過該臨界風速后，位移呈現(xiàn)發(fā)散突增，冷卻塔塔筒發(fā)生大變形，進入塑性連續(xù)倒塌階段。研究表明：數(shù)值仿真結(jié)果的臨界風壓為18.62 kPa，規(guī)范[19]考慮整體穩(wěn)定性時的臨界風壓為17.67 kPa，兩者相差5.10%，規(guī)范相對于數(shù)值模擬結(jié)果較為保守；該大型冷卻塔在風荷載與自重荷載耦合作用下，臨界倒塌風速為80 m/s。

圖11 冷卻塔喉部最大位移隨風速變化示意圖Fig.11 Schematic diagram of maximum displacement of cooling tower throat as a function of wind speed

圖12和圖13分別給出大型冷卻塔在臨界風速80 m/s的數(shù)值仿真倒塌形態(tài)示意圖與有效塑性應變示意云圖。由圖分析發(fā)現(xiàn)：在風致冷卻塔倒塌全過程中，首先是冷卻塔喉部凹陷，俯視冷卻塔喉部變成橢圓型，此時，冷卻塔有效塑性應變在迎風面由喉部延伸到塔底；然后在喉部繼續(xù)凹陷發(fā)展中，而頂部剛性環(huán)凹陷較緩慢，牽扯住喉部，形成兩道從冷卻塔喉部迎風點斜向上的褶痕，同時，有效塑性應變也隨褶橫斜向上發(fā)展；接著隨著冷卻塔頂部剛性環(huán)繼續(xù)凹陷，達到失效應變發(fā)生斷裂崩塌，剛性環(huán)從迎風點向內(nèi)折陷，同時喉部褶痕處剪切撕扯破壞，有效塑性應變從喉部沿子午向向上發(fā)展至剛性環(huán)；最后，頂部剛性環(huán)碎塊沖擊背面冷卻塔，由重力與風荷載共同作用，由于殼體牽連破壞，發(fā)生連續(xù)倒塌。

圖12 冷卻塔倒塌形態(tài)示意圖Fig.12 Schematic diagram of cooling tower collapse

圖13 冷卻塔倒塌過程有效塑性應變示意圖Fig.13 Schematic diagram of effective plastic strain during cooling tower collapse

3.2 倒塌初始(t=0 s)應力響應

圖14給出了大型冷卻塔在臨界風速下倒塌初始(t=0 s)的應力分布圖，圖中正值表示拉應力。由圖14發(fā)現(xiàn)：冷卻塔達到承載極限開始倒塌時，子午向三維應力以 0°子午軸左右對稱，子午向應力從塔頂往下遞增至喉部附近出現(xiàn)最大值，在喉部至塔筒底部的子午向應力相近；大型冷卻塔在臨界風速下倒塌初始(t=0 s)時，喉部環(huán)向應力以0°子午線為軸左右對稱，最大 Mises應力并非在正中間而是出現(xiàn)在子午向15°與-15°上，從而形成喉部左右的褶皺區(qū)。

圖14 塔筒倒塌初始(t=0 s)應力分布圖Fig.14 Initial (t=0 s)stress distribution of tower collapse

3.3 倒塌受力特點

大型冷卻塔的塔筒是典型的雙曲線薄殼系統(tǒng)，在倒塌過程中經(jīng)歷了小變形和大變形兩個階段。圖15給出冷卻塔倒塌受力特點與相應的力學模型示意圖。大型冷卻塔失穩(wěn)前變形很小且處于彈性階段，其幾何非線性和材料非線性的影響可以忽略不計，在小變形下，大型冷卻塔主要通過薄殼的受彎承載力抵抗結(jié)構(gòu)倒塌。冷卻塔塔筒的倒塌是突然發(fā)生的，首先，在迎風面的喉部附近發(fā)生局部破壞，隨后碎裂范圍沿高度上下及環(huán)向不斷擴展，在大變形下，塔筒喉部迎風點產(chǎn)生塑性鉸，進入塑性變形而發(fā)揮懸鏈線機制。

圖15 冷卻塔倒塌受力特點與力學模型示意圖Fig.15 Mechanical characteristics and mechanical model of cooling tower collapse

3.4 單元失效參數(shù)

倒塌過程中的構(gòu)件破壞會導致結(jié)構(gòu)內(nèi)力重分布，剩余結(jié)構(gòu)構(gòu)件發(fā)生卸載或過載，引起結(jié)構(gòu)破壞的連鎖反應。由于本數(shù)值仿真方法采用具有失效準則的雙線性隨動強化模型(PLAW)，其中雙線性隨動強化模型的屈服應力反映殼體開裂時的應力，而失效應變反映殼體斷開時的應變。為了考慮不同單元失效參數(shù)對結(jié)構(gòu)倒塌的影響，在彈性模量、切線模量、硬化參數(shù)、應變率參數(shù)保持不變的前提下，采取屈服應力和失效應變不同取值的三種組合，基于增量動力分析(IDA)方法，對結(jié)構(gòu)的這三種工況進行非線性分析提取倒塌風速，表3給出不同單元失效參數(shù)的倒塌風速列表。

表3 不同單元失效參數(shù)的倒塌風速Table 3 Collapsed wind speeds with different failure parameters

對比發(fā)現(xiàn)，雙線性隨動強化模型的屈服應力參數(shù)對倒塌風速影響較大，而失效應變參數(shù)對倒塌風速影響較小，對于典型的雙曲線大型冷卻塔，倒塌主要是因為整體失穩(wěn)，局部殼體的開裂破壞，剩余構(gòu)件的內(nèi)力重分布不能達到穩(wěn)定，進而連續(xù)倒塌。

4 結(jié)論

本文提出了基于CFD與LS-DYNA技術(shù)的大型冷卻塔風致倒塌全過程數(shù)值仿真模擬方法，數(shù)值再現(xiàn)了世界第一高度220 m超大型冷卻塔風致倒塌全過程。在此基礎上，對比研究了此類超大型冷卻塔風致倒塌的受力特征與作用機理，為基于結(jié)構(gòu)強健性的大型冷卻塔抗風設計提供理論依據(jù)。

研究表明：本文數(shù)值方法可以有效模擬超大型冷卻塔風致倒塌全過程；大型冷卻塔風致倒塌過程始于塔筒喉部大變形，并在兩側(cè) 30°范圍內(nèi)呈現(xiàn)褶皺現(xiàn)象，最終因變形不協(xié)調(diào)而相互牽扯垮塌；其倒塌受力機制可分為彎拱機制與懸絞線機制兩種，大型冷卻塔失穩(wěn)前變形很小且處于彈性階段，主要以圓拱形薄殼的受彎承載力抵抗結(jié)構(gòu)倒塌，大型冷卻塔失穩(wěn)時變形突增進入塑性階段，在塔頂迎風點形成塑性鉸抵抗結(jié)構(gòu)連續(xù)倒塌；材料模型的單元失效參數(shù)對超大型冷卻塔風致倒塌的影響不可忽略，雙線性隨動強化模型的屈服應力參數(shù)對倒塌風速影響較大，而失效應變參數(shù)對倒塌風速影響較小。