楊其文，張煒，李曉俊，陳峰，朱祖超

(1.中石化洛陽工程有限公司，河南 洛陽 471003；2.浙江理工大學浙江省流體傳輸技術(shù)研究重點實驗室，浙江 杭州 310018；3.浙江機電職業(yè)技術(shù)學院機械技術(shù)系，浙江 杭州 310053)

空化通常出現(xiàn)在液體局部靜壓低于當?shù)責崃W狀態(tài)下的飽和蒸汽壓而發(fā)生氣化的過程，是一種復雜的相變現(xiàn)象[1].通常對于熱力學不敏感的流體(例如常溫下的水)，空化被假設(shè)為等溫過程[2].然而，對于熱力學敏感流體(例如低溫流體，氟化酮以及高溫下的水)，這種等溫假設(shè)是無效的[3].由于氣液密度較小，當氣化形成同樣大小的空腔時，熱力學敏感流體需要比熱力學不敏感流體氣化更多的液體.

低溫離心泵通常承擔所在系統(tǒng)“心臟”的作用，涉及的空化問題直接決定低溫流體裝備及系統(tǒng)的穩(wěn)定性.由于流道結(jié)構(gòu)復雜、影響因素繁多，低溫離心泵的空化是困擾所在行業(yè)發(fā)展的“卡脖子”難題，它不僅影響低溫泵的水動力性能，同時造成高幅度的結(jié)構(gòu)振蕩.美國Fastrac火箭和日本H-ⅡA火箭[4]在研制過程中均出現(xiàn)了因誘導輪低溫空化造成的破壞.

與常溫水介質(zhì)空化相比，熱敏流體空化的本質(zhì)區(qū)別是流體介質(zhì)的物理屬性隨溫度變化敏感.SARSDY等[5]通過試驗首次發(fā)現(xiàn)水和氟利昂中空化的差異.發(fā)現(xiàn)當空化在水中發(fā)生時，氣液界面清晰透明，然而氟利昂發(fā)生空化時，空泡界面模糊，呈泡霧狀.HORD[6]在NASA的資助下進行了一系列有關(guān)水翼、文丘里管和回轉(zhuǎn)體的熱敏空化試驗.試驗數(shù)據(jù)常常被用來驗證數(shù)值模擬的有效性.KELLY等[7]通過高速攝影技術(shù)研究了不同空化數(shù)、不同入口溫度以及不同攻角下氟化酮繞NACA 0015水翼的非穩(wěn)態(tài)空化特征，并記錄了水翼表面的壓力數(shù)據(jù)，但并未給出溫度數(shù)據(jù).LE等[8]使用修正的Merkle空化模型數(shù)值研究了液氫和液氮的空化流動，并考慮了可壓縮性.LIANG等[9]基于輸運方程的修正空化模型數(shù)值研究了不同溫度的水繞NACA 0015水翼的空化流動.SUN等[10]數(shù)值研究了氟化酮和液氮繞NACA 0015水翼空化流動.發(fā)現(xiàn)溫度分布與空腔的演化有關(guān)，而且對于同樣的空化數(shù)，液氮的溫降更明顯.朱佳凱等[11]通過試驗和數(shù)值模擬觀察了液氮在文氏管和水翼中的空化特性，測量了空化過程的流量、溫度及壓力的動態(tài)變化，證實了壓比大于2.23時，空化潰滅產(chǎn)生的壓力波為空化脫落的主要機理.

文中目的是驗證修正的ZGB空化模型在熱敏空泡流中的有效性，研究氟化酮空化流動的熱力學效應以及動態(tài)脫落過程，并以此揭示熱力學空化流動特征.

1 數(shù)學模型及數(shù)值方法

1.1 控制方程

熱力學空化流動是一種復雜的多相流問題.考慮到計算資源和精度，均勻混合流模型是一種比較合理的選擇.均勻混合模型是將氣相和液相視為一種均勻混合物.兩相共用一套控制方程，為了降低計算資源，兩相之間的相互作用被忽略.質(zhì)量方程、動量方程、能量方程和輸運方程的計算公式為

(1)

(2)

(3)

(4)

式中：ρm，u，p和μm分別為混合密度、速度、壓強和混合層流黏度；i，j，k分別為直角坐標系中的x，y，z方向；T，C，K和L分別為溫度、比熱容、熱導率和氣化潛熱.

它們的定義式如下：

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

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

(5)

式中：下標v和l分別表示氣相和液相.

1.2 基于輸運方程的空化模型

輸運方程的源項由Rayleigh-Plesset方程推導得出，即

(6)

式中：RB為單個球形氣泡半徑；S0為表面張力；pv為當?shù)販囟认碌娘柡驼羝麎?忽略黏性項、二階導數(shù)項和表面張力項，可以得到氣泡半徑對時間求導項：

(7)

式(7)考慮了壓力變化對氣泡大小的影響，但是忽略了溫度變化對氣泡大小的影響.對于熱力學敏感流體(如液氫、液氮和氟化酮等)而言，其物質(zhì)屬性對溫度的變化非常敏感.氣泡壁面的自然對流引起的溫差計算公式為

