基于嵌套网格的对称水翼流致振动数值模拟

姚志峰，邹圣可，曾永顺，刘宝熙，张佳云

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

可逆式水泵水轮机是抽水蓄能电站的核心装备，运行时需要频繁切换工况，不可避免在非设计工况运行，导叶或叶片将承受较大幅值的动态水力载荷.如何在设计阶段预测导叶或叶片等结构振动特性是目前水泵水轮机设计的难点问题[1-3].水翼可近似为导叶或叶片的核心工作单元.水翼绕流过程会在尾部产生交替脱落的卡门涡街[4-6]，水翼结构受脱落频率影响产生受迫振动，其脱落频率和流速通常满足斯特劳哈尔定律.当脱落频率接近结构固有频率时，水翼振动频率在某段流速范围内不再跟随旋涡脱落频率，而是保持其固有频率振动，即为“锁频”现象.此时，水翼振动幅值急剧增强，可达非“锁频”时的数百倍量级[5-6].

已有不少学者研究了水力机械叶片、导叶等部件的涡激振动现象[7-9].基于高速水洞实验，瑞士洛桑联邦理工大学AUSONi 等[4]获得了5～30 m/s 流速下钝型尾部对称水翼的结构固有频率和旋涡脱落频率，并发现当旋涡脱落频率接近水翼的第一阶扭转模态频率时，噪声和振动水平显著增加，旋涡脱落频率锁定在第一阶扭转模态频率上，流速范围为11～13 m/s.基于数值模拟方法，挪威科技大学TENGS 等[10]采用准双向流固耦合方法获得了水翼5～25 m/s 流速下的振动频率曲线和振动幅值曲线，并发现当流速范围为10～12 m/s 时，结构振动水平加剧，振动频率锁定在一个固有频率上，而不是随着流速线性变化.

已有研究表明，借助分离式双向迭代流固耦合方法可模拟水翼结构动力响应.LIAGHAT[11]等利用该方法，通过对水翼施加瞬态激励力，结合频谱分析和自由振动衰减法，获得了水翼结构振动固有频率和旋涡脱落频率，得到水力阻尼比与流速之间的线性关系.曾永顺等[12-14]采用类似方法系统分析了网格、时间步长、近壁区处理模式、激励方式和数值阻尼等关键参数对计算结果的影响，显著提高了振动水翼水力阻尼比的计算精度.

但受限于湍流模型和网格尺度等影响，以上分离式流固耦合数值计算都面临结构变形过大导致出现负体积网格的问题，因此研究只限于小变形的非锁频流速.

嵌套网格技术可以解除物体与网格之间在拓扑结构上的约束关系，能够让各物体在计算域内自由地运动，不受网格结构的影响[15-16].与传统动网格划分方式相比，嵌套网格对计算模型进行分块处理，各子区域独立生成网格，降低了计算网格的生成难度；而且嵌套网格在计算过程中，不会产生网格缠绕、畸变和负体积等问题，能始终保持计算网格的质量，从而保证计算结果的准确性[15].

本文将基于嵌套网格方法对NACA 0009 对称尾部形状水翼进行双向流固耦合数值模拟，得到其在锁频区和非锁频区流速下的振动响应.验证了基于声学单元的模态分析方法，验证嵌套网格方法在水翼流场计算中的适用性，分析了锁频工况和非锁频工况下水翼流场特征.

1 数值模拟理论

1.1 控制方程

流体运动由连续性方程和动量方程控制[17]：

式中：i，j为张量角标（取1，2，3）；ui为流体的速度；ρ为流体密度；t为时间；p为压力； µ为动力黏性系数；fi为作用于流体上的质量力.

结构动力学方程如下[18]：

式 中：y为结构振动位 移；y˙ 为 结 构 振 动速度；y¨为结构振动加速度； ζ为阻尼比； ωn为角频率；fn为固有频率；F(t)为外载荷；M、C和K分别为质量、阻尼和刚度系数.

在流体与结构交界面上，应满足动力学和运动学条件[19]：

式中： τ为应力；n为外法线单位矢量的分量；d为位移；下标f 和s 分别表示流场和结构场.

1.2 湍流模型

本文考虑水翼表面的边界层转捩.采用 γ-Reθ转捩SST 模型，能有效地捕捉边界层转捩等细致流场结构.对比传统SSTk-ω模型，该模型在k方程中耦合 γ-Reθ转捩模型进行修改，方程如下[20-21]：

式中：k为湍动能； ω为比耗散率； ρ为流体密度； τij为壁 面 剪切 应 力； µ为 动力 黏 度 系数； µt为 涡 黏系 数；γeff为有效间歇因子； β*和 σk为常系数；t为时间；u为速度.

