劉宏升，馬杰，李亮，吳丹，解茂昭

(1. 大連理工大學海洋能源利用與節(jié)能教育部重點實驗室，遼寧大連，116024；2. 濰柴動力股份有限公司，山東濰坊，261000)

面對節(jié)能與環(huán)保的雙重挑戰(zhàn)，目前許多能量轉換裝置的工作壓力已經(jīng)超過了所涉及流體的臨界壓力，例如，液體火箭發(fā)動機、汽輪機和壓燃式發(fā)動機等，這使噴霧燃燒過程常常在接近或超過燃料臨界壓力的環(huán)境下進行[1]。

在發(fā)動機領域，高壓液體噴射系統(tǒng)通常存在2種極端情況：一是在亞臨界環(huán)境下，射流存在一個明確的氣液分界面，該界面在氣動不穩(wěn)定性作用下會發(fā)生初次及二次破碎霧化，形成噴霧；二是在跨/超臨界環(huán)境下，伴隨高梯度傳熱傳質和臨界相變等復雜工況的出現(xiàn)，氣液界面將不復存在，表面張力趨近于零，因缺乏分子內(nèi)聚力導致液滴形成的可能性減小甚至完全消除，此時燃料噴霧混合過程由分子和湍流擴散所主導，燃料變成一種介于液體與氣體之間的高密度流體，即所謂的超臨界流體[2-3]?？梢姡谶@種環(huán)境下，燃料的噴射與混合氣的形成機理極其復雜，與亞臨界工況存在明顯差異。因此，發(fā)動機領域的燃油噴霧混合必然涉及跨/超臨界工況，掌握跨/超臨界噴射混合的機理及其控制是至關重要的[4]。鑒于此，眾多學者開始對跨/超臨界燃料噴霧混合過程進行實驗和數(shù)值模擬方面的研究。

MAYER 等[5]通過液氧/氣氫的射流實驗發(fā)現(xiàn)，超臨界環(huán)境下的射流表面與亞臨界環(huán)境下的情形截然不同，具體表現(xiàn)為從液核分離出極細的流體線，進而發(fā)展成各種形狀的高密度氣體團。最近幾年，出現(xiàn)了以內(nèi)燃機為背景的超臨界噴霧實驗研究。MANIN 等[6]研究了正十二烷噴霧在噴嘴附近的形態(tài)變化，發(fā)現(xiàn)隨缸內(nèi)壓力和溫度升高，傳統(tǒng)噴霧中陡峭的兩相界面轉變?yōu)闆]有明確液滴的擴散層。然而，由于缸內(nèi)氣體劇烈密度波動的影響，實驗的圖像質量和加速環(huán)境受到了很大限制。POURSADEGH 等[7]采用燃油直噴技術將丙烷噴射到氣態(tài)氮內(nèi)進行超臨界噴霧實驗研究，他們將液滴形成時間和液滴蒸發(fā)時間的比值定義為量綱一的時間τ，以此來分析燃料液滴出現(xiàn)的條件，發(fā)現(xiàn)在較低溫度下，當τ＜1時為傳統(tǒng)的噴霧結構，有液滴存在；當τ＞1時噴霧呈稠密的混合形態(tài)，沒有明顯液滴存在。

數(shù)值模擬方面，SELLE 等[8]利用立方型PR 狀態(tài)方程對MAYER 的實驗進行大渦模擬(LES)，跨臨界射流下他們的研究結果與實驗值吻合良好；KIM 等[9]采用RANS 方法對比研究了PR 和SRK 這2個立方型狀態(tài)方程，他們發(fā)現(xiàn)狀態(tài)方程的選擇不僅會影響熱力學狀態(tài)變量，還會影響湍流強度。然而2種狀態(tài)方程在超臨界區(qū)域都表現(xiàn)出不可忽視的誤差；基于LES方法，ZONG等[10]對液氮的超臨界混合進行了模擬，結果表明射流表面存在一個高密度梯度區(qū)域，具有類似于固體壁面的作用，能夠有效抑制噴霧表面正常的速度波動，從而減緩了噴霧的破碎過程；MüLLER 等[11]對超臨界壓力下液氮和預熱氫的同軸噴射進行了大渦模擬，介紹了基于立方型PR 狀態(tài)方程的體積平移方法，用于多組分流體的大渦模擬，并對跨臨界和超臨界射流條件進行了測試；PARK 等[12-13]對MAYER等的氮氣噴霧實驗進行了模擬，對比研究了多個SGS 模型 ， 發(fā)現(xiàn) VREMAN 模型[14]和SMAGORINSKY 模型[15]都得到比較理想的結果；解茂昭等[4,16]提出一種跨臨界/超臨界環(huán)境下液體燃料噴射與混合氣形成的數(shù)理模型，即“混合層回縮”模型，他們采用LES 方法，對跨臨界和超臨界條件下的低溫液氮射流進行了數(shù)值研究，并討論了偽沸騰現(xiàn)象對射流演化過程的影響。

