汪 泉,王馮云,胡 聰,柳柏楊

(湖北工業(yè)大學(xué)機械工程學(xué)院,湖北 武漢430068)

1 引言

風(fēng)能是一種綠色可再生能源,越來越受到世界各國的重視。風(fēng)力發(fā)電機是一種將風(fēng)能轉(zhuǎn)化成電能的重要機器,而葉片是風(fēng)力機中最關(guān)鍵的零部件之一,占整機價格的20%左右。目前,風(fēng)力發(fā)電機葉片的設(shè)計主要包括葉片氣動外形設(shè)計和內(nèi)部結(jié)構(gòu)鋪層設(shè)計[1]。Wang等[2]提出了一種將葉片氣動外形與結(jié)構(gòu)拓?fù)鋬?yōu)化相耦合的方法,利用葉片外形的非均勻有理B樣條曲線(NURBS)控制點位置作為形狀設(shè)計變量參數(shù),葉片橫截面的有限元的相對密度作為拓?fù)湓O(shè)計變量的特征參數(shù),使水平軸風(fēng)力機具有更高功率輸出和更優(yōu)的結(jié)構(gòu)性能。Akram,MT等[3]提出了NREL S809 翼型的空氣動力學(xué)優(yōu)化設(shè)計方法,以提高風(fēng)力機葉片設(shè)計的性能,同時開發(fā)了基于遺傳算法的集成代碼,利用類形狀變換(CST)和參數(shù)截面(PARSEC)對NREL S809翼型進行優(yōu)化,分析其氣動特性。T. Ashuri等[4]集成風(fēng)輪及塔架耦合特性,以最小成本為目標(biāo)函數(shù),對某5MW海上風(fēng)力機葉片進行多學(xué)科優(yōu)化設(shè)計,優(yōu)化結(jié)果顯示其發(fā)電量成本減少了2.3%。王瓏,王同光[5]等提出基于均勻分解與差分進化算法的風(fēng)力機葉片多目標(biāo)優(yōu)化設(shè)計方法,針對某1.5MW風(fēng)力機葉片,分別以兩個目標(biāo)函數(shù)、三個目標(biāo)函數(shù)及四個目標(biāo)函數(shù)進行優(yōu)化設(shè)計,從而驗證了該改進算法具有均勻分布性及有效收斂性等優(yōu)點。楊陽,李春[6]等以年發(fā)電量最大和葉片質(zhì)量最輕為優(yōu)化目標(biāo)函數(shù),通過多目標(biāo)遺傳算法設(shè)計某5MW風(fēng)力機葉片,得到Pareto分布優(yōu)化解集。并與NREL 5MW風(fēng)力機葉片比較,結(jié)果表明年發(fā)電量提高了3.3%,質(zhì)量減少了8.7%。此外,還有很多學(xué)者[7-10]在風(fēng)電葉片的設(shè)計方面做出了卓越的貢獻,均具有一定的價值。

近年來,逐漸有學(xué)者將等幾何分析方法[11]應(yīng)用到葉片結(jié)構(gòu)有限元分析中。Y. Bazilevs,M. C. Hsu等[12]將基于NURBS(非均勻有理B樣條)曲面的等幾何分析思想引入到風(fēng)力機葉片曲面建模及流體網(wǎng)格劃分中,構(gòu)建葉片無主梁流固耦合模型,并通過文獻資料驗證了該流固耦合方法的正確性。EtanaFerede[13]將風(fēng)力機葉片鋪層材料特性、截面剛度特性等參數(shù)簡化到彈性梁模型中去,利用等幾何分析方法通過NURBS幾何控制點來研究葉片結(jié)構(gòu)變形。陳龍等[14]在等幾何中加入靈敏度分析,通過懸臂梁形狀優(yōu)化驗證了方法可行性以及收斂速度快的特點。Austin J. Herrema等[15]提出基于Kirchoff-Love殼理論的葉片位移與扭轉(zhuǎn)連續(xù)性罰函數(shù)表達式,以解決葉片等幾何分析過程中NURBS曲面分塊非匹配性問題,通過葉片變形、應(yīng)變及屈曲分析來確定罰函數(shù)因子的范圍。上述文獻對風(fēng)力機葉片等幾何分析作了初步的研究,但并未對復(fù)雜葉片曲面建模過程進行深入研究,大多對葉片進行了一定的簡化,并沒有將葉片氣動外形與等幾何分析模型有機統(tǒng)一起來。

