基于分離渦方法的臺北101大廈流固耦合風致響應分析

盧春玲，劉宇杰，陳錦焜，李秋勝

(1.桂林理工大學 土木與建筑工程學院，桂林 541004；2.廣西巖土力學與工程重點實驗室，桂林 541004；3.桂林理工大學 廣西有色金屬隱伏礦床勘查及材料開發(fā)協(xié)同創(chuàng)新中心，桂林 541004；4.中南大學 土木工程學院，長沙 410082)

流固耦合指的是可變形體在風流場中，由于風載荷作用產(chǎn)生變形，而結構變形又反過來作用于風流場，使風流場產(chǎn)生重分布，導致流體對固體作用的荷載大小和分布發(fā)生變化。它研究的是流體與固體之間的相互作用。對于阻尼低，自振頻率低的超高層建筑物而言，這種流固耦合效應更加明顯。近年來，隨著計算技術的不斷提高，流固耦合相關理論與應用也得到了極大的發(fā)展，國內外學術界已經(jīng)意識到流固耦合效應對建筑結構的抗風穩(wěn)定性有不可忽視的影響[1-5]。在對建筑結構進行流固耦合相關研究中，對進行流場分析時采用分離渦(Detached Eddy Simulation，DES)算法是一種較為新型的數(shù)值模擬方法，相比于目前較常用的大渦模擬方法(Large Eddy Simulation, LES)，分離渦模型能在保持精度的前提下有效降低計算量。目前，采用分離渦算法對高層建筑進行流固耦合的研究在國外已有進行[6-7]。本文應用分離渦模擬并結合一種新的可滿足大氣邊界層中風場特性的湍流脈動速度生成方法——離散再合成的隨機湍流生成法(DSRFG)[8]模擬非穩(wěn)態(tài)邊界層湍流風場。以復雜超高層建筑臺北101大廈為對象，對該超高層建筑進行建筑表面風壓分布與周圍流場分布的數(shù)值模擬，模擬得到臺北101大廈剛體模型周圍的風流場及作用于其上的風荷載時程數(shù)據(jù)，然后將模型表面風荷載時程施加在超高層建筑有限元模型上，獲得未考慮流固耦合效應的結構風致響應分析結果。同時基于WORKBENCH平臺構建有限元模型，采用考慮雙向流固耦合的方法對超高層建筑模型進行模擬計算，將計算結果與現(xiàn)場實測以及風洞試驗的相應數(shù)據(jù)進行對比，以對該數(shù)值模擬方法的實用性和準確性進行檢驗。同時將考慮與未考慮流固耦合的計算結果進行對比，分析流固耦合效應對101大廈風致響應及結構周圍流場分布特性的影響。

1 風洞試驗

臺北101大廈高頻天平測力風洞試驗是由Rowan Williams Davies&Irwin Inc.(RWDI)于安大略圭爾夫市的風洞實驗室進行的。臺北101大廈模型采用與建筑物的幾何外形相似的剛性模型，模型縮尺比為1∶500，如圖1所示。在RWDI 4.9 m×2.4 m的邊界層風洞中，考慮主體建筑600 m半徑范圍內所有主要建筑物，設定風速采用經(jīng)過實際風速折減后的試驗風速，采用高頻天平對臺北101大廈的剛性縮尺模型進行測力試驗。風洞試驗測得湍流強度剖面，見圖2。

圖1 臺北101大廈風洞試驗模型

圖2 風洞試驗湍流強度剖面

2 大廈振動監(jiān)控系統(tǒng)

臺北101大廈動力特性，如自振頻率、模態(tài)及阻尼比等數(shù)據(jù)，均來源于在臺北101大廈上安裝的由傳感器和數(shù)據(jù)采集處理系統(tǒng)組成的振動監(jiān)測系統(tǒng)。該傳感器系統(tǒng)包括30個加速傳感器，分別安裝在大廈的-5 th(地下室)，1 st，36 th，60 th，89 th和101 st層，如圖3(a)所示。加速度傳感器在大廈平面上的布置如圖3(b)所示。用于對比分析的現(xiàn)場實測數(shù)據(jù)是由采樣頻率為200 Hz的數(shù)據(jù)采集系統(tǒng)從2005年8月到2008年5月期間，三次臺風作用和一次地震作用下采集得到。