基于以上研究，本文作者采用自主構建的超臨界射流模型對低溫液氮噴入高溫氮氣環(huán)境進行RANS 模擬，研究和分析跨/超臨界射流的特點及混合過程。

1 理論模型

針對超臨界流體的流動特性及其在臨界點附近的熱物理性質的非理想性和輸運屬性的瞬變特性，本文結合真實流體狀態(tài)方程、熱力學和CFD的相關理論構建了一個新的CFD求解程序。其中，采用RANS 模型對湍流進行封閉，對經(jīng)典的PISO算法進行修正，以實現(xiàn)壓力和速度的耦合迭代計算。

1.1 控制方程

式中：ρ為流體密度；u為速度；h為焓；p為壓力；μ為動力黏度；D 為物質導數(shù)符號；下標i，j和k代表笛卡爾坐標方向；αeff為有效熱擴散系數(shù)；τij為黏性應力張量，其與流體的應變率張量Sij呈線性關系：

式中：Skk為流體散度divu；δij為二階單位張量，當i=j時，δij=1，否則δij=0。

1.2 真實氣體狀態(tài)方程

在跨/超臨界條件下，分子間的作用力不能忽略，理想氣體狀態(tài)方程不再適用，因此需要采用真實氣體狀態(tài)方程來分析流體的熱物理性質?？紤]到立方型狀態(tài)方程在保證高計算效率的前提下，仍能得到較好的預測結果。本研究中選用PR狀態(tài)方程來模擬跨/超臨界射流噴霧，其具體形式如下：

式中：T為溫度；p為壓力；Tc為臨界溫度；pc為臨界壓力；ρ為密度；Ru為通用氣體常數(shù)；Mw為物質的量；a和b分別表示分子間引力和斥力的作用；α(T)為交叉系數(shù)；SPR為中間計算變量；ω為偏心因子。

1.3 熱力學和輸運屬性

1.3.1 熱力學屬性

在跨臨界和超臨界條件下，單組分流體熱力學屬性，例如焓h(T,p)和質量定壓熱容cp(T,ρ)，可由理想狀態(tài)的參考值(p0=0.1 MPa)以及偏離函數(shù)共同決定，如下式所示：

其中，以PR 方程為例的壓力偏導數(shù)項的推導如下：

式中：cp0(T)和h0(T)均為理想狀態(tài)參考值；可由NASA多項式方程求解。

式中：a1～a6為多項式系數(shù)。

1.3.2 輸運屬性

單組分流體的運動黏度可以通過摩擦理論求解。本研究中選用的組分為氮氣，下面以氮氣為例進行描述。

根據(jù)摩擦理論[17]，總黏度η包括稀氣體黏度η0和殘余摩擦ηf：

式中：稀氣體黏度η0定義為零密度極限下的黏度，而殘余摩擦項ηf與經(jīng)典力學的摩擦概念有關。根據(jù)摩擦理論，殘余摩擦項可以表示如下：

稀氣體黏度可表示為

式中：d1，d2，d3，κr，κa和κrr為模型系數(shù)；pr和pa分別為范德華吸引力和排斥力項，根據(jù)PR方程推出的表達式如下：

氮氣的詳細模型系數(shù)如表1所示。

表1 摩擦理論中氮氣模型系數(shù)Table 1 Nitrogen coefficient in friction theory

單組分導熱系數(shù)通過VASSERMAN 等[18]推出的標準方程進行求解：

式中：ζ和τ分別為對比密度和對比溫度，即

