董依帆， 柯世堂， 楊 青

(南京航空航天大學 土木工程系，南京 210016)

近年來隨著國內(nèi)能源結(jié)構(gòu)轉(zhuǎn)型步伐的加快，風力發(fā)電設(shè)施陸續(xù)興建。其中我國西部風能資源充沛[1]，為更好的捕捉風能，風力機必然朝著高大化方向發(fā)展。然而此地區(qū)戈壁沙漠面積廣大，沙粒受到近地邊界層風場的驅(qū)動，在一定高度范圍內(nèi)隨機運動，兩者之間相互耦合并發(fā)生能量遷移，使得沙粒以較大速度沖擊至風力機表面，形成附加荷載及風蝕效應(yīng)，對風力機的安全性能和使用壽命產(chǎn)生顯著影響[2-3]。因此，開展氣-固兩相流理論研究，探討不同風速和沙粒粒徑組合對風力機氣動載荷的作用機理，具有重要的工程意義和理論價值。

當前國內(nèi)外對風-沙共同作用的研究多集中在風沙現(xiàn)象和風沙運動本身，包括沙粒運動力學過程[4]、輸沙量的垂直分布[5]、輸沙率模型[6]、臨界起沙風速[7]和風速廓線規(guī)律[8]。相較于風荷載作用，風-沙共同作用下建筑物或構(gòu)筑物荷載效應(yīng)更加顯著[9]。文獻[10]通過風洞測力試驗，研究了不同沙濃度和風速條件下風沙對低矮建筑整體受力影響，主要考察風沙環(huán)境中建筑物表面的風壓特性；文獻[11]模擬了大氣邊界層中單體及兩相鄰對稱高層建筑的風沙流場，并指出隨著沙粒相體積分數(shù)的增大，風壓系數(shù)不斷增大；文獻[12]對風沙環(huán)境下水平軸風力機進行數(shù)值模擬，研究了在某一來流風速和不同沙粒直徑組合工況下的風力機風輪轉(zhuǎn)矩。已有研究均未基于風力機整機和沙粒真實環(huán)境分布特性展開多參數(shù)條件下大型風力機風沙特性的系統(tǒng)研究。

鑒此，本文以南京航空航天大學自主研發(fā)的5 MW水平軸三葉片風力機[13]為研究對象，采用改進的k-ε湍流模型和粉塵釋放模型(Dust Production Model，DPM)，設(shè)置不同風速和不同沙粒粒徑形成9種組合計算工況，對處于典型停機狀態(tài)的塔架-葉片耦合體系展開數(shù)值模擬，以獲取沙粒在風力機不同部位的分布密度及撞擊速度，并對比分析風沙共同作用下風力機體系繞流特性、表面壓力分布隨風沙流速度和沙粒濃度的變化規(guī)律。

1 風-沙雙相耦合算法

1.1 沙塵暴等級劃分

我國西北地區(qū)沙塵暴天氣頻發(fā)，主要集中在春季4～5月，其發(fā)生時間集中、頻率高、強度大。根據(jù)探空資料，沙塵暴發(fā)生時形成的沙塵壁高達300 m，沙塵暴影響高度在2 100 m以內(nèi)[14-15]。沙塵暴天氣根據(jù)最大風速與最小能見度劃分為4個等級，如表1所示。本文選取中、強、特強沙塵暴天氣典型風速18 m/s，22 m/s，25 m/s作為風場計算風速。

表1 沙塵暴天氣強度劃分標準

1.2 風沙流密度

風沙流密度是指空氣運動速度在大于起沙風速情況下，單位體積空氣氣流中所含沙的質(zhì)量，單位為kg/m3。在近地層中，風受地面摩擦阻力的影響而速度降低，一般而言，摩擦力隨著高度的增加而減小，故風速隨著高度的增加而增大，風沙流密度與高度、風力等均有密切關(guān)系。

據(jù)文獻[16]統(tǒng)計，當風速在40 m/s以下時，風沙流密度的數(shù)量級一般約為10-5，且隨高度的增加而下降，基本以3 m高度作為分界點。在3 m高度以下，風沙流密度隨高度的增加而大幅下降，在3 m高度以上風沙流密度隨高度增加逐步下降，且沙粒運動以飛揚為主?；诠こ贪踩钥紤]，采取極端工況進行數(shù)值模擬，故風場入口處3 m以下高度風沙流密度取為10-4kg/m3，速度與風速相同；3 m以上高度風沙流密度取為10-5kg/m3，速度為0。

