非对称尾部形状水翼水力阻尼识别方法研究

曾永顺，姚志峰，2，杨正军，王福军，2

（1.中国农业大学 水利与土木工程学院，北京 100083；2.北京市供水管网系统安全与节能工程技术研究中心，北京 100083；3.中交天津航道局有限公司，天津 300461；4.天津市疏浚工程技术企业重点实验室，天津 300461）

1 研究背景

由于风能、太阳能等新能源调蓄的要求，抽水蓄能水利工程越来越多。水泵水轮机作为核心水力机械，通常需要频繁切换工况，尤其是在非设计工况运行，增加了叶轮承受高幅值动态结构应力的危险［1］。这就需要在设计阶段准确预测叶轮等结构疲劳寿命。水力阻尼是预测结构疲劳寿命的关键参数之一，直接决定了流激振动幅值精度［2］。水翼作为水力机械叶轮和导叶的核心工作单元，当流体绕过时，会在尾部形成交替脱落的卡门涡街［3-4］。对于非对称尾部形状水翼，不仅受周期性卡门涡街的流激振动影响，同时受到较大升力作用［5-6］。这种复杂非定常外力作用，导致水翼动态响应信号的平衡位置具有时变特性，采用传统的自由振动衰减法将大幅度增加水力阻尼比的识别误差，在很多工况下，该方法甚至无法使用。

对于水力阻尼参数的识别，近几年受到国内外学者越来越多的关注。加拿大蒙特利尔大学Liaghat［7］等通过双向流固耦合数值计算发现对于对称尾部形状水翼，在较低流速下通过频谱分析和自由振动衰减法，获得了水翼结构振动固有频率和旋涡脱落频率，后者振幅约为前者振幅的15%，所识别的水力阻尼比与实验值误差达20%。基于实验方法，瑞士洛桑联邦理工大学Roth 等［8］应用巴特沃斯低通滤波器和带通滤波器处理振动响应信号，分别获得了弯曲模态和扭转模态的振动响应，表明滤波处理可有效消除涡激振动干扰。挪威科技大学Bergan［9］和中国农业大学姚志峰［10］等分别在零攻角条件下开展了高流速水翼流激振动测试，实验都表明自由振动衰减法在非对称尾部形状水翼的场合几乎无法使用。只能通过共振放大法或半功率（宽带）法，但这些方法需要在共振区域以较小的频率间隔持续激励水翼，实验难度大，耗费时间长。

已有研究表明，借助双向流固耦合数值计算方法可模拟水翼结构受短时冲击后的动力响应信号［6，11］。为了消除平衡点位置变化对自由振动衰减信号的影响，需要首先拟合平衡点位置的变化曲线，并依据该曲线做平衡位置的修正。常用的曲线拟合方法包括线性拟合法、多项式拟合法和样条曲线拟合法等。文献［12］基于谐波信号对插值算法的测量不确定度进行理论分析，通过仿真实验证明了提高采样频率并对噪声进行滤波处理可提高测量精度，同时验证了线性插值方法在谐波分析中的可靠性；分段最小二乘的多项式拟合方法在水泵马鞍形特性曲线拟合方面具有一定的优势，能够保证拟合的连续性且精度高［13］；样条曲线能够实现在控制点处一阶和二阶导数连续［14］，采用该方法拟合所有平衡点，可得到光滑拟合曲线且曲率过渡均匀。

本文将基于双向流固耦合数值模拟方法，得到NACA 0009 对称和非对称尾部形状水翼在动水中的振动响应，采用带通滤波处理振动信号，消除旋涡脱落频率对振动响应的影响。针对非对称尾部形状水翼，分析升力作用对振动信号的干扰，比较三种不同振动响应校准方法对水力阻尼比识别精度的影响，并给出工程建议。

2 理论背景

2.1 流激振动计算理论基于Ansys workbench 平台的双向流固耦合数值模拟方法，流场和结构场单独求解，通过耦合交界面传递数据。流动计算基于Navier-Stokes 方程，采用雷诺时均方法（RANS），湍流模型采用转捩SST 模型，由SST k-ω 模型耦合γ-Reθt转捩模型形成，可在近壁区完全求解，理论上对转捩点位置及转捩区长度预测更加精确［15-16］，具体公式见文献［17］。

