逆風(fēng)條件下柔性風(fēng)帆性能的流固耦合研究

方澤江，譚俊哲*,2，紀(jì)光英，袁鵬,2，孫澤，王樹杰,2

1 中國海洋大學(xué) 工程學(xué)院，山東 青島 266100

2 青島市海洋可再生能源重點(diǎn)實(shí)驗(yàn)室，山東 青島 266100

0 引 言

風(fēng)帆性能與帆船動力性能直接相關(guān)。柔性風(fēng)帆由薄織物制成，航行時呈現(xiàn)為帶有拱度的形狀，捕獲風(fēng)能，以產(chǎn)生帆船航行的驅(qū)動力。柔性風(fēng)帆在風(fēng)載荷下易發(fā)生非線性大變形，它的初始形狀與航行時的真實(shí)形狀存在很大差異。柔性風(fēng)帆航行時的真實(shí)形狀可簡稱為飛行形狀[1]，風(fēng)帆的幾何形狀影響著周圍的空氣流動和產(chǎn)生的力。有必要研究柔性風(fēng)帆的變形并采集其飛行形狀[2]，更準(zhǔn)確地預(yù)測柔性風(fēng)帆的空氣動力學(xué)性能。

計算流體力學(xué)（CFD）能夠幫助工程師完成最佳風(fēng)帆設(shè)計[3]。用CFD 預(yù)測風(fēng)帆性能，絕大多數(shù)以風(fēng)帆初始形狀為形狀模型[4]，風(fēng)帆形狀不隨風(fēng)況變化而改變[5-6]。柔性風(fēng)帆變形是非線性大變形問題，與風(fēng)載荷之間的相互作用屬于復(fù)雜的流固耦合現(xiàn)象。Lee 等[7]采用滑動網(wǎng)格方法對游艇風(fēng)帆的二維截面進(jìn)行流固耦合分析，研究風(fēng)帆變形對性能的影響。Bak 等[8]采用單向流固耦合方法研究了無拱度的織物材料的簡易風(fēng)帆變形的非線性行為，如褶皺現(xiàn)象。Bak 等[9-10]將風(fēng)帆視為殼體結(jié)構(gòu)，基于集成的CFX和LS-DYNA對30 ft（9.144 m）游艇風(fēng)帆進(jìn)行部分雙向流固耦合分析，比較了初始風(fēng)帆和變形風(fēng)帆的空氣動力學(xué)性能；還研究了風(fēng)帆和桅桿變形對風(fēng)帆空氣動力學(xué)性能的影響。柔性風(fēng)帆變形的流固耦合試驗(yàn)研究也取得了許多進(jìn)展[11-14]。

柔性風(fēng)帆性能試驗(yàn)成本高昂、試驗(yàn)條件受限；而CFD 數(shù)值模擬成本低、操作方便，但難以準(zhǔn)確解決柔性風(fēng)帆非線性大變形問題。為此，本文將對逆風(fēng)條件下30 ft 游艇主帆的變形進(jìn)行流固耦合研究。采用顯式有限元分析程序AUTODYN解決柔性風(fēng)帆非線性大變形問題，構(gòu)建變形后的柔性風(fēng)帆形狀。采用STAR-CCM+軟件預(yù)測柔性風(fēng)帆空氣動力學(xué)性能，比較柔性風(fēng)帆初始形狀和變形形狀表面的壓力分布、周圍的空氣流動情況、產(chǎn)生的升力和阻力，以及受力中心位置之間的差異，討論柔性風(fēng)帆變形對空氣動力學(xué)性能的影響。

1 數(shù)值模擬方法

1.1 流固耦合方法

采用Euler-Lagrange 作用的流固耦合方法。在全耦合中，拉格朗日網(wǎng)格和歐拉網(wǎng)格進(jìn)行動態(tài)相互作用。流固耦合的結(jié)構(gòu)模型使用需定義人工厚度的殼單元，人工厚度至少是其周圍歐拉單元最小尺寸的2 倍。對于基于有限元方法（FEM）的柔性風(fēng)帆非線性大變形分析，在ANSYS Workbench 中通過集成的顯式動力學(xué)模塊Explicit Dynamics 和AUTODYN 進(jìn)行流固耦合研究。

1.2 流體運(yùn)動控制

