汪 泉， 楊書益， 朱曉霜

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

0 引言

目前，有限元分析是研究風(fēng)力機(jī)葉片空氣動(dòng)力學(xué)性能最常見的計(jì)算方法，但在風(fēng)力機(jī)葉片有限元求解中會(huì)遇到幾個(gè)問題。 一是有限元求解精度與網(wǎng)格的細(xì)化程度有關(guān)，網(wǎng)格越細(xì)，計(jì)算結(jié)果精度越高，但計(jì)算時(shí)間和成本也隨之增加。 二是網(wǎng)格劃分工具對(duì)風(fēng)力機(jī)葉片等高級(jí)復(fù)雜曲面識(shí)別精度較低，無法精確劃分，會(huì)導(dǎo)致計(jì)算精度的缺失。三是當(dāng)風(fēng)力機(jī)葉片模型發(fā)生改變，必須重新劃分網(wǎng)格，耗費(fèi)大量時(shí)間。 基于以上原因，Hughes T J R 提出等幾何分析（IGA）的方法，該方法直接結(jié)合了CAD 中的幾何模型，將其中的幾何信息作為有限元分析的輸入信息，大大節(jié)省了劃分網(wǎng)格的時(shí)間[1]。IGA 是基于有限元分析方法的等參單元思想，是將CAD 中用于表達(dá)幾何模型的非均勻有理B 樣條的基函數(shù)作為形函數(shù)，實(shí)現(xiàn)了CAD 到CAE的無縫結(jié)合。

許多學(xué)者對(duì)IGA 進(jìn)行了研究。 陳明飛基于一階剪切變形理論，并采用等幾何有限元方法對(duì)任意曲率的功能梯度曲梁進(jìn)行了自由振動(dòng)分析[2]。王悅提出了一種基于T 樣條曲面的變網(wǎng)格柔性系統(tǒng)等幾何分析方法，形成了系統(tǒng)的T 樣條單元局部細(xì)化算法[3]。 榮吉利研究了IGA 中的HB 樣條基函數(shù)構(gòu)造以及多層貝塞爾提取方法問題， 并通過帶圓孔平板靜力學(xué)算例驗(yàn)證了基于HB 樣條IGA 方法的有效性和準(zhǔn)確性[4]。劉宏亮針對(duì)結(jié)構(gòu)等幾何優(yōu)化設(shè)計(jì)， 系統(tǒng)梳理和綜述了主要的等幾何優(yōu)化方法及其在結(jié)構(gòu)優(yōu)化設(shè)計(jì)中的應(yīng)用[5]。 Deng X W 將Kirchhoff-Love 殼單元用于3D 葉片建模，在基于Rhino 的Grasshopper 中的集成框架進(jìn)行建模和分析[6]。 Herrema A J 利用IGA 創(chuàng)建了一個(gè)優(yōu)化設(shè)計(jì)的框架， 用壓管為例證實(shí)了該框架的實(shí)用性[7]。 Aung T L 開發(fā)了集成IGA 的優(yōu)化設(shè)計(jì)框架， 并將其應(yīng)用于優(yōu)化兩個(gè)焊縫的形狀以減少應(yīng)力集中[8]。

然而以上研究大多基于簡(jiǎn)單模型的IGA，風(fēng)力機(jī)葉片作為一種高級(jí)復(fù)雜的曲面， 在等幾何方面的應(yīng)用還很少， 且目前的研究大多是借助三維軟件平臺(tái)的二次開發(fā)， 沒有實(shí)現(xiàn)真正意義上的設(shè)計(jì)和分析的統(tǒng)一。 為了實(shí)現(xiàn)復(fù)雜三維風(fēng)力機(jī)葉片設(shè)計(jì)與分析的完全統(tǒng)一， 避免復(fù)雜的網(wǎng)格劃分過程，提高其求解精度與效率，本文開展了風(fēng)力機(jī)葉片的IGA 研究。首先構(gòu)建適合IGA 的風(fēng)力機(jī)葉片NURBS 模型， 然后推導(dǎo)三維靜力學(xué)IGA 相關(guān)計(jì)算公式，實(shí)現(xiàn)三維復(fù)雜風(fēng)力機(jī)葉片的IGA。一方面通過簡(jiǎn)單算例精確解驗(yàn)證IGA 算法正確性，另一方面通過與商業(yè)軟件對(duì)比驗(yàn)證其優(yōu)越性。

1 風(fēng)力機(jī)葉片模型的構(gòu)建

