李樹民，孫 壯，李澤祥，李鑫洋，羅朝勻，李 輝，孫鑫宇，劉偉毅

(成都流體動力創(chuàng)新中心，成都 610072)

0 引言

風(fēng)能資源作為最具開發(fā)潛力的可再生能源之一，其精細(xì)化評估直接關(guān)系到風(fēng)能利用效率和極端風(fēng)速下的風(fēng)電機(jī)組安全性。近年來，風(fēng)電機(jī)組不斷向大型化方向發(fā)展，葉尖高度超過200 m、風(fēng)輪直徑超過150 m 的大型風(fēng)電機(jī)組屢見不鮮，這種發(fā)展趨勢給風(fēng)能資源評估帶來了新的挑戰(zhàn)。現(xiàn)有的風(fēng)資源評估技術(shù)建立在經(jīng)典的近地層相似理論基礎(chǔ)上，僅適用于高度100 m 左右的近地層[1]。在100～300 m 高度的大氣邊界層（atmospheric boundary layer,ABL）下層不再是地表摩擦力和氣壓梯度力二力平衡的大氣運(yùn)動，而是加入了地轉(zhuǎn)偏向力的三力平衡情況。因此重新認(rèn)識風(fēng)特性，拓展或重構(gòu)風(fēng)資源理論和方法，揭示典型地形影響下300 m 高度內(nèi)的風(fēng)和湍流特性及其形成機(jī)理，已成為國內(nèi)外風(fēng)電領(lǐng)域的新的研究重點(diǎn)。計(jì)算流體力學(xué)（CFD）方法可以獲得精細(xì)化的流場細(xì)微結(jié)構(gòu)，已成功應(yīng)用于很多工業(yè)領(lǐng)域，在精細(xì)化風(fēng)資源評估中具有一定潛力。由于目前傳統(tǒng)工業(yè)用CFD 軟件尚不具備大氣效應(yīng)的求解能力，而主流的風(fēng)資源評估軟件均被國外壟斷，使得我國風(fēng)資源精細(xì)化評估工作受到嚴(yán)重制約，因此建立一套能夠適用于我國典型地形風(fēng)場特征的大氣模式的仿真模擬方法，使其具備中性、非中性（穩(wěn)定、不穩(wěn)定等）、地轉(zhuǎn)效應(yīng)等大氣風(fēng)場仿真能力，將中尺度大氣模式模擬方法移植到傳統(tǒng)CFD 中（本文中稱之為大氣CFD 模式），實(shí)現(xiàn)天氣背景風(fēng)場的大氣模式與傳統(tǒng)CFD 模擬的有效銜接，將對我國風(fēng)資源精細(xì)化評估、微觀選址、風(fēng)電預(yù)測起到重要支撐作用。

風(fēng)資源評估經(jīng)歷了從無熱力變化的中性大氣到考慮熱力變化的非中性大氣的研究歷程。在中性大氣方面，熱力作用不明顯，風(fēng)速較為強(qiáng)勁，研究者多集中在計(jì)算風(fēng)工程領(lǐng)域。Richards等[2]最早提出了中性大氣邊界層流場模擬方法及邊界條件，并在計(jì)算風(fēng)工程領(lǐng)域得到廣泛應(yīng)用。Hargreaves等[3]通過改變壁面邊界條件，并在計(jì)算域上部考慮大氣的剪切作用，實(shí)現(xiàn)了中性風(fēng)速廓線在較長計(jì)算域中的保持性研究。楊易等[4]從邊界條件應(yīng)滿足湍流模型的角度出發(fā)，引入了新的入流邊界條件，在原有速度入口廓線的基礎(chǔ)上，給出了入口湍動能邊界條件隨高度變化關(guān)系式。張曉東[5]針對壁面函數(shù)模型與入流剖面不一致會導(dǎo)致流向上出現(xiàn)不合理的梯度變化情況，提出了一種壁面函數(shù)模型，緩解了入口速度廓線在計(jì)算域中的梯度變化，并應(yīng)用于Bolund Hill 和Askervein Hill兩個(gè)典型山地風(fēng)場的仿真模擬與驗(yàn)證。遆子龍等[6]基于管道流動的粗糙度理論研究了地表粗糙度對山區(qū)峽谷橋梁區(qū)域風(fēng)場的影響。

