超臨界噴霧中渦的演化過程數(shù)值研究

李 瀟，王佔元，劉宗寬，周 磊，衛(wèi)海橋

（天津大學(xué)內(nèi)燃機燃燒學(xué)國家重點實驗室，天津 300072）

0 概述

隨著環(huán)境污染問題的日益嚴(yán)重和化石能源短缺問題的不斷加劇，內(nèi)燃機領(lǐng)域的排放法規(guī)也逐漸嚴(yán)格。優(yōu)化燃燒過程和提高內(nèi)燃機的熱效率是降低污染物排放的有效手段，而燃油的霧化及混合氣的形成直接影響燃燒的進行，對燃油在缸內(nèi)的擴散混合過程進行研究有重要意義。

在現(xiàn)代發(fā)動機上，先進燃燒模式和增壓技術(shù)的使用提高了熱效率，也使得缸內(nèi)溫度和壓力不斷增高。在某些高增壓發(fā)動機上，燃油在噴射初始時刻缸內(nèi)壓力和溫度可達2.5 MPa 和700 K，已經(jīng)超過大部分烴類燃油的臨界壓力和溫度，例如正十二烷的臨界壓力和溫度分別為1.82 MPa 和658 K。由此可見，缸內(nèi)壓力和溫度超過燃料臨界點的狀況會大概率發(fā)生，對超臨界環(huán)境下的噴霧特性進行研究顯得尤為重要。

超臨界流體具有獨特性質(zhì)，文獻［1］中通過試驗證明超臨界流體密度類似液體，黏度類似氣體，而擴散性質(zhì)介于液體和氣體之間。由于超臨界流體的特殊物性，燃油在超臨界環(huán)境下的噴霧也表現(xiàn)出不同于亞臨界環(huán)境的特有屬性。亞臨界環(huán)境下由于表面張力的存在，射流與周圍氣體之間存在明確的氣液界面，而在超臨界環(huán)境中的噴霧過程中，表面張力趨近于零導(dǎo)致氣液界面消失，擴散過程靠湍流擴散作用控制，射流表面各參數(shù)存在較大的梯度變化。關(guān)于對超臨界流體熱物性的預(yù)測，文獻［2］中指出超臨界流體分子間的相互作用不能被忽略，需使用真實流體狀態(tài)方程。超臨界流體的輸運系數(shù)（如黏度系數(shù)和導(dǎo)熱系數(shù)）的計算也要考慮流體超臨界狀態(tài)的影響。

超臨界噴霧的試驗研究已經(jīng)取得了較大進展，文獻［3］中通過試驗發(fā)現(xiàn)在超臨界環(huán)境下射流表面不會產(chǎn)生亞臨界環(huán)境下出現(xiàn)的液滴和液絲，而由液態(tài)條狀組成的梳子狀結(jié)構(gòu)代替。文獻［4］中利用平面激光誘導(dǎo)熒光（planar laser induced fluoresence，PLIF）技術(shù)，實現(xiàn)了在熱力學(xué)大范圍內(nèi)對流體密度和密度梯度的測量。許多學(xué)者也開展了關(guān)于超臨界噴霧的數(shù)值模擬研究。文獻［5］中基于大渦模擬，利用描述真實流體行為的計算模型，對正庚烷噴入超臨界環(huán)境下的行為進行了分析，發(fā)現(xiàn)密度梯度在噴霧過程中起著重要作用。文獻［6］中研究了不同超臨界壓力對低溫液氮噴霧的影響，指出超臨界下產(chǎn)生的密度梯度層會抑制射流徑向速度的振蕩，而促進軸向速度的振蕩，且環(huán)境壓力的增加將導(dǎo)致射流提前過渡到自相似狀態(tài)。文獻［7］中利用直接數(shù)值仿真方法（direct numerical simulation，DNS）對超臨界射流中混合層的發(fā)展進行了分析，證明密度梯度對混合層的穩(wěn)定性和湍流特性起關(guān)鍵作用。

從對超臨界噴霧的研究來看，目前的試驗研究證明了在超臨界環(huán)境下噴霧形態(tài)的特殊性。由于超臨界流體的特殊性，亞臨界環(huán)境下的數(shù)值模擬策略不再適用，超臨界環(huán)境下的噴霧模型構(gòu)建具有挑戰(zhàn)性。關(guān)于密度梯度在超臨界噴霧中的作用原理及超臨界環(huán)境下射流擴散演化過程的影響因素研究更是少見。

