張亮亮，吳波，楊陽，劉連杰

(重慶大學 土木工程學院，重慶 400045)

山區(qū)橋址處CFD計算域的選取方法

張亮亮，吳波，楊陽，劉連杰

(重慶大學 土木工程學院，重慶 400045)

山區(qū)橋址處風場具有較強的隨機性與不確定性，在選取其CFD計算域時，盲目參考已有的工程可能會造成較大的模型誤差，或增加大量的計算開銷。為解決該問題，提出了山區(qū)復雜地形CFD計算域選取的方法與步驟，驗證了其準確性。具體方法為：設置一個大范圍的基準計算域進行初算，通過后處理軟件繪制平均風壓系數(shù)極差隨高度的變化曲線、壁面附近水平面上的靜壓偏差等值云圖以及橫風向各平面靜壓值與對應邊界面的均方根差值曲線，分別篩選出基準計算域在高度方向、順風向及橫風向上對整體風場貢獻可忽略的區(qū)域，余下部分則可用于該區(qū)風場的求解。

計算域設置；山區(qū)地形；計算流體力學；模型誤差；求解效率

在山區(qū)建設大跨度橋梁時，風環(huán)境參數(shù)對橋梁結(jié)構(gòu)的設計與安全評價具有重要的參考價值，是必須首要解決的問題之一[1-2]?，F(xiàn)有規(guī)范對山區(qū)風特性的描述較少，風洞實驗難以模擬大范圍區(qū)域的風場，現(xiàn)場實測又難以捕捉大范圍的風特性數(shù)據(jù)，而計算風工程(CFD)具有周期短、數(shù)據(jù)全面、費用低等優(yōu)點[3-4]。隨著計算機技術的發(fā)展，CFD在實際工程中的應用越來越廣泛，其模擬精度也得到了進一步的提高[5]。

利用CFD技術進行風場模擬時，首先要進行計算域的設置。CFD數(shù)值模擬的誤差主要包括模型誤差、離散誤差與迭代誤差[6]。其中，計算域設置不合理所引起的模型誤差是影響計算精度的首要原因，應盡量消除。若計算域設置過小，則模型誤差較大，即使將離散誤差、迭代誤差控制到最低，也無法準確反映實際的風場特性[7]。Fujiwara等[8]發(fā)現(xiàn)，不同大小的計算域?qū)ν晃恢玫挠嬎憬Y(jié)果差異很大，認為模型邊界應遠離計算點足夠遠。但是，若將計算域設置太大，雖可降低模型誤差，但卻大大增加了計算開銷，降低了求解效率。尤其在進行高精度、多工況的實際模擬時，過大的計算域?qū)τ谟嬎阒芷诘暮馁M是難以估量的[8]。如何設置大小合適的計算域，兼顧計算精度與求解效率，是筆者主要解決的問題。目前已有一些關于CFD計算域設置的研究或建議[9-12]。崔利民等[9]針對一個雙向?qū)ΨQ、孤立的正弦山丘提出了山體地形下低矮房屋數(shù)值風洞模擬的計算域設定方法；Franke等[12]給出了群體建筑物CFD風場的豎向、側(cè)向與順風向長度建議值。然而，這些研究均是針對孤立山丘或規(guī)則建筑物的風場，其風場分布對稱、較為規(guī)則且有跡可循。而山地地形在大范圍區(qū)域內(nèi)包含不計其數(shù)的山丘、溝壑、河流、房屋等，風場具有較強的隨機性與不確定性，不能簡單地參照已有工程的經(jīng)驗進行設置，針對規(guī)則風場所提出的方法或建議值也不適用。有經(jīng)驗的學者在進行山區(qū)CFD模擬時，往往先進行試算[13-14]。筆者以地處長江交匯口山地地貌的大寧河特大橋為依托，提出山區(qū)橋址處復雜地形CFD計算域設置的一般方法與步驟。

1 工程概況

重慶大寧河特大橋是G42重慶段二期高速公路的重要工程，橋型為上承式無鉸鋼桁拱橋，跨度400 m，矢高80 m，橋軸線垂直穿過大寧河和兩岸山坡。橋址區(qū)位于大寧河與長江交匯口附近，地表為緩、陡相間的折線型斜坡，坡度55°～75°，屬于典型的山地地貌。橋位范圍內(nèi)最大地面標高為518 m左右，河底標高僅為90 m，其相對高差達428 m，切割深度大。