(a)傳感器立面布置圖

3 數(shù)值模擬方法

本文對臺北101大樓進行考慮與未考慮流固耦合數(shù)值風洞模擬。對未考慮流固耦合的數(shù)值風洞，首先建立剛性模型與ANSYS有限元模型，運用fluent軟件計算并提取剛性模型表面風壓時程，施加在ANSYS有限元模型中獲取其風致響應。對考慮流固耦合的數(shù)值風洞，則基于結構動力學特性相似原理和workbench平臺，建立外形與剛性模型相同，并具有與ANSYS有限元模型相同質量、體積、振型與阻尼比的workbench有限元模型，依靠流固耦合模塊直接計算得到其風致響應。

3.1 計算模型及網(wǎng)格劃分

數(shù)值風洞計算域見圖4(a)，在順風向(X軸方向)的長度為32Db(-9

(1)

經(jīng)計算流域與模型滿足阻塞率小于3%的要求。

本文對計算域中的流體和固體區(qū)域網(wǎng)格分別進行劃分。對于流體區(qū)域，不考慮流固耦合的剛體模型流場采用混合網(wǎng)格劃分，網(wǎng)格總數(shù)量320萬左右，網(wǎng)格劃分方式見圖4(b)。考慮流固耦合的有限元模型流場采用非結構化網(wǎng)格進行劃分，網(wǎng)格數(shù)量313萬左右。

對于建筑實體部分，未考慮流固耦合的臺北101大樓有限元模型通過有限元軟件ANSYS建立，分別使用梁單元BEAM188、殼單元SHELL181、桿單元LINK8以及質量單元MASS21構建全尺寸主體結構，并忽略建筑圍護結構，大廈的三維有限元模型共計20 532個梁單元，24 048個板單元, 以及3 496個桿單元。ANSYS有限元模型如圖4(d)所示?？紤]流固耦合的建筑有限元模型基于Workbench平臺建立，為便于流體與固體之間耦合面的建立，并在保留建筑模態(tài)特性的同時提高計算效率，將建筑模型構建為均質等效模型。建筑單元類型為solid186，單元總數(shù)約32萬。建筑物整體質量為2.391×108kg，模型體積為1.039×106m3，故材料密度取為230.12 kg/m3；彈性模量取為1.3×109Pa，泊松比取為0.2。

(a)數(shù)值風洞計算域整體尺寸

3.2 邊界條件

對于流體計算域，剛性模型與彈性模型均采用分離渦模擬的數(shù)值模擬方法并采用相同的邊界條件。采用作者提出的一種新的湍流入口生成方法-DSRFG方法。速度入口處平均風速沿高度的變化服從指數(shù)率，計算公式為

(2)

式中:v、v0為任意高度與標準參考高度處的平均風速(m/s)；z、z0為任意高度與標準參考高度(m)；α為地面粗糙度指數(shù)，取值0.15。文中z0=10 m，10 m高度處50年重現(xiàn)期平均風速v0為43.27 m/s。

基于Shiau[9]在臺灣基隆港的觀測塔測(距離臺北北部約20 km)得的26 m高度處的強風(臺風)數(shù)據(jù)的結果確定臺北101大廈入口處的湍流強度剖面及湍流積分尺度的分布。湍流強度與風洞試驗結果進行對比如圖5所示，從圖5可看出數(shù)值風洞與風洞試驗的湍流強度剖面趨勢相同，但在數(shù)值上有差異。采用DSRFG方法生成的入口處的瞬時風速分布圖，如圖6所示。

圖5 入口湍流強度及積分尺度剖面

