基于单/双向流固耦合的水翼水力阻尼数值预测方法研究

曾永顺，姚志峰，2，洪益平，2，王福军，2

（1.中国农业大学 水利与土木工程学院，北京 100083；2.北京市供水管网系统安全与节能工程技术研究中心，北京 100083）

1 问题提出

水力机械运行时承受瞬变的水力载荷，满足共振条件时，引起剧烈振动[1-2]。在共振预测时，通常需要已知流固耦合过程中的附加质量和水力阻尼。前者主要影响水下固有频率，且研究较为充分，后者决定了共振幅值大小，但对其了解不够深入[1，3]。主要制约因素是高速流动条件下的水下结构水力阻尼参数测定非常困难[4-5]。近年来，数值计算水平已得到很大的提升，但工程计算依然只能经验性地考虑系统总阻尼，对水力阻尼的定量预测方法还缺少系统的研究。

分离式单/双向流固耦合主要基于计算流体力学（CFD）、计算固体力学（CSD）和计算网格动力学（CMD）等较为成熟的算法，在流体域和固体域分别独立计算，只在流体与固体交界面传递载荷，这已在工程上得到较多应用[6-7]。采用该方法如何预测水力阻尼比是目前国际上水力机械领域研究的热点问题之一。

单向流固耦合方法通过求解振动周期内振动结构的能量耗散，间接计算得到水力阻尼比。Monette 等[6]较早开展了流固耦合条件下水翼与动水的能量交换机理研究，提出了有限元法、模态做功法和单自由度法3种水力阻尼比预测方法。其中，模态做功法得到较多的运用[8-11]。Tengs 等[8]和Bergan 等[9]利用模态做功法得到水翼的第1阶弯曲模态在2.5～45 m/s之间的水力阻尼比，但在低速流动条件下会产生过高预测。Nennemann 等[10]和Gauthier 等[11]提出一种水体附加刚度的预测方法，对模态做功法进行了完善，并建立了水翼第1阶弯曲模态水力阻尼比与流速之间的线性关系，相对偏差在10%以内，但是该模拟没有考虑振型幅值对计算结果的影响。Tengs 等[8]计算得到了水翼在不同振型幅值下的水力阻尼比，结果表明当振型幅值从0.05 mm 增大到2.5 mm时，水力阻尼比显著增大，在10 m/s 流速下相对偏差可达300%。

双向流固耦合方法考虑了流固耦合交界面变形的影响，通过流场非定常计算（CFD）和结构瞬态动力学（CSD）的迭代计算，借助动网格变形，直接获取特定激励下自由振动幅值[7]。该方法相对单向耦合需要更多的计算资源，但在水力阻尼比求解的理论上更加完备。基于该方法，Liaghat 等[12]得到水力阻尼比与流速之间的线性关系，但没有对边界层速度和尾迹区涡激振动进行细致分析，与实验结果相比较，平均相对偏差在12%左右。曾永顺等[13]较为系统地分析了网格、时间步长、近壁区处理模式、激励方式和数值阻尼等关键参数对计算结果的影响，显著提高了水翼水力阻尼比的计算精度，在2～15 m/s 流速下，模拟结果相对实验结果的最大偏差为8.82%。另外，涡激振动[14]和升力作用[15]对水力阻尼比识别精度有重要的影响。针对非对称尾部形状水翼振动响应的时变特性，曾永顺等[16]提出了带通滤波加振动响应校准，对传统自由振动衰减法进行了改进。

在水力阻尼比预测时，目前缺乏对单/双向流固耦合方法的系统比较。本文采用上述两种方法，以NACA0009 钝形尾部形状水翼为分析对象，预测不同流速下的水力阻尼比，分析两种预测方法的可靠性，比较计算精度和耗时，讨论了实际工程中预测方法的选择策略。

2 理论背景