在非中性大氣風(fēng)場模擬方面，研究主要集中在風(fēng)能領(lǐng)域，以丹麥科技大學(xué)（DTU）的研究為代表[7]，主要根據(jù)莫寧-奧布霍夫相似理論，將莫寧-奧布霍夫長度、熱浮力、科氏力等大氣特征參數(shù)引入湍流模型，實(shí)現(xiàn)了非中性大氣的仿真模擬，在這一思想指導(dǎo)下開發(fā)出了EllipSys3D軟件，并衍生出了風(fēng)資源評估軟件WAsP[8]。法國美迪公司的Meteodyn WT 風(fēng)資源評估軟件采用了同樣的思想，基于混合長度理論將莫寧-奧布霍夫長度引入一方程湍流模型中，并通過入口速度和溫度廓線的分層理論，實(shí)現(xiàn)了非中性大氣的仿真[9]。在復(fù)雜地形非中性大氣仿真方面，比較有代表性的是采用域前模擬方法來開展[7]：第一步，在地形計(jì)算域前建立一個(gè)空的計(jì)算域，利用循環(huán)計(jì)算的方法，配合非中性的湍流模型和邊界條件，發(fā)展出平坦地形條件下的非中性風(fēng)場狀態(tài)；第二步，以第一步發(fā)展出的風(fēng)速為邊界條件，施加進(jìn)帶有復(fù)雜地形的計(jì)算域中，實(shí)現(xiàn)非中性風(fēng)場的仿真模擬。雖然經(jīng)過了多年的發(fā)展，但采用CFD 方法模擬大氣風(fēng)場依然是個(gè)非常有挑戰(zhàn)性的難題[10-11]。

由于國外商業(yè)軟件的通用性、便捷性和強(qiáng)大的市場運(yùn)作能力，嚴(yán)重壓縮了我國自主軟件的生存空間。目前在風(fēng)電領(lǐng)域，國內(nèi)已有的風(fēng)電場數(shù)值模擬軟件主要依靠國外。國產(chǎn)自主軟件大多是為航空、航天、航海等問題計(jì)算而開發(fā)，對10 km 級及以上的大范圍風(fēng)場仿真的工程應(yīng)用還比較少見。我國依托數(shù)值計(jì)算開展的風(fēng)力發(fā)電技術(shù)創(chuàng)新研究也面臨著CFD 軟件使用受限的風(fēng)險(xiǎn)。鑒于此，本文避開通用軟件市場，在專業(yè)軟件方面精準(zhǔn)發(fā)力，開發(fā)了能夠模擬仿真我國不同地形、不同大氣穩(wěn)定度特征的自主軟件LSWS，并開展了典型地形的仿真驗(yàn)證工作。本文所述典型地形是指我國具有較強(qiáng)風(fēng)資源開發(fā)潛力區(qū)域所包含的主要地形，包括：內(nèi)蒙古草原、丘陵、長江中下游平原、山地、秦嶺山脈、高原山地等。在研發(fā)大氣CFD模式的過程中，選取了典型地形中的兩種地形開展研究，包括內(nèi)蒙草原錫林浩特平坦草原地形和山西丘陵地形。在這兩種地形下，我國建有大量的風(fēng)電場，對于研究風(fēng)特性、風(fēng)能特性有重要的參考意義。文章第一節(jié)為大氣CFD 模式及邊界條件，詳細(xì)介紹了大氣CFD 模式理論，基于錫林浩特風(fēng)場實(shí)測數(shù)據(jù)提出的修正Gryning 風(fēng)廓線模型和擬合廓線模型，以及所用計(jì)算方法和自主軟件架構(gòu)；第二節(jié)為我國典型錫林浩特平坦草原地形、山西丘陵地形風(fēng)場計(jì)算結(jié)果，并與商業(yè)軟件Fluent 的計(jì)算結(jié)果進(jìn)行了對比；第三節(jié)給出了結(jié)論及建議。

1 理論及方法

1.1 大氣CFD 模式及邊界條件

1.1.1 大氣CFD 模式

為了實(shí)現(xiàn)對大氣邊界層風(fēng)場的模擬，考慮不同大氣穩(wěn)定度的影響，在質(zhì)量守恒方程和動量守恒方程的基礎(chǔ)上引入了位溫方程?？刂品匠探M可以寫為：

其中：ui為i方向的平均速度分量；μ為空氣動力黏度；μt為湍流渦黏；P為壓力；ρ 為密度；σθ為位勢溫度θ的施密特?cái)?shù)，θ的表達(dá)式為 θ=T(P0/P)R/cp；T為大氣溫度；P0≈1×105Pa 為標(biāo)準(zhǔn)大氣壓；R=8.313 J/(mol·K)為普適氣體常數(shù)；cp≈1 J/(g·K)為干空氣在恒定壓力下的比熱容。

大氣邊界層內(nèi)部的大尺度運(yùn)動受到熱分層和科里奧利效應(yīng)的影響，分別由地球自轉(zhuǎn)和慣性質(zhì)量引起。當(dāng)對整個(gè)ABL 進(jìn)行建模時(shí)，這些效應(yīng)可以通過附加的源項(xiàng)（Sv）引入RANS 方程組：

其中：gi為重力加速度，gi=(0,0,-g)T；ρ0為參考密度；εi=(-1,1,0)T；fc=2ΩEsinλ為科里奧力參數(shù)，ΩE為地球自轉(zhuǎn)速度，λ為緯度。

1.1.1.1 中性大氣湍流閉合方程

使用k-ε湍流模型封閉上述控制方程組。適用于中性ABL 的k-ε模型表達(dá)式為：

其中：σk和 σε分別為k和 ε 的施密特?cái)?shù)；Cε1和Cε2為模型系數(shù)；Pk為湍動能的生成項(xiàng)。通過求解湍動能k和耗散率 ε的兩個(gè)偏微分輸運(yùn)方程，得到湍流渦黏。由此得到的混合長度lt和渦流黏度 μt用k和 ε可表示為：