本研究基于上述超臨界噴霧的研究現(xiàn)狀，以開源軟件OpenFoam 為計算平臺，構(gòu)建了適合超臨界環(huán)境的噴霧模型，分析了超臨界條件下渦的演化過程。模型采用真實流體的狀態(tài)方程及熱物性和輸運系數(shù)的計算方法，結(jié)合大渦模擬通過對比超臨界和亞臨界環(huán)境下噴霧過程，重點研究了超臨界和亞臨界環(huán)境下燃油擴散混合過程的區(qū)別與聯(lián)系，從根本上揭示了超臨界環(huán)境下射流擴散的特殊性質(zhì)。

1 理論模型

1.1 大渦模擬基本方程

本文中采用大渦模擬的方法對超臨界噴霧過程進行了模擬，大渦模擬是將流場在空間上進行過濾，通過濾波函數(shù)使大渦和小渦分離，大尺度渦可以直接通過數(shù)值模擬求解，而小尺度渦（亞網(wǎng)格尺度）需要亞網(wǎng)格模型模化求解，因此大渦模擬能捕捉到雷諾平均無法捕捉到的流場細(xì)微結(jié)構(gòu)，得到的計算結(jié)果更接近真實噴霧形態(tài)。

控制方程為可壓縮形式的Navier-Stokes（NS）方程，包括連續(xù)方程、動量方程、能量方程和組分方程，見式（1）～式（4）。黏性應(yīng)力張量表示見式（5）。

式中，ρ為密度；v為速度；p為壓力；hs為比焓；Yi為組分質(zhì)量分?jǐn)?shù)；τ為黏性應(yīng)力張量；Sρ、Sv、Sh和SY分別為密度、速度、焓和組分源項；μe和αe分別為有效動力黏度和有效熱擴散系數(shù)；(?v)T為速度梯度的轉(zhuǎn)置；I為單位矩陣。

液相與氣相間的傳質(zhì)及其他作用通過在控制方程中添加作用項實現(xiàn)，因噴霧過程中不考慮燃燒反應(yīng)，所以控制方程中未求解由于化學(xué)反應(yīng)產(chǎn)生的源項。

在大渦模擬中，上述控制方程經(jīng)過空間濾波將會產(chǎn)生亞網(wǎng)格尺度（subgrid scale，SGS）項，例如亞網(wǎng)格應(yīng)力張量τSGS、亞網(wǎng)格熱流量hSGS等無法直接求解的量，因此需要亞網(wǎng)格模型來封閉控制方程。本文選取的亞網(wǎng)格模型為動態(tài)單方程模型［8］，其結(jié)構(gòu)簡單且有足夠的精度，得到廣泛的應(yīng)用。

1.2 狀態(tài)方程及熱動力學(xué)屬性

超臨界流體的部分參數(shù)會偏離亞臨界環(huán)境下的變化規(guī)律，為了準(zhǔn)確描述流體的物性參數(shù)變化，需要采用真實流體的狀態(tài)方程及熱物性和輸運系數(shù)的計算方法。由于立方型狀態(tài)方程具有形式簡單、精度高、計算成本低等優(yōu)點被大量使用，常用立方型狀態(tài)方程有Soave-Redlich-Kwong（SRK）和Peng-Robinson（PR）。文獻［9］中對比了SRK 和PR 兩種狀態(tài)方程對超臨界環(huán)境下噴霧預(yù)測的差異性，根據(jù)其研究結(jié)果，本文中選用的是PR 狀態(tài)方程，具體形式見式（6）～式（10）。

式中，T為流體的溫度；Mw為流體的物質(zhì)的量；pc和Tc分別為臨界壓力和臨界溫度；ω為流體的偏心因子；Ru為氣體常數(shù)；Spr、a和b為物質(zhì)性質(zhì)決定的常數(shù)。

超臨界流體黏度系數(shù)和導(dǎo)熱系數(shù)的計算選用文獻［10］中提出的計算方法，基本思想為利用特定的混合準(zhǔn)則計算多組分混合物在低壓下的黏度系數(shù)和導(dǎo)熱系數(shù)，之后利用修正函數(shù)得到高溫高壓狀態(tài)下的黏度系數(shù)和導(dǎo)熱系數(shù)。以黏度系數(shù)為例給出了計算過程。

