非對稱尾部形狀水翼水力阻尼識別方法研究

曾永順，姚志峰，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）平衡點數量過多導致多項式擬合階次過高，微小平衡點位置變化會很大程度地影響擬合函數，因此水力阻尼比識別精度不高。工程上，若可接受一定的誤差，為了節省數據處理時間可采用線性插值法。若對水力阻尼比精度有更高的要求，則需采用光滑樣條曲線法。