1 引言
近年来，随着船舶行业的蓬勃发展，舰艇的航速大大提升，舰艇水下结构的水弹性振动得到了造船界的重视。若舰艇的航速大，舵翼的结构强度不高，舵翼在水下运动时有可能会产生颤振，造成结构破坏，威胁航行安全。而在舰艇航速较低时，会产生一种容易被忽视的振动现象“涡振”，涡振即结构在流场中产生涡的频率与结构固有频率一致时发生的共振现象。涡振虽然是限幅运动，但会造成结构疲劳，长此以往，造成结构破坏。对于潜艇来说，发生涡振时会产生显著的噪声，破坏隐身性能。
针对水翼的振动现象，国内外的众多学者展开了研究。Chen等[1]采用SST k-ω湍流模型对平行槽水翼叶片进行数值分析，研究了水翼位置和间距对涡产生的抑制效果和流动特性的影响。结果表明，叶片间距越小，越容易产生涡，并根据分析结果选取翼型。Huang等[2]选用NACA0012翼型进行了试验，着重探究了翼梢对涡激振动的影响，发现翼梢处容易产生涡量大的涡。Huerre等[3]以及Oertel[4]通过试验分别给出了流动不稳定性时涡发放特点和钝体尾流形态的概述。Chae[5-6]通过对二自由度下二维水翼的研究，讨论了水动力阻尼、非线性流固耦合响应等对水翼振动特性的影响。余志兴等[7-8]采用N-S方程以弹簧系统代替实际固体变形产生恢复力，发现要构成翼型涡振，攻角要足够大。侯磊等[9]采用离散涡方法，分析高雷诺数下来流速度及攻角对水翼涡发放频率影响规律。刘胡涛等[10]对二维水翼结构进行流固耦合运动分析，探讨了初始攻角、刚心位置以及来流速度对水翼振动的影响。郭春雨等[11] 采用DEM模型对三维水翼梢涡处流场和气核运动特性的问题进行数值分析，为梢涡空化初生的机理研究建立基础。贾文隽等[12]分别采用传统CFD方法和格子玻尔兹曼流固耦合方法对平直机翼和仿生机翼进行数值模拟，并通过气结构响应对比找出了仿生机翼的优点及不足。赵宇等[13]采用数值方法对绕水翼的非定常空化流动进行了流动计算分析，结果表明空化现象对水气两相转换起到了促进作用，空化现象影响涡结构空间分布。陈衍茂等[14-15]采用同伦分析法、增量谐波法等研究含有立方非线性的二自由度机翼的极限环颤振，并对极限环的幅值和频率进行求解。
综上所述，对于水翼涡振特性的研究主要集中在单一翼型涡振机理等方面的研究，对于相似翼型的涡振特性差异方面存在空白。本文首先将所研究的翼型设置为刚体，探究其泄涡特性，然后通过双向流固耦合的方法，研究其涡振特性，并对几种相似翼型进行对比分析，为舰艇舵翼选型提供参考。
2 控制方程及流固耦合求解方法
2.1 流体控制方程
由于舵翼为流线型，后缘产生的涡多为细涡，雷诺平均湍流模型难以捕捉，故采用大涡模拟。大涡模拟综合采用了直接模拟和雷诺平均模拟的思想与方法。将湍流中的脉动划分尺度，大尺度湍流直接采用数值求解。而对小于计算网格尺度的湍流模型建立模型求解，其控制方程为：
(1)
(2)
其中，τij为亚网格应力，
为流体密度， μ为动力粘性系数，
均为滤波后的场变量。
2.2 结构振动控制方程
在本研究中，将三维水翼简化为二维截面，并引入弹簧阻尼系统，将水翼的振动简化为竖向和绕旋转轴扭转的两自由度的运动。
(3)
(4)
式(3)和式(4)分别是结构竖向和绕旋转轴扭转的振动控制方程。
和y(t)分别为t时刻结构竖向振动加速度、速度和位移；
和α(t)分别为t时刻结构的扭转振动加速度、速度和位移；ξh是结构竖向运动的阻尼比，ξα是结构做扭转运动的阻尼比。ωh和ωα分别为结构的竖向振动圆频率和扭转振动圆频率；FL(t)、Mα(t)分别为t时刻结构受到的流体自激力和力矩；m、Im分别为结构的单位展长质量和结构的扭转惯性矩。
2.3 双向流固耦合求解方法
双向流固耦合即流场计算结果压力载荷通过流固耦合界面传递给结构场，同时结构场的位移变形结果通过流固耦合界面传递给流场，双向流固耦合分析流固耦合交界面处的数据传递是双向的。流场和结构场之间是相互影响的，考虑结构的变形对流场的影响。
本文在计算水翼涡激振动响应的每一个时间步长内，先求解流体控制方程，得到作用在水翼上的水动力FL(t)和力矩Mα(t)，将得到的水动力和力矩代入结构振动方程式。采用Newmark-β法，求解结构的振动响应，利用动网格技术，将得到的速度赋予结构及结构壁面附近的网格，使其随结构一起做刚体运动，以实现水翼与流场的双向耦合。
3 物理模型及网格
计算选取了对称翼型NACA0009、NACA0012和NACA0015，其中重点介绍了NACA0012翼型。NACA翼型末尾两位数字代表了相对厚度的大小，NACA0015翼型最厚。水翼弦长为0.1 m长，旋转轴设置在距前缘1/3弦长处，来流攻角为0°。
采用的网格策略为“刚性边界层运动区域+动网格变形区域+外流场”的划分策略，翼型表面附有贴体网格，计算Y+值，保证近壁面处网格的质量。翼型外圆形区域为随动网格，随翼型一起运动，以保证翼型计算的收敛性及精度。圆环区域为动网格区采用三角形网格，以保证具有较好的运动性。最外围为流场网格。图1为所使用网格和局部放大图。
边界条件设置：左侧弧形边界以及上下边界设置为速度入口(velocity-inlet)，右侧边界设置为压力出口(pressure-outlet)，翼型采用无滑移壁面边界(wall)。
图1 网格示意图
Fig.1 Mesh
湍流模型选用大涡模拟。计算算法采用“Coupled”， 该算法同时求解流体的连续方程、动量方程和能量方程，适合求解本文关于不可压缩及等温流动的耦合问题。采用二阶迎风格式，以提高计算精度。水翼由于其结构特点，泄涡频率较高，时间步的选取影响计算精度，本文经过时间步无关性验证，选取时间步长为0.000 5 s，限于篇幅不再赘述。
本文二维水翼的结构参数取值来自于工程实例。水翼弦长为0.1 m，展长为0.2 m，材料为POM(聚甲醛树脂)，采用一端固支的方法，基于ANSYS平台对其提取质量、模态等结构参数。竖弯刚度由式(5)计算得出:
(5)
式中: f为水翼一阶竖弯固有频率，m为水翼单位展长质量。
在对水翼进行模态求解时发现：水翼的扭转固有频率要远大于竖弯固有频率。竖弯方向的涡振在低速下更容易发生，本文只分析了竖弯方向单自由度的一阶涡振现象.本文使用的翼型结构参数取值如表1。
表1 水翼结构参数
Table 1 Structural parameters of hydrofoil
NACA翼型单位展长质量/(kg·m-1)扭转惯性矩/(kg·m2·m-1)竖弯固有频率/Hz竖弯刚度/(N·m-1)000923.164.8742.291.63×106001230.886.5355.153.71×106001538.608.2464.186.24×106
4 刚性水翼泄涡特性分析
水翼设置为刚体，分别计算了NACA0012水翼在流速0.5 m/s、1.0 m/s、1.5 m/s、2.0 m/s下的升力系数。图2、图3分别为水翼在流速0.5 m/s、1.0 m/s下的升力系数随时间变化曲线，可以看出，升力系数变化的范围不大，在计算稳定后，升力系数曲线呈规律的正弦波形，反映了水翼表面受到了周期性变化的垂向载荷。随着流速的增加，水翼受到的垂向载荷频率变大。
图2 0.5 m/s时升力系数曲线
Fig.2 Lift coefficient curve at 0.5 m/s
图3 1.0 m/s时升力系数曲线
Fig.3 Lift coefficient curve at 1.0 m/s
对4种流速下的水翼的升力曲线进行快速傅里叶变换(FFT)，得到水翼受到的垂向载荷的频率，即水翼的泄涡频率，表2为提取到的各个流速下的泄涡频率。图4表示流速与泄涡频率的关系，可直观看出泄涡频率与流速呈现线性关系。
表2 泄涡频率
Table 1 Vortex shedding frequency
流速/(m·s-1)0.51.01.52.0泄涡频率/Hz23.3350.3777.30103.22
图4 流速与泄涡频率关系曲线
Fig.4 Relationship between velocity and frequency of vortex shedding
以上为NACA0012翼型的分析，采用同样的方法，计算了NACA0009以及NACA0015翼型在同流速下的泄涡频率，如图5所示。
图5 各翼型泄涡频率曲线
Fig.5 Comparison of vortex shedding frequency of each airfoil
从图5可以看出，本文选取的3种对称翼型的泄涡频率都随着流速的增加而增加。在同流速的比较下，水翼的厚度对泄涡频率产生影响，厚度越大的水翼，泄涡频率越低。
5 弹性水翼涡振特性分析
前文对刚性水翼型的泄涡特性进行了分析，只考虑了流场对水翼的作用，得到的规律较为简单。本节将水翼设置为弹性水翼，使用双向流固耦合的方法，对水翼的涡振特性进行分析。
在计算刚性水翼NACA0012时，水翼在1 m/s时泄涡频率达到了50.37 Hz，与该翼型的自振频率55.15 Hz十分接近，现在1 m/s左右的区间进行细化，提取水翼在不同流速下的竖向位移图。
图6、图7和图8分别给出了水翼在0.5 m/s、1.0 m/s、1.5 m/s时的竖向振动时域图。在低流速时水翼振动较为规律，0.5 m/s时水翼振幅达到最大数值后振幅略有缩小，而后以稳定振幅振动。1 m/s时水翼振幅呈现“喇叭形”，推测在该流速下发生涡振。1.5 m/s时水翼振幅不规则，反映流场在较高流速下湍流度变大，水翼受力情况更复杂。
图6 0.5 m/s竖向位移曲线
Fig.6 0.5 m/s vertical displacement diagram
图7 1.0 m/s竖向位移曲线
Fig.7 1.0 m/s vertical displacement diagram
图8 1.5 m/s竖向位移曲线
Fig.8 1.5 m/s vertical displacement diagram
图9—图11为水翼的振动频谱曲线。在0.5 m/s的流速下，水翼的振动频率较低，与前文计算刚性水翼在该流速下的泄涡频率23.33 Hz相近。在1 m/s时，水翼振动频率与自振频率相近，符合发生涡振的特征。随着流速增加，频谱图出现双峰特征，左侧峰集中在水翼的自振频率附近，呈现明显宽谱特性，右侧峰为水翼的泄涡频率，较为集中。
图9 0.5 m/s振动频谱曲线
Fig.9 0.5 m/s vibration spectrum
图10 1.0 m/s振动频谱曲线
Fig.10 1.0 m/s vibration spectrum
图11 1.5 m/s振动频谱曲线
Fig.11 1.5 m/s vibration spectrum
对水翼在各个流速下的振幅进行无量纲化
(y为水翼的振幅，D为水翼的弦长)，得到水翼无量纲振幅随流速变化的趋势曲线，如图12。水翼振幅随流速变化图像呈现倒“V”形。在低流速下，水翼振幅随流速增加而变大。在1 m/s的流速下水翼发生共振，振幅达到最大。流速继续增加，振幅反而减小。
图12 无量纲振幅随流速变化曲线
Fig.12 The variation of dimensionless amplitude with velocity
图13—图16是不同流速下水翼升力和振动位移时程曲线。在流速为0.5 m/s下，水翼的升力曲线与位移曲线的相位相同，即相位差β=0，流场对水翼做正功，但频率低于水翼的自振频率，此时水翼做小幅振动，以耗散流场所输入的能量。当流速增加到1 m/s时，水翼发生涡激共振，此时水翼的振幅在与自振频率同频率的涡激力的作用下达到最大。当流速小幅增加至1.1 m/s时候，升力和振动幅度产生了约
的相位差。涡激力所做的功要比发生涡振时所做的功小，振幅比发生共振时的振幅要小。当来流速度继续增大到1.5 m/s时，水翼的升力曲线与竖向位移曲线相位差正好为π，此时系统脱离了涡振区间。表明：随着流速的增加，发生涡振时的位移响应与涡激力相位逐渐增加。随着涡振区间的发展，水流从由结构振动不断提供能量转化成由结构振动不断消耗能量。
图13 0.5 m/s时升力位移曲线
Fig.13 Lift displacement diagram at 0.5 m/s
图14 1.0 m/s时升力位移曲线
Fig.14 Lift displacement diagram at 1.0 m/s
图15 1.1 m/s时升力位移曲线
Fig.15 Lift displacement diagram at 1.1 m/s
图16 1.5 m/s时升力位移曲线
Fig.16 Lift displacement diagram at 1.5 m/s
以相同方法计算了NACA0009和NACA0015，限于篇幅，不再详细介绍。表3列出了各翼型发生涡振时的流速及无量纲振幅。从表中可以发现如下规律：对于同一阶的涡振现象，厚度越厚的翼型，发生涡振的区间越延后。翼型越薄，发生涡振时的无量纲振幅越小，对于本文计算的二维翼型，振幅甚至出现了数量级的差异。
表3 各翼型涡振流速及无量纲振幅
Table 3 Vortex vibration velocity and dimensionless amplitude of each airfoil
翼型涡振速度/(m·s-1)无量纲振幅00090.80.600121.08.900151.639
6 结论
1) 刚性水翼的泄涡频率与流速呈正比。水翼的形状对泄涡频率的影响较大，在同弦长同流速的前提下，翼型厚度越大，泄涡频率越低。
2) 水翼发生涡振时在区间内振幅最大，此时泄涡频率与结构的固有频率一致，来流流速的变化使涡激力和振动位移产生相位转换现象，决定了涡激振动是限幅振动。
3) 在保证结构强度的前提下，舰艇的舵结构应尽量选取薄翼型。较薄的翼型虽然发生涡振的流速低，但振幅小.
[1] CHEN Y，ZHANG W，FANG D H，et al.Vortex Suppression and Flow Pattern Analysis of a Hydrofoil with Parallel Grooves[J].Processes，2021，9(05):816-816.
[2] HUANG R F，LIN C L.Vortex shedding and shear-layer instability of wing at low-Reynolds numbers[J].AIAA Journal，2010，33(8):1398-1403.
[3] HUERRE P，MONKEWITZ P A.Local and global instabilities in spatially developing flows[J].Annual Review of Fluid Mechanics，1990，22:473-537.
[4] OERTEL H.Jr.Wakes behind blunt bodies[J].Annual Review of Fluid Mechanics，1990，22:539-564.
[5] CHAE E J.Dynamic response and stability of flexible hydrofoils in incompressible and viscous flow[D].Michigan:University of Michigan，2015.
[6] CHAE E J，AKCABAY D T，YOUNG Y L.Dynamic response and stability of a flapping foil in a dense and viscous fluid[J].Physics of Fluids，2013，25(10):104106.
[7] Yu Z X.Hydroelastic calculation in viscous flow field[D].Shanghai:Shanghai Jiao Tong University，1999.
[8] 余志兴，刘应中，缪国平.二维机翼弹簧系统的涡激振动[J].船舶力学，2002(05):25-32.
Yu Z X，Liu Y Z，Mou G P.Vortex induced vibration of two dimensional wing spring system[J].Journal of Ship Mechanics，2002，6(05):25-32.
[9] 侯磊，丁云峰，王晴，等.高雷诺数下水翼涡发放频率预报方法[J].中国舰船研究，2019，14(06):88-97.
Hou L，Ding Y F，Wang Q，et al.Prediction method of vortex shedding frequency of hydrofoil at high Reynolds number[J].Research on Chinese warships，2019，14(06):88-97.
[10] 刘胡涛，张怀新，姚慧岚.水翼涡激振动的数值模拟研究[J].舰船科学技术，2016，38(11):7-13.
Liu H T，Zhang H X，Yao H L.Numerical simulation of vortex induced vibration of hydrofoil[J].Ship science and technology，2016，38(11):7-13.
[11] Guo C Y，Xue R，Hu J，et al.Numerical simulation of tip vortex core motion characteristics of three dimensional hydrofoil[J].Journal of Harbin Engineering University，2021，42(05):618-625.
[12] 贾文隽，华欣，史继拓，等.基于流固耦合的仿生海鸥机翼气动弹性分析[J].兵器装备工程学报，2020，41(07):106-111.
Jia W J，Hua X，Shi J T，et al.Aeroelastic analysis of bionic seagull wing based on fluid structure interaction[J].Journal of Ordnance Equipment Engineering，2020，41(07):106-111.
[13] Zhao Y，Wang G Y，Huang B，et al.Eddy dynamics analysis of unsteady cavitating flow around hydrofoil[J].Ship Mechanics，2015，19(08):896-904.
[14] Chen Y M，Liu J K.Complex normal form method for stability of nonlinear flutter limit cycle[J].Journal of Aerodynamics，2004，22(4):614-617.
[15] Chen Y M，Liu J K.Study on Hopf bifurcation points of both supercritical and subcritical nonlinear flutter systems[J].Applied Mathematics and Mechanics，2008，29(2):181-187.