2.1 流场控制方程为了对边界层和尾迹区流动进行精确模拟，流场计算采用转捩SST 模型[17-18]。相对于传统SST k-ω模型，该模型在k 方程中耦合γ-Reθt转捩模型进行修正，控制方程[17]如下：

式中：k为湍动能，m2·s-2；ω为比耗散率，s-1；ρ为流体密度，kg·m-3；τij为壁面剪切应力，N·m-2；μ为动力黏度系数，N·s·m-2；μt为涡黏系数，N·s·m-2；γeff为有效间歇因子（-）；β*和σk为常系数；t为时间，s；u为速度，m·s-1。

2.2 单向耦合水力阻尼比识别方法假设水翼的振型及固有频率不随流速发生改变。获得水翼对流体的模态做功后，根据水翼在静水中的模态质量、固有频率和振型幅值能够计算得到水力阻尼 比[19]：

式中： ζi为第i阶水中模态的水力阻尼比；Wi为第i阶模态一个振动周期内水翼对流体的模态做功，J；mf，i为第i阶水中模态的模态质量，kg；ωnf，i为第i阶水中模态的角固有频率，rad·s-1；ωnf，i=2πfnf，i，fnf，i为第i阶水中模态固有频率，Hz；A0，i为振型的幅值，m；Wn第n个单位时间步内的模态做功，J；N为一个振动周期内的时间步数；P为水体压力，N·m-2；u(t)为振型运动速度，m·s-1；S为水翼表面面积，m2；∆t为时间步长，s。第i阶水中模态的模态质量计算公式为：

式中：ma，i、 fna，i和Ei分别为第i阶空气中模态的模态质量（kg）、固有频率（Hz）和网格单元储存的总动能（J）。

2.3 双向耦合水力阻尼比识别方法双向耦合采用迭代求解方法，直接求解振动衰减曲线，进而识别水力阻尼比。水下结构动力学方程[20]可表示为：

对于自由振动，F (t)=0，求解方程（7）后，结构运动方程[8]可表示为：

式中： A为幅值，m；ωdf，i为第i阶水中模态的阻尼角固有频率，为初相位，rad。

3 计算设置

3.1 物理模型研究对象为NACA0009 钝形尾部形状水翼，物理模型及计算域如图1所示。水翼弦长、翼展和尾部厚度分别为100、150和3.22 mm。水翼材料为不锈钢，密度、弹性模量和泊松比分别为7700 kg/m3、215 GPa和0.3。计算域为750 mm×150 mm×150 mm的方形水洞段，水翼放置攻角为0°。

3.2 边界条件水翼一端被固定约束，另一端自由。边界条件为5～20 m/s的速度进口，2.5 Bar的静压出口。为了让水翼自由端的约束方式和绕流情况与真实实验相符合，将自由端对应的边界设置为对称面，其他边界为无滑移壁面。流场非定常计算的收敛标准为均方根残差小于1×10-5，最大迭代步数为10。非定常计算尾部出现周期性脱落涡后，将该结果作为流固耦合计算的初始文件。基于动网格变形技术，实施单/双向流固耦合计算。

对于单向流固耦合，将水翼水中第1阶弯曲模态的振型导入到流场计算域中，如图1所示。采用差值法将振型赋给流场的网格节点，并根据振型对应的固有频率做周期性振动。将水翼表面单位面积的模态做功进行面积分，累计10个振动周期的结果，取平均后得到1个周期的模态做功。

对于双向流固耦合，为激励第1阶弯曲模态，在水翼表面中心线上沿y 轴正方向施加200 N的瞬态激励，如图1所示，激励时间小于1/4个振动周期。耦合迭代的收敛标准为均方根残差小于1×10-4，最大迭代步数为30。监测流场的网格位移获得结构的振动响应，监测点在尾部中心。

图1 计算域

图2 网格

