随着能源和调水需求的增加，水力机械向着大功率、大尺寸、高转速的方向发展，其运行稳定性问题更加突出^[1]。水力机械运行时，受到内部复杂流场激振力作用，不可避免地发生一定程度的振动^[2]。当激振力频率与结构固有频率相近时，可能发生共振，造成重大损失^[3]。水下结构振动特性取决于流场-结构的耦合特性，特别是附加质量和阻尼效应^[4]。水利机械偏离设计工况运行时，在叶片流线型表面的前缘位置可能出现附着型空化区域，即为“前缘空化”。前缘空化是水力机械的主要空化类型之一，它显著降低了结构部件周围流体介质密度，使水体附加质量和附加阻尼发生变化，从而对结构场振型产生影响。因此，探究并预测空化条件下的结构动力响应具有重要意义。

国内外不少学者通过实验和数值模拟方法研究了水翼前缘空化的发展形态和流场特性^[5-7]。Dupont^[5]测量了NACA0009水翼在不同空化数和Reynold数下的水翼流场分布、空泡发生位置及形态。De La Torre等^[6]使用激光测振仪测量了相同水翼结构在不同空化数下的振型及固有频率。本文作者及其同事^[7]针对相同水翼，测试了无空化条件下不同流速时的附加阻尼比。以上实验研究为本文的数值模拟提供了验证数据。

当结构浸没在水下时，结构运动将受到周围介质的阻碍。相比空气，水中结构的动态响应会发生较大改变^{[1, 4]}。周围流体的阻碍作用主要体现在使振动系统增加了附加质量和阻尼。在附加质量的作用下，水中结构振动的固有频率降低^[8]；在附加阻尼的作用下，振动幅值衰减加快^[1]。但是，目前对于附加质量和阻尼的研究主要集中在无空化条件下。在附加质量研究方面，朱文若等^[8]运用流固耦合的有限元方法对离心式叶轮水中固有频率经验下降系数进行了分析及优化；Lindholm等^[9]测量了悬臂梁结构在空气中、完全以及部分浸没在水中的固有频率，并建立了附加质量效应的经验公式；Blevins等^[10]的研究表明，水中振动结构体的附加质量与结构体的外形、结构体与边界条件的相对位置、振动幅值和振动方向以及量纲为1的参数fD²/ν (f为结构场振动频率，D为特征参数，ν为流体运动黏度)4个因素相关。在附加阻尼研究方面，近十年的研究重点主要在流速与附加阻尼的关系上：Seeley等^[11]实验测定了附加阻尼随流速增加线性增大；本文作者及其同事^[7]的实验研究表明，水翼附加阻尼在低速区和高速区的表现并不一致，分界点为共振流速点，在共振流速点以下，附加阻尼随流速几乎保持不变，在共振流速点以上，附加阻尼随流速线性增加；Bergan等^[12]的实验也验证了上述结果。

到目前为止，对于空化发生后如何影响附加质量和阻尼的问题尚不明确。特别是，由于空化本身强非线性、宽频域和局部瞬态冲击的激振源特性，使得水中结构动态响应实验异常困难。因此，不少学者通过数值模拟来寻求解决这类问题的途径。Liu等^[13]拓展了声固耦合方法，数值模拟结果表明水翼固有频率随着空泡面积的增加而增大。Hu等^[14]采用紧耦合的数值模拟研究方法表明，水翼微小变形对局部空化的影响有限。曾永顺等^[15]通过双向流固耦合数值模拟方法揭示了NACA0009水翼在非空化条件下的固有频率和附加阻尼随流速的变化规律，但并未涉及空化场的计算。

基于上述研究，本文以NACA0009弹性水翼为研究对象，利用分离式双向迭代的流固耦合方法，开展不同前缘空化程度下的水翼弯曲和扭转模态的激振响应数值模拟，获得弹性振动水翼的固有频率与附加阻尼比，以揭示前缘空化对弹性水翼振动特性的影响。

1 数值计算理论 1.1 流场控制方程

为了兼顾计算精度和所需计算资源，本文选择基于剪应力输运(shear stress transport, SST) k -ω湍流模型的Reynold平均方法(RANS)；为了模拟空化多相流场，采用Euler-Euler方法中的均相流模型。流场计算的连续性方程和动量方程分别为^[16]：

