曾永順，姚志峰，2，洪益平，2，王福軍，2

（1.中國農(nóng)業(yè)大學 水利與土木工程學院，北京 100083；2.北京市供水管網(wǎng)系統(tǒng)安全與節(jié)能工程技術研究中心，北京 100083）

1 問題提出

水力機械運行時承受瞬變的水力載荷，滿足共振條件時，引起劇烈振動[1-2]。在共振預測時，通常需要已知流固耦合過程中的附加質(zhì)量和水力阻尼。前者主要影響水下固有頻率，且研究較為充分，后者決定了共振幅值大小，但對其了解不夠深入[1，3]。主要制約因素是高速流動條件下的水下結構水力阻尼參數(shù)測定非常困難[4-5]。近年來，數(shù)值計算水平已得到很大的提升，但工程計算依然只能經(jīng)驗性地考慮系統(tǒng)總阻尼，對水力阻尼的定量預測方法還缺少系統(tǒng)的研究。

分離式單/雙向流固耦合主要基于計算流體力學（CFD）、計算固體力學（CSD）和計算網(wǎng)格動力學（CMD）等較為成熟的算法，在流體域和固體域分別獨立計算，只在流體與固體交界面?zhèn)鬟f載荷，這已在工程上得到較多應用[6-7]。采用該方法如何預測水力阻尼比是目前國際上水力機械領域研究的熱點問題之一。

單向流固耦合方法通過求解振動周期內(nèi)振動結構的能量耗散，間接計算得到水力阻尼比。Monette 等[6]較早開展了流固耦合條件下水翼與動水的能量交換機理研究，提出了有限元法、模態(tài)做功法和單自由度法3種水力阻尼比預測方法。其中，模態(tài)做功法得到較多的運用[8-11]。Tengs 等[8]和Bergan 等[9]利用模態(tài)做功法得到水翼的第1階彎曲模態(tài)在2.5～45 m/s之間的水力阻尼比，但在低速流動條件下會產(chǎn)生過高預測。Nennemann 等[10]和Gauthier 等[11]提出一種水體附加剛度的預測方法，對模態(tài)做功法進行了完善，并建立了水翼第1階彎曲模態(tài)水力阻尼比與流速之間的線性關系，相對偏差在10%以內(nèi)，但是該模擬沒有考慮振型幅值對計算結果的影響。Tengs 等[8]計算得到了水翼在不同振型幅值下的水力阻尼比，結果表明當振型幅值從0.05 mm 增大到2.5 mm時，水力阻尼比顯著增大，在10 m/s 流速下相對偏差可達300%。

雙向流固耦合方法考慮了流固耦合交界面變形的影響，通過流場非定常計算（CFD）和結構瞬態(tài)動力學（CSD）的迭代計算，借助動網(wǎng)格變形，直接獲取特定激勵下自由振動幅值[7]。該方法相對單向耦合需要更多的計算資源，但在水力阻尼比求解的理論上更加完備?；谠摲椒ǎ琇iaghat 等[12]得到水力阻尼比與流速之間的線性關系，但沒有對邊界層速度和尾跡區(qū)渦激振動進行細致分析，與實驗結果相比較，平均相對偏差在12%左右。曾永順等[13]較為系統(tǒng)地分析了網(wǎng)格、時間步長、近壁區(qū)處理模式、激勵方式和數(shù)值阻尼等關鍵參數(shù)對計算結果的影響，顯著提高了水翼水力阻尼比的計算精度，在2～15 m/s 流速下，模擬結果相對實驗結果的最大偏差為8.82%。另外，渦激振動[14]和升力作用[15]對水力阻尼比識別精度有重要的影響。針對非對稱尾部形狀水翼振動響應的時變特性，曾永順等[16]提出了帶通濾波加振動響應校準，對傳統(tǒng)自由振動衰減法進行了改進。

在水力阻尼比預測時，目前缺乏對單/雙向流固耦合方法的系統(tǒng)比較。本文采用上述兩種方法，以NACA0009 鈍形尾部形狀水翼為分析對象，預測不同流速下的水力阻尼比，分析兩種預測方法的可靠性，比較計算精度和耗時，討論了實際工程中預測方法的選擇策略。

2 理論背景

2.1 流場控制方程為了對邊界層和尾跡區(qū)流動進行精確模擬，流場計算采用轉(zhuǎn)捩SST 模型[17-18]。相對于傳統(tǒng)SST k-ω模型，該模型在k 方程中耦合γ-Reθt轉(zhuǎn)捩模型進行修正，控制方程[17]如下：

