空化模型在低溫流體空化流動三維計算中的應用與評價

孫鐵志，魏英杰，王 聰，路中磊

（1.大連理工大學 船舶工程學院，遼寧 大連 116024；2.哈爾濱工業大學 航天學院，哈爾濱 150001）

0 引 言

當流場中液體的壓強降到當地溫度下的飽和蒸汽壓強以下時，將會發生空化現象[1]。空化多存在于渦輪泵、螺旋槳等流體機械中，空化的發生會導致機械振動、噪聲以及機械性能下降等影響，所以空化問題引起了國內外學者開展相關研究。低溫流體例如液氫、液氧等介質的空化常存在于液體火箭發動機推進系統中，由于低溫流體空化過程中熱力學效應的影響以及低溫環境條件等因素的限制，增加了開展相關研究的難度，但國內外學者一直致力于這方面問題的研究。

Stahl和Stepanoff[2]最早提出了B因子理論，分析了熱力學效應對泵的揚程的影響；后來Sarosdy和Acosta[3]開展了水和氟利昂空化的對比試驗，發現了氟利昂液體空化時空泡呈泡霧狀；二十世紀七十年代Hord[4-5]在NASA的資助下，系統地開展了液氮的低溫空化試驗，從而加深了對低溫流體空化特性的直觀認識。隨著計算機技術的發展，近些年通過數值計算也開展了低溫流體空化特性的研究，Hosangadi等人[6]采用Merkle空化模型對液氫和液氮流體空化流動進行了計算，指出了熱力學效應對空化特性具有顯著影響；Tseng等人[7]對質量傳輸模型的參數進行了敏感性分析。國內的馬相孚、張小斌和黃彪等人[8-10]對低溫流體空化流動問題開展了二維數值模擬研究，初步建立了二維數值計算方法。

本文在計算過程中通過CFX Expression Language語言將不同空化模型嵌入到CFX軟件中，并將液氮、液氫的密度、比熱容、熱傳導系數以及飽和蒸汽壓強等隨溫度變化的物性參數引入到求解代碼中，各個參數隨著溫度的改變實時更新；同時考慮空化過程汽化潛熱的影響，并將其以源相的形式添加到能量方程中，從而形成了一套計算低溫流體空化的三維數值方法。

通過將數值計算結果與試驗數據進行對比和分析，評價不同空化模型的有效性，為今后深入準確全面地預測低溫流體空化特性奠定基礎。

1 數值計算方法

1.1 基本控制方程

對于考慮熱力學效應的低溫液體空化流動問題計算的控制方程，除了連續性方程和動量方程外，還包括包含能量源項的能量方程，依次如下：

其中：ρm=αlρl+αvρv為混合相的密度；下標v和l分別代表汽相和液相；下標i和j代表坐標方向；u為速度；μ為流體的動力粘度，μt是湍流粘度；Keff為熱傳導系數，L為汽化潛熱。

1.2 湍流模型

選取的SST（Shear Stress Transport）湍流模型通過在BSL（Baseline k-ω）湍流模型添加渦粘度限制方程后可有效預測流動分離問題，其中BSL模型方程如下：

其中：μ為流體的動力粘度；μt為湍流粘度；Pkb表示浮力引起的湍動能的產生項；Pk表示由于粘性力引起的湍動能的產生項；常數 β′=0.09、α1=5/9；F1和F2為混合函數；S為應變率的不變測度。

1.3 空化模型

在物理上空化過程是由熱力學和動力學約束的相變過程，通過建立液—汽兩相之間質量傳輸的蒸發源相（）和凝結源相（）的方程來表示兩相之間的轉換過程，通過求解傳輸方程以獲得體積分數（αl）或質量分數（fv）的分布，兩種形式方程如下：

SST模型中添加的渦粘度限制方程形式為：

目前空化模型的建立主要是基于Rayleigh-Plesset汽（氣）泡動力學方程（9）的推導，得到描述汽（氣）泡生長和潰滅的基本方程。

式中：pv為泡內壓強，p為氣泡周圍壓強，R為汽泡半徑，ρl為液體密度，T為液體表面張力。

