曾永順，姚志峰，2，楊正軍，王福軍，2

（1.中國農(nóng)業(yè)大學(xué) 水利與土木工程學(xué)院，北京 100083；2.北京市供水管網(wǎng)系統(tǒng)安全與節(jié)能工程技術(shù)研究中心，北京 100083；3.中交天津航道局有限公司，天津 300461；4.天津市疏浚工程技術(shù)企業(yè)重點(diǎn)實(shí)驗(yàn)室，天津 300461）

1 研究背景

由于風(fēng)能、太陽能等新能源調(diào)蓄的要求，抽水蓄能水利工程越來越多。水泵水輪機(jī)作為核心水力機(jī)械，通常需要頻繁切換工況，尤其是在非設(shè)計(jì)工況運(yùn)行，增加了葉輪承受高幅值動(dòng)態(tài)結(jié)構(gòu)應(yīng)力的危險(xiǎn)［1］。這就需要在設(shè)計(jì)階段準(zhǔn)確預(yù)測(cè)葉輪等結(jié)構(gòu)疲勞壽命。水力阻尼是預(yù)測(cè)結(jié)構(gòu)疲勞壽命的關(guān)鍵參數(shù)之一，直接決定了流激振動(dòng)幅值精度［2］。水翼作為水力機(jī)械葉輪和導(dǎo)葉的核心工作單元，當(dāng)流體繞過時(shí)，會(huì)在尾部形成交替脫落的卡門渦街［3-4］。對(duì)于非對(duì)稱尾部形狀水翼，不僅受周期性卡門渦街的流激振動(dòng)影響，同時(shí)受到較大升力作用［5-6］。這種復(fù)雜非定常外力作用，導(dǎo)致水翼動(dòng)態(tài)響應(yīng)信號(hào)的平衡位置具有時(shí)變特性，采用傳統(tǒng)的自由振動(dòng)衰減法將大幅度增加水力阻尼比的識(shí)別誤差，在很多工況下，該方法甚至無法使用。

對(duì)于水力阻尼參數(shù)的識(shí)別，近幾年受到國內(nèi)外學(xué)者越來越多的關(guān)注。加拿大蒙特利爾大學(xué)Liaghat［7］等通過雙向流固耦合數(shù)值計(jì)算發(fā)現(xiàn)對(duì)于對(duì)稱尾部形狀水翼，在較低流速下通過頻譜分析和自由振動(dòng)衰減法，獲得了水翼結(jié)構(gòu)振動(dòng)固有頻率和旋渦脫落頻率，后者振幅約為前者振幅的15%，所識(shí)別的水力阻尼比與實(shí)驗(yàn)值誤差達(dá)20%?；趯?shí)驗(yàn)方法，瑞士洛桑聯(lián)邦理工大學(xué)Roth 等［8］應(yīng)用巴特沃斯低通濾波器和帶通濾波器處理振動(dòng)響應(yīng)信號(hào)，分別獲得了彎曲模態(tài)和扭轉(zhuǎn)模態(tài)的振動(dòng)響應(yīng)，表明濾波處理可有效消除渦激振動(dòng)干擾。挪威科技大學(xué)Bergan［9］和中國農(nóng)業(yè)大學(xué)姚志峰［10］等分別在零攻角條件下開展了高流速水翼流激振動(dòng)測(cè)試，實(shí)驗(yàn)都表明自由振動(dòng)衰減法在非對(duì)稱尾部形狀水翼的場(chǎng)合幾乎無法使用。只能通過共振放大法或半功率（寬帶）法，但這些方法需要在共振區(qū)域以較小的頻率間隔持續(xù)激勵(lì)水翼，實(shí)驗(yàn)難度大，耗費(fèi)時(shí)間長(zhǎng)。

