計算風工程中k－ε模型的邊界條件研究

郭 威，張 悅

（1.福建省交通規(guī)劃設計院，福建福州 350004；2.長安大學公路學院，陜西西安 710064）

計算風工程中k－ε模型的邊界條件研究

郭 威1，＊，張 悅1，2

（1.福建省交通規(guī)劃設計院，福建福州 350004；2.長安大學公路學院，陜西西安 710064）

從k－ε湍流模型控制方程入手，對兩類滿足平衡大氣邊界層理論要求的入口邊界條件進行了詳細分析，比較了二者的理論差異。然后以我國公路橋梁抗風規(guī)范中建議的A類風場為來流條件，定義了相應的k－ε模型常數(shù)與壁面條件，使用兩類邊界條件進行數(shù)值模擬。計算結果表明：兩類邊界條件結合相應的模型常數(shù)與壁面條件，均能在數(shù)值模擬中構建基本滿足規(guī)范要求的平衡大氣邊界層，但各有其適用范圍。在此基礎上，對模型常數(shù)C1和湍動能耗散率與數(shù)值模擬結果之間的關系進行了研究，采用修改模型常數(shù)與邊界條件的方法，改善了數(shù)值模擬結果。

k－ε湍流模型；平衡大氣邊界層；邊界條件

0 引 言

正如風洞試驗需要在風洞內重現(xiàn)充分發(fā)展的大氣邊界層，在計算風工程中，構建水平均一的（horizontally homogenous）大氣邊界層是保證數(shù)值模擬結果準確可靠的關鍵前提。所謂水平均一，是指來流物理量沿流向梯度為0，即計算域入口處的來流在無干擾條件下到達出口保持不變。這樣的邊界層可稱為平衡邊界層。構建平衡邊界層，需要來流邊界、湍流模型和壁面條件三者相互協(xié)調。由于問題十分復雜，通常難以給定一個既能滿足理論要求，還能與場地實測數(shù)據(jù)相符的來流邊界條件，因此在數(shù)值模擬中實現(xiàn)平衡大氣邊界層十分不易。

邊界層在計算域入口段產生的物理量變化，將與結構物對來流的擾動作用疊加，共同體現(xiàn)在計算結果中。同時，邊界層的物理量變化，將使得作用在結構物上的風剖面與入口處的預設條件不符。兩種效應降低了數(shù)值模擬的可靠性與準確性。

對于如何構建平衡大氣邊界層的問題，諸多學者進行了研究。Blocken［1］等指出入口風速剖面中包含的粗糙長度應該與壁面函數(shù)相匹配，否則來流在計算域內將發(fā)展成一個新的平衡邊界。然而僅按照其建議的方法修改壁面參數(shù)，來流在計算域內依然發(fā)生變化。Richards和Hoxey［2］（1993）從標準k－ε模型方程出發(fā)，基于對平衡大氣邊界層的一系列假設，推導出一組邊界條件及其相應的模型常數(shù)，但是其給出的湍動能剖面只能是常數(shù)，與實際不符。在此基礎上，Yang［3］等推導出標準k－ε模型方程的通解，給出了一組符合實測的邊界條件，但Yang沒有給出模型常數(shù)設置方法，給出的邊界方程不能直接得到相應的模型常數(shù)。

本文首先對Richards與Yang提出的兩類邊界條件進行理論對比。然后應用兩種方法，對《公路橋梁抗風規(guī)范》（JTG－TD60－01－2004）中建議的A類風場進行足尺數(shù)值模擬［4］。通過對兩個算例結果的比較，分析了兩類邊界條件的差異及各自的適用條件。根據(jù)模擬結果，本文對模型常數(shù)和邊界條件做了調整，改善了平衡大氣邊界層的模擬結果。

1 平衡大氣邊界層的邊界分析

1.1 第一類邊界條件

定常流的標準k－ε模型控制方程為：

式中，湍流粘性的定義為：

Richards與Hoxey（1993）首先對平衡大氣邊界層做了以下假設：（1）邊界層內豎向速度為0，對于二維流場，有邊界層內水平切應力沿高度方向相等，即：

其中，τ為大氣內部切應力，u＊為剪切速度。在大氣邊界層局部區(qū)域內，該假設近似成立。