1.3 嵌套网格概述

嵌套网格由背景网格和组件网格组成，通过挖洞、嵌套区域最小化和供体网格搜寻等步骤建立起两套网格之间的插值关系，形成插值区域，两套网格在插值区域内通过插值完成数据传递.由上述步骤产生死亡单元、受体单元、供体单元和求解单元，在插值区域中，受体单元从多个供体单元接收信息[15-16].

网格嵌套处理时，如果嵌套区域过于小，在网格尺度一定的条件下，嵌套区域内没有足够的网格，则会出现一套网格的受体单元与另一套网格的受体单元在几何位置上相对应的情况，此类网格称为孤儿网格，会对计算精度产生影响；两套网格在嵌套交界面处的网格尺寸差异较大也会产生孤儿网格.

1.4 模态分析理论

采用完全流固耦合方法计算水中水翼湿模态，不考虑水的黏性作用和湍流流动对水翼的影响，该方法将流场控制方程简化为声波方程，并采用声学单元代替流场[22-23]：

式中：p为水体压力；v为水中声速；t为时间.

将结构动力学方程与方程（9）耦合，可表示为：

式中：Ms、Cs和Ks分别为结构的质量、阻尼和刚度矩阵；Mf、Cf和Kf分别为水体附加的质量、阻尼和刚度矩阵；Mfs和Kfs分别为耦合系统的质量和刚度矩阵；下标f 和s 分别表示流场和结构场.

2 数值模拟方案

2.1 计算模型

数值模拟对象见图1，NACA0009 水翼以0°攻角放置在 150 mm×150 mm×750 mm的矩形水槽中，水翼弦长L=100 mm，水翼最大厚度h1=10 mm，尾部厚度h2=3.22 mm，展向宽度w=150 mm.该水翼材料为不锈钢，材料密度 ρ=7 900 kg/m3，弹性模量E=200 GPa，泊松比v=0.3.水翼一侧为固定支座，另一侧为枢轴嵌入.在水翼上表面设置振动监测点1，监测点位置与实验[4]保持一致，监测点坐标为（x/L，z/w）=（0.80，0.75）；在流场中设置压力脉动监测点2 获取旋涡脱落频率，监测点坐标为（x/L，z/w）=（1.1，0.5）.

图1 几何模型Fig.1 Geometricl model

2.2 模态分析设置

为了确定水体环境中水翼的固有频率和振型，采用ANSYS Modal Acoustics 开展水翼湿模态分析，采用Full Damped 方法进行求解，得到水翼在水中的各阶模态.在湿模态分析的计算域尺寸为750 mm（长）×150 mm（宽）×150 mm（高），水翼攻角为0°，计算域出口和进口设置为全吸收面，其他壁面则为全反射面.

2.3 双向流固耦合设置

采用ANSYS Workbench 平台耦合Fluent 和Transient Structural 实现双向流固耦合计算.进口边界条件采用速度进口，出口边界条件为压力出口，数值为250 kPa，流体与水翼相接触的面设置为流固耦合交界面，其余面均为无滑移壁面.对于结构场，采用有限元法进行瞬态动力学分析.对于流场采用有限体积法求解，以流场非定常计算结果为流固耦合计算的流场初始文件.流固耦合迭代收敛标准为 1×10-4，每个时间步长最大迭代次数为20，时间步长Δt=2×10-5s.通过记录流场动网格变形或者结构场网格变形来获取水翼的振动响应.

流场网格见图2，分为背景网格和组件网格，网格类型均采用正六面体结构化网格.在背景网格中，对组件区域所在几何位置处进行局部加密.为了有效捕捉层流和转捩边界层，必须保证流场水翼表面y+≈1，经过计算后得到水翼表面第一层网格厚度为2×10-3mm，为保证流场中组件网格的质量，对水翼前缘和尾部进行局部加密.

图2 流场网格Fig.2 Mesh for fluid field

3 关键计算参数分析

3.1 嵌套区域大小分析

为排除嵌套区域大小对数值计算结果的影响，比较了6 种不同大小的部件网格区域的计算效果.部件网格区域为矩形，考虑到计算效率和硬件限制的情况下选取最优部件网格区域大小.不同嵌套区域大小参数见表1，嵌套区域相对水翼位置及大小如图3 所示.

表1 嵌套区域尺寸及网格数Tab.1 Overlapping area size and its grid number

图3 不同嵌套区域大小示意图Fig.3 Diagram of different overlapping area sizes

在进行嵌套区域大小验证时，进口流速v=20 m/s.取相同时间内水翼尾部监测点的绝对压力进行快速傅里叶变换，得到该监测点旋涡脱落频率.以旋涡脱落频率为关键参数进行验证，验证结果如图4 所示.