風(fēng)帆周圍的湍流使用STAR-CCM+通過求解三維Navier-Stokes (N-S)方程組模擬。質(zhì)量守恒方程、動量守恒方程和能量守恒方程分別為[15]：

對N-S 方程進(jìn)行時間積分，獲得雷諾平均Navier-Stokes (RANS)方程。湍流運(yùn)動計算時，還需附加湍流方程。本文研究采用SSTk-ω 湍流模型，因?yàn)镾STk-ω 考慮了湍流中剪切應(yīng)力的輸運(yùn)，可以精確計算在逆壓梯度下從光滑表面流動分離的時刻和數(shù)量，能準(zhǔn)確模擬風(fēng)帆周圍的氣流分離，計算成本比更高階湍流模型低[16]。

2 柔性風(fēng)帆模型

2.1 形狀模型

對Yoo 等[17]設(shè)計的30 ft 游艇主帆進(jìn)行數(shù)值模擬。風(fēng)帆的拱度形狀通過零厚度NACA 6 位數(shù)機(jī)翼截面a=0.8 的中弧線表示，中弧線通過載荷分布線得到，從x/c=0 至x/c=a載荷均勻分布，從x/c=a至x/c=1 載荷線性減小至0，函數(shù)表達(dá)式如式(4)所示[18]，Lee 等[7]，Bak 等[9]，Yoo 等[19]，Ciortan等[20]， Bak 等[10,21]都采用此風(fēng)帆模型進(jìn)行了數(shù)值模擬或者比例模型試驗(yàn)研究。

式中：y為中弧線的縱坐標(biāo)；c為弦長；x為沿弦長方向的距離；a=0.8，為中弧線類型的標(biāo)記；cli=1，為設(shè)計升力系數(shù)；g和l的表達(dá)式為

NACAa=0.8 的中弧線如圖1 所示，采用5 個機(jī)翼截面的中弧線通過放樣構(gòu)建出風(fēng)帆形狀，5 個截面從風(fēng)帆底邊沿高度方向依次命名為foot，25%，50%，75%，top，風(fēng)帆幾何參數(shù)如表1 所示，幾何形狀如圖2 所示，帆高11.9 m，面積24.7 m2。

圖1 NACA a=0.8 的中弧線Fig. 1 Mean camber line with NACA a=0.8

表1 風(fēng)帆幾何參數(shù)Table 1 Geometrical parameters of the sail

圖2 風(fēng)帆形狀模型Fig. 2 Shape model of the sail

2.2 材料模型

常用薄織物作為柔性風(fēng)帆材料，假設(shè)材料為各向同性(Isotropic)彈性材料，材料模型采用密度ρ、彈性模量E和泊松比μ來表征。Kevlar?49 織物的材料特性如表2 所示[9]。

3 柔性風(fēng)帆的流固耦合分析

3.1 流固耦合分析的條件設(shè)置

基于AUTODYN 軟件采用Euler-Lagrange 作用的全耦合方法，設(shè)置以下條件進(jìn)行柔性風(fēng)帆非線性大變形分析。

1) 以圖2 所示風(fēng)帆形狀作為柔性風(fēng)帆的初始形狀，表2 所示Kevlar?49 織物作為柔性風(fēng)帆材料。

2) 為了保證較高的流固耦合精度，流體域的網(wǎng)格尺寸應(yīng)不大于結(jié)構(gòu)的網(wǎng)格尺寸。結(jié)構(gòu)的基本網(wǎng)格尺寸設(shè)置為60 mm，網(wǎng)格單元數(shù)量為6 988，計算域的網(wǎng)格數(shù)量為511 632。

3) 相對風(fēng)向角為20°，風(fēng)速為10.3 m/s。

4) 終止時間設(shè)置：流體域長度為4.3 m，風(fēng)速10.3 m/s，時間t=4.3/10.3=0.417 s，考慮風(fēng)和帆的耦合作用需要耗費(fèi)時間，所以取終止時間為0.5 s。

5) 柔性風(fēng)帆的固定方案對變形具有很大影響，實(shí)際上主前緣固定在桅桿上，后角由橫桿末端的拉力支撐，因此在前緣和后角添加固定邊界條件約束所有移動，以近似代替桅桿和橫桿對主帆的固定和支撐，如圖3 所示。

圖3 固定邊界條件Fig. 3 Fixed support conditions

6) 風(fēng)帆被視為殼體結(jié)構(gòu)，定義人工殼單元厚度為0.68 mm。