式(23)中相關系數(shù)的計算公式如下:

1.4 湍流模型

眾所周知，基于渦流黏度概念，學者們提出了諸多RANS 模型，其中求解湍流動能(k)及其耗散率(ε)的k-ε模型應用最為廣泛。根據(jù)KIM等[9]和PARK 等[12]對跨/超臨界射流模擬的結果，湍流模型對射流特性的影響并不像狀態(tài)方程那樣重要。因此，為了降低計算成本，本研究使用了RANS模型中的可壓縮k-ε模型來封閉湍流問題，該模型的控制方程如下：

式中：k和ε分別為湍流動能和湍流能量的耗散率。紊流能量的產(chǎn)生速率，Pk表示為

以符合Boussinesq 渦黏度近似的方式評估Pk，本研究中使用的模型常數(shù)如表2所示。

表2 k-ε湍流模型系數(shù)Table 2 Coefficient in k-ε turbulence model

2 數(shù)值方法和計算模型

根據(jù)MAYER 的跨/超臨界實驗[19]，液氮射流從直徑為2.2 mm且長徑比大于40的完全發(fā)展湍流管中噴出，從而確保了出口處湍流的充分發(fā)展。實驗中采用的是一個截面邊長為60 mm 的方管，長度約為1 m，并在超臨界壓力下充滿了高溫氣態(tài)氮氣。由于計算腔室內(nèi)的整個流場需要大量的計算成本，在當前的模擬研究中，為了降低計算消耗，在OpenFOAM 中采用具有周向圓周角為5°的2維軸對稱結構，噴嘴內(nèi)徑為1.1 mm，混合腔室管徑為61 mm、長度為350 mm，以確保出口處湍流的充分發(fā)展。計算域均采用六面體結構網(wǎng)格，噴嘴內(nèi)徑使用22 個網(wǎng)格點離散，腔室采用350×160個網(wǎng)格點離散。根據(jù)PARK等[12,20]的網(wǎng)格無關性驗證，當前網(wǎng)格系統(tǒng)足以滿足網(wǎng)格獨立性。

本文模型在一系列液氮的跨/超臨界實驗中選取3個算例(Case)進行研究，環(huán)境壓力分別為4 MPa和6 MPa，均高于氮氣的臨界壓力(pc=3.34 MPa)，環(huán)境溫度均為298 K。三者的入射溫度都高于氮氣的臨界溫度(Tc=126.2 K)，Case 1 和Case 3 均低于各自環(huán)境壓力下的偽沸點溫度Tpb(Case 1 的射流入口初始溫度Tinj=126.9 K，Tpb=129.8 K(4 MPa)；Case 3 的Tinj=135.4 K，Tpb=139.1 K(6 MPa))。表3所示為3 個Case 的邊界條件和內(nèi)部場初始條件。其中：vinj和Tinj分別為射流入口初始速度和溫度；p∞為腔室環(huán)境初始壓力；T∞為腔室環(huán)境初始溫度；Re為雷諾數(shù)。

表3 跨/超臨界射流初始及邊界條件Table 3 Initial and boundary conditions of trans/supercritical jets

3 結果與討論

由表3可知：Case 1和Case 3的射流中將會跨越比熱容峰值點而出現(xiàn)偽沸騰現(xiàn)象，因此，稱為跨臨界射流。Case 2 的射流入口初始溫度為137.0 K，高于4 MPa下的臨界溫度和偽沸騰溫度，稱為嚴格意義上的超臨界射流。LI 等[21]發(fā)現(xiàn)：在初始速度相差較小時，跨臨界射流與超臨界射流僅在湍流屬性上有細微差別。

