引言

復合材料葉片是風機設備中將風能轉化為機械能的關鍵部件[1]。目前，葉片尺寸正在朝著大型化的方向發展，而其結構性能試驗的成本也隨之增加，因此，找到一種有效的結構計算分析方法對于節約成本以及結構校核和開發新型葉片就顯得尤為重要。

隨著計算機技術的發展，有限元法在結構分析中得到了廣泛的應用。有限元強大的建模和結構分析功能適用于復合材料葉片的應力、變形、頻率、屈曲、疲勞及葉根強度分析。ANSYS是一款著名的商業化大型通用有限元軟件，廣泛應用于航空航天、機械制造等領域。ANSYS多物理場仿真及耦合的獨特功能，以及200多種單元類型，可以對復合材料結構設計、材料研究及制造工藝提供完整的解決方案。總之，對于復合材料結構計算分析，完全可以通過ANSYS程序來實現[2]。

但是，因復合材料葉片結構的特殊性，例如：①形狀不規則（每個截面都不同）；②鋪層復雜，過渡層很多；③大量夾層結構（“三明治”結構）；④大量粘結區域。所以，其有限元模型的建立是葉片結構有限元分析中的一大難題，而單元類型的選擇又決定著建立有限元模型的難易。目前，復合材料風電葉片有限元模型在單元類型的選擇上主要采用三種單元類型：shell99殼單元、shell91殼單元、Solid46實體單元。選擇實體單元，雖然能提高有限元的計算精度，但是建立葉片的有限元模型會花費大量的工作時間，且很難定義單元坐標，這非常不利于工程上葉片的結構校核及分析；選擇殼單元，可以方便地設置和修改鋪層厚度，單元坐標的設置容易實現，建模和計算時間比采用實體單元少，這極大地提高了工作效率，而且其計算精度完全可以滿足工程需要。

因此，本文使用shell99殼單元，通過三維建模，建立了葉片的有限元模型，并以懸臂梁的方式，對葉片的模態和靜力變形進行了計算分析，通過計算，得出了葉片的重量、振型及最大變形，并與試驗數據進行了對比。 

1 有限元模型的建立

通常，在整個有限元求解過程中最重要的環節是有限元前處理模型的建立。一般包括幾何建模、定義材料屬性和實常數（要根據單元的幾何特性來設置，有些單元沒有實常數）、定義單元類型，網格劃分、添加約束與載荷等。

由于葉片形狀復雜，而一般有限元軟件所提供的幾何建模工具功能相當有限，所以在ANSYS中難以快速方便地對其建模。因此，針對較復雜的結構，可以先在三維CAD軟件（如在PROE中）建立幾何模型，然后在有限元分析軟件ANSYS中通過輸入接口讀入實體模型，最后，在ANSYS環境下，通過幾何修補和簡化、板殼中面抽取、節點偏置、網格自動劃分等技術對葉片模型進行處理，并形成高效準確的有限元模型，使之適用于CAE分析。

1.1 單元設置與材料屬性

針對葉片中的梁、殼等復合材料層合結構，ANSYS提供了一系列的特殊單元——結構多層復合材料單元，以模擬各種復合材料。鋪層單元中可以考慮復合材料特有的鋪層特性和各向異性特性。

本計算采用的是相對簡單的線性鋪層單元Shell99。該單元是一種八節點3D殼單元，每個節點有六個自由度，主要適用于薄到中等厚度的板和殼結構，一般要求寬厚比應大于10。Shell99可實現多達250層的等厚材料層，或者125層厚度在單元面內呈現雙線性變化的不等厚材料層。如果材料層大于250，用戶可通過輸入自己的材料矩陣形式來建立模型，還可以通過一個選項將單元節點偏置到結構的表層或底層。

單元鋪層主要是確定纖維方向和纖維量，是復合材料風電葉片結構設計的一個重要環節。鋪層設計的優劣在很大程度上決定著結構設計的成敗[3]。

本計算的鋪層完全按照工藝鋪層進行設計。在ANSYS環境下，針對Shell99單元，通常有兩種方法來定義材料層的配置：①通過定義各層材料的性質；②通過定義表示宏觀力、力矩與宏觀應變、曲率之間相互關系的本構矩陣。第一種方法是由下到上一層一層定義材料層的配置，底層為第一層，后續的層沿單元坐標系的Z軸正方向自底向上疊加，對于每一層材料，由單元實常數表來定義材料性質、鋪層方向角、厚度，如圖1所示為葉片某部分的單元鋪層；第二種方法是定義各層材料性質的另一種方式，矩陣表示了單元的力-力矩與應變-曲率的關系，必須在ANSYS外進行計算。

 

圖1.單元鋪層圖

Fig.1 element laminated diagram