圖6 DSRFG方法產(chǎn)生的入口瞬時風速場

通過DSRFG方法生成的入口處脈動風速時程的功率譜與目標譜卡曼譜的比較如圖7所示。由圖7可知，大渦模擬所生成的脈動風速時程功率譜與目標譜一致，故其功率譜特性符合大氣邊界層湍流。其他邊界條件，見表1。

表1 邊界條件設置

圖7 DSRFG方法生成的脈動風速時程功率譜與目標譜比較

在結構動力響應計算中，設置大樓在基礎處的約束為3向固定約束(Fixed Support)。在workbench有限元模型中，設置建筑物除去底面的其他面為流固耦合面(Fluid Solid Interface)。

3.3 湍流模型

對未考慮流固耦合的剛性模型與考慮流固耦合的模型流體部分均采用分離渦模型作為湍流模型。分離渦模型是Spalart于1997年提出的一種三維非穩(wěn)態(tài)數(shù)值模型，其原理是在模型近壁面采用RANS模型進行模擬，而在遠端采用LES方法模擬，在保證相對高的計算精度的同時，有效降低了計算量。為獲得較高的分離渦模擬精度，本文選用延遲分離渦模擬[10-11](Delayed Detached Eddy Simulation，DDES)進行分離渦的計算。該方法構造了一個轉換函數(shù)fd，表達式為：

fd=1-tanh[(8rd)3]

(3)

其中，

(4)

式中：νt為運動學渦黏性；ν為分子黏度；Uij為速度梯度；k為卡門常數(shù)；d為距壁面距離。

利用fd對經(jīng)典分離渦模型進行修正，能夠重新定義分離渦模型的長度尺度，延遲了RANS模型的作用范圍，有效解決由于邊界層網(wǎng)格加密產(chǎn)生的模型應力損耗問題和網(wǎng)格誘導分離問題，提高模擬的準確度。為使分離渦模型的方程封閉，需要設定亞格子模型和RANS模型。分離渦模型選擇湍動能運輸亞格子模型(Kinetic-Engergy Transport，KET)，由Kim等[12]提出，這種亞格子模型對不均勻湍流的復雜流場有著較好的模擬效果。RANS模型則設定為Realizablek-ε兩方程模型。對于考慮雙向流固耦合的模型，采用分離渦模型進行流場的數(shù)值模擬，亞格子模型設置與剛性模型相同。

3.4 求解參數(shù)

對流體的求解基于流體計算軟件Fluent，采用分離渦法對臺北101大廈進行模擬計算。對風場的計算采用3D單精度、分離式求解器，空氣模型選用不可壓縮的常密度空氣模型，基于壓力求解器(Pressure-Based)，時間類型選擇瞬態(tài)(Transient)，求解算法采用SIMPLEC算法，離散格式選用二階迎風格式。時間步長設置為0.05 s，時間步數(shù)設置為12 000。為獲得較穩(wěn)定的入口湍流，將定常Standardk-ε模型計算結果通過瞬態(tài)化處理作為分離渦模型計算的初始場。

在workbench有限元模型中，不僅要設置流體的求解參數(shù)，還要對耦合求解模塊進行設置。在Workbench平臺耦合計算模塊中設置求解類型為瞬態(tài)(Transient)，迭代結束時間設置為600 s，時間步長為0.3 s，每步最大迭代子步采用默認的5步，將流體部分的建筑物表面與固體部分的流固耦合面進行數(shù)據(jù)傳輸，并設置計算順序，即流體計算-固體計算依次迭代，最終計算完成后在分析模塊中查看計算結果。

4 數(shù)值模擬結果分析

4.1 流場分析

本文對考慮流固耦合與未考慮流固耦合的數(shù)值模擬模型進行流場分析，分離渦模型在考慮流固耦合與未考慮流固耦合情況下風場中平均流速分布及模型y=0平面渦量分布見圖8、9。

