范珍涔，范 瑋，靳 樂(lè)

（1.中航工業(yè)燃?xì)鉁u輪研究院，成都610500；2.西北工業(yè)大學(xué)動(dòng)力與能源學(xué)院，西安710072）

0 引言

隨著高超聲速飛行器的發(fā)展，對(duì)高性能動(dòng)力裝置的需求越來(lái)越強(qiáng)烈。隨著飛行速度的提高，發(fā)動(dòng)機(jī)的工作參數(shù)將會(huì)發(fā)生一定改變，其燃燒室中燃燒的初始環(huán)境溫度和壓力有望超過(guò)航空煤油的臨界值。另外，高超聲速飛行器對(duì)冷卻系統(tǒng)的要求不斷提高，作為飛行器的惟一自帶介質(zhì)的燃油將作為首選的冷卻介質(zhì)，燃油在冷卻機(jī)體和發(fā)動(dòng)機(jī)熱端部件的同時(shí)自身溫度也被提升，逐漸達(dá)到并超過(guò)燃油的臨界點(diǎn)溫度。因此，未來(lái)高速飛行器動(dòng)力裝置的燃燒室中必定會(huì)出現(xiàn)超臨界工況，而超臨界噴射、蒸發(fā)和燃燒過(guò)程與常規(guī)亞臨界工況差別極大。探索液態(tài)碳?xì)淙剂系某R界噴射和燃燒機(jī)理成為高性能空天動(dòng)力領(lǐng)域亟需解決的關(guān)鍵問(wèn)題之一。

對(duì)超臨界流體的基礎(chǔ)研究表明[1]，當(dāng)?shù)湫鸵簯B(tài)碳?xì)淙剂稀娇彰河吞幱诔R界狀態(tài)時(shí)，其既不屬于液相也不屬于氣相，煤油的物性表現(xiàn)為既具有液態(tài)煤油的密度又具有氣體的性質(zhì)，具有等價(jià)的氣相和液相密度值，兼具氣體的低黏度和液體的高密度以及介于氣體和液體之間的高擴(kuò)散系數(shù)等特征，并且其物性參數(shù)隨著溫度和壓力的改變而顯得非常敏感。這些物性的特殊性是造成超臨界噴射及其混合燃燒過(guò)程與亞臨界情況下迥異的直接原因，國(guó)內(nèi)外學(xué)者針對(duì)超臨界碳?xì)淙剂系膰娚溥M(jìn)行了相關(guān)研究并得到了一定成果。Newman[2]在1971年就考慮到未來(lái)先進(jìn)發(fā)動(dòng)機(jī)的發(fā)展趨勢(shì)，率先對(duì)近臨界熱力學(xué)條件下的液態(tài)噴射進(jìn)行了研究，發(fā)現(xiàn)了表面張力大大減小的現(xiàn)象，拉開(kāi)了超臨界噴射研究的帷幕；Mayer等[3]試驗(yàn)研究了液氮和液氮/氦氣混合物噴射到氮?dú)夂蜌鈶B(tài)氮/氦混合物超臨界環(huán)境中的噴射過(guò)程；Wu等[4-5]對(duì)超臨界乙烯噴射到靜止的亞臨界氮?dú)猸h(huán)境中的現(xiàn)象，試驗(yàn)研究了臨界點(diǎn)附近輸運(yùn)參數(shù)變化對(duì)激波結(jié)構(gòu)的影響；Oscbwald[6]利用拉曼光譜法對(duì)超臨界氮?dú)庾杂蓢娚涞脑囼?yàn)研究發(fā)現(xiàn)，在臨界點(diǎn)附近的噴射所吸收的熱量不會(huì)提高流體的溫度而使得噴射更為擴(kuò)展，同時(shí)表面張力消失，噴射速度和動(dòng)量對(duì)于噴射的影響很小，噴射的主要控制因素是環(huán)境的熱力學(xué)參數(shù)；Chehroudi等[7]分別對(duì)亞臨界流體噴射到超臨界環(huán)境中噴射與可變密度的湍流氣體噴射的噴射擴(kuò)張角進(jìn)行了對(duì)比試驗(yàn)研究，分析發(fā)現(xiàn)2種噴射在噴射擴(kuò)張角上存在定量的相似；Barata等[8]利用可變密度噴射的數(shù)值模型對(duì)試驗(yàn)數(shù)據(jù)進(jìn)行了驗(yàn)證計(jì)算，得到了較好的吻合；Sui等[9]在非常接近臨界點(diǎn)的情況下進(jìn)行了超臨界SF6流體噴射到靜止的N2和CO2高壓環(huán)境的試驗(yàn)研究；Chen[10]研究發(fā)現(xiàn)超臨界噴射的射流穿透長(zhǎng)度的主要影響因素是噴射溫度和腔室的壓力；Douthip[11]通過(guò)試驗(yàn)和數(shù)值計(jì)算對(duì)比研究了超臨界噴射長(zhǎng)度和噴射擴(kuò)張角隨著噴射流量和溫度的變化規(guī)律，并與相同流量和出口壓力下的亞臨界噴射進(jìn)行了對(duì)比分析，發(fā)現(xiàn)超臨界的噴射長(zhǎng)度要小于亞臨界的。上述研究中，絕大部分都以小分子、單組分的燃料為研究對(duì)象，使用實(shí)際液態(tài)碳?xì)淙剂系臉O少。同時(shí)，大部分研究針對(duì)跨臨界噴射，對(duì)于超臨界流體噴射到超臨界環(huán)境中的研究很少。