式中：k為湍動能，m2·s-2；ω為比耗散率，s-1；ρ為流體密度，kg·m-3；τij為壁面剪切應力，N·m-2；μ為動力黏度系數(shù)，N·s·m-2；μt為渦黏系數(shù)，N·s·m-2；γeff為有效間歇因子（-）；β*和σk為常系數(shù)；t為時間，s；u為速度，m·s-1。

2.2 單向耦合水力阻尼比識別方法假設水翼的振型及固有頻率不隨流速發(fā)生改變。獲得水翼對流體的模態(tài)做功后，根據(jù)水翼在靜水中的模態(tài)質(zhì)量、固有頻率和振型幅值能夠計算得到水力阻尼 比[19]：

式中： ζi為第i階水中模態(tài)的水力阻尼比；Wi為第i階模態(tài)一個振動周期內(nèi)水翼對流體的模態(tài)做功，J；mf，i為第i階水中模態(tài)的模態(tài)質(zhì)量，kg；ωnf，i為第i階水中模態(tài)的角固有頻率，rad·s-1；ωnf，i=2πfnf，i，fnf，i為第i階水中模態(tài)固有頻率，Hz；A0，i為振型的幅值，m；Wn第n個單位時間步內(nèi)的模態(tài)做功，J；N為一個振動周期內(nèi)的時間步數(shù)；P為水體壓力，N·m-2；u(t)為振型運動速度，m·s-1；S為水翼表面面積，m2；?t為時間步長，s。第i階水中模態(tài)的模態(tài)質(zhì)量計算公式為：

式中：ma，i、 fna，i和Ei分別為第i階空氣中模態(tài)的模態(tài)質(zhì)量（kg）、固有頻率（Hz）和網(wǎng)格單元儲存的總動能（J）。

2.3 雙向耦合水力阻尼比識別方法雙向耦合采用迭代求解方法，直接求解振動衰減曲線，進而識別水力阻尼比。水下結構動力學方程[20]可表示為：

對于自由振動，F(xiàn) (t)=0，求解方程（7）后，結構運動方程[8]可表示為：

式中： A為幅值，m；ωdf，i為第i階水中模態(tài)的阻尼角固有頻率，為初相位，rad。

3 計算設置

3.1 物理模型研究對象為NACA0009 鈍形尾部形狀水翼，物理模型及計算域如圖1所示。水翼弦長、翼展和尾部厚度分別為100、150和3.22 mm。水翼材料為不銹鋼，密度、彈性模量和泊松比分別為7700 kg/m3、215 GPa和0.3。計算域為750 mm×150 mm×150 mm的方形水洞段，水翼放置攻角為0°。

3.2 邊界條件水翼一端被固定約束，另一端自由。邊界條件為5～20 m/s的速度進口，2.5 Bar的靜壓出口。為了讓水翼自由端的約束方式和繞流情況與真實實驗相符合，將自由端對應的邊界設置為對稱面，其他邊界為無滑移壁面。流場非定常計算的收斂標準為均方根殘差小于1×10-5，最大迭代步數(shù)為10。非定常計算尾部出現(xiàn)周期性脫落渦后，將該結果作為流固耦合計算的初始文件?；趧泳W(wǎng)格變形技術，實施單/雙向流固耦合計算。

對于單向流固耦合，將水翼水中第1階彎曲模態(tài)的振型導入到流場計算域中，如圖1所示。采用差值法將振型賦給流場的網(wǎng)格節(jié)點，并根據(jù)振型對應的固有頻率做周期性振動。將水翼表面單位面積的模態(tài)做功進行面積分，累計10個振動周期的結果，取平均后得到1個周期的模態(tài)做功。

對于雙向流固耦合，為激勵第1階彎曲模態(tài)，在水翼表面中心線上沿y 軸正方向施加200 N的瞬態(tài)激勵，如圖1所示，激勵時間小于1/4個振動周期。耦合迭代的收斂標準為均方根殘差小于1×10-4，最大迭代步數(shù)為30。監(jiān)測流場的網(wǎng)格位移獲得結構的振動響應，監(jiān)測點在尾部中心。

圖1 計算域

圖2 網(wǎng)格