1.3 計算模型與顆粒運動方程

標準k-ε模型在科學研究及工程實際中得到了廣泛的應(yīng)用，但用于強旋流、彎曲壁面流動或彎曲流線流動時，會產(chǎn)生一定的失真。為彌補其缺陷，本文采用Realizablek-ε模型，Realizablek-ε模型中k與ε的模數(shù)化輸運方程組為

(1)

(2)

與標準k-ε模型比較發(fā)現(xiàn)，Realizablek-ε模型主要變化是：湍流黏度計算公式發(fā)生了變化，引入了與旋轉(zhuǎn)和曲率有關(guān)的內(nèi)容；ε方程發(fā)生了很大變化，方程中的產(chǎn)生項不再包含有k方程中的產(chǎn)生項Gk；ε的倒數(shù)第二項不具有任何奇異性，即使值很小或為零，也不會為零。

相比于視為連續(xù)相的空氣來流，沙粒在空氣中所占比例較小，即使在特強沙塵暴(風速≥25 m/s)時，沙粒占空氣體積分數(shù)仍遠小于10%。基于此類實際情況，DPM模型被廣泛應(yīng)用于多相流聯(lián)合研究。本文在風場計算穩(wěn)定后，將離散相模型作為第二相，介入至連續(xù)相中展開風-沙雙向耦合運算。沙粒在風場中運動平衡方程為

(3)

(4)

式中：μ為流體黏性系數(shù)；dp為顆粒直徑；Re為相對雷諾數(shù)； 可表示為

(5)

考慮沙粒離散相影響后，風連續(xù)相基本控制方程可表示為

(6)

(7)

(8)

式中：I為單位張量； 等式右邊第二項為體積膨脹作用。

1.4 顆粒碰撞模型

沙粒沖擊到風力機表面過程服從動量守恒定律，求解沖擊力的關(guān)鍵在于碰撞時間。計算中忽略沙粒在沖擊過程中可能發(fā)生的破裂現(xiàn)象，認為沙粒與結(jié)構(gòu)間相互作用遵循牛頓第二定律并假設(shè)反彈后的速度與撞擊前的速度一致。由動量定理

(9)

式中：f(t)為單個沙粒沖擊力矢量， N；vs為沙粒速度矢量。

沙粒在單位時間內(nèi)對結(jié)構(gòu)的沖擊力F(τ)為

(10)

將撞擊風力機的沙粒近似看作球體，則

(11)

由于沙粒直徑在1 mm及以下，撞擊前水平末速度相對較大，可將風力機表面簡化為無限大平面[17]，將碰撞時間τ取為

(12)

(13)

式中：m為沙粒的質(zhì)量；φ為速度恢復系數(shù)取1；k1，k2取值與沙粒、風力機材料楊氏模量、泊松比有關(guān)。

則沙粒對結(jié)構(gòu)的沖擊力可簡化為

(14)

2 數(shù)值模擬

2.1 工程簡介

表2給出了該5 MW風力機的主要結(jié)構(gòu)尺寸及示意圖，葉片長度為60 m，機艙尺寸18 m×6 m×6 m，各葉片沿周向成120°夾角分布，塔架通長變厚，高為124 m，塔底半徑3.5 m，塔頂半徑3.0 m，塔頂厚度90 mm，塔底厚度200 mm。

表2 5 MW風力發(fā)電結(jié)構(gòu)體系主要參數(shù)列表

2.2 工況設(shè)置

該風力機位于B類地貌，綜合考慮我國沙塵暴天氣等級劃分，對比研究3種風速和3種沙粒粒徑組合下極端工況對風力機整機的作用規(guī)律。其中，小風、中風和大風分別以沙塵暴等級中、強和特強的基本風速進行劃分；沙粒粒徑分別取為最大統(tǒng)計概率區(qū)間內(nèi)的0.10 mm，0.25 mm以及極端工況下的1 mm，3種風速和沙粒粒徑組合共計9種工況，如表3所示。

表3 風力機風-沙耦合計算工況劃分表

2.3 計算域設(shè)置

