基于改進剪切應(yīng)力輸運k-ω湍流模型對考慮風向隨高度偏轉(zhuǎn)大氣邊界層的數(shù)值模擬

馮成棟，顧 明

（同濟大學(xué)土木工程防災(zāi)國家重點實驗室，上海 200092）

當前，計算風工程（computational wind engineering，CWE）在結(jié)構(gòu)風工程領(lǐng)域的研究與應(yīng)用越來越廣泛。在計算風工程中，正確模擬不可壓縮中性層結(jié)水平均勻穩(wěn)態(tài)正壓的大氣邊界層（以下簡稱“大氣邊界層”）是非常重要的一項工作，它直接關(guān)系到后續(xù)建筑繞流模擬的準確性。所謂“水平均勻”，是指各流動物理量的水平梯度為零，大氣邊界層風特性在入口到模型之間不發(fā)生明顯的變化，也即大氣邊界層的“自保持”。為了滿足這一要求，諸多學(xué)者從不同方面進行了研究。Richards 和Hoxey［1］基于標準k-ε模型，提出了一組適用于近地層的邊界條件。這組邊界條件已得到廣泛應(yīng)用［2-3］。Blocken等［4］、Hargreaves和Wright［5］、方平治等［6］據(jù)此對壁面函數(shù)以及粗糙度參數(shù)進行修正，以實現(xiàn)近地面風場的準確模擬和自保持。Yang等［7-8］推導(dǎo)了近似滿足大氣邊界層自保持要求的湍動能表達式，建議了一組新的入流邊界條件，將其成功應(yīng)用于標準k-ε模型和剪切應(yīng)力輸運（shearstress transport，SST）k-ω模型。唐煜等［9］基于SSTk-ω模型和平衡湍流假設(shè)，推導(dǎo)了適用于一般風工程計算的入口來流邊界條件表達式，并給出了適用于地形尺度風場計算的湍流模型常數(shù)建議值。Parente等［10-11］、Longo 等［12］則對k-ε兩方程分別添加了源項，以實現(xiàn)入流邊界條件和湍流模型方程的相容。Hargreaves 和Wright［5］、O’Sullivan 等［13］、Richards 和Norris［14］特別強調(diào)了流域頂部邊界條件對于實現(xiàn)自保持的重要性。

然而，以上研究僅著眼于近地層特性的模擬。實際上，大氣邊界層由近地層和Ekman層組成，其中的氣流可近似看作是水平氣壓梯度力、科里奧利力和湍流摩擦力平衡的結(jié)果?！叭ζ胶狻辈粌H導(dǎo)致了平均風速隨高度的增加而增加，也導(dǎo)致了平均風向隨高度的增加而沿順時針方向偏轉(zhuǎn)（北半球）。圖1給出了大氣邊界層中風向隨高度的變化示意圖（θ表示風偏角，Vgr表示梯度風高度處的風速），稱為“Ekman 螺線”［15］。一般認為，近地層高度數(shù)值在20～200 m 左右［16］，此范圍內(nèi)風向基本不隨高度變化，而在約占大氣邊界層總高度90%的Ekman 層中，風向隨高度的偏轉(zhuǎn)表現(xiàn)得尤為顯著［17-18］。

近年來，全球出現(xiàn)了多座超過500 m 高度的超高層建筑，如上海中心大廈（Shanghai Tower，632 m，2015）、哈利法塔（Burj Khalifa Tower，828 m，2010）等，而在規(guī)劃中的國王塔（Kingdom Tower）更是超過1 000 m。對這些千米級超高層建筑而言，高空風特性將會對風荷載和風致響應(yīng)產(chǎn)生很大影響，甚至可能是決定性的。我國現(xiàn)行建筑結(jié)構(gòu)荷載規(guī)范［19］認為各類地貌梯度風高度以上平均風速大小保持不變，但是近來有研究表明，在千米量級附近，平均風速隨高度單調(diào)增加［20］。另外，荷載規(guī)范也沒有考慮風向隨高度偏轉(zhuǎn)。風向隨高度的偏轉(zhuǎn)將使得千米級建筑的風力沿高度分布和響應(yīng)更加復(fù)雜［21］，在結(jié)構(gòu)設(shè)計時應(yīng)該引起高度重視。從基本方程出發(fā)，對考慮風向隨高度偏轉(zhuǎn)的大氣邊界層進行計算流體動力學(xué)（computational fluid dynamics，CFD）數(shù)值模擬和自保持實現(xiàn)方法的探討，是研究千米級建筑風荷載和風致響應(yīng)特性的基礎(chǔ)。

