風力機復合材料葉片流固耦合分析*

張 旭， 劉安宇

(天津工業大學 機械工程學院， 天津 300387)

風力發電已經成為世界上最具發展前景的能源開發方式.目前，利用風能最為行之有效的裝置是風力發電機，其中最具重要性的部件即為葉片，因為其性能與風力機運行現狀及穩定性存在直接關聯.在風力機的運行過程中，葉片周圍空氣的流動對葉片產生的作用力促使葉片出現一定程度的形變，這一變化會導致周邊空氣的流向發生變化.為了得到精準的葉片動力學響應，首要要考慮葉片及流場間存在的作用，即必須考慮葉片與其周圍流場的耦合作用.

國外研究人員對流體和固體風力發電機組進行了大量研究.Tongchitpakdee等對偏航條件下的大型風機進行了數值模擬計算，并對其氣動性能進行了系統分析[1]；Rachid等在定常情況及非定常情況下對NACA4415翼型的葉片彎曲及扭轉效應進行分析研究，延長了風力機的使用壽命[2]；Kim等考慮了氣動彈性對葉片的影響，并對大型風力機葉片性能進行了預測[3].

自80年代中期至今，國內研究人員對實體研究流固耦合問題與風力發電機整機的流固耦合進行了大量的研究.徐建源等利用ADAMS軟件建立了風波聯合作用海上風力機模型，并對其進行流固耦合分析，得到了復雜工況下的風力機葉片的動態響應參數[4]；李少華等利用Fluent軟件并結合SSTk-ω模型對1.2 MW風力機旋轉風輪流場進行模擬，但是并沒有考慮風力機塔筒與流場間的相互作用對數據造成的影響[5]；胡丹梅等在水平軸風力機模型不同尖速比條件下對風輪下游流場速度進行測量，獲得風力機葉片尾跡的流場定量信息，由于該測量方法容易受外界因素的影響，因此，流動情況與實際情況有一定誤差[6-7]；田琳琳等結合制動盤理論與CFD方法，并通過Fluent軟件對9臺風力機的尾流互相干擾情況進行模擬，得到風力機尾流相互影響的結果[8]；陳海萍等研究了葉片周圍流場對葉片結構特性的影響[9]；潘萍萍等通過Fluent軟件對1.5 MW風力發電機塔筒進行流固耦合分析，得到了風力機塔筒的動態評估方法[10].

本文應用ANSYS軟件對改進后的2 MW風力機葉片模型進行流固耦合分析，通過對傳統翼型NACA4412的改進，應用鋪層設計把葉梢、葉根、展向等厚度改成漸變式.首先得到風力機葉片及風輪的流場計算結果，再將風力機流場的計算結果施加到固體計算中，得到風力機的應力分布，進而得出應力最大值及變形量最大位置，所得結果可以為風力機的設計和安全運行提供指導[11].

1 流固耦合方法

流固耦合問題從數據傳輸的角度來看，其計算可以分為兩種：單向流固耦合和雙向流固耦合.單向流固耦合是指通過流體計算獲得的壓力、速度等數據可以通過連接界面導入固體結構，或者將從固體計算獲得的網格位移轉移到流體的計算.

計算流體力學簡稱CFD，主要使用離散點上的一系列變量值代替時域和空間域中的一系列物理量.通過建立可以表示變量之間關系的代數方程，將基本方程離散化，通過CFD的數值模擬，可以顯示流動領域的運動.在通過CFD分析的工程應用中，可以深入研究耦合問題以實現實驗指南，達到節約成本的目的.流固耦合分析流程圖如圖1所示.

圖1 耦合系統流程示意圖Fig.1 Schematic diagram of coupling system

2 三維模型建立

本文利用NASYS APDL建立風力機葉片的三維實體模型.采用傳統翼型NACA4412應用復合材料建立葉片有限元模型如圖2a所示，該模型總長44.5 m，翼型最大弦長3.1 m，根部為直徑2.34 m的圓形，質量為10.6 t.為了加強葉片穩定性，在葉片中部添加兩個厚為1.6 cm的加強筋，葉根部也進行了加厚，具體模型如圖2b所示.葉片的建模工作主要在傳統翼型NACA4412模型下進行改進，應用鋪層設計把葉梢、葉根、展向等厚度改成漸變式，這樣不僅更加接近實體模型，在強度、剛度上也有所加強，風輪模型如圖2c所示.

3 靜止單葉片流固耦合分析

3.1 流場計算分析

根據葉片結構在ANSYS軟件中繪得實體模型，有針對地對流體實際范圍進行了設定，同時進行了網格劃分，不過因常規風力機體積過大，因此在制定計算域的過程中需要將其設置為葉片的數倍以上.設置的流體域范圍為202 m×103 m×63 m，所劃分的網格數為4 601 159，計算域及計算域網格模型如圖3所示.