在風(fēng)力機(jī)葉片設(shè)計(jì)和分析過程中， 設(shè)計(jì)師通常會(huì)先設(shè)計(jì)其CAD 模型，然后通過有限元分析軟件對(duì)其進(jìn)行網(wǎng)格劃分。 網(wǎng)格劃分的質(zhì)量決定了仿真分析的精度與效率， 特別對(duì)于風(fēng)力機(jī)葉片這種復(fù)雜的曲面，容易劃分低質(zhì)量的網(wǎng)格，從而影響分析結(jié)果的精度與效率。 為此本文基于NURBS 構(gòu)建了適合IGA 的風(fēng)力機(jī)葉片二維三維模型，模型設(shè)計(jì)和分析均采用同一的基函數(shù)， 不僅避免了復(fù)雜的網(wǎng)格劃分過程， 同時(shí)也消除了模型設(shè)計(jì)與分析之間的誤差。

1.1 風(fēng)力機(jī)葉片二維翼型的構(gòu)建

本文將NURBS 作為設(shè)計(jì)和分析基函數(shù)，連接了設(shè)計(jì)和分析。 一條p 次NURBS 曲線定義為

式中：c（u）為NURBS 樣條函數(shù)；pi為控制點(diǎn)；wi為權(quán)因子；Ni，p（u）為定義在非均勻節(jié)點(diǎn)u 上第i 個(gè)P次B 樣條基函數(shù)；a，b 為節(jié)點(diǎn)矢量u 的取值范圍。

NURBS 曲線具有幾何不變性、凸包性和局部支撐性等性質(zhì)， 是計(jì)算機(jī)輔助設(shè)計(jì)系統(tǒng)中最常見的表示方法。 采用NURBS 構(gòu)建風(fēng)力機(jī)葉片二維翼型可以提高翼型的光滑連續(xù)性，可以更準(zhǔn)確、更光滑地構(gòu)建風(fēng)力機(jī)葉片的翼型，為后續(xù)的IGA 提供準(zhǔn)確的模型支撐。 利用NURBS 來實(shí)現(xiàn)風(fēng)力機(jī)葉片翼型的參數(shù)化是翼型設(shè)計(jì)中至關(guān)重要的一步，也是IGA 的基礎(chǔ)。

在Matlab 中， 先編寫翼型反求控制點(diǎn)程序，得到二維翼型的控制點(diǎn)， 然后基于NURBS 理論編寫構(gòu)建翼型程序，其流程如圖1 所示。

圖1 NURBS 流程圖Fig.1 Flow chart of NURBS

首先選擇翼型的型值點(diǎn)， 然后通過參數(shù)化處理確定樣條的階數(shù)，求得節(jié)點(diǎn)矢量信息，再通過切矢邊界條件，求解線性方程組A*D=E，獲得控制頂點(diǎn)向量D， 接著將得到的控制點(diǎn)和節(jié)點(diǎn)矢量及權(quán)因子作為輸入?yún)?shù)，得到NURBS 構(gòu)建的翼型。圖2 為DU97-W-300 翼型分布圖。

圖2 DU97-W-300 翼型分布圖Fig.2 DU97-W-300 airfoil

選取了部分型值點(diǎn)，求得其控制點(diǎn)為

1.2 風(fēng)力機(jī)葉片三維模型的構(gòu)建

葉片是以氣動(dòng)中心軸進(jìn)行扭轉(zhuǎn)的， 翼型的氣動(dòng)中心為1/4 弦長(zhǎng)位置， 一般變弦長(zhǎng)變扭角的三維葉片集成表達(dá)式為

式中：θ 為幅角；u 為葉片展向位置；ρ（θ）為翼型分布函數(shù)；a 為擬合系數(shù)；xM為翼型氣動(dòng)中心展向位置，一般取0.25；c（u）為弦長(zhǎng)分布函數(shù)；β（u）為扭角分布函數(shù)。

根據(jù)三維葉片的形函數(shù)/分布函數(shù)公式將其耦合到NURBS 建模中。 分別設(shè)置21 個(gè)翼型的數(shù)據(jù)，每個(gè)翼型由24 個(gè)控制點(diǎn)構(gòu)成，其中12 個(gè)控制上翼型，12 個(gè)控制下翼型。

圖3 為基于NURBS 構(gòu)建的三維風(fēng)力機(jī)葉片模型。模型由21 個(gè)翼型組成，每個(gè)翼型由24 個(gè)控制點(diǎn)組成。

