考慮葉片停機位置大型風力機塔架風-沙致結構響應分析

柯世堂　董依帆

摘要： 強風停機狀態下葉片位置會顯著影響風力機塔架的繞流及穩定性能，尤其在沙塵暴極端天氣條件下，沙粒的附著也會影響風的湍流特征并對塔架產生附加沖擊力，現有工作均缺乏風?沙耦合作用對大型風力機體系氣動性能及結構響應的探索。以南京航空航天大學自主研發的5 MW水平軸風力機體系為對象，以風?沙雙向耦合算法為核心，基于CFD技術分別采用連續相和離散相模型進行典型風場和沙粒組合的同步迭代模擬。考慮葉片單個旋轉周期內8個停機位置，對比分析了風力機塔架表面等效壓力系數和氣動載荷分布特性;結合有限元方法系統探討了不同停機位置下風力機體系動力特性、風?沙致結構響應和屈曲穩定性能。研究表明：風?沙耦合環境下沙粒在風力機塔架迎風面底部產生的荷載效應最為顯著，部分區域沙致壓力系數可達0.55，沙粒沖擊荷載與風荷載的比值最高可達23.7%;不同停機位下風荷載均對塔架結構響應起主導作用，但沙粒沖擊作用在塔架底部的放大效應不可忽視。

關鍵詞： 大型風力機體系; 風?沙耦合環境; 結構響應; CFD數值模擬; 停機位置

中圖分類號： TK83??? 文獻標志碼： A??? 文章編號： 1004-4523（2021）01-0060-12

DOI：10.16385/j.cnki.issn.1004-4523.2021.01.007

引? 言

中國西部風能資源充沛，然而此地區戈壁沙漠面積廣大，沙粒受到近地邊界層風場的驅動，在一定高度范圍內隨機運動^[1?2];同時，沙粒的附著也會影響風的湍流特征并對風力機體系產生附加沖擊力，兩者之間相互耦合并發生能量遷移^[3]，使得沙粒以較大速度沖擊至風力機表面，形成附加荷載及風蝕效應^[4]。在沙塵暴極端天氣條件下，風力機組一般處于停機狀態，此時葉片不同停機位置將顯著影響整個風力機塔架?葉片耦合體系的氣動力分布和結構穩定性能^[5]。因此，對風?沙耦合環境下大型風力機體系氣動性能及結構響應開展探索性研究具有重要的工程與理論意義。

大型風力機體系屬于典型的風敏感結構，現有研究主要集中在流場作用機理、結構氣動力分布^[6]、氣彈耦合效應^[7]和體系風致響應特性^[8]等內容。針對風?沙環境下風力機的氣動性能研究較少，文獻[9]通過低角度加速沖蝕磨損試驗研究了風沙環境下不同速度、沖角對風力機葉片的沖刷磨損行為特性，分析發現葉片涂層的沖蝕磨損量隨著沖蝕速度的增大而增加;文獻[10]采用虛擬風洞的方法研究了風沙環境下模型風輪的轉矩特性，獲得了風力機風輪轉矩隨沙塵顆粒粒徑和體積分數變化的規律，分析發現，風輪轉矩隨沙塵體積分數的增大先增大后減小，隨沙塵顆粒直徑的增大，先增大后減小，最后趨于平穩;綜上所述，風?沙環境對于風力機氣動性能的影響不可忽略，尤其對于沙塵暴極端天氣下大型風力機體系的風?沙耦合作用效應缺乏定量和定性的研究。此外，國內外風電機組設計規范和標準^[11?15]也尚未考慮此類極端環境下風力機的氣動和結構性能評估。

鑒于此，本文以南京航空航天大學自主研發的5 MW水平軸三葉片風力機為研究對象，基于風?沙雙向耦合算法，分別采用連續相和離散相模型進行典型風場和沙粒組合的同步迭代模擬。在此基礎上主要研究：（1）不同停機位置下風力機塔架表面沙致壓力系數和氣動載荷分布特性;（2）不同停機位置下風力機體系動力特性、風?沙致結構響應和屈曲穩定性能。研究結論可為此類極端天氣條件下大型風力機風?沙沖擊荷載取值和結構設計提供參考。

