空化模型在低溫流體空化流動(dòng)三維計(jì)算中的應(yīng)用與評(píng)價(jià)

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

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

0 引 言

當(dāng)流場(chǎng)中液體的壓強(qiáng)降到當(dāng)?shù)販囟认碌娘柡驼羝麎簭?qiáng)以下時(shí)，將會(huì)發(fā)生空化現(xiàn)象[1]?？栈啻嬖谟跍u輪泵、螺旋槳等流體機(jī)械中，空化的發(fā)生會(huì)導(dǎo)致機(jī)械振動(dòng)、噪聲以及機(jī)械性能下降等影響，所以空化問題引起了國內(nèi)外學(xué)者開展相關(guān)研究。低溫流體例如液氫、液氧等介質(zhì)的空化常存在于液體火箭發(fā)動(dòng)機(jī)推進(jìn)系統(tǒng)中，由于低溫流體空化過程中熱力學(xué)效應(yīng)的影響以及低溫環(huán)境條件等因素的限制，增加了開展相關(guān)研究的難度，但國內(nèi)外學(xué)者一直致力于這方面問題的研究。

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

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

通過將數(shù)值計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比和分析，評(píng)價(jià)不同空化模型的有效性，為今后深入準(zhǔn)確全面地預(yù)測(cè)低溫流體空化特性奠定基礎(chǔ)。

1 數(shù)值計(jì)算方法

1.1 基本控制方程

對(duì)于考慮熱力學(xué)效應(yīng)的低溫液體空化流動(dòng)問題計(jì)算的控制方程，除了連續(xù)性方程和動(dòng)量方程外，還包括包含能量源項(xiàng)的能量方程，依次如下：

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

1.2 湍流模型

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

其中：μ為流體的動(dòng)力粘度；μt為湍流粘度；Pkb表示浮力引起的湍動(dòng)能的產(chǎn)生項(xiàng)；Pk表示由于粘性力引起的湍動(dòng)能的產(chǎn)生項(xiàng)；常數(shù) β′=0.09、α1=5/9；F1和F2為混合函數(shù)；S為應(yīng)變率的不變測(cè)度。

1.3 空化模型

在物理上空化過程是由熱力學(xué)和動(dòng)力學(xué)約束的相變過程，通過建立液—汽兩相之間質(zhì)量傳輸?shù)恼舭l(fā)源相（）和凝結(jié)源相（）的方程來表示兩相之間的轉(zhuǎn)換過程，通過求解傳輸方程以獲得體積分?jǐn)?shù)（αl）或質(zhì)量分?jǐn)?shù)（fv）的分布，兩種形式方程如下：

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

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

式中：pv為泡內(nèi)壓強(qiáng)，p為氣泡周圍壓強(qiáng)，R為汽泡半徑，ρl為液體密度，T為液體表面張力。

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

表1 不同空化模型方程蒸發(fā)源相與凝結(jié)源相表達(dá)式Tab.1 Evaporation and condensation source term expression of different cavitation models

其中：Cdest和Cprod為蒸發(fā)系數(shù)和凝結(jié)系數(shù)的經(jīng)驗(yàn)常數(shù)；αug為非凝結(jié)氣體的體積分?jǐn)?shù)；U∞為無窮遠(yuǎn)處流場(chǎng)速度；t∞為特征時(shí)間。

2 計(jì)算結(jié)果與分析

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

2.1 液氮繞水翼空化流動(dòng)計(jì)算結(jié)果與分析

2.1.1 計(jì)算模型和網(wǎng)格劃分

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

圖1 計(jì)算模型結(jié)構(gòu)Fig.1 Structure of computational model

圖2 計(jì)算域模型網(wǎng)格劃分Fig.2 Grids of computational model

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

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

2.1.2 空化區(qū)域相分布特性

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

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

2.1.3 溫度場(chǎng)和壓力場(chǎng)分布對(duì)比

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

圖4 三種空化模型下壓力和溫度計(jì)算結(jié)果與試驗(yàn)對(duì)比Fig.4 Comparison of the numerical results with experimental data by three cavitation models

2.2 液氫空化流動(dòng)計(jì)算結(jié)果與分析

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

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

同樣，根據(jù)Hord試驗(yàn)隨機(jī)選取三種試驗(yàn)工況作為計(jì)算對(duì)比驗(yàn)證，如表3所示。計(jì)算邊界條件與

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

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

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

圖6 液氫中三種空化模型計(jì)算的壓力和溫度結(jié)果與試驗(yàn)數(shù)據(jù)對(duì)比Fig.6 Computed pressure and temperature profiles in liquid hydrogen compared with experimental data

2.3 不同空化模型機(jī)理分析與評(píng)價(jià)

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

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

3 結(jié) 論

為評(píng)價(jià)Kubota、Merkle和Kunz三種空化模型對(duì)低溫流體空化過程流場(chǎng)特性預(yù)測(cè)的適用性，通過對(duì)CFX軟件的二次開發(fā)，將空化模型、能量源相和液氮、液氫物性參數(shù)嵌入到CFX求解代碼中，建立了考慮熱力學(xué)效應(yīng)下液氮繞水翼空化流動(dòng)和尖頂拱在液氫中空化流動(dòng)三維數(shù)值的計(jì)算方法，通過將計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)對(duì)比分析，得到主要結(jié)論如下：

（1）不同空化模型計(jì)算得到的空化區(qū)域空化特性差異明顯，Kubota模型計(jì)算的空泡厚度最大、液相體積分?jǐn)?shù)沿徑向變化緩慢，Merkle模型次之；Kunz模型計(jì)算的空泡長(zhǎng)度最大，空化閉合區(qū)域回射流現(xiàn)象明顯。

（2）空化模型質(zhì)量傳輸機(jī)制顯著影響著空化區(qū)域流場(chǎng)特性，空化區(qū)域壓差、液/汽相體積分?jǐn)?shù)值以及流體介質(zhì)密度是決定著空化模型方程中蒸發(fā)源相和凝結(jié)源相質(zhì)量傳輸大小的主要因素，進(jìn)而影響著空化流場(chǎng)特性，從而導(dǎo)致Kubota模型計(jì)算的空化區(qū)域壓降較大，空泡閉合位置壓力變化梯度最大，溫度變化迅速；而Kunz模型在空化區(qū)域從低壓恢復(fù)到無窮遠(yuǎn)處壓強(qiáng)緩慢。

（3）三種空化模型對(duì)低溫流體空化流動(dòng)流場(chǎng)結(jié)果預(yù)測(cè)存在差異，Kubota空化模型可有效計(jì)算流場(chǎng)壓力分布，Merkle模型可較好地反映熱力學(xué)效應(yīng)下空化區(qū)域溫度降，而Kunz模型計(jì)算的壓力、溫度結(jié)果與試驗(yàn)數(shù)據(jù)偏差最大，總體上Kubota模型和Merkle模型適用性較強(qiáng)。

[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.