摘 要：針對(duì)風(fēng)資源精細(xì)化評(píng)估需求，提出采用實(shí)際測(cè)風(fēng)數(shù)據(jù)修正k-l 湍流模型的關(guān)鍵參數(shù)，以提高數(shù)值仿真對(duì)大尺度復(fù)雜地形流動(dòng)的準(zhǔn)確性和適用性。在此過程中，首先介紹k-l 模型的求解過程并開展網(wǎng)格分辨率研究；其次利用實(shí)際測(cè)風(fēng)數(shù)據(jù)修正模型參數(shù)并使用典型案例進(jìn)行驗(yàn)證；最后以中國(guó)華南地區(qū)某實(shí)際風(fēng)電場(chǎng)為對(duì)象，將其與k-ε 模型進(jìn)行對(duì)比分析以驗(yàn)證其有效性和優(yōu)越性。結(jié)果表明，相比于k-ε 模型，原始k-l 模型由于高估了尾流影響，導(dǎo)致預(yù)測(cè)誤差偏大；而修正參數(shù)后的k-l 模型能有效提升精度，風(fēng)加速因子和發(fā)電量等效小時(shí)數(shù)的平均相對(duì)誤差分別降低5.1%和7%。因此，在k-l 模型的實(shí)際應(yīng)用過程中，需利用實(shí)際測(cè)風(fēng)數(shù)據(jù)修正參數(shù)，或?qū)?shù)B1的量級(jí)調(diào)整為102。

關(guān)鍵詞：風(fēng)力發(fā)電；湍流模型；數(shù)值方法；復(fù)雜地形；測(cè)風(fēng)數(shù)據(jù)

中圖分類號(hào)：TK81 文獻(xiàn)標(biāo)志碼：A

0 引 言

風(fēng)資源評(píng)估是風(fēng)能資源開發(fā)規(guī)劃的關(guān)鍵步驟［1-2］，亦是保障風(fēng)電場(chǎng)投資效益的重要手段。隨著風(fēng)電產(chǎn)業(yè)規(guī)模的不斷加大，風(fēng)電場(chǎng)地形日趨復(fù)雜，對(duì)風(fēng)資源評(píng)估的準(zhǔn)確性和適應(yīng)性提出了更高要求?；诰€性求解器的傳統(tǒng)方法（如WAsP等）無法滿足復(fù)雜山地的評(píng)估需求，基于雷諾平均（Reynoldsaverage Navier-Stockes， RANS）的計(jì)算流體力學(xué)方法（computational fluid dynamics， CFD）憑借其計(jì)算精度高、地形適應(yīng)性好，同時(shí)兼顧計(jì)算資源等優(yōu)勢(shì)，近年來務(wù)受關(guān)注［3-7］。湍流模型作為RANS 方法的重要組成部分，直接決定了計(jì)算結(jié)果的準(zhǔn)確性，因此研究者們嘗試了多種手段以提升湍流模型的有效性和準(zhǔn)確性。

k-ε 湍流模型憑借其泛用性在各類型流動(dòng)仿真中得到廣泛使用，許多研究者嘗試將其應(yīng)用于風(fēng)電場(chǎng)尺度的流動(dòng)仿真［8-10］。Beljaars 等［8］以Askervein 山為對(duì)象使用k-ε 模型進(jìn)行了計(jì)算，得到的切應(yīng)力在邊界層區(qū)域是可以接受的，但在主流區(qū)域的預(yù)測(cè)值偏大。為使k-ε 模型適用于復(fù)雜地形流動(dòng)，Richards 等［9］從大氣邊界條件出發(fā)，對(duì)k-ε 模型參數(shù)進(jìn)行了調(diào)整，結(jié)果表明修正參數(shù)后可得到更接近實(shí)測(cè)數(shù)據(jù)的結(jié)果；劉文等［10］以湍動(dòng)能標(biāo)準(zhǔn)誤差為優(yōu)化目標(biāo)，通過正交試驗(yàn)優(yōu)選k-ε 模型參數(shù)，對(duì)比其他組參數(shù)結(jié)果，風(fēng)速、湍動(dòng)能和耗散率的誤差分別可降低1%、3% 和1.6%。