已有研究表明，借助雙向流固耦合數(shù)值計(jì)算方法可模擬水翼結(jié)構(gòu)受短時(shí)沖擊后的動(dòng)力響應(yīng)信號(hào)［6，11］。為了消除平衡點(diǎn)位置變化對(duì)自由振動(dòng)衰減信號(hào)的影響，需要首先擬合平衡點(diǎn)位置的變化曲線，并依據(jù)該曲線做平衡位置的修正。常用的曲線擬合方法包括線性擬合法、多項(xiàng)式擬合法和樣條曲線擬合法等。文獻(xiàn)［12］基于諧波信號(hào)對(duì)插值算法的測(cè)量不確定度進(jìn)行理論分析，通過仿真實(shí)驗(yàn)證明了提高采樣頻率并對(duì)噪聲進(jìn)行濾波處理可提高測(cè)量精度，同時(shí)驗(yàn)證了線性插值方法在諧波分析中的可靠性；分段最小二乘的多項(xiàng)式擬合方法在水泵馬鞍形特性曲線擬合方面具有一定的優(yōu)勢(shì)，能夠保證擬合的連續(xù)性且精度高［13］；樣條曲線能夠?qū)崿F(xiàn)在控制點(diǎn)處一階和二階導(dǎo)數(shù)連續(xù)［14］，采用該方法擬合所有平衡點(diǎn)，可得到光滑擬合曲線且曲率過渡均勻。

本文將基于雙向流固耦合數(shù)值模擬方法，得到NACA 0009 對(duì)稱和非對(duì)稱尾部形狀水翼在動(dòng)水中的振動(dòng)響應(yīng)，采用帶通濾波處理振動(dòng)信號(hào)，消除旋渦脫落頻率對(duì)振動(dòng)響應(yīng)的影響。針對(duì)非對(duì)稱尾部形狀水翼，分析升力作用對(duì)振動(dòng)信號(hào)的干擾，比較三種不同振動(dòng)響應(yīng)校準(zhǔn)方法對(duì)水力阻尼比識(shí)別精度的影響，并給出工程建議。

2 理論背景

2.1 流激振動(dòng)計(jì)算理論基于Ansys workbench 平臺(tái)的雙向流固耦合數(shù)值模擬方法，流場(chǎng)和結(jié)構(gòu)場(chǎng)單獨(dú)求解，通過耦合交界面?zhèn)鬟f數(shù)據(jù)。流動(dòng)計(jì)算基于Navier-Stokes 方程，采用雷諾時(shí)均方法（RANS），湍流模型采用轉(zhuǎn)捩SST 模型，由SST k-ω 模型耦合γ-Reθt轉(zhuǎn)捩模型形成，可在近壁區(qū)完全求解，理論上對(duì)轉(zhuǎn)捩點(diǎn)位置及轉(zhuǎn)捩區(qū)長(zhǎng)度預(yù)測(cè)更加精確［15-16］，具體公式見文獻(xiàn)［17］。

結(jié)構(gòu)計(jì)算基于結(jié)構(gòu)動(dòng)力學(xué)方程：

2.2 水力阻尼表征對(duì)于單自由度系統(tǒng)，振動(dòng)的時(shí)間響應(yīng)函數(shù)可以寫成：

2.3 數(shù)字信號(hào)處理方法利用Matlab 軟件，對(duì)水下結(jié)構(gòu)振動(dòng)響應(yīng)時(shí)域信號(hào)，采用傅里葉變換得到頻域信息，采用巴特沃斯帶通數(shù)值濾波器進(jìn)行濾波處理，對(duì)濾波后數(shù)據(jù)的平衡點(diǎn)進(jìn)行擬合修正。

離散數(shù)字信號(hào)的傅里葉變換可表示為［19］：

式中： xn為時(shí)域輸入序列； Xk為頻域輸出序列； N 為正整數(shù)。采用快速傅里葉變換（FFT）實(shí)現(xiàn)時(shí)頻域轉(zhuǎn)換。巴特沃斯數(shù)字濾波器對(duì)振動(dòng)信號(hào)的處理方式可表示為：

式中： |H (jω)|2=A(ω2)為振幅平方； jω 為頻率（Hz）； ωc為截止頻率（Hz）；M 為濾波器的階數(shù)。使用Matlab 中巴特沃斯帶通數(shù)字濾波函數(shù)實(shí)現(xiàn)數(shù)字信號(hào)濾波。具體公式見文獻(xiàn)［18］。

3 計(jì)算設(shè)置