3.3 計算參數(shù)無關性分析采用結構化網(wǎng)格離散流場和結構場，如圖2所示。采用O 形拓撲結構劃分邊界層，第一層網(wǎng)格厚度為2×10-6m，網(wǎng)格拓展比為1.05。流場和結構場網(wǎng)格質(zhì)量分別在0.35和0.15 以上，且有85%以上的網(wǎng)格質(zhì)量大于0.8。前期研究[13]中對該計算域進行網(wǎng)格收斂性分析和時間步長無關性檢查，結果表明當網(wǎng)格單元數(shù)大于246.5 萬，時間步長小于2×10-5時對脫落渦頻率沒有影響；且當結構場網(wǎng)格單元數(shù)為1.4 萬時，模態(tài)分析結果與實驗結果相對偏差在0.5%以內(nèi)。本研究采用上述參數(shù)進行流固耦合計算。在20 m/s 流速下，y+平均值為1.29，滿足轉(zhuǎn)捩SST模型在近壁區(qū)直接求解的要求。

4 結果與討論

4.1 單向耦合計算結果

4.1.1 模態(tài)分析 水翼低階模態(tài)振型及固有頻率如圖3所示。對于相同模態(tài)，空氣中振型與水中振型一致，但是在水體附加質(zhì)量作用下，水中固有頻率顯著低于空氣中。對于不同模態(tài)，圖3（a）和圖3（c）為彎曲變形，前者沿流向只有一個節(jié)線，后者有兩個節(jié)線，兩個模態(tài)分別命名為第1和第2階彎曲模態(tài)。圖3（b）為扭轉(zhuǎn)變形，沿展向只有一個節(jié)線，命名為第1階扭轉(zhuǎn)模態(tài)。將模擬結果與實驗結果[5]相比較，第1階彎曲和扭轉(zhuǎn)模態(tài)固有頻率的相對誤差分別為0.44%和0.58%。根據(jù)式（5）和式（6），計算得到第1階彎曲模態(tài)在空氣中的模態(tài)質(zhì)量為0.2486 kg，在水中的模態(tài)質(zhì)量為0.5815 kg。

圖3 低階模態(tài)振型及固有頻率

4.1.2 振幅無關性驗證 定義相對幅值A?為振型幅值與第1層網(wǎng)格厚度的比值。10 m/s 流速下，A?=10和1000時的模態(tài)做功特性如圖4所示。對于時域，相對振幅越大，水翼對流體的模態(tài)做功越大，能量耗散的越多。對于頻域，當A?=1000時頻率僅有2倍的第1階彎曲模態(tài)固有頻率；當A?=10時，除該頻率為主頻外，還有第1階彎曲模態(tài)固有頻率和脫落渦頻率。結果表明，當相對振幅過大時，脫落渦對水力阻尼比的貢獻被忽略。將脫落渦頻率模擬結果與實驗結果[21]相比較，相對偏差為2.52%。

圖4 不同相對幅值的模態(tài)做功特性

在10 m/s 流速下，不同相對幅值對應的水力阻尼比如圖5所示。隨著相對幅值的增大，水力阻尼比逐漸增大，這一趨勢與Tengs 等[8]的結果相一致。當相對幅值小于50時對水力阻尼比計算結果沒有影響，并將該相對幅值用于不同流速下水力阻尼比的計算。

圖5 不同相對幅值的水力阻尼比

4.1.3 單向耦合水力阻尼比 在5～20 m/s 流速下，基于單向流固耦合計算得到的第1階彎曲模態(tài)的水力阻尼比，如表1所示。結果表明隨著流速增大，水力阻尼比逐漸增大，這一趨勢與文獻[5]的實驗結果相一致。在模擬流速范圍內(nèi)，相對文獻[5]的實驗平均偏差為11.42%。

表1 單向流固耦合得到的不同流速下的水力阻尼比

4.2 雙向耦合計算結果

4.2.1 振動特性為了消除渦激振動和其他模態(tài)的振動對第一階彎曲模態(tài)水力阻尼比識別的影響，采用帶通濾波將第一階彎曲模態(tài)的位移從原始位移中分離出來。20 m/s 流速下，水翼振動位移在濾波前后的時域及頻域如圖6所示。如圖6（b）所示，濾波前振動響應有4個頻率成分，包括第1階彎曲、第1階扭轉(zhuǎn)、第2階彎曲模態(tài)固有頻率和脫落渦頻率。在外部激勵作用下，第1階彎曲模態(tài)的幅值顯著大于其他頻率。經(jīng)過濾波處理后，振動響應的頻率成分僅剩第1階彎曲模態(tài)固有頻率。

圖6 濾波前后振動響應（v=20m/s）

