胡俊， 侯夏伊， 于勇

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

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

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

基于這一思路，在RANS框架下模擬空化流時，也可以引入平均剪切應(yīng)變率修正空化模型. 本文針對考慮非冷凝氣體的Schnerr-Sauer空化模型，將修正后的空化模型通過UDF嵌入Fluent18.2商業(yè)軟件，采用Mixture多相流模型和Realizablek-ε湍流模型對攻角為2.5°，Re=2×106的二維NACA0009 MOD水翼的定?？栈鲌鲞M(jìn)行數(shù)值模擬. 通過與已有的實驗結(jié)果對比，分析了經(jīng)驗參數(shù)的取值對采用兩種空化模型的計算結(jié)果的影響，驗證了引入當(dāng)?shù)亓鲃犹卣餍拚栈Ｐ偷姆椒ㄔ诙S定?？栈饔嬎阒械目尚行?

1 數(shù)學(xué)模型和數(shù)值計算方法

1.1 控制方程

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

(1)

(2)

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

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

(3)

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

(4)

式中：αv為氣相體積分?jǐn)?shù)；μl和μv分別為液相動力黏度及氣相動力黏度.

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

(5)

1.2 空化模型

1.2.1氣泡動力學(xué)基本方程

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

(6)

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

(7)

1.2.2Schnerr-Sauer空化模型

Schnerr-Sauer模型基于Rayleigh-Plesset單氣泡動力學(xué)方程，將式(1)與式(5)聯(lián)立，得到氣相質(zhì)量變化率與氣相體積分?jǐn)?shù)變化率之間的關(guān)系為

(8)

其中氣相體積分?jǐn)?shù)αv與氣泡數(shù)密度n及氣泡半徑RB的關(guān)系為

(9)

由式(9)可得到氣相體積分?jǐn)?shù)變化率與氣泡半徑變化率之間的關(guān)系

(10)

Schnerr-Sauer空化模型的最終表達(dá)式為

(11)

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

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

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

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

(12)

(13)

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

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

(15)

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

(16)

式中：dnuc、fnuc、αnuc和nnuc分別為非冷凝氣體直徑、體積、體積分?jǐn)?shù)以及數(shù)密度.

1.2.4引入剪切應(yīng)變率修正模型系數(shù)

t∞=c/U∞

(18)

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

2 計算網(wǎng)格及邊界條件

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

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

圖2 計算網(wǎng)格Fig.2 Computational meshes

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

3 結(jié)果與討論

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)束語

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

引入當(dāng)?shù)亓鲃犹卣?，采用剪切?yīng)變率修正考慮非冷凝氣體的Schnerr-Sauer空化模型后，參數(shù)非冷凝氣體數(shù)密度n的選取對數(shù)值模擬的結(jié)果不再有顯著影響，同時，該修正方法提高了考慮非冷凝氣體的Schnerr-Sauer空化模型的準(zhǔn)確性，準(zhǔn)確地預(yù)測了空化造成壓力梯度突變的位置.