基于單/雙向流固耦合的水翼水力阻尼數(shù)值預(yù)測(cè)方法研究

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

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

1 問(wèn)題提出

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

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

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

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

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

2 理論背景

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

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

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

式中： ζi為第i階水中模態(tài)的水力阻尼比；Wi為第i階模態(tài)一個(gè)振動(dòng)周期內(nèi)水翼對(duì)流體的模態(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個(gè)單位時(shí)間步內(nèi)的模態(tài)做功，J；N為一個(gè)振動(dòng)周期內(nèi)的時(shí)間步數(shù)；P為水體壓力，N·m-2；u(t)為振型運(yùn)動(dòng)速度，m·s-1；S為水翼表面面積，m2；?t為時(shí)間步長(zhǎng)，s。第i階水中模態(tài)的模態(tài)質(zhì)量計(jì)算公式為：

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

2.3 雙向耦合水力阻尼比識(shí)別方法雙向耦合采用迭代求解方法，直接求解振動(dòng)衰減曲線，進(jìn)而識(shí)別水力阻尼比。水下結(jié)構(gòu)動(dòng)力學(xué)方程[20]可表示為：

對(duì)于自由振動(dòng)，F(xiàn) (t)=0，求解方程（7）后，結(jié)構(gòu)運(yùn)動(dòng)方程[8]可表示為：

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

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

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

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

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

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

圖1 計(jì)算域

圖2 網(wǎng)格

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

4 結(jié)果與討論

4.1 單向耦合計(jì)算結(jié)果

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

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

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

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

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

圖5 不同相對(duì)幅值的水力阻尼比

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

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

4.2 雙向耦合計(jì)算結(jié)果

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

圖6 濾波前后振動(dòng)響應(yīng)（v=20m/s）

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

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

圖7 不同流速下振動(dòng)響應(yīng)的頻率

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

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

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

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

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

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

5 結(jié)論

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