基于SST k－ω 湍流模型的平衡大氣邊界層模擬

胡 朋，李永樂，廖海黎

（西南交通大學橋梁工程系，四川 成都 610031）

0 引言

計算風工程（Computing Wind Engineering，簡稱CWE）是結(jié)構(gòu)風工程中極具前景的一個方向，也是當前國際風工程研究的熱點。CWE的核心內(nèi)容是計算流體動力學（Computational Fluid Dynamics，簡稱CFD），與傳統(tǒng)的風洞試驗相比，CFD 方法具有費用低、周期短、便于模擬真實環(huán)境、描述流場細節(jié)及給出流場定量結(jié)果等優(yōu)點。自20世紀60年代興起以來，隨著數(shù)值求解方法的發(fā)展和計算機硬件技術的進步，CFD 在結(jié)構(gòu)風工程中的應用越來越廣泛。

由于結(jié)構(gòu)風工程研究的對象位于大氣邊界層內(nèi)，研究的重點是鈍體繞流，其流動可視為低速、不可壓縮粘性流體的湍流運動，因此準確模擬大氣邊界層內(nèi)的湍流流動是精確預測結(jié)構(gòu)風荷載的前提，這在CFD 數(shù)值模擬中就首先要求其滿足水平均勻性［1］的條件，即在沒有任何障礙物的簡單邊界層流動下，數(shù)值模擬區(qū)域內(nèi)的平均風剖面和湍流風剖面自入口到出口保持不變，這種滿足水平均勻性條件的邊界層稱為平衡大氣邊界層，邊界層內(nèi)流場的水平均勻性程度直接影響數(shù)值模擬的精度。Blocken等［2］采用k－ε湍流模型模擬了兩平行建筑物內(nèi)的風速分布情況并與相應的風洞試驗作對比，結(jié)果表明兩建筑物中心線上行人高度處風速放大系數(shù)的計算值與試驗值相差達15%左右，其原因主要是由于入口處的風速剖面在流動過程中發(fā)生了改變，即不滿足水平均勻性條件，從而使計算風速放大系數(shù)的參考風速發(fā)生了變化，當根據(jù)風速剖面變化對參考風速進行修正后，風速放大系數(shù)的計算值與試驗值吻合較好，文中還建議在進行數(shù)值模擬之前，應首先對無模型的空流場進行水平均勻性檢驗。嚴格說來，當且僅當計算區(qū)域入口處的來流邊界條件與湍流模型（包括壁面函數(shù)）及地表粗糙特性相協(xié)調(diào)時，才能得到水平均勻性的平衡大氣邊界層［3］，但由于湍流的復雜性，在數(shù)值模擬中精確實現(xiàn)平衡大氣邊界層常存在困難。Richards和Hoxey［1］基于k－ε湍流模型提出了滿足平衡大氣邊界層的幾個原則，給出了一組來流邊界條件，并對湍流模型參數(shù)的取值進行了探討，但文中所給出的湍動能邊界條件為一常數(shù)，這與實測結(jié)果及試驗結(jié)果不符。張建和楊慶山［4］基于標準k－ε湍流模型，將近壁區(qū)第一層網(wǎng)格的湍流量進行平均處理并施加剪切應力，從而使邊界層內(nèi)沿流場不同位置處的流動特征量吻合程度較高。方平治等［5］通過修正的壁面函數(shù)并采用k－ε湍流模型對平衡大氣邊界層進行了研究，研究結(jié)果表明修正的壁面函數(shù)對粗糙度較大的風場效果較好。Yang等［6］從標準k－ε湍流模型方程出發(fā)，假設湍流在整個邊界層內(nèi)處于平衡態(tài)，推導出一組隨高度變化的來流邊界條件，并以風洞試驗數(shù)據(jù)為參考，通過改變湍流模型參數(shù)，并考慮壁面粗糙度，較好地實現(xiàn)了邊界層流場的平衡，并將上述方法應用于SSTk－ω湍流模型［7］，但已有試驗表明，湍流能量僅在邊界層的內(nèi)層范圍內(nèi)才處于平衡［8－9］，并且湍流模型中的參數(shù)也是根據(jù)漸近性原則通過試驗確定［10］。已有研究多數(shù)是針對兩方程湍流模型中的k－ε模型，結(jié)構(gòu)風工程研究中SSTk－ω兩方程湍流模型對鈍體分離流適用性較好，應用范圍較廣，但針對基于SSTk－ω湍流模型的平衡大氣邊界層的研究較為少見。

