水平軸風(fēng)力機(jī)葉片曲面建模與等幾何分析研究

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

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

1 引言

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

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

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

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

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

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

(1)

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

xM—?dú)鈩?dòng)中心展向位置

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

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

表1 翼型位置分布

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

該葉片的弦長(zhǎng)及扭角分布如圖2和3所示。

圖2 葉片弦長(zhǎng)分布

圖3 葉片扭角分布

根據(jù)葉片弦長(zhǎng)、扭角分布數(shù)據(jù),可擬合出葉片弦長(zhǎng)及扭角高階多項(xiàng)式系數(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)力機(jī)葉片復(fù)雜曲面模型構(gòu)建

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

圖4 葉片三維曲面圖

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

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

3.1 NURBS曲線

P次B樣條曲線的定義可以用如下公式來(lái)進(jìn)行定義

(4)

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

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

{pi}—控制點(diǎn)

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

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

(5)

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

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

表2 翼型控制點(diǎn)

圖5 NACA63415翼型廓線

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

圖6 翼型基函數(shù)

圖7 B樣條曲線翼型

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

(6)

式中:{pi}—控制點(diǎn);

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

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

同時(shí),為了簡(jiǎn)化式(6),令:

(7)

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

(8)

3.2 NURBS參數(shù)化

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

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

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

表3 葉素1的控制點(diǎn)

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

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

圖10 葉素控制點(diǎn)曲線(綠),原始曲線(藍(lán)),NURBS曲線(紅)(一個(gè)葉素)

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

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

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

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

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

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

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

圖13 生成相交曲線

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

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

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

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

4.2.1 等幾何分析平臺(tái)

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

圖16 Rhino分析平臺(tái)各功能按鍵

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

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

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

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

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

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

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

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

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

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

圖20 風(fēng)力機(jī)葉片等幾何分析應(yīng)變?cè)茍D

5 結(jié)果驗(yàn)證

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

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

在SolidWorks完成葉片模型構(gòu)建后,保存通用格式導(dǎo)入到傳統(tǒng)的有限元分析軟件Abaqus中,對(duì)該模型施加相同的材料屬性及邊界條件。利用Abaqus求出風(fēng)力機(jī)葉片位移云圖和應(yīng)變?cè)茍D,如圖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。通過(guò)對(duì)比,兩種方法得到的應(yīng)變最大值偏差為6.277%,產(chǎn)生偏差的原因在于:在進(jìn)行限元分析時(shí),會(huì)進(jìn)行網(wǎng)格劃分,此時(shí)會(huì)對(duì)特征進(jìn)行簡(jiǎn)化,使分析模型與實(shí)體模型產(chǎn)生差異,而等幾何分析模型與實(shí)體模型一致,因此兩種方法得出結(jié)論會(huì)存在較小的偏差,從而證明了將等幾何法應(yīng)用到風(fēng)力機(jī)葉片靜力學(xué)分析中的可行性與先進(jìn)性。

圖22 葉片位移云圖

圖23 葉片應(yīng)變?cè)茍D

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

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

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

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

6 結(jié)論

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

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

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