其中Cμ為常數(shù)。

1.1.1.2 非中性大氣湍流閉合方程

在非中性條件下，為了考慮熱穩(wěn)定性的影響，本文參考Koblitz[7]論文中的理論公式，對k-ε模型進(jìn)行修正，最大的特色是通過經(jīng)驗(yàn)的湍動能理論估計(jì)對全局最大混合長度進(jìn)行顯式積分計(jì)算，避免了原文中繁瑣且沒有理論依據(jù)的預(yù)前模擬才能求解的限制。

湍動能方程中增加的浮力源項(xiàng)B與系數(shù)Cε3一起出現(xiàn)在耗散方程。浮力源項(xiàng)B的表達(dá)式為：

在熱不穩(wěn)定條件下，B為正，支持湍動能的產(chǎn)生；而在穩(wěn)定條件下，B為負(fù)，抑制湍動能的產(chǎn)生。

附加擴(kuò)散項(xiàng)D的表達(dá)式為：

其中，系數(shù)Cμ為常數(shù)。

其中，le為全局最大混合長度。在中性條件下的平坦草原地形上，le=l0=0.000 27G/fc，G為地轉(zhuǎn)風(fēng)速。對于非中性條件和復(fù)雜地形，有：

通常取系數(shù)α=0.075，用于保證在中性條件下l0和lMY兩/個(gè)長度尺度完全相同。對湍動能k取理論估計(jì)并忽略大氣邊界層高度以外的湍動能，這里u*為大氣邊界層的摩擦速度。

式（15）可近似表達(dá)為lMY=αhg/2 。hg=u*/(6fc)為大氣邊界層高度的估算。αB的表達(dá)式為：

其中，理查德森數(shù)Rig=B/Pk，B依賴于渦黏 μt、重力加速度gi以及密度 ρ。

在數(shù)值模擬過程中對多種湍動能入口廓線進(jìn)行了對比。對湍動能廓線的一種擬合形式為：

其中，C1和C2為擬合系數(shù)。針對這一湍動能廓線形式，對lMY的積分表達(dá)式進(jìn)行積分可得：

數(shù)值模擬系數(shù)取值見表1。

表1 數(shù)值模擬系數(shù)取值Table 1 Values of the numerical simulation coefficients

1.1.2 壁面函數(shù)

在地形表面附近，通過壁面函數(shù)來計(jì)算地面上方第一層網(wǎng)格中心高度的函數(shù)值。在本文計(jì)算中，熱穩(wěn)定度的影響通過入口廓線來體現(xiàn)，壁面函數(shù)的形式基于平衡假設(shè)來推導(dǎo)，唯一的湍流能量產(chǎn)生機(jī)制是剪切作用，不考慮浮力等因素的影響。對數(shù)速度剖面的表達(dá)形式為：

其中，E為對數(shù)分布律中的常數(shù)，通常取8.143。由式（19）可推導(dǎo)出相關(guān)物理量的壁面函數(shù)。速度的壁面函數(shù)通過壁面切應(yīng)力 τ0來實(shí)現(xiàn)，

其中：kp為網(wǎng)格中心處的湍動能；up為網(wǎng)格中心處的速度；為地面到網(wǎng)格中心zp的無量綱距離，

第一層網(wǎng)格中心的溫度通過壁面熱流密度q0來計(jì)算，

其中：Tp為第一層網(wǎng)格中心處的溫度；T0為壁面處的溫度；λt為第一層網(wǎng)格中心與壁面之間區(qū)域的當(dāng)量導(dǎo)熱系數(shù)，且有關(guān)系式：

其中：σl為分子Prandtl 數(shù)（即Pr）；σT為湍流Prandtl 數(shù)。

1.2 錫林浩特風(fēng)速測量實(shí)驗(yàn)

入口條件是影響大氣CFD 模式準(zhǔn)確性的重要因素。以內(nèi)蒙古自治區(qū)錫林浩特大氣邊界層湍流觀測實(shí)驗(yàn)場的實(shí)測風(fēng)數(shù)據(jù)對大氣CFD 模式的入口廓線進(jìn)行標(biāo)定。該實(shí)驗(yàn)場位于錫林浩特東北方向，距離錫林浩特約30 km，為高原草場，地形平坦，零星分布一些低山丘陵，可以代表內(nèi)陸平坦地形。前期朱蓉等[12]已經(jīng)對該地區(qū)的風(fēng)能資源和氣候特征做了詳細(xì)研究。該地區(qū)設(shè)立了100 m 觀測塔和70 m 觀測塔各1 個(gè)，本文所用數(shù)據(jù)的觀測時(shí)段為2009 年6 月至2010年7 月。測風(fēng)塔外觀如圖1 所示。對100 m 高度以下的風(fēng)速，采用測風(fēng)塔實(shí)測風(fēng)數(shù)據(jù)。對于100 m 以上高度的風(fēng)速，測風(fēng)塔難以測量，采用激光雷達(dá)測風(fēng)儀進(jìn)行測量。

