胡俊， 侯夏伊， 于勇

(北京理工大學 宇航學院，北京 100081)

空化是指液體中的壓力降低到飽和蒸汽壓時，產生充滿氣體空穴[1]，是一種因流體因素作用而在波體內部或在波體與固體界面上發(fā)生的液體與其蒸氣的相變過程與現象[2-3]. 空化涉及非定常、可壓縮、相間傳質和湍流脈動等許多復雜的流動現象，并且會引起水下翼型、水中機械的材料損傷，產生空泡噪聲[4-8]. 研究空化現象具有很現實的工程意義，發(fā)展空泡流的數值模擬技術是目前空化研究領域中的熱點.

研究空化問題的數值模擬方法可分為狀態(tài)方程法(barotropic model)和輸運方程法(transport equation model,TEM). 其中，輸運方程法通常不考慮流場的可壓縮性，假設氣液兩相速度相等，需要添加額外的質量輸運模型即空化模型來描述氣液兩相間的相變規(guī)律. 目前常用的空化模型均涉及經驗參數，如氣泡的數密度、非冷凝氣體質量分數和氣泡半徑等，這些參數的取值很大程度上影響了數值模擬的準確性. Asnaghi等[9]基于大渦模擬(LES)的方法，引入當地瞬時剪切應變率修正了考慮非冷凝氣體的Schnerr-Sauer空化模型，對二維NACA0009 MOD水翼的定常空化流場計算. 結果表明該修正方法可以有效消除考慮非冷凝氣體的Schnerr-Sauer空化模型中經驗參數的取值對數值模擬結果的影響.

基于這一思路，在RANS框架下模擬空化流時，也可以引入平均剪切應變率修正空化模型. 本文針對考慮非冷凝氣體的Schnerr-Sauer空化模型，將修正后的空化模型通過UDF嵌入Fluent18.2商業(yè)軟件，采用Mixture多相流模型和Realizablek-ε湍流模型對攻角為2.5°，Re=2×106的二維NACA0009 MOD水翼的定?？栈鲌鲞M行數值模擬. 通過與已有的實驗結果對比，分析了經驗參數的取值對采用兩種空化模型的計算結果的影響，驗證了引入當地流動特征修正空化模型的方法在二維定?？栈饔嬎阒械目尚行?

1 數學模型和數值計算方法

1.1 控制方程

Mixtrue模型是一種簡化的多相流模型，將流體中各相視為相互混合的單一流體，通過求解混合物的連續(xù)方程，動量方程及氣相的體積分數方程來模擬氣液兩相流體的運動. 其中，混合相連續(xù)方程及動量方程為

(1)

(2)

式中：ρm、ρl和ρv為混合相密度、液相密度和氣相密度；μm為混合相動力黏度；Ui、Uj為混合相速度. 其中，混合相密度及動力黏度為

ρm=ρvαv+ρl(1-αv)

(3)

μm=μvαv+μl(1-αv)

(4)

式中：αv為氣相體積分數；μl和μv分別為液相動力黏度及氣相動力黏度.

除混合相連續(xù)方程和動量方程外，還需要額外的氣相體積分數方程來描述氣液兩相間的質量傳遞

(5)

1.2 空化模型

1.2.1氣泡動力學基本方程

Rayleigh-Plesset單氣泡動力學模型假設液體區(qū)域內均為球形氣泡，氣泡間無相互作用:

(6)

式中：RB為氣泡半徑；PB為氣泡表面壓力；P為環(huán)境壓力；S為氣液分界面上的表面張力系數；?l為液相運動黏度；t為時間. 忽略單氣泡動力學方程中的黏性項，表面張力項及高階導數項，可得到單個氣泡半徑變化率為

(7)

1.2.2Schnerr-Sauer空化模型

Schnerr-Sauer模型基于Rayleigh-Plesset單氣泡動力學方程，將式(1)與式(5)聯(lián)立，得到氣相質量變化率與氣相體積分數變化率之間的關系為

(8)

其中氣相體積分數αv與氣泡數密度n及氣泡半徑RB的關系為

(9)

由式(9)可得到氣相體積分數變化率與氣泡半徑變化率之間的關系

(10)

Schnerr-Sauer空化模型的最終表達式為

(11)