(8)

式中：hb為對流換熱系數(shù).

進一步推導得出：

(9)

其中，

(10)

式中：C0為經(jīng)驗系數(shù)；熱擴散系數(shù)λl=Kl/ρlCl.

將式(10)代入式(9)中可得：

(11)

重新修正ZGB空化模型，可以得出考慮熱力學效應的空化模型為

(12)

(13)

式中：RB=1.0×10-6m，經(jīng)驗系數(shù)Fvap和Fcond對熱力學空化流動數(shù)值模擬的結(jié)果產(chǎn)生重要影響.SUN等[12]重新修正了該經(jīng)驗系數(shù)，發(fā)現(xiàn)Fvap=5.0和Fcond=0.001能夠較好地預測溫度和壓強，文中數(shù)值模擬將使用該經(jīng)驗系數(shù).

試驗表明湍流對空化流動產(chǎn)生重要影響，通常使用湍動能k和混合密度ρm來修正飽和蒸汽壓，修正的飽和蒸汽壓公式為

ptur=0.39kρm,

(14)

(15)

式中：pv(T)為飽和蒸汽壓，是溫度的函數(shù).

2 數(shù)值計算與驗證

2.1 繞二維水翼穩(wěn)態(tài)空化流動的數(shù)值驗證

采用文獻[6]在NASA的資助下進行的熱敏空化試驗來驗證數(shù)值計算方法的有效性.試驗模型是一個安裝在截面為方形的水洞中且前緣半徑為3.96 mm的半球頭型水翼.5個溫度傳感器和壓力傳感器分別安裝在水翼的上表面.由于水翼的上下表面對稱，所以只選取其中的一半作為計算域.二維的計算域邊界條件如圖1所示.

圖1 計算域和邊界條件

水翼表面設(shè)置為絕熱，自由滑移.入口的速度和溫度根據(jù)試驗數(shù)據(jù)設(shè)定，出口壓強需要進行多次調(diào)整，使計算得到的入口壓強與試驗數(shù)據(jù)保持一致.對流項設(shè)置為高分辨率格式.

在水翼表面附近進行網(wǎng)格加密，并設(shè)計6套網(wǎng)格驗證網(wǎng)格無關(guān)性，網(wǎng)格1的單元總數(shù)為5 641，網(wǎng)格2的單元總數(shù)為10 205，網(wǎng)格3的單元總數(shù)為25 325，網(wǎng)格4的單元總數(shù)為50 423，網(wǎng)格5的單元總數(shù)為83 505，網(wǎng)格6的單元總數(shù)為128 561.圖2為6套網(wǎng)格的無關(guān)性分析，基于6套網(wǎng)格的計算結(jié)果差別不大，綜合考慮計算精度和計算時間，后續(xù)計算采用網(wǎng)格3(單元總數(shù)為25 325).

圖2 模擬壓力系數(shù)與試驗數(shù)據(jù)對比

文獻[6]試驗中選擇的293F工況作為參考數(shù)據(jù)，根據(jù)試驗參數(shù)，設(shè)置進口速度為23.9 m/s，溫度為77.9 K.圖3為ZGB模型和文中修正的ZGB空化模型(MZGB模型)的計算結(jié)果對比.

圖3 數(shù)值計算得到的壓力和溫度與試驗數(shù)據(jù)對比

由圖3可知，在水翼隨邊附近出現(xiàn)了較大的溫降和壓降.隨著流動往下游移動，由于氣體凝結(jié)的緣故，溫降和壓降逐漸減小.這表明水翼導邊附近氣化吸熱導致了明顯的熱效應出現(xiàn).從圖中對比結(jié)果可知，MZGB計算得到的結(jié)果與試驗數(shù)據(jù)更加接近.隨著入口溫度的降低，空化強度逐漸增加，空腔尾部的閉合區(qū)域沿著水翼表面向下游延伸.空化的傳質(zhì)過程能夠影響局部流場和溫度場.由該圖可以發(fā)現(xiàn)MZGB模型在處理熱敏空化流動方面具有優(yōu)勢，能夠較好地預測水翼表面溫度和壓力的分布.需要注意的是在空腔尾部的閉合區(qū)域預測的壓力與試驗數(shù)據(jù)存在一定程度的差別，這可能是出口處的阻塞效應引起.

2.2 繞三維水翼非穩(wěn)態(tài)空化流動的數(shù)值驗證

熱敏流體空化流動的計算域與文獻[7]的試驗裝置一致，如圖4所示.

圖4 幾何模型及邊界條件

圖5 水翼表面及附近的網(wǎng)格

網(wǎng)格節(jié)點總數(shù)為2 412 720.對流項采用高分辨率格式，瞬態(tài)項采用二階向后歐拉格式.為了平衡計算資源和計算精度，將RMS的殘差設(shè)置為10-4.研究中，將穩(wěn)態(tài)的無空化流動計算結(jié)果作為非穩(wěn)態(tài)空化流動的初始計算文件，時間步長設(shè)置為Δt=Tr/100=6.8×10-5s(Tr=c/u∞).