$ \frac{\partial \rho_{\mathrm{m}}}{\partial t}+\frac{\partial\left(\rho_{\mathrm{m}} u_{j}\right)}{\partial x_{j}}=0, $

(1)

$ \begin{gathered} \frac{\partial\left(\rho_{\mathrm{m}} u_{i}\right)}{\partial t}+\frac{\partial\left(\rho_{\mathrm{m}} u_{i} u_{j}\right)}{\partial x_{j}}=-\frac{\partial p}{\partial x_{i}}+ \\ \frac{\partial}{\partial x_{j}}\left[\mu_{\mathrm{m}}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}-\frac{2}{3} \frac{\partial\left(u_{i}\right)}{\partial x_{j}} \delta_{i j}\right)\right]+S_{\mathrm{u}}. \end{gathered} $

(2)

其中：t为时间，s；x为空间坐标，m；下标i和j代表坐标方向；u为速度，m/s；p为压强，Pa；ρ_m为混合介质密度，m³/s；μ_m为混合湍流黏性系数，Pa·s；δ_ij为判定符号，当i=j时，δ_ij=1，当i≠j时，δ_ij=0；S_u为源项。

1.2 空化模型

采用Zwart-Gerber-Belamri(ZGB)空化模型模拟相间质量传输，相间质量传输率可表示为^[17]：

$ \dot{m}^{-}=F_{\text {vap }} \frac{3 \alpha_{\text {nuc }}\left(1-\alpha_{\mathrm{v}}\right) \rho_{\mathrm{v}}}{R_{\text {nuc }}} \sqrt{\frac{2}{3} \frac{\left|p_{\mathrm{v}}-p\right|}{\rho_{\mathrm{f}}}}, $

(3)

$ \dot{m}^{+}=F_{\text {cond }} \frac{3 \alpha_{\mathrm{v}} \rho_{\mathrm{v}}}{R_{\mathrm{B}}} \sqrt{\frac{2}{3} \frac{\left|p_{\mathrm{v}}-p\right|}{\rho_{\mathrm{f}}}}. $

(4)

其中：F_vap为蒸发项系数，默认值为50；凝结项系数F_cond=0.01；α_v为蒸汽体积分数；p_v为环境温度下液体汽化压强；ρ_v为蒸汽密度，m³/s；ρ_f为液体密度，m³/s；核子半径R_nuc=1 μm；核子位置体积分数为α_nuc=0.000 5；R_B为空泡半径，m。

采用空化数σ作为定量表征空化环境的参数，

$ \sigma=\left(p-p_{\mathrm{v}}\right) /\left(0.5 \rho_{\mathrm{f}} u^{2}\right). $

(5)

1.3 声固耦合方法

对于水中模态分析，将Navier-Stokes(N-S)方程简化为Helmholtz方程(声波方程)^[18]，

$ \frac{1}{v_{\mathrm{w}}} \frac{\partial^{2} p}{\partial t^{2}}=\nabla^{2} p. $

(6)

其中v_w为水中声速，m/s。

将方程(5)离散，并耦合空气中结构动力学方程，可得声固耦合方法的有限元离散方程，

$ \begin{array}{c} &{\left[\begin{array}{cc} \boldsymbol{M}_{\mathrm{s}} & {\bf{0}} \\ \boldsymbol{M}_{\mathrm{fs}} & \boldsymbol{M}_{\mathrm{f}} \end{array}\right]\left[\begin{array}{c} \ddot{\boldsymbol{y}} \\ \ddot{\boldsymbol{p}} \end{array}\right]+\left[\begin{array}{cc} \boldsymbol{C}_{\mathrm{s}} & {\bf{0}} \\ {\bf{0}} & \boldsymbol{C}_{\mathrm{f}} \end{array}\right]\left[\begin{array}{l} \dot{\boldsymbol{y}} \\ \dot{\boldsymbol{p}} \end{array}\right]+} \\ &{\left[\begin{array}{cc} \boldsymbol{K}_{\mathrm{s}} & \boldsymbol{K}_{\mathrm{fs}} \\ {\bf{0}} & \boldsymbol{K}_{\mathrm{f}} \end{array}\right]\left[\begin{array}{c} \boldsymbol{y} \\ \boldsymbol{p} \end{array}\right]=\left[\begin{array}{c} \boldsymbol{F} \\ {\bf{0}} \end{array}\right]}. \end{array} $