式中：Cevap=1.0和Ccond=-1.0分別為蒸發(fā)項及凝結項系數；Pv為飽和蒸氣壓. 由于Schnerr-Sauer空化模型沒有考慮到實際流體中存在的非冷凝氣體對空化初生及發(fā)展的影響，且需要確定模型中參數氣泡數密度n的值，在模擬空化流場時存在一定的局限性.

1.2.3考慮了非冷凝氣體的Schnerr-Sauer空化模型

經研究表明，實際流體中存在的非冷凝氣體對空化的初生與發(fā)展起到了重要作用，該部分氣體在空化發(fā)生的過程中充當空化核，使液體抗拉強度下降，當液體壓力下降到一定值時，相變首先從空化核處發(fā)生. 由于原始的Schnerr-Sauer空化模型沒有考慮非冷凝氣體對空化過程的影響，在模擬空化流時存在一定缺陷. 為了提高Schnerr-Sauer空化模型模擬空化流的能力，本文將流場視為氣體、液體及非冷凝氣體三相混合物，推導出一種考慮了非冷凝氣體的Schnerr-Sauer空化模型.

Asnaghi等假設空化流場中的單個氣泡由非冷凝氣體發(fā)展而來，故流場中氣泡數密度n與非冷凝氣體數密度nnuc相等，并將非冷凝氣體的體積fnuc及體積分數αnuc定義為[8]：

(12)

(13)

參考原始的Schnerr-Sauer空化模型對氣相體積分數的定義即式(9)，將修正后的Schnerr-Sauer空化模型的氣相體積分數與非冷凝氣體體積分數之和定義為[8]

由式(14)可得到單個氣泡半徑的表達式

(15)

由于非冷凝氣體對蒸發(fā)過程的影響遠大于空泡潰滅過程，本文僅在Schnerr-Sauer空化模型的蒸發(fā)項中引入非冷凝氣體，最終得到考慮非冷凝氣體的Schnerr-Sauer空化模型：

(16)

式中：dnuc、fnuc、αnuc和nnuc分別為非冷凝氣體直徑、體積、體積分數以及數密度.

1.2.4引入剪切應變率修正模型系數

t∞=c/U∞

(18)

式中：Ce-mod為修正后的蒸發(fā)系數；t∞為主流時間尺度；c為翼型弦長；U∞為來流速度.

2 計算網格及邊界條件

計算采用了和試驗相同的二維NACA0009 MOD型水翼及流動條件，翼型攻角為2.5°，翼型弦長c=100 mm，對應雷諾數Re=2×106. 圖1給出了計算域及其邊界條件，計算區(qū)域入口距翼型前緣2.0c，出口距翼型尾緣4.0c. 如圖2所示，整個計算域被劃分為9個塊，總網格數為34 720. 翼型附近的區(qū)域采用C型結構化網格劃分，可以較好地匹配翼型的形狀，并在翼型周圍的近壁區(qū)域進行網格加密，近壁面y+值為 30～200之間，滿足壁面函數要求，計算殘差為10-4.

圖1 計算域及邊界條件Fig.1 Computational domain and boundary conditions

圖2 計算網格Fig.2 Computational meshes

數值計算邊界條件采用速度入口及壓力出口，上下邊界均為對稱條件，翼型表面采用絕熱、無滑移固壁條件. 為了保證與實驗條件一致，入口速度為U∞=20 m/s，計算了4種空化數下的流場.

3 結果與討論

3.1 非冷凝氣體數密度n對考慮非冷凝氣體的Schnerr-Sauer空化模型的影響

為了探究非冷凝氣體數密度n的取值對考慮非冷凝氣體的Schnerr-Sauer空化模型的影響，在空化數σ分別為0.75、0.80、0.85、0.90時，取dnuc為定值10-5，計算了n取106、108、1010和1012時的空化流場.

圖3所示的計算結果表明：隨著空化數σ由0.75增大到0.90，空化流場的流動狀態(tài)更加穩(wěn)定，計算出的空穴長度也有所減小. 同時，參數n的取值對空穴長度的影響隨空化數的減小而減小. 4種空化數下計算出的空穴長度均隨n的增加而增加，其中，n=1012時，空穴長度達到最大，且壓力系數的分布最接近于實驗值. 但由于引入當地流動特征修正前的空化模型的蒸發(fā)項系數為1，即使n取較大值，即n=1012時，在4種空化數下計算出的空穴長度都明顯小于實驗值. 當σ=0.75、n=1012時，采用該空化模型能夠捕捉到空穴閉合區(qū)域的壓力峰值及逆壓梯度；當σ=0.80、n=1012時，這種逆壓梯度較小且與實驗值相差較大；當σ=0.85、0.90時，即使非冷凝氣體數密度n取較大值，采用該空化模型也難以捕捉到空穴閉合區(qū)域的逆壓梯度及壓力峰值.