為了進(jìn)一步提高復(fù)雜地形流動(dòng)求解的準(zhǔn)確性和適用性，Mellor 等［11］發(fā)展了k-l 湍流模型，推導(dǎo)了該模型參數(shù)的物理意義，并根據(jù)典型流動(dòng)的實(shí)驗(yàn)結(jié)果確定了關(guān)鍵參數(shù)取值。在此基礎(chǔ)上，Yamada［12］對(duì)k-l 模型進(jìn)行了簡(jiǎn)化，并指出可利用代數(shù)表達(dá)式來確定l 取值，以提高模型計(jì)算效率。Arritt［13］使用簡(jiǎn)化后的k-l 湍流模型對(duì)湖陸風(fēng)現(xiàn)象進(jìn)行了仿真，結(jié)果表明模型預(yù)測(cè)值與實(shí)測(cè)結(jié)果一致性良好。此外，在風(fēng)資源評(píng)估領(lǐng)域得到廣泛使用的商業(yè)軟件Meteodyn WT 亦采用k-l 模型［14］作為默認(rèn)湍流模型。

k-l 模型自提出之初根據(jù)典型的小尺度流動(dòng)（尺度為0.1～1.0 m），如管道流［15］、槽道流［16］、平板流［17］和剪切流［18］等實(shí)驗(yàn)結(jié)果確定了模型參數(shù)，然而對(duì)于風(fēng)電場(chǎng)尺度（可達(dá)102～104 m）的流動(dòng)仿真而言，兩者存在量級(jí)上的差異，這可能導(dǎo)致根據(jù)前者確定的模型參數(shù)并不適用于大尺度的復(fù)雜地形流動(dòng)。因此，本文首先開發(fā)k-l 模型的相關(guān)代碼；然后根據(jù)實(shí)際測(cè)風(fēng)數(shù)據(jù)對(duì)k-l 模型的關(guān)鍵參數(shù)進(jìn)行修正；最后利用典型案例和實(shí)際案例對(duì)比不同湍流模型的性能，并分析流場(chǎng)差異，以驗(yàn)證k-l 模型參數(shù)修正的有效性和準(zhǔn)確性。

1 方 法

1.1 數(shù)值方法

本文利用開源軟件OpenFOAM（Open Source FieldOperation and Manipulation）求解復(fù)雜地形的流場(chǎng)結(jié)構(gòu)，數(shù)值方法中的控制方程為三維、穩(wěn)態(tài)、不可壓縮的雷諾平均Navier-Stokes 方程，采用SIMPLE 算法［19］求解控制方程。速度對(duì)流項(xiàng)使用二階迎風(fēng)格式離散，擴(kuò)散項(xiàng)使用二階中心格式離散。在湍流模型方面，本文采用k-l 湍流模型［12］，其中湍動(dòng)能k 的輸運(yùn)方程為：

基于文獻(xiàn)［12，20-22］，在OpenFOAM 中實(shí)現(xiàn)k-l 湍流模型的相關(guān)代碼，其求解流程采用偽代碼進(jìn)行說明，如表1 所示。在完成當(dāng)前迭代步的流場(chǎng)求解后，根據(jù)邊界條件和流場(chǎng)對(duì)湍動(dòng)能方程進(jìn)行迭代求解，直至滿足殘差要求，并用于下一迭代步的求解。

1.2 算例設(shè)置

為保持各風(fēng)向下計(jì)算域的一致性，本文中算例的物理模型均設(shè)置為圓柱體（如圖1 所示）其中區(qū)域1 為內(nèi)部域，根據(jù)實(shí)際地形生成；區(qū)域2 為過渡域，實(shí)現(xiàn)由實(shí)際地形向平地過渡；區(qū)域3 為外部域，地形為平地。計(jì)算域的生成步驟為：1）根據(jù)實(shí)際地形文件生成三維物理模型；2）在物理模型基礎(chǔ)上利用開源軟件cfMesh 生成網(wǎng)格；3）設(shè)置邊界條件。其中，頂面采用對(duì)稱面邊界；底面采用絕熱無滑移壁面；圓柱環(huán)面采用入流出流邊界條件，其特征是根據(jù)風(fēng)向確定當(dāng)?shù)剡吔珙愋?，即?dāng)環(huán)面段處空氣流入時(shí)，邊界條件類型為速度入口，此時(shí)采用Richards 等［9］提出的入口邊界條件，如式（10）～式（12）所示。