(7)

其中：M、C和K分别为质量、阻尼和刚度矩阵；下标s和f分别代表结构和流体；M_fs和K_fs分别为耦合系统等效质量和刚度矩阵；y为位移，m；F为外载荷，N。

1.4 双向流固耦合方法

耦合计算流体动力学(CFD)和计算结构动力学(CSD)，对瞬态结构振动响应进行分离式求解。流固耦合交界面同时满足运动学条件和动力学条件^[19-20]：

$ \int(\boldsymbol{p n}+\tau) \varphi(y) \mathrm{d} s=F_{\mathrm{m}}. $

(8)

其中：τ为黏性应力，Pa；φ (y)为模态振型；s为面积，m²；F_m为模态力，N。

1.5 附加质量和附加阻尼识别方法

水下结构动力学方程可表示为^[21]

$ \left(m_{\mathrm{s}}+m_{\mathrm{f}}\right) \ddot{y}+\left(c_{\mathrm{s}}+c_{\mathrm{f}}\right) \dot{y}+\left(k_{\mathrm{s}}+k_{\mathrm{f}}\right) y=F(t). $

(9)

其中：m、c和k分别代表质量、阻尼和刚度系数。

水下圆周固有频率和附加阻尼比ζ可分别表示为：

$ \omega_{\mathrm{nw}}=\sqrt{\frac{k_{\mathrm{s}}+k_{\mathrm{f}}}{m_{\mathrm{s}}+m_{\mathrm{f}}}}, $

(10)

$ \zeta=\frac{c_{\mathrm{s}}+c_{\mathrm{f}}}{2 \sqrt{\left(k_{\mathrm{s}}+k_{\mathrm{f}}\right)\left(m_{\mathrm{s}}+m_{\mathrm{f}}\right)}}. $

(11)

将附加质量系数定义为水体附加质量与结构质量的比值，若忽略水体附加刚度，可表示为^[22]

$ C_{\mathrm{M}}=\frac{m_{\mathrm{f}}}{m_{\mathrm{s}}}=\frac{\omega_{\mathrm{na}}^{2}}{\omega_{\mathrm{nw}}^{2}}-1. $

(12)

其中：ω_na为空气中固有角频率，rad/s；ω_n=2πf_n，f_n为固有频率，Hz。

将式(10)、(11)代入式(9)，求解后可得结构运动方程^[21]，

$ y(t)=A \mathrm{e}^{-\zeta \omega_{\mathrm{n}}} \sin \left(\omega_{\mathrm{d}}+\varphi\right). $

(13)

其中：A为幅值，m；ω_d为有阻尼固有角频率， ${\omega _{\rm{d}}} = {\omega _{\rm{n}}}\sqrt {(1 - {\zeta ^2})} $，rad/s；φ为初相位，rad。

2 计算方案 2.1 物理模型

NACA0009弹性水翼弦长L=100 mm，展向长度w=150 mm，尾部厚度h=3.22 mm，材料为铝合金，密度为2 770 kg/m³，弹性模量为7.1×10¹⁰ Pa，Poisson比为0.33。物理模型及计算域如图 1所示。水翼攻角为2°，计算域进口距水翼前缘2.5L，出口距尾部4L，底部距水翼0.8L，均与空化水洞测试段尺寸保持一致^[6]。本文的双向流固耦合计算未考虑间隙对水翼动力响应特性影响，水翼展向长度w与计算域宽度相等，均为150 mm。

2.2 网格与边界条件

对水翼绕流流体域和结构体进行结构化网格划分，在水翼周围采用O-Block，并对水翼壁面及测试段上下壁面进行网格加密，网格近壁区扩展比为1.05。为消除计算域网格对计算结果的影响，基于Richardson外推加速法对网格收敛性指数^[23]进行定量分析，结果表明当流场网格数大于246万个时对流场内部流动预测没有影响。对于结构场进行网格无关性检查发现，当网格数大于1.4万个时对固有频率预测没有影响。基于相同水翼，该网格尺度与Zeng等^[21]的一致，此时流场壁面平均y⁺≈1，符合SST k-ω湍流模型计算要求。采用六面体网格，通过计算每一个六面体的Jacobi行列式值并利用标准化行列式的矩阵来表征单元变形，其值在0.3以上可满足大多数流动现象。本算例中水翼头部和尾部网格质量见图 2，网格质量都在0.8以上。

