董依帆， 柯世堂， 楊 青

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

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

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

鑒此，本文以南京航空航天大學自主研發(fā)的5 MW水平軸三葉片風力機[13]為研究對象，采用改進的k-ε湍流模型和粉塵釋放模型(Dust Production Model，DPM)，設置不同風速和不同沙粒粒徑形成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 風沙流密度

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

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

(1)

(2)

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

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

(3)

(4)

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

(5)

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

(6)

(7)

(8)

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

1.4 顆粒碰撞模型

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

(9)

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

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

(10)

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

(11)

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

(12)

(13)

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

則沙粒對結構的沖擊力可簡化為

(14)

2 數(shù)值模擬

2.1 工程簡介

表2給出了該5 MW風力機的主要結構尺寸及示意圖，葉片長度為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ā)電結構體系主要參數(shù)列表

2.2 工況設置

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

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

2.3 計算域設置

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

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

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

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

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

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

2.4 有效性驗證

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

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

3.1 風場分布特性

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

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

圖5和圖6分別給出了未干擾及干擾區(qū)段塔架渦量云圖。對比分析可知：未干擾區(qū)段氣流在風力機塔架迎風面前緣發(fā)生流動分離，并在側邊出現(xiàn)加速效應；而在干擾區(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)生極端荷載效應。為清晰展示沙粒在風力機不同部位的分布密度，基于顆粒水平向速度對沙粒軌跡進行追蹤，同時定義單位時間內(nèi)沖擊塔架沙粒的位置分布為沙粒著點，圖7給出了9種工況風沙場中沙粒著點圖，并對沙粒的密集程度進行了等比例粗化處理。

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

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

3.3 風沙荷載分布

3.3.1 沙粒空間分布規(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) 撞擊風力機塔架的沙粒，隨著風速增大，小粒徑和大粒徑沙粒的平均速度均增大，而中等粒徑的沙粒平均速度隨風速的增大而減??；風力機葉片上追蹤的沙粒平均速度特性也呈現(xiàn)相似規(guī)律，不同在于風力機葉片上的大粒徑沙粒平均速度隨風速的增大而減小。

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

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

定義與塔架重合的葉片為A、迎風面左側葉片為B、右側葉片為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的遮擋效應，塔架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)控點沙粒撞擊荷載轉化成沙致壓強；②計算監(jiān)控點沙致壓強與對應參考高度處風壓比值，即沙壓系數(shù)。圖12給出了各工況下風力機塔架沙壓系數(shù)三維分布示意圖。由圖12可知：

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

(2) 沙壓系數(shù)分布呈現(xiàn)明顯三維效應，其隨塔架高度增加而迅速減小，塔架與葉片重合部分沙壓系數(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ù)作為沙致效應的衡量參數(shù)。

考慮沙粒著點分布特性及葉片與塔筒之間的氣動干擾效應，選取塔架未干擾區(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高度塔筒斷面迎風面兩側0～30°內(nèi)壓力系數(shù)差異明顯，工況3、工況6、工況9大于其余工況下等效壓力系數(shù)，且最大壓力系數(shù)隨風速的增大而減小，背風面呈負壓；

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

圖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 結 論

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

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

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

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