研究對(duì)象為NACA 0009 水翼，計(jì)算域?yàn)?50 mm×150 mm×750 mm 的矩形通道，水翼在通道中攻角為0°，具體尺寸如圖1所示。對(duì)于對(duì)稱尾部形狀水翼，弦長(zhǎng)L=100 mm，水翼尾部厚度h=3.22 mm，展向?qū)挾葁=150 mm。對(duì)稱尾部形狀水翼材料為不銹鋼，密度為 ρ =7700 kg/m3，彈性模量為E=215 GPa，泊松比為υ =0.3。

圖1 計(jì)算模型（單位：mm）

進(jìn)口邊界條件采用速度進(jìn)口，v=15 m/s，出口邊界條件為平均靜壓，P=2.5Bar，流體與水翼相接觸的面設(shè)置為交界面，遠(yuǎn)離支座的面采用對(duì)稱面，其余面為無滑移壁面。對(duì)于結(jié)構(gòu)場(chǎng)，采用有限元法進(jìn)行瞬態(tài)動(dòng)力學(xué)分析。對(duì)于流場(chǎng)采用有限體積法，瞬態(tài)項(xiàng)采用二階歐拉后差分格式。采用雙向流固耦合數(shù)值模擬方法對(duì)水翼進(jìn)行水力阻尼特性研究，以流場(chǎng)非定常計(jì)算結(jié)果為流固耦合計(jì)算的流場(chǎng)初始文件。流固耦合迭代收斂標(biāo)準(zhǔn)為1×10-4，迭代步數(shù)為30。計(jì)算50 個(gè)時(shí)間步長(zhǎng)后，沿y 軸正方向在水翼表面中線上施加一個(gè)激勵(lì)使其產(chǎn)生彎曲變形。為了避免激勵(lì)太大超過許用應(yīng)力，或者激勵(lì)太小導(dǎo)致激勵(lì)引起的振動(dòng)淹沒在渦激振動(dòng)中，激勵(lì)大小選擇為200 N。該激勵(lì)方法與Zeng 等［20］的計(jì)算設(shè)置相一致。通過記錄流場(chǎng)動(dòng)網(wǎng)格變形獲取水翼的振動(dòng)響應(yīng)，記錄點(diǎn)具體位置為水翼尾緣的幾何中心，如圖1所示。

流場(chǎng)及結(jié)構(gòu)場(chǎng)網(wǎng)格如圖2所示，網(wǎng)格類型均采用正六面體結(jié)構(gòu)化網(wǎng)格，流場(chǎng)水翼表面y+控制在1左右。圖中區(qū)域A、B 和C 分別代表水翼的前緣、近壁區(qū)和尾部，水翼表面第一層網(wǎng)格厚度僅有2×10-6m。為保證非對(duì)稱尾部形狀水翼的網(wǎng)格質(zhì)量，在尖角處網(wǎng)格局部加密。

如圖3所示，在10 m/s 流速下，以旋渦脫落頻率為關(guān)鍵參數(shù)進(jìn)行時(shí)間步長(zhǎng)無關(guān)性驗(yàn)證。結(jié)果表明當(dāng)時(shí)間步長(zhǎng)Δt 從2×10-5s 下降到1×10-5s 時(shí)，頻率相對(duì)變化量在1%以內(nèi)?？烧J(rèn)為Δt=2×10-5s 時(shí)，時(shí)間步長(zhǎng)對(duì)計(jì)算結(jié)果沒有影響。將該時(shí)間步長(zhǎng)下得到模擬結(jié)果與實(shí)驗(yàn)結(jié)果［4］比較，誤差為2.67%。

圖2 網(wǎng)格

圖3 時(shí)間步長(zhǎng)無關(guān)性檢查

4 結(jié)果與討論

4.1 對(duì)稱尾部形狀水翼水力阻尼比識(shí)別流場(chǎng)尾跡區(qū)速度分布如圖4所示，由速度流線可以看出在水翼尾部有旋渦產(chǎn)生；由速度云圖可以看出旋渦交替對(duì)稱分布，脫落渦引起的渦激振動(dòng)將會(huì)對(duì)振動(dòng)響應(yīng)造成干擾；由局部放大的速度矢量圖可以看出近壁區(qū)速度分布情況，流動(dòng)速度從水翼表面到主流區(qū)逐漸增大，這一變化規(guī)律與真實(shí)流動(dòng)情況相一致。

