曾永順，齊 鑫，姚志峰，3，王福軍，3

(1.中國農業(yè)大學 水利與土木工程學院，北京 100083；2.北京航天石化技術裝備工程有限公司，北京 100176；3.北京市供水管網(wǎng)系統(tǒng)安全與節(jié)能工程技術研究中心，北京 100083)

1 研究背景

現(xiàn)代抽水蓄能機組的能源轉化效率可達到90%以上，在新能源并入電網(wǎng)、調峰填谷和調頻調相等方面有重要作用[1]。根據(jù)國家能源局發(fā)布的抽水蓄能中長期發(fā)展規(guī)劃，我國抽水蓄能機組的投產(chǎn)總規(guī)模將從2021年的3249萬千瓦提高到2030年的1.2億千瓦。為滿足未來電力系統(tǒng)需求，抽水蓄能機組有向大流量和高揚程發(fā)展的趨勢，但存在流激振動引起的穩(wěn)定性問題。如機組需頻繁轉換工況，此時水力激振的頻率范圍加寬，葉片結構存在共振風險[2]。特別是大流量工況，空化可能引起額外激振，加大結構振動水平的評估難度[3]。

葉片結構的固有頻率和水力阻尼分別是評估共振工況點和共振幅值的關鍵參數(shù)[4-5]?；谀B(tài)分析方法，可在非空化工況預測水翼和葉輪等結構的固有頻率，相對于實驗結果的誤差在10%以內[6-7]。空化發(fā)生后，Torre等[8]和曾永順[9]基于實驗觀測到水翼固有頻率隨空化數(shù)的降低逐漸增大。但傳統(tǒng)模態(tài)分析方法在理論上將流動方程簡化為聲波方程，難以模擬局部壓強降低引發(fā)的空化現(xiàn)象，導致空化時的結構固有頻率預測存在困難。對于水力阻尼，已發(fā)表文獻表明流動環(huán)境對其影響顯著。如水翼在高流速下，水力阻尼隨流速線性增長[10-14]，但其隨葉頂間隙和攻角的增大而逐漸減小[15-16]。在空化時，曾永順[9]基于實驗得到水翼低階模態(tài)的水力阻尼隨空化數(shù)降低而減小的規(guī)律，但目前缺乏數(shù)值模擬工作對其內部機理進行揭示。

水力阻尼數(shù)值計算方法包括非定常渦格法、能量法和雙向流固耦合法等[4，17-19]。非定常渦格法基于勢流理論，有較高的計算效率，但精度有限，如Kohtanen等[17]對平板進行水力阻尼計算時的誤差有58%。能量法假設振型和固有頻率不隨流速改變，通過計算一個振動周期內的能量耗散獲取水力阻尼，在計算時兼顧了精度和效率，與實驗結果相比較誤差在11%以內[6，19]。但傳統(tǒng)能量法只能求解空化周期內某一個瞬時的水力阻尼，難以捕捉非定常空化的演化過程，計算結果存在隨機性。目前雙向流固耦合法理論上最完善，能夠實時進行流場和結構場的數(shù)據(jù)交換，在水翼的水力阻尼計算時也表現(xiàn)出較高精度[19]。特別是在空化時，雙向流固耦合能描述空化形態(tài)實時變化對結構振動的影響。但由于前期缺乏空化時的水力阻尼實驗數(shù)據(jù)，導致雙向流固耦合法的各關鍵計算參數(shù)取值尚不明確。

本文將葉片簡化為NACA 0009鈍形尾緣水翼，對其在空化時的水力阻尼特性開展研究。實現(xiàn)以下目標：(1)確定適用于雙向流固耦合的各關鍵計算參數(shù)；(2)實現(xiàn)空化時的固有頻率和水力阻尼定量計算；(3)揭示空化對水力阻尼影響的內部機理。

2 理論背景

采用迭代式雙向流固耦合法進行水力阻尼數(shù)值模擬，流場采用有限體積法，結構場采用有限元法，通過流固耦合交界面進行數(shù)值傳輸，基于商業(yè)軟件ANSYS Workbench完成計算。

(1)

Chybrid=χ(ρm/ρf)CFBM+1-χ(ρm/ρf)CDCM

(2)

(3)

式中：μt為渦黏系數(shù)；ρm為汽液混合密度；k為湍動能；ε為湍流耗散率；ρf為流體密度；CFBM為濾波函數(shù)；CDCM為密度修正函數(shù)；系數(shù)C2=0.2，Cμ=0.09。