因此,該文以風(fēng)力機葉片復(fù)雜曲面為研究對象,在風(fēng)力機葉片三維曲面集成表達理論的基礎(chǔ)上構(gòu)建變弦長、變扭角的風(fēng)力機葉片通用設(shè)計平臺,并以9個翼型為原始數(shù)據(jù),建立風(fēng)力機葉片復(fù)雜曲面模型,然后結(jié)合NURBS和T-spline曲面技術(shù),將曲面導(dǎo)入Rhino軟件平臺進行等幾何分析,并解決葉片邊界條件離散化繁瑣、計算精度低等問題。同時,將等幾何分析計算結(jié)果與通用有限元軟件ABAQUS計算出來的結(jié)果進行對比分析,從而驗證了本文提出的方法的可行性。

2 風(fēng)力機葉片復(fù)雜曲面模型構(gòu)建

2.1 葉片統(tǒng)一三維集成表達公式

變弦長,變扭角風(fēng)力葉片的統(tǒng)一的三維集成表達[1]形式如式(1)所示

(1)

式中:ρ(θ)—翼型形狀函數(shù);

xM—氣動中心展向位置

2.2 風(fēng)力機葉片復(fù)雜曲面模型構(gòu)建

本文選用某2WM風(fēng)力機葉片進行復(fù)雜曲面建模,其翼型類型如圖1所示,位置分布如表1所示。

表1 翼型位置分布

圖1 某2MW風(fēng)力機葉片翼型類型

該葉片的弦長及扭角分布如圖2和3所示。

圖2 葉片弦長分布

圖3 葉片扭角分布

根據(jù)葉片弦長、扭角分布數(shù)據(jù),可擬合出葉片弦長及扭角高階多項式系數(shù):

C(u)=5.052u4+10.869u3+5.105u2+3.803u+3.985

(2)

β(u)=17.353u4-49.655u3+57.352u2-42.234u+17.181

(3)

2.3 風(fēng)力機葉片復(fù)雜曲面模型構(gòu)建

基于風(fēng)力機葉片的設(shè)計參數(shù)并結(jié)合復(fù)雜三維葉片曲面集成表達式(1),利用Matlab編制程序,輸出葉片截面數(shù)據(jù)點,構(gòu)造風(fēng)力機葉片曲面模型,如圖4所示。

圖4 葉片三維曲面圖

3 風(fēng)力機葉片NURBS參數(shù)化

為了對風(fēng)力機葉片進行等幾何分析,需將上述葉片曲面數(shù)據(jù)需進行NURBS參數(shù)化。

3.1 NURBS曲線

P次B樣條曲線的定義可以用如下公式來進行定義

(4)

式中:U—節(jié)點失量

U={a,…,a,up+1,…,um-p-1,b,…b}

{pi}—控制點

{Ni,p(u)}—節(jié)點矢量

B樣條曲面是二維方向控制的點的網(wǎng)格,因此需要用兩個互相垂直的基函數(shù)進行定義

(5)

U={0,…,0,up+1,…ur-p-1,1,…,1}和V={0,…,0,vp+1,…vs-p-1,1,…,1}為節(jié)點矢量。

翼型NACA63415形狀如圖5所示,當(dāng)節(jié)點矢量U={0,0,0,0,1/5,1/2,3/5,4/5,1,1,1,1},節(jié)點個數(shù)m+1=12,用3次B樣條繪制NACA63415翼型,即p=3,根據(jù)公式m=n+p+1,可求得控制點個數(shù)n+1=m-p=8,即表2。

表2 翼型控制點

圖5 NACA63415翼型廓線

3次基函數(shù)利用三角形陣列進行求解,可求得翼型NACA63415的基函數(shù),以及由B樣條曲線生成的翼型,如圖6、圖7所示。

圖6 翼型基函數(shù)

圖7 B樣條曲線翼型

按照定義,p次NURBS曲線可以用如下公式進行表示

(6)

式中:{pi}—控制點;

{wi}—權(quán)因子;

{Ni,p(u)}—定義在節(jié)點矢量U上的p次B樣條基函數(shù);

同時,為了簡化式(6),令:

(7)

于是,式(6)可以寫成如下形式

(8)

3.2 NURBS參數(shù)化

風(fēng)力機葉片等幾何分基礎(chǔ)是將葉片NURBS化,利用Matlab完成風(fēng)力機葉片NURBS參數(shù)化。NURBS曲線擬合程序邏輯框圖如圖8所示。

圖8 葉片NURBS參數(shù)化邏輯框圖

在風(fēng)力機葉片設(shè)計中,一共得到41葉素的數(shù)據(jù)。因此,在求NURBS的控制點的時候,也得出了41個葉素的數(shù)據(jù)的控制點,限于篇幅有限,這里展示一個葉素的控制點部分?jǐn)?shù)據(jù),如表3所示。