式中，c、d為模型常數(shù)；η0為低壓黏度系數(shù)；η為高壓黏度系數(shù)；ηp為模型中間變量，無實際物理意義；Mwm、Vcm和Tcm分別為混合物的物質(zhì)的量、臨界體積和臨界溫度；k和?分別為勢能和玻爾茲曼常數(shù)；A6～A10為偏心因子和偶極矩的線性函數(shù)；Ω*、Fcm、T*和G2的具體計算方法見文獻［10］。

質(zhì)量擴散系數(shù)采用文獻［11］中提出的高溫高壓修正方法，其利用Tr構(gòu)造了超臨界狀態(tài)下的修正函數(shù)，具體形式見式（16）。

式中，DH和DL分別為高壓和低壓下的質(zhì)量擴散系數(shù)；Tr為溫度與臨界溫度的比值；A、B、C和E是隨壓力與臨界壓力的比值pr變化的模型經(jīng)驗系數(shù)。

上述超臨界環(huán)境下物性參數(shù)計算模型的準(zhǔn)確性已從模型自身［10-11］和其在超臨界環(huán)境下的噴霧模擬中得到試驗數(shù)據(jù)的驗證，大量關(guān)于超臨界狀態(tài)下的模擬研究［1，12］均使用了此計算模型。

2 數(shù)值模擬

針對超臨界噴霧現(xiàn)象，本文中選取的模擬對象為美國Sandia 國家實驗室關(guān)于正十二烷的噴射試驗數(shù)據(jù)，該試驗被大量用來進行超臨界噴霧的研究。文獻［13］中的研究證明了此試驗中存在超臨界狀態(tài)的油氣混合過程。根據(jù)已有試驗數(shù)據(jù)建立噴霧模型，該試驗將正十二烷通過直徑為0.09 mm 的噴孔噴入定容彈內(nèi)，由于不考慮燃燒過程，定容彈內(nèi)的環(huán)境氣體中不含O2，由N2、CO2和H2O 組成，密度約為22.8 kg/m3。初始溫度為363 K 的正十二烷在150 MPa 的壓力下噴入溫度900 K、壓力6.0 MPa 的環(huán)境中，而正十二烷的臨界溫度和臨界壓力分別為658 K 和1.82 MPa，在噴霧過程中會跨越燃油臨界點，使得燃油在超臨界環(huán)境中進行霧化和混合。表1 給出了模擬中具體的邊界條件，其中N2、CO2和H2O 的物質(zhì)的量分?jǐn)?shù)分別為89.71 %、6.52 %和3.77 %。

表1 邊界條件

為探究超臨界和亞臨界環(huán)境下燃油擴散混合的異同點，在上述模擬條件的基礎(chǔ)上，更改環(huán)境氣體的溫度和壓力，使整個噴霧過程都處于亞臨界環(huán)境中。本文中根據(jù)文獻［14］中的研究，環(huán)境溫度和壓力設(shè)為440 K 和2.9 MPa，并保證了環(huán)境氣體密度不變，均為22.8 kg/m3，表2 給出了兩種環(huán)境下的算例設(shè)置。

表2 不同環(huán)境下算例設(shè)置

本研究所采用的計算流體力學(xué)（computational fluid dynamics，CFD）軟件平臺為OpenFoam，該軟件平臺具有開源、免費等優(yōu)點，其自帶求解器類型齊全，泛用性強并可基于用戶需求自行擴展，得到廣泛的應(yīng)用。本文中使用的SprayFoam 求解器基于歐拉-拉格朗日噴霧模型，求解器耦合PR 狀態(tài)方程和真實流體輸運參數(shù)的計算方法，并結(jié)合大渦模擬策略，構(gòu)建了適用于超臨界環(huán)境下的噴霧模型。