图4 不同嵌套区域大小验证结果Fig.4 Verification results of different overlapping area sizes

可见，随着嵌套区域扩大，旋涡脱落频率逐渐趋于稳定且接近于实验结果.考虑到计算结果准确性，确定嵌套区域长度为L=1.4l，嵌套区域宽度H=4h1为最优方案,与实验结果相对偏差为1.18%.

3.2 背景网格和组件网格尺寸差异分析

为了保证水翼近壁区计算的准确性，组件网格的尺寸保持不变，改变背景网格的尺寸来验证两套网格尺寸差异对计算结果的影响, 如图5 所示.

图5 嵌套网格尺寸差异Fig.5 Overlapping grid size difference

采用上一节确定的嵌套区域，在嵌套交界面处通过调整背景网格尺寸来控制两套网格的尺寸差异.嵌套交界面处网格尺寸差异见表2.

表2 嵌套网格尺寸差异Tab.2 Overlapping grid size difference

取相同时间内水翼尾部监测点的绝对压力值进行快速傅里叶变换，图6 为背景网格和部件网格尺寸差异对旋涡脱落频率的影响.

图6 背景网格、部件网格尺寸差异对计算结果影响Fig.6 The effect of mesh size differences on the calculation results

当背景网格和部件网格尺寸差异过大，旋涡脱落频率计算不准确，与实验值对比有一定偏差.随着两套网格尺寸差异减小，计算结果与实验结果接近.

在确定背景网格和组件网格的网格尺寸大小时，应保证组件区域外层网格和对应几何位置处的背景网格尺寸大小一致，以确保计算结果的准确性.

4 结果与讨论

4.1 模态分析

本文实验结果来自瑞士洛桑联邦理工大学空化实验水洞，采用激光测振仪测量水翼特定位置点的振动[4]，振动监测点如图7 所示.

图7 实验振动监测点布置[4]Fig.7 Arrangement of experimental vibration monitoring points[4]

模态分析结果如图8～图10 所示.图8 是水翼第一阶弯曲模态，振幅随着从固定端到枢轴嵌入侧的跨度增加而增加，在靠近枢轴嵌入侧又逐渐减小.图9是水翼第一阶扭转模态，在水翼弦长的中间截线区域振幅约为0，最大振幅出现在水翼前缘和尾缘.图10是水翼第二阶弯曲模态，在靠近固定端和z/w=0.66处振幅约为0，随着远离这两个位置的跨度增加，振幅逐渐增加.

图8 第一阶弯曲模态实验、模拟对比Fig.8 Experimental and simulation comparison of first-order bending mode

图9 第一阶扭转模态实验、模拟对比Fig.9 Experimental and simulation comparison of first-order torsional mode

图10 第二阶弯曲模态实验、模拟对比Fig.10 Experimental and simulation comparison of second-order bending mode

对比模态分析结果和实验结果[4]，振型基本一致.表3 将数值模拟结果和实验结果作比较.对比前三阶模态，模拟结果和实验结果相对偏差均在4%以内，验证了模态分析方法的可靠性.

表3 模拟结果和实验结果比较Tab.3 Comparison of simulation results and experimental results

4.2 流致振动分析

4.2.1 水翼边界层速度分布

对来流速度v=20 m/s 和v=12 m/s 进行边界层速度分布分析.图11 为在弦长相对位置x/L=0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9 和0.99 处的边界层速度分布.y表示空间点到水翼表面的垂直距离，将y除以水翼尾部厚度h2进行量纲一化.Uxmean表示流场平均速度，Uref表示自由剪切流速度，将Uxmean/Uref进行流向速度量纲一化.将x/L=0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9 和0.99 处的速度分布分别向后移动0.5、1.0、1.5、2.0、2.5、3.0、3.5、4.0 和4.5 个单位.

图11 不同流速下边界层速度分布Fig.11 Velocity distribution in boundary layer at different velocities

在水翼表面速度为0，离水翼表面越远，边界层速度越大直到达到湍流核心区的流速.图11（a）为20 m/s来流速度下的边界层速度分布，模拟结果与实验结果吻合较好.空间上某点速度达到0.99Uref时，将该点到水翼表面的垂直距离定义为边界层厚度，图12为不同流速下的边界层厚度.

图12 不同流速下边界层厚度Fig.12 Thickness of boundary layer at different velocities

边界层厚度沿来流方向逐渐增大，在99%弦长位置边界层厚度达到最大.随着流速增大，边界层厚度越小，这种现象在与文献[24]一致.但在90%～99%弦长范围内，锁频区12 m/s 流速下边界层厚度小于非锁频区20 m/s 流速下边界层厚度，这可能是由于锁频区水翼振幅较大，最终导致边界层厚度减小.