本文以液態(tài)碳?xì)淙剂系某R界噴射為背景，考慮超臨界流體物性的特殊性，對(duì)比PR方法與理想氣體方法應(yīng)用于超臨界噴射的可行性，得出了能夠較為準(zhǔn)確模擬超臨界噴射的數(shù)值模型。

1 PR狀態(tài)方程及其熱力學(xué)關(guān)系

在超臨界噴射中應(yīng)該考慮到實(shí)際氣體效應(yīng)，壓縮因子z 的定義為

式中：T 和P 分別為溫度和壓力；Vm為摩爾體積；R 為氣體常數(shù)。

對(duì)于理想氣體來(lái)說(shuō)z=1，而對(duì)于超臨界流體來(lái)說(shuō)，通常z<1。

Peng-Robinson[12]考慮到實(shí)際氣體的效應(yīng)，建立了2參數(shù)的狀態(tài)方程，即

式中：參數(shù)a 和b 的值可以由物質(zhì)的臨界性質(zhì)求得，其與臨界壓縮因子zc的具體取值見(jiàn)表1。

表1 PR狀態(tài)方程中參數(shù)的取值

從表中可見(jiàn)，fω是1個(gè)關(guān)于偏心因子ω 的函數(shù)，可以由式（3）求出。Pc、Tc分別為物質(zhì)的臨界壓力、臨界溫度。典型物質(zhì)的臨界狀態(tài)參數(shù)見(jiàn)表2。從表中可見(jiàn)，一些典型物質(zhì)的臨界參數(shù)、壓縮因子和偏心因子的值，大部分物質(zhì)的實(shí)際臨界壓縮因子值與上述方法計(jì)算得到的值都非常接近。

表2 典型物質(zhì)的臨界狀態(tài)參數(shù)

進(jìn)行超臨界流體熱力學(xué)參數(shù)的計(jì)算可以采用相對(duì)偏差法：選定參考狀態(tài)，將熱力學(xué)參數(shù)通過(guò)狀態(tài)方程與P、V、T 相關(guān)聯(lián)。各種熱力學(xué)參數(shù)的解析表達(dá)形式為

式中：S 和H 分別為熵和焓；S0和H0分別為參考狀態(tài)下的熵和焓值。

對(duì)于定壓比熱和定容比熱，同樣可以采用相對(duì)偏差法來(lái)求得。真實(shí)氣體比熱可以表示為

式中：ΔCp和ΔCv稱(chēng)為剩余比熱容，表達(dá)式為

式（8）、（9）定義了真實(shí)氣體比熱和理想氣體比熱之間的關(guān)系，通過(guò)將PR狀態(tài)方程代入上述公式中，即可以求出所有的熱力學(xué)參數(shù)值。

2 數(shù)值模型的建立

2.1 物理模型

本文主要對(duì)航空煤油的常用替代燃料——C10H22（正癸烷）噴射到超臨界N2環(huán)境的超臨界噴射過(guò)程進(jìn)行了數(shù)值模擬。計(jì)算的物理模型如圖1所示。從圖中可見(jiàn)，燃料進(jìn)口直徑為0.254mm，N2進(jìn)口寬度為3.968mm，噴射腔室充滿N2。具體邊界和工況條件的設(shè)置與Douthip[11]的試驗(yàn)中的保持一致，見(jiàn)表3。

圖1 物理模型

表3 數(shù)值計(jì)算模型的邊界和工況條件參數(shù)

2.2 數(shù)值模型

將控制方程微分化之后，可以得到