式中：uτ——摩擦速度，m/s；z0——粗糙長(zhǎng)度，m。

當(dāng)環(huán)面段空氣流出時(shí)，邊界條件為壓力出口；采用該邊界條件的優(yōu)勢(shì)是：能夠使用連續(xù)不間斷的環(huán)面作為計(jì)算域邊界，從而實(shí)現(xiàn)在不同工況下復(fù)用同一套網(wǎng)格。

1.3 網(wǎng)格分辨率

為減少網(wǎng)格尺度對(duì)數(shù)值模擬結(jié)果的影響，本節(jié)以典型案例Bolund 島［23-24］為研究對(duì)象開展網(wǎng)格分辨率研究。圖2 顯示了Bolund 島及其附近的地形分布，其中白色圓點(diǎn)為測(cè)風(fēng)塔，來流主風(fēng)向?yàn)?39°。利用2.2 節(jié)中的方法生成計(jì)算網(wǎng)格，計(jì)算域中內(nèi)部域的水平網(wǎng)格分辨率分別設(shè)置為20、10、5 m，對(duì)應(yīng)粗、中、細(xì)尺度網(wǎng)格，其中中尺度網(wǎng)格如圖3 所示。

同一工況下，不同網(wǎng)格各測(cè)風(fēng)塔點(diǎn)位的風(fēng)速分布如圖4所示，其中風(fēng)速采用來流風(fēng)速進(jìn)行無量綱化?？梢钥吹?，粗尺度網(wǎng)格相比中、細(xì)尺度網(wǎng)格的風(fēng)速誤差較大，平均誤差可達(dá)3.4%；而中尺度網(wǎng)格的結(jié)果與細(xì)尺度網(wǎng)格基本一致。因此，本文后續(xù)采用中尺度網(wǎng)格進(jìn)行計(jì)算。

2 結(jié)果與討論

為提高k-l 湍流模型對(duì)大尺度的復(fù)雜地形流動(dòng)的適用性，本節(jié)首先利用實(shí)際測(cè)風(fēng)數(shù)據(jù)修正模型參數(shù)，然后分別以Bolund島和某實(shí)際風(fēng)電場(chǎng)為對(duì)象開展研究，分析不同參數(shù)下k-l 湍流模型和k-ε 湍流模型的性能差異，并與實(shí)測(cè)結(jié)果進(jìn)行對(duì)比。

2.1 參數(shù)修正及驗(yàn)證

2.1.1 參數(shù)修正

根據(jù)式（1）可以看到B1 是影響湍動(dòng)能方程求解的關(guān)鍵參數(shù)，文獻(xiàn)［11］指出B1 的計(jì)算應(yīng)遵循以下關(guān)系：

式中：uτ——摩擦速度，m/s；τ——切應(yīng)力，Pa；ρ——密度，kg/m3。

在此基礎(chǔ)上，文獻(xiàn)［11］根據(jù)典型小尺度流動(dòng)的實(shí)驗(yàn)結(jié)果（表2）計(jì)算了B1 的取值。然而，基于小尺度流動(dòng)計(jì)算得到的B1 可能是不適用于大尺度的復(fù)雜地形流動(dòng)。因此，本節(jié)基于Bolund 島的實(shí)測(cè)數(shù)據(jù)［23-24］對(duì)B1 進(jìn)行修正，并驗(yàn)證修正參數(shù)后k-l 模型的有效性。

根據(jù)Bolund 島的實(shí)際測(cè)風(fēng)數(shù)據(jù)計(jì)算得到的B1 分布如圖5 所示，其中圓點(diǎn)為根據(jù)不同工況下各測(cè)點(diǎn)的實(shí)測(cè)數(shù)據(jù)計(jì)算得到的B1 值，點(diǎn)劃線為基于篩除異常值后的實(shí)測(cè)結(jié)果的擬合平均值?？梢钥吹?，對(duì)于該案例，B1 應(yīng)取119.8，相比原取值（16.6）存在量級(jí)上的差異。進(jìn)一步地，對(duì)于大尺度的復(fù)雜地形流動(dòng)可根據(jù)實(shí)測(cè)數(shù)據(jù)對(duì)湍流模型的參數(shù)進(jìn)行修正，數(shù)量級(jí)范圍為B1～102。

