基于单向流固耦合的不同攻角下奥运会帆板帆翼空气动力特性数值模拟

2021-03-30贺阳映蔺世杰

中国体育科技 2021年1期

贺阳映，马勇，张松，蔺世杰

帆翼是帆船帆板前进的动力源，其性能直接影响运动帆船帆板航行的速度（Parolini et al.，2005；Persson et al.，2017；Viola，2011）。帆板运动由冲浪运动演变而来（马勇等，2013a；Ma et al.，2016），其比赛的基本航线可以分为梯形外绕航线和梯形内绕航线，整个航线涵盖了迎风、横风、顺风等多个航段（蔺世杰等，2017；马勇等，2013c），运动员需要进行调帆、换舷等操作来完成超越、绕标等动作（马勇等，2013a，2013c，2016）。在比赛中，不同航段中的每一次操作都会改变帆翼的攻角，而帆翼攻角的改变会影响帆翼的升力与阻力，因此，需要对不同攻角下帆板帆翼的空气动力性能进行研究（马勇等，2016；Ma et al.，2016）。对于奥运会级别的 Neil Pryde RS：X帆板，比赛中帆翼的攻角会主动调整或者受风作用改变，帆翼会发生形变，而帆翼外形的变化又影响帆翼空气动力性能，也就是帆翼与其周围流场会发生流固耦合（Fluid-Structure Interaction，FSI）作用（雷晓珊等，2019；Augier，2012；Bak et al.，2013；Durand et al.，2014），因此，对于帆翼空气动力性能研究需要考虑FSI。

越来越多的学者利用计算流体动力学（Computation‐al Fluid Dynamic，CFD）的方法研究运动帆翼的气动特性（Parolini et al.，2005；Viola 2011；Lee et al.，2016；Persson et al.，2017）。在试验研究和求解Navier-Stokes方程数值模拟方面，Giacobone（2017）对“美洲杯”帆船在不同湍流模型（Spalart-Allmaras、k-ε、SSTk-ω）和不同网格下的升阻力系数、力矩系数进行了计算，并将计算结果与风洞试验结果进行对比。Nava等（2017）基于雷诺时均方程（Reynolds Averaged Navier-Stokes，RANS）和大涡模拟（Large Eddy Simulation，LES），对迎风条件下速度为7.4m/s、攻角为18.7°的美洲杯帆船前帆和主帆的空气动力性能进行了数值模拟，所得结果与已发表的试验数据吻合良好。在Neil Pryde RS：X级别帆板帆翼空气动力性能研究方面，在不考虑流固耦合情况下，何海峰（2012）和罗晓川（2012）基于RANS方法进行了不同攻角和不同风速下的数值模拟，发现帆翼的失速角为20°～30°之间，随着攻角的增加，帆面压力中心向前移动至帆翼前方约1/3处。

随着计算机水平的发展与提升，考虑流固耦合效应的帆翼空气动力性能研究越来越多。Lee等（2016）对二维的三面水翼帆0°～360°（间隔为15°）攻角范围空气动力性能进行数值模拟，发现 45°、90°和 135°3 个常见视风角下三面水翼帆的最佳攻角范围。Sacher等（2017）考虑到FSI的IMOCA主帆空气动力性能数值模拟与试验较吻合，认为通过高斯方程可以将帆翼基本参数和性能参数进行快速优化。Deparday等（2018）对风角为50°～140°下的三角帆帆翼进行了流固耦合研究，发现帆翼气动弹性相对较弱，非稳定性振动发生在前帆。雷晓珊等（2019）应用弹簧光顺和局部网格重构的动网格技术，进行了同一攻角时不同风速下Neil Pryde RS：X级别帆板帆翼的单/双向流固耦合数值模拟，得到了帆翼气动力、表面压力、变形及应变特征。