1 風-沙雙相耦合算法

1.1 沙塵暴等級劃分與風沙流密度

中國西北地區沙塵暴天氣頻發，主要集中在春季4至5月，其發生時間集中、頻率高、強度大。根據探空資料，沙塵暴發生時形成的沙塵壁高達300 m，沙塵暴影響高度在2100 m以內^[16?17]。沙塵暴天氣根據最大風速和最小能見度劃分為四個等級，如表1所示。本文選取中級沙塵暴天氣典型風速18 m/s風場計算風速，雙相耦合迭代算法中沙粒粒徑取為1 mm^[18]。

2 數值模擬

2.1 工況設置

表2給出了南京航空航天大學自主研發的5 MW上風型水平軸三葉片風力機主要結構設計參數及模型示意圖，其中塔架為通長變厚度結構，頂壁厚40 mm，底壁厚90 mm;葉片傾角為5°，切出風速25 m/s，各葉片之間成120°夾角，沿周向均勻分布，葉片長度為60 m，沿翼展各葉素截面的詳細參數如表3所示;機艙尺寸為18 m×6 m×6 m，分別對應長、寬、高方向。根據以上設計參數依次建立塔架、葉片和機艙等部件，通過布爾運算形成大型風力機三維實體模型。

該風力機位于B類地貌，根據葉片與塔架的相對位置，考慮到三葉片體系旋轉過程中存在的周期性，定義葉片與塔架完全重合時為工況1，順時針旋轉15°為一個工況共計8個。此外定義工況1中與塔架完全重合的葉片編號為Ⅰ，順時針旋轉分別為葉片Ⅱ和葉片Ⅲ。不同工況下風力機葉片單周期旋轉方向如圖1所示。

2.2 計算域設置

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

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

計算域左側和頂部邊界條件為速度入口，右側為壓力出口，數值計算采用3D單精度、分離式求解器，流場流速為絕對速度，空氣模型等效為理想不可壓縮流體，計算域入口采用冪指數為0.15的風廓線模型，離地10 m高度處的風速設置18 m/s，采用SIMPLEC算法實現速度與壓力之間的耦合，對流項求解格式為二階，計算過程中設置了網格傾斜校正以提高混合網格計算效果，控制方程的能量計算殘差設置為1×10^-6，最后初始化風場進行迭代計算。圖3給出了平均風速、湍流度剖面模擬結果與理論值的對比曲線，結果表明平均風速和湍流度剖面均與理論值吻合良好。圖4給出了8種工況下塔架30 m高度環向壓力系數分布圖。可知各工況平均風壓系數分布曲線的負壓極值點和分離點對應角度與規范曲線^[11]一致，迎風和側風區域風壓系數數值吻合較好，僅在背風區負壓系數略大于規范值，對比結果驗證了風場數值模擬的有效性。

2.3 沙粒沖擊荷載

采用式（6）分別進行8種工況下風力機塔架、葉片沙荷載計算，圖5?6分別給出了各葉片、塔架沙荷載。由圖可知：風力機不同葉片沙荷載分布呈現不同的規律，葉片Ⅰ沙荷載隨葉片高度增加先減小后增大;葉片Ⅱ沙荷載隨高度增加先增加后減小，且在工況5（60°）達到極大值;葉片Ⅲ沙荷載隨高度減小先減小后增大。

式中n為單位時間內撞擊某區域的沙粒數量，F（τ）為單個沙粒對結構的沖擊力，F_W為該區域風荷載。表5給出了沿塔架子午向高度風、沙荷載比值，可以發現：不同葉片停機位置對風力機塔架表面沙荷載分布影響較大，但均在（0?0.05）H高度范圍內產生該工況下極大沙荷載，工況1，2沙荷載沿高度分布呈現較為明顯的變化趨勢，其隨高度的變化先減小后增大再減小;其余工況沙荷載在塔架0.1H高度以上分布較為一致。