2.2 結構振動方程水翼受脈沖激勵后，多階模態(tài)同時被激發(fā)，其振動方程為[9]：

(4)

為對總阻尼進行識別，采用巴特沃斯帶通濾波器將單階模態(tài)的振動響應提?。?/p>

y(t)=Ae-ωnζtsin(ωdt+φ)

(5)

假設總阻尼比是其他類型阻尼比的線性疊加，則水力阻尼比為：

ζh=ζt-ζa

(6)

式中ζa為空氣中阻尼比。

2.3 交界面數(shù)據(jù)傳輸方程交界面同時滿足運動學和動力學條件：

選取我院選出接受心臟起搏器植入術患者26例，隨機將患者分成對照組和觀察組，各13例，本次研究是經(jīng)過醫(yī)院倫理委員會同意后實施的，所有患者的精神狀態(tài)良好，認知功能正常。其中對照組男6例，女7例，年齡53～86歲，平均（65.1±4.3）歲，觀察組男7例，女6例，年齡51～85歲，平均（66.2±4.1）歲，通過對兩組患者一般資料進行比較后，差異無統(tǒng)計學意義（P＞0.05）。

(7)

(8)

2.4 空化參數(shù)定義采用無量綱空化數(shù)表征空化程度[20]：

σ=(p-pv)/(0.5ρfv2)

(9)

式中：σ為空化數(shù)；pv為飽和蒸氣壓；v為來流速度。

3 計算設置

3.1 物理模型以NACA 0009鈍尾緣水翼為研究對象，在方形計算域中進行雙向流固耦合數(shù)值模擬，見圖1。水翼弦長、翼展和尾緣厚度分別為70、67和2.25 mm，計算域長、寬和高分別為700、70和190 mm。水翼在計算域中攻角為10°，一端約束，另一端自由，葉頂間隙3 mm。水翼材料為鋁合金，密度、彈性模量和泊松比分別為2700 kg/m3、69 GPa和0.334。

圖1 計算模型

3.2 邊界條件采用速度進口(v=6.64 m/s)，壓力出口，其他邊界為無滑移壁面?？栈瘮?shù)1.04～2.02，根據(jù)出口壓力計算得到。根據(jù)前期文獻調研[21]，水中振動結構的低階模態(tài)更容易被激發(fā)。因此，本文以水翼第一階彎曲模態(tài)為例進行水力阻尼特性研究。為激發(fā)該模態(tài)，在水翼表面施加脈沖激振，見圖1(b)。

4 數(shù)值計算可靠性分析

4.1 計算參數(shù)可靠性分析流場和結構場都采用結構化網(wǎng)格進行劃分，展向投影見圖2。對于流場，水翼表面第一層厚度為0.01 mm，網(wǎng)格拓展比為1.1。由于流固耦合計算通過交界面進行數(shù)據(jù)傳輸，流場和結構場的網(wǎng)格節(jié)點無需嚴格對應。因此，為提高計算效率，結構場網(wǎng)格數(shù)顯著低于流場。根據(jù)文獻[22]對相同水翼的網(wǎng)格收斂性分析，當流場和結構場的網(wǎng)格單元數(shù)分別大于247萬和1.4萬時，對流場繞流特性和結構場振動特性沒有影響。本文也采用相同網(wǎng)格單元數(shù)，以網(wǎng)格質量著色，網(wǎng)格質量在0.4以上，見圖2。

圖2 計算域網(wǎng)格劃分

時間步長影響雙向流固耦合計算效率，在保證精度的前提下其值越大越適合工程應用。對于簡單平板和水翼，Liaghat等[19]認為時間步長需小于1/33個振動周期。曾永順[9]的研究表明一個振動周期內的時間步數(shù)大于25時，水力阻尼的變化量在3%。本文參考前期工作，一個振動周期內的時間步數(shù)取值45，絕對時間步長為1×10-4s。

數(shù)值阻尼在結構瞬態(tài)動力學計算時引入，穩(wěn)定了時間積分格式并抑制了高頻虛假數(shù)值噪聲。其在商業(yè)軟件的默認值為0.1，但針對具體物理模型、模態(tài)階次和時間步長，該值的適用性還需驗證。為消除數(shù)值阻尼對水力阻尼計算結果的影響，根據(jù)試錯法進行數(shù)值阻尼等效。分三步：