从国内外研究进展来看，目前对于运动帆翼的研究主要是不考虑流固耦合时帆翼气动特性的数值模拟与试验，而对Neil Pryde RS：X级别帆板帆翼流固耦合的问题也仅讨论了风速对于帆翼气动特性的影响，忽略了帆翼在不同攻角时风致涡激振动对帆翼结构的影响。因此，本文拟通过求解RANSE，基于单向流固耦合，对风速为6m/s时迎风、横风、顺风航段常见攻角下Neil Pryde RS：X帆板帆翼的流场及结构场进行数值模拟，分析攻角变化对帆翼周围流场及帆翼结构的影响。

1 研究对象与方法

1.1 研究对象

研究对象是奥运会级别Neil Pryde RS：X帆板帆翼。首先，利用Pro/ENGINEER Wildfire对测量得到的帆翼外型数据进行逆向简化建模，忽略索具、加强筋等部分，只考虑桅杆、帆面。帆翼模型如图1和图2所示，NP帆板帆翼的基本参数如表1所示。Neil Pryde RS：X帆板帆翼的失速角处于20°～30°之间（何海峰，2012；罗晓川，2012），但是由于所选攻角间隔较大，计算得到的帆板帆翼失速角范围也过于宽泛，因此，本文对该帆翼在风速为6 m/s、攻角为 20°～50°及80°～100°、间隔5°时的工况进行数值仿真。

图1 Neil Pryde RS：X帆板帆翼模型正视图Figure 1.Front View of a Sail Wing Model

图2 Neil Pryde RS：X帆板帆翼模型俯视图Figure 2.Top View of a Sail Wing Model

表1 Neil Pryde RS：X帆板帆翼参数Table 1 Parameters of the Sail Wing for Neil Pryde RS：X Class

本文在对帆板帆翼进行单向流固耦合数值模拟时，根据实际需要假定帆翼为：帆面由均质材料组成，帆面材料属于各项同性材料，帆翼所发生的变形均在其屈服强度之内。帆翼的物理参数：帆面弹性模量为14.4 GPa、泊松比为0.33、密度为100.86 kg/m3，桅杆弹性模量为161.7 GPa、泊松比为0.23、密度为100.86 kg/m3（雷晓珊等，2019）。

1.2 计算域设置及网格划分

在对帆翼进行数值模拟时，为了保证计算域既可以满足计算精度又不浪费计算资源（马勇等，2013b），根据帆板帆翼尺寸，计算域和边界条件为：1）前端（x轴正方向）为12 m，边界条件为速度入口；2）后端（x轴负方向）为20 m，边界条件为自由出流；3）侧边界（y轴方向左右）为12 m，边界条件为对称平面；4）上边界（z轴正方向）为12 m，边界条件为对称平面；5）下边界（z轴负方向）为1 m，边界条件为对称平面；6）帆翼表面边界条件为wall壁面，如图3所示。在数值模拟时为了更加真实地模拟帆翼的运动，将桅杆、桅杆底以及帆翼底部设置为固定约束（图4）。将流场计算所得帆翼表面压力加载至结构场的帆翼模型，此处流固耦合面的选择应与流场中保持一致，施加载荷后帆翼表面如图5所示。

图3 计算域及边界条件Figure 3.Computational Domain and Boundary Conditions

图4 帆翼约束Figure 4.Constraint of the Sail Wing

图5 帆翼载荷Figure 5.Load of the Sail Wing

流体域及结构域进行计算前需要对其进行离散，即把原有的计算区域划分成若干个子计算区域，并且确定区域中的各个节点，进而生成网格。本文在数值模拟过程中需要计算多种攻角，且当风吹过帆翼时，帆翼会产生运动或者发生变形。综合考虑计算工况及计算过程中模型的变化，本文对计算域及帆翼模型进行非结构化网格划分。为了更加准确地描述帆翼周围流场，对帆翼周围进行适当加密，但由于桅杆的存在导致空气绕流不顺畅，使得整个帆翼背风面产生漩涡（马勇，2009），因此，在对帆翼进行网格划分时，对桅杆处的网格进行加密。网格划分完成后帆翼模型如图6所示，计算域网格情况如图7所示。