2.4 沙致壓力系數

為定量比較不同停機位置下風力機塔筒表面壓力分布，定義斷面風壓系數、沙致壓力系數、等效壓力系數，如下式所示：

式中C_pi為i點的平均風壓系數;P_i為測點平均壓力（Pa）;P_H為參考高度處遠前方的靜壓（Pa）;ρ為空氣密度，本文數值模擬過程中均取1.225 kg/m³;V_H為參考高度處遠前方的平均風速;C_psi為i點的沙致壓力系數;n為單位時間內撞擊某區域的沙粒數量;F（τ）為單個沙粒對結構的沖擊力;S為計算區域面積;C_pei為i點的等效壓力系數。

其等效思路為：①將表面各監控點沙粒撞擊荷載轉化成沙致壓強;②計算監控點沙致壓強與對應參考高度處風壓比值，即沙壓系數;③將風壓系數與轉化后的沙致壓力系數相加，得到等效壓力系數。等效壓力系數可作為等效目標，定量比較不同停機位下風力機塔架所受風?沙耦合作用大小。

考慮沙粒著點分布特性及葉片與塔架之間的氣動干擾效應，選取塔架未干擾區段沙粒分布密度差異較大，干擾區段塔架與葉片重合面積不同的2個斷面為典型斷面。圖7?9給出各工況下風力機塔架4個典型斷面風壓系數、沙致壓力系數、等效壓力系數對比曲線，分析可知：1）未干擾區（0.025H，0.30H）塔架斷面呈現良好的對稱性;顯著干擾區（0.70H，0.90H）塔架斷面工況1，2，3，6等效壓力系數、風壓系數分布曲線不再保持對稱，塔架迎風面中心點處呈現負壓。2）不同工況下塔架0.025H高度塔架斷面迎風面兩側0°?60°范圍內等效壓力系數與風壓系數差異明顯，沙致壓力系數最高可達0.549;其余典型斷面（0.30H，0.70H，0.90H）等效壓力系數與風壓系數的數值和分布規律基本一致。

3 響應分析

3.1 有限元建模及動力特性

以圖1顯示的8個風力機不同停機位的計算工況為例，本文基于大型通用有限元分析軟件ANSYS建立大型風力機塔架?葉片一體化有限元模型如圖10所示，塔架和葉片采用Shell63單元，機艙采用Beam189單元，環基采用Solid65單元。建立模型后采用Block Lanczos方法求解風力機自振頻率和振型^[22?23];同時，為確保有限元分析的可信度，本文增加了網格無關性驗證，如表6給出了工況1下、不同網格方案下風力機基頻。由表可知，隨著網格總數的增加，風力機基頻呈現逐漸減小的趨勢，而4680網格數和6240網格數的計算結果無明顯差異，綜合計算精度和效率，本文選取網格總數為4680的方案。

圖11給出了不同停機位置下風力機模型前100階模態固有頻率分布曲線及工況2?8相較于工況1風力機有限元模型基頻的增/減量，可以看出大型風力機葉片?塔架耦合模型基頻均較小，僅為0.138 Hz左右，且各模態之間頻率間隔很小。葉片不同停機位置主要影響風力機體系的低階固有頻率，各工況低階頻率出現較小的差別，高階頻率基本一致。可見停機狀態下葉片不同位置對風力機體系的頻率和振型影響微弱，低階振型主要以葉片帶動機艙及塔架進行前后揮舞和左右擺動，高階模態出現塔架及葉片本身的結構變形和失穩形態，如圖12所示。

3.2 塔架響應

圖13?14給出了不同工況下塔架風致徑向位移響應、風?沙致徑向位移響應分布圖。由圖可得：葉片不同停機位置塔架主要影響塔架中上部位移，沙粒沖擊作用則主要影響塔架中下部位移。徑向位移隨著塔架高度的增加而逐漸增大，最大正負位移均出現在塔頂處0°和180°位置;不同停機位置下塔架位移極值在迎風面和背風面的分布范圍不同，極值區域隨葉片與塔架位置相對干涉作用增大而增大;工況1葉片完全遮擋塔架，塔架頂部位移出現“逆向效應”，沙粒沖擊作用放大，塔架頂部風?沙致位移為負。