計算域為直徑16 mm、高50 mm 的圓柱型區(qū)域，采用O-Block 對圓柱型計算域進行劃分，并對局部網(wǎng)格進行加密，總網(wǎng)格數(shù)量約為330 萬，最小網(wǎng)格尺寸為50 μm。文獻［15］中研究了250.00 μm、125.00 μm、62.50 μm 和41.67 μm 網(wǎng)格對正十二烷噴射試驗?zāi)M結(jié)果的影響，其指出大渦模擬采用62.50 μm 網(wǎng)格便能準(zhǔn)確捕捉噴霧特征，由此可知本文中使用的網(wǎng)格可以保證計算結(jié)果的準(zhǔn)確性。時間步長由最大庫朗數(shù)控制，最大庫朗數(shù)為0.33，時間步長約為3×10-8s。

3 結(jié)果分析

圖1 為超臨界和亞臨界環(huán)境下試驗和計算的貫穿距隨噴射時間變化的對比圖。其中氣相貫穿距定義為燃油質(zhì)量分?jǐn)?shù)為0.1% 時對應(yīng)的位置與燃油噴嘴之間的最遠(yuǎn)距離，液相貫穿距定義為軸線方向上燃油質(zhì)量為已噴燃油90% 區(qū)域所占的軸向距離。從圖1 中可以看出計算模型能準(zhǔn)確捕捉不同環(huán)境下噴霧的液相和氣相貫穿距。在超臨界環(huán)境下的噴霧初期計算的液相貫穿距稍高于計算值，這也導(dǎo)致了氣相貫穿距的增加，亞臨界環(huán)境下也有類似的情況，但整體上計算結(jié)果和試驗值基本吻合。本文中構(gòu)建的噴霧求解器可以在寬范圍的溫度和壓力下準(zhǔn)確預(yù)測噴霧特征。

圖1 不同環(huán)境下試驗和計算噴霧貫穿距對比

為更好地證明超臨界環(huán)境下計算結(jié)果的準(zhǔn)確性，圖2 給出了超臨界環(huán)境下距離噴油位置處軸向距離X=17.85 mm 和25.00 mm 處沿徑向燃油質(zhì)量分?jǐn)?shù)分布的試驗和計算結(jié)果對比圖。從試驗和計算結(jié)果的對照中可以看出計算的燃油質(zhì)量分?jǐn)?shù)分布基本上與試驗結(jié)果吻合，在噴霧邊緣的燃油質(zhì)量分?jǐn)?shù)稍微高出了試驗值，這是由噴霧形態(tài)分布的不確定性所致。

圖2 不同軸向距離處沿徑向燃油質(zhì)量分?jǐn)?shù)分布

3.1 射流擴散率

超臨界環(huán)境下燃油液滴的蒸發(fā)、熱物性參數(shù)和油氣混合層的發(fā)展演化與亞臨界環(huán)境下相比都表現(xiàn)出明顯的特殊性，這也使超臨界環(huán)境下燃油的射流擴散過程不同于亞臨界。下面從燃油的分布入手研究兩種環(huán)境下的射流擴散的差異。

圖3 和圖4 分別給出了0.6 ms 時刻兩種環(huán)境下燃油蒸氣的密度分布圖，后面的分析都基于此時刻進行。從圖中可以明顯地看出亞臨界環(huán)境下噴霧有更長的貫穿距，而超臨界環(huán)境下噴霧沿徑向充分發(fā)展。從燃油蒸氣密度的大小可以看出，超臨界環(huán)境下的燃油充分蒸發(fā)，霧化效果更好。

圖3 超臨界環(huán)境下密度分布

為了更細(xì)致地描述圖3 和圖4 中的密度分布現(xiàn)象，本文中利用射流擴散率來量化燃油和環(huán)境氣體混合的程度。射流擴散率是一個能衡量射流擴散混合效率的指標(biāo)。射流擴散率可以利用射流的速度場和密度場計算，需要根據(jù)應(yīng)用場景進行選擇，本文中基于密度場的半寬半峰（half-width，half-maximum）規(guī)則［16］來計算射流擴散率和擴散角度。

圖4 亞臨界環(huán)境下密度分布

圖5 展示了基于射流密度場的半寬半峰規(guī)則計算射流擴散率的具體過程，其中Δρ為密度增加量，Lρ為半寬半峰長度。擴散率定義為Lρ沿射流軸線的線性回歸斜率。為了保證計算的準(zhǔn)確性，線性回歸必須在Lρ線性增加的區(qū)域進行，在射流中通常是自相擬區(qū)域。根據(jù)半寬半峰規(guī)則得到射流擴散角度，如式（17）所示。