與已有的平衡大氣邊界層模擬方法不同，本文通過在湍流模型輸送方程中添加源項的方法使來流邊界條件與湍流模型相協(xié)調(diào)，從而實現(xiàn)CFD 分析中空流場大氣邊界層的平衡。研究中首先基于SSTk－ω湍流模型依次確定出平均速度剖面、湍動能剖面，并根據(jù)剪切應力為常量確定出比耗散率剖面，再通過湍動能與比耗散率輸送方程的平衡推導出源項的表達形式，最后通過相應的CFD 數(shù)值算例分別考察了不考慮源項與考慮源項時沿流場不同位置處的平均速度、湍動能及比耗散率剖面的一致性，驗證了本文方法的有效性。

1 SST k－ω 湍流模型

剪切應力輸送（Shear－Stress Transport，簡稱SST）k－ω湍流模型由Menter提出［11］，該模型綜合了k－ω模型在近壁區(qū)計算的優(yōu)點和k－ε模型在遠場計算的優(yōu)點，增加了交叉擴散項，并在湍流粘性系數(shù)的定義中考慮了湍流剪切應力的輸送過程，從而使SSTk－ω湍流模型應用范圍更廣。已有研究表明［12－13］，SSTk－ω湍流模型較傳統(tǒng)的k－ε湍流模型要更適用于具有逆壓梯度流動或分離流動的計算，因而前者更廣泛地應用于大氣邊界層鈍體繞流的計算中。以張量形式表達的SSTk－ω湍流模型的流場輸送方程為：

式中Γ、Γk、Γω為速度u（v或w）、湍動能k及比耗散率ω的有效擴散系數(shù)，其各自定義為：

μt為湍流粘性系數(shù)，其定義為：為k、ω的產(chǎn)生項，Yk、Yω為k、ω的耗散項，Dω為交叉擴散項，各項分別定義為：

式（1）～式（10）中，ρ為空氣密度，p為壓強；μ為分子粘性系數(shù)，S為平均應變率的張量模量；F1、F2為混合函數(shù)；σk、σω、α＊、α∞、β＊及βi為湍流模型中的系數(shù)；σω，2、α1為湍流模型常數(shù)；Si、Sk及Sω為各輸送方程的自定義源項，各項定義詳見文獻［14］。

對于SSTk－ω湍流模型模擬二維、定常、不可壓縮流動的平衡大氣邊界層，一般要滿足以下條件［1］：

（1）豎向速度v為零；

（2）壓強p為常數(shù)；

（3）剪切應力τ為常數(shù)；

（4）湍動能k與比耗散率ω滿足各自的輸送方程。

由以上條件并考慮到μ?μt［1，14］，由此可將流場輸送方程式（1）－（3）化簡可得：

2 邊界條件的確定

2.1 平均風速邊界條件［15］

若將大氣邊界層分為滿足“壁面定律”的近地層及滿足“速度虧損律”的外層兩部分，當存在某一重合區(qū)域使兩種定律都成立時，可得到對數(shù)律的風剖面表達式：

式中u（z）為不同高度處的平均風速，u＊為摩擦速度，Κ為卡門常數(shù)，取Κ＝0.42，z0為粗糙長度。一般認為，對數(shù)律表示大氣底層的風剖面比較理想，當為強風時適用范圍可更高。為考慮地面處z＝0m 時的速度，可將上式改寫為：

2.2 湍流邊界條件

FLUENT［14］為SSTk－ω湍流模型提供了多組湍流邊界條件形式，本文采用經(jīng)常使用的湍動能k與比耗散率ω的一組湍流參數(shù)。對于湍動能k剖面目前沒有統(tǒng)一的表達形式，較多地采用包含湍流強度Iu（z）與速度u（z）形如k＝Ck［Iu（z）·u（z）］2的湍動能剖面形 式（其 中Ck為一常數(shù)）［2，16－18］，部分采用湍動能平衡假設而確定的湍動能剖面［4，6－7］，還有采用沿高 度 為 常數(shù)的湍動能剖面［1，19－20］等?？紤]到湍動能始終大于零，且湍動能的表達式常含u（z）2的因子，并注意到式（15），由此本文將湍動能剖面近似地表示為：

其中Cu1、Cu2為湍流邊界條件參數(shù)，可由試驗數(shù)據(jù)擬合得到，其它參數(shù)意義同前。

