NACA 0012摆动水翼水动力特性的二维数值模拟

吕元博 ，田新亮 ，李欣 ，宋春辉

1上海交通大学海洋工程国家重点实验室，上海200240

2高新船舶与深海开发装备协同创新中心，上海200240

0 引言

自Hine等［1］率先研制出波浪滑翔机以来，作为一种直接将波浪能转化为前进动能的海洋装置，波浪滑翔机被广泛应用于海上长时间观测作业。由于波浪滑翔机具有续航能力强、噪声小等优点［2-3］，其在检测海洋环境要素、海洋灾害预报、海洋科学研究等方面呈现出广阔的应用前景。

波浪滑翔机的水翼是影响其航行性能的重要因素。针对水翼的水动力性能，国内外学者进行了大量研究。Kraus［4］通过建立模型，对波浪滑翔机各部位进行了六自由度模拟，确定水翼最大摆角为20°时最优，模拟得出滑翔机在不同海况下的各项水动力学参数。贾立娟［5］利用FLUENT软件研究了翼型、摆角等对水翼水动力学特性的影响，发现摆角为18°时水翼的水动力性能较好。Sun等［6］针对果蝇翅膀的主动摆动进行了数值模拟，发现动态平均升力系数可以达到准静态升力系数的2倍。Andersen等［7］对某一对称翼型进行了主动摆动和主动垂荡运动的数值与实验研究，得出2种运动模型尾流场相似，但对于垂荡运动在低频高幅值时，在1个振荡周期内将产生2对对称的尾窝。胡合文［8］通过势流与粘流结合的方法对滑翔机整体及水翼部分进行了计算，得到水翼的最佳摆角在 15°附近。张晓庆等［9］和刘焕兴等［10］以鱼类尾鳍摆动推进为背景，利用计算流体动力学（Computational Fluid Dynamics，CFD）软件，模拟了二维刚性水翼摆动的非定常水动力性能。

本文选用NACA 0012作为波浪滑翔机的翼型。针对该翼型，Ohtake等［11］和 Yonemoto等［12］均进行过大量研究，获得了定常流场下机翼的性能参数。Read等［13］和 Triantafyllou等［14］通过一套水翼主动摆动装置进行实验，分析了NACA 0012型水翼在主动摆动和升沉运动下的水动力性能，得出有效攻角的不平滑变化会使水翼性能降低。Schouveiler等［15］在麻省理工学院（MIT）拖曳水池通过实验给定NACA 0012翼型正弦规律的升沉运动和主动摆动，得出最佳推进效率的斯特哈尔数St值为0.25～0.35，St较大时水动力性能会降低。由于波浪滑翔机是通过垂荡方向的运动引起机翼的被动旋转，因此，通过水翼的主动运动模拟机翼的实际运动可能存在较大误差。被动旋转则需要考虑水翼的摆动运动与流场运动的瞬态耦合，且要考虑水翼在流场中运动时所受力矩以及限位角的约束条件，存在一定的仿真难度，但对水翼摆动过程有更好的模拟，而目前针对非定常流场中的水翼被动摆动的研究少有涉及。

本文将基于STAR-CCM+软件中的DFBI模块，以NACA 0012翼型为研究对象，对波浪滑翔机水翼做垂荡运动时引起的绕自身中心轴的被动旋转运动进行模拟，对水翼主动旋转与被动旋转2种情况下的推力系数进行对比，并分析限位角、波浪参数等因素对水翼性能的影响，以期为波浪滑翔机的设计与研究提供参考。

1 理论概述

1.1 波浪滑翔机运动原理

波浪滑翔机的运动原理如图 1（a）［16］所示。波浪滑翔机由水面母船和水下滑翔机2部分组成。当水面母船随波浪上升时，缆绳会拉紧，带动水下滑翔机向上运动，由于水的作用力，会使水下滑翔机的水翼逆时针翻转。当水面母船随波浪下降时，缆绳松弛，水下滑翔机自由下降，水翼受到水的反作用力而顺时针摆动。水翼在不断的上、下摆动过程中会产生向前的推力，拖拽水面母船不断前进。

图1（b）所示为波浪滑翔机水翼在波浪起伏过程中的受力分析结果。水翼在随波浪上升的过程中，受垂直于翼面向下的水动力作用，会产生水平方向的分力Fx和垂直向下的分力Fy。水翼在随波浪下降的过程中，受垂直于翼面向上的水动力作用，同样会分解为水平方向的分力Fx和垂直向上的分力Fy。两个过程中的水平分力Fx成为使波浪滑翔机向前运动的推进力。在此，将整个研究对象简化为图1（b）中的水翼模型。