本文計算采用的空化模型是近些年提出并在常溫流體空化計算中具有廣泛應用的模型，如表1所示。其中Kubota模型可有效解釋渦流和空化氣泡之間粘性效應的相互作用，并重點考慮了空化發生時初生和生長過程空泡體積變化的影響，對于非定常空化流動具有較強適用性；Merkle模型的推導建立過程是用大氣泡群代替單個氣泡，蒸發源相和凝結源相表達式直接與壓力差相關，這區別于Kubota模型壓力差的平方根；Kunz模型在蒸發和凝結過程中起主導因素的分別是壓力差和體積分數。而Kunz模型中的質量傳遞建立是基于兩種不同策略，從液相到汽相的轉換主要取決于空化區域壓力低于飽和蒸汽壓力的差值，而從汽相到液相轉換時，質量傳輸是體積分數的三階多項式，所以體積分數的改變對Kunz模型中影響顯著[11-13]。

表1 不同空化模型方程蒸發源相與凝結源相表達式Tab.1 Evaporation and condensation source term expression of different cavitation models

其中：Cdest和Cprod為蒸發系數和凝結系數的經驗常數；αug為非凝結氣體的體積分數；U∞為無窮遠處流場速度；t∞為特征時間。

2 計算結果與分析

Hord在二十世紀七十年代系統地開展了低溫條件下液氫和液氮空化流動的試驗研究，其獲得的試驗數據被國內外學者廣泛地用于數值計算結果的驗證。本文基于Hord試驗，開展了液氮繞水翼空化流動和尖頂拱在液氫中空化流動的三維計算研究，通過與試驗數據對比，綜合了評價Kubota，Kunz以及Merkle三種空化模型在不同低溫介質、不同幾何模型條件下空化流動計算中的適用性。

2.1 液氮繞水翼空化流動計算結果與分析

2.1.1 計算模型和網格劃分

液氮繞水翼空化計算模型與Hord試驗一致，試件結構及尺寸如圖1所示，根據試驗條件，計算域入口采用速度入口，出口為壓力開口，流域邊界及水翼設置為絕熱不可滑移壁面，并在翼型表面分別設置五個溫度監測點和五個壓力監測點，以便將數值計算結果與試驗數據進行對比。圖2給出了計算模型網格劃分，計算域內采用H-型和C-型網格，同時在水翼附近網格進行加密處理，以提高計算效率，有效捕捉翼型周圍流場參數，且在計算過程中對網格無關性進行了驗證。

圖1 計算模型結構Fig.1 Structure of computational model

圖2 計算域模型網格劃分Fig.2 Grids of computational model

為有效評價不同空化模型的適用性，根據試驗數據隨機選取三種不同工況，如表2所示。計算過程中保持計算網格、湍流模型、邊界條件以及試驗參數等條件一致。通過對比不同空化模型下空化區域兩相分布特性、翼型周圍壓力和溫度數據等綜合因素，對三種不同空化模型的適用性進行評價。

表2 計算工況（液氮）Tab.2 Computational cases(Liqiud nitrogen)

2.1.2 空化區域相分布特性

為分析不同空化模型下液氮繞水翼空化流場兩相分布特性，首先給出290C工況三種模型的空化特性對比，如圖1所示。從圖1（a）中三維空泡形態的對比可以看出，空化發生起于翼型頭部，并沿翼型向下游發展，發展過程中空泡厚度增加，且在空泡閉合位置出現回射流。從圖中也可以看出不同空化模型下空化特性的差異：Kubota空化模型計算得到的空泡厚度最大，并且最大厚度位置接近空泡尾部；Merkle模型計算的空泡形態類似于橢球形，整體汽—液交界面輪廓過渡圓滑；而Kunz模型計算的空泡在尾部閉合區域表現尖銳，回射現象更加明顯，并且空泡厚度最小。同時Kubota模型計算得到的空泡長度為23.0 mm，Merkle模型為23.6 mm，Kunz模型為27.4 mm，試驗值為22.96 mm。雖然三種模型計算的結果都與試驗值接近，但Kubota模型計算的空泡長度與試驗結果吻合最好。

圖3 290C工況空化區域相分布特性對比Fig.3 Phase distribution of 290C case at the cavitation region

2.1.3 溫度場和壓力場分布對比