如圖1所示，本文研究的風力機均處于0°來流風向角。為保證風力機尾流能夠充分發(fā)展，計算域尺寸設(shè)置為12D×5D×5D(流向X×展向Y×豎向Z，D為風輪直徑)，風力機置于距離計算域入口3D處。為兼顧計算效率與精度，同時考慮到葉片表面扭曲復雜，網(wǎng)格劃分采用混合網(wǎng)格離散形式，將整個計算域劃分為局部加密區(qū)域和外圍區(qū)域。局部加密區(qū)域內(nèi)含風力機模型，采用非結(jié)構(gòu)化網(wǎng)格進行劃分，外圍區(qū)域形狀規(guī)整，采用高質(zhì)量的結(jié)構(gòu)化網(wǎng)格進行劃分。

為確保數(shù)值模擬結(jié)果的可信度，本文增加了網(wǎng)格無關(guān)性驗證，表4給出了不同網(wǎng)格方案下網(wǎng)格質(zhì)量和迎風面壓力系數(shù)。由表4可知，隨著網(wǎng)格總數(shù)的增加，網(wǎng)格質(zhì)量逐漸增加，網(wǎng)格歪斜度和迎風面風壓系數(shù)呈現(xiàn)逐漸減小的趨勢，而840萬網(wǎng)格數(shù)和1 100萬網(wǎng)格數(shù)的網(wǎng)格質(zhì)量和計算結(jié)果無明顯差異，綜合計算精度和效率，本文選取840萬網(wǎng)格總數(shù)的方案。計算域及具體的網(wǎng)格劃分如圖1所示。

表4 不同網(wǎng)格方案下網(wǎng)格質(zhì)量和迎風面壓力系數(shù)

圖1 計算域及網(wǎng)格劃分示意圖Fig.1 Schematic diagram of computational domain and encrypted mesh

計算域左側(cè)和頂部邊界條件為速度入口，右側(cè)為壓力出口，數(shù)值計算采用3D單精度、分離式求解器，流場流速為絕對速度，空氣模型等效為理想不可壓縮流體，計算域入口采用冪指數(shù)為0.15的風廓線模型，離地10 m高度處的風速分別設(shè)置為“2.2”節(jié)中3種基準風速。流場求解采用SIMPLEC算法實現(xiàn)速度與壓力之間的耦合，對流項求解格式為二階，計算過程中設(shè)置了網(wǎng)格傾斜校正以提高混合網(wǎng)格計算效果，控制方程的計算殘差設(shè)置為1×10-6，最后初始化風場進行迭代計算。圖2給出了平均風速、湍流度剖面模擬結(jié)果與理論值的對比曲線，結(jié)果表明平均風速和湍流度剖面均與理論值吻合良好，風場模擬標準滿足工程要求。

圖2 不同風速下風場速度及湍流度剖面示意圖Fig.2 Velocity and turbulence profiles under different wind speeds

2.4 有效性驗證

圖3給出了3種風速下塔架30 m高度環(huán)向壓力系數(shù)分布圖。對比可知各工況平均風壓系數(shù)分布曲線的負壓極值點和分離點對應(yīng)角度與規(guī)范曲線[18]一致，迎風和側(cè)風區(qū)域風壓系數(shù)數(shù)值吻合較好，僅在背風區(qū)負壓系數(shù)略大于規(guī)范值，對比結(jié)果驗證了風場數(shù)值模擬的有效性。

圖3 3種風速下風力機塔架30 m高度壓力系數(shù)與規(guī)范對比曲線Fig.3 Pressure coefficients on the 30 m height tower under three wind speeds and domestic codes

3 結(jié)果分析

為研究風-沙荷載特征值隨環(huán)境參數(shù)變化的作用規(guī)律，首先分析風-沙雙向耦合作用下風場分布特性、沙粒著點及其分布規(guī)律；在此基礎(chǔ)上，通過計算沙粒沖擊荷載、等效外壓系數(shù)定量分析不同工況下沙致效應(yīng)差異。

3.1 風場分布特性

圖4給出了不同風速下大型風力機體系風速流線圖。由圖4可知，由于葉片的遮擋效應(yīng)，塔架迎風面壓力明顯減小，葉片與塔架重合部位內(nèi)側(cè)出現(xiàn)附著于塔架迎風面的自由剪切層，并在塔架背風面形成成熟的旋渦脫落。