圖6以氣體體積分數(shù)10%的等值面表示空腔形狀，并選擇4個典型時刻的圖片與Kelly試驗觀察的照片做比較.從中可以看出該模型可以合理地預測到試驗過程中三維水翼附近空腔的脫落及其演化過程.這表明該數(shù)值方法可以模擬出熱力學空化流動特征.

圖6 數(shù)值模擬得到的空腔形狀與試驗對比

圖7 水翼吸力面壓力系數(shù)分布的對比

2.3 熱敏流體空化脫落的動態(tài)特征

由圖6，7可知，文中三維水翼非穩(wěn)態(tài)空化流動的仿真結(jié)果與試驗數(shù)據(jù)基本一致.為了描述熱力學空化流動的動態(tài)特征，圖8顯示了空腔總體積隨時間的變化過程，其計算公式為

(16)

式中：N為計算域中控制單元的總數(shù)目；αi為每個控制單元中氣體的體積分數(shù)；Vi為每個控制單元的體積.

圖8說明了氣體總體積隨時間的變化過程是準周期的.選取一個典型周期中8個時刻來顯示空腔的形成、發(fā)展、脫落和潰滅過程，并以氣體體積分數(shù)為10%的等值面.圖9同時顯示了各時刻中截面的溫度云圖.從圖中可知，空腔區(qū)域的溫度要比周圍液體的溫度低，這是空化時出現(xiàn)的氣化吸熱導致的.若將Tref定義為圖9中典型周期的時間，則可對非穩(wěn)態(tài)熱力學空化流動進行以下描述：① 從0.125Tref到0.375Tref時刻，由于回射流的作用，水翼吸力面上的附著空腔脫落，并被主流輸運到下游變成流動更加紊亂的云空化；② 與此同時，水翼前緣的附著空腔開始形成，并向水翼尾部生長；③ 從0.500Tref到0.625Tref時刻，脫落的云空化潰滅，產(chǎn)生瞬間高壓，并以壓力波的形式向四周傳播；④ 從0.750Tref到1.000Tref時刻，回射流開始形成，并開始誘導附著空腔脫落.

圖8 空腔總體積的變化

該8個時刻基于Ω方法識別的旋渦等值面，其定義式如下

(17)

(18)

(19)

(20)

式中：A為對稱矩陣；B為反對稱矩陣；ζ為一個避免分母為0而設(shè)置的很小的正數(shù).

LIU等[13]認為Ω=0.52的等值面能夠較好地識別旋渦的表面.從圖9中看到Ω準則可以清楚地顯示出許多發(fā)夾渦結(jié)構(gòu)，Ω的等值面與脫落的空腔形狀相似，這表明脫落的空腔區(qū)域中呈現(xiàn)大規(guī)模的旋渦運動，這與早期試驗觀察的結(jié)果一致[14].值得注意的是，在水翼的尾緣附近，Ω的等值面比空腔形狀更復雜.這是由于附著的空腔脫落形成云空化后，隨即被卷向下游.當它被輸運到高壓區(qū)域，其迅速潰滅形成大量的微型氣泡.在此過程中，Ω的等值面不但沒有變小，反而變更大更加復雜了，并且被主流輸運到下游.這說明空腔的脫落和潰滅促進了旋渦結(jié)構(gòu)的形成.

圖9 非定常熱力學空化周期性發(fā)展過程

圖10展現(xiàn)了一個典型周期內(nèi)的水翼截面處預測溫降與B因子預測溫降的對比.

圖10 一個典型周期內(nèi)水翼中截面的溫降

利用B因子對熱敏空化流的熱效應進行分析，其計算公式為

(21)

從圖10中可以看出，B因子能夠有效地預測大部分區(qū)域的溫降，與預測值十分符合.但圖中仍有2個地方存在明顯的誤差，其中一個是B因子無法預測由于液化放熱導致的溫升，已在圖中用綠線標出.水翼前緣由于空腔的生長，存在一個明顯的溫降，但B因子理論過度預測了該溫降，已在圖中用粉線標出.整體上，B因子還是能夠有效地預測空化區(qū)域內(nèi)的溫降.

3 結(jié) 論

1)基于輸運方程的修正空化模型結(jié)合混合湍流模型能夠較合理的預測流動細節(jié)，模擬結(jié)果與試驗數(shù)據(jù)基本吻合.

2)回射流、渦旋與流動分離之間有很強的關(guān)聯(lián)性.黏性和逆壓梯度誘導回射流的形成.當回射流與主流相互作用時，將誘導渦旋的形成和脫落.形成的渦旋反過來又影響流場的結(jié)構(gòu)，并常常伴隨著流動分離的出現(xiàn).

3)Ω的等值面與脫落的空腔形狀相似，這表明脫落的空腔區(qū)域中呈現(xiàn)大規(guī)模的旋渦運動.此外，在水翼的尾緣附近，Ω的等值面比空腔形狀更復雜，這說明空腔的脫落和潰滅促進了旋渦結(jié)構(gòu)的形成.

4)在空化區(qū)域，水翼壁面存在明顯的溫降，B因子能夠有效地預測腔內(nèi)溫降，但仍存在過度預測.