式中：q 向量為相關(guān)的獨(dú)立守恒變量；f、g 和h 分別為在3個(gè)空間方向上的通量；S 為源項(xiàng)。下標(biāo)i 和v 分別為無(wú)黏和黏性流動(dòng)項(xiàng)?；诶字Z平均的Navier-Stokes方程，獨(dú)立變量和無(wú)黏通量可以表示為

式中：e 為單位體積內(nèi)的總能量；u、v 和w 分別為各方向的速度；σi為在湍流模型中的一些湍流傳輸特性，例如湍流動(dòng)力學(xué)能等。

上述結(jié)構(gòu)形式適用于多相流計(jì)算。式中前5 行表示標(biāo)準(zhǔn)的歐拉方程，分別為能量方程、連續(xù)方程和3個(gè)動(dòng)量方程。將實(shí)際氣體的狀態(tài)方程（P=zρRT）用守恒變量表示為

式中：γ 為比熱比：z 可通過(guò)狀態(tài)方程得出。黏性項(xiàng)定義為

式中：K 為熱傳導(dǎo)系數(shù)；D 為擴(kuò)散系數(shù)；τij為黏性壓力，定義為

式（14）首先假定了斯托克斯法在氣體計(jì)算中的正確性，因此衍生出第2黏性系數(shù)κ，其值為動(dòng)力黏度μ 的2/3。

溫度同樣可以通過(guò)狀態(tài)方程由守恒變量求出，即

源項(xiàng)則可表示為

式中：gx、gy和gz為體作用力；Ωi為源項(xiàng)。

3 結(jié)果分析

利用上文的物理和數(shù)值模型對(duì)采用了PR狀態(tài)方程的實(shí)際氣體模型和理想氣體模型進(jìn)行了對(duì)比數(shù)值研究。2種模型對(duì)進(jìn)口參數(shù)預(yù)測(cè)的相對(duì)誤差（理想氣體模型預(yù)測(cè)值－實(shí)際氣體模型預(yù)測(cè)值）/ 理想氣體模型預(yù)測(cè)值×100%）見(jiàn)表4。從表中可見(jiàn)，2種模型對(duì)氮?dú)鈪?shù)的預(yù)測(cè)差距很小，而對(duì)燃油參數(shù)的預(yù)測(cè)差距較大。

表4 2種計(jì)算模型對(duì)進(jìn)口參數(shù)預(yù)測(cè)的相對(duì)差異

表中Re 的數(shù)值可根據(jù)式（17）由進(jìn)口流速V、進(jìn)口處動(dòng)力黏度μ、密度ρ 和進(jìn)口直徑Dinj求出。

為了評(píng)估這種影響，基于數(shù)值計(jì)算結(jié)果，對(duì)2種模型所求得的噴射穩(wěn)定后的各物性參數(shù)的分布進(jìn)行了對(duì)比。

3.1 密度的分布

2種模型求得的噴射穩(wěn)定后的密度分布和不同徑向距離上密度沿軸向的分布情況分別如圖2、3所示。圖中，橫坐標(biāo)X/Dinj為軸向距離X 與噴嘴直徑Dinj之間的比值，縱坐標(biāo)Y/Dinj為徑向距離Y 與噴嘴直徑Dinj之間的比值。通過(guò)計(jì)算結(jié)果可知，圖2中的實(shí)際氣體模型模擬結(jié)果與文獻(xiàn)[11]中的數(shù)值模擬和試驗(yàn)結(jié)果較為吻合，而理想氣體模型模擬結(jié)果的誤差很大。理想氣體的密度預(yù)測(cè)偏低，速度預(yù)測(cè)值偏高，動(dòng)量預(yù)測(cè)偏大。因此，造成密度較大的區(qū)域分布長(zhǎng)度預(yù)測(cè)偏高。

圖2 2種模型求得的密度分布

圖3 不同徑向位置的密度沿軸向的分布

從圖3中可見(jiàn)，距離噴嘴中軸線越遠(yuǎn)的位置，2種模型的密度預(yù)測(cè)值越接近。越遠(yuǎn)離中軸線，氮?dú)饬吭酱?，理想氣體模型預(yù)測(cè)氣態(tài)密度的準(zhǔn)確性較高，因而越遠(yuǎn)離中軸線，理想氣體模型密度預(yù)測(cè)越接近真實(shí)值。同時(shí)，在理想氣體模型預(yù)測(cè)結(jié)果中，沿軸向方向，在X/Dinj=45～50處出現(xiàn)陡降，而在實(shí)際氣體模型預(yù)測(cè)結(jié)果中，該現(xiàn)象發(fā)生于X/Dinj=60～70處。這是由于理想氣體模型預(yù)測(cè)超臨界燃油密度值偏低，與氮?dú)獾拿芏炔罹喔≡斐傻?。燃油和氮?dú)饷芏炔罹嘣叫。诨旌线^(guò)程中，在相同氮?dú)赓|(zhì)量分?jǐn)?shù)下，混合體系密度越接近氮?dú)?，出現(xiàn)明顯下降的位置越靠前。