現(xiàn)場為期2 a的實測數(shù)據(jù)表明，橋址處的主導風向為東北風和東南風，100 a重現(xiàn)期的基本風速為26.8 m/s。

2 計算域的設置方法

2.1 基準計算域

2.1.1 計算域大小 定義如下參數(shù)：H為壁面最高點離計算域頂面的高度，LE、LW、LS、LN分別為橋址中心距計算域東、西、南、北邊界面的距離(圖1)。

山區(qū)橋址處的CFD計算域一般為長/寬5～15 km、高1～3 km[13-15]的長方體區(qū)域。先設置一個較大范圍的基準計算域進行初算，將橋位附近有代表性的地貌全部包含在內(nèi)，其大小為25 km(橋軸方向)×20 km(橋軸法向)×6 km(豎直方向)，即LE=LW=12.5 km，LS=LN=10 km，H=6.0 km?；鶞视嬎阌虻淖饔脼椋簭钠溆嬎憬Y(jié)果中提取相關指標，作為設置H、LE、LW、LS、LN的依據(jù)；將其解作為標準解，計算各檢驗計算域的模型誤差。

圖1 基準計算域及定義的參數(shù)Fig.1 Reference computational domain and its parameters

使用Google Earth獲取地形底面的高程數(shù)據(jù)，取樣間隔30 m，共計獲得584 714個離散高程點。將其導入逆向工程軟件Imageware，擬合四階地形曲面。然后將曲面導入網(wǎng)格劃分平臺Gambit，形成計算域。設置邊界條件：入流面取速度入口Velocity-inlet，出流面為自由流Outflow，地形底面為Wall壁面，其余為對稱邊界Symmetry[13]。

2.1.2 網(wǎng)格劃分與求解 為便于基準計算域與后續(xù)工作中檢驗計算域的對比，所有計算域均分塊為內(nèi)部區(qū)域與外部區(qū)域。內(nèi)部區(qū)域是風場計算所重點關心的區(qū)域，覆蓋橋址中心附近5 km×4 km的范圍，壁面劃分60 m尺度的三角形網(wǎng)格；外部區(qū)域的壁面劃分90 m尺度的非結(jié)構(gòu)化網(wǎng)格。高度方向上，靠近壁面的第一層網(wǎng)格厚度10 m，逐漸加大，增長因子1.1，最大尺度150 m?；鶞视嬎阌蚬矂澐掷庵w網(wǎng)格單元3 542 798個。

將網(wǎng)格文件導入Ansys Fluent 15.0，選用全隱式分離求解器，時間、空間離散均采用二階差分格式，壓力與速度耦合選用SIMPLE算法。本工作的目的在于選擇計算域大小，尚未進入風場特性的高精度模擬的階段。因此，為提高工作效率，采用層流(laminar)模型，來流設置為均勻流，風速取本區(qū)的基本風速26.8 m/s。

2.2 高度H的設置

在模擬孤立建筑物的風場時，一些學者[10-11]將平均風壓系數(shù)極差ΔCp(z)作為設定計算域高度H的依據(jù)，認為ΔCp(z)越大，該高度的水平平面對建筑物周圍風場的影響越大；反之，影響越小，即該平面已不在建筑物的風場區(qū)。山地地貌的風場雖復雜得多，但在豎直方向上具有和孤立建筑物風場同樣的特征[10]。因此，將ΔCp(z)引入到為山地風場，作為計算域高度設置的依據(jù)?；鶞视嬎阌虻摩p(z)計算式為

(1)

式中：Cp為測點平均風壓系數(shù);pi為測點靜壓；p∞為參考點(計算域頂面)靜壓；vH為參考點風速；ρ為空氣密度；Cp,max(z),Cp,min(z)分別為基準計算域Z高度水平面上平均風壓系數(shù)的最大值、最小值，二者差值即該平面的平均風壓系數(shù)極差ΔCp(z)。