為對比葉片位置停機位置、沙粒沖擊作用引起塔架底部彎矩值的差異，分別以工況1在凈風、風?沙共同作用下響應值作為初始狀態，將各工況下塔底內力響應值與之作差，得到塔架風致內力響應特征值、風?沙致內力響應特征值;再將風致、風?沙致內力響應值作差，得到沙致響應特征值。

如圖15所示：塔底徑向彎矩只在0°，10°，70°以及320°位置出現差異，工況2在0°位置數值相差最大達到146.82 N·m;風致塔底徑向彎矩與風?沙致環向彎矩在迎風面差異顯著，沙致塔底徑向彎矩最高達12.78 N·m;背風面環向彎矩呈左右對稱，風致響應與風?沙致響應分布基本一致，沙粒主要沖擊塔架迎風面左右60°，沙致彎矩響應極值均出現在0°位置，背風面及側風面可忽略沙粒沖擊作用。

3.3 屈曲穩定性分析

圖16給出了不同工況下風力機塔架?葉片耦合體系屈曲位移和臨界風速對比示意圖。對比可知：1）葉片停機位置對整體結構屈曲失穩時最大位移影響顯著，隨葉片對塔架遮擋面積的減小最大位移逐漸增大，在工況5時達到峰值;2）風與風?沙耦合作用下臨界風速、屈曲最大位移分布曲線規律基本一致，風?沙耦合作用使得臨界風速略微增大;工況3對應的臨界風速最低，工況5對應的臨界風速最高。

4 結? 論

本文系統探討了風?沙耦合環境下大型風力機體系氣動力和結構響應特性，主要涉及風?沙雙向耦合數值模擬、風荷載分布與風致響應特性、沙荷載分布與沙致響應特性、風?沙致等效荷載及其效應等。得到如下主要結論：

1）不同葉片停機位置對風力機塔架表面沙荷載分布影響較大，但均在（0?0.05）H高度范圍內產生該工況下極大沙荷載，工況6（75°）下該高度范圍內沙荷載與風荷載比值達23.70%;且隨高度增加，沙荷載與風荷載比值逐漸減小。

2）不同工況下沙粒撞擊位置主要集中在塔架0.6H（未干擾段）高度以下，集中于迎風區域兩側各60°范圍，沙致壓力系數最大值為0.549，發生在工況4的塔架底部迎風面。考慮葉片停機位置影響時，葉片遮擋效應對塔架中上部位移影響顯著，沙粒沖擊作用則加大塔架中下部位移。

3）塔架徑向位移分布規律差異較大，工況1（葉片完全遮擋）徑向位移對沙粒沖擊作用敏感，徑向位移隨著塔架高度的增加而逐漸增大，最大正負位移均出現在塔頂處0°和180°位置。沙致彎矩響應極值均出現在0°位置，總體分布大致呈左右對稱，背風面及側風面可忽略沙粒沖擊作用。

參考文獻：

[1]??????? Bagnold R A. The movement of desert sand[J]. Geographical Journal， 1935，85（4）：342-365.

[2]??????? Bagnold R A. The transport of sand by wind[J]. Geographical Journal， 1937，89（5）：409-438.

[3]??????? Cheng J， Jiang F， Xue C， et al. Computational method for maximum sediment discharge and sand-carrying wind load in the prevention and treatment of wind drift sand for railway in strong wind area[J]. China Railway Science， 2012，33（1）：1-5.

[4]??????? 蔣富強，李? 熒，李凱崇，等.蘭新鐵路百里風區風沙流結構特性研究[J].鐵道學報，2010，32（3）：105-110.

Jiang Fuqiang， Li Ying， Li Chongkai， et al. Study on structural characteristics of Gobi wind sand flow in 100 km wind area along Lan-Xin Railway[J]. Journal of the China Railway Society， 2010，32（3）：105-110.