其次，研究者對計算域入口處的平均風速剖面使用對數(shù)律表達式擬合：

式中K為馮卡門常數(shù)，通常取值為0.4～0.42。

在風速u、湍動能k與湍流耗散率ε水平方向梯度為0的條件下，由方程（1）可以得到μtGk＝ρε，即：

由于μ遠小于μt，可將其忽略。將方程（3）～（5）聯(lián)立求解，不難得出Richards等建議的入口湍動能k及其耗散率ε表達式：

這組由k－ε湍流模型方程（1）推導出的邊界條件表達式（5）、（7）和（8）自然能夠滿足方程（1）。再將其代入模型方程（2）中，要使得方程成立，模型常數(shù)需滿足關系式：

1.2 第二類邊界條件

第二類邊界條件由Yang等推導，沒有使用方程（4），而是將方程（3）、（5）和（6）直接代入到方程（1）中得到：

分離變量后得該微分方程通解：

代入方程（6）得到：

方程（5）、（12）和（13）就是Yang等建議的入口邊界條件。

將方程（5）、（6）和（8）代入到方程（2）中，方程左邊為0，而方程的右邊：

可見Yang推導的邊界條件沒有直接對標準k－ε湍流模型中的常數(shù)關系進行限制，但在入口邊界條件確定后應驗算F1是否為0。若F1不等于0，則模型方程（2）不能滿足，來流將在計算域內發(fā)展成一個滿足模型方程的新邊界層。

1.3 兩類邊界條件對比

從以上兩類邊界條件的推導過程可以看出，二者的最大差異在于是否假設邊界層內切應力相等。該假設僅在大氣邊界層的小范圍內近似成立，與真實情況略有出入。使用該假設，引入方程（4），使得第一類邊界條件中的湍動能k沿高度方向保持不變，與實際大氣邊界層不符。但第一類邊界條件中，平均風速剖面中的剪切速度u＊物理意義明確，湍動能耗散率ε表達式簡單，同時標準k－ε模型中的常數(shù)也有確定的關系。

第二類邊界條件沒有使用方程（4），而是使用與第一類邊界條件相同的平均風速剖面表達式，推導出標準k－ε模型控制方程（1）的通解（9）。該解能夠通過改變算式中的參數(shù)A和B，擬合湍動能k隨高度變化的數(shù)值，使得計算域入口的湍動能設置更接近實際情況。然而這樣的改進也引來新問題：將方程（2）整理后得到的等式

并不能得出確定的模型常數(shù)關系。此外，將方程（3）、（6）代入切應力表達式（4），得到：

由此可以看出，隨著湍動能k沿高度的變化，大氣邊界層內的水平方向切應力沿豎直方向的分布是變化的。因此，平均風速剖面表達式中的u＊值也應沿高度變化。

2 算例驗證

2.1 風場參數(shù)確定

按照我國《公路橋梁抗風規(guī)范》推薦的A類風場設置計算域入口平均風速剖面與湍動能剖面。

首先假定風速的參考高度為10m，設該高度處的風速為10m／s，則A類風場的平均風速剖面為：

“規(guī)范”中只給出A類風場不同高度的流向湍流強度，需要換算才能得到相應的湍動能剖面。在只關心流向湍流強度的前提下，可進一步假定豎向脈動風速為0，按關系式進行換算。結果如表1所示。

表1 A類風場湍動能換算值Table 1 Turbulent kinetic energy conversion value of class A wind field

2.2 入口邊界條件設置

取馮卡門常數(shù)K＝0.42，則剪切速度u＊＝平均風速剖面表達式為

第一類邊界條件中的湍動能k值只能為常數(shù)，這里取把K與Cμ代入等式（9）中，得到σε＝3.68。

圖1 湍流強度擬合結果Fig.1 Turbulence intensity fitting results

模型常數(shù)取值與第一類邊界條件相同，得到湍動能耗散率剖面為

兩類邊界條件設置如表2所示。

2.3 計算模型

采用二維模型驗證兩類邊界條件的適用性，矩形計算域大小為400m×5000m（y×x），計算域中不設置任何結構物。采用結構化網格對計算域進行離散，水平向網格間距為10m，豎向首層網格高度取0.32m，漸變率1.08，首層網格中心高度zp＝0.16m。網格總數(shù)為60 000。采用SIMPILEC算法求解，對流項的離散格式選用QUICK。以各物理量的無量綱殘差低于10－5為計算收斂標準。計算平臺為Fluent6.2。