1.2 水翼被动运动方程

波浪滑翔机水翼的运动分为2部分：一部分是在波浪驱动下的y方向的平移运动；另一部分是在平移过程中受水动力作用的摆角运动。

水翼的垂向运动h(t)为

式中：t为运动时间；f为运动频率；h0为垂向运动幅值。

水翼摆角的被动运动控制方程为

式中：J为绕重心的转动惯量；̈为水翼被动旋转的角加速度；M为作用在水翼上的水动力力矩，可以根据水翼表面压强积分得到。对式（2）进行积分，可以得到角速度̇和摆角θ。

在设计波浪滑翔机的水翼时，当水翼摆角θ达到限位角θ0时，摆角达到最大值，直到力矩反向，水翼开始反向摆动。即

1.3 水翼主动运动方程

对于水翼摆动问题，如果在垂直方向和摆动方向均采用主动运动模型，则其摆动方向的运动规律为

式中，ψ为平移运动和摆动运动的相位差。

对于主动旋转工况，有效攻角α(t)为水平来流速度和平移速度的合速度与水翼的摆角之间形成的实际攻角，其表达式为

式中：h'(t)为平移速度；U为来流速度，即水翼前进速度。

在此，引入流体力学中的相似准则斯特哈尔数St：

式中，A为尾涡区域在平移方向的宽度，近似表示为平移幅度的2倍，即A=2h0。代入式（6），得St=2fh0/U。与式（5）联立，可得

式中，αmax为最大有效攻角。

1.4 推力系数

水翼的瞬时推力系数为

式中：Fx(t)为水翼的瞬时推力；ρ为水的密度；b为水翼弦长；s为水翼展长。

式中，T为水翼旋转周期。

平均推力系数为

2 数值模型及验证

2.1 模型参数

水翼采用NACA 0012翼型剖面，最大厚度位于弦长1/3处，如图2所示。本文取水翼特征弦长b=0.1 m，摆动轴位于最大厚度处。

2.2 数值模型

如图3所示，二维摆动水翼的计算区域为50b×30b。翼表面为圆形计算域，其直径为5b，用于翼的旋转运动。对于边界条件的设定，左边界和上、下边界均为速度入口，右边界为压力出口。

计算采用SST k-ω湍流模型，利用STAR-CCM+软件的切割体网格对流域进行分割。如图4所示，在水翼表面，采用了网格无相对变形且数值交换较好的重叠网格，并在计算域及尾流区域进行局部加密。边界层取为30层，控制y+严格小于1。在水翼运动方面，采用STAR-CCM+软件中的流体与固体相互作用（Dynamic Fluid Body Interaction，DFBI）模块，通过输入上、下表面边界的正弦边界条件，模拟波浪起伏过程，然后通过DFBI模块计算出水翼在升沉运动下的被动旋转运动情况，并通过单自由度旋转模块控制水翼可旋转的最大角度，即限位角θ0。

2.3 数值模型验证

为验证所建立的数值模型是否满足计算要求，对网格尺寸、时间步长、计算结果可靠性进行验证。

2.3.1 网格尺寸验证

考虑到计算的准确性及高效性，计算分析了4种网格尺寸固定角度下的升力系数CL，并进行比较。表1为在固定角度5°时4种网格尺寸A1～A4对应的CL值及相对误差（指计算值与实验值的相对误差），从中可观察到从A2开始计算值趋于稳定。根据对翼型的经验分析，CL实验值一般略小于CFD计算值，且可观察到计算值与理论值十分接近。

表1 网格尺寸验证Table 1 The mesh size validation

图5为4种网格尺寸的水翼上、下表面压力系数Cp分布曲线。从中可观察到：A2，A3和A4均重合较好；A1在上表面约1/3处出现了波动，结合该处的流场进行分析，发现此处的边界层出现了明显分离；A1边界层网格尺寸较大，无法精确捕捉。综合考虑计算量及计算准确性，选取A2作为网格模型。

2.3.2 时间步长验证

在A2网格模型下，通过STAR-CCM+软件中的Motion模块，对水翼施加主动升沉运动叠加绕旋转中心的旋转运动，时间步长Δt分别取T0/100，T0/200，T0/300，和T0/400共4种工况，计算升力系数峰值CLmax（表2）。由表2可知，除T0/100时间步长以外，其他工况的CLmax均较接近。