第一步，在空氣中進行結構瞬態(tài)動力學計算，根據(jù)瑞利阻尼矩陣輸入阻尼，可表示為[23]：

[C]=a0[M]+a1[K]

(10)

(11)

式中：[C]、[M]和[K]分別為阻尼矩陣、質量矩陣和剛度矩陣；a0和a1分別為質量和剛度矩陣系數(shù)，工程中一般忽略質量矩陣系數(shù)，即a0=0。

假設結構振動是多階模態(tài)的疊加，則質量矩陣和剛度矩陣可分別定義為多階模態(tài)質量和模態(tài)剛度構成的對角矩陣。若在水中進行結構動力學計算，則需在式(10)中考慮水體附加的質量矩陣、剛度矩陣和阻尼矩陣。

第二步，通過調節(jié)數(shù)值阻尼，并基于數(shù)值計算得到不同的輸出阻尼。

第三步，當輸入阻尼與輸出阻尼一致時，認為數(shù)值阻尼取值可靠。

在空氣中以水翼第一階彎曲模態(tài)為例進行數(shù)值阻尼等效，調節(jié)a1令輸入阻尼為0.002。當數(shù)值阻尼為0.3時，結構瞬態(tài)動力學計算得到的阻尼比為0.0057，顯著大于輸入阻尼，見圖3(a)。說明數(shù)值阻尼的存在會導致水力阻尼的過預測。在不同數(shù)值阻尼下，輸入阻尼與輸出阻尼的相對誤差見圖3(b)。當數(shù)值阻尼為默認值(0.1)時，相對誤差可達42%。當數(shù)值阻尼取0.01時，輸入阻尼與輸出阻尼的相對誤差在4%以內，可認為消除了數(shù)值阻尼的影響。當數(shù)值阻尼進一步減小到0.001時，振動響應發(fā)散。

圖3 數(shù)值阻尼可靠性驗證

4.2 計算結果可靠性分析基于模態(tài)分析得到水翼第一階彎曲模態(tài)振型，見圖4，其在空氣中和靜水中一致。但在水體附加質量的作用下，水中固有頻率(230.1 Hz)顯著低于空氣中(405.3 Hz)。根據(jù)實驗測量[9]，該水翼在空氣中和靜水中的第一階彎曲模態(tài)固有頻率分別為425.7 Hz和214.8 Hz。將固有頻率模態(tài)分析結果與實驗結果[9]相比較，誤差在7.12%以內，驗證了計算設置和計算參數(shù)的可靠性。

圖4 水翼第一階彎曲模態(tài)振型

為驗證空化時流場非定常計算的可靠性，定量獲取前緣空化的周期性脫落頻率，并將其與實驗結果相比較，見表1。當空化數(shù)σ在1.04到2.02之間時，模擬結果與實驗結果[9]吻合良好，相對偏差在4.98%到8.76%之間。

表1 空化脫落頻率實驗結果[9]與模擬分析對比

5 結果與討論

5.1 結構振動特性典型空化條件(σ=1.39)的結構振動見圖5(a)。在預設脈沖激振作用下，水翼第一階彎曲模態(tài)被激發(fā)。對比圖4和圖5(a)，位移都是從固定端到自由端逐漸增大，說明空化本質上未改變水翼振動方式?；趧泳W(wǎng)格技術，結構場將位移傳遞給流場，而經(jīng)過流場非定常計算后再將力傳遞回結構場，形成循環(huán)。通過監(jiān)測流場的網(wǎng)格變形可等效獲取結構的振動特性，測點見圖5(b)。

圖5 雙向流固耦合獲取的水翼振動和流場變形

σ=1.39時，水翼的振動響應時域見圖6(a)。脈沖激勵后，振幅在阻尼作用下逐漸衰減。將振動響應進行快速傅里葉變換(FFT)，得到頻域特性，見圖6(b)。主要頻率包括空化脫落渦頻率(fca)和水翼固有頻率(fn)?；诎吞匚炙箮V波，將水翼第一階彎曲模態(tài)的振動響應提取。濾波前后，固有頻率對應幅值的變化量在3%以內，證明了濾波器的可靠性。將濾波后的振動響應無量綱化，并和相同工況下的實驗結果進行比較，見圖6(c)。就振幅隨時間衰減的規(guī)律而言，模擬和實驗結果吻合良好，水力阻尼相對誤差約15%。