流场瞬态项采用二阶Euler后差分格式，对流项及湍流模型相关量均采用高阶精度求解。流体域进口采用速度进口，流速为14 m/s，与文[6]的实验条件一致；出口采用静压边界条件，压强值通过空化数计算得到；计算域上下壁面采用无滑移边界条件，水翼固定端和自由端采用的边界条件分别为无滑移壁面和对称面。根据Zeng等^[21]以旋涡脱落频率为关键参数进行的时间步长无关性检验结果，本研究中流场和结构场时间步长均取2×10^-5 s。单步迭代次数最高为30次，数值阻尼为0.01。

2.3 激励方式与振动监测点

为获得水翼弯曲模态、扭转模态下固有频率、附加质量与阻尼，在双向流固耦合计算中对水翼施加激振力。对于弯曲模态，沿水翼表面中心位置施加y轴正方向的激振力，大小为1 200 N，其方向与水翼所受升力方向相同，如图 3a所示。对于扭转模态，在水翼前缘及尾部施加垂直于翼弦且方向相反的激振力，大小为1 200 N。由于本文水翼存在2°攻角，在流固耦合计算初始时在升力作用下发生变形，为排除水翼所受升力对施加激振力的影响，在水翼尾部监测点振动稳定后，再对水翼施加激振力，并对施加激振力后的振动响应进行分析，激振力施加时刻及持续时间如图 3b所示。施加激振力的时间需小于该模态下振动周期的1/4，即1/(4f)，f为振动频率，Hz。对于弯曲模态，激振力持续时间为1.8×10^-3 s；对于扭转模态，激振力持续时间为3.8×10^-4 s。在水翼尾部布置等距离的7个监测点，用于获取水翼的振动响应。另外，测点8为流场压强脉动监测点。

2.4 计算工况

为探究空化数对水翼动力响应特性的影响，使用双向流固耦合方法计算了静水中水翼的弯曲模态、扭转模态激振响应，并与后续流动和空化条件下的计算结果进行对比。选取4个不同空化数的工况进行双向流固耦合计算，在空化数σ分别为5.07(即无空化发生)和0.80、0.62、0.40(有空化发生)工况下，对弯曲和扭转模态进行模拟，得到其振动响应。

3 结果与讨论 3.1 流场结果分析 3.1.1 前缘空化形态分析

图 4为空化数分别为5.07、0.80、0.62和0.40时的空化区域及绕流流场。可以看出，在空化数σ=5.07时，水翼绕流流场无空化现象；随着空化数的降低，水翼前缘出现附着在上壁面的空泡，空泡长度随空化数的降低而增加，且空化区域内气体体积分数也随空化数的降低而增大。当σ=0.40时，空泡长度超过水翼的尾部，在x/L=0.7至x/L=1.1区域存在大尺度旋涡。

计算得到不同空化数下水翼前缘空化形态定量结果，将空泡长度l_cavity及高度h_cavity与水翼长度L作比，进行量纲归一化，得到相对长度及高度，如表 1所示。前缘空化长度和高度均随空化数的降低而增加。其中，空化数σ=0.62时，文[5]实验得到的相对空泡长度为0.754，本文模拟值为0.79，相对误差为4.8%。

3.1.2 尾迹区涡量分布

在水翼展向z=0.05 m处作截面，得到该截面z方向涡量云图及水翼绕流流线，如图 5所示。可以看出，当空化数σ为5.07、0.80和0.62时，水翼尾部从上壁面脱落负向旋涡，从水翼下壁面脱落正向旋涡，交替形成Karman涡街；当空化数σ=0.40时，由于空泡长度超过水翼尾部，水翼上侧负涡与壁面分离，沿流动方向并未形成Karman涡街。水翼上壁面负涡区域面积随空化数的降低而增大，下壁面正涡区域则无明显变化。随着空化数的降低，空泡长度增大，水翼尾迹区旋涡的长度也同时增大。

3.2 附加质量和附加阻尼 3.2.1 振型

通过在动水环境中施加2种类型激振力，获得了所研究水翼的弯曲和扭转振型，见图 6。

