逆風條件下柔性風帆性能的流固耦合研究

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

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

2 青島市海洋可再生能源重點實驗室，山東 青島 266100

0 引 言

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

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

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

1 數值模擬方法

1.1 流固耦合方法

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

1.2 流體運動控制

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

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

2 柔性風帆模型

2.1 形狀模型

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

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

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

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

表1 風帆幾何參數Table 1 Geometrical parameters of the sail

圖2 風帆形狀模型Fig. 2 Shape model of the sail

2.2 材料模型

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

3 柔性風帆的流固耦合分析

3.1 流固耦合分析的條件設置

基于AUTODYN 軟件采用Euler-Lagrange 作用的全耦合方法，設置以下條件進行柔性風帆非線性大變形分析。

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

2) 為了保證較高的流固耦合精度，流體域的網格尺寸應不大于結構的網格尺寸。結構的基本網格尺寸設置為60 mm，網格單元數量為6 988，計算域的網格數量為511 632。

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

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

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

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

6) 風帆被視為殼體結構，定義人工殼單元厚度為0.68 mm。

7) 在風帆的5 個截面上添加探測點，如圖4所示，通過輸出探測點的空間坐標數據，近似構建風帆變形的形狀。

圖4 帆面上的探測點Fig. 4 Gauges on sail surface

3.2 流固耦合結果分析

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

圖5 風帆變形的位移云圖Fig. 5 Displacement contours of sail deformation

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

4 柔性風帆的性能預測

4.1 STAR-CCM+計算的條件設置

為了說明柔性風帆變形對其空氣動力學性能的影響，采用STAR-CCM+在相同條件下對初始風帆和變形風帆進行性能預測。

1) 風帆被視為計算域內零厚度的內部面。

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

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

圖7 風帆周圍網格加密和邊界層Fig. 7 Refinements of the meshes around the sail and prism layers

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

4.2 性能預測結果分析

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

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

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

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

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

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

表3 初始風帆與變形風帆的升、阻力系數Table 3 Lift and drag coefficients of the initial and deformed sails

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

5 結 論

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

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

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

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

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

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

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