结构计算基于结构动力学方程：

2.2 水力阻尼表征对于单自由度系统，振动的时间响应函数可以写成：

2.3 数字信号处理方法利用Matlab 软件，对水下结构振动响应时域信号，采用傅里叶变换得到频域信息，采用巴特沃斯带通数值滤波器进行滤波处理，对滤波后数据的平衡点进行拟合修正。

离散数字信号的傅里叶变换可表示为［19］：

式中： xn为时域输入序列； Xk为频域输出序列； N 为正整数。采用快速傅里叶变换（FFT）实现时频域转换。巴特沃斯数字滤波器对振动信号的处理方式可表示为：

式中： |H (jω)|2=A(ω2)为振幅平方； jω 为频率（Hz）； ωc为截止频率（Hz）；M 为滤波器的阶数。使用Matlab 中巴特沃斯带通数字滤波函数实现数字信号滤波。具体公式见文献［18］。

3 计算设置

研究对象为NACA 0009 水翼，计算域为150 mm×150 mm×750 mm 的矩形通道，水翼在通道中攻角为0°，具体尺寸如图1所示。对于对称尾部形状水翼，弦长L=100 mm，水翼尾部厚度h=3.22 mm，展向宽度w=150 mm。对称尾部形状水翼材料为不锈钢，密度为 ρ =7700 kg/m3，弹性模量为E=215 GPa，泊松比为υ =0.3。

图1 计算模型（单位：mm）

进口边界条件采用速度进口，v=15 m/s，出口边界条件为平均静压，P=2.5Bar，流体与水翼相接触的面设置为交界面，远离支座的面采用对称面，其余面为无滑移壁面。对于结构场，采用有限元法进行瞬态动力学分析。对于流场采用有限体积法，瞬态项采用二阶欧拉后差分格式。采用双向流固耦合数值模拟方法对水翼进行水力阻尼特性研究，以流场非定常计算结果为流固耦合计算的流场初始文件。流固耦合迭代收敛标准为1×10-4，迭代步数为30。计算50 个时间步长后，沿y 轴正方向在水翼表面中线上施加一个激励使其产生弯曲变形。为了避免激励太大超过许用应力，或者激励太小导致激励引起的振动淹没在涡激振动中，激励大小选择为200 N。该激励方法与Zeng 等［20］的计算设置相一致。通过记录流场动网格变形获取水翼的振动响应，记录点具体位置为水翼尾缘的几何中心，如图1所示。

流场及结构场网格如图2所示，网格类型均采用正六面体结构化网格，流场水翼表面y+控制在1左右。图中区域A、B 和C 分别代表水翼的前缘、近壁区和尾部，水翼表面第一层网格厚度仅有2×10-6m。为保证非对称尾部形状水翼的网格质量，在尖角处网格局部加密。

如图3所示，在10 m/s 流速下，以旋涡脱落频率为关键参数进行时间步长无关性验证。结果表明当时间步长Δt 从2×10-5s 下降到1×10-5s 时，频率相对变化量在1%以内。可认为Δt=2×10-5s 时，时间步长对计算结果没有影响。将该时间步长下得到模拟结果与实验结果［4］比较，误差为2.67%。

图2 网格

图3 时间步长无关性检查

4 结果与讨论

4.1 对称尾部形状水翼水力阻尼比识别流场尾迹区速度分布如图4所示，由速度流线可以看出在水翼尾部有旋涡产生；由速度云图可以看出旋涡交替对称分布，脱落涡引起的涡激振动将会对振动响应造成干扰；由局部放大的速度矢量图可以看出近壁区速度分布情况，流动速度从水翼表面到主流区逐渐增大，这一变化规律与真实流动情况相一致。

如图5（a）所示，在v=15 m/s 流速下通过记录动网格变形得到水翼的振动响应。如图5（b）所示，对原始振动信号进行快速傅里叶变换，可知振动响应包括三个频率成分。这三个频率成分分别为第一阶模态固有频率195.31 Hz、第二阶模态固有频率854.49 Hz 和旋涡脱落频率1074.21 Hz。