為進一步對比三種空化模型計算結果的差異和驗證數值計算結果的準確性，圖4給出了三種空化模型在不同工況下數值計算結果與Hord試驗數據的對比，計算結果與試驗數據具有較好的一致性，從而驗證了計算方法的有效性。對比不同空化模型的計算結果可以看出，不同工況下同一空化模型計算的壓力和溫度整體變化趨勢一致，在空化區域壓力和溫度值降低，且壓力和溫度值并不保持恒定，這主要是由于液氮等低溫液體空化時由于熱力學效應會引起空化區域溫度降低，而液氮飽和蒸汽壓強是溫度的函數，且隨溫度變化敏感，從而導致空化區域壓強值呈梯度變化，這一特性與常溫水有明顯區別。三種空化模型計算的空化區域的最大壓降、最大溫降有差異，且在空化區域壓力和溫度曲線變化的梯度不同。壓力方面：在空化區域Kubota模型計算得到的最大壓降（①位置）最大，在空化區域壓力值相對穩定，恢復到無窮遠處值壓力曲線起點（②位置）滯后，且變化梯度最大，Kubota模型計算的空泡閉合區域（③位置）會有壓力峰值；Kunz模型計算的壓力曲線梯度總體最小，Merkle模型計算的壓力梯度值介于兩者之間；溫度方面：總體而言Kubota空化模型計算得到的溫度曲線變化梯度最大，空化區域最大溫降值較大，Kunz模型得到的溫度曲線變化最緩慢，且最大溫降最小。綜合對比可以看出Kubota模型對液氮繞水翼空化流動預測較好，Merkle模型可較好地反映出空化區域溫度場分布。

圖4 三種空化模型下壓力和溫度計算結果與試驗對比Fig.4 Comparison of the numerical results with experimental data by three cavitation models

2.2 液氫空化流動計算結果與分析

液氫空化計算模型與Hord[5]試驗一致，如圖5所示。計算幾何模型為尖頂拱形回轉體，且在尾部擴張，計算域上游為圓柱形回轉體，并在尖頂拱肩部位置開始擴張，如圖5（a）所示。計算模型具體尺寸如圖5（b）所示，其中尖形拱直徑為9.1 mm，肩部倒角半徑為2.3 mm。計算模型網格劃分方法與2.1節中水翼模型一致，在尖形拱頭部空化區域網格精細劃分以有效捕獲空化流動特性。

圖5 液氫空化計算模型Fig.5 Computational model of cavitating flow in liquid hydrogen

同樣，根據Hord試驗隨機選取三種試驗工況作為計算對比驗證，如表3所示。計算邊界條件與

2.1 節中液氮繞水翼空化設置一致，尖形拱和水洞邊界都設置為絕熱壁面。

圖6給出了三種空化模型下尖頂拱在液氫中空化流動數值計算結果，沿尖頂拱周圍沿軸向方向分布的壓力和溫度變化趨勢與2.1節中液氮繞水翼空化結果整體變化趨勢一致，數值計算結果與試驗數據吻合較好。從圖6（a）中可以看出，三種空化模型計算的347B和366B工況壓力結果差異不明顯，在513B工況中Kubota模型計算的壓力結果與試驗數據最接近，其次是Merkle模型。從圖6（b）給出的三種空化模型計算得到的溫度變化與試驗對比可以看出，Kubota空化模型計算得到的最大溫降值最高，Kubota和Merkle兩種空化模型對溫度場的預測要優于Kunz空化模型。綜合對比可以看出Kubota和Merkle模型對尖頂拱在液氫中空化流動的溫度變化預測較準確，Kubota模型可較好地計算出尖頂拱周圍壓力場分布。

表3 計算工況（液氫）Tab.3 Computational cases(Liqiud hyrogen)

圖6 液氫中三種空化模型計算的壓力和溫度結果與試驗數據對比Fig.6 Computed pressure and temperature profiles in liquid hydrogen compared with experimental data

2.3 不同空化模型機理分析與評價

為說明不同空化模型得到的空化特性、壓力分布和溫度降低的差別，以290C（液氮）和366B（液氫）工況為例，從兩種流體介質不同空化模型下空化區域蒸發源相和凝結源相分布特性進行分析，如圖7所示。從圖中蒸發源相分布對比可以看出，兩種介質中分布趨勢一致，其中Kubota模型的值最小，Merkle模型和Kunz模型計算的結果整體接近；從表1中不同空化模型源相公式可知，蒸發源相值主要取決于液相體積分數αl和pv-p的大小，Kubota空化模型中蒸發源相是pv-p的二分之一次方，所以Kubota模型需要較大的壓降以驅動相同的空化強度。