如圖5（a）所示，在v=15 m/s 流速下通過記錄動(dòng)網(wǎng)格變形得到水翼的振動(dòng)響應(yīng)。如圖5（b）所示，對(duì)原始振動(dòng)信號(hào)進(jìn)行快速傅里葉變換，可知振動(dòng)響應(yīng)包括三個(gè)頻率成分。這三個(gè)頻率成分分別為第一階模態(tài)固有頻率195.31 Hz、第二階模態(tài)固有頻率854.49 Hz 和旋渦脫落頻率1074.21 Hz。

圖4 對(duì)稱尾部形狀水翼尾跡區(qū)速度分布（v=15m/s）

圖5 濾波前振動(dòng)響應(yīng)時(shí)域圖及頻域圖（v=15m/s）

如圖6所示，水中第一階模態(tài)和第二階模態(tài)的振型分別為彎曲變形和扭轉(zhuǎn)變形，圖5（b）中振動(dòng)響應(yīng)中的主頻為水中彎曲模態(tài)對(duì)應(yīng)的頻率。瑞士洛桑聯(lián)邦理工大學(xué)高速水洞中得到的固有頻率結(jié)果［8］和旋渦脫落頻率結(jié)果［4］，在本文中作為數(shù)值模擬方法可靠性的驗(yàn)證。由表1可知，雙向流固耦合模擬結(jié)果與實(shí)驗(yàn)結(jié)果很好地吻合，第一階模態(tài)固有頻率、第二階模態(tài)固有頻率和旋渦脫落頻率與實(shí)驗(yàn)結(jié)果偏差分別為0.68%、2.90%和1.12%。

圖6 對(duì)稱尾部形狀水翼低階模態(tài)對(duì)應(yīng)的振型

采用巴特沃斯帶通數(shù)值濾波器對(duì)振動(dòng)信號(hào)進(jìn)行濾波處理，通帶頻率為100～200 Hz，圖7（a）為濾波處理后的振動(dòng)響應(yīng)。將濾波后的振動(dòng)響應(yīng)經(jīng)過快速傅里葉變換后，得到振動(dòng)響應(yīng)頻域圖。如圖7（b）所示，濾波處理后頻率成分僅為第一階模態(tài)固有頻率。

表1 對(duì)稱尾部形狀水翼固有頻率及旋渦脫落頻率模擬結(jié)果與實(shí)驗(yàn)結(jié)果比較

圖7 濾波前后振動(dòng)響應(yīng)時(shí)域圖及頻域圖（v=15m/s）

根據(jù)式（2），采用自由振動(dòng)衰減法進(jìn)行函數(shù)擬合。擬合算法為Trust-Region，其中，DiffMin-Change 參數(shù)設(shè)置為1×10-20，DiffMaxChange 參數(shù)設(shè)置為1×10-12，MaxFunEvals 參數(shù)設(shè)置為6×107，MaxIter 參數(shù)設(shè)置為4×107，TolFun 參數(shù)設(shè)置為1×10-30，TolX 參數(shù)設(shè)置為1×10-30。函數(shù)擬合后，得到的R2及Adjusted R2均在0.99 以上。

對(duì)稱尾部形狀水翼在零攻角條件下不受升力作用，且激勵(lì)水翼的時(shí)間小于1/4 個(gè)振動(dòng)周期，可認(rèn)為在第一個(gè)峰值點(diǎn)后的振動(dòng)響應(yīng)不受激勵(lì)的影響。經(jīng)過濾波處理后，前1/2 個(gè)周期的振動(dòng)響應(yīng)有較大變化，故在水力阻尼識(shí)別時(shí)舍棄第一個(gè)振動(dòng)響應(yīng)的峰值點(diǎn)。圖8為振動(dòng)信號(hào)濾波處理前后，上下峰值點(diǎn)擬合的水力阻尼比。不受外力作用下，理論上基于上下峰值點(diǎn)識(shí)別的水力阻尼比一致。

圖8 水力阻尼比識(shí)別