2.1.2 典型案例驗(yàn)證

為驗(yàn)證k-l 湍流模型參數(shù)修正的有效性和準(zhǔn)確性，使用不同參數(shù)的k-l 模型和k-ε 模型對(duì)Bolund 島進(jìn)行計(jì)算并對(duì)比分析。模型參數(shù)設(shè)置如表3 所示，其中G1 采用Richards 等［9］提出的參數(shù)，G2 采用劉文等［10］優(yōu)化后得到的參數(shù)，G3 為k-l模型的原始參數(shù)［12］，G4 為本文采用的k-l 模型參數(shù)。

圖6 和7 顯示了各模型的點(diǎn)位處風(fēng)速分布和點(diǎn)位M3 處風(fēng)廓線與實(shí)測(cè)結(jié)果的對(duì)比。可以看到，對(duì)于該案例，不同模型以及不同的模型參數(shù)并未對(duì)結(jié)果產(chǎn)生顯著影響，各組仿真得到的風(fēng)速分布較為接近；同時(shí)，風(fēng)廓線對(duì)比結(jié)果顯示不同模型及參數(shù)得到風(fēng)速隨高度的變化趨勢(shì)基本一致，G4 結(jié)果整體上更接近實(shí)測(cè)結(jié)果。

采用平均絕對(duì)誤差（mean absolute error， MAE）和平均相對(duì)誤差（mean relative error， MRE）評(píng)估模型精度，其定義式分別為：

式中：X——仿真值；X?——實(shí)測(cè)值。

表4 顯示了各組的誤差結(jié)果。從平均結(jié)果上來看，G4 的風(fēng)加速因子的MAE 和MRE 均低于其他組，僅為0.059 和8.8%。進(jìn)一步從整體分布上對(duì)風(fēng)加速因子的誤差進(jìn)行分析，其箱線圖如圖8 所示?？梢钥吹?，G4 的誤差的最大值和最小值均低于其他組，同時(shí)超過半數(shù)點(diǎn)位的誤差低于0.05。因此，參數(shù)修正后的k-l 模型在Bolund 島案例中預(yù)測(cè)精度高于其他模型，驗(yàn)證了修正參數(shù)的k-l 湍流模型的有效性。

2.2 實(shí)際案例對(duì)比

進(jìn)一步以中國(guó)華南地區(qū)某復(fù)雜山地風(fēng)電場(chǎng)為例，使用不同湍流模型進(jìn)行計(jì)算并分析結(jié)果。圖9 顯示了該風(fēng)電場(chǎng)地形，其高差可達(dá)1000 m，其中圓點(diǎn)為測(cè)風(fēng)塔，方點(diǎn)為機(jī)位點(diǎn)。

基于同一套的計(jì)算網(wǎng)格，分別采用不同湍流模型及模型參數(shù)開展16 扇區(qū)的定向計(jì)算（即對(duì)0～360°方向16 等分），分別進(jìn)行對(duì)應(yīng)風(fēng)向下的CFD 計(jì)算?；诙ㄏ蛴?jì)算結(jié)果、測(cè)風(fēng)數(shù)據(jù)和機(jī)型的功率曲線，采用文獻(xiàn)［25］中的發(fā)電量評(píng)估方法計(jì)算各機(jī)位點(diǎn)的風(fēng)速及發(fā)電量。在此基礎(chǔ)上，利用實(shí)際機(jī)型對(duì)應(yīng)的機(jī)艙傳遞函數(shù)反推實(shí)際來流風(fēng)速，并考慮發(fā)電量折減，本文采用發(fā)電量等效小時(shí)數(shù)He（下文簡(jiǎn)稱小時(shí)數(shù)）來評(píng)估發(fā)電量。

He = E／Wr（17）

式中：E——實(shí)際發(fā)電量，MWh；Wr——機(jī)型的額定功率，MW。