圖1 錫林浩特風(fēng)速測量實(shí)驗(yàn)Fig.1 Xilinhot wind measurement experiment

1.3 入口邊界條件

本節(jié)給出了中性、非中性風(fēng)廓線的經(jīng)典模型及修正方法，對比了兩種經(jīng)典模型描述我國典型地形風(fēng)特性的優(yōu)缺點(diǎn)，基于錫林浩特實(shí)測風(fēng)數(shù)據(jù)對經(jīng)典模型進(jìn)行了修正，提出了一種修正模型，同時(shí)基于分段擬合方法提出了一種擬合廓線模型，這兩種模型均可作為大氣CFD 模式入口邊界條件。

1.3.1 廓線修正方法

大氣邊界層可以垂直分為層流底層、地表層、?？寺鼘? 個(gè)部分：最低部分稱為層流底層，其高度等于氣動粗糙度高度z0；z0以上是湍流充分發(fā)展的地表層，其垂直范圍可能從20 m 到100 m 不等[5]，具體取決于空氣的熱分層，通常為整個(gè)大氣邊界的10%；地表層之上是?？寺鼘樱涓叨瓤赡艹^1 000 m，具體取決于空氣的穩(wěn)定性、科里奧利參數(shù)和地表的粗糙度高度。

風(fēng)速隨著高度的垂直升高而升高，是最重要的大氣特征。最經(jīng)典的方法是使用對數(shù)函數(shù)描述。在地表層：

式（25）只適用于中性條件。對于非中性條件，可采用莫寧-奧布霍夫相似理論對式（25）進(jìn)行修正，

對于不穩(wěn)定條件，修正方程[13-14]為：

對于穩(wěn)定條件，對數(shù)風(fēng)廓線修正方程為[15-17]：

其中：a=5；A=1；B=2/3；C=5；D=0.35。

當(dāng)從地表層到達(dá)艾克曼層時(shí)，地球自轉(zhuǎn)產(chǎn)生的科氏力不可忽略，上述方法不能準(zhǔn)確描述速度隨高度的變化，需要對整個(gè)大氣邊界層的風(fēng)廓線進(jìn)行統(tǒng)一描述。目前常用的描述方法可分為兩種。

1.3.2 EKM 模型

第一個(gè)方法是采用分段擬合，并確保在兩個(gè)區(qū)域內(nèi)能夠光滑過渡[18]。該模型在下文中稱為EKM 模型：

其中：ug為地轉(zhuǎn)風(fēng)速；zp為地表層厚度；α0為表面層風(fēng)與地轉(zhuǎn)風(fēng)的夾角。對中性條件，

對非中性條件，

ψ(zp/L)在穩(wěn)定和不穩(wěn)定條件下不同。對于穩(wěn)定條件，

式中，a的取值在4.7～5 范圍內(nèi)。對于不穩(wěn)定條件，

1.3.3 修正Gryning 模型

第二個(gè)方法是修正混合長度對高度的依賴性，其由Gryning等[19]提出。他們用公式表達(dá)了與高度相關(guān)的混合長度，以限定其隨著高度的生長，從而把對數(shù)法則式的有效性延伸到了表面層以上。對于穩(wěn)定條件：

對中性和不穩(wěn)定條件：

式中，zi為大氣邊界層高度。

其他幾項(xiàng)可由以下公式獲得：

式（36）、式（37）在下文中稱為Gryning 模型。不同穩(wěn)定度下，風(fēng)速隨高度變化的規(guī)律如圖2 所示。圖中方塊為基于測風(fēng)塔和激光雷達(dá)的實(shí)測數(shù)據(jù)，實(shí)線為EKM 模型和Gryning 模型的擬合曲線?？梢钥闯?，在100 m 高度以下時(shí)，兩種模型的擬合曲線與實(shí)測值均能良好吻合，其中Gryning 模型稍優(yōu)于EKM 模型。隨著高度增加，尤其是在穩(wěn)定條件下，擬合曲線與實(shí)測值之間的偏差增大。這是由于理論模型高估了速度隨高度增長的趨勢，導(dǎo)致其在更高的高度下的速度估值偏高。

圖2 風(fēng)速隨高度變化規(guī)律Fig.2 Variation of wind speed with height

相比于EKM 模型，Gryning 模型結(jié)構(gòu)更加簡單。因此，基于實(shí)測數(shù)據(jù)，對方程（36）、方程（37）右側(cè)的第二項(xiàng)和最后一項(xiàng)進(jìn)行重新擬合修正，獲得Gryning修正模型。對于穩(wěn)定條件：

對中性和不穩(wěn)定條件：

Gryning 修正模型與實(shí)測數(shù)據(jù)的對比如圖3 所示，修正模型與實(shí)測數(shù)據(jù)吻合良好。修正模型較好地修正了原模型與實(shí)測值隨高度增加而存在的偏差。因此，在保證Gryning 模型整體結(jié)構(gòu)不變的情況下，僅需要對其相關(guān)常數(shù)項(xiàng)進(jìn)行簡單的修正，即可與實(shí)測值吻合。