對于比耗散率ω的剖面形式目前研究較少。由上文的條件（3）：“剪切應力τ為常數(shù)”，則要求τ滿足：

其中τw為壁面剪切應力，由此可知當且僅當Su＝0時式（11）才 成立。對于式（5）μt表達式中的可視為一系數(shù)，為簡化μt的形式并便于后續(xù)公式的推導，在此先近似地定義湍流粘性系數(shù)為：

將式（15）、式（18）代入式（17）可得：

注意到式（16），由此可得比耗散率ω的剖面形式為：

3 源項的確定

與已有的平衡大氣邊界層模擬方法不同，本文通過在流場輸送方程式（11）～（13）中添加源項的方法使來流邊界條件與湍流模型相協(xié)調(diào)。由于Su＝0，以下著重考慮源項Sk及Sω的表達形式。

3.1 源項Sk 的表達形式

對于輸送方程式（12）左端第一項，由式（16）、式（18）及式（19）并根據(jù)函數(shù)積的求導法則可知：

對于式（12）左端第三項Yk，由式（8）、式（16）及式（20）可知：

其表達式形式與式（21）相似，而10ρβ＊kω的表達式與式（22）相似。由此可知，當項相疊加時，可將合并到這兩項之中，考慮方程式（12）成立，由此可得源項Sk的表達形式為：

式中引入兩個綜合系數(shù)C1k、C2k，其中C1k為考慮系數(shù)β＊／Κ及－10β＊／Κ，C2k為考慮系數(shù)及1／Κ，通過對項的合并處理大大簡化了Sk的形式。

3.2 源項Sω 的表達形式

源項Sω的確定方法與源項Sk的方法類似。對于方程式（13）左端第一項由式（18）、（19）及（20）并根據(jù)函數(shù)積的求導法則可知：

對于式（13）左端第二項Gω，由（7）、式（15）可知：

對于式（13）左端第三項Yω，由式（9）、式（20）可知：

對于式（13）左端第四項Dω，由式（10）、式（16）及式（20），并注意到k與ω僅為z的函數(shù)可知：

由方程式（13）成立，由此可得源項Sω的表達形式為：

同理上式中也引入了四個綜合系數(shù)C1ω、C2ω、C3ω及C4ω，這四個綜合系數(shù)分別與原位置未簡化的系數(shù)項相對應。

由于SSTk－ω湍流模型中各系數(shù)定義的復雜性，加之源項Sk、Sω均是基于簡化的湍流粘性系數(shù)推導而成，因此上述綜合系數(shù)C1k、C2k與C1ω、C2ω、C3ω及C4ω的表達形式可能比式（24）、式（29）中對應的未簡化的系數(shù)項要復雜得多。為弱化湍流模型系數(shù)及對各綜合系數(shù)具體形式的影響，并方便源項的應用，在此假設上述6個綜合系數(shù)均為常數(shù)。需要指出的是，由于SSTk－ω湍流模型的復雜性，即使在假設上述綜合系數(shù)均為常數(shù)的情況下也很難確定各系數(shù)的值，但注意到式（29）中簡化后的源項Sω中第一項為正，第二項與第四項均為負，而第三項的正負號則與參數(shù)Cu1及Cu2的取值有關，可認為第二項與第四項對流場的影響相似，若第三項也為正，同理可認為第一項與第三項對流場的影響也相似；另一方面，考慮到湍流封閉方程的構(gòu)造中通常假設湍動能耗散相關項與湍動能方程中的對應項有相同的機制和公式［10］，對于SSTk－ω湍流模型則有：

其中f（ω）是比耗散率ω的相關項（如產(chǎn)生項、耗散項等；對于k－ε湍流模型，式中比耗散率ω則變?yōu)楹纳⒙师牛?，而f（k）是對應的湍動能相關項，Cω是系數(shù)。若將上述思想引入到源項Sk及Sω的表達式中，則由式（24）可知：式（29）中Sω的表達式中沒有第二項及第三項?；谝陨戏治?，假設式（29）中第二項與第三項對流場的作用可被其它兩項代替，于是可令簡化后的Sω表達式中第二項與第三項的系數(shù)C2ω＝C3ω＝0，由此就減少了未知參數(shù)，同時也簡化了Sω的表達形式。而對于Sk、Sω中剩余的綜合系數(shù)C1k、C2k與C1ω、C4ω，通過初步試算發(fā)現(xiàn)C1k與C1ω對流場水平均勻性的敏感性較強，而C2k與C4ω對流場水平均勻性的敏感性相對較弱，由此可先大致確定出C1k與C1ω的值，然后再確定出C2k與C4ω的值，最后再對各系數(shù)進行局部調(diào)整就可得到基本滿足水平均勻性的平衡大氣邊界層。