首先利用風(fēng)場(chǎng)內(nèi)測(cè)風(fēng)塔的實(shí)測(cè)數(shù)據(jù)驗(yàn)證不同模型得到的流場(chǎng)的合理性。測(cè)風(fēng)塔A 處的風(fēng)廓線對(duì)比如圖10 所示。可以看到，在高度較低時(shí)不同組的結(jié)果基本一致且與實(shí)測(cè)結(jié)果吻合較好；當(dāng)高度增加，不同組計(jì)算得到的風(fēng)廓線出現(xiàn)差異，其中G3 呈隨高度的增加風(fēng)速變小的趨勢(shì)，其他組基本保持不變。不同組的風(fēng)廓線相對(duì)誤差如表5 所示。k-l 模型經(jīng)參數(shù)修正后（G4）能有效提升對(duì)風(fēng)場(chǎng)中風(fēng)廓線預(yù)測(cè)的可靠性，相比原參數(shù)（G3）可降低誤差7.41%。

進(jìn)一步利用風(fēng)場(chǎng)運(yùn)行過程中機(jī)位點(diǎn)的實(shí)測(cè)數(shù)據(jù)，包括風(fēng)速和發(fā)電量，對(duì)比不同的計(jì)算模型并結(jié)合流場(chǎng)進(jìn)行分析。圖11 所示為不同湍流模型得到的風(fēng)加速因子和小時(shí)數(shù)分布與運(yùn)行實(shí)測(cè)結(jié)果的對(duì)比。可以看到風(fēng)加速因子與小時(shí)數(shù)的變化趨勢(shì)基本一致。然而，G3 的結(jié)果在部分點(diǎn)位（如L3、L22和L23 等）處顯著偏離實(shí)測(cè)結(jié)果，同時(shí)各點(diǎn)位之間的結(jié)果波動(dòng)較為明顯，而實(shí)測(cè)結(jié)果和其他模型均未出現(xiàn)該變化趨勢(shì)，說明k-l 模型的原始參數(shù)可能并不適用于大尺度的復(fù)雜地形流動(dòng)，會(huì)導(dǎo)致誤差增大。

為分析造成結(jié)果差異的原因，對(duì)比了不同模型獲得的流場(chǎng)。該風(fēng)電場(chǎng)的主風(fēng)向處于157.5°扇區(qū)內(nèi)，以L3 點(diǎn)位為例進(jìn)行分析。圖12 顯示了157.5°扇區(qū)下L3 點(diǎn)位附近的流場(chǎng)分布?？梢钥吹剑珿3 的流場(chǎng)與其他組存在較大差異，流經(jīng)第一個(gè)山坡（圖12c 中紅圈）后的尾流區(qū)域的范圍顯著大于其他組，因此L3 點(diǎn)位受尾流影響導(dǎo)致風(fēng)速降低，預(yù)測(cè)的小時(shí)數(shù)顯著低于其他模型和實(shí)測(cè)結(jié)果。這意味著采用原始參數(shù)的k-l 模型（G3）高估了山坡后尾流的影響，導(dǎo)致發(fā)電量預(yù)測(cè)偏低。除G3外，其他模型的流場(chǎng)基本一致，采用修正參數(shù)后的k-l 模型（G4）的流場(chǎng)的尾流范圍略小于k-ε 模型的結(jié)果（G1 和G2）。

表6 所示為不同模型的平均誤差，圖13 和圖14 分別展示了整體誤差分布?？梢钥吹剑瑹o論從平均結(jié)果還是整體分布上，G4 均表現(xiàn)出最低的誤差水平，其風(fēng)加速因子和小時(shí)數(shù)的平均誤差分別為6.3% 和10.9%。此外，G3 的誤差最大，分別達(dá)11.4% 和17.9%，說明原始參數(shù)的k-l 模型高估了尾流影響，導(dǎo)致預(yù)測(cè)誤差增大。上述結(jié)果表明，利用實(shí)際測(cè)風(fēng)數(shù)據(jù)修正參數(shù)能有效降低誤差，基于實(shí)測(cè)數(shù)據(jù)修正參數(shù)的k-l模型更適用于大尺度的復(fù)雜地形流動(dòng)。

3 結(jié) 論