圖5 基于密度場的半寬半峰規(guī)則

基于上述射流擴散率的概念，表3 給出了兩種環(huán)境下的射流擴散率和擴散角度。超臨界環(huán)境下的擴散率和擴散角度均大于亞臨界，其燃油擴散混合效果也優(yōu)于亞臨界，利用此方法得的超臨界環(huán)境下擴散角度與文獻［17］中的試驗結(jié)果一致性較好。

表3 擴散率和擴散角度

3.2 渦的分布

為了更加深入地研究兩種環(huán)境下擴散率差異的原因，需要從更微觀的層面進行分析，從根本上探究不同環(huán)境下射流擴散混合過程的影響因素。在高速射流噴霧過程中，射流與周圍環(huán)境流體間的速度梯度非常大，初始噴油區(qū)域附近會產(chǎn)生剪切層，剪切層受Kelvin-Helmholtz 不穩(wěn)定性的影響，會發(fā)生渦的卷起、配對和分裂等級聯(lián)過程。隨著渦的不斷產(chǎn)生并向下游運動，周圍氣體會不斷地被卷吸到射流中，從而使燃油與環(huán)境氣體混合。由此可見渦的分布對噴霧過程中油氣混合有著關(guān)鍵的影響，根據(jù)渦的分布可以分析兩種環(huán)境下油氣混合狀態(tài)的差異。

圖6 和圖7 給出了兩種環(huán)境下噴油位置附近區(qū)域內(nèi)渦量的分布圖。超臨界和亞臨界環(huán)境下渦的分裂位置分別約為距離噴油位置軸向距離X=6.50 mm 和10.00 mm 處?？梢钥闯龀R界環(huán)境下渦開始產(chǎn)生的位置更靠近射流上游，超臨界環(huán)境下渦更早產(chǎn)生，加速了射流的破碎，使得更多的周圍氣體被卷吸，這種作用造成噴霧初始階段超臨界環(huán)境下渦量分布表面多褶皺，而亞臨界環(huán)境下表面較光滑，超臨界環(huán)境下射流的破碎趨勢更加明顯。超臨界下渦更早地產(chǎn)生促進了油氣混合，同時也增強噴霧沿徑向的發(fā)展，擴展率增大，相應(yīng)造成貫穿距降低，此結(jié)論與圖3 和圖4 的密度分布相一致。

圖6 超臨界環(huán)境下渦量分布

圖7 亞臨界環(huán)境下渦量分布

超臨界環(huán)境下渦更早產(chǎn)生的原因有待探究，需要從更深層次的渦的演化角度分析。

3.3 渦的演化

渦的分布直接影響著射流的破碎及其后續(xù)混合過程的進行，為了進一步分析渦的產(chǎn)生與發(fā)展，需要對渦的演化及其影響因素進行探究，并量化各種影響因素所占的比例，從根本上分析超臨界和亞臨界環(huán)境下渦的分布不同的原因。根據(jù)文獻［18］中的研究，在不考慮其他外力的情況下，渦的演化方程為：

式中，d/dt代表物質(zhì)導(dǎo)數(shù)；右側(cè)第1 項表示渦的伸展，代表了速度梯度對渦的影響，它是產(chǎn)生復(fù)雜渦結(jié)構(gòu)的關(guān)鍵原因，并影響著湍流中動能級聯(lián)過程；右側(cè)第2 項為渦的膨脹，這是由于流體體積的收縮和膨脹造成的，代表著流體的可壓縮性對渦的影響；右側(cè)第3 項為壓力梯度與密度梯度錯位所產(chǎn)生的斜壓轉(zhuǎn)矩，是Rayleigh-Taylor 不穩(wěn)定性形成的原因；最后一項來自于由黏性應(yīng)力引起的渦的擴散。