3.2 溫度

2種模型求得的噴射穩(wěn)定后的溫度分布如圖4所示。從圖中可見(jiàn)，由實(shí)際氣體模型求出的溫度分布場(chǎng)在軸向溫度逐漸升高，在X/Dinj≈65時(shí)與環(huán)境氣體的溫度一致，而由理想氣體模型求出的溫度分布場(chǎng)則在軸向在計(jì)算域內(nèi)一直未能與環(huán)境氣體的溫度一致。

圖4 2種模型求得的溫度分布

由于理想氣體的定壓比熱在溫度較高的情況下，預(yù)測(cè)值比實(shí)際氣體的大得多，偏離真實(shí)情況量太大，因此升溫過(guò)程的模擬變得很緩慢，造成了圖中這種明顯的區(qū)別。從圖中還可見(jiàn)，在噴射溫度較低的核心區(qū)域長(zhǎng)度內(nèi)，理想氣體模擬得到的結(jié)果比實(shí)際氣體模擬得到的結(jié)果要低一些，這同樣是由比熱預(yù)測(cè)的差異造成的，在溫度較低時(shí)，采用理想氣體方法求出的比熱值比實(shí)際情況的低，造成升溫過(guò)程較快，低溫核心區(qū)域長(zhǎng)度預(yù)測(cè)偏小。2種模型求得的燃油和環(huán)境的交界面情況有所不同，由理想氣體模型求得的溫度分布云圖出現(xiàn)了湍流特征，這是由于理想氣體模型求得的噴射雷諾數(shù)過(guò)大所造成的。

不同徑向位置上的溫度沿軸向的分布如圖5所示。從圖中可見(jiàn)，在燃油溫度較低的情況下（X/Dinj<40），理想氣體模型求得的溫度隨軸向距離的增大變化相對(duì)劇烈；而在燃油溫度較高的情況下（X/Dinj>40），理想氣體求得的溫度隨軸向距離增大的變化相對(duì)緩慢。這同樣是由于理想氣體熱力學(xué)參數(shù)預(yù)測(cè)的巨大誤差造成的，在燃油溫度較低的情況下，理想氣體模型預(yù)測(cè)的定壓比熱比實(shí)際值低，而在燃油溫度較高的情況下，理想氣體模型預(yù)測(cè)的定壓比熱比實(shí)際值高。

圖5 不同徑向位置上的溫度沿軸向的分布

3.3 質(zhì)量分?jǐn)?shù)分布的預(yù)測(cè)

圖6 不同軸向位置上的燃油質(zhì)量分?jǐn)?shù)沿軸向的分布

圖7 不同徑向位置上的燃油質(zhì)量分?jǐn)?shù)沿徑向的分布

不同軸向和徑向位置上燃油質(zhì)量分?jǐn)?shù)的分布分別如圖6、7所示。從圖6中可見(jiàn)，理想氣體模型計(jì)算得出的燃油質(zhì)量分?jǐn)?shù)減小為零的軸向位置大概在X/Dinj=50處，而實(shí)際氣體模型得出的則約在X/Dinj=70處。這說(shuō)明理想氣體模型得出的燃油軸向擴(kuò)散區(qū)比實(shí)際氣體模型模擬得出的低很多，這是由于理想氣體模型只適用于中低壓條件下的流體噴射擴(kuò)散過(guò)程的模擬，且由于理想氣體模型在預(yù)測(cè)密度、溫度等物理量時(shí)均有一定誤差，而這些物理量都是影響擴(kuò)散過(guò)程的重要因素。同樣，圖7中的沿徑向燃油質(zhì)量分?jǐn)?shù)的分布預(yù)測(cè)的差異也是上述原因造成的。

3.4 噴射長(zhǎng)度和噴射擴(kuò)張角的預(yù)測(cè)