图6 帆翼网格分布Figure 6.Mesh of the Sail Wing

图7 计算域网格Figure 7.Mesh of Computational Domain

1.3 控制方程

帆翼周围空气在流动过程中遵循质量守恒、动量守恒以及能量守恒定律，但因其周围风速远小于声速的1/3，认为帆板帆翼周围的空气为不可压缩，且帆翼周围空气流动过程中热交换量较小，所以计算时未考虑能量守恒方程。为了避免计算量大等问题，流体域基于雷诺平均法对质量守恒方程和动量守恒方程中各量进行瞬时和脉动叠加处理，其形式如下（林虹兆，2016）：

其中，ui=(u，v，w)是速度在xi=(x，y，z)各方向上的分量，和分别代表时均速度和脉动速度；ρ为流体密度，ν为流体的运动粘度系数，为湍流产生的雷诺应力张量。

增加的雷诺应力项导致雷诺方程组本身不封闭，为了使方程组封闭，本文引入Realizablek-ε模型建立由湍流脉动引起的附加应力与时均应变率之间的联系，从而使雷诺方程封闭实现求解。Realizablek-ε的湍流动能方程（k方程）和耗散方程（ε方程）的最后形式（Shih et al.，1995）：

其中，k为湍流动能，ε为湍流动能耗散率，μt为湍流粘度，μ为湍流粘度系数，Gk为由平均速度梯度而产生的湍流动能项，σk=1.0，σε=1.2，C2=1.9，

1.4 升阻力系数分析

升力系数和阻力系数公式如下：

其中，L为升力，D为阻力，ρ为空气密度，U为风速，S为帆翼面积。

1.5 数值方法验证

为了验证数值模拟方法的准确性，在西北工业大学低湍流风洞对缩尺比为1∶15的刚性帆翼进行试验，帆翼上的力学数据通过ATI nano17系列六分量天平测得。采用Fluent软件对试验模型进行数值模拟，计算结果如图8所示。数值模拟结果与试验结果对比发现：在风速为5m/s和6m/s、攻角为10°～20°时，试验得到的帆翼上的升阻力与数值模拟计算结果显示出良好的一致性，验证了本文采用的数值方法的可靠性。

图8 数值模拟与试验结果对比Figure 8.Comparison of the Numerical Simulation and Experimental Results

2 结果

2.1 升阻力系数

升力系数随攻角的变化分为3个阶段，当攻角处于20°～25°范围内时，升力系数随着攻角的增加逐渐增大，25°攻角时升力系数达到最大值；当攻角处于25°～90°范围时，帆翼的升力系数随着攻角的增大而减小，所以对于女子Neil Pryde RS：X级别帆板帆翼而言，其失速角位于25°附近；当攻角大于90°时，升力系数变为负值，导致升力系数变为反向，此时升力变成使帆板发生横移的力，随着攻角的继续增加，升力系数逐渐增大，其阻碍帆板前进的作用显著增加（图9）。

图9 不同攻角下帆翼升阻力系数Figure 9.Lift and Drag Coefficient at Different Attack Angles

2.2 迎风面与背风面压力

所有航段帆翼背风面均为负压；迎风及横风航段攻角从20°增加到50°时，压力波动较为明显；攻角为20°时，帆翼背风面帆尾角处以及桅杆上部压力最大，负压最低点位于桅杆中下部，形成狭长区域；攻角为25°时，帆尾角处压力减小，但桅杆上部高压区域增大，且后帆边上部压力也开始增大，桅杆中下部最低负压区域逐渐减小；攻角为30°～50°时，帆尾角处、桅杆上部以及后帆边上部压力减小，但桅杆上部与后帆边上部所包围的正压区域面积逐渐扩大（图10）。

图10 不同攻角下帆翼（左）背风面、（右）迎风面压力云图Figure 10.Pressure Distribution of the Sail Wing（Left）Leeward Side and（Right）Windward Side at Different Attack Angles