图4 对称尾部形状水翼尾迹区速度分布（v=15m/s）

图5 滤波前振动响应时域图及频域图（v=15m/s）

如图6所示，水中第一阶模态和第二阶模态的振型分别为弯曲变形和扭转变形，图5（b）中振动响应中的主频为水中弯曲模态对应的频率。瑞士洛桑联邦理工大学高速水洞中得到的固有频率结果［8］和旋涡脱落频率结果［4］，在本文中作为数值模拟方法可靠性的验证。由表1可知，双向流固耦合模拟结果与实验结果很好地吻合，第一阶模态固有频率、第二阶模态固有频率和旋涡脱落频率与实验结果偏差分别为0.68%、2.90%和1.12%。

图6 对称尾部形状水翼低阶模态对应的振型

采用巴特沃斯带通数值滤波器对振动信号进行滤波处理，通带频率为100～200 Hz，图7（a）为滤波处理后的振动响应。将滤波后的振动响应经过快速傅里叶变换后，得到振动响应频域图。如图7（b）所示，滤波处理后频率成分仅为第一阶模态固有频率。

表1 对称尾部形状水翼固有频率及旋涡脱落频率模拟结果与实验结果比较

图7 滤波前后振动响应时域图及频域图（v=15m/s）

根据式（2），采用自由振动衰减法进行函数拟合。拟合算法为Trust-Region，其中，DiffMin-Change 参数设置为1×10-20，DiffMaxChange 参数设置为1×10-12，MaxFunEvals 参数设置为6×107，MaxIter 参数设置为4×107，TolFun 参数设置为1×10-30，TolX 参数设置为1×10-30。函数拟合后，得到的R2及Adjusted R2均在0.99 以上。

对称尾部形状水翼在零攻角条件下不受升力作用，且激励水翼的时间小于1/4 个振动周期，可认为在第一个峰值点后的振动响应不受激励的影响。经过滤波处理后，前1/2 个周期的振动响应有较大变化，故在水力阻尼识别时舍弃第一个振动响应的峰值点。图8为振动信号滤波处理前后，上下峰值点拟合的水力阻尼比。不受外力作用下，理论上基于上下峰值点识别的水力阻尼比一致。

图8 水力阻尼比识别

将基于上下峰值点识别的水力阻尼比之间的相对偏差，定义为水力阻尼比识别偏差度，偏差度越小，说明识别方法的精度越高。如表2所示，对于双向流固耦合数值模拟结果，滤波前后水力阻尼比识别的偏差度分别为7.51%和1.92%。滤波处理后，基本消除了涡激振动对振动响应信号的干扰。将滤波处理后的模拟结果与实验结果［8］相比较，最大相对偏差为8.44%。验证了在识别水力阻尼比时，自由振动衰减法的可靠性。

4.2 非对称尾部形状水翼水力阻尼比识别对于非对称尾部形状水翼，进行双向流流固耦合数值模拟时，计算设置与对称尾部形状水翼相一致。流场尾迹压力分布如图9（a）所示，水翼上表面压力明显高于下表面，且由速度流线可以看出在水翼尾部有旋涡脱落，此时水翼振动响应受到涡激振动和升力作用共同干扰。

表2 滤波前后基于上下峰值点识别的水力阻尼比

图9 升力作用（v=15m/s）

流场非定常计算结果表明，出现交替脱落的卡门涡后，升力出现明显的周期性变化。图9（b）为15 m/s 流速下，非对称尾部形状水翼受到的升力作用随时间的变化，升力平均值Clave=69.81 N，升力方向指向y 轴负方向。

非对称尾部形状水翼低阶模态固振型如图10所示，尾部修型后第一阶和第二阶模态振型分别为弯曲变形和扭转变形，振型与对称尾部形状水翼相一致。

图10 非对称尾部形状水翼低阶模态对应的振型