表3 葉素1的控制點

圖9為葉片NURBS曲線與41個葉素的控制點曲線、原始曲線對比圖。有圖可知:三者的差別用肉眼幾乎看不出來,基本重合。圖10為三種曲線的細節(jié)信息。在圖10中,紅色表示NURBS曲線擬合成的葉素,綠色表示NURBS的控制曲線,藍色表示直接擬合成的風(fēng)力機葉片。因為 NURBS 基函數(shù)具有高階連續(xù)性,克服了有限元分析通常僅有C0連續(xù)性的弊端,所以用NURBS曲線擬合的葉素更加光滑,精確度更高。

圖9 葉素控制點曲線,原始曲線,NURBS曲線(整體展示)

圖10 葉素控制點曲線(綠),原始曲線(藍),NURBS曲線(紅)(一個葉素)

4 風(fēng)力機葉片等幾何分析

4.1 風(fēng)力機葉片T-spline模型

將控制點導(dǎo)入Rhino中,生成風(fēng)力機葉片NURBS模型,如圖11所示:

圖11 風(fēng)力機葉片NURBS模型

這種方法生成的模型具有很多的局限性,既不能實時保存成.iga文件讀取有關(guān)信息,同時也不能對網(wǎng)格進行細化。為了克服此問題,在Rhino中安裝了T-splines插件,實現(xiàn)CAD與CAE之間有機結(jié)合。力機葉片T-splines建模具體思路如下框圖12所示:

圖12 風(fēng)力機葉片T-splines模型創(chuàng)建邏輯框圖

采用相交曲線的方法,每個葉素的垂直距離為0.016mm,一共1000個平面與NURBS風(fēng)力機模型相交,如圖13所示?？梢缘玫?000多條的翼型曲線,其平面的數(shù)目和間距共同控制T-spline的縱向網(wǎng)格密度。

圖13 生成相交曲線

通過多次Loft T-spline surface操作,克服了網(wǎng)格偏轉(zhuǎn),網(wǎng)格錯亂等諸多網(wǎng)格生成過程中所遇到的困難,完成了風(fēng)力機葉片T-spline的模型,如圖14所示。與NURBS曲線生成的模型對比,葉片T-spline模型更加具體,細節(jié)描述的更加完整。同時得到了風(fēng)力機葉片的控制點,整個模型一共有3萬多個控制點和單元,如圖15所示。

圖14 風(fēng)機機葉片T-spline模型網(wǎng)格

圖15 風(fēng)力機葉片T-spline模型控制點

4.2 風(fēng)力葉片等幾何分析

4.2.1 等幾何分析平臺

Rhino中提供了各種“插件”滿足用戶的需求,它提供了一種圖形編程設(shè)計工具,用于參數(shù)設(shè)計,即為Grasshopper,同時可以利用免費和開源軟件開發(fā)包(SDK)對Rhino進行開發(fā)。Hsu[16]創(chuàng)建的平臺將等幾何分析的八大步驟分別集中于各個功能按鍵中,在利用每個功能的時候只需點擊對應(yīng)的圖標(biāo),如圖16所示。

圖16 Rhino分析平臺各功能按鍵

通過去除不重要的特征,曲面重建和不需要的結(jié),筆者完成了適合等幾何分析的風(fēng)力機T-spline的模型網(wǎng)格。圖17展示了風(fēng)力機葉片的T樣條表面網(wǎng)格,模型的自由度更高。

圖17 適合分析的風(fēng)力機葉片T-spline模型網(wǎng)格

4.2.2 增添材料屬性,邊界條件和載荷

完成了適合等幾何分析模型后,下一步需要增添風(fēng)力機葉片材料特性,用圖16 (b)中的“TsSelet”命令來選擇T-spline樣條模型的不同曲面元素,為不同的元素組定義不同的集合,在本文中將所有的曲面元素放在同一個集合中。在完成集合定義后,可以利用圖16(c)中的“Material setup”命令對集合賦予材料屬性。該功能可以允許選擇定義各向同性或者正交異性材料的復(fù)合層。邊界條件選擇在風(fēng)力機葉片的大端施加完全固定,如圖18所示,至于載荷,為了簡化計算,這里只考慮自身重力的影響,即g=9.8N/kg。

圖18 風(fēng)力機葉片的邊界條件(圖中紅色部分所示)

4.2.3 等幾何分析計算結(jié)果

將本文所設(shè)計的風(fēng)力機葉片T-spline模型導(dǎo)入IGA shell solver平臺中,進行動態(tài)仿真,其時間步長為0.001s,計算時間為1s,可以得到風(fēng)力機葉片計算前后對比。