3.3 计算参数无关性分析采用结构化网格离散流场和结构场，如图2所示。采用O 形拓扑结构划分边界层，第一层网格厚度为2×10-6m，网格拓展比为1.05。流场和结构场网格质量分别在0.35和0.15 以上，且有85%以上的网格质量大于0.8。前期研究[13]中对该计算域进行网格收敛性分析和时间步长无关性检查，结果表明当网格单元数大于246.5 万，时间步长小于2×10-5时对脱落涡频率没有影响；且当结构场网格单元数为1.4 万时，模态分析结果与实验结果相对偏差在0.5%以内。本研究采用上述参数进行流固耦合计算。在20 m/s 流速下，y+平均值为1.29，满足转捩SST模型在近壁区直接求解的要求。

4 结果与讨论

4.1 单向耦合计算结果

4.1.1 模态分析 水翼低阶模态振型及固有频率如图3所示。对于相同模态，空气中振型与水中振型一致，但是在水体附加质量作用下，水中固有频率显著低于空气中。对于不同模态，图3（a）和图3（c）为弯曲变形，前者沿流向只有一个节线，后者有两个节线，两个模态分别命名为第1和第2阶弯曲模态。图3（b）为扭转变形，沿展向只有一个节线，命名为第1阶扭转模态。将模拟结果与实验结果[5]相比较，第1阶弯曲和扭转模态固有频率的相对误差分别为0.44%和0.58%。根据式（5）和式（6），计算得到第1阶弯曲模态在空气中的模态质量为0.2486 kg，在水中的模态质量为0.5815 kg。

图3 低阶模态振型及固有频率

4.1.2 振幅无关性验证 定义相对幅值A∗为振型幅值与第1层网格厚度的比值。10 m/s 流速下，A∗=10和1000时的模态做功特性如图4所示。对于时域，相对振幅越大，水翼对流体的模态做功越大，能量耗散的越多。对于频域，当A∗=1000时频率仅有2倍的第1阶弯曲模态固有频率；当A∗=10时，除该频率为主频外，还有第1阶弯曲模态固有频率和脱落涡频率。结果表明，当相对振幅过大时，脱落涡对水力阻尼比的贡献被忽略。将脱落涡频率模拟结果与实验结果[21]相比较，相对偏差为2.52%。

图4 不同相对幅值的模态做功特性

在10 m/s 流速下，不同相对幅值对应的水力阻尼比如图5所示。随着相对幅值的增大，水力阻尼比逐渐增大，这一趋势与Tengs 等[8]的结果相一致。当相对幅值小于50时对水力阻尼比计算结果没有影响，并将该相对幅值用于不同流速下水力阻尼比的计算。

图5 不同相对幅值的水力阻尼比

4.1.3 单向耦合水力阻尼比 在5～20 m/s 流速下，基于单向流固耦合计算得到的第1阶弯曲模态的水力阻尼比，如表1所示。结果表明随着流速增大，水力阻尼比逐渐增大，这一趋势与文献[5]的实验结果相一致。在模拟流速范围内，相对文献[5]的实验平均偏差为11.42%。

表1 单向流固耦合得到的不同流速下的水力阻尼比

4.2 双向耦合计算结果

4.2.1 振动特性为了消除涡激振动和其他模态的振动对第一阶弯曲模态水力阻尼比识别的影响，采用带通滤波将第一阶弯曲模态的位移从原始位移中分离出来。20 m/s 流速下，水翼振动位移在滤波前后的时域及频域如图6所示。如图6（b）所示，滤波前振动响应有4个频率成分，包括第1阶弯曲、第1阶扭转、第2阶弯曲模态固有频率和脱落涡频率。在外部激励作用下，第1阶弯曲模态的幅值显著大于其他频率。经过滤波处理后，振动响应的频率成分仅剩第1阶弯曲模态固有频率。

图6 滤波前后振动响应（v=20m/s）