圖2給出了基準計算域的ΔCp(z)曲線?？梢姡p隨離地高度Z的增加而逐漸減小，即，離壁面越高的平面受地形的影響越小。對于Z=3.5 km以上的平面，其ΔCp均低于0.01，遠小于Z=3.5 km以下的部分，而且其曲線斜率接近于0，幾乎不再變化。說明基準計算域Z=3.5 km以上的區(qū)域?qū)φw風場的貢獻微弱，已不在壁面風場區(qū)范圍內(nèi)。

圖2 基準計算域的平均風壓系數(shù)極差ΔCp(z)曲線Fig.2 ΔCp(z)curve of reference computational domain

2.3 順風向長度LE、LW的設置

崔利民等[9]以Y軸中心線(順風向)上的相對靜壓變化曲線作為設置計算域長度的指標。但其計算模型是一個雙向?qū)ΨQ的、孤立的正弦山丘，其靜壓分布關于Y軸對稱。然而對于山地地貌的風場而言，Y軸中心線附近的風場與邊界附近風場相去甚遠，崔利民的方法難以全面反應水平面上風場的分布。

CFD模擬風場時，順風向的兩個邊界(x方向，本例對應LE、LW)分別為速度入口與壓力出流。為評價計算域東、西面部分對整體風場的貢獻，經(jīng)試算與篩選，發(fā)現(xiàn)靜壓偏差δP可綜合反映風場在x方向上的分布，其計算式為

(2)

為保持出流面的一致，避免額外誤差，在評價東面部分對整體的影響時，以東面邊界為入流面；同樣，在評價西面部分的影響時，以西面邊界為入流面。

由圖2可知，離壁面越高的平面對風場的影響越小。偏于保守地考慮，選擇壁面最高點附近的水平面(Z=0.5 km)為代表，在后處理軟件Tecplot中繪出其δp(z)的等值云圖(圖3)。

圖3 基準計算域Z=0.5 km平面的靜壓偏差δp等值云圖(大于均方根值S(δp)者)Fig.3 δp-contour of reference computational domain at Z=0.5km plane

圖3給出的是大于均方根值S(δp)的部分——貢獻大的區(qū)域。就圖(a)中x=8.5～12.5 km、圖(b)中x=-8～-12.5 km的區(qū)域而言，其δp幾乎均低于S(δp)，說明基準計算域以東8.5～12.5 km、以西8～12.5 km的區(qū)域?qū)φw的貢獻很小。

2.4 橫風向?qū)挾萀N、LS的設置

CFD模擬風場時，橫風向的兩個邊界面通常設置為對稱邊界Symmetry(y方向，本例對應LN、LS)。其上所有物理量的梯度(grad(φn))均為0。若要以其他面代替現(xiàn)有邊界，則該面上各物理量的解必須與現(xiàn)有邊界的解十分接近。

以0.5 km的增量，分別在基準計算域y=2.0～9.5 km的各平面上以200 m的間距設置監(jiān)控點，每個面共計2 500個；考察這些監(jiān)控點的靜壓與北邊界面(y=10.0 km)上對應點的差值，差值的均方根計算式為

式中：pi,y為y平面上第i個監(jiān)控點的靜壓；pi,0為北邊界面上對應點的靜壓。

同樣地，在y=-2.0～-9.5 km的各平面上設置監(jiān)控點，考察其與南邊界面(y=-10.0 km)上對應點的差值，均方根差值的計算方法同式(3)，不過pi,0對應于南邊界面。

Sp(y)曲線見圖4。由圖4(a)可見，y=-8.0～-9.5 km各平面的Sp均低于0.5%，曲線斜率接近于0。說明將南面邊界移動至計算域中心以南8.0～9.5 km的區(qū)域所造成的誤差是可以忽略不計的。同樣，由圖4(b)，在不顯著增加計算誤差的情況下，可將北面邊界移動至計算域中心以北6.0～9.5 km的區(qū)域。

圖4 基準計算域各y平面與南/北邊界的靜壓均方根差值Sp(y)Fig.4 Sp(y) curve of reference computational domain

實際上，y=-6.5～-8.0 km的各平面的Sp也低于0.5%，在實際工程中基本上可忽略不計。但偏于保守地認為，則以接近于0的斜率趨近于0的曲線段才是對整體貢獻可以忽略的部分。顯然，y=-6.5～-8.0 km曲線段的斜率并不接近于0，因此，不認為該部分對整體的貢獻可以忽略。