圖6 基于雙向流固耦合的水力阻尼識別(α=10°，σ=1.39)

不同空化數(shù)下的固有頻率和水力阻尼見圖7。隨空化數(shù)降低，空化區(qū)域與結構相接觸的面積增大，水體附加質量降低，導致第一階彎曲模態(tài)的固有頻率逐漸增大，見圖7(a)。將雙向流固耦合法計算的固有頻率與實驗結果[9]相比較，相對誤差在3.04%～9.69%之間，平均誤差為7.34%。對于第一階彎曲模態(tài)的水力阻尼，實驗[9]和雙向流固耦合模擬都得到其隨空化數(shù)降低而逐漸減小的規(guī)律，見圖7(b)。說明在小空化數(shù)工況，共振發(fā)生時更有可能引發(fā)高幅值振動。將水力阻尼模擬結果與實驗結果[9]相比較，相對誤差在6.96%～22.48%之間，平均誤差為13.56%。且模擬得到的水力阻尼基本在實驗的誤差棒范圍內，驗證了雙向流固耦合法在空化時進行水力阻尼計算的可靠性。

圖7 不同空化數(shù)下的模態(tài)參數(shù)

5.2 流動機理討論在z/w=0.5的截面上，σ=1.39時，一個空化發(fā)展周期內的氣體體積分數(shù)分布見圖8。其中，前緣附著空化剛開始脫落的時刻定義為t=0T。T為空化的發(fā)展周期，通過圖6中的空化脫落渦頻率(fca)計算得到，即T=1/fca。當t=0T～2/3T時，前緣附著空化的長度逐漸增大。當t=5/6T時，存在明顯的反向射流現(xiàn)象。反向射流從水翼尾緣向前緣發(fā)展，并在靠近前緣的位置切斷前緣附著空化，導致前緣附著空化的長度驟減。在反向射流機制下，前緣附著空化周期性的發(fā)展、回縮和斷裂，進而對水翼施加周期性的激振。

圖8 空化發(fā)展周期內的氣體體積分布(σ=1.39，z=0.5w)

在z/w=0.5的截面上，t=1/3T時，流場壓力分布見圖9(a)。在前緣附著空化和空化脫落云所在區(qū)域，絕對壓力低于飽和蒸汽壓。在逆壓梯度作用下，靠近壁面的空化區(qū)域出現(xiàn)負向速度。該現(xiàn)象可能是導致圖8中反向射流的主要原因。該時刻水翼繞流速度分布見圖9(b)，主要分為兩個區(qū)域。一是水翼壓力面以下的區(qū)域，該區(qū)域內的速度為來流速度。二是在水翼吸力面靠近壁面區(qū)域，該區(qū)域內的速度降低，甚至出現(xiàn)局部負向速度。根據(jù)前人實驗和數(shù)值結果[4，6，10-14]，水翼一階彎曲模態(tài)的水力阻尼隨流速線性增長。因此，對于水力阻尼隨空化數(shù)降低而下降的現(xiàn)象，可解釋為空化作用下的水翼吸力面速度下降。

圖9 水翼繞流特性(σ=1.39，z=0.5w，t=1/3T)

6 結論

基于雙向流固耦合數(shù)值模擬方法，在攻角α=10°、流速v=6.64 m/s和空化數(shù)σ=1.04～2.02時，對NACA 0009鈍形尾緣水翼進行水力阻尼特性研究。主要結論如下所示：(1)明確了水力阻尼雙向流固耦合計算時的網(wǎng)格尺度、時間步長和數(shù)值阻尼的定量選取原則。基于試錯法調整數(shù)值阻尼，定量確定該值為0.01時，空氣中輸入阻尼與輸出阻尼的相對誤差在4%以內。(2)基于雙向流固耦合法實現(xiàn)了空化時的水翼固有頻率和水力阻尼定量計算，得到固有頻率和水力阻尼隨空化數(shù)降低分別增大和減小的規(guī)律，兩者相較于實驗的平均誤差分別為9.69%和13.65%。(3)空化發(fā)生后在水翼吸力面可觀測到逆壓梯度，反向射流在其作用下產(chǎn)生并導致前緣空化周期性脫落；且水翼吸力面的繞流速度在空化發(fā)生后顯著下降，可能是導致水力阻尼在空化時下降的主要原因。