表2 兩類邊界條件對比Table 2 Comparison on two kinds of boundary conditions

2.4 壁面函數(shù)及其他參數(shù)設置

Blocken［1］與Richards［5］均指出，壁面函數(shù)中粗糙高度Ks與邊界層粗糙長度z0并不相等。二者之間的關系需由平均風速剖面，湍動能剖面與壁面函數(shù)三者共同確定。

Fluent中標準壁面函數(shù)［6］為：

式中，Up為首層網格中心高度處平均風速；yp是首層網格中心高度；v是動力粘度，取值為為壁面粗糙常數(shù)；u＊為剪切速度，等于為經驗常數(shù)，約取9.793。K＋s＝u＊Ks／v，K＋s＞90時，

壁面剪切速度：

在平衡大氣邊界層中u＊＝uτ。此時壁面函數(shù)表達式為：

設Cs＝1，可以得到Ks＝0.092。

第二類邊界條件中的平均風速剖面的形式雖然與第一類邊界條件相同，但剪切速度u＊在不同高度取值不同。首層網格中心高度處的實際摩擦速度應為代入方程（20），整理得：

設Cs＝1，解得Ks＝0.122，即為第二類邊界條件的壁面參數(shù)。

此外，應該注意到空氣分子粘性μ＝1.71× 10－5Ns／m2，空氣湍流粘性為：

按照第一類邊界條件，在計算域頂部，μt＝ 125Ns／m2。湍流粘性比應修改Fluent軟件中湍流粘性比的默認限制值1×105，否則湍動能將在入口附近的計算域內迅速衰減。

其他模型常數(shù)為C1＝1.5，C2＝1.92，σk＝1。

2.5 計算結果

使用兩類邊界條件計算得到的入口與出口平均風速和湍動能的對比如圖2所示。

圖2（a）表明，使用兩類邊界條件計算得到的出口與入口平均風速剖面相差較小。在計算域頂部，出口平均風速低于入口設定，這主要是受計算域頂面的邊界條件影響，可以通過在頂面施加剪應力而改善［78］。在計算域底部，出口處平均風速均略高于入口設定，且第二類邊界條件計算得到的平均風速改變量小于第一類邊界條件。這是由于第二類邊界條件在地面附近設置的湍動能值高于第一類邊界條件，由此推算出的壁面Ks值更高導致的。方平治［9］討論了Ks對平均風速剖面的影響，從中可以看出，增大Ks值能夠減小計算域底部平均風速剖面沿高度方向的梯度。圖2（b）表明使用兩類邊界條件計算得到的出口處湍動能要小于出口處的設置值。同樣，使用第二類邊界條件湍動能的改變量略小于第一類邊界條件。

圖2 兩類邊界條件計算結果對比Fig.2 Comparison of two kinds of boundary conditions calculation results

比較兩個算例可以看出，使用兩類邊界條件得到的計算結果都能近似滿足構建平衡大氣邊界層的要求。第一類邊界條件的設置方法較簡單便捷，但湍動能的設置與實測結果有偏差，僅適宜在湍動能沿高度變化小的情況下使用。第二類邊界條件可以設置湍動能沿高度的變化，更接近實際風場。而且地面附近區(qū)域設置的湍動能值更大，也使得壁面Ks增大，從而減小了計算域底部平均風速剖面沿高度方向的梯度，與入口平均風速剖面更加吻合。因此，第二類邊界條件適宜在計算準確性要求更高及湍動能強度隨高度變化大的情況下使用。

3 對第二類邊界條件的改善

為得到更加平衡的大氣邊界層，針對第二類邊界條件得到的湍動能偏小的情況，本文首先從滿足模型控制方程（2）的要求入手，對F1是否為0進行驗算。從算式（14）可以看出，F(xiàn)1在z較小的壁面附近數(shù)值較大，減小C1將減小F1值，因此可以適當修改C1改善計算結果。

Durbin［10］指出C1合理取值應為：

其有效取值范圍在1.3～1.55之間。在不改變邊界條件與其他模型常數(shù)的前提下，本文嘗試修改C1＝1.45與C1＝1.40。C1不同取值時F1的數(shù)值如圖3所示。