將基于上下峰值點(diǎn)識(shí)別的水力阻尼比之間的相對(duì)偏差，定義為水力阻尼比識(shí)別偏差度，偏差度越小，說明識(shí)別方法的精度越高。如表2所示，對(duì)于雙向流固耦合數(shù)值模擬結(jié)果，濾波前后水力阻尼比識(shí)別的偏差度分別為7.51%和1.92%。濾波處理后，基本消除了渦激振動(dòng)對(duì)振動(dòng)響應(yīng)信號(hào)的干擾。將濾波處理后的模擬結(jié)果與實(shí)驗(yàn)結(jié)果［8］相比較，最大相對(duì)偏差為8.44%。驗(yàn)證了在識(shí)別水力阻尼比時(shí)，自由振動(dòng)衰減法的可靠性。

4.2 非對(duì)稱尾部形狀水翼水力阻尼比識(shí)別對(duì)于非對(duì)稱尾部形狀水翼，進(jìn)行雙向流流固耦合數(shù)值模擬時(shí)，計(jì)算設(shè)置與對(duì)稱尾部形狀水翼相一致。流場(chǎng)尾跡壓力分布如圖9（a）所示，水翼上表面壓力明顯高于下表面，且由速度流線可以看出在水翼尾部有旋渦脫落，此時(shí)水翼振動(dòng)響應(yīng)受到渦激振動(dòng)和升力作用共同干擾。

表2 濾波前后基于上下峰值點(diǎn)識(shí)別的水力阻尼比

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

流場(chǎng)非定常計(jì)算結(jié)果表明，出現(xiàn)交替脫落的卡門渦后，升力出現(xiàn)明顯的周期性變化。圖9（b）為15 m/s 流速下，非對(duì)稱尾部形狀水翼受到的升力作用隨時(shí)間的變化，升力平均值Clave=69.81 N，升力方向指向y 軸負(fù)方向。

非對(duì)稱尾部形狀水翼低階模態(tài)固振型如圖10所示，尾部修型后第一階和第二階模態(tài)振型分別為彎曲變形和扭轉(zhuǎn)變形，振型與對(duì)稱尾部形狀水翼相一致。

圖10 非對(duì)稱尾部形狀水翼低階模態(tài)對(duì)應(yīng)的振型

圖11（a）為15 m/s 流速下非對(duì)稱尾部形狀水翼的振動(dòng)響應(yīng)時(shí)域圖，在渦激振動(dòng)作用下，振動(dòng)響應(yīng)出現(xiàn)上下波動(dòng)，在升力作用下，平衡位置移向y 軸負(fù)方向。由圖11（b）可知振動(dòng)響應(yīng)頻率成分只包括第一階模態(tài)固有頻率146.11 Hz 和旋渦脫落頻率964.34 Hz，與實(shí)驗(yàn)結(jié)果［10］比較相對(duì)誤差分別為10.89%和1.42%。圖11（c）采用巴特沃斯帶通數(shù)值濾波器進(jìn)行濾波，通帶頻率為100～200 Hz。圖11（d）為濾波后的頻域圖，振動(dòng)響應(yīng)頻率成分只保留第一階模態(tài)固有頻率。

對(duì)于周期恒定簡(jiǎn)諧運(yùn)動(dòng)，從平衡位置分別到波峰和波谷經(jīng)歷的時(shí)間一致。若將升力作用下的振動(dòng)響應(yīng)波峰和波谷找出，將兩者中間時(shí)刻所在的位置定義為振動(dòng)響應(yīng)平衡點(diǎn)。用擬合曲線將這些點(diǎn)連接，再用濾波后的振動(dòng)響應(yīng)減去該曲線，可對(duì)原始振動(dòng)信號(hào)進(jìn)行校準(zhǔn)。如圖12所示，平衡點(diǎn)的連接方式可以是多項(xiàng)式擬合法、線性插值法或者光滑樣條曲線法。

圖11 渦激振動(dòng)和升力共同作用下的振動(dòng)響應(yīng)（v=15m/s）

對(duì)于多項(xiàng)式擬合法，采用N-1 次多項(xiàng)式擬合曲線可過所有的平衡點(diǎn)，其中N 為平衡點(diǎn)個(gè)數(shù)，可直接得到擬合函數(shù)。對(duì)于線性插值法，采用分段線性函數(shù)將所有平衡點(diǎn)連接，需求出N-1 個(gè)線性函數(shù)。對(duì)于光滑樣條曲線法，只能獲得時(shí)間序列所對(duì)應(yīng)的函數(shù)值，無法實(shí)現(xiàn)擬合函數(shù)表達(dá)。將振動(dòng)響應(yīng)原始信號(hào)曲線分別減去三種擬合曲線后，可得校準(zhǔn)后的振動(dòng)響應(yīng)，校準(zhǔn)后的固有頻率都為 fn=141.24 Hz。校準(zhǔn)前固有頻率為146.11 Hz，相對(duì)變化量為3.33%。將校準(zhǔn)后的第一階模態(tài)固有頻率與實(shí)驗(yàn)值［10］比較，相對(duì)偏差為7.58%。