表2 网格尺寸A2时间步长验证Table 2 Time step validation of A2

图6为不同时间步长下，升力系数CL与时间T的关系曲线。与其他时间步长相比，T0/100在波峰波动处出现了较大的相位偏差，且波峰数量减少。综合考虑计算量及计算准确性，最终选用Δt/T0=1/200为本计算模型的时间步长。

2.3.3 可靠性验证

为验证模型的准确性，在固定摆角、相同雷诺数下，将计算值与文献［12］中的2个实验值进行对比，攻角α分别在-5°～5°之间每隔1°选取，Re=105。图7为固定摆角下水翼的升力及阻力系数。由图7（a）可见计算值与实验值吻合较好，且CL比Yonemoto的计算值更精确、数据点范围更广，图中的直线为斜率a=2π的理论值。而从图7（b）可知，阻力系数CD和Yonemoto计算值均出现了略低于实验值的情况，且计算结果基本一致。因此，本文数值模型可准确用于计算水翼的水动力特性问题。

3 计算结果与分析

针对水翼主动与被动旋转2种工况，在STAR-CCM+软件中分别对水动力性能进行仿真，并对比仿真结果，着重讨论限位角、波高、波浪频率等对水翼推力系数的影响。主动与被动旋转工况分别如表3和表4所示。

表3 水翼主动摆动计算工况Table 3 The calculation conditions for active oscillating wing

表4 水翼被动摆动计算工况Table 4 The calculation conditions for passive oscillating wing

3.1 水翼主动旋转与被动旋转对推进性能的影响

为了比较水翼主动旋转与被动旋转对推进性能的影响，选取进流速度U=0.4 m/s，最大有效攻角αmax=15°，h0=0.8 m，对不同St数下的主动与被动工况进行分析，结果如图8所示。

从图8可以发现，不论是主动摆动运动还是被动摆动运动，在一个周期内水翼的平均推力系数均随St数的增大而增大，且主动摆动的平均推力系数大于被动摆动，采取主动摆动方法时其平均推力系数将比被动摆动高出约30%。分析原因，主要在于2种方式水翼摆动的旋转速率不同。采取被动摆动方式模拟时，水翼在水中受到水反作用力的力矩，在升沉时瞬间完成翻转，旋转速率在瞬间达到最大，达到最大限位角后变为0，之后以此限位角做平移运动。而传统的主动摆动方法给定水翼的正弦摆动速率，水翼在靠近平衡位置处旋转速率最快，从平衡位置到波峰的过程中水翼并非瞬间完成翻转，而是以正弦速率周期性摆动，所以此正弦速率的摆动过程会对水产生向后的额外推力，从而使水翼获取额外的向前推进力，使平均推力系数增大。从能量的角度分析，被动摆动水翼的输入能量仅仅为垂荡运动所需能量；而主动摆动的输入能量则为主动垂荡运动与主动旋转所需的能量之和，相较于被动摆动，增加了主动旋转部分的能量，故也将产生更大的推力。

由图8可以发现，在St数较小时，两者的平均推力系数差值较小，随着St数的逐渐增大，两者的误差逐渐增加并趋于稳定，二者平均推力系数的差别可达30%。可见，主动旋转与被动旋转对推进性能的影响较大，对于波浪滑翔机水翼的水动力性能研究，采用传统的主动摆动运动模型会增加计算误差，不能准确模拟波浪滑翔机水翼的实际运动状态。

3.2 水翼摆动限位角θ对推进性能的影响

为了研究水翼摆动限位角对推进性能的影响，选取5种限位角进行分析，其他参数设置如下：进流速度U=0.4 m/s，波浪频率f=0.2 Hz，波高h=0.8 m，推力系数时历如图9所示。

从图9（a）可以分析得到，当水翼被动摆动的最大限位角θ0在 10°，15°和 20°时，推力系数随时间做周期性的类正弦变化，且随着最大限位角增加 C（t）也逐渐增加，在限位角θ0=20°时，推力系数峰值达到最大。

图9（b）为水翼被动摆动的最大限位角θ0在20°，25°和30°工况下对应的C（t）-T曲线。由图可知，在前半个周期，推力系数的峰值仍随限位角的增加而增大，但在后半周期，在限位角θ0=25°和θ0=30°时 C（t）峰值出现了突变和尖点，没有呈现很好的周期性。分析突变产生的原因，当摆角大于25°时，水翼开始逐渐发生失速，会受远处来流与水翼上表面流动分离的相互影响，对推力性能产生较大影响。因此当摆角超过某一临界值时，水翼推力系数将不再继续升高，会出现明显的波动。