由圖8～9可以看出，對于考慮流固耦合的氣彈模型與未考慮流固耦合的剛性模型，來流在建筑周圍表現(xiàn)出相似的典型鈍體擾流的流場特點。上部來流越過建筑物頂端后向下分離，在建筑頂端背風處形成駐渦區(qū)。來流通過建筑物頂端與兩側繞流，不同頻率的風速脈動與分離區(qū)域中的渦旋結構相互作用，在建筑側面與背風面形成密集的渦旋脫落，這也是高層建筑橫風向脈動風力產(chǎn)生的主要原因。

(a)剛性模型y=0平面

由圖9(a)、(c)可以看出，相比于有限元模型，剛性模型周圍渦量分布更為密集，尾流區(qū)域也更廣泛，這使得不考慮流固耦合的模型順風向加速度與位移響應均大于考慮流固耦合的模型。由圖9(b)、(d)可以看出，雙流固耦合作用下建筑背面尾流區(qū)域較剛性模型存在明顯的旋渦脫落現(xiàn)象，旋渦更多且渦量更大。但在有限元模型中，建筑側面的分離泡在形成初期即與建筑分離，形成脫落渦旋時也與建筑距離較遠，而對于剛性模型，兩側的分離泡與建筑有較大的附著面積，建筑物兩側渦量明顯大于氣彈模型，這使得剛性模型中建筑受到更大的橫向力影響，在風致響應中則表現(xiàn)為未考慮流固耦合的模型擁有更大的橫風向位移與加速度。在有限元模型中，漩渦脫落時其渦量較大，可能會對其下游的毗鄰建筑造成影響。

(a)剛性模型y=0

4.2 模態(tài)分析

對臺北101大廈ANSYS有限元模型與workbench有限元模型進行模態(tài)分析，將前五階自振頻率與實測數(shù)據(jù)進行對比，將結果列于表2。從表2可看出，ANSYS模型模態(tài)計算值與實測值差距都在15%以內，而workbench模型模態(tài)計算值與實測差距在5%以內，具有較高的精確度與可信度。

表2 前五階振型自振頻率計算值與實測值對比

4.3 等效靜力風荷載

本文采用規(guī)范計算以及時域分析兩種方法，對臺北101大廈等效靜力風荷載進行計算。

規(guī)范中順風向風荷載標準值ωk按下式計算：

ωk=βzμsμzω0

(5)

式中：βz為高度z處的風振系數(shù)；μz為風壓高度變化系數(shù)，根據(jù)建筑結構荷載規(guī)范GB—50009—2012[13]中B類風場要求分層取值；μs為風荷載體形系數(shù)，對該建筑以高度超過45 m，長寬比為1的矩形截面高層建筑取值，取1.4；ω0為基本風壓(kN/m2)，臺北地區(qū)50年重現(xiàn)期下基本風壓取值0.7 kN/m2。橫風向風荷載標準值按下式計算：

(6)

時域分析法是對結構的動力響應時程，采用慣性風荷載方法計算建筑等效靜力風荷載(ESWLS)，在未考慮流固耦合的數(shù)值風洞中，在剛性模型計算中的每個時間步對大廈的表面風壓進行積分得到大廈每層的風荷載時程數(shù)據(jù)，并將其施加在臺北101大廈三維有限元模型上，在時域內用逐步積分法求解微分方程，得到結構的動力響應時程，在考慮流固耦合的情況下，結構的動力響應時程數(shù)據(jù)可直接由氣彈模型計算得出。等效靜力風荷載可表示為：

(7)

兩種計算方法與風洞試驗數(shù)據(jù)對比見表3。從數(shù)據(jù)對比結果可以看出，順風向上考慮雙向流固耦合與未考慮流固耦合的建筑等效靜力風荷載較為一致，二者與風洞試驗結果相差不大；橫風向和扭轉方向上，未考慮流固耦合的等效風荷載略小于高頻天平測力風洞試驗結果。這是因為數(shù)值模擬中未對建筑周邊建筑進行考慮，風洞試驗中考慮了主體建筑600 m半徑范圍內所有主要建筑物，且風洞試驗中風場的湍流特性也與數(shù)值模擬存在一定差別；值得注意的得是考慮雙向流固耦合后，橫風向等效靜力風荷載則大大減小。此外，數(shù)值模擬計算結果與規(guī)范計算結果有一些差別，這是因為規(guī)范是對結構進行簡化計算，且只考慮結構一階響應，對建筑層間結構引起的風場湍流特性與高階振型影響無法充分考慮。