4 平衡大氣邊界層數(shù)值模擬驗證

根據(jù)上述分析，以下從數(shù)值模擬的角度來說明在湍流模型輸送方程中添加源項以實現(xiàn)平衡大氣邊界層的有效性。

4.1 數(shù)值計算模型

采用CFD 商業(yè)軟件FLUENT［14］對二維的空流場進行數(shù)值模擬分析，計算區(qū)域尺寸設定為12.0m×1.8m（x×z）（對于三維區(qū)域流動，由于側(cè)面常設置為對稱邊界條件，相當于二維區(qū)域流動），區(qū)域中沒有放置任何物體。采用結(jié)構(gòu)化網(wǎng)格對計算區(qū)域進行離散，最小網(wǎng)格尺寸zmin為0.01m，最終劃分的四邊形網(wǎng)格總數(shù)為8000個，計算區(qū)域網(wǎng)格劃分如圖1所示（縱向為x軸，豎向為z軸）。流動設為二維不可縮的定常流動，計算模型采用SSTk－ω湍流模型，計算區(qū)域各邊界條件的定義如表1所示，其中入口處邊界條件表達式中的各參數(shù)均通過風洞試驗數(shù)據(jù)擬合得到［6］。此外，為說明源項Sk、Sω對實現(xiàn)平衡大氣邊界層的有效性，在流場的k、ω輸送方程中分不考慮源項和考慮源項兩種情況，不考慮源項時各綜合系數(shù)均取零值；考慮源項時，源項中各綜合系數(shù)值根據(jù)上文方法最終確定情況如表2所示。

表1 計算區(qū)域邊界條件Table 1 Boundary conditions of the computational domain

表2 輸送方程源項Table 2 Source terms of transport equations

對于求解控制參數(shù)的設置，壓力與速度耦合采用SIMPLEC算法，動量方程、湍動能及比耗散率輸送方程均采用二階離散格式，使用入口處來流邊界條件各特征量對流場進行初始化，殘差收斂精度均設置為1.0×10－6。

圖1 計算區(qū)域網(wǎng)格劃分示意Fig.1 Schematic diagram of mesh arrangement

4.2 數(shù)值模擬結(jié)果

當不考慮源項Sk、Sω時，流場收斂后提取沿流場縱向x＝－6m（入口處）、x＝－2m（靠近入口的1／3 總長處，也是模型常放置的地方）、x＝2m 及x＝6m（出口處）處的平均速度u、湍動能k及比耗散率ω的剖面如圖2所示，由圖可知，各處的速度剖面及比耗散率剖面總體上能保持一致，但入口處的湍動能剖面與其它位置的湍動能剖面差異明顯，表明此時的流場尚未達到平衡。而考慮源項時對應位置的各特征量剖面如圖3所示，由圖可知，各處的速度剖面、湍動能剖面及比耗散率剖面均有較高程度的一致性，表明此時的流場基本達到了平衡，較好地實現(xiàn)了平衡大氣邊界層。此外，對比圖2（a）與圖3（a）中的速度剖面可知，圖2（a）頂端處速度剖面的一致性較圖3（a）中的要稍差一些，同樣對于圖2（c）中比耗散率剖面在高度較低時其剖面的一致性也稍差于圖3（c）的。通過上述無源項與有源項時流場特征量剖面的對比表明，在湍流模型輸送方程中添加源項能較有效地實現(xiàn)大氣邊界層的平衡。

圖2 不考慮源項時不同位置的平均速度剖面及湍流特征剖面Fig.2 Profiles of average wind and turbulence wind at different positions without considering the source terms

圖3 考慮源項時不同位置的平均速度剖面及湍流特征剖面Fig.3 Profiles of average wind and turbulence wind at different positions with considering the source terms

5 結(jié)論