7) 在風(fēng)帆的5 個截面上添加探測點(diǎn)，如圖4所示，通過輸出探測點(diǎn)的空間坐標(biāo)數(shù)據(jù)，近似構(gòu)建風(fēng)帆變形的形狀。

圖4 帆面上的探測點(diǎn)Fig. 4 Gauges on sail surface

3.2 流固耦合結(jié)果分析

圖5 所示為風(fēng)帆變形的位移云圖，最大變形約為26 mm，在沒有固定邊界條件的風(fēng)帆邊緣變形較大。圖6 所示為風(fēng)帆變形前、后截面的拱度曲線，可見風(fēng)帆前緣附近拱度向背風(fēng)面變形更深，風(fēng)帆后緣附近帆高較低截面的拱度向背風(fēng)面靠近，帆高較高截面的拱度向迎風(fēng)面靠近，風(fēng)帆在風(fēng)載荷作用下發(fā)生明顯扭曲現(xiàn)象。風(fēng)帆初始形狀面積為24.7 m2，變形形狀面積為24.3 m2，面積減少約1.6%，說明風(fēng)帆面積變化很小，風(fēng)帆在風(fēng)載荷下的拉伸對風(fēng)帆變形的影響可以忽略不計。

圖5 風(fēng)帆變形的位移云圖Fig. 5 Displacement contours of sail deformation

圖6 初始風(fēng)帆和變形風(fēng)帆拱度曲線比較Fig. 6 A comparison on camber lines of the initial and deformed sails

4 柔性風(fēng)帆的性能預(yù)測

4.1 STAR-CCM+計算的條件設(shè)置

為了說明柔性風(fēng)帆變形對其空氣動力學(xué)性能的影響，采用STAR-CCM+在相同條件下對初始風(fēng)帆和變形風(fēng)帆進(jìn)行性能預(yù)測。

1) 風(fēng)帆被視為計算域內(nèi)零厚度的內(nèi)部面。

2) 假設(shè)溫度為25℃，空氣密度為1.169 kg/m3，湍流邊界層中第1 層網(wǎng)格與帆面之間的距離y+=118。風(fēng)為縱向風(fēng)，改變風(fēng)帆位置，使相對風(fēng)向角為20°，風(fēng)速設(shè)置為10.3 m/s。

3) 采用多面體網(wǎng)格，加密風(fēng)帆附近的網(wǎng)格并設(shè)置邊界層網(wǎng)格(圖7)，初始風(fēng)帆和變形風(fēng)帆的網(wǎng)格數(shù)量分別為36 893 和35 671，計算域的網(wǎng)格數(shù)量分別為3 557 850 和3 535 506。

圖7 風(fēng)帆周圍網(wǎng)格加密和邊界層Fig. 7 Refinements of the meshes around the sail and prism layers

4) 采用RANS 方程和k-ω SST 湍流模型。

4.2 性能預(yù)測結(jié)果分析

比較初始風(fēng)帆和變形風(fēng)帆表面的壓力分布、周圍的空氣流動情況、產(chǎn)生的升力和阻力、受力中心位置的變化。

圖8 所示為初始風(fēng)帆和變形風(fēng)帆表面的壓力分布，變形風(fēng)帆迎風(fēng)面壓力值有所增大，前緣附近正壓值的范圍也有所擴(kuò)大，這解釋了如圖6 所示的風(fēng)帆前緣附近拱度向背風(fēng)面靠近。圖9 所示為初始風(fēng)帆在帆高25%截面處的壓力分布和速度流線，在背風(fēng)面后緣附近出現(xiàn)流動分離現(xiàn)象，流動分離導(dǎo)致失速使得背風(fēng)面后緣附近壓力值升高，約為?30 Pa，而迎風(fēng)面后緣附近壓力值約為5 Pa，因此帆高25%截面后緣附近的壓差為正，這解釋了如圖6 所示的帆高25% 截面后緣附近的拱度仍然向背風(fēng)面靠近。圖10 所示為初始風(fēng)帆和變形風(fēng)帆在帆高25%截面的速度分布，變形風(fēng)帆較初始風(fēng)帆在背風(fēng)面后緣附近產(chǎn)生的流動分離現(xiàn)象更嚴(yán)重。