表3 各工況等效靜力風荷載(結構阻尼比=2%)

4.4 位移和加速度響應

計算與現(xiàn)場實測得到的大樓最高居住層第89層處，X及Y軸方向上的加速度功率譜對比如圖10所示。從圖10中可看出，計算的加速度功率譜與實測譜吻合較好。

(a)x方向加速度功率譜

通過有限元瞬態(tài)動力分析得到臺北101大廈最高居住層(第89層)在50年重現(xiàn)期[14]風速條件下的位移時程與加速度時程曲線，50年重現(xiàn)期風速條件下建筑物最高居住層加速度與位移時程曲線，如圖11所示。

(a)順風向加速度

從圖11可以看出，計算穩(wěn)定后大廈最高居住層的橫風向加速度響應占據(jù)主導地位?？紤]流固耦合作用時順風向最大加速度為0.11 m/s2，最大位移為0.56 m，橫風向最大加速度為0.15 m/s2，最大位移為0.32 m；而未考慮流固耦合時順風向最大加速度為0.18 m/s2，最大位移為0.59 m，橫風向最大加速度為0.86 m/s2，最大位移為0.67 m。在順風向上，ANSYS有限元模型的加速度與位移響應均略大于workbench有限元模型，這是由于有限元模型背風面分布密集的分離渦對建筑產(chǎn)生了負風壓。在橫風向上，ANSYS有限元模型在加速度和位移響應均遠大于workbench有限元模型，這是由于有限元模型中緊密附著在建筑兩側的分離渦對建筑造成的巨大橫向力。對比兩組模型的加速度時程數(shù)據(jù)可以看出，總體上看，ANSYS有限元模型最大位移、最大加速度比workbench有限元模型結果大30%以上。由計算可知，臺北101大廈最高居住層居民在50年重現(xiàn)期風速條件下會感到不適，但在可接受范圍之內，滿足在工程實踐中的要求。

5 結 論

本文以臺北101大廈為研究對象，基于結構動力學特性相似原理，利用 Workbench平臺建立大廈的有限元模型并進行流固耦合數(shù)值模擬，并將計算結果與現(xiàn)場實測以及風洞試驗的相應數(shù)據(jù)進行了對比，以驗證數(shù)值風洞的有效性。監(jiān)測振動穩(wěn)定后400 s內建筑物最高居住層的位移與加速度時程等風致響應，并將其與未考慮流固耦合的ANSYS三維有限元模型進行對比，主要結論如下：

(1)由剛性模型的數(shù)值風洞與實體風洞測得等效風荷載對比可知，采用DSRFG湍流生成方法的DDES模型具有較高的準確性與可信度。在非定風場的模型中，DDES模型能夠較好的預測出建筑物側面與背后的旋渦分布現(xiàn)象，對復雜風場有足夠強的預測能力。

(2)臺北101大廈超高層建筑振動以橫風向為主。在考慮流固耦合效應后，建筑順風向位移及加速度響應變化較小，而建筑物最高居住層(89層)橫風向位移響應與加速度響應均明顯減小，這表明考慮流固耦合效應在對該棟建筑進行數(shù)值風洞模擬時是必要的。

(3)考慮流固耦合情況下，建筑模型橫風向等效風荷載顯著降低，順風向等效風荷載則保持不變。在雙向流固耦合作用下，來流將在建筑背風處形成渦量較大的脫落旋渦，可能對下游毗鄰建筑風穩(wěn)定性造成影響。