SSTk-ω模型能夠有效計算湍流切應(yīng)力在逆壓梯度邊界層的輸運，被視為目前RANS（Reynoldsaveraged Navier-Stokes）方法中最適用于鈍體分離流模擬的兩方程湍流模型［22］。為了與后續(xù)建筑繞流模擬工作銜接，本文基于SSTk-ω模型，通過修改模型參數(shù)、方程源項和湍動粘度，提出了“改進SSTk-ω模型”，并采用該模型對考慮風向隨高度偏轉(zhuǎn)的大氣邊界層進行了CFD數(shù)值模擬。通過預(yù)前模擬和主模擬兩個步驟，對風場自保持的實現(xiàn)方法進行了探討。通過對比CFD模擬結(jié)果和實測結(jié)果，驗證了改進SSTk-ω模型模擬結(jié)果的合理性；通過出流面和入流面各物理量的對比，驗證了自保持方法的有效性。

1 預(yù)前模擬

所謂預(yù)前模擬，是指在模擬主要研究對象之前，設(shè)立單獨的計算域，對入口湍流進行模擬，而后將模擬得到的入口湍流施加于主模擬計算域的方法。由于本文研究的大氣邊界層滿足穩(wěn)態(tài)假定和水平均勻性假定，因此并不涉及時間存儲的問題，且預(yù)前模擬只需滿足豎向的網(wǎng)格劃分與隨后主模擬所用網(wǎng)格的豎向劃分一致即可，在水平方向可以采用數(shù)目較少的網(wǎng)格。這種方法可在模擬精度達到一定要求的前提下，節(jié)省計算資源，因此在氣象學(xué)的相關(guān)模擬、地形風場模擬和建筑繞流的前期風場模擬中被廣泛采用［2，13，23］。

基于SSTk-ω模型，通過修改模型參數(shù)、方程源項和湍動粘度，對考慮風向隨高度偏轉(zhuǎn)的大氣邊界層進行預(yù)前模擬，為后續(xù)主模擬提供入口信息。

1.1 基本方程

引入不可壓縮中性層結(jié)水平均勻穩(wěn)態(tài)正壓大氣邊界層假定，忽略分子粘度，并假定垂直方向運動與水平方向運動相比可以忽略，則連續(xù)性方程自動滿足，水平方向的動量方程可以簡化為

式中：（u，v）為水平平均風速分量；p為壓力；ρ為空氣密度；f為科里奧利參數(shù)（f=2ωsinφ，φ為緯度，ω為地球自轉(zhuǎn)角速度）；′和′表示雷諾應(yīng)力（湍流切應(yīng)力）。式（1）和式（2）等號左端從左到右依次為水平氣壓梯度力、科里奧利力和湍流摩擦力。上述方程稱為經(jīng)典Ekman 模型（也稱為“三力平衡”模型）。

為了使動量方程閉合，引入Boussinesq 假定，將雷諾應(yīng)力與平均速度梯度聯(lián)系起來：

式中：μt為湍動粘度；k為湍動能。

馮成棟和顧明［24］已通過修改標準k-ε模型的基本參數(shù)，對考慮風向偏轉(zhuǎn)的大氣邊界層進行了模擬。為了與后續(xù)建筑繞流模擬工作銜接，本文將大氣邊界層模擬方法推廣至SSTk-ω模型，即從k-ε模型方程出發(fā)，根據(jù)相關(guān)數(shù)學(xué)關(guān)系，推導(dǎo)得到適用于SSTkω模型的、可實現(xiàn)大氣邊界層內(nèi)風向偏轉(zhuǎn)的改進方案。推導(dǎo)過程如下。

引入比耗散率ω=ε/（Cμk）（其中ε為湍動能耗散率，Cμ為模型常數(shù)），將其代入簡化的標準k-ε模型方程中（參照文獻［24］中公式（4）—（6）），可得（忽略源項Sk和Sε）：