圖8 初始風(fēng)帆和變形風(fēng)帆的帆面壓力分布Fig. 8 Pressure distribution on the surfaces of the initial and deformed sails

圖9 初始風(fēng)帆在帆高25%截面的壓力分布和速度流線Fig. 9 Pressure distribution and velocity lines around the initial sail on the section of 25% sail height

圖10 初始風(fēng)帆和變形風(fēng)帆在帆高25%截面的速度分布Fig. 10 Velocity distribution around the initial and deformed sails on the section of 25% sail height

通過比較初始風(fēng)帆和變形風(fēng)帆的升力系數(shù)（CL）和阻力系數(shù)（CD），分析風(fēng)帆變形對空氣動力學(xué)性能的影響。由表3 可知，與初始風(fēng)帆相比，變形風(fēng)帆的升力系數(shù)減小約0.98%，阻力系數(shù)增大約2.7%，表征空氣動力學(xué)性能的升阻比減小約3.6%。該結(jié)果與Bak 等[9]通過流固耦合研究的三維風(fēng)帆變形對性能的影響效果一致。分析表明，產(chǎn)生上述現(xiàn)象的主要原因在于：帆高較低的帆面所占面積較大，是決定風(fēng)帆性能的主要部分，在相對風(fēng)向角為20°的逆風(fēng)條件下，變形風(fēng)帆較低帆面的背風(fēng)面后緣附近產(chǎn)生了更嚴(yán)重的流動分離現(xiàn)象(圖10)，較低帆面后緣附近的拱度向背風(fēng)側(cè)靠近(圖6)，導(dǎo)致有效攻角減小。因此在該條件下拱度增大帶來的升力增大效應(yīng)相對較小，總體升力減小。

表3 初始風(fēng)帆與變形風(fēng)帆的升、阻力系數(shù)Table 3 Lift and drag coefficients of the initial and deformed sails

帆船在穩(wěn)定航行時，橫傾力矩等于扶正力矩，風(fēng)帆變形時其受力中心位置會發(fā)生變化，從而使帆船失去平衡。假設(shè)風(fēng)帆橫向投影面的質(zhì)心為風(fēng)帆受力中心[22]，采用SolidWorks 評估風(fēng)帆初始形狀和變形形狀模型的質(zhì)心，變形風(fēng)帆的受力中心向右下方輕微移動。

5 結(jié) 論

本文采用非線性有限元方法計算柔性風(fēng)帆在風(fēng)載荷作用下的變形，以固定邊界條件代替桅桿和橫桿對風(fēng)帆的支撐和固定作用，在帆面上添加探測點(diǎn)，通過輸出探測點(diǎn)的空間坐標(biāo)數(shù)據(jù)構(gòu)建變形風(fēng)帆的形狀，比較初始風(fēng)帆和變形風(fēng)帆的形狀特性；采用STAR-CCM+在相同條件下預(yù)測初始風(fēng)帆和變形風(fēng)帆的性能。主要結(jié)論如下：

1) 變形風(fēng)帆較初始風(fēng)帆面積減小約1.6%，說明材料本身的拉伸對拱度變化的影響可忽略不計。

2) 變形風(fēng)帆在帆面前緣附近拱度向背風(fēng)面靠近；在后緣附近帆高較低截面的拱度向背風(fēng)面靠近，帆高較高截面的拱度向迎風(fēng)面靠近，有明顯扭曲現(xiàn)象。

3) 與初始風(fēng)帆相比，變形風(fēng)帆的迎風(fēng)面整體壓力值增大。

4) 與初始風(fēng)帆相比，變形風(fēng)帆的升力系數(shù)減小約0.98%，阻力系數(shù)增大約2.7%，升阻比減小約3.6%。主要原因在于，變形風(fēng)帆較低帆面的背風(fēng)面后緣附近產(chǎn)生了更嚴(yán)重的流動分離現(xiàn)象，導(dǎo)致失速，使得分離區(qū)域的壓強(qiáng)增加，從而使風(fēng)帆兩側(cè)的壓強(qiáng)差降低，進(jìn)而引起升力減小；較低帆面后緣附近的拱度向背風(fēng)面靠近，造成有效攻角減小,因此風(fēng)帆的整體升力減小。

5) 變形風(fēng)帆的受力中心向右下方輕微移動。

6) 柔性風(fēng)帆變形的流固耦合研究適用于精確評估或者優(yōu)化風(fēng)帆性能。