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

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

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

1 風(fēng)洞試驗

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

圖1 臺北101大廈風(fēng)洞試驗?zāi)Ｐ?/p>

圖2 風(fēng)洞試驗湍流強(qiáng)度剖面

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月期間，三次臺風(fēng)作用和一次地震作用下采集得到。

(a)傳感器立面布置圖

3 數(shù)值模擬方法

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

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

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

(1)

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

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

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

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

3.2 邊界條件

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

(2)

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

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

圖5 入口湍流強(qiáng)度及積分尺度剖面

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

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

表1 邊界條件設(shè)置

圖7 DSRFG方法生成的脈動風(fēng)速時程功率譜與目標(biāo)譜比較

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

3.3 湍流模型

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

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

(3)

其中，

(4)

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

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

3.4 求解參數(shù)

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

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

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

4.1 流場分析

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

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

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

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

(a)剛性模型y=0

4.2 模態(tài)分析

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

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

4.3 等效靜力風(fēng)荷載

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

規(guī)范中順風(fēng)向風(fēng)荷載標(biāo)準(zhǔn)值ωk按下式計算：

ωk=βzμsμzω0

(5)

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

(6)

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

(7)

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

表3 各工況等效靜力風(fēng)荷載(結(jié)構(gòu)阻尼比=2%)

4.4 位移和加速度響應(yīng)

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

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

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

(a)順風(fēng)向加速度

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

5 結(jié) 論

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

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

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

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