在進行流場分析時應先抑制葉片固體域，以實現對流體域的保護，設置流體控制溫度為25 ℃.依照入口實際條件對其速度進行控制，并將這一設定速度稱為恒速，流體的方向則需要垂直于入口界面，結合實際情況設置流體速度大小為13 m/s，而流體出口處設定為自由流出，風力機葉片邊界條件設置為光滑無滑移壁面.風場流速圖如圖4a所示，葉片在迎風面和背風面的壓力分布如圖4b、c所示，在背風面葉片受到較小的壓力，而迎風面葉片受到壓力最大.

圖2 葉片及風輪模型Fig.2 Models for blade and wind wheel

圖3 計算域及計算域網格模型(靜止狀態)Fig.3 Models for computational domain and computationaldomain grid (static state)

圖4 風場流速及迎風面、背風面壓力分布Fig.4 Wind field velocity and pressure distributionof windward and leeward surfaces

3.2 葉片靜力分析

葉片翼面材料為玻璃鋼，該材料為線性材料，用MP命令定義其材料特性.其各層的線性材料特性為正交異性，具體參數為：展向模量42.6×109Pa，徑向模量16.5×109Pa，剪切模量5.5×109Pa，泊松比0.22，密度1 950 kg/m3.

通過Static Structural接收風場的分析結果數據，并實現結構靜力的相關研究.首先單獨保留葉片，確保流程部分處于制約狀態.由于葉片實際網格劃分直接影響最后的研究結果，因此在劃分過程中要確保其規則性.在Static Structural模塊下選擇Engineering Data，定義材料參數密度、彈性模量、泊松比和剪切模量，導入固體模型，定義網格劃分參數劃分網格，設置的網格數為980 916，定義固定約束.劃分網格后的模型如圖5所示.

圖5 葉片網格模型Fig.5 Model for blade grid

4 旋轉狀態風輪流固耦合分析

4.1 流場計算分析

將Fluent模塊的Solution與Static Structural模塊的Setup建立連接進行單向流固耦合分析，導入風輪分別建立旋轉流體域和外部空氣流體域，旋轉流體域半徑為48 m，高度為6 m.外部空氣流體域尺寸為606 m×296 m×296 m，設置網格劃分參數及膨脹層，計算域及計算域網格模型如圖6所示.

圖6 計算域及計算域網格模型(旋轉狀態)Fig.6 Models for computational domain and computationaldomain grid (rotation state)

設置室溫為25 ℃，入口速度為13 m/s，定義旋轉軸為y軸，旋轉速度為1.894 49 rad/s，風場流速圖如圖7a所示，通過Contour云圖顯示的壓力云圖如圖7b所示，可以看出壓力從葉根到葉尖逐漸減少，葉根最大，而在葉尖處最小.

圖7 風場流速及風輪壓力分布Fig.7 Wind field velocity and pressuredistribution of wind wheel

4.2 風輪結構靜力分析

在Static Structural模塊下選擇Engineering Data進行材料參數定義：密度為1 950 kg/m3，彈性模量為42.6×109Pa，泊松比為0.22，剪切模量為5.5×109Pa，導入固體模型，劃分網格，設置葉輪網格數為1 393 827，定義固定約束.劃分網格后的模型如圖8所示.

圖8 風輪網格模型Fig.8 Model for wind wheel grid

4.3 結果輸出

導入流體分析結果即壓力，選擇導入壓力的面為葉片外表面，得到應變分布如圖9所示.由圖9可知，在額定速度下風力機運行時，風力機結構受風場的影響出現變形，變形的趨勢與在風場中所受到的葉片表面壓力分布是基本吻合的，此外在形變量最大的葉尖處，相對于葉根處的約束部分最遠，并且此處葉片材料是最薄的，與文獻[11]實驗結果基本吻合.

圖9 風輪應變分布Fig.9 Stain distribution of wind wheel

5 結 論

本文合理應用ANSYS Workbench實現風力系統中風輪及葉片的單向流固耦合分析，得到兩者在運行數據上存在的關聯.實驗中得到的數據為進行葉片疲勞及共振頻率測試提供了一定的參考，為之后的研究提供基礎.在實際運行工作中，葉片的固有頻率為耦合之后整個系統的固有頻率，這個固有頻率包括了風力機塔架在系統中對葉片的影響，系統中旋轉離心力的作用，風場在系統中對葉片的影響，系統中重力作用等一系列的作用，這些影響與單個葉片在靜止狀態下的模態頻率有一定的區別，但一般而言頻率的變化很小.在今后的葉片優化設計過程中，同樣需要關注該部分問題，以便得到的模型與實際運作更為接近，減小并優化設計周期.