图10为在不同限位角下水翼在T/4时刻被动摆动运动时的涡量图。由涡量图可知，在限位角θ0=10°～20°时，水翼上表面分离区较小，呈现出流线型的绕流。但随着摆角逐渐增大，从25°开始，水翼的尾缘产生明显的涡旋，翼型上表面的分离区面积逐渐增加，后缘附近开始有翼尖涡，一定程度上降低了水翼的水动力性能。如图10所示，水翼的脱落涡均呈逗号形式，尾部产生的射流增加了水翼的推力，随着限位角的增大，脱落涡逐渐减小。

综合以上分析，选取θ0=20°作为水翼被动摆动的最大限位角。

3.3 波高h和波浪频率f对水翼推进性能的影响

确定限位角θ=20°后，在来流速度U=0.4 m/s，波浪频率f=0.2 Hz的条件下，选取波高h=0.1，0.4和0.8 m共3种工况，得到了对应的推力系数C（t）关于周期T的曲线（图11）。

由图11可见，C（t）随波高h的增大而增大，且在波高h=0.1和0.4 m时，推力系数曲线未呈现出明显的周期性，C（t）的大部分小于0，说明流场并没有趋于稳定。当h=0.8 m时，曲线开始出现明显的周期性变化。此外，在波浪频率f=0.2 Hz时，随着波高h的增加，推力系数逐渐增大。且在一定来流速度和波浪频率下，波高需达到一定的高度才可产生稳定向前的推力，流场和推力系数也可趋于稳定。

上述对于波高h的讨论限定在单一频率，即f=0.2 Hz条件下，规律不具有普遍性。为找到波高h对水翼推进性能的影响规律，分别选取多组工况进行了计算对比。

由于波高h=0.8 m时推力系数周期化变化，故取波高h在0.8 m附近变化，波浪频率分别取f=0.1，0.2和0.4 Hz共3种情况，研究不同波高h对C（t）的影响规律（图12）。由图12可见，在不同波浪频率下，推力系数C（t）均随波高h的升高而增大。

由上述计算结果可知，在其他条件一定时，从波高h=0.8 m开始，推力系数C（t）趋于稳定。故在波高h=0.8 m，来流速度U=0.4 m/s时，选取不同波浪频率进行计算。

图13是h=0.8 m时不同波浪频率下被动摆动水翼的C（t）-T曲线。由图可知，当波浪频率f=0.2 Hz时，推力系数在C（t）=0附近波动，波浪频率f每增加一倍，C（t）增幅远高于一倍。当f=0.8 Hz时，C（t）峰值达到20。可见，波浪频率f对推力系数C（t）的影响更为明显。在波高h=0.8 m时，推力系数C（t）随波浪频率f的增大而增加。

同样，上述对于波浪频率的讨论也限定在单一波高，即h=0.8 m条件下，不具有普遍性，为找到波浪频率f对水翼推进性能的影响规律，分别取多组工况进行了计算对比。

图14所示为不同波浪频率f对C（t）的影响规律。由图可知，在不同波高下，推力系数C（t）随f增加而增大，当波浪频率从f=0.2 Hz扩大一倍至f=0.4 Hz时，C（t）曲线的峰值增加了 5～10倍。在波浪频率较低时，即f=0.1 Hz和f=0.2 Hz时，推力系数C（t）峰值相差不大。由此可知，在波浪频率较小时，波浪频率f的改变对推力系数的影响较小。

4 结 论

本文以波浪滑翔机的水下滑翔部分水翼为研究对象，对NACA 0012型水翼在波浪中的被动摆动情况进行了水动力学分析及计算，得到如下主要结论：

1）通过自主建立的STAR-CCM+软件的重叠网格与DFBI模块相结合的模型，实现了对水翼垂直方向主动平移状态下的被动旋转的模拟。

2）相同条件下，主动摆动和被动摆动运动对水翼的推进性能有较大影响，二者的平均推力系数可相差30%。

3）通过分析计算结果与涡量图，得到滑翔机水翼的最大限位角在20°附近时滑翔机的水动力性能达到最优。

4）在其他条件一定时，滑翔机水翼的推进性能受波浪参数影响较大，且随波高h和波浪频率f的增加而增加，在波浪频率较小时，波浪频率f的改变对推力系数影响不大。

参考文献：