圖3 采用系數修正前模型計算出的翼面壓力分布Fig.3 Pressure distribution of hydrofoil calculated by using the original model

3.2 非冷凝氣體直徑dnuc對考慮非冷凝氣體的Schnerr-Sauer空化模型的影響

為了研究非冷凝氣體直徑dnuc取不同值時，考慮非冷凝氣體的Schnerr-Sauer空化模型的表現，在空化數σ=0.75，n=108和1010時，改變非冷凝氣體直徑dnuc，使dnuc=10-3、10-4、10-5、10-6和10-7m，對二維水翼的定?？栈鲌鲞M行數值模擬，如圖4所示.

圖4 采用系數修正前模型計算出的翼面壓力分布Fig.4 Pressure distribution of hydrofoil calculated by using the original model

圖4表明dnuc的取值對數值模擬結果的影響較小. 當n=108，dnuc取10-3～10-5時，以及n=1010，dnuc取10-3～10-6時，計算出的翼型表面空穴長度及壓力系數分布接近. 但即使dnuc取值較大，即dnuc=10-3時，壓力系數分布與實驗值差距依然較大，翼型表面的空穴長度明顯小于實驗值，且難以捕捉到空穴閉合區(qū)域的壓力峰值及逆壓梯度. 當n=108，dnuc=10-6、10-7以及n=1010，dnuc=10-7時，由于非冷凝氣體直徑dnuc取值過小，導致非冷凝氣體體積分數過小，空穴長度也會明顯縮短. 相比非冷凝氣體數密度n，考慮非冷凝氣體的Schnerr-Sauer空化模型對非冷凝氣體直徑dnuc的取值并不敏感.

3.3 引入平均剪切應變率修正后的考慮非冷凝氣體的Schnerr-Sauer空化模型

由于考慮非冷凝氣體的Schnerr-Sauer空化模型中的參數非冷凝氣體數密度n的取值對數值模擬結果的影響較大，非冷凝氣體直徑dnuc的取值對數值模擬的結果影響較小，本文僅探究了當參數n取不同值時，引入剪切應變率修正后的考慮非冷凝氣體的Schnerr-Sauer空化模型的表現. 在空化數σ分別為0.75、0.80、0.85和0.90時，取dnuc為定值10-5，計算了n=106、108、1010和1012時的空化流場，如圖5所示.

當σ分別為0.75、0.80和0.85，非冷凝氣體數密度n取106、108、1010和1012時，采用系數修正后的空化模型清晰地捕捉到了空穴閉合區(qū)域的壓力峰值及逆壓梯度；當σ=0.90，n=106時,難以捕捉到這種現象. 當n=106時，由于非冷凝氣體密度n的取值過小，采用修正蒸發(fā)項系數后的空化模型計算出的壓力分布雖比采用修正前的模型計算出的結果更接近實驗值，但空穴長度仍明顯小于實驗值.

當非冷凝氣體數密度n取108、1010和1012時，引入當地流動特征，采用平均剪切應變率修正考慮非冷凝氣體的Schnerr-Sauer空化模型的蒸發(fā)項系數可以有效地消除參數非冷凝氣體數密度n對空化流數值模擬的顯著影響，準確預測了空穴長度及翼面壓力分布，提高了該空化模型的適用性與準確性.

圖5 采用系數修正后模型計算出的翼面壓力分布Fig.5 Pressure distribution of hydrofoil calculated by using the modified model

4 結束語

在考慮非冷凝氣體的Schnerr-Sauer空化模型中，經驗參數非冷凝氣體數密度n的取值對數值模擬結果影響較大，非冷凝氣體直徑dnuc的取值對數值模擬結果的影響較小.

引入當地流動特征，采用剪切應變率修正考慮非冷凝氣體的Schnerr-Sauer空化模型后，參數非冷凝氣體數密度n的選取對數值模擬的結果不再有顯著影響，同時，該修正方法提高了考慮非冷凝氣體的Schnerr-Sauer空化模型的準確性，準確地預測了空化造成壓力梯度突變的位置.