式中：Gk為湍動能生成項；C1ε、C2ε、σk和σε均為模型常數(shù)。

化簡式（6），并用式（6）-式（5）乘以ω可得：

式（7）兩端同除以k，可得：

結(jié)合比耗散率定義和式（4），進一步簡化式（8）得：

若令σk=σε=σω，則式（9）可進一步簡化如下：

令C1ε-1=C1ω，C2ε-1=C2ω，為了 將Apsley 和Castro［25］提出的方法推廣至SSTk-ω模型，參照文獻［24］中的公式（7），將式（10）等號左端的C1ω（=C1ε-1）用C*1ω替代，如式（11）：

式中：lmax表示邊界層最大混合長度；lm為局地混合長度（即湍流長度尺度），用式（12）計算：

綜上，為了使大氣邊界層模擬結(jié)果更加接近實測，基于SSTk-ω模型的改進方案如下（簡化起見，取σk=σε=σω）：

式（13）—（15）表達的模型在本文中被稱為“改進SSTk-ω模型”。需要注意的是，該模型方程與ANSYS Fluent 15.0 內(nèi)置的SSTk-ω模型方程有所區(qū)別，為此，應(yīng)當參照ANSYS Fluent Theory Guide［26］中SSTk-ω湍流模型方程，通過UDF（userdefined function）中的源項宏和湍動粘度宏，對內(nèi)置的SSTk-ω模型方程進行修改，使之與式（13）—（15）吻合。

1.2 模型參數(shù)

湍流模型常數(shù)的取值，隨研究問題性質(zhì)的不同而異。諸多研究表明，Cμ在大氣邊界層的模擬中，取值當與默認推薦值不同，常取Cμ≈0.03［7，23］。需要注意的是，Cμ在SSTk-ω模型中對應(yīng)β*∞，故其值也應(yīng)取為0.03。C1ε和C2ε的取值參考文獻［23］，分別取為1.52 和1.833，故C1ω和C2ω應(yīng) 分 別 取 為0.52 和0.833。關(guān)于σε，應(yīng)由以下約束關(guān)系確定，以使在近地層中能夠再現(xiàn)對數(shù)律［23］：

式中：κ為馮·卡門常數(shù)κ，可取0.42。σk和σω的取值則沒有嚴格規(guī)定，簡化起見，可取σk=σω=σε。另外，對比式（15）與ANSYS Fluent Theory Guide［26］中SSTk-ω模型方程可以發(fā)現(xiàn)，C2ω×Cμ對應(yīng)SSTk-ω模型常數(shù)βi。

除了本身的模型常數(shù)外，式（11）還引入了參數(shù)lmax。Koblitz 等［23］指出，lmax取28 m（實尺）時，模擬結(jié)果和實測結(jié)果［27］吻合得很完美。本文采用此值。

表1 給出了ANSYS Fluent 15.0 內(nèi)置SSTk-ω模型常數(shù)的推薦取值和本文大氣邊界層模擬所用改進SSTk-ω模型常數(shù)取值（僅列出有改動的模型常數(shù)）。

表1 模型常數(shù)取值Tab.1 Model constant values

1.3 目標風場和相似準則

常規(guī)大氣邊界層風洞的風場模擬中，往往只關(guān)注平均風剖面、湍流度剖面、脈動風速功率譜和湍流積分尺度等。然而，在真實的大氣環(huán)境中，風場是水平氣壓梯度力、科里奧利力、湍流摩擦力等相互作用的結(jié)果，Ekman 層內(nèi)的風速矢量呈現(xiàn)出隨高度偏轉(zhuǎn)的普遍特征［28-29］。本文擬采用著名的Leipzig實測風剖面作為大氣邊界層數(shù)值模擬的目標風場。Leipzig風剖面于1931 年由探空測風氣球測得，并經(jīng)Lettau于1950 年重新復(fù)核［27］。表2 給出了Leipzig 實測風剖面的風場參數(shù)，G為地轉(zhuǎn)風速，φ為緯度，u*為摩擦速度，z0為空氣動力學(xué)粗糙長度，f為科里奧利參數(shù)，θ0表示地面附近水平風速矢量到水平地轉(zhuǎn)風速矢量的總風偏角。關(guān)于目標風場的湍動能剖面，有大量現(xiàn)場實測和數(shù)值模擬資料可供參考［30-32］。

