逆风条件下柔性风帆性能的流固耦合研究

方泽江，谭俊哲*,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) 柔性风帆变形的流固耦合研究适用于精确评估或者优化风帆性能。