圖3 風(fēng)力機(jī)葉片NURBS 模型Fig.3 NURBS model of wind turbine blade

2 風(fēng)力機(jī)葉片的IGA

2.1 IGA 方法

設(shè)定風(fēng)力機(jī)葉片的材料為各向同性材料，風(fēng)力機(jī)葉片的IGA 屬于線彈性問題，風(fēng)力機(jī)葉片在載荷和邊界條件下的位移可以表示為邊界值問題，其位移ui為

式中：Γ 為域邊界；σ 為應(yīng)力。

式（7）為狄利克雷邊界條件，式（8）為紐曼邊界條件。

IGA 是對(duì)上述式進(jìn)行變分后得到弱形式，對(duì)于邊界條件來說，u∈s，v∈V 滿足以下條件：

式中：K 為剛度矩陣；P 為載荷向量。

采用高斯積分方法，在得到了位移值之后，可以求得單元的應(yīng)變?chǔ)?為

式中：L 為微分算子。

由此可得應(yīng)力σ 為

式中：D 為彈性矩陣，與材料的泊松比和彈性模量有關(guān)。

IGA 流程如圖4 所示。

圖4 IGA 流程圖Fig.4 IGA flow chart

2.2 NURBS 算例

為了驗(yàn)證NURBS 三維實(shí)體IGA 的正確性，首先對(duì)簡(jiǎn)單的算例進(jìn)行IGA 并與其精確解對(duì)比?？紤]到受均勻拉力作用下的球形孔洞問題， 算例模型結(jié)構(gòu)尺寸及相關(guān)參數(shù)如圖5 所示。

圖5 算例模型結(jié)構(gòu)尺寸Fig.5 Model structure size of the calculation example

圖中各參數(shù)均為無量綱參數(shù)， 其在球坐標(biāo)系（r，?，θ）的精確解[9]為

式中：a 為球體的半徑；S 為施加的拉力；ν 為泊松比。

由于該算例的對(duì)稱性，本文考慮對(duì)模型的1/8 進(jìn)行IGA，并將計(jì)算得到的結(jié)果與精確解進(jìn)行對(duì)比。 當(dāng)S=1，ν=0.3 時(shí)，σmax=3S（9-5ν）/（14-10ν），其最大應(yīng)力精確解為2.045 45。 圖6 為球形孔洞的應(yīng)力云圖。 圖中IGA 計(jì)算得到最大應(yīng)力值為2.05， 而其精確解為2.045 45， 其相對(duì)誤差為0.22%。 求解結(jié)果與精確解相差無幾，通過此算例驗(yàn)證了本文構(gòu)造的IGA 程序的正確性，為后文風(fēng)力機(jī)葉片的IGA 提供了有效的理論依據(jù)。

圖6 球形孔洞的應(yīng)力云圖Fig.6 Stress cloud diagram of spherical holes

2.3 風(fēng)力機(jī)葉片的IGA

在1.2 節(jié)中構(gòu)建了適合IGA 的三維風(fēng)力機(jī)葉片模型，考慮了葉片的彎扭特性，根據(jù)文獻(xiàn)[10]，為了方便計(jì)算， 將實(shí)際風(fēng)力機(jī)葉片模型進(jìn)行適當(dāng)簡(jiǎn)化，將其材料定義為各向同性材料，材料參數(shù)設(shè)置為彈性模量為1×105MPa，泊松比為0.3，密度為7.8×103kg/m3。

風(fēng)力機(jī)葉片在正常工作時(shí)， 其實(shí)際受載情況為揮舞方向、拍打方向和扭矩的耦合。本文考慮在風(fēng)力機(jī)葉片揮舞方向受載這一工況下的三維風(fēng)力機(jī)葉片的等幾何靜力學(xué)分析， 根據(jù)揮舞方向?qū)嶋H受載情況，其具體約束設(shè)置為每個(gè)節(jié)點(diǎn)的Y，Z 方向固定，X 方向受位移約束-1 mm （圖7 中載荷B），將風(fēng)力機(jī)葉片的葉根位置固定。

圖7 風(fēng)力機(jī)葉片的模型Fig.7 Model of wind turbine blade

將求解結(jié)果輸出為.vtu 文件格式， 然后通過將風(fēng)力機(jī)葉片IGA 結(jié)果導(dǎo)入開源軟件Paraview中進(jìn)行可視化，得到了風(fēng)力機(jī)葉片的位移云圖（圖8），在該載荷作用下，其最大位移為2.929 9 mm，最小位移為0.027 5 mm。