建筑結(jié)構(gòu)的風洞試驗和數(shù)值模擬通常采用縮尺模型。為了與后續(xù)建筑繞流模擬工作銜接，可根據(jù)常規(guī)風洞試驗的風場參數(shù)和相似準則，確定各物理量的縮尺比。本文在數(shù)值模擬中，采用的幾何縮尺比λL為1：1333、風速縮尺比λV為1：1.12，密度縮尺比λρ為1：1，在此縮尺比下得到的Leipzig 風剖面參數(shù)u*和z0與TJ-2 某風洞試驗［33］目標風剖面取值相同。此外，由于動量方程式（1）和式（2）引入了科氏力，故必須考慮動力氣象學(xué)中羅斯貝數(shù)Ro 相似，羅斯貝數(shù)Ro定義如下：

式中：V表示特征速度；L表示特征長度；f為科里奧利參數(shù)。由羅斯貝數(shù)相似準則可以導(dǎo)出科里奧利參數(shù)f的縮尺比：

縮尺后用于數(shù)值模擬的Leipzig 風剖面的風場參數(shù)一并示于表2。

表2 Leipzig風剖面參數(shù)Tab.2 Leipzig wind profile parameters

1.4 網(wǎng)格劃分

為了驗證本文大氣邊界層預(yù)前模擬結(jié)果的網(wǎng)格無關(guān)性，采用表3中的三組網(wǎng)格方案，水平向網(wǎng)格均勻劃分，高度方向網(wǎng)格采用增長比率，由密到疏。

表3 預(yù)前模擬所用網(wǎng)格Tab.3 Grids for precursor simulations

1.5 邊界條件和求解設(shè)置

大氣邊界層預(yù)前模擬采用的邊界條件設(shè)置見圖2和表4。

表4 預(yù)前模擬邊界條件Tab.4 Boundary conditions for precursor simulations

圖2 預(yù)前模擬邊界條件Fig.2 Boundary conditions for precursor simulations

假定地轉(zhuǎn)風方向平行于x軸，為了驅(qū)動空氣形成風場，需要施加壓力梯度。考慮到在正壓大氣邊界層中，水平氣壓梯度力不隨高度變化，則其值可由地轉(zhuǎn)平衡關(guān)系推導(dǎo)得到：

式中：ρ為空氣密度；f為科里奧利參數(shù)；（ug，vg）為地轉(zhuǎn)風速分量。

為了使風向隨高度形成偏轉(zhuǎn)，需要在動量方程中以源項形式加入科氏力，可通過UDF（user-defined function）源項宏實現(xiàn)。為了使模擬結(jié)果更加接近實測，“改進SSTk-ω模型”需借助UDF源項宏和湍動粘度宏實現(xiàn)。

本文數(shù)值模擬采用ANSYS Fluent 15.0模擬平臺，穩(wěn)態(tài)求解，設(shè)置壓力速度耦合方式為SIMPLEC，動量方程和湍流模型方程非線性對流項采用二階迎風格式離散，壓力插值格式采用Standard，梯度插值方法采用Least Squares Cell Based。所有變量和連續(xù)性方程的殘差收斂標準設(shè)置為10-6。

1.6 模擬結(jié)果分析

圖3～圖7 分別給出了預(yù)前模擬所得平均風速u、平均風速v、平均合風速Uavg、平均風偏角θ隨高度的變化曲線，以及量綱一化湍動能隨量綱一化高度（大氣邊界層高度h取Detering 和Etling［34］給出的1600 m，縮尺后為1.2 m）的變化曲線，并將Leipzig實測結(jié)果［27］、Ekman 理論解［28］、同濟大學(xué)TJ-2 某風洞試驗?zāi)繕孙L剖面［33］及其對數(shù)律擬合曲線、ESDU（engineering sciences data unit，英國工程科學(xué)數(shù)據(jù)庫）規(guī)范平均風剖面［35］，以及由Detering 和Etling［34］推演得到Leipzig 實測中的湍動能、Brost 等［30］和Grant［31］的湍動能實測結(jié)果、Esau［32］針對不同大氣邊界層高度h的LES模擬得到的湍動能結(jié)果一并示于圖中，以資比較。