圖3 Gryning 修正模型與實(shí)測值的對比Fig.3 Comparison between the modified Gryning model and the measured value

1.3.4 擬合廓線模型

除使用理論模型外，基于測風(fēng)塔和激光雷達(dá)的實(shí)測數(shù)據(jù)，結(jié)合分段函數(shù)擬合方法和理論公式，本文構(gòu)建了不同熱穩(wěn)定度下的風(fēng)速u隨高度z變化的入口廓線模型，描述如下。

當(dāng)熱穩(wěn)定度分級小于或等于3 時(shí)：

當(dāng)熱穩(wěn)定度分級為4 或5 時(shí)：

其中：u*為通過實(shí)測獲得的地表摩擦速度；u4為通過實(shí)驗(yàn)測量獲得的4 m 處風(fēng)速；z0為地表摩擦度；L為莫寧-奧布霍夫長度；κ=0.41為馮·卡門常數(shù)；ug為地轉(zhuǎn)風(fēng)速，

A、B與穩(wěn)定度有關(guān)（Garratt，1994）：

c和c′為與熱穩(wěn)定度有關(guān)的常數(shù)，由表2 可獲得。

表2 不同熱穩(wěn)定度下擬合參數(shù)Table 2 Fitting parameters under different thermal stability degrees

擬合曲線與實(shí)測數(shù)據(jù)對比如圖4 所示，模型擬合曲線在不同穩(wěn)定度下均能與實(shí)測數(shù)據(jù)良好吻合，該模型能夠?yàn)榈湫偷匦畏侵行源髿釩FD 模式計(jì)算提供準(zhǔn)確的速度入口條件。

圖4 擬合模型曲線與實(shí)測數(shù)據(jù)對比Fig.4 Comparison between the fitting model curve and the measured data

1.3.5 入口湍動能

在大氣邊界層流動中，被廣泛采用的湍流動能k和動能耗散率ε的廓線可由以下兩個(gè)方程給出：

其中，Cμ一般取0.03。本文基于Zhang[5]的工作，引入一隨高度變化的修正項(xiàng)，修正上述湍動能方程：

圖5（a）給出了實(shí)測數(shù)據(jù)與方程（53）曲線的對比。由圖可見，僅在穩(wěn)定條件下，方程（53）的預(yù)測值與實(shí)測值吻合良好；在弱穩(wěn)定條件和不穩(wěn)定條件下，方程（53）低估了湍動能；在中性和弱不穩(wěn)定條件下，方程（53）高估了湍動能。基于此，對方程（53）進(jìn)行重新擬合修正，新的湍動能廓線由下式給出：

圖5 湍動能隨高度變化規(guī)律Fig.5 Variation of turbulent kinetic energy with height

式中，C1和C2為通過擬合獲得的常數(shù)，其在不同穩(wěn)定條件下的值由表3 獲得。圖5（b）給出了修正后的曲線與實(shí)測值的對比。從圖中可以看出，修正后的曲線與實(shí)測值吻合良好。

表3 不同熱穩(wěn)定度下湍動能方程擬合參數(shù)Table 3 Fitting parameters for the turbulent kinetic energy equation under different thermal stability degrees

1.4 軟件設(shè)計(jì)及計(jì)算方法

1.4.1 軟件設(shè)計(jì)及功能描述

近地面風(fēng)場特性受地表起伏影響較大。在復(fù)雜地形的近地面形成精細(xì)湍流結(jié)構(gòu)進(jìn)而對復(fù)雜地形開展CFD 模擬的精度要求很高，導(dǎo)致網(wǎng)格數(shù)量龐大、計(jì)算耗時(shí)長。傳統(tǒng)的中尺度大氣模擬方法雖然計(jì)算速度快，但分辨率在10 km 量級，無法得到精細(xì)化的精確風(fēng)場數(shù)據(jù)。為了實(shí)現(xiàn)精細(xì)化大氣風(fēng)場模擬，本文在上述理論基礎(chǔ)上，開發(fā)了一款具有自主知識產(chǎn)權(quán)的風(fēng)場CFD 模擬軟件，能夠快速、準(zhǔn)確地得到典型地形的風(fēng)場結(jié)果，實(shí)現(xiàn)不同大氣穩(wěn)定度條件下的風(fēng)場分類模擬，可為風(fēng)資源評估、風(fēng)電場微觀選址和進(jìn)一步風(fēng)能預(yù)測提供技術(shù)支撐，并可節(jié)省大量實(shí)驗(yàn)成本和建設(shè)成本。