經(jīng)過IGA可視化,得到如圖19所示的風(fēng)力葉片表面的位移大小,位移最大的位置出現(xiàn)在葉片的最右端,云圖最大值0.2413mm。

圖19 風(fēng)力機葉片等幾何分析表面位移云圖

同理可得圖20所示的風(fēng)力葉片應(yīng)變云圖,位移最大的位置出現(xiàn)在葉片彎折最嚴(yán)重出,應(yīng)力最大值為0.000282Mpa。

圖20 風(fēng)力機葉片等幾何分析應(yīng)變云圖

5 結(jié)果驗證

為了驗證基于該平臺算出的風(fēng)力機葉片結(jié)果的正確性,本文將41個葉素數(shù)據(jù)導(dǎo)入Solidworks中,利用放樣生成風(fēng)力機葉片的曲面模型,如圖21所示。

圖21 風(fēng)力機葉片曲面模型

在SolidWorks完成葉片模型構(gòu)建后,保存通用格式導(dǎo)入到傳統(tǒng)的有限元分析軟件Abaqus中,對該模型施加相同的材料屬性及邊界條件。利用Abaqus求出風(fēng)力機葉片位移云圖和應(yīng)變云圖,如圖22和圖23所示。在圖22中,葉片在在重力作用下的位移范圍為0mm～0.2561mm,利用等幾何方法算出的葉片位移范圍為0mm～0.2413mm,兩種方法得到的位移最大值偏差為6.133%。在圖23中,葉片在重力作用下應(yīng)變范圍為1.707e-8～2.997e-4Mpa,利用等幾何方法算出的葉片應(yīng)變范圍為0～2.82e-4Mpa。通過對比,兩種方法得到的應(yīng)變最大值偏差為6.277%,產(chǎn)生偏差的原因在于:在進行限元分析時,會進行網(wǎng)格劃分,此時會對特征進行簡化,使分析模型與實體模型產(chǎn)生差異,而等幾何分析模型與實體模型一致,因此兩種方法得出結(jié)論會存在較小的偏差,從而證明了將等幾何法應(yīng)用到風(fēng)力機葉片靜力學(xué)分析中的可行性與先進性。

圖22 葉片位移云圖

圖23 葉片應(yīng)變云圖

為了進一步驗證等幾何分析方法的正確性,在有限元模型中求出位移與應(yīng)變的最大值點隨時間變化的規(guī)律,與T-spline模型中的最大值點進行對比分析。導(dǎo)出這兩個點位移和應(yīng)變隨時間變化的數(shù)據(jù),如表4所示。

表4 時間與位移、應(yīng)變

將兩種不同的方法計算出的最大值的位移和應(yīng)變導(dǎo)入到一個圖中,如圖24所示。從圖中可知,兩種方法計算出的最大值點位移和應(yīng)變隨時間的變化規(guī)律近似為一條直線,變化規(guī)律相同。其中利用有限元計算出的位移曲線斜率k1=0.2561,等幾何斜率k2=0.2413,偏差在6.1%左右;利用有限元算出的應(yīng)變曲線斜率k1=0.000300,;利用等幾何方法計算出來的斜率k2=0.000282,偏差在6.4%左右。綜上所述,這兩種方法算出最大值點隨時間變化的趨勢相同,K值偏差保持在較小的范圍內(nèi),表明將等幾何方法應(yīng)用在風(fēng)力機葉片有限元分析上具有一定的可行性與先進性。

圖24 兩種方法的應(yīng)力及位移最大值隨時間變化對比

6 結(jié)論

1)利用三維葉片曲面集成表達式完成了風(fēng)力機葉片建模,獲得了41個葉素數(shù)據(jù)。

2)根據(jù)NURBS的性質(zhì),完成了NURBS擬合程序,結(jié)合41個葉素數(shù)據(jù),求出了風(fēng)力機葉片曲面的控制點,完成了風(fēng)力機葉片NURBS構(gòu)建。同時借用T-splines插件,完成風(fēng)力機葉片曲面的T-spline建模。

3)將風(fēng)力機葉片T-spline模型導(dǎo)入等幾何平臺中進行靜力學(xué)分析,得到葉片的位移與應(yīng)變云圖。同時將幾何模型導(dǎo)入Abaqus中進行有限元靜力學(xué)計算。將兩種方法得到的結(jié)果進行對比,最大應(yīng)力及最大位移偏差控制在6%左右。說明將等幾何分析應(yīng)用到風(fēng)力機葉片有限元分析具有一定的可行性與實用性。