圖8 風(fēng)力機(jī)葉片位移云圖Fig.8 Wind turbine blade displacement cloud diagram

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

為了驗(yàn)證風(fēng)力機(jī)葉片IGA 結(jié)果的正確性，將計(jì)算所得風(fēng)力機(jī)葉片模型控制點(diǎn)信息通過插入曲線命令導(dǎo)入三維建模軟件SolidWorks 中進(jìn)行建模，然后將模型導(dǎo)入軟件ANSYS 中構(gòu)建其有限元模型（圖9）。 將IGA 的結(jié)果和與有限元分析的結(jié)果進(jìn)行對(duì)比，以驗(yàn)證風(fēng)力機(jī)葉片IGA 的正確性。

圖9 有限元模型Fig.9 Finite element model

在商業(yè)軟件ANSYS 中施加與IGA 相同的材料屬性、 邊界條件和載荷約束， 并進(jìn)行求解（圖10）。

圖10 有限元結(jié)果Fig.10 Finite element results

為了驗(yàn)證IGA 的正確性，將風(fēng)力機(jī)葉片在相同的條件下進(jìn)行了有限元分析，葉片在載荷的作用下最大位移為3.012 4 mm，最小位移為0.028 7 mm，IGA 與有限元分析的最大誤差為5.1%。 為了進(jìn)一步驗(yàn)證IGA 的正確性，分別取固定位置處的節(jié)點(diǎn)，分析對(duì)比IGA 和有限元分析的結(jié)果（表1）。從不同位置節(jié)點(diǎn)位移的對(duì)比結(jié)果可以看出，IGA與有限元分析的結(jié)果誤差均在5.1%以內(nèi)，進(jìn)一步驗(yàn)證了IGA 的準(zhǔn)確性。

表1 載荷下的位移Table 1 Displacement under load mm

為了驗(yàn)證IGA 的優(yōu)越性，將不同的控制點(diǎn)數(shù)目情況下IGA 的結(jié)果與ANSYS 結(jié)果進(jìn)行了對(duì)比（圖11）。

圖11 兩種方法不同節(jié)點(diǎn)數(shù)目的位移值對(duì)比Fig.11 Comparison of displacement values of the two methods with different numbers of nodes

當(dāng)節(jié)點(diǎn)數(shù)為600 時(shí)，IGA 的計(jì)算結(jié)果為2.8831，計(jì)算時(shí)間為42 s； 當(dāng)ANSYS 節(jié)點(diǎn)數(shù)量為2 300時(shí)，計(jì)算時(shí)間為151 s，其計(jì)算結(jié)果與IGA 計(jì)算結(jié)果基本一致。 這也證明了在不損失計(jì)算結(jié)果精度的條件下，IGA 方法較有限元方法節(jié)省了大量的節(jié)點(diǎn)，提高了求解的效率。

4 結(jié)論

本文針對(duì)復(fù)雜三維風(fēng)力機(jī)葉片， 首先構(gòu)建了適合IGA 的風(fēng)力機(jī)葉片NURBS 模型， 然后推導(dǎo)三維靜力學(xué)IGA 相關(guān)計(jì)算公式，并編寫了相關(guān)程序，實(shí)現(xiàn)了三維復(fù)雜風(fēng)力機(jī)葉片的IGA。并通過簡(jiǎn)單算例精確解和商業(yè)軟件對(duì)其進(jìn)行驗(yàn)證， 得到以下結(jié)論。

①構(gòu)建IGA 框架， 并編寫相關(guān)計(jì)算程序，通過簡(jiǎn)單算例的精確解驗(yàn)證了其正確性。

②將風(fēng)力機(jī)葉片變弦長(zhǎng)變扭角公式與NURBS 相耦合， 構(gòu)建了適合IGA 的三維風(fēng)力機(jī)葉片模型， 并對(duì)構(gòu)建的模型進(jìn)行IGA， 并將IGA結(jié)果與商業(yè)軟件進(jìn)行對(duì)比。對(duì)比結(jié)果表明：位移最大誤差為5.1%， 驗(yàn)證了風(fēng)力機(jī)葉片IGA 的正確性； 對(duì)兩種方法在不同節(jié)點(diǎn)數(shù)情況下的計(jì)算結(jié)果進(jìn)行對(duì)比，在不損失求解精度的情況下，IGA 較有限元分析節(jié)省了大量的節(jié)點(diǎn)，提高了求解的效率。