為了對跨/超臨界射流模型的計算結果進行驗證，將模擬結果與實驗數(shù)據(jù)進行對比。圖1～2所示分別為Case 1和Case 2工況下模擬與實驗的軸心線密度對比，x/d為軸心距離與直徑的比值。由圖1～2可以看出：在軸心線上流體密度分布的模擬值與實驗結果基本一致。對于Case 1 而言，在x/d＜7的區(qū)域，PR 狀態(tài)方程的模擬結果(PR EOS)過高地預測了噴嘴附近的高密度分布，而當x/d＞7時，特別是噴嘴下游區(qū)域，當前模型結果與實驗值吻合較好，沿軸心線的密度衰減趨勢與實驗完全一致，證明模型具有較高的預測精度。對于Case 2 工況下的超臨界射流，無論是在噴嘴附近的高密度區(qū)還是在遠離噴嘴的下游區(qū)域，模型的模擬結果均與實驗值吻合良好。對比2種工況的軸線密度分布還可以看出，噴嘴附近(x/d＜7)存在明顯的高密度核心，可以視為稠密液核區(qū)，之后，由于噴霧的擴散以及周圍環(huán)境中高溫氣體的加熱作用使得密度迅速下降。在密度下降區(qū)，當前模型的結果則顯示了更高的精度，原因在于采用真實氣體狀態(tài)方程對壓力泊松方程進行修正，考慮了壓力波動的影響。

圖3所示為Case 1工況下軸距分別為5d和25d處密度的徑向分布情況，并與理想氣體狀態(tài)方程的計算結果(Ideal EOS)以及實驗值進行了對比。由圖3可以看出：在x/d=5處，本模型中使用PR狀態(tài)方程的計算結果與實驗值基本吻合，但仍有一些偏差?？拷S心區(qū)域(徑向距離與直徑的比值r/d＜0.5)，本模型的模擬值預測的密度峰值偏低；然而，在r/d＞0.5 的區(qū)域，模擬值與實驗值吻合更好，特別是在遠離軸心區(qū)(r/d＞1)的位置，兩者完全吻合。KIM等[9]的研究證實：流場中的湍流擾動可導致密度峰值的偏差，同時光學測量中的反射與折射也會引起這種偏差。在x/d=25 處，靠近射流中心線處的模擬值與實驗值吻合較好，準確地預測了密度峰值；然而在遠離軸心區(qū)的位置，本模型預測的密度衰減較實驗而言相對緩慢，但整體的密度衰減趨勢與實驗一致。同時可以看出，無論在x/d=5處，還是x/d=25處，利用理想氣體狀態(tài)方程構建模型的模擬結果預測的核心區(qū)域的密度分布均明顯偏低。

可見，利用PR 狀態(tài)方程構建的新模型在噴嘴附近及射流下游的預測能力較好，噴嘴附近的恒值密度長度與實驗也基本相符。尤其在超臨界射流中，噴嘴附近表現(xiàn)為射流液體進入腔室后密度立即下降。在跨/超臨界條件下，無法觀察到明顯的氣液界面，反而表現(xiàn)為一種拓寬的氣-液混合層。為理解跨/超臨界的射流機理，有必要對其混合熱動力學屬性及混合層進行分析。

圖4 所示為超臨界工況Case 2 和跨臨界工況Case 3在5個不同軸距處密度的徑向分布情況。由圖4 可以看出：跨/超臨界的混合層中均存在一個密度峰值。對比相同位置處，可以發(fā)現(xiàn)Case 3 的密度峰值遠比Case 2的高，并且Case 3的密度峰值沿軸向的下降趨勢明顯比Case 2 的快，但跨臨界的密度沿徑向的擴散速度卻明顯比超臨界射流的速度低。其原因在于：在跨臨界射流中，偽沸騰發(fā)生之前，低溫液核周圍的氣液混合層內(nèi)具有高的密度梯度分層，類似于“固體壁效應”，從而使得低溫液氮在軸線方向擴散更遠，徑向擴散受到明顯抑制；雖然Case 2 的環(huán)境壓力(4 MPa)低于Case 3(6 MPa)，但因為其初始噴射溫度高于偽沸騰溫度，射流的初始狀態(tài)即為超臨界態(tài)，因此，射流直接進入混合狀態(tài)，并且混合層中的密度梯度小，使超臨界射流產(chǎn)生更寬的混合層。

對于超臨界工況Case 2，當軸距從5d增至10d時，密度迅速下降并在15d之后曲線接近水平。這也再次證明射流流體與環(huán)境流體之間的密度差異已十分有限，流體性質因湍流混合而趨于一致。反之，跨臨界射流Case 3在20d時徑向密度依然存在較大差異。這表明超臨界射流比跨臨界射流更早進入混合自相似狀態(tài)[22]。