圖7 三種空化模型下蒸發源相和凝結源相對比Fig.7 Comparison of the evaporation and condensation source term by three cavitation models

3 結 論

為評價Kubota、Merkle和Kunz三種空化模型對低溫流體空化過程流場特性預測的適用性，通過對CFX軟件的二次開發，將空化模型、能量源相和液氮、液氫物性參數嵌入到CFX求解代碼中，建立了考慮熱力學效應下液氮繞水翼空化流動和尖頂拱在液氫中空化流動三維數值的計算方法，通過將計算結果與試驗數據對比分析，得到主要結論如下：

（1）不同空化模型計算得到的空化區域空化特性差異明顯，Kubota模型計算的空泡厚度最大、液相體積分數沿徑向變化緩慢，Merkle模型次之；Kunz模型計算的空泡長度最大，空化閉合區域回射流現象明顯。

（2）空化模型質量傳輸機制顯著影響著空化區域流場特性，空化區域壓差、液/汽相體積分數值以及流體介質密度是決定著空化模型方程中蒸發源相和凝結源相質量傳輸大小的主要因素，進而影響著空化流場特性，從而導致Kubota模型計算的空化區域壓降較大，空泡閉合位置壓力變化梯度最大，溫度變化迅速；而Kunz模型在空化區域從低壓恢復到無窮遠處壓強緩慢。

（3）三種空化模型對低溫流體空化流動流場結果預測存在差異，Kubota空化模型可有效計算流場壓力分布，Merkle模型可較好地反映熱力學效應下空化區域溫度降，而Kunz模型計算的壓力、溫度結果與試驗數據偏差最大，總體上Kubota模型和Merkle模型適用性較強。

[1]Joseph D D.Cavitation in a flowing liquid[J].Physical Review E,1995,51(3):R1649.

[2]Stahl H A,Stepanoff A J.Thermodynamic aspects of cavitation in centrifugal pumps[J].Trans.ASME,1956,78(8):1691-1693.

[3]Sarosdy L R,Acosta A J.Note on observations of cavitation in different fluids[J].Journal of Basic Engineering,1961,83:399-400.

[4]Hord J.Cavitation in liquid cryogens II-Hydrofoil,1973[R].NASA CR-2156.

[5]Hord J.Cavitation in liquid cryogens III-Ogives,1973[R].NASA CR-2242.

[6]Hosangadi A,Ahuja V.Numerical study of cavitation in cryogenic fluids[J].Journal of Fluids Engineering,2005,127(2):267-281.

[7]Tseng C C,Shyy W.Modeling for isothermal and cryogenic cavitation[J].International Journal of Heat and Mass Transfer,2010,53(1):513-525.

[8]Ma Xiangfu,Wei Yingjie,Wang Cong,et al.Simulation of liquid nitrogen around hydrofoil cavitation flow[J].Journal of Ship Mechanics,2012(12):1345-1352.

[9]Zhang X B,Wu Z,Xiang S J,et al.Modeling cavitation flow of cryogenic fluids with thermodynamic phase-change theory[J].Chinese Science Bulletin,2013,58(4-5):567-574.

[10]Huang B,Wu Q,Wang G.Numerical investigation of cavitating flow in liquid hydrogen[J].International Journal of Hydrogen Energy,2014,39(4):1698-1709.

[11]Kubota A,Kato H,Yamaguchi H,et al.Unsteady structure measurement of cloud cavitation on a foil section using conditional sampling technique[J].Journal of Fluids Engineering,1989,111(2):204-210.

[12]Merkle C L,Feng J,Buelow P E O.Computational modeling of dynamics of sheet cavitation[C]//Proceedings of the 3rd International Symposium on Cavitation.Grenoble,France,1998.

[13]Kunz R F,Boger D A,Stinebring D R,et al.A preconditioned Navier-Stokes method for two-phase flows with application to cavitation prediction[J].Computers&Fluids,2000,29(8):849-875.