葉片的材料體系為玻纖/環氧，葉片制作采用真空灌注工藝，所用復合材料有：三軸向玻璃布、雙軸向玻璃布、單軸向玻璃布、PVC泡沫、Balsa木、氈等。玻璃鋼復合材料與泡沫材料的主要力學性能見表1、表2所示，其中，Ex為材料的纖維方向，玻璃鋼復合材料的密度取ρ=1888kg/m3，Balsa木密度取ρ=150kg/m3，PVC密度取ρ=80kg/m3。

表1.玻璃鋼復合材料力學性能

Table 1. Mechanical properties of FRP materials

名稱	符號	單位	UD	Biaxial	Triaxial
玻纖/環氧	Ex	Mpa	39000	11400	28500
	Ey	Mpa	8920	11400	13500
	Ez	Mpa	8920	8920	8920

表2.泡沫材料力學性能

Table 2. Mechanical properties of foam materials

名稱	符號	單位	模量值
Balsa木	Ex	Mpa	1000
	Ey	Mpa	35
PVC	Ex	Mpa	65

1.2 模型建立與網格劃分

首先，依據三維坐標變換原理求解出葉片空間截面翼型的實際位置，然后以大型三維軟件PROE為工作平臺，通過導入空間坐標點，生成B樣條曲線，如圖2所示為本計算模型的三維線框圖。

圖2.葉片線框投影圖

Fig.2 Blade frame drawing

　　最后，由曲面掃掠命令生成葉片三維外形圖，再結合曲面曲線分析命令對所生成的曲線、曲面進行檢驗和修改，直至生成符合要求的葉片三維外形圖，如圖3所示。將生成的三維模型轉化為IGES格式文件，為后續建立有限元模型做準備。

 

圖3.葉片外形圖

Fig.3 Blade outline diagram

將PROE導出的IGES格式文件，輸入到ANSYS系統中，得到了ANSYS環境下的葉片三維模型。采用Shell99單元對葉片殼體、梁、腹板進行網格劃分，有限元模型單元數為29914，節點數為88680，如圖4所示。

 

圖4.網格劃分圖

Fig.4 Meshing diagram

1.3 約束與載荷

葉片根部采用剛性固定的約束形式，即根部所在節點的6個自由度被固定，整個葉片簡化為懸臂梁模型。加載方式與試驗加載方法保持一致，第一個工況，即在flapwise方向，選擇單點加載，施加集中力39KN，如圖5所示；第二個工況，即在edgewise方向，選擇四點加載，從左至右分別施加集中力57.5KN、21.4KN、20.4KN、36.2KN，如圖6所示。

 

2 計算結果與分析

2.1 質量計算結果

表3即為ANSYS輸出的風電葉片質量計算結果，重心位置與實測值基本吻合，葉片總質量低于實際值。產生葉片質量計算值比實測值小的主要原因是葉片灌膠后，其泡沫的密度應該大于實際泡沫的密度，以及建立的葉片有限元模型沒有考慮實際葉片中的附件（如接閃器等金屬件）重量等。

表3.風電葉片質量計算結果

Table 3. The mass result of the blade

	質量[kg]	重心（距葉根）[m]
實際值	5950	12.1
有限元計算值	5571	11.9
誤差	6.4%	1.65%

2.2 模態分析結果

表4即為ANSYS輸出的風電葉片一階固有頻率計算結果，并且提取了葉片的前五階振型，如圖7所示。

表4.風電葉片一階固有頻率計算結果

Table 4. The first natural frequency of the blade

	頻率值[Hz]
	一階flapwise方向	一階edgewise方向
試驗值	0.8	1.46
計算值	0.89	1.61
誤差	11.25%	10.27%

由圖7可知，一階頻率為flapwise方向一階固有頻率，二階頻率為edgewise方向一階固有頻率。由表4知葉片一階固有頻率的計算值比實測值大，造成計算值偏大的主要原因是葉根約束方式與試驗（通過螺栓固定）不一致，以及計算質量小于實際質量等。

 

圖7.振型圖

Fig.7 Vibration mode diagram

2.3 靜力分析結果

　　表5即為兩種工況下計算出的風電葉片最大撓度值，計算結果與實測值吻合較好。

表５.兩種載荷工況下風電葉片的最大撓度

Table 5. The most displacement of the blade in the two load case

	葉片最大撓度[m]
	flapwise方向	edgewise方向
試驗值	5.36	1.11
計算值	4.83	1.01
誤差	9.9%	9%

葉片變形如圖8、圖9所示。

 

圖9.edgewise方向變形圖

Fig.9 Deformation diagram in edgewise

3 結論

(1) 采用殼單元模擬風電葉片計算出葉片總質量、撓度變形能與實測結果相對誤差小于10%，證明了該方法在工程應用上的可行性和可靠性。

(2) 由于葉根約束方式與試驗（通過螺栓固定）不一致以及計算質量小于實際質量等原因，葉片固有頻率的計算值略高于實測結果。

(3) 采用殼單元計算風電葉片剛度，既可保證計算結果的可靠性又可縮短建模時間提高工作效率，對風電葉片結構分析的實際工程應用具有重要價值。