圖3 不同C1取值時F1的數(shù)值Fig.3F1values with diffirentC1

圖3表明，隨著C1的減小近壁面位置F1更接近0。不同C1取值得到的計算結果如圖4所示。

從計算結果可以看出，隨著常數(shù)C1的降低，平均風速的取值幾乎沒有變化。而湍動能的變化有其特點：高空位置（＞200m）的湍動能差異微小，而200m以下位置的湍動能值隨著C1的降低而增大，越接近壁面，增大的數(shù)值越多。

為減少高空位置湍動能的衰減，本文嘗試修改其他邊界條件。注意到即便是第一類邊界條件，當標準k－ε湍流模型兩個控制方程完全滿足平衡大氣邊界層條件時，來流物理量經過計算域后依然發(fā)生變化。這是由于邊界條件的前提假定與計算域內的實際流場的微小差異引起。這樣的誤差難以預估，因此，鑒于以上模擬結果顯示湍動能在出口處小于入口值，本文嘗試將入口處的湍動能耗散率剖面表達式上增加一個系數(shù)，使得

D分別取0.9、0.8和0.7，得到的計算結果如圖5所示。

圖4 不同C1值計算結果Fig.4 Calculation results with diffirent C1

從圖5可以看出，計算域入口處設置湍動能耗散率的變化同樣不影響出口位置平均風速值，而且對出口位置全部高度范圍內的湍動能值都有影響。具體表現(xiàn)為：在入口處湍動能耗散率減小的同時，出口處的湍動能值在全部高度范圍內都有增加，但增加的數(shù)值并不相同。高處的湍動能增加值較大，越接近壁面增加值越小。D取0.7時，高空處的湍動能與入口條件吻合度較好，但在200m以下低風速位置湍動能值高于入口值。D取1與0.9時，出口處湍動能仍小于入口值。取D＝0.8時，計算域出入口的湍動能吻合度較好。

圖5 不同D值計算結果Fig.5 Calculation results with diffirentD

4 結 論

本文首先分析了兩類平衡大氣邊界層的入口邊界條件，并通過算例驗證了其準確性與適用性。進而對邊界條件中的參數(shù)進行了研究，以期改善數(shù)值模擬結果?？偨Y全文，得出如下結論：

1）第一類邊界條件的設置方法相對簡單，但湍動能為定值，只適宜在湍動能沿高度分布均勻的大氣邊界層模擬中使用。第二類邊界條件的設置方法較復雜，但允許湍動能沿高度變化分布，在大氣邊界層模擬中的適用范圍廣泛。

2）減小C1會增大計算域內壁面附近位置的湍動能，但不影響高空處湍動能。同時C1對平均風速剖面影響微小。

3）減小ε能增大計算域內的湍動能值，效果隨高度的增大而增加。ε對平均風速剖面影響微小。

［1］Blocken B，Stathopoulos T，Carmeliet J.CFD simulation of the atmospheric boundary layer：wall function problems［J］.Atmos.Environ.，2007，41（2）：238－252.

［2］Richards P J，Hoxey R P.Appropriate boundary conditions for computational wind engineering models using thek－εturbulence model［J］.Journal of Wind Engineering and Industrial，1993，46－47：145－153.

［3］Yang Wei，Gu Ming，Chen Suqin.A set ofturbulence boundary condition ofk－εmodel for CWE［J］.Acta Aerodynamica Sinica，2005，23（1）：97－102.（in Chinese）楊偉，顧明，陳素琴.計算風工程中k－ε模型的一類邊界條件［J］.空氣動力學學報，2005，23（1）：97－102.

［4］JTG／TG－T D 60－01－2004，Wind－resistant design specification for highway bridges［S］.JTG－／T D 60－01－2004，公路橋梁抗風設計規(guī)范［S］.

［5］Richards P J，Norris S E.Appropriate boundary conditions for computational wind engineering models revisited［J］.J.Wind Eng.Ind.Aerodyn.，2011，99：257－266.

［6］Fluent 6.1Documentation［R］.Fluent Inc.2003.