3.2.2 弯曲模态

图 7给出了静水中及空化数σ分别为5.07、0.80和0.62工况下的水翼激励响应结果。以水翼施加激振力的时刻为时间起点0时刻，采用带通滤波将固有频率振动响应分离出来，并拟合振动响应的上下峰值点，从而得到附加阻尼比，如图 7所示。

表 2给出了弯曲模态下的固有频率和附加阻尼比，并根据固有频率计算得到弯曲模态的附加质量系数，并与文[6]的实验值进行比较。由表 2可知，静水中附加质量系数模拟值和实验值^[6]的相对误差为6.0%，空化数σ=5.07时，相对误差为2.2%，当空化数σ=0.62时，相对误差为16.4%，均在可接受范围内。空气中水翼固有频率为283.9 Hz，当计算域介质为静水时，由于水体附加质量的影响，水翼固有频率降低，为133.7 Hz。当计算域进口流速为14 m/s时，水翼弯曲模态固有频率较静水中增大，附加质量系数由静水时的3.31降低至3.19。随着空化数的降低，附着在水翼表面的空化区域逐渐增大，显著降低了围绕在水翼周围的介质密度，水翼固有频率进一步增大，附加质量系数减小。具体地，在空化出现后，两次空化数的降低分别使水翼固有频率增加了2.25%和0.21%；附加阻尼比分别降低了2.34%和9.00%。

3.2.3 扭转模态

采用与弯曲模态相同的附加阻尼识别方法，对扭转模态监测点振动曲线进行分析，如图 8所示。

表 3给出了水翼扭转模态的固有频率、附加阻尼比和附加质量系数，并与文[6]的实验值进行对比。静水中附加质量系数相对误差为4.0%，空化数σ=5.07时相对误差为0%，空化数σ=0.62时相对误差为32%。水翼在空气中的扭转模态固有频率为1 053.9 Hz；当介质为静水时，由于水体附加质量的影响，固有频率降低至645.2 Hz。当计算域中流速增至14 m/s时，固有频率较静水中有所增大，随着空化数的降低，前缘空化覆盖水翼面积增大，水翼固有频率明显增大，附加质量系数由无空化时的1.58降低至0.68，在静水中，水翼附加阻尼比为0.010 82，流速增大至14 m/s时，附加阻尼比也较静水中增大。当空化数降低、水翼前缘有空化产生时，监测点振动信号在激振力撤销后迅速衰减，并恢复至施加激振力前的振动状态，附加阻尼比分别为0.081 81和0.091 63，较无空化产生时显著增大。具体地，在空化出现后，两次空化数的降低分别使水翼固有频率增加了20.12%和3.25%；对应的附加阻尼比的增幅分别为165.70%和12.00%。

总体而言，前缘空化对水翼固有频率和附加阻尼比具有显著影响，其影响方式与振型有关。双向流固耦合方法已经能够较好地模拟振动水翼的空化流场和激振特性，为下一步深入研究多叶片结构的真实水力机械叶轮提供了有效手段。

4 结论

为了研究水力机械叶片在前缘空化条件下的振动特性，将叶片简化为水翼开展研究。基于SST k-ω湍流模型和ZGB空化模型的RANS方法预测了前缘空化流场。利用双向流固耦合方法获得了弯曲和扭转振型激励下的水翼自由振动衰减振动信号，分析得到振动水翼的固有频率与附加阻尼比，主要结论如下：

1) 随着空化数的降低，水翼前缘空化长度、厚度和边界层厚度均增大。随着空泡长度的增加，水翼尾部脱落涡频率降低，当空泡长度超过水翼尾部时，不再出现交替脱落的旋涡。相比静水环境，动水中水翼的固有频率稍有增大。

2) 空化系数变化对弯曲振型影响较小，但对扭转振型的影响较大。随着前缘空化区域的增大，弹性水翼弯曲模态与扭转模态的固有频率都逐渐增大，但附加阻尼比的变化依赖于振型。

3) 当弹性水翼在攻角2°、流速14 m/s的工况运行时，前缘空化出现后，且空化区域长度未超过水翼长度，弯曲振型下水翼附加阻尼比可减少9.00%，而扭转振型下附加阻尼比的最大增幅为165.70%。