胡俊, 侯夏伊, 于勇
(北京理工大學(xué) 宇航學(xué)院,北京 100081)
空化是指液體中的壓力降低到飽和蒸汽壓時,產(chǎn)生充滿氣體空穴[1],是一種因流體因素作用而在波體內(nèi)部或在波體與固體界面上發(fā)生的液體與其蒸氣的相變過程與現(xiàn)象[2-3]. 空化涉及非定常、可壓縮、相間傳質(zhì)和湍流脈動等許多復(fù)雜的流動現(xiàn)象,并且會引起水下翼型、水中機械的材料損傷,產(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行數(shù)值模擬. 通過與已有的實驗結(jié)果對比,分析了經(jīng)驗參數(shù)的取值對采用兩種空化模型的計算結(jié)果的影響,驗證了引入當(dāng)?shù)亓鲃犹卣餍拚栈P偷姆椒ㄔ诙S定??栈饔嬎阒械目尚行?
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氣泡動力學(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空化模型的最終表達式為
(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)空化核,使液體抗拉強度下降,當(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)可得到單個氣泡半徑的表達式
(15)
由于非冷凝氣體對蒸發(fā)過程的影響遠大于空泡潰滅過程,本文僅在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∞為來流速度.
計算采用了和試驗相同的二維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ū)域進行網(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ù)下的流場.
為了探究非冷凝氣體數(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時,空穴長度達到最大,且壓力系數(shù)的分布最接近于實驗值. 但由于引入當(dāng)?shù)亓鲃犹卣餍拚暗目栈P偷恼舭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
為了研究非冷凝氣體直徑dnuc取不同值時,考慮非冷凝氣體的Schnerr-Sauer空化模型的表現(xiàn),在空化數(shù)σ=0.75,n=108和1010時,改變非冷凝氣體直徑dnuc,使dnuc=10-3、10-4、10-5、10-6和10-7m,對二維水翼的定??栈鲌鲞M行數(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的取值并不敏感.
由于考慮非冷凝氣體的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
在考慮非冷凝氣體的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ù)測了空化造成壓力梯度突變的位置.