［1］HINE R，WILLCOX S，HINE G，et al.The wave glid⁃er：a wave-powered autonomous marine vehicle［C］//OCEANS 2009，MTS/IEEE Biloxi-Marine Technology for Our Future：Global and Local Challenges.Missis⁃sippi：IEEE，2009：1-6.

［2］徐春莺，陈家旺，郑炳焕.波浪驱动的水面波力滑翔机研究现状及应用［J］.海洋技术学报，2014，33（2）：111-117.XU C Y，CHEN J W，ZHENG B H.Research status and applications of wave gliders［J］.Journal of Ocean Technology，2014，33（2）：111-117（in Chinese）.

［3］杨海，刘雁集，张凯.实验尺度下无人水下滑翔机设计与试验［J］.中国舰船研究，2016，11（1）：102-107，120.YANG H，LIU Y J，ZHANG K.Design and experi⁃ment for laboratory-scale autonomous underwater glid⁃ers［J］.Chinese Journal of Ship Research，2016，11（1）：102-107，120（in Chinese）.

［4］KRAUS N D.Wave glider dynamic modeling parame⁃ter identification and simulation［D］.Hawaii：Universi⁃ty of Hawaii，2012.

［5］贾立娟.波浪动力滑翔机双体结构工作机理与动力学行为研究［D］.天津：国家海洋技术中心，2014.JIA L J.Study of operation principal of two-part archi⁃tecture and dynamic behavior of wave glider［D］.Tian⁃jin：National Ocean Technology Center，2014（in Chi⁃nese）.

［6］SUN M，TANG J.Unsteady aerodynamic force genera⁃tion by a model fruit fly wing in flapping motion［J］.Journal of Experimental Biology，2002，205（1）：55-70.

［7］ANDERSEN A，BOHR T，SCHNIPPER T，et al.Wake structure and thrust generation of a flapping foil in two-dimensional flow［J］.Journal of Fluid Mechan⁃ics，2017，812：R4-1-R4-12.

［8］胡合文.波浪滑翔机的水动力分析［D］.哈尔滨：哈尔滨工程大学，2015.HU H W.Hydrodynamic analysis of wave glider［D］.Harbin：Harbin Engineering University，2015（in Chi⁃nese）.

［9］张晓庆，王志东，张振山.二维摆动水翼仿生推进水动力性能研究［J］.水动力学研究与进展（A辑），2006，21（5）：632-639.ZHANG X Q，WANG Z D，ZHANG Z S.Hydrodynam⁃ic study of bionic propulsion for 2D flapping foil［J］.Journal of Hydrodynamics（Ser.A），2006，21（5）：632-639（in Chinese）.

［10］刘焕兴，苏玉民，庞永杰.非正弦摆动对水翼水动力性能的影响［J］.华中科技大学学报（自然科学版），2016（2）：15-20.LIU H X，SU Y M，PANG Y J.Effects of the nonsinu⁃soidal oscillating on the hydrodynamic performance of the hydrofoil［J］.Journal of Huazhong University of Science and Technology（Natural Science Edition），2016（2）：15-20（in Chinese）.

［11］OHTAKE T，NAKAE Y，MOTOHASHI T.Nonlinear⁃ity of the aerodynamic characteristics of NACA0012 aerofoil at low Reynolds numbers［J］.Journal of the Japan Society for Aeronautical and Space Sciences，2007，55（644）：439-445.

［12］YONEMOTO K，TAKATO K，OCHI H，et al.Kutta condition violation in two-dimensional NACA0012 airfoil at low Reynolds number［C］//26th AIAA Ap⁃plied AerodynamicsConference.Hawaii： AIAA，2008：6399.

［13］READ D A，HOVER F S，TRIANTAFYLLOU M S.Forces on oscillating foils for propulsion and maneu⁃vering［J］.Journal of Fluids and Structures，2003，17（1）：163-183.

［14］TRIANTAFYLLOU G S，TRIANTAFYLLOU M S，GROSENBAUGH M A.Optimal thrust development in oscillating foils with application to fish propulsion［J］.Journal of Fluids and Structures，1993，7（2）：205-224.

［15］SCHOUVEILER L，HOVER F S，TRIANTAFYL⁃LOU M S.Performance of flapping foil propulsion［J］.Journal of Fluids and Structures，2005，20（7）：949-959.

［16］Liquid Robotics.Converting wave motion into propul⁃sion［EB/OL］.（2014-09）［2017-06］.https：//www.liq⁃uid-robotics.com/platform/how-it-works/.