在5～20 m/s 流速下，固有频率及脱落涡频率随流速的变化如图7所示。第1阶弯曲、扭转及第2阶弯曲模态固有频率基本不随流速发生改变，相对变化范围分别在0.43%、2.61%和3.95%以内。对于第1阶弯曲和扭转模态的固有频率，与文献[5]的实验结果相比较，平均相对偏差分别为4.36%和4.07%。脱落涡频率随流速线性增长，相对文献[21]的实验结果平均偏差为4.24%。模拟流速对应的脱落涡频率均大于第1阶弯曲模态固有频率，错开了共振工况。

4.2.2 双向耦合水力阻尼比 20 m/s 流速下，水力阻尼比的识别如图8所示。将滤波后振动响应的上下峰值点取出，根据自由振动衰减法进行函数拟合。基于上下峰值点的水力阻尼比分别为0.076 78和0.079 22，二者的相对偏差为3.18%，将二者取平均后得到最终的水力阻尼比为0.078。为避免滤波对振动幅值的影响，在水力阻尼比识别时舍弃第一个振动响应的峰值点。

图7 不同流速下振动响应的频率

在5～20 m/s 流速下，水力阻尼比随流速的变化如图9所示。数值模拟计算得到的水力阻尼比随流速线性增长的趋势与文献[5]的实验结果相一致，双向流固结果与实验结果相比较，平均相对偏差为4.95%。在10 m/s 流速下，单/双向流固耦合计算得到水力阻尼比与实验结果相比较，相对偏差分别为17.07%和0.64%。原因是该流速接近共振工况，此时结构振动对流场流动有较为强烈的反作用影响，而双向流固耦合数值模拟与真实流动更加吻合，因此模拟精度显著高于单向流固耦合。

图8 水力阻尼比的识别方法

图9 不同流速下的水力阻尼比

4.3 讨论对比单/双向流固耦合数值模拟方法，前者假设振型和固有频率不随流速发生改变[11]，且不考虑结构变形的反作用影响[8]，而后者理论上更加完善。在20 m/s 流速下，基于双向流固耦合计算得到的第一阶弯曲模态的振型及固有频率如图10所示。与静水中模态分析结果相比较，振型相一致，固有频率相对偏差为3.65%。在低阶模态下，水翼固有频率在5～20 m/s 流速下的相对变化量在3.95%以内，这也证明了单向流固耦合方法的假设的合理性。但是，在大变形情况下，单向流固耦合忽略了涡激振动对水翼能量耗散的贡献，会导致水力阻尼比的预测误差过大，此时需要双向流固耦合进行求解。

图10 20 m/s 流速下第1阶弯曲模态振型

基于相同的计算设置，在10 m/s 流速下采用48 线程进行第1阶弯曲模态的水力阻尼比预测，单/双向流固耦合计算一个振动周期分别需要2.1 h和31.3 h，后者需要的计算资源为前者的15倍。实际工程中，水力机械结构较为复杂，优先推荐单向流固耦合方法预测水力阻尼比；若在大变形工况下，或需要高精度的预测结果，适当简化工程对象，可采用双向流固耦合方法。

5 结论

采用单/双向流固耦合数值模拟方法，在5～20 m/s 流速下对NACA0009 钝形尾部水翼的水力阻尼比进行了预测，验证了两种计算方法的可靠性，对比分析了两种方法的计算精度和计算时间，并对实际工程中水力阻尼比预测方法的选择进行了讨论。主要结论如下所示：（1）静水中水翼第1阶弯曲模态与20 m/s 流速下对应振型一致，且低阶模态固有频率在5～20 m/s之间的相对变化在3.95%以内。据此，可假设振型及固有频率与流速无关，开展单向流固耦合的水力阻尼比预测。（2）单/双向流固耦合数值模拟计算得到的第1阶弯曲模态水力阻尼比随流速线性增长的趋势与实验一致，平均相对偏差分别为11.42%和4.95%，后者计算精度高于前者，但计算资源需求量是前者的15倍。（3）工程上，预测复杂水力机械的水力阻尼比时，优先推荐单向流固耦合方法，若需要高精度或结构产生大变形时，则需要简化分析对象，并采用双向流固耦合方法。