2.3 流线图与速度云图

从图11可以发现，攻角为20°～25°时，流线沿帆翼弦向顺着帆翼表面流向后缘，流动顺畅且未发生流动分离；经过失速角之后，帆翼背风面流动分离明显，出现多个回流，并且随着攻角的增大流线分析区域加大。攻角为40°和100°时，帆翼背风面不仅产生回流，且伴随明显的漩涡，流动更加复杂。

图11 不同攻角下帆翼XY（Z=1.45m）截面上的速度云图及流线图Figure 11.Velocity Distribution and Streamlines on XY（Z=1.45 m）Section of the Sail Wing at Different Attack Angles

2.4 帆翼结构变形

从图12中可以看出，不同攻角下帆翼发生形变的位置基本一致，但最大变形值略有差异。综合来看，最大变形均发生在帆顶角处，这是由于帆顶角处与流体相对运动速度较高，并且到桅杆以及帆翼底面距离最大，所以最容易发生变形，且变形最大（图13）。

图12 不同攻角下帆板帆翼变形云图Figure 12.Deformation Distribution of the Sail Wing at Different Attack Angles

图13 帆翼最大变形曲线图Figure 13.Maximum Deformation of the Sail Wing

2.5 等效应力

等效应力表示各个方向上的材料应力差值，对帆翼表面进行等效应力变化分析，可以快速确定模型中受力最大的区域。根据单向流固耦合计算结果，得到不同攻角下柔性帆翼等效应力云图（图14）。从图14中可以看出，帆翼表面的等效应力分布不均匀，存在着应力集中的现象，桅杆顶端与帆上角连接处等效应力的数值最大，应力分布最为集中。此外，帆翼后帆边也存在较大的等效应力，而桅杆中下部以及帆翼底部等效应力较小。可见，在本文所涉及的攻角范围内，帆板帆翼的主要受风区域位于桅杆中上部和后帆边的中下部，以及两者之间的区域。

图14 不同攻角下帆板帆翼等效应力云图Figure 14.Equivalent Stress Distribution of the Sail Wing at Different Attack Angles

图15 帆翼最大等效应力曲线图Figure 15.Maximum Equivalent Stress of the Sail Wing

3 讨论

3.1 帆翼空气动力特性分析

Neil Pryde RS：X帆板比赛有外绕比赛航线（Ma et al.，2016）（图16），及与外绕航线对应的内绕航线（马勇等，2013a，2013c），其不同点在于内绕航线在航标1和航标4之间多进行一次绕标。但无论是内绕航线还是外绕航线，在运动过程中都要完成迎风航段、横风航段、顺风航段以及绕标，航线中涉及到的核心操作就是帆翼攻角调整（马勇等，2013a，2013c）。其中迎风航段、横风航段以及绕从1标、绕3标向2标航行、绕4标向1标航行时，帆翼攻角都在20°～50°之间；顺风航段帆翼攻角在80°～100°之间；而外绕航线绕2标向3标航行时，帆翼攻角则从20°～50°之间转向80°～100°之间。在不同航段以及绕航标1、2、3、4时，帆船运动员需要根据外界风况（风速、风向角）变化以及对手所在位置调整帆翼攻角，从而使帆板按照自己的意图航线。

图16 奥运会比赛外绕航线图（Ma et al.，2016）Figure 16.Olympic Outer Circuit in the Sailing Regatta

在迎风、横风航行阶段航行时，帆翼迎流攻角调整范围一般为20°～50°，结合图9（a）的计算结果，攻角增加使帆翼流动分离导致帆翼压差阻力迅速增加；而在顺风航段中，从图9（b）结果可知，阻力系数变化曲线较为平缓，升力最大时攻角为25°。从迎风、横风航行阶段属于典型的升力型航段分析，为获取更快的航行速度，需要使帆翼上的升力尽可能大而阻力尽量小，因此，攻角选择为20°～25°。本文考虑单向流固耦合影响并进行了攻角每隔5°的计算，确定最佳攻角为25°左右，这与已有研究最佳攻角在20°～30°之间（何海峰，2012；罗晓川，2012）范围符合，并进一步精确。所以，帆板运动员在迎风航段、横风航段以及绕从1标、绕3标向2标航行、绕4标向1标航行时，将帆翼攻角调整到25°左右可以使帆板航速更快。