自主軟件LSWS，即考慮大氣穩(wěn)定度的典型地形風(fēng)場模擬軟件，版本號V1.0。LSWS 軟件工作流程及軟件構(gòu)架如圖6 所示。軟件將大氣模式與CFD 結(jié)合起來，為典型地形風(fēng)場模擬提供了大氣CFD 模式及模擬方法。軟件提供了完備的各種典型地形風(fēng)場求解方案，入口邊界條件可采用多種方案（包括經(jīng)典的傳統(tǒng)大氣穩(wěn)定度廓線、自主修正的Gryning 風(fēng)廓線等模型），并以實(shí)測數(shù)據(jù)來驅(qū)動模式運(yùn)轉(zhuǎn)，并可以根據(jù)使用者的需求來設(shè)置入口參數(shù)。軟件嵌入了適用于我國典型地形、具有不同穩(wěn)定度修正的廓線模型和相關(guān)參數(shù)，在CFD 通用方程的基礎(chǔ)上增加了熱浮力、科氏力等大氣運(yùn)動對流場及湍流的影響模型，具備復(fù)雜地形風(fēng)場大分離、非線性的求解能力，可模擬風(fēng)場瞬態(tài)和穩(wěn)態(tài)流動特征，可有效求解典型地形的包含時(shí)間項(xiàng)在內(nèi)的風(fēng)場結(jié)構(gòu)。軟件包含適用于典型地形風(fēng)場模擬的網(wǎng)格預(yù)處理、幾何量處理、求解器與后處理模塊，簡單易用，且具有跨平臺使用能力。

圖6 LSWS 軟件工作流程及軟件架構(gòu)圖Fig.6 Workflow and software architecture of LSWS

1.4.2 計(jì)算方法

控制方程的離散采用隱式有限體積法。根據(jù)導(dǎo)數(shù)項(xiàng)離散采用的不同方法和適用的網(wǎng)格結(jié)構(gòu)類型，離散方法又分為采用局部坐標(biāo)變換的方法和無坐標(biāo)變換的方法。本軟件的一個(gè)特色是給出了無坐標(biāo)變換方法在斜網(wǎng)格情況下的修正技術(shù)，使用該技術(shù)分別給出了對流項(xiàng)、黏性項(xiàng)和源項(xiàng)的斜網(wǎng)格修正離散方法。對修正QUICK 類迎風(fēng)格式，利用邊界零體積和同位網(wǎng)格的處理技術(shù)，給出了邊界無需降階和外插的處理方案。同時(shí)，給出了斜網(wǎng)格修正的邊界條件處理方法。軟件還開發(fā)了每個(gè)單元的通量守恒自檢功能，可以通過該功能相應(yīng)地調(diào)整網(wǎng)格和計(jì)算方法，來保證計(jì)算過程中的通量守恒性、收斂性和計(jì)算結(jié)果的準(zhǔn)確性。具體細(xì)節(jié)可參閱文獻(xiàn)[20]。

2 典型地形算例及分析

本文選取錫林浩特和山西南樺山兩處典型地形進(jìn)行大氣CFD 模式驗(yàn)證。錫林浩特為平坦的草原地形，山西南樺山為丘陵地形，兩者具有代表性，且前期在這兩地開展了長期風(fēng)速測量實(shí)驗(yàn)，數(shù)據(jù)較為完備。圖7（a）為錫林浩特平坦草原地形，圖7（b）為山西南樺山丘陵地形，圖7（c）為南樺山測風(fēng)點(diǎn)位的山前、山脊和山后分布情況。經(jīng)過調(diào)研和數(shù)值實(shí)驗(yàn)，本文建立了錫林浩特平坦草原地形計(jì)算域模型和山西南樺山丘陵地形計(jì)算域模型，如圖8 所示。計(jì)算域整體采用圓形，便于不同風(fēng)向的批量計(jì)算。在計(jì)算域建模方面，將計(jì)算域分成平坦區(qū)、過渡區(qū)和核心區(qū)三部分，可使本文提出的風(fēng)廓線模型既能應(yīng)用于平坦地形又能應(yīng)用于復(fù)雜地形，同時(shí)從計(jì)算流體力學(xué)角度可使計(jì)算域整體滿足守恒性，加速收斂，且能避免復(fù)雜地形“人造斷崖”帶來的誤差干擾。

圖7 錫林浩特和山西南樺山風(fēng)場位置Fig.7 Locations of the Xilinhot and Nanhuashan wind fields

圖8 計(jì)算域及網(wǎng)格拓?fù)銯ig.8 Computational domain and grid topology

由于錫林浩特為平坦草原地形，關(guān)注區(qū)域?yàn)? km范圍內(nèi)的核心區(qū)，通過一定的平滑過渡和拉平得到過渡區(qū)和平坦區(qū)，最終圓形計(jì)算域半徑為8 km，計(jì)算域高度為0.4 km。山西南樺山計(jì)算域，考慮到丘陵的相互影響，關(guān)注區(qū)域在15 km 范圍內(nèi)的核心區(qū)，山的當(dāng)?shù)叵鄬Ω叨葹?50 m，在山的西側(cè)有1 號和2 號測風(fēng)點(diǎn)，山的東側(cè)有3 號測風(fēng)點(diǎn)。同樣通過平滑過渡到平坦區(qū)，最終山西南樺山圓形計(jì)算域半徑為18 km，計(jì)算域高度為3 km。因兩處的實(shí)際主導(dǎo)風(fēng)向均為西南風(fēng)，因此本文以西南風(fēng)方向?yàn)槔M(jìn)行仿真模擬，使用的入口廓線為1.3.4 節(jié)提出的廓線模型。按照計(jì)算域中的結(jié)構(gòu)拓?fù)?，利用商業(yè)軟件進(jìn)行結(jié)構(gòu)網(wǎng)格劃分，首層網(wǎng)格高度為3 m，核心區(qū)水平分辨率控制在30 m左右，兩款軟件采用相同網(wǎng)格和邊界條件。