在5～20 m/s 流速下，固有頻率及脫落渦頻率隨流速的變化如圖7所示。第1階彎曲、扭轉(zhuǎn)及第2階彎曲模態(tài)固有頻率基本不隨流速發(fā)生改變，相對變化范圍分別在0.43%、2.61%和3.95%以內(nèi)。對于第1階彎曲和扭轉(zhuǎn)模態(tài)的固有頻率，與文獻[5]的實驗結果相比較，平均相對偏差分別為4.36%和4.07%。脫落渦頻率隨流速線性增長，相對文獻[21]的實驗結果平均偏差為4.24%。模擬流速對應的脫落渦頻率均大于第1階彎曲模態(tài)固有頻率，錯開了共振工況。

4.2.2 雙向耦合水力阻尼比 20 m/s 流速下，水力阻尼比的識別如圖8所示。將濾波后振動響應的上下峰值點取出，根據(jù)自由振動衰減法進行函數(shù)擬合?；谏舷路逯迭c的水力阻尼比分別為0.076 78和0.079 22，二者的相對偏差為3.18%，將二者取平均后得到最終的水力阻尼比為0.078。為避免濾波對振動幅值的影響，在水力阻尼比識別時舍棄第一個振動響應的峰值點。

圖7 不同流速下振動響應的頻率

在5～20 m/s 流速下，水力阻尼比隨流速的變化如圖9所示。數(shù)值模擬計算得到的水力阻尼比隨流速線性增長的趨勢與文獻[5]的實驗結果相一致，雙向流固結果與實驗結果相比較，平均相對偏差為4.95%。在10 m/s 流速下，單/雙向流固耦合計算得到水力阻尼比與實驗結果相比較，相對偏差分別為17.07%和0.64%。原因是該流速接近共振工況，此時結構振動對流場流動有較為強烈的反作用影響，而雙向流固耦合數(shù)值模擬與真實流動更加吻合，因此模擬精度顯著高于單向流固耦合。

圖8 水力阻尼比的識別方法

圖9 不同流速下的水力阻尼比

4.3 討論對比單/雙向流固耦合數(shù)值模擬方法，前者假設振型和固有頻率不隨流速發(fā)生改變[11]，且不考慮結構變形的反作用影響[8]，而后者理論上更加完善。在20 m/s 流速下，基于雙向流固耦合計算得到的第一階彎曲模態(tài)的振型及固有頻率如圖10所示。與靜水中模態(tài)分析結果相比較，振型相一致，固有頻率相對偏差為3.65%。在低階模態(tài)下，水翼固有頻率在5～20 m/s 流速下的相對變化量在3.95%以內(nèi)，這也證明了單向流固耦合方法的假設的合理性。但是，在大變形情況下，單向流固耦合忽略了渦激振動對水翼能量耗散的貢獻，會導致水力阻尼比的預測誤差過大，此時需要雙向流固耦合進行求解。

圖10 20 m/s 流速下第1階彎曲模態(tài)振型

基于相同的計算設置，在10 m/s 流速下采用48 線程進行第1階彎曲模態(tài)的水力阻尼比預測，單/雙向流固耦合計算一個振動周期分別需要2.1 h和31.3 h，后者需要的計算資源為前者的15倍。實際工程中，水力機械結構較為復雜，優(yōu)先推薦單向流固耦合方法預測水力阻尼比；若在大變形工況下，或需要高精度的預測結果，適當簡化工程對象，可采用雙向流固耦合方法。

5 結論

采用單/雙向流固耦合數(shù)值模擬方法，在5～20 m/s 流速下對NACA0009 鈍形尾部水翼的水力阻尼比進行了預測，驗證了兩種計算方法的可靠性，對比分析了兩種方法的計算精度和計算時間，并對實際工程中水力阻尼比預測方法的選擇進行了討論。主要結論如下所示：（1）靜水中水翼第1階彎曲模態(tài)與20 m/s 流速下對應振型一致，且低階模態(tài)固有頻率在5～20 m/s之間的相對變化在3.95%以內(nèi)。據(jù)此，可假設振型及固有頻率與流速無關，開展單向流固耦合的水力阻尼比預測。（2）單/雙向流固耦合數(shù)值模擬計算得到的第1階彎曲模態(tài)水力阻尼比隨流速線性增長的趨勢與實驗一致，平均相對偏差分別為11.42%和4.95%，后者計算精度高于前者，但計算資源需求量是前者的15倍。（3）工程上，預測復雜水力機械的水力阻尼比時，優(yōu)先推薦單向流固耦合方法，若需要高精度或結構產(chǎn)生大變形時，則需要簡化分析對象，并采用雙向流固耦合方法。