本文利用實(shí)際測(cè)風(fēng)數(shù)據(jù)對(duì)k-l 湍流模型的參數(shù)進(jìn)行修正，使用Bolund 島和某實(shí)際風(fēng)電場(chǎng)對(duì)比了不同參數(shù)不同模型的影響，分析了流場(chǎng)差異。參數(shù)修正后的k-l 湍流模型能有效提升精度，說明采用k-l 湍流模型對(duì)實(shí)際風(fēng)電場(chǎng)進(jìn)行仿真時(shí)需通過該方法修正模型參數(shù)，主要結(jié)論如下：

1）相比參數(shù)優(yōu)化后的k-ε 模型，采用原始參數(shù)的k-l 模型的誤差偏大，風(fēng)加速因子和發(fā)電量等效小時(shí)數(shù)的平均相對(duì)誤差分別可達(dá)11.4% 和17.9%；基于實(shí)際測(cè)風(fēng)數(shù)據(jù)修正參數(shù)后能有效提升k-l 模型的仿真精度，風(fēng)加速因子和發(fā)電量等效小時(shí)數(shù)的平均誤差分別可降低5.1% 和7%。

2）采用原始參數(shù)的k-l 湍流模型會(huì)高估山坡后的尾流影響，從而導(dǎo)致預(yù)測(cè)風(fēng)速偏低，同時(shí)低估發(fā)電量；利用實(shí)際測(cè)風(fēng)數(shù)據(jù)修正模型參數(shù)后，可有效改善流場(chǎng)結(jié)構(gòu)、提高準(zhǔn)確性。

3） k-l 湍流模型的原始參數(shù)并不適用于大尺度的復(fù)雜地形流動(dòng)，在實(shí)際應(yīng)用過程中，應(yīng)采用當(dāng)?shù)販y(cè)風(fēng)數(shù)據(jù)修正模型參數(shù)，或考慮將參數(shù)B1 的量級(jí)調(diào)整為102。

［參考文獻(xiàn)］

［1］ BURTON T， SHARPE D， JENKINS N， et al. Windenergy handbook［M］. New York： John Wiley amp; SonsLtd.， 2001

［2］ 賀德馨. 風(fēng)工程與工業(yè)空氣動(dòng)力學(xué)［M］. 北京： 國(guó)防工業(yè)出版社， 2006： 221-225.

HE D X. Wind engineering and industrial aerodynamics［M］. Beijing： National Defense Industry Press， 2006：221-225.

［3］ DHUNNY A Z， LOLLCHUND M R， RUGHOOPUTH S DD V. Wind energy evaluation for a highly complex terrainusing Computational Fluid Dynamics（CFD）［J］.Renewable energy， 2017， 101： 1-9.

［4］ HAN X X， LIU D Y， XU C， et al. Similarity functions anda new k- ε closure for predicting stratified atmosphericsurface layer flows in complex terrain［J］. Renewableenergy， 2020， 150： 907-917.

［5］ 余文林， 柯世堂. 基于WRF與CFD嵌套的臺(tái)風(fēng)下大型風(fēng)力機(jī)流場(chǎng)作用與氣動(dòng)力分布［J］. 太陽能學(xué)報(bào)， 2020，41（12）： 260-269.

YU W L， KE S T. Flow field action and aerodynamic loadsdistribution for large-scale wind turbine under typhoonbased on nesting of WRF and CFD［J］. Acta energiaesolaris sinica， 2020， 41（12）： 260-269.

［6］ 趙子涵， 李朝， 肖儀清， 等. 基于NWP/CFD嵌套的復(fù)雜地形風(fēng)場(chǎng)模擬研究［J］. 太陽能學(xué)報(bào)， 2021， 42（2）： 205-210.

ZHAO Z H， LI C， XIAO Y Q， et al. Wind field simulation over complex terrain by coupling NWP/CFD approach［J］.Acta energiae solaris sinica， 2021， 42（2）： 205-210.

［7］ 侯亞麗， 呂愛靜， 邸建琛， 等. 建筑密度對(duì)建筑物群內(nèi)風(fēng)能利用的影響［J］. 太陽能學(xué)報(bào)， 2022， 43（5）： 336-342.

HOU Y L， LYU A J， DI J C， et al. Influence of buildingdensity on wind energy utilization in buildings［J］. Actaenergiae solaris sinica， 2022， 43（5）： 336-342.