2.1 典型錫林浩特平坦草原地形仿真

分別利用自主軟件LSWS 和商業(yè)軟件Fluent，對錫林浩特典型平坦草原地形開展CFD 仿真。利用錫林浩特實(shí)測風(fēng)數(shù)據(jù)標(biāo)定修正廓線模型，選取不穩(wěn)定、中性和穩(wěn)定三種穩(wěn)定度廓線作為入口條件。選取與測風(fēng)設(shè)備點(diǎn)位一致的3 個(gè)測點(diǎn)的風(fēng)速廓線，進(jìn)行歸一化處理，可得到2 個(gè)測點(diǎn)位置的風(fēng)廓線信息，如圖9所示。圖9（d）為錫林浩特平坦草原地形的海拔高度及測風(fēng)設(shè)備的空間位置。從圖9 中可以看出，LSWS和Fluent 計(jì)算的結(jié)果總體趨勢一致，與實(shí)測風(fēng)數(shù)據(jù)吻合度較高，兩軟件均能很好地反映出不同穩(wěn)定度廓線的垂直分布情況，地形對廓線的影響較小。對于不穩(wěn)定狀態(tài)（圖9a），風(fēng)速分布斜率較大；中性狀態(tài)（圖9c），斜率接近對數(shù)率分布特征；穩(wěn)定狀態(tài)（圖9e），斜率相對平緩。對于穩(wěn)定狀態(tài)風(fēng)速廓線的計(jì)算而言，LSWS相對Fluent，在捕捉廓線廓形上的精度相對更好。

圖9 測點(diǎn)位置速度廓線Fig.9 Velocity profiles at the testing locations

測點(diǎn)1 和測點(diǎn)2 在70 m 高度處實(shí)測風(fēng)數(shù)據(jù)的誤差對比如表4 所示。自主軟件LSWS 商業(yè)軟件Fluent誤差減小量的計(jì)算方式為：Fluent 與實(shí)測風(fēng)速之間的絕對誤差E1與LSWS 與實(shí)測風(fēng)速之間的絕對誤差E2的差值的百分比，即(E1-E2) × 100%?？傮w而言，在平坦地形條件下，LSWS 與Fluent 總體精度相當(dāng)。所有穩(wěn)定度下，測點(diǎn)1 的誤差在1%以內(nèi)；測點(diǎn)2 在穩(wěn)定條件下，LSWS 計(jì)算精度比Fluent 提升6.33%。

表4 70 m 高度測風(fēng)數(shù)據(jù)誤差對比Table 4 Error comparison of the wind measurement data at the height of 70 m

2.2 典型山西南樺山丘陵地形仿真

利用自主軟件LSWS 和商業(yè)軟件Fluent，對山西南樺山丘陵地形進(jìn)行仿真計(jì)算對比研究，同時(shí)進(jìn)一步比較LSWS 求解復(fù)雜地形風(fēng)場的能力。同樣采用本文的修正廓線模型作為入口邊界條件。通過比較CFD 模擬計(jì)算結(jié)果與實(shí)測結(jié)果，考察大氣CFD 模式及模擬方法對山區(qū)復(fù)雜地形的影響、風(fēng)場的空間分布、加減速效應(yīng)的模擬精度和計(jì)算能力。

圖10 顯示了三種穩(wěn)定度下、距地面10 m 高度處的速度矢量圖。自主軟件LSWS 和商業(yè)軟件Fluent計(jì)算的風(fēng)場分布整體相似。相比之下，F(xiàn)luent 計(jì)算的山區(qū)復(fù)雜地形迎風(fēng)側(cè)的風(fēng)速偏高、背風(fēng)側(cè)山谷處的風(fēng)速偏低。

圖10 三種穩(wěn)定度下距地面10 m 高度處速度矢量圖Fig.10 Velocity vectors at 10 m above the ground under three stability degrees

自主軟件LSWS 和商業(yè)軟件Fluent 計(jì)算得到的三種穩(wěn)定度下的風(fēng)場云圖如圖11 所示。從圖中可以看出，對于不穩(wěn)定狀態(tài)，兩軟件的計(jì)算結(jié)果總體趨勢相同；而對于中性狀態(tài)和穩(wěn)定狀態(tài)的風(fēng)速分布，LSWS 得到的邊界層分布更符合真實(shí)大氣邊界層狀態(tài)，F(xiàn)luent 計(jì)算得到的邊界層形成過快、速度梯度較大且對地形的起伏變化不敏感。

圖11 三種穩(wěn)定度下的風(fēng)場云圖Fig.11 Wind field contours under three stability degrees