如圖13所示，采用自由振動(dòng)衰減法識(shí)別振動(dòng)響應(yīng)上下峰值點(diǎn)的水力阻尼比?；谌N振動(dòng)響應(yīng)校準(zhǔn)方法的水力阻尼比識(shí)別結(jié)果如表3所示。

由表3可以看出，在升力作用下，不同水力阻尼比識(shí)別方法得到的結(jié)果從0.046 89 到0.1441，需要消除識(shí)別方法對(duì)水力阻尼比的影響。對(duì)比三種水力阻尼比識(shí)別方法的精度，多項(xiàng)式擬合法、線性插值法和光滑樣條曲線法水力阻尼比識(shí)別的偏差度分別為34.93%、3.53%和0.16%。對(duì)于多項(xiàng)式擬合法，偏差度較大。其原因是由于多項(xiàng)式擬合法完全擬合所有的數(shù)據(jù)點(diǎn)形成高階多向式（6 階），若非定常計(jì)算得到平衡點(diǎn)出現(xiàn)微小變化，將會(huì)引起高階多項(xiàng)式的很大改變。線性插值法和光滑樣條曲線法都有較高的識(shí)別精度，與實(shí)驗(yàn)值比較最大誤差分別為7.62%和4.48%。若在工程中需要極其精準(zhǔn)的水力阻尼參數(shù)，可采用光滑樣條曲線法識(shí)別水力阻尼比。若在工程中可接受一定的誤差，為了節(jié)省時(shí)間可采用線性插值法。

5 結(jié)論

圖12 振動(dòng)響應(yīng)校準(zhǔn)方式

圖13 水力阻尼比識(shí)別

表3 振動(dòng)響應(yīng)校準(zhǔn)后基于上下峰值點(diǎn)識(shí)別的水力阻尼比

采用雙向流固耦合數(shù)值模擬方法，激勵(lì)水翼后通過記錄動(dòng)網(wǎng)格變形獲取振動(dòng)響應(yīng)。兩種尾緣形狀水翼的第一階模態(tài)固有頻率、第二階模態(tài)固有頻率和旋渦脫落頻率預(yù)測(cè)結(jié)果與實(shí)驗(yàn)吻合良好，在15 m/s 流速下最大偏差分別7.58%、2.90%和1.42%。采用自由振動(dòng)衰減法識(shí)別水力阻尼比，分析了渦激振動(dòng)作用和升力作用下的水力阻尼比識(shí)別方法，主要結(jié)論如下：（1）在渦激振動(dòng)的干擾下，振動(dòng)響應(yīng)包含脫落渦的頻率成分。15 m/s 流速下，帶通濾波前后對(duì)稱尾部形狀水翼的水力阻尼比識(shí)別偏差度分別7.51%和1.92%。濾波處理可顯著提高水力阻尼比的識(shí)別精度。（2）對(duì)于非對(duì)稱尾部形狀水翼，振動(dòng)響應(yīng)受到渦激振動(dòng)和升力作用共同干擾。濾波處理后找到平衡點(diǎn)，采用多項(xiàng)式擬合法、線性插值法和光滑樣條曲線法校準(zhǔn)振動(dòng)響應(yīng)，三種方法水力阻尼比識(shí)別的偏差度分別為34.93%、3.53%和0.16%。（3）平衡點(diǎn)數(shù)量過多導(dǎo)致多項(xiàng)式擬合階次過高，微小平衡點(diǎn)位置變化會(huì)很大程度地影響擬合函數(shù)，因此水力阻尼比識(shí)別精度不高。工程上，若可接受一定的誤差，為了節(jié)省數(shù)據(jù)處理時(shí)間可采用線性插值法。若對(duì)水力阻尼比精度有更高的要求，則需采用光滑樣條曲線法。