［7］Yang Wei，Jin Xinyang，Gu Ming，et al.A study on the selfsustaining equilibrium atmosphere boundary layer in computational wind engineering and its application［J］.China Civil Engineering Journal，2007，40（12）：1－5.（in Chinese）楊偉，金新陽，顧明，等.風工程數(shù)值模擬中平衡大氣邊界層的研究與應用［J］.土木工程學報，2007，40（2）：1－5.

［8］Hargreavesa，D M，Wright N G.The use of commercial CFD software to model the atmospheric boundary layer［C］／／The Fourth International Symposium on Computational Wind Engineering，Yokohama，Japan，2006：797－800.

［9］Fang Pingzhi，Gu Ming，Tan Jianguo.Numerical wind fields based on thek－εturbulent models in computational wind engineering［J］.Chinese Journal of Hydroynamics，2010，25（4）：475－483.（in Chinese）方平治，顧明，談建國.計算風工程中基于k－ε系列湍流模型的數(shù)值風場［J］.水動力學研究與進展，2010，25（4）：475－483.

［10］Durbin P A.Separated flow computations with thek－ε－v2model［J］.AIAA Journal，1995，33（4）：659－664.

［11］Svetlana Poroseva，Gianluca Iaccarino.Simulating separated flow usingk－εmodel［R］.Annual Research Briefs 2001，Center for Turbulence Research，Standford University，2001：375－383.

［12］Yang Yi，Gu Ming，Chen Suqin，et al.New inflow boundary conditions for modelling the neutral equilibrium atmospheric boundary layer in computational wind engineering［J］.J.Wind Eng.Ind.Aerodyn.，2009，97（2）：88－95.

［13］Zeng Kai，Wang Congjun，Huang Bencai，et al.Suggestion and analysis of several key factors in computational wind engineering［J］.Acta Aerodynamica Sinica，2007，25（4）：504－508.（in Chinese）曾鍇，汪叢軍，黃本才，等.計算風工程中幾個關鍵影響因素的分析與建議［J］.空氣動力學學報，2007，25（4）：504－508.

［14］Launder B E，Spalding D B.The numerical computation of turbulent flows［J］.Computer Methods in Applied Mechanics and Engineering，1974，3（2）：269－289.

［15］O’Sullivan J P，Archer R A，F(xiàn)lay R G J.Consistent boundary conditions for flows within the atmospheric boundary layer［J］.J.Wind Eng.Ind.Aerodyn，2011，99（99）：65－77.

A study on boundary conditions of k－εmodel in Computational Wind Engineering

Guo Wei1，＊，Zhang Yue1，2

（1.Fujian Communications Planning ＆Designing Institute，F(xiàn)uzhou 350004，China；2.School of High way，Chang’an University，Xi’an 710064，China）

Proceeded from thek－εturbulence model control equation，two kinds of boundary conditions meet the theoretical conditions of horizontally homogenous Atmosphere Boundary Layer（ABL）were studied in detail，and their theoretical differences were compared.And then the class A wind field，which is suggested by the Chinese Wind－resistant Design Specification for Highway Bridges，was simulated as the inlet flow conditions in full scale using these two kinds of boundary conditions，with adaptedk－εmodel constants and wall conditions well defined.The results show that，in numerical simulation，both two kinds of boundary conditions，combined with adapted model constants and wall conditions，are able to construct horizontally homogenous ABL，which basically satisfy the conditions required by the code，but each has its own applicable scope.Based on these works，the relation between model constantC1，turbulent dissipation rateεand the numerical simulation results was studied，and a method to improve the results by modifying model constant andεwas introduced.

k－εturbulence model；homogenous atmosphere boundary layer；boundary conditions

V211.3

Adoi：10.7638／kqdlxxb－2014.0075

0258－1825（2016）04－0530－06

2014－07－11；

2015－01－28

國家自然科學基金（51078038）

郭威＊（1988－），男，福建永泰人，碩士，助理工程師，專業(yè)：橋梁工程.E－mail：guowei－14＠163.com

郭威，張悅.計算風工程中k－ε模型的邊界條件研究［J］.空氣動力學學報，2016，34（4）：530－535.

10.7638／kqdlxxb－2014.0075 Guo W，Zhang Y.A study on boundary conditions ofk－εmodel in Computational Wind Engineering［J］.Acta Aerodynamica Sinica，2016，34（4）：530－535.