圖4 風力機體系風速流線圖Fig.4 Wind speed streamline of wind turbine

圖5和圖6分別給出了未干擾及干擾區(qū)段塔架渦量云圖。對比分析可知：未干擾區(qū)段氣流在風力機塔架迎風面前緣發(fā)生流動分離，并在側(cè)邊出現(xiàn)加速效應(yīng)；而在干擾區(qū)段由于葉片對塔架的干擾作用，塔架和葉片之間的區(qū)域隨著風速的增加出現(xiàn)明顯的能量積聚，持續(xù)發(fā)展后在背風面形成尾流渦旋以及回流，且隨風速的增加而加強。

圖5 風力機塔架未干擾段湍動能圖Fig.5 Wind turbine tower undisturbed turbulent kinetic energy

圖6 風力機塔架干擾段湍動能圖Fig.6 Turbulence kinetic energy of wind turbine tower interference section

3.2 風沙運動特征

風沙流場中沙粒沖擊風力機表面，易產(chǎn)生極端荷載效應(yīng)。為清晰展示沙粒在風力機不同部位的分布密度，基于顆粒水平向速度對沙粒軌跡進行追蹤，同時定義單位時間內(nèi)沖擊塔架沙粒的位置分布為沙粒著點，圖7給出了9種工況風沙場中沙粒著點圖，并對沙粒的密集程度進行了等比例粗化處理。

由圖7可知，由于葉片的對塔架的遮擋效應(yīng)，各工況沙粒撞擊位置多集中分布風力機塔架迎風區(qū)域0～0.6H高度范圍內(nèi)，H為風力機塔架高度。受氣流旋渦驅(qū)動作用，塔架遮擋區(qū)和背風區(qū)壁面僅有少量沙粒附著；塔架、葉片收集到的沙粒均以工況7最多，且隨風速的增加而增加，隨沙粒粒徑的增大逐漸變小。

圖7 風沙場中沙粒著點示意圖Fig.7 Sand-reflect motion trajectory in sand and wind field

3.3 風沙荷載分布

3.3.1 沙?？臻g分布規(guī)律

圖8和圖9分別給出了9種工況下風力機塔架和葉片收集到的沙粒數(shù)量、撞擊速度及速度占有率(各速度分布區(qū)間內(nèi)沙粒數(shù)量與總數(shù)量的比值)對比曲線。由圖8和圖9對比可知：

(1) 當沙粒粒徑較小時，風力機塔架收集到的沙粒數(shù)量隨風速的增大而增加，當粒徑達到中等以上時，風力機塔架收集到的沙粒數(shù)量隨風速的增大而減少；然而，粒徑和風速對風力機葉片的沙粒著點影響并不顯著，這是由于葉片位置較高，處于風沙流密度影響較小區(qū)域。

圖8 風力機塔架沙粒數(shù)量與水平末速度分布曲線Fig.8 Wind turbine tower sand number and horizontal velocity distribution curve

圖9 風力機葉片沙粒數(shù)量與水平末速度分布曲線Fig.9 Wind turbine blade grit and horizontal velocity distribution curve

(2) 撞擊風力機塔架的沙粒，隨著風速增大，小粒徑和大粒徑沙粒的平均速度均增大，而中等粒徑的沙粒平均速度隨風速的增大而減?。伙L力機葉片上追蹤的沙粒平均速度特性也呈現(xiàn)相似規(guī)律，不同在于風力機葉片上的大粒徑沙粒平均速度隨風速的增大而減小。

(3) 風力機塔架、葉片收集到的沙粒水平速度占有率分布規(guī)律基本一致，速度占有率隨著水平末速度(順風向平均速度)的增大先增加后減小。此外，收集到的沙粒平均水平速度均遠小于基準風速，且沙粒水平撞擊速度均隨著沙粒直徑的增加而增大。

3.3.2 葉片、塔架沙粒沖擊荷載

定義與塔架重合的葉片為A、迎風面左側(cè)葉片為B、右側(cè)葉片為C。采用式(12)分別進行9種工況下風力機塔架、葉片沙荷載計算，同時給出了塔架、葉片不同高度(葉展長度)范圍沙荷載以及風力機不同高度范圍沙荷載與該區(qū)域風荷載的比值，如圖10、圖11和表5所示。經(jīng)對比發(fā)現(xiàn)：