實現(xiàn)平衡大氣邊界層是進行數(shù)值模擬的前提，滿足平衡大氣邊界層的流場能顯著提高數(shù)值計算結(jié)果的精度［7］。本文基于SSTk－ω湍流模型，提出通過在湍動能與比耗散率輸送方程中添加源項的方法使來流邊界條件與湍流模型相協(xié)調(diào)，從而實現(xiàn)CFD 分析中空流場大氣邊界層的平衡。CFD 數(shù)值算例表明在湍流模型輸送方程中考慮源項能較好地實現(xiàn)大氣邊界層的平衡，驗證了本文方法的有效性。本文提出的方法雖然僅針對SSTk－ω湍流模型，但從方法實現(xiàn)的過程來看，同樣適用于標準k－ω及k－ε系列等兩方程湍流模型，具有一定的普遍性。本研究為平衡大氣邊界層的研究提供了新的思路和研究方法，對結(jié)構(gòu)風工程中CFD 數(shù)值模擬有一定的理論及應用價值。

［1］RICHARDS P J，HOXEY R P.Appropriate boundary conditions for computational wind engineering models using thek－εmodel［J］.JournalofWindEngineering andIndustrialAerodynamics，1993，46－47：145－153.

［2］BLOCKEN B，CARMELIET J，STATHOPOULOS T.CFD evaluation of wind speed conditions in passages between parallel buildings－effect of wall－function roughness modifications for the atmospheric boundary layer flow［J］.JournalofWindEngineeringandIndustrial Aerodynamics，2007，95（9－11）：941－962.

［3］BLOCKEN B，STATHOPOULOS T，CARMELIET J.CFD simulation of the atmospheric boundary layer：wall function problems［J］.AtmosphericEnvironment，2007，41（2）：238－252.

［4］張建，楊慶山.基于標準k－ε模型的平衡大氣邊界層模擬［J］.空氣動力學學報，2009，27（6）：729－735.

［5］方平治，顧明，談建國.計算風工程中基于k－ε系列湍流模型的數(shù)值風場［J］.水動力學研究與進展，2010，25（4）：475－483.

［6］YANG Y，GU M，CHEN S Q，et al.New inflow boundary conditions for modelling the neutral equilibrium atmospheric boundary layer in computational wind engineering［J］.JournalofWindEngineeringandIndustrial Aerodynamics，2009，97（2）：88－95.

［7］YANG W，QUAN Y，JIN X Y，et al.Influences of equilibrium atmosphere boundary layer and turbulence parameter on wind loads of low－rise building［J］.JournalofWindEngineeringandIndustrialAerodynamics，2008，96（10－11）：2080－2092.

［8］陳懋章.粘性流體動力學基礎［M］.北京：高等教育出版社，2002.

［9］是勛剛.湍流［M］.天津：天津大學出版社，1994.

［10］張兆順.湍流［M］.北京：國防工業(yè)出版社，2002.

［11］MENTER F R.Two－equation eddy－viscosity turbulence models for engineering applications［J］.AIAAJournal，1994，32（8）：1598－1605.

［12］LOUREIRO J B R，ALHO A T P，SILVA FREIRE A P.The numerical computation of near－wall turbulent flow over a steep hill［J］.JournalofWindEngineering andIndustrialAerodynamics，2008，96（5）：540－561.

［13］EI－BEHERY S M，HAMED M H.A comparative study of turbulence models performance for separating flow in a planar asymmetric diffuser［J］.Computers＆ Fluids，2011，44（1）：248－257.

［14］FLUENT.Fluent 6.3user's guide［M］.Labnon，New Hampshire.Fluent Inc.，2006.

［15］SIMIU E，SCANLAN R H.Wind effects on structures：fundamentals and applications to design［M］.USA：Wiley，1996.

［16］LAKEHAL D.Application of thek－εmodel to flow over a building placed in different roughness sublayers［J］.JournalofWindEngineeringandIndustrialAerodynamics，1998，73（1）：59－77.

［17］SAGRADO A P G，BEECK J V，RAMBAUD P，et al.Numerical and experimental modelling of pollutant dispersion in a street canyon［J］.JournalofWindEngineeringandIndustrialAerodynamics，2002，90（4－5）：321－339.

［18］GAO Y，CHOW W K.Numerical studies on air flow around a cube［J］.JournalofWindEngineeringand IndustrialAerodynamics，1995，93（2）：115－135.

［19］PASLEY H，CLARK C.Computational fluid dynamics study of flow around floating－roof oil storage tanks［J］.JournalofWindEngineeringandIndustrialAerodynamics，2000，86（1）：37－54.

［20］HARGREAVES D M，WRIGHT N G.On the use of thek－εmodel in commercial CFD software to model the neutral atmospheric boundary layer［J］.Journalof WindEngineeringandIndustrialAerodynamics，2007，95（5）：355－369.