由圖3～圖7可見，無論是平均風速、風偏角，還是湍動能，模擬結(jié)果均與實測結(jié)果吻合得很好。由于經(jīng)典Ekman 解的前提是湍動粘性系數(shù)為常量（縮尺后此處取為0.0053 m2·s-1），所以其值與模擬結(jié)果有一定的偏差，但是趨勢一致。對于量綱一化湍動能，本文模擬結(jié)果與實測結(jié)果較為吻合，且基本落在LES 模擬結(jié)果范圍內(nèi)?？傮w來講，本文大氣邊界層模擬結(jié)果合理，可以滿足后續(xù)自保持研究和建筑繞流模擬的需要。

圖3 平均風速u隨高度的變化Fig.3 Profiles of mean velocity component u with height

圖7 量綱一化湍動能隨量綱一化高度的變化Fig.7 Profiles of non-dimensional turbulent kinetic energy with non-dimensional height

值得注意的是，對數(shù)律平均風剖面對真實大氣邊界層風速存在嚴重低估。具體說來，圖5 中的“TJ-2 目標風剖面對數(shù)律擬合”得到的u*和z0與表2縮尺之后的Leipzig 風剖面對應(yīng)參數(shù)（“模擬”欄）相同，但由圖5可見，對數(shù)律平均風剖面僅在近地層范圍內(nèi)有效，其適用高度為100 m左右，這與Li等［36］和Drew 等［37］通過分析實測結(jié)果得到的對數(shù)律在高空風的預(yù)測中對實際風速低估的結(jié)論相吻合。此外，ESDU規(guī)范平均風剖面［35］與實測平均合風速結(jié)果總體上較為接近（在較高位置處存在低估，但相比對數(shù)律風剖面已有明顯改善），然而尚不能準確反映風向沿高度偏轉(zhuǎn)的規(guī)律。

圖5 平均合風速隨高度的變化Fig.5 Profiles of mean velocity magnitude with height

2 主模擬

大氣邊界層風場自保持的實現(xiàn)至少需要處理以下三方面的相容問題：①來流邊界條件和湍流模型的相容問題；②壁面函數(shù)和湍流模型的相容問題；③來流邊界條件和壁面函數(shù)的相容問題［6］。預(yù)前模擬得到的大氣邊界層風場可以很好地與壁面函數(shù)和湍流模型相適應(yīng)［13］，而壁面函數(shù)和湍流模型的相容問題已由CFD計算軟件（如ANSYS Fluent 15.0）內(nèi)置解決。因此，采用預(yù)前模擬和主模擬兩個步驟可使以上三方面問題基本解決，從而較好地實現(xiàn)大氣邊界層風場的自保持。

圖4 平均風速v隨高度的變化Fig.4 Profiles of mean velocity component v with height

圖6 平均風偏角隨高度的變化Fig.6 Profiles of mean wind veering angle with height

2.1 網(wǎng)格劃分

為了驗證本文大氣邊界層主模擬結(jié)果的網(wǎng)格無關(guān)性，采用表5中的三組網(wǎng)格方案，水平向網(wǎng)格均勻劃分，高度方向網(wǎng)格采用增長比率，由密到疏。

表5 主模擬所用網(wǎng)格Tab.5 Grids for main simulations

2.2 邊界條件和求解設(shè)置