按這3種方法分別進行高度、順風向長度、橫風向?qū)挾鹊脑O置，舍棄掉對整體貢獻可被忽略的區(qū)域，所選定的用于實際求解的計算域為：H=3.5 km、LE=8.5 km、LW=8.0 km、LS=8.0 km、LN=6.0 km，其體積僅為基準計算域的0.27倍。

2.5 計算域設置步驟

山區(qū)橋址處CFD計算域的設置步驟為：

1)設置一個大范圍的基準計算域，將本區(qū)域內(nèi)所有具代表性的地貌包含在內(nèi)；獲取本區(qū)氣象數(shù)據(jù)，以本區(qū)的盛行風向為入口，以相對于實際求解較疏的網(wǎng)格、較低階的計算方法進行快速初算；

2)繪制基準計算域的平均風壓系數(shù)極差隨離地高度Z的變化曲線ΔCp(z)，以曲線上趨近于零值的拐點高度作為實際計算域的高度；

3)分別以順風向(盛行風)上的兩個邊界面為速度入口進行求解，繪制其壁面附近水平面的靜壓偏差云圖δp，舍棄掉低于其均方根值S(δp)的部分，余下部分的長度即為實際計算域的順風向長度；

4)在橫風向上的各平面設置監(jiān)控點，考察各監(jiān)控點靜壓與對應兩個邊界面的差值，分別繪制其均方根誤差曲線Sp(y)，以趨近于零值的拐點作為實際計算域的橫風向長度。

3 算法驗證

計算域高度、順風向長度、橫風向?qū)挾鹊脑O置方法是在不斷試算、反復篩選的基礎上而提出的，必須對其進一步的驗證。驗證思路為：按上述方法已經(jīng)評價出對整體風場貢獻可忽略不計的區(qū)域，分別舍棄掉這些區(qū)域，建立若干個檢驗計算域。對比檢驗計算域的解與標準解(基準計算域的解)的差異。若差異可忽略不計，說明舍棄掉區(qū)域不會增強顯著的模型誤差，即驗證了上述方法的準確性；若差異不可忽略，則說明上述方法并不適用。

3.1 檢驗計算域

以0.5 km為增量，分別將H、LE、LW、LS、LN作為唯一變量，設置H=1.0～5.5 km的10個、LE=2.5～12.0 km的20個、LW=2.5～12.0 km的20個、LN=2.0～9.5 km的18個、LS=2.0～9.5 km的18個(共計86個)檢驗計算域。

實測數(shù)據(jù)表明，該區(qū)的主導風向為東北風與東南風。因此，除LE=2.5～12.0 km的20個檢驗計算域外，其余檢驗計算域均以西面邊界入流，形成由西向東的流向。與之對比的基準計算域也取西面邊界為入流面進行求解，其網(wǎng)格劃分、求解方法均與各檢驗計算域完全一致，從而保證了H、LW、LS、LN分別作為變量的唯一性；對于LE=2.5～12.0 km的20個檢驗計算域，若以西面邊界入流，則各計算域與基準計算域的出流面不一致，可能導致額外的誤差。因此，這20個檢驗計算域均以東面邊界入流，形成由東向西的流向，與之對比的基準計算域也取東面邊界入流。

3.2 五參數(shù)對計算誤差的影響

在各檢驗計算域內(nèi)部區(qū)域Z=1/4Hj、2/4Hj、3/4Hj、Hj(Hj為檢驗計算域的高度)處的水平面上以100 m的間隔設置2 000個監(jiān)控點，每個檢驗計算域共計4×2 000=8 000個監(jiān)控點。對比檢驗計算域各監(jiān)控點的順風向、橫風向、豎向風速與基準計算域?qū)c的差值，其均方根誤差的計算式為

(4)

式中：vi,j為檢驗計算域第i個監(jiān)控點的風速；vi,0為基準計算域?qū)c的風速。若某檢驗計算域的順風向速度vx、橫風向速度vy、豎向速度vz的誤差(Svx、Svy、Svz)均趨近于0，且所在曲線段斜率接近于0，則認為其模型誤差可忽略不計。