[5]??????? 樓文娟，余? 江，潘小濤.風力機葉片揮舞擺振氣彈失穩分析[J].工程力學，2015，32（11）：236-242.

Lou Wenjun， Yu Jiang， Pan Xiaotao. Calculating for aerodynamic stability response of wind turbine blade in flapwise and edgewise direction?[J]. Engineering Mechanics， 2015，32（11）：236-242.

[6]??????? 柯世堂，余? 瑋，王同光.停機狀態葉片位置對風力機體系氣動性能影響[J].浙江大學學報（工學版），2016，50（7）：1230-1238.

Ke Shitang， Yu Wei， Wang Tongguang. Impact for blade position on aerodynamic performance of wind turbine system under stopped status[J]. Journal of Zhejiang University （Engineering Science）， 2016，50（7）：1230-1238.

[7]??????? 李德源，莫文威，夏鴻建，等.水平軸風力機柔性葉片氣彈耦合分析[J].太陽能學報，2015，36（3）：734-742.

Li Deyuan， Mo Wenwei， Xia Hongjian， et al. Aeroelastic coupling analysis of flexible blades for horizontal axis wind turbines[J]. Chinese Journal of Solar Energy， 2015，36（3）：734-742.

[8]??????? 柯世堂，余? 瑋，王同光.基于大渦模擬考慮葉片停機位置大型風力機風振響應分析[J].振動與沖擊，2017，36（7）：92-98.

Ke Shitang， Yu Wei， Wang Tongguang. Analysis of wind-induced vibration response of large wind turbine based on large eddy simulation[J]. Journal of Vibration and Shock， 2017，36（7）：92-98.

[9]??????? 張? 永，黃? 超，劉? 召，等.挾沙風作用下風力機葉片涂層沖蝕特性研究[J].材料導報， 2016，30（10）：95-99.

Zhang Yong， Huang Chao， Liu Zhao， et al. Study on erosion characteristics of wind turbine blade coating under the effect of wind and sand[J]. Material Review， 2016，30（10）：95-99.

[10]????? 李仁年，范向增，李德順，等.風沙環境下水平軸風力機的轉矩特性[J].蘭州理工大學學報，2015，41（4）：55-59.

Li Rennian， Fan Xiangzeng， Li Deshun， et al. Torque characteristics of horizontal-axis wind turbine in sand-blowing environment[J]. Journal of Lanzhou University of Technology， 2015， 41（4）： 55-59.

[11]????? 中華人民共和國住房和城鄉建設部.GB50009-2012，建筑結構荷載規范[S]. 北京： 中國建筑工業出版社， 2012.

MOHURD. GB50009-2012， Building structure load specification [S]. Beijing： China Building Industry Press， 2012.

[12]????? 中華人民共和國國家發展和改革委員會.DL/T 5383-2007， 風力機發電場設計技術規范[S].北京：中國電力出版社， 2007.

National Development and Reform Commision. DL/T 5383-2007， Technical specification of wind power plant design[S]. Beijing： China Electric Power Press， 2007.

[13]????? 中華人民共和國國家質量監督檢驗檢疫總局，中國國家標準化管理委員會.GB/T 20319-2006， 風力發電機組驗收規范[S].北京：中國機械工業出版社， 2006.

General Administration of Quality Supervision， Inspection and Quarantine of the People's Republic of China， Standardization Administration of the People's Republic of China. GB/T 20319-2006. Code for acceptance of wind turbines generator systems[S]. Beijing： Chinese Machinery Industry Press， 2006.

[14]????? IEC 61400-21：2001. Wind turbine generator systems-Part21： Measurement and assessment of power quality characteristics of grid connected wind turbines[S]. International Elector-technical Commission， 2001.

[15]????? IEC 61400-23：2001. Wind turbine generator systems-Part23：Full-scale structural testing of rotor blades[S]. International Elector-technical Commission， 2001.

[16]????? 唐國利，巢清塵.近48年中國沙塵暴的時空分布特征及其變化[J].應用氣象學報， 2005，16（b03）：128-132.