圖8 和圖9 分別為噴油位置附近區(qū)域內(nèi)兩種環(huán)境下不同項對渦的演化影響分布圖，自上而下依次為渦的伸展項、膨脹項、斜壓轉(zhuǎn)矩項和黏性應(yīng)力項。為便于分析，圖中僅給出了不同項垂直于平面方向的大小，在其他方向上各項具有同樣的趨勢。根據(jù)圖中各項的大小可以看出在兩種環(huán)境下渦的伸展的作用近乎相同，渦的膨脹作用也有同樣的趨勢，并且渦的伸展對渦的演化影響最大（達到1012水平）。在亞臨界環(huán)境下的斜壓轉(zhuǎn)矩和黏性應(yīng)力的影響較小，而斜壓轉(zhuǎn)矩的影響稍高于黏性應(yīng)力。在超臨界環(huán)境下斜壓轉(zhuǎn)矩的影響增強，與渦的膨脹作用相當(dāng)（均為1011），此時黏性應(yīng)力造成的影響遠(yuǎn)低于斜壓轉(zhuǎn)矩。

圖8 超臨界環(huán)境中不同項的影響

圖9 亞臨界環(huán)境中不同項的影響

整體對照兩種環(huán)境下不同項在渦的演化過程中的作用，可以明顯地看出斜壓轉(zhuǎn)矩的作用在超臨界環(huán)境中會大幅增強，其他3 項的影響基本保持不變，且3 項中速度梯度的作用最強。在高速射流中（例如本研究中的射流速度最高可達600 m/s），射流表面和周圍氣體之間會形成巨大的速度梯度，成為渦的演化過程中的決定因素，而渦的膨脹和黏性應(yīng)力的影響在溫度和壓力劇烈變化后仍保持不變，說明周圍環(huán)境氣體的影響在高速射流中會相對降低，該現(xiàn)象是由于射流高速運動使得其與環(huán)境的相互作用時間縮短所致。

超臨界環(huán)境中斜壓轉(zhuǎn)矩的作用增強與其形成較大的密度梯度有關(guān)，圖10 和圖11 給出了兩種環(huán)境下的密度梯度分布。超臨界環(huán)境下的密度梯度明顯高于亞臨界，超臨界狀態(tài)下的大密度梯度使密度梯度和壓力梯度的錯位所產(chǎn)生的斜壓轉(zhuǎn)矩增大，強化了渦的發(fā)展，進而對混合層的演化過程造成重要的影響。由此可見，超臨界環(huán)境下產(chǎn)生的較大密度梯度是其提前出現(xiàn)渦的分離的關(guān)鍵影響因素，說明密度梯度對射流混合過程起著重要作用，密度梯度的大小直接影響著超臨界環(huán)境噴霧的發(fā)展。

圖10 超臨界環(huán)境下密度梯度分布

圖11 亞臨界環(huán)境下密度梯度分布

由此可見超臨界條件下渦更早地產(chǎn)生的原因是超臨界下形成的較大密度梯度層，密度梯度的增大加大了壓力梯度與密度梯度錯位所產(chǎn)生的斜壓轉(zhuǎn)矩，即渦的演化公式中?(1p) ×?p項增大，加速了渦的產(chǎn)生。超臨界環(huán)境下射流密度梯度的增大是由于射流進入超臨界環(huán)境時其與周圍環(huán)境的質(zhì)量擴散率接近于零［19］，射流與環(huán)境氣體的物質(zhì)交換受到抑制，使其形成一個高密度的混合層，混合層與環(huán)境氣體間產(chǎn)生較大的密度梯度。

綜上，欲強化渦的產(chǎn)生與演化以增強卷吸作用來強化燃料的噴霧混合：首先需要在噴油裝置允許的壓力條件下盡可能提高噴油壓力，提高射流與環(huán)境氣體間的速度梯度，加速渦的伸展；再者要提高環(huán)境溫度和壓力，使射流進入超臨界狀態(tài)，以提高射流與環(huán)境間壓力梯度和密度梯度錯位所產(chǎn)生的斜壓轉(zhuǎn)矩，增大Rayleigh-Taylor 不穩(wěn)定性。通過上述噴射策略，可大幅提高燃油霧化效率，使之與環(huán)境氣體充分混合，利于后續(xù)燃燒過程的進行。

4 結(jié)論

（1）基于射流密度分布，超臨界環(huán)境下的射流擴散率和擴散角度均大于亞臨界。

（2）超臨界環(huán)境會促進渦的形成，使周圍更多的氣體被卷吸入射流中，加速射流的破碎及后續(xù)的油氣混合過程。

（3）在高速射流中，速度梯度對渦的演化過程起主要作用，超臨界環(huán)境下由于密度梯度的增大會使斜壓轉(zhuǎn)矩的作用增強，加速渦的演化。