Svx、Svy、Svz隨H、LE、LW、LS、LN的變化曲線如圖5。由圖5可見，模型誤差可忽略不計的檢驗計算域有：H=3.5～5.5 km的5個，LE=8.5～12.0 km的8個，LW=8～12.0 km的9個，LN=6.0～9.5 km的8個以及LS=7.5～9.5 km的5個。說明H=3.5～6.0 km、LE=8.5～12.5 km、LW=8～12.5 km、LN=6.0～10.0 km、LS=7.5～10.0 km這些區(qū)域是可以被舍棄的。

而按照該設置方法，選定的計算域為：H=3.5 km，LE=8.5 km，LW=8.0 km，LS=8.0 km，LN=6.0 km?？梢姡齃W偏于保守外(多估計了0.5 km)，其余參數(shù)均與與Svx、Svy、Svz誤差曲線的結(jié)果準確吻合，從而驗證了該方法的準確性。

圖5 高度(圖(a))、順風向(圖(b)、(c))及橫風向長度(圖(d)、(e))對計算誤差Svx、Svy、Svz的影響曲線Fig.5 Impact of H、LE、LW、LSand LNon calculation error

4 結(jié) 論

1)山區(qū)風場具有較強的隨機性與不確定性，在選取其計算域時，不可盲目參照已有的工程經(jīng)驗。

2)提出了山區(qū)橋址處CFD計算域的選取方法，驗證了其準確性。

3)進行了一些偏于保守的處理：認為曲線上以近于0斜率趨近于零值者才是對整體貢獻可忽略的部分；工程運用時可視精度需要選擇貢獻大小的分界；以受壁面影響最大的水平面代表計算域在順風向的分布，工程運用時可視精度需要選擇離壁面稍遠一些的平面。

4)計算域的選取工作尚未涉及實際求解的高精度需求，因此,可以較疏的網(wǎng)格、較低階算法進行快速初算；再通過簡單處理即可準確、直觀地篩選出實際求解的計算域，并不會增加額外的計算負擔。

5)本例所選定的實際計算域體積僅為基準計算域的0.27倍，在進行密網(wǎng)格、高階算法、復雜湍流模型、多工況的實際求解時，可大大減少工作量，縮短計算周期。

[1] 周志勇，肖亮，丁泉順，等.大范圍區(qū)域復雜地形風場數(shù)值模擬研究[J].力學季刊，2010，31(1)：101-107. Zhou Z Y,Xiao L,Ding Q S,et al. Numerical simulation study of wind environment for the flow around large region with complex terrain [J]. Chinese Quarterly of Mechanics,2010,31(1): 101-107. (in Chinese)

[2] 賀德馨.我國風工程研究現(xiàn)狀與展望[J].力學與實踐，2002，24(4)：10-19. He D X. The status and vista of wind engineering studies in China [J]. Mechanics in Engineering,2002,24(4): 10-19. (in Chinese)

[3] 項海帆.進入21世紀的橋梁風工程研究[J].同濟大學學報:自然科學版,2002,30(5): 529-532. Xiang H F. Study on bridge wind engineering into the 21st century [J]. Journal of Tongji University: Natural Science,2002,30(5): 529-532. (in Chinese)

[4] Ding Y L,Zhou G D,Li A Q,et al. Statistical characteristics of sustained wind environment for a long-span bridge based on long-term field measurement data [J]. Wind and Structures,2013,17(1): 43-68.

[5] Blocken B,Stathopoulos T,Carmeliet J,et al. Application of computational fluid dynamics in building performance simulation for the outdoor environment: An overview[J]. Journal of Building Performance Simulation,2011,4(2): 157-184.

[6] 曾鍇，汪叢軍，黃本才，等.計算風工程中幾個關鍵影響因素的分析與建議[J].空氣動力學學報，2007，25(4)：504-508. Zeng K,Wang C J,Huang B C,et al.Suggestion and analysis of several key factors in computational wind engineering [J]. Acta Aerodynamica Sinica,2007,25(4): 504-508. (in Chinese)