在主模擬中，采用UDF 將1.6 小節(jié)中改進SSTk-ω模型預(yù)前模擬所得平均風速剖面、湍動能剖面和比耗散率剖面以入流邊界條件的形式施加于主模擬區(qū)域入口。需要特別說明的是，由1.6小節(jié)圖3～圖5 可見，在相當一部分高度范圍內(nèi)，平均風速大小超越了地轉(zhuǎn)風速G（表2），且在平均風向與地轉(zhuǎn)風向第一次平行后，平均風速的v分量出現(xiàn)了負值，這一現(xiàn)象已為諸多實測和理論分析所證實［16］。這種“超地轉(zhuǎn)”現(xiàn)象導(dǎo)致了主模擬區(qū)域邊界條件指定的復(fù)雜性，特別是給側(cè)向邊界的指定增加了難度；側(cè)向邊界的回流對應(yīng)于平均風速v分量的負值，這意味著側(cè)向邊界的入流和出流情況很難明確地予以指定［38］。為了克服這一困難，可以將1.6 小節(jié)模擬所得來流風場整體旋轉(zhuǎn)一定的角度，而流域外輪廓保持不變。方便起見，可將來流風場整體順時針旋轉(zhuǎn)26.1o（此角度為來流地面風與地轉(zhuǎn)風的夾角，參見表2），即將“原地面風”和“原地轉(zhuǎn)風”旋轉(zhuǎn)至“新地面風”和“新地轉(zhuǎn)風”的位置，如圖8 所示。在新的來流風場下，流域邊界的鑒定比較明確，在出流邊界不會出現(xiàn)回流現(xiàn)象，因而更加合理。圖9和表6給出了大氣邊界層主模擬邊界條件設(shè)置。

表6 主模擬邊界條件.Tab.6 Boundary conditions for main simulations

圖8 來流風場旋轉(zhuǎn)示意圖Fig.8 Overall inflow rotation diagram

此外，由于以速度入口邊界給定的風場風向已隨高度偏轉(zhuǎn)，且速度入口邊界可作為“風場驅(qū)動源”［38］，故可在主模擬中去掉動量方程中的科氏力源項和水平壓力梯度驅(qū)動項。為了實現(xiàn)風場自保持，預(yù)前模擬中改進SSTk-ω模型添加的源項宏和湍動粘度宏仍需在主模擬中添加。其余參數(shù)設(shè)置和求解設(shè)置同前。

2.3 模擬結(jié)果分析

圖10～圖12分別給出了入流面（inlet）和出流面（outlet）的平均合風速Uavg、平均風偏角θ和湍動能k隨高度的變化曲線。由圖可見，雖然在地面附近出流面的各物理量相比入流面有所變化，但是總體而 言，大氣邊界層風場在流域中得到了較好的保持。

圖10 出流面與入流面的平均合風速對比Fig.10 Comparison of mean velocity magnitude of inlet and outlet

圖12 出流面與入流面的湍動能對比Fig.12 Comparison of turbulent kinetic energy of inlet and outlet

3 結(jié)論

本文基于SSTk-ω模型，通過修改模型參數(shù)、方程源項和湍動粘度，提出了“改進SSTk-ω模型”，并采用該模型對考慮風向隨高度偏轉(zhuǎn)的大氣邊界層進行了CFD 數(shù)值模擬。通過預(yù)前模擬和主模擬兩個步驟，對風場自保持的實現(xiàn)方法進行了探討。得到如下結(jié)論：

（1）本文提出的“改進SSTk-ω模型”引入了比耗散率生成項和湍流長度尺度之間的反饋機制，綜合考慮了模型參數(shù)取值的合理性，可取得與大氣邊界層實測較一致的模擬結(jié)果。對數(shù)律平均風剖面僅適用于近地層，對高空風速存在嚴重低估的現(xiàn)象。ESDU規(guī)范風剖面對整個大氣邊界層平均風速預(yù)測較好，但尚不能準確反映風向沿高度偏轉(zhuǎn)的規(guī)律。

圖11 出流面與入流面的平均風偏角對比Fig.11 Comparison of mean wind veering angle of inlet and outlet

（2）通過預(yù)前模擬和主模擬兩個步驟，即把“改進SSTk-ω模型”預(yù)前模擬得到的風剖面以入流邊界條件的形式施加于主模擬入口，平均風速、平均風偏角和湍動能均可以在流域中得到較好的保持。

（3）本文提出的“改進SSTk-ω模型”，是對已有修正k-ε模型的推廣，豐富和完善了考慮風向偏轉(zhuǎn)的大氣邊界層RANS模擬方法?；诖诉M行的大氣邊界層風場模擬及自保持方法研究，可為后續(xù)建筑繞流模擬提供前提和基礎(chǔ)。由于SSTk-ω模型在鈍體分離流模擬中更具優(yōu)勢，因而“改進SSTk-ω模型”可在實現(xiàn)偏轉(zhuǎn)風場準確模擬的基礎(chǔ)上，對建筑平均風壓和繞流進行較為精準的預(yù)測。