［8］ BELJAARS A C M， WALMSLEY J L， TAYLOR P A. Amixed spectral finite-difference model for neutrallystratified boundary-layer flow over roughness changes andtopography［J］. Boundary-layer meteorology， 1987， 38（3）： 273-303.

［9］ RICHARDS P J， HOXEY R P. Appropriate boundaryconditions for computational wind engineering modelsusing the k-ε turbulence model［J］. Journal of windengineering and industrial aerodynamics， 1993， 46/47：145-153.

［10］ 劉文， 王曉東， 閆姝， 等. 基于校準(zhǔn)k-ε 模型的復(fù)雜地形流動(dòng)數(shù)值模擬［J］. 工程熱物理學(xué)報(bào)， 2021， 42（2）： 377-385.

LIU W， WANG X D， YAN S， et al. Numerical simulations offlow within complex terrain based on calibrated k- ε model［J］. Journal of engineering thermophysics， 2021， 42（2）：377-385.

［11］ MELLOR G L， YAMADA T. Development of a turbulenceclosure model for geophysical flow problems［J］. Review ofgeophysics and space physics， 1982， 20：831-875.

［12］ YAMADA T. Simulations of nocturnal drainage flows by aq2l turbulence model［J］. Journal of the atmosphericsciences， 1983， 40：91-106.

［13］ ARRITT R W. The effect of water surface temperature onlake breezes and thermal internal boundary layers［J］.Boundary-layer meteorology， 1987， 40（1）： 101-125.

［14］ MANNING J， HANCOCK P， WHITING R. A study of theability of Meteodyn WT to replicate measurements aroundsteep hills using wind tunnel data from the ‘RUSHIL’experiment［J］. Wind engineering， 2010， 34（5）： 477-499.

［15］ LAUFER J. The structure of turbulence in fully developedpipe flow［ R］. NACA-TR-1174， 1954.

［16］ LAUFER J. Investigation of turbulent flow in a twodimensionalchannel［R］. NACA-TR-1053， 1951.

［17］ KLEBANOFF P S. Characteristics of turbulence in aboundary layer with zero pressure gradient［ R］. NACA-TN-3178， 1954.

［18］ ROSE W G. Results of an attempt to generate ahomogeneous turbulent shear flow［J］. Journal of fluidmechanics， 1966， 25： 97-120.

［19］ PATANKAR S V， SPALDING D B. A calculationprocedure for heat， mass and momentum transfer in threedimensionalparabolic flows［J］. International journal ofheat and mass transfer， 1972， 15（10）： 1787-1806.

［20］ BLACKADAR A K. The vertical distribution of wind andturbulent exchange in a neutral atmosphere［J］. Journal ofgeophysical research， 1962， 67（8）： 3095-3102.

［21］ LOUIS J F. A parametric model of vertical eddy fluxes inthe atmosphere［J］. Boundary-layer meteorology， 1979， 17（2）： 187-202.

［22］ YAMADA T. An application of a three-dimensional，simplified second-moment closure numerical model tostudy atmospheric effects of a large cooling-pond［J］.Atmospheric environment， 1967， 13（5）： 693-704.

［23］ BERG J， MANN J， BECHMANN A， et al. The bolundexperiment， part I： flow over a steep， three-dimensionalhill［J］. Boundary-layer meteorology， 2011， 141（2）： 219-243.

［24］ BECHMANN A， S?RENSEN N N， BERG J， et al. Thebolund experiment， part II： blind comparison ofmicroscale flow models［J］. Boundary-layer meteorology，2011， 141（2）： 245-271.

［25］ 姜婷婷， 葉杭冶， 申新賀， 等. 基于風(fēng)輪面等效風(fēng)速的風(fēng)電場(chǎng)發(fā)電量評(píng)估方法研究［J］. 太陽能學(xué)報(bào)， 2021， 42（9）： 244-249.

JIANG T T， YE H Y， SHEN X H， et al. Research on windfarm power generation assessment based on rotorequivalent wind speed［J］. Acta energiae solaris sinica，2021， 42（9）： 244-249.

基金項(xiàng)目：國(guó)家重點(diǎn)研發(fā)計(jì)劃（2018YFB1501100）