基于Mixture模型對翼形空化的數值模擬

肖鑫明，桂 林

(四川大學水利水電學院，成都 610065)

0 引 言

在水力機械中翼形空化是由于水流繞流葉片引起壓力降低而產生的[1]，當流場中某點的壓力低于液體介質工作溫度所對應的汽化壓力時，該處液體發生相變，產生氣泡，氣泡生長直至潰滅的現象即為空化現象[2]。當潰滅發生在水力機械固體表面附近時，流體中連續破滅的空泡所產生的高壓力不斷作用于固體表面，會造成流體機械表面的損壞，并兼有振動、噪聲及性能下降等現象[3]。葉片的進口處截面便是一個簡單翼形，可以通過研究翼形來研究這種汽蝕性能的變化。于鳳榮等人通過在翼形上打孔作了相關研究，利用單相模型證實了打孔可以破壞翼形上方區的低壓流場[4]。王鑫等人進一步對孔的布置方式和位置、孔的數量和深度等參數對翼形的影響進行了具體研究[5]。戴月進等人研究了翼形表面粗糙帶對空化效果的抑制[6]。2018年李雅琴等人利用Mixture多相流模型模擬了二維翼形的空化情況[7]。但以上翼形空化的研究要么限于單相流的研究，要么僅僅是利用多相流模型進行初步的模擬，隨著CFD技術的不斷發展，現在的多相流模型已能達到足夠的精度來勝任進一步的模擬，本文將基于Mixture模型對微孔及翼形表面粗糙度對三維翼形的空化影響作進一步的探究。

1 翼形流動特性

1.1 微孔和翼形表面粗糙度影響空化原理

微孔改善翼形性能，主要是因為翼形表面微孔提前將固體表面的層流邊界層轉換為湍流邊界層，后者有更好的附著性，從而推遲了固體表面邊界層的脫落與分離，減小了翼形后漩渦區的面積與強度，大幅降低了翼形的壓差阻力[5]。而在翼形表面施加粗糙帶可有效增加翼形表面附近流場的湍動能，改變壓力分布，從而影響水翼初生空化的發生[6]。

1.2 升力、阻力特性

翼形升力公式為：

(1)

翼形阻力公式為：

(2)

式中：Cl為升力系數；Cd為阻力系數；ρ為流體密度，kg/m3；v∞為來流速度，m/s；S為翼形平面面積，m2。

在改善翼形空化性能的同時，還要注意翼形的水動力性能，往往我們會通過升阻比這一新的特性系數來確定翼形的升阻性能以及來改善優化它的升力性能[8,9]。升阻比一般用L/D表示，升阻比又稱“舉阻比”。翼形在運行過程中，在同一迎角的升力與阻力的比值與翼形的迎角、流體的速度等參數有關，此值愈大說明該翼形工作面的水力動力性能愈好[10]。因此在后文比較完幾種方案的空化情況后，還會對升阻比進行分析。

2 數值模擬過程

2.1 計算區域與網格劃分

(1)模型建立。在Solideworks中導入翼形NACA-0016的數據，翼形長為100 mm，位于一個長方體的流道內，x=0～300,y=-40～40,z=0～240,翼形頂部距離入口80 mm。為了更好地觀察升力系數和阻力系數的變化，將翼形的攻角設置為8°，見圖1。在微孔試驗中還會為翼形低壓區添加微孔，孔徑為1 mm。為了更好地觀察添加微孔后的變化，在y方向上添加等距的3個微孔。

圖1 模型尺寸Fig.1 Model size drawing

(2)空化區的確定。圖2為數值計算得到的無孔翼形工作面空泡份額分布曲線，可以確定開孔的空化區，選擇開孔位置位于距離翼形前緣10 mm的位置進行開孔。

圖2 翼形空泡分布曲線Fig.2 Airfoil bubble distribution curve