图11（a）为15 m/s 流速下非对称尾部形状水翼的振动响应时域图，在涡激振动作用下，振动响应出现上下波动，在升力作用下，平衡位置移向y 轴负方向。由图11（b）可知振动响应频率成分只包括第一阶模态固有频率146.11 Hz 和旋涡脱落频率964.34 Hz，与实验结果［10］比较相对误差分别为10.89%和1.42%。图11（c）采用巴特沃斯带通数值滤波器进行滤波，通带频率为100～200 Hz。图11（d）为滤波后的频域图，振动响应频率成分只保留第一阶模态固有频率。

对于周期恒定简谐运动，从平衡位置分别到波峰和波谷经历的时间一致。若将升力作用下的振动响应波峰和波谷找出，将两者中间时刻所在的位置定义为振动响应平衡点。用拟合曲线将这些点连接，再用滤波后的振动响应减去该曲线，可对原始振动信号进行校准。如图12所示，平衡点的连接方式可以是多项式拟合法、线性插值法或者光滑样条曲线法。

图11 涡激振动和升力共同作用下的振动响应（v=15m/s）

对于多项式拟合法，采用N-1 次多项式拟合曲线可过所有的平衡点，其中N 为平衡点个数，可直接得到拟合函数。对于线性插值法，采用分段线性函数将所有平衡点连接，需求出N-1 个线性函数。对于光滑样条曲线法，只能获得时间序列所对应的函数值，无法实现拟合函数表达。将振动响应原始信号曲线分别减去三种拟合曲线后，可得校准后的振动响应，校准后的固有频率都为 fn=141.24 Hz。校准前固有频率为146.11 Hz，相对变化量为3.33%。将校准后的第一阶模态固有频率与实验值［10］比较，相对偏差为7.58%。

如图13所示，采用自由振动衰减法识别振动响应上下峰值点的水力阻尼比。基于三种振动响应校准方法的水力阻尼比识别结果如表3所示。

由表3可以看出，在升力作用下，不同水力阻尼比识别方法得到的结果从0.046 89 到0.1441，需要消除识别方法对水力阻尼比的影响。对比三种水力阻尼比识别方法的精度，多项式拟合法、线性插值法和光滑样条曲线法水力阻尼比识别的偏差度分别为34.93%、3.53%和0.16%。对于多项式拟合法，偏差度较大。其原因是由于多项式拟合法完全拟合所有的数据点形成高阶多向式（6 阶），若非定常计算得到平衡点出现微小变化，将会引起高阶多项式的很大改变。线性插值法和光滑样条曲线法都有较高的识别精度，与实验值比较最大误差分别为7.62%和4.48%。若在工程中需要极其精准的水力阻尼参数，可采用光滑样条曲线法识别水力阻尼比。若在工程中可接受一定的误差，为了节省时间可采用线性插值法。

5 结论

图12 振动响应校准方式

图13 水力阻尼比识别

表3 振动响应校准后基于上下峰值点识别的水力阻尼比

采用双向流固耦合数值模拟方法，激励水翼后通过记录动网格变形获取振动响应。两种尾缘形状水翼的第一阶模态固有频率、第二阶模态固有频率和旋涡脱落频率预测结果与实验吻合良好，在15 m/s 流速下最大偏差分别7.58%、2.90%和1.42%。采用自由振动衰减法识别水力阻尼比，分析了涡激振动作用和升力作用下的水力阻尼比识别方法，主要结论如下：（1）在涡激振动的干扰下，振动响应包含脱落涡的频率成分。15 m/s 流速下，带通滤波前后对称尾部形状水翼的水力阻尼比识别偏差度分别7.51%和1.92%。滤波处理可显著提高水力阻尼比的识别精度。（2）对于非对称尾部形状水翼，振动响应受到涡激振动和升力作用共同干扰。滤波处理后找到平衡点，采用多项式拟合法、线性插值法和光滑样条曲线法校准振动响应，三种方法水力阻尼比识别的偏差度分别为34.93%、3.53%和0.16%。（3）平衡点数量过多导致多项式拟合阶次过高，微小平衡点位置变化会很大程度地影响拟合函数，因此水力阻尼比识别精度不高。工程上，若可接受一定的误差，为了节省数据处理时间可采用线性插值法。若对水力阻尼比精度有更高的要求，则需采用光滑样条曲线法。