4.2.2 水翼尾迹区流场分析

水翼尾迹区不同时刻瞬时速度场如图13 和图14 所示，图13(a)、图14(a)均为实验结果[4]，图13(b)、图14(b)为模拟结果.

图13 非锁频区20 m/s 流速水翼尾迹区不同时刻瞬时速度场实验、模拟对比Fig.13 Typical instantaneous velocity fields in the hydrofoil wake for lock-off conditions at different moments (20 m/s)

图14 锁频区12 m/s 流速水翼尾迹区不同时刻瞬时速度场实验、模拟对比Fig.14 Typical instantaneous velocity fields in the hydrofoil wake for lock-in condition at different moments (12 m/s)

在非锁频区20 m/s 流速工况下，尾流结构表现出三维涡脱落模式，脱落过程沿水翼展向上不同步，速度场显示为不规则平行涡脱落，如图13 所示.对比实验结果和模拟结果，实验和模拟中都出现了较为明显的涡旋破裂.

在锁频区12 m/s 流速工况下，水翼扭转模态被激发，水翼尾部的横向运动显著增强了脱落涡展向结构.速度场显示为平行涡脱落，与非锁频工况相比，锁频工况下水翼尾迹区涡脱落呈展向同相位平行脱落，涡街整体强度更大，如图13 所示.在实验和模拟中，涡结构破裂都较少.

4.2.3 流致振动分析

图15（a）为实验中非锁频区20 m/s 流速下完整振动响应时域图.通过记录结构场网格变形获取水翼振动响应，水翼以拍频形式振动，最大振幅约为10 mm/s，与实验偏差较小，见图15（b）.对振动响应信号进行行快速傅里叶变换，见图15（c），可知振动响应主要包括2 个频率成分，分别是水翼第一阶扭转模态1 023 Hz 和旋涡脱落频率1 464 Hz.由于水翼第一阶弯曲模态对水翼振动响应影响较小，在模拟结果中未能识别出水翼第一阶弯曲模态.旋涡脱落频率和第一节扭转模态固有频率与实验结果对比，偏差分别为14.94%和1.12%.

图15 非锁频区20 m/s 流速振动响应时域图和频域图Fig.15 Time signal and spectrum of the vortex-induced vibration signal for lock-off condition (20 m/s)

图16（a）为实验中锁频区12 m/s 流速下完整振 动响应时域图.在锁频区12 m/s 流速下，水翼以第一阶扭转模态形式振动，见图16（b）.对振动响应信号进行快速傅里叶变换，见图16（c），可知振动响应包括一个频率成分，为水翼第一阶扭转模态固有频率878 Hz.

图16 锁频区12 m/s 流速振动响应时域图和频域图Fig.16 Time signal and spectrum of the vortex-induced vibration signal for lock-in condition (12 m/s)

在锁频区12 m/s 流速情况下，水翼振幅约为160 mm/s，与实验结果600 mm/s 相比有偏差，但振动响应频率与实验偏差仅为1.34%.曾永顺[25]在实验中得到在锁频区存在负水力阻尼现象，振动结构从水体中吸收能量，导致结构振幅增大.对于水翼振幅，锁频区结果远大于非锁频区，符合实验趋势.表明采用分离式双向流固耦合方法已能预测锁频与非锁频的振幅相对差异.关于锁频区绝对振动幅值预测不精确问题，需要从数值阻尼、动网格、迭代策略等方面开展进一步研究.

5 结 论

基于嵌套网格方法，对水力机械领域常用典型对称水翼进行了双向流固耦合计算，主要对锁频区12 m/s 和非锁频区20 m/s 两个典型流速进行了数值计算.分析了嵌套网格的适用性、水翼的固有频率和振型、以及水翼的振动特性和流场特性，主要结论如下：

①基于声学单元的湿模态分析方法，可较准确得到水下的水翼模态参数，与实验结果相比，前三阶预测的固有频率偏差均在4%以内.

②分析了嵌套区域范围、背景网格和组件网格尺寸差异对计算结果的影响.当嵌套区域越大时，旋涡脱落频率逐渐趋于稳定且接近于实验结果，与实验结果相对偏差为1.18%；两套网格在交界面处应保持网格尺寸一致.

③对水翼进行了非锁频区和锁频区的双向流固耦合计算，对比了振动特性和流场特性.与非锁频工况相比，锁频工况下水翼尾部边界层变薄，尾迹区涡脱落呈平行脱落，涡结构破坏较少；在锁频区，水翼以扭转模态振动，振动频率为878 Hz，与实验相对偏差为1.34%，振动幅值较非锁频工况明显增大.