圖12 為3 個(gè)典型測點(diǎn)位置的風(fēng)速廓線分布對比。整體而言，自主軟件LSWS 模擬的不同穩(wěn)定度狀態(tài)下的風(fēng)速廓線均與實(shí)驗(yàn)值接近，可以較為準(zhǔn)確地刻畫山地地形動力強(qiáng)迫作用對山前、山脊和山后風(fēng)速廓線的影響。不同穩(wěn)定度入口廓線對山地風(fēng)速廓線的分布影響較小，地形的作用占主導(dǎo)作用。從廓線的總體垂直分布來看，F(xiàn)luent 在模擬山地風(fēng)速分布上，速度垂直梯度變化較快，與山地邊界層風(fēng)場特性有一定的偏差。

圖12 測點(diǎn)位置風(fēng)廓線Fig.12 Velocity profiles at the testing locations

1 號、2 號和3 號測風(fēng)塔在100 m 高度處實(shí)測風(fēng)數(shù)據(jù)的誤差對比如表5 所示?？傮w而言，在山區(qū)復(fù)雜地形條件下，自主軟件LSWS 的計(jì)算精度明顯優(yōu)于商業(yè)軟件Fluent，所有穩(wěn)定度下的計(jì)算精度均取得提升。對于山前1 號測風(fēng)塔，在不穩(wěn)定、中性和穩(wěn)定條件下，LSWS 與Fluent 計(jì)算精度相當(dāng)。隨著氣流繞過山脊，對于山脊2 號測風(fēng)塔，在不穩(wěn)定、中性和穩(wěn)定條件下，LSWS 與Fluent 對比，計(jì)算精度最大可分別提升5.42%、5.34%、40.77%。就山后3 號測風(fēng)塔而言，在不穩(wěn)定、中性和穩(wěn)定條件下，LSWS 與Fluent 對比，計(jì)算精度最大可分別提升38.39%、39.57%和41.80%。由此可見，LSWS 對復(fù)雜山區(qū)地形影響風(fēng)場結(jié)果的模擬能力有顯著的提升。

表5 測風(fēng)塔100 m 高度測風(fēng)數(shù)據(jù)誤差對比Table 5 Error comparison of the wind data measured at the height of 100 m

3 結(jié)論及建議

本文根據(jù)風(fēng)場實(shí)測數(shù)據(jù)，提出了兩種能夠反映不同大氣穩(wěn)定度特征的風(fēng)廓線模型，建立了能夠求解不同大氣穩(wěn)定度狀態(tài)的大氣CFD 模式，改進(jìn)了非中性大氣湍流閉合方程及求解方法，提高了可實(shí)現(xiàn)性和計(jì)算效率，開發(fā)了自主可控的CFD 軟件，并利用錫林浩特平坦草原地形和山西丘陵地形實(shí)測風(fēng)數(shù)據(jù)進(jìn)行了驗(yàn)證。研究得出如下結(jié)論：

1）基于實(shí)測風(fēng)數(shù)據(jù)提出的符合我國地形與風(fēng)特性的修正Gryning 廓線模型和新的擬合廓線模型，可應(yīng)用于平坦和復(fù)雜地形風(fēng)場CFD 模擬。

2）通過建立的大氣CFD 模式及方法，對錫林浩特平坦草原地形風(fēng)場進(jìn)行仿真計(jì)算，對比了自主軟件LSWS 與商業(yè)軟件Fluent 的計(jì)算精度。在平坦地形下兩款軟件的計(jì)算精度總體相當(dāng)。對于穩(wěn)定大氣狀態(tài)，自主軟件LSWS 的計(jì)算精度較商業(yè)軟件Fluent 最高可提升6.33%。

3）對山西丘陵地形風(fēng)場進(jìn)行仿真計(jì)算，總體上自主軟件LSWS 的計(jì)算精度明顯高于商業(yè)軟件Fluent。在所有大氣條件下，自主軟件LSWS 較商業(yè)軟件Fluent 計(jì)算精度最高可提升41.80%。特別是對于繞山后流場，本文方法得到的結(jié)果更加合理準(zhǔn)確，對解決復(fù)雜山區(qū)風(fēng)特性評估等問題有積極的意義和工程應(yīng)用價(jià)值。

綜合多年CFD 軟件開發(fā)和商業(yè)軟件Fluent 使用經(jīng)驗(yàn)，作者認(rèn)為本文開發(fā)的軟件與目前求解域模型和網(wǎng)格更加適配，在計(jì)算過程中，還可以查看每一個(gè)單元網(wǎng)格的通量守恒性并進(jìn)行修正，防止了流場的數(shù)值擴(kuò)散，使得計(jì)算結(jié)果有更好精度。目前自主軟件LSWS 還處于起步階段，操作還不夠便捷，需進(jìn)一步優(yōu)化改進(jìn)交互界面，提升用戶體驗(yàn)感，進(jìn)而提升市場競爭力。

致謝：在大氣邊界層風(fēng)特性理論、風(fēng)場實(shí)測數(shù)據(jù)、實(shí)測數(shù)據(jù)處理方法等方面，得到了國家氣候中心朱蓉研究員、中國科學(xué)院大氣物理研究所程雪玲研究員、華北電力大學(xué)李莉副教授的支持和指導(dǎo)，在此表示感謝。