(1) 大粒徑沙粒與不同風速的組合工況(工況3、工況6、工況9)在風力機表面產(chǎn)生顯著荷載，由于位置和形狀的不同，風力機塔架、葉片沙荷載分布呈現(xiàn)不同的規(guī)律；

(2) 葉片A、葉片B、葉片C沙荷載均沿葉展方向呈減小趨勢，由于所處高度不同，風沙流密度影響程度不同，葉片A承受最大沙荷載約為葉片B、葉片C的兩倍；

(3) 不同工況對風力機塔架表面沙荷載分布影響較大，但均在0～0.05H高度范圍內(nèi)產(chǎn)生該工況下極大沙荷載，工況3下該高度范圍內(nèi)沙荷載與風荷載比值達22.57%；

(4) 顯著工況(工況3、工況6、工況9)下風力機塔架極大沙荷載隨風速的變大而增大，沙荷載沿塔架高度先減小后增大再減小，0.2H～0.6H高度范圍內(nèi)沙荷載極大值也隨風速的增大而增大。由于葉片A的遮擋效應(yīng)，塔架0.5H高度以上沙荷載減小。

圖10 各工況風力機葉片沙荷載示意圖Fig.10 Sand load distribution on different height of wind turbine blades

圖11 風力機塔筒沿子午線高度荷載示意圖Fig.11 Comparison curves of equivalent pressure in typical meridian of wind turbine tower

表5 工況1～工況9塔架不同高度范圍風沙荷載比值列表

3.3.3 沙壓系數(shù)

沙致壓力系數(shù)(簡稱沙壓系數(shù))計算方法為：①將表面各監(jiān)控點沙粒撞擊荷載轉(zhuǎn)化成沙致壓強；②計算監(jiān)控點沙致壓強與對應(yīng)參考高度處風壓比值，即沙壓系數(shù)。圖12給出了各工況下風力機塔架沙壓系數(shù)三維分布示意圖。由圖12可知：

(1) 各工況風力機塔架沙壓系數(shù)均主要集中于迎風面兩側(cè)各60°范圍內(nèi)，其余范圍數(shù)值基本為0，各工況沙壓系數(shù)最大值均發(fā)生在0～0.05H高度范圍內(nèi)，各工況沙壓系數(shù)最大值為0.517；

(2) 沙壓系數(shù)分布呈現(xiàn)明顯三維效應(yīng)，其隨塔架高度增加而迅速減小，塔架與葉片重合部分沙壓系數(shù)趨于0；

(3) 不同組合工況下沙壓系數(shù)大小差異明顯，且沙粒粒徑的增大而顯著增大，此外，其最大值隨風速的增大呈減小趨勢，而其總體分布更趨于平均化。

圖12 各工況風力機塔架沙致壓力系數(shù)三維分布示意圖Fig.12 Three-dimensional distribution of sand-reflect and sand pressure coefficient on surface of wind turbine tower

3.3.4 風-沙共同作用等效壓力系數(shù)

為定量比較不同組合工況下的塔架風-沙致壓力分布，將沙壓系數(shù)與風壓系數(shù)矢量加和，定義等效壓力系數(shù)作為沙致效應(yīng)的衡量參數(shù)。

考慮沙粒著點分布特性及葉片與塔筒之間的氣動干擾效應(yīng)，選取塔架未干擾區(qū)段沙粒分布密度差異較大，干擾區(qū)段塔架與葉片重合面積不同的兩個斷面為典型斷面。圖13給出各工況下風力機塔架4個典型斷面等效壓力系數(shù)對比曲線，分析可知：

圖13 風力機塔架典型斷面環(huán)向等效壓力系數(shù)對比曲線Fig.13 Comparison curves of equivalent pressure coefficient in typical section of wind turbine tower

(1) 不同工況下相同高度截面處等效壓力系數(shù)分布規(guī)律及數(shù)值基本一致，位于未干擾區(qū)(0.025H，0.30H)塔筒斷面環(huán)向風壓系數(shù)呈現(xiàn)良好的對稱性；而位于顯著干擾區(qū)(0.70H，0.90H)塔筒斷面，其環(huán)向表面壓力系數(shù)分布曲線不再保持對稱，塔筒迎風面中心點處呈現(xiàn)負壓；