噴射長(zhǎng)度和噴射擴(kuò)張角的計(jì)算方法與文獻(xiàn)[11]中一致，即：噴射長(zhǎng)度定義為噴射穩(wěn)定后噴嘴出口到中心軸線上C10H22質(zhì)量分?jǐn)?shù)為0.2位置的距離；噴射擴(kuò)張角定義為噴嘴出口處密度明顯分界面的切線夾角。2 種不同模型計(jì)算得出的噴射長(zhǎng)度和噴射擴(kuò)張角的結(jié)果分別如圖8、9所示。計(jì)算結(jié)果與Doungthip試驗(yàn)值[11]的對(duì)比見(jiàn)表5。

圖9 2種模型對(duì)噴射擴(kuò)張角的計(jì)算結(jié)果

表5 噴射長(zhǎng)度和擴(kuò)張角計(jì)算值與試驗(yàn)值的對(duì)比

從表中可見(jiàn)，采用實(shí)際氣體模型能夠較為精確地預(yù)測(cè)出噴射長(zhǎng)度和擴(kuò)張角的值，而采用理想氣體方法得到的結(jié)果與實(shí)際情況偏差非常大。這是由于噴射長(zhǎng)度和擴(kuò)張角的定義是基于密度或者質(zhì)量分?jǐn)?shù)的，而理想氣體模型預(yù)測(cè)密度和質(zhì)量分?jǐn)?shù)都有較大誤差，從而造成預(yù)測(cè)噴射特性產(chǎn)生巨大誤差。

4 結(jié)論

以未來(lái)高性能航空發(fā)動(dòng)機(jī)燃燒室中可能出現(xiàn)的超臨界工況為研究背景，基于PR狀態(tài)方程法，建立了考慮超臨界噴射物性特點(diǎn)的數(shù)值計(jì)算模型。并將該計(jì)算模型與基于理想氣體的計(jì)算模型進(jìn)行了對(duì)比研究，分別得到了不同噴射模型的密度、溫度、質(zhì)量分?jǐn)?shù)分布以及超臨界噴射長(zhǎng)度和擴(kuò)張角的變化規(guī)律和差異性。通過(guò)與已有的試驗(yàn)結(jié)果對(duì)比，驗(yàn)證了基于PR方法模型的有效性和可靠性。該數(shù)值模型的建立和驗(yàn)證將為研究碳?xì)淙剂系某R界噴射現(xiàn)象提供1種有效的數(shù)值工具。

[1]Holland M，Eaton E，Hanley H.A correlation of the viscosity and thermal conductivity data of gaseous and liquid ethylene[J].Journal of Physical and Chemical Reference Data，1983，12（4）：917-932.

[2]Newman J，Brzustowski T.Behavior of a liquid jet near the thermodynamic critical region[J].AIAA Journal，1971，9（8）：1595-1602.

[3]Mayer W，Schik A，Vielle B，et al.Atomization and breakup of cryogenic propellants under high-pressure subcritical and supercritical conditions [J].Journal of Propulsion and Power，1998，14（5）：835-842.

[4]Wu P，Chen T，Abdollah S，et al.Injection of supercritical ethylene in Nitrogen[J].Journal of Propulsion and Power，2006，22（4）：809-819.

[5]Wu P，Shahnam M，Kirkendall K，et al.Expansion and mixing processes of under expanded supercritical fuel jets injection into superheated conditions[R].AIAA-97-2852.

[6]Oscbwald M，Schik A.Supercritical Nitrogen free jet investigation by spontaneous Raman scattering [J].Experiments in Fluids，1999，27（6）：497-506.

[7]Chehroudi B，Talley D，Coy E.Visual characteristics and initial growth rates of round cryogenic jets at subcritical and supercritical pressures[J].Physical Fluids，2002，14（2）：850-861.

[8]Barata J，Gokalp I，Silva A.Numerical study of cryogenic jets under supercritical conditions[J].Journal of Propulsion and Power，2003，19（1）：142-147.

[9]Chen L，Sui P.Atomization during the injection of supercritical fluid into high-pressure environment[C]//International Union of Theoretical and Applied Mechanics(IUTAM)Symposia on Droplets and Sprays，Taiwan：Springer，1994：6-10.

[10]Chen L.Heat transfer，fouling and combustion of supercritical fuels[R].AFOSR-TR-940321，1994.

[11]Douthip T，Ervin J，Williams T，et al.Studies of injection of jet fuel at supercritical conditions [J].Industrial and Engineering Chemistry Research，2002，41（23）：5856-5866.

[12]Peng D Y，Robinson D B.A new two-constant equation of state[J].Industrial and Engineering Chemistry Research Fundamentals，1976，15（1）：59-64.