(3)網格劃分。利用ANSYS Meshing軟件對模型計算區域采用非結構網格進行劃分(見圖3、圖4)，并在翼形壁面進行邊界層網格的劃分，且在翼形工作面邊界、前緣、后緣都進行加密設置，所得網格數量在100 萬個左右。接下來對網格進行疏密的敏感性分析。由于計算機計算能力可較快計算出本模型100 萬個網格的結果，故在敏感性分析中只對網格進行再加密至200 萬個，所得的翼形空化情況及生升阻力系數的值差別非常細微，但計算收斂時間卻增加了3倍左右。平衡計算精度與計算時間，最后選取了網格在100 萬個左右的方案。

圖3 網格整體劃分Fig.3 The overall grid division diagram

圖4 網格細部劃分Fig.4 The detail grid division diagram

2.2 湍流模型

SSTk-ε模型是用Menter[11]發展而來的，以便其在廣泛的領域中可以優于k-ε模型，而且在帶有壓力梯度的流動中其計算結果更加準確。占梁梁[12]通過對比幾個方程得出從模擬空化狀態下翼形能量特性的角度來看，優先選擇SST模型較好的結論，故本文的模擬選用SSTk-ε模型。該模型將k-ω模型和k-ε模型以加權平均的方式結合起來，兼具k-ε模型對遠場條件依賴性較小和k-ω模型在近壁面模擬準確度較高的特點，方程表達式為：

(3)

(4)

式中:方程右側的前3項分別為湍流生成項、耗散項和擴散項；ω方程右側的C項代表交叉擴散項；各常、系數的取值詳見文獻[13]。

2.3 多相流Mixture模型

該模型是一種簡化后的多相流模型，主要用于模擬各相存在不同速度的多相流動計算。由于假定了流動在短空間尺度上局部的平衡，相與相間的耦合是很強的[14]。因此可以應用于模擬存在強烈耦合的各相同性多相流和各相似或相同速度運動的多相流。本文研究選用此模型基于空化兩相流來進行數值計算。其主要方程如下。

氣相輸運方程：

(5)

氣泡動力學方程：

(6)

2.4 空化模型

Singhal完全空化模型[15]綜合考慮了湍動能以及水中存在的氣核的影響，因此得到了廣泛的使用。其質量輸運方程中2個參數計算方法如下：

P≤Pv

(7)

(8)

式中：Fvap為蒸發系數；Fcond為凝結系數；fv為氣相質量分數；fg為不可凝氣體質量分數；k為湍動能；Pv為飽和蒸汽壓。

2.5 邊界條件

使用Fluent軟件對三維翼形空化進行數值模擬時，邊界條件設置如下:①入口邊界條件水的速度設置為10 m/s，空泡份額為0，根據來流條件計算設置湍流強度為5.11%，水力直徑為120 mm。②出口邊界為了防止回流設置成壓力出口，壓力大小與大氣壓一致，空泡份額設置為0，湍流強度和水力直徑設置與入口相同。③沿流向兩側邊界采用對稱邊界條件，各變量沿法向的分量為零。④上下兩邊面及微型孔壁采用無滑移光滑固壁邊界條件。⑤翼形四周的壁面在表面粗糙影響試驗中采用粗糙高度0.3 mm，其余試驗均采用無滑移光滑固壁邊界條件。

2.6 離散方法

采用Fluent求解器中基于壓力的求解器進行計算，采用Coupled算法實現對壓力場和速度場的迭代求解，對壓力方程采用Presto格式，其余采用Quick格式進行離散。

2.7 試驗方案

本文總共采取4種方案，第1種方案不作任何處理(A方案)，第2種方案給翼形表面添加0.3 mm粗糙高度(B方案)，第3種方案為在翼形空化區開微孔(C方案)，第4種方案在翼形開微孔方案上給翼形表面添加0.3 mm粗糙高度(D方案)。

3 計算結果及分析

3.1 模擬結果的收斂性

按照給定條件對相關參數設置后，分別對4種方案進行模擬計算。由圖5的殘差監測曲線可以看出，中途出現了幾次波動，是由于計算過程為了促進收斂不斷調整殘差而產生的，結果除Continuity之外其余的殘差收斂值皆小于1e-04。判斷Continuity是否收斂還可以通過計算進口與出口質量流量差來判斷，如果差值小于總流量的0.5%，也可以判定為收斂。由圖6可知，差值遠小于總流量的0.5%，質量流量達到了穩定值。由此說明通過區域的質量流量是連續的，滿足質量守恒定律。