随着攻角的增大，帆翼背风面负压最低点沿桅杆逐渐下移，最低负压区域也逐渐减小，帆翼迎风面与来流夹角较小。由于桅杆与来流垂直，所以其受压最大，而后帆边受压最小是桅杆和帆翼相互影响造成的（马勇，2009）。此外，帆翼迎风面桅杆处速度最大，均大于来流速度（6 m/s），后帆边处速度最小，背风面后帆边尾流处速度最大，帆翼拱度最大处速度最小，这与其他级别帆板帆翼的空气动力性能结果类似（马勇等，2016；Ma et al.，2016）。由于帆翼拱度特点，顺风航段帆翼背风面压力变化较小，正压区域位于帆翼中部，随着攻角的增大，正压中心逐渐向桅杆方向移动；帆翼迎风面与来流几乎垂直，导致帆翼整个迎风面均受正压。

帆船帆板运动员在比赛中要随时注意风况（风速或者风向角）变化，变化后的风速和风向角度可能使得帆翼产生明显的涡脱现象，从而引起帆翼的变形和振动或者失速，此时，帆船帆板运动员要及时调整帆翼参数，使帆翼处于较优攻角状态（雷晓珊等，2019）。

3.2 帆翼结构变形分析

有关奥运会Neil Pryde RS：X级别帆板帆翼与其周围流场的流固耦合作用的结论与其他帆翼流固耦合作用类似（Augier，2012；Bak et al.，2013；Durand et al.，2014）。从结构域来看，帆翼发生形变的位置基本一致，最大变形基本均发生在帆顶角处，且迎风及横风航段、顺风航段帆翼最大变形与帆翼的升力、阻力变化有关（图9、图13）。在迎风航段，当攻角小于失速角（25°附近）时，随着攻角的增加升力迅速增大，而此时阻力较小，因此，在攻角为20°～25°时，最大变形量随攻角逐渐增大，且增大幅度较为明显；当攻角达到失速角之后，升力突然减小，阻力继续增大，此时最大变形迅速减小；当攻角继续增大，升力缓慢减小，阻力小幅度增大，此时作用于帆翼面上的升力仍大于阻力，所以帆翼最大变形仍继续增加，但幅度变小。在顺风航段，当攻角＜90°时，作用于帆翼上的升力逐渐减小，阻力基本无变化，但此时阻力大于升力，所以帆翼的最大变形量逐渐减小；当攻角＞90°时，帆翼拱度发生转向，阻力成为构成帆板前进的推进力的一部分，升力成为阻碍帆板前进的力，此时尽管阻力远大于升力，但是升力随攻角不断增加，所以最大变形与之前相比有明显增大，但其趋势为逐渐减小。

此外，帆翼表面应力分布并不均匀，后帆边处存在较大的应力集中，桅杆顶端与帆上角连接处应力分布最为集中，但桅杆中下部以及帆翼底部等效应力较小（图14）。由计算结果对比可知，与顺风航段相比，迎风及横风航段各个攻角下帆翼所受最大等效应力值偏小（图15）。迎风、横风航段，20°～30°攻角范围内，最大等效应力值相差较小；攻角为45°时等效应力最大，约为4.5×105 Pa。顺风航段下，85°攻角时最大等效应力最小，约为5.2×105 Pa，90°时最大，约为6.3×105 Pa；95°和100°攻角时，最大等效应力几乎相等。帆船帆板运动员在比赛中要根据风况（风速和风向角）变化，随时观察帆翼应力集中位置与大小，及时调整帆翼攻角，避免应力集中与失速导致涡激振动等叠加作用。