Tang Guoli， Chao Qingchen. The temporal and spatial distribution characteristics and variation of sandstorm in China in recent 48 years[J]. Journal of Applied Meteorology， 2005，16（b03）：128-132.

[17]????? 戴雪榮，師育新，薛? 濱.蘭州現代特大塵暴沉積物粒度特征及其意義[J].蘭州大學學報（自科版）， 1995，31（4）：168-174.

Dai Xuerong， Shi Yuxin， Xue Bin. Particle size characteristics and its significance of modern duststorms in lanzhou[J]. Journal of Lanzhou University（Natural Science）， 1995，31（4）：168-174.

[18]????? 董依帆，柯世堂，楊? 青.考慮風-沙雙向耦合作用大型風力機體系氣動力分布研究[J].振動與沖擊， 2019，38（16）：84-92.

Dong Yifan， Ke Shitang， Yang Qing. Wind-sand movement characteristics and aerodynamic distribution of large wind turbine systems considering bidirectional coupling effect?[J]. Journal of Vibration and Shock， 2019，38（16）：84-92.

[19]????? 李天涯.沖量及動量定理的示例探討[J].物理教學探討， 2014，32（12）：36-36.

Li Tianya. An example of the impulse and momentum theorem[J]. Discussion on Physics Teaching， 2014，32（12）：36-36.

[20]????? 楊慶峰.用牛頓第二定律解決問題的幾種方法[J].數理化解題研究，2018，（13）： 78-79.

[21]????? 羅貝格.彈性及彈-塑性介質中的沖擊波[M].北京：科學出版社， 1965.

Robberg. Shock Waves in Elastic and Elastic-Plastic Media [M]. Beijing： Science Press， 1965.

[22]????? Ke Shitang， Wang Tongguang， Ge Yaojun， et al. Wind-induced responses and equivalent static wind loads of tower-blade coupled large wind turbine system[J]. Structural Engineering and Mechanics， 2014，52（3）：485-505.

[23]????? 鄧? 露，肖志穎，黃民希，等.考慮流固耦合的近海風機動力響應數值計算[J].湖南大學學報（自然科學版）， 2015，42（7）：1-8.

DENG Lu， XIAO Zhiying， HUANG Mingxi， et al. Numerical simulation of dynamic response for offshore wind turbines including fluid-structure interaction[J]. Journal of Hunan University （Natural Sciences）， 2015，42（7）：1-8.

Abstract： The position of the blades in a strong wind-off state will significantly affect the flow around wind turbine towers and stability of them. Especially in the extreme weather conditions of sandstorms， the adhesion of sand will also affect the turbulence characteristics of the wind and generate additional impact force on the tower. The existing researches lack the exploration of wind-sand coupling effect to explore the aerodynamic performance and structural response of large-scale wind turbines. Taking the 5 MW horizontal-axis wind turbine system independently developed by Nanjing University of Aeronautics and Astronautics as the object and wind-sand-wave two-way coupling algorithm as the core， the simultaneous iterative simulation of the typical wind field and sand combination is performed on the continuous phase and discrete phase model respectively by CFD technology. Considering the eight stop positions in a single rotation cycle of the blade， the equivalent pressure coefficient and the aerodynamic load distribution characteristics are compared and analyzed. Then combined with the finite element method， the dynamic characteristics， wind-sand induced structural response and buckling stability of the wind turbine system at different stopping positions are discussed. The results show that the load effect of sand particles on the windward surface of wind turbine tower is the most significant under the wind-sand coupling environment. The pressure coefficient of sand in some regions can reach 0.55， and the ratio of sand impact load to wind load can reach 23.7%. The wind load under different stand positions plays a dominant role in the response of the tower structure， but the effect of sand impact on the bottom of the tower cannot be ignored.

Key words： large wind turbine system; wind-sand coupling environment; structural response; CFD numerical simulation; stop position

作者簡介： 柯世堂（1982-），男，博士，教授。電話：13621581707; E-mail： keshitang@163.com