圖5 殘差監測曲線Fig.5 Residual monitoring curve

圖6 質量流量連續性檢查Fig.6 Quality flow continuity inspection

3.2 翼形表面壓力及空化分布

結合圖7的壓力分布、圖8的速度分布、圖9的翼形表面空化分布可以看出，翼形吸力面前段流體速度增大導致這個區域壓力降低，當壓力下降至飽和蒸汽壓之下時，這個局部區域中的水發生相變形成氣相空腔，可以看到，先是一部分氣相空腔潰滅，再達到一個完全氣相空腔潰滅，空化的前段為氣液兩相混合區，中段充滿氣相，后半段是氣液兩相混合區，也就是本文主要研究的翼形空化現象，另外從速度分布圖中還可看出在翼形的尾部形成回射流漩渦。

圖7 方案A側面壓力分布Fig.7 The lateral pressure distribution of case A

圖8 方案A側面速度分布Fig.8 The lateral velocity distribution of case A

圖9 方案A翼形表面空化分布Fig.9 Airfoil surface cavitation distribution of case A

從圖8、圖9對比圖7來看，無論是在翼形表面增加微孔還是增加翼形表面的粗糙程度，均對翼形空化有一定的改善效果，增加翼形表面的粗糙程度可以將整體的空化情況降低，而增加微孔僅可以降低微孔及微孔之后區域的空化情況，這與文獻[5]和文獻[6]中用單相流模型模擬所得出的結論一致，可以判斷多相流模型模擬翼形空化的正確性。從圖10可以看出，當在翼形表面增加微孔并且增加翼形表面的粗糙程度時，翼形空化的情況改善得最好。方案B、C、D的翼形表面空化分布分別見圖10、圖11、圖12。

圖10 方案B翼形表面空化分布Fig.10 Airfoil surface cavitation distribution of case B

圖11 方案C翼形表面空化分布Fig.11 Airfoil surface cavitation distribution of case C

圖12 方案D翼形表面空化分布Fig.12 Airfoil surface cavitation distribution of case D

3.3 3種情況對水力特性的影響

通過表1和圖13可以看出，方案B增加了翼形表面粗糙度后，翼形的升力系數減少了7%左右，阻力系數增加了5%左右，升阻比下降11%左右。增加微孔方案C在升力系數和阻力系數上均與未增加微孔和表面粗糙度的方案A持平。方案D表明增加微孔在一定程度上提升了僅增加表面粗糙度的升力系數，但同時也增加了阻力系數，導致最后升阻比和方案B基本一致。基于以上的對比可知，增加微孔對翼形的主要水力特性升阻力基本不會造成什么影響，而增加翼形表面的粗糙度會對翼形的主要水力特性升阻力造成一定的削弱影響。

表1 翼形表面升阻力系數對比Tab.1 Comparison table of lift drag coefficient of airfoil surface

圖13 4種方案對應的變化曲線Fig.13 Corresponding variation

4 結 論

根據以上的數值模擬及分析可以得出以下結論。

(1)本文對三維翼形進行了空化模擬，初步表明現在的多相流模型及空化模型能夠較好地模擬出三維翼形空化的情況，而不僅僅局限于二維或者單相流的模擬。

(2)本文對增加翼形表面粗糙度及增加微孔對翼形空化的影響進行模擬，模擬結果表明2者均能在一定程度上改善翼形空化的性能，但是增加翼形表面粗糙度對翼形的水力特性有一定程度上的減弱影響。至于孔數的多少、位置及粗糙度的大小等因素對翼形空化的影響還有待進一步研究。

(3)將本文結果初步應用到流體機械上，可優先考慮給葉片空化區增加微孔，再考慮增加葉片空化區的粗糙程度，這樣可以減少這些措施對葉片水動力性能的影響，但具體實施還需視具體情況而定。