(2) 同一工況下4個典型斷面等效壓力系數(shù)數(shù)值均略有差異，其中最大負壓隨著高度的增加先增大后減小，在葉片遮擋區(qū)域附近約為-1.5；

(3) 不同工況下0.025H高度塔筒斷面迎風面兩側(cè)0～30°內(nèi)壓力系數(shù)差異明顯，工況3、工況6、工況9大于其余工況下等效壓力系數(shù)，且最大壓力系數(shù)隨風速的增大而減小，背風面呈負壓；

(4) 位于干擾區(qū)段塔架斷面，隨高度增加迎風面、側(cè)風面和背風面負壓分別呈增大、減小和增大的趨勢。

圖14給出了風力機塔架0°，90°，180°及270°四條典型子午線的等效壓力系數(shù)對比曲線，分析可得：

(1) 風力機塔架90°與270°子午線等效壓力系數(shù)隨高度先增大后減小再增大，均在0.5H，0.8H高度處出現(xiàn)拐點；0°和180°子午線迎風面等效壓力系數(shù)發(fā)生數(shù)值分離多集中在沙粒著點數(shù)目較多的0.5H高度以下的未干擾區(qū)段，0.1H～0.6H高度范圍內(nèi)背風面等效壓力系數(shù)出現(xiàn)明顯差異；

(2) 90°和270°子午線等效壓力系數(shù)最大值均產(chǎn)生在0.5H高度處，最大值約為為-1.890；受風沙流影響，0°子午線最大等效壓力系數(shù)發(fā)生0.005H高度處，約為1.092；270°子午線方向最大等效壓力系數(shù)發(fā)生在塔架頂部附近，最大值為-0.282；

(3) 0°子午線底部等效壓力系數(shù)受風沙流影響顯著，大粒徑沙粒對塔架的正面沖擊作用加強，沙粒著點隨風速的增大而減少，故等效壓力系數(shù)隨著風速的增大而減小。

圖14 風力機塔筒典型子午線等效壓力系數(shù)對比曲線Fig.14 Comparison curves of equivalent pressure coefficient in typical meridians of wind turbine tower

4 結(jié) 論

(1) 僅考慮單向風作用時風力機體系風場和風壓系數(shù)均呈現(xiàn)明顯的三維效應(yīng)，考慮風-沙雙向耦合后沙粒主要沖擊塔架迎風面中下部，僅有少量沙粒隨氣流進入塔筒背風面；沙粒水平方向作用力與分布范圍大小受風速和粒徑的影響，在塔架和葉片等不同部位呈現(xiàn)不同的規(guī)律。

(2) 風力機塔架、葉片沙粒沖擊數(shù)量隨風速的增大而增加，隨粒徑的增大而減少，沙粒統(tǒng)計數(shù)量最多的為工況7(特強沙塵暴發(fā)生風速+典型工況下較小沙粒粒徑)。塔架上的沙粒速度占有率最大值隨風速、粒徑的增大而降低，速度分布更加均勻，而沖擊葉片沙粒速度占有率分布規(guī)律不受風速、粒徑等因素干擾，主要表現(xiàn)為峰值右移。

(3) 不同工況下沙粒撞擊位置主要集中在塔架0.6H高度以下，受兩側(cè)流動分離渦裹挾，集中于迎風區(qū)域兩側(cè)各60°范圍，且大粒徑沙粒沖擊塔架表面更易形成顯著荷載，沙致壓力系數(shù)最大值為0.517，發(fā)生在工況3(中等沙塵暴發(fā)生風速+極端工況下沙粒粒徑)的塔架底部迎風面。

(4) 不同工況下風力機塔架最大負壓值隨著高度的增加先增大后減小，背風區(qū)域負壓值隨著高度的增加先減小后增大，在葉片遮擋區(qū)域附近約為-1.5。各子午線等效壓力系數(shù)隨著高度的增加變化趨勢不同，0°子午線底部等效壓力系數(shù)受風沙流影響顯著，大粒徑沙粒對塔架的正面沖擊作用加強，最大等效壓力系數(shù)發(fā)生0.005H高度處，數(shù)值為1.092。