圖5所示為3種工況下的質量定壓熱容和溫度沿軸心線的分布情況。由圖5可看出：對于跨臨界射流Case 1 而言，質量定壓熱容在低溫液核末端(7d)開始迅速增大，并在12d處達到峰值，即跨臨界射流中存在明顯的質量定壓熱容峰值。但從Case 1的溫度曲線可以看出，溫度在該范圍內(nèi)并沒有顯著增加。而對于跨臨界射流Case 3 而言，環(huán)境壓力增大反而導致質量定壓熱容在偽臨界點附近的值減小，且溫升速率明顯增加，這是因為氮氣的臨界壓力為3.34 MPa，當環(huán)境壓力與其接近時，低溫射流從高溫氣體中吸收的能量主要用于跨越偽沸騰線引起的體積膨脹，而非提高自身溫度，同時也進一步解釋了Case 1 中密度曲線的迅速下降與偽沸騰現(xiàn)象有關；而在環(huán)境壓力更高的跨臨界射流中，其質量定壓熱容峰值會明顯下降，即偽沸騰現(xiàn)象顯著減弱。對于超臨界射流Case 2，由于不存在偽沸騰現(xiàn)象，其與跨臨界射流相比，不存在質量定壓熱容峰值，且溫升起點更早，溫升速率更高。即由于不存在偽沸騰現(xiàn)象，從高溫流體吸收的熱量并不會導致明顯的體積膨脹，而是完全用于自身升溫。

圖6～7所示分別為在x/d＜30區(qū)域內(nèi)，Case 1和Case 2工況下射流的質量定壓熱容、密度、溫度及湍流動能的瞬時云圖。由圖6和圖7可知跨臨界射流與超臨界射流的射流輪廓存在顯著差異。

由圖6可以看出：對于跨臨界射流Case 1，由于跨臨界射流穿越偽沸點導致出現(xiàn)比熱峰值，從而形成一個過渡混合層，且比熱峰值等值面可視為氣液界面。然而，對于超臨界射流Case 2，質量定壓熱容在射流進入腔室后便不斷減小。對于密度而言，跨臨界射流存在更長的低溫液核區(qū)。超臨界工況下，盡管噴嘴附近存在較短的液核，但密度值較低，且不同于亞臨界射流，是一種超臨界狀態(tài)。從Case 2 下游區(qū)域可見，超臨界射流的密度擴散區(qū)更寬。

由圖7可以看出：跨臨界射流Case 1和超臨界射流Case 2 在軸向和徑向的溫度分布上，跨臨界射流比超臨界射流具有更寬的徑向分布，同時軸向同長度范圍內(nèi)溫升十分有限，而超臨界射流的溫度在軸向呈現(xiàn)連續(xù)梯度變化，溫升速率明顯加快。這是由于超臨界工況下不存在偽沸騰，射流從周圍高溫流體中吸收的熱量主要用于提高溫度，使其溫升速率更快。在湍動能方面，跨臨界射流的湍流混合主要發(fā)生在下游區(qū)域，且分布范圍更廣。在跨臨界工況下，存在一包含了高密度射流核高比熱層，該層延遲了噴霧的非穩(wěn)態(tài)發(fā)展，抑制了徑向波動，導致湍流動能從流動徑向向軸向重新分配，從而使湍流動能從徑向向軸向轉移。

4 結論

1)對于液氮的跨/超臨界射流噴射，利用PR狀態(tài)方程構建的新模型在噴嘴附近及噴嘴下游都表現(xiàn)出良好的預測能力，噴嘴附近的致密液核長度與實驗也基本相符。

2)跨臨界射流中，在偽沸騰發(fā)生之前，低溫液核周圍的氣液混合層內(nèi)存在強密度分層，從而明顯抑制了射流的徑向擴散，使低溫液氮在軸線方向擴散更遠。超臨界射流中，射流進入混合狀態(tài)更早，并且由于混合層中的密度梯度較小，超臨界射流的混合層更寬。

3)跨臨界射流吸收的能量主要用于自身的體積膨脹，而非提高自身溫度，混合層中的密度變化也與其有關；并且在高環(huán)境壓力下，跨臨界射流的偽沸騰現(xiàn)象明顯減弱。對于超臨界射流，偽沸騰現(xiàn)象的缺失使射流從高溫流體吸收的熱量完全用于自身的升溫。