[7] Zhu Z W,Liu Z Q. CFD prediction of local scour hole around bridge piers [J]. Journal of Central South University of Technology: English Edition,2012,19(1): 273-281.

[8] Fujiwara T. Application of statistical system identification method to the analysis of wind-induced current and wind set-up [J]. Coastal Engineering Journal,1985,28: 117-123.

[9] 崔利民，彭興黔，時凌琳，等.山體地形下低矮房屋數(shù)值風洞模擬的計算域設定[J].華僑大學學報:自然科學版，2010，31(4)：463-467. Cui L M,Peng X Q,Shi L L,et al. Computational domain setting about numerical wind tunnel in the simulation of the low-rise housing in the mountain terrain [J]. Journal of Huaqiao University: Natural Science,2010,31(4): 463-467. (in Chinese)

[10] Jakob N,Gurevych I. Extracting opinion targets in asingle-and cross-domain setting with conditional random fields [C]// Conference on Empirical Methods in Natural Language Processing,Cambridge,MA,United States,2010: 1035-1045.

[11] Hu W L,Zhang Y H,Wang L B. Numerical simulation on turbulent fluid flowand heat transfer enhancement of a tube bank fin heat exchanger with mounted vortex generators on the fins [J]. Journal of Enhanced Heat Transfer,2010,18(5): 361-374.

[12] Franke J,Hellsten A,Heinke K,et al. The COST 732 best practice guideline for CFD simulation of flows in the urban environment: A summary [J]. International Journal of Environment and Pollution,2011,44: 419-427.

[13] 祝志文，張士寧，劉震卿，等.橋址峽谷地貌風場特性的CFD模擬[J].湖南大學學報:自然科學版，2011，38(10)：13-17. Zhu Z W,Zhang S N,Liu Z Q,et al. CFD simulation of wind field at bridge site on gorge terrain [J]. Journal of Hunan University: Natural Science,2011,38(10): 13-17. (in Chinese)

[14] 李永樂，蔡憲棠，唐康，等.深切峽谷橋址區(qū)風場空間分布特性的數(shù)值模擬研究[J].土木工程學報，2011，44(2)：116-122. Li Y L,Cai X T,Tang K,et al. Study of spatial distribution feature of wind fields over bridge site with a deep-cutting gorge using numerical simulation [J]. China Civil Engineering Journal,2011,44(2): 116-122. (in Chinese)

[15] Huang W F,Xu Y L. Prediction of typhoon design wind speed and profile over complex terrain [J]. Structural Engineering and Mechanics,2013,45(1): 1-18.

(編輯 胡英奎)

Method of setting up the wind field of a mountain bridge site

ZhangLiangliang,WuBo,YangYang,LiuLianjie

(School of Civil Engineering,Chongqing University,Chongqing400045,P. R. China)

The wind field of a mountain bridge site always shows a strong randomness and uncertainty. As a result,setting up its computational domain in CFD simulation by a simply reference to the existing experience causes a large amount of model errors and huge additional computational overhead. To solve this problem,a method to the selection of CFD computational domain of complex terrains is explored and verified. To start the procedure,a large reference computational domain should be chosen and preliminary solved. Then,its solution should be post-processed to draw the following curves or contour images: range of mean pressure coefficient align with height,static pressure deviation contour of the horizontal plane near the top of the bottom surface,the root-mean-square error of the static pressure of the crosswind planes compared to the corresponding edge surface. According to these curves and contour images,areas parts which have slight contribution on the overall wind field are figured out and abandoned,and the selected computational domain is constitutive of the rest parts of the reference domain.

computational domain settings; mountain terrain; CFD; model error;computational efficiency

10.11835/j.issn.1674-4764.2015.05.002

2015-07-05 基金項目：國家自然科學基金(51578098)

張亮亮(1956-)，教授，博士生導師，主要從事橋梁力學性能分析、橋梁抗風研究，(E-mail)zll200510@126.com。

Foundation item：National Natural Science Foundation of China(No.51578098)

V211.3;O368

A

1674-4764(2015)05-0011-07

Received:2015-07-05

Author brief：Zhang Liangliang (1956-),professor,doctoral supervisor,main research interests:bridge mechanics performance analysis & wind resistance of bridge,(E-mail) zll200510@126.com.