三維直接數(shù)值模擬部分預(yù)混發(fā)動(dòng)機(jī)內(nèi)的燃燒過(guò)程

張　帆　堯命發(fā)

(天津大學(xué)，內(nèi)燃機(jī)燃燒學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津300072)

?

三維直接數(shù)值模擬部分預(yù)混發(fā)動(dòng)機(jī)內(nèi)的燃燒過(guò)程

張帆堯命發(fā)*

(天津大學(xué)，內(nèi)燃機(jī)燃燒學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津300072)

采用三維直接數(shù)值模擬方法研究了一個(gè)類似于部分預(yù)混燃燒(PPC)發(fā)動(dòng)機(jī)條件下高辛烷值燃料PRF70的著火過(guò)程。文章采用了簡(jiǎn)化的PRF化學(xué)動(dòng)力學(xué)機(jī)理，包含33個(gè)組分和38步基元反應(yīng)。計(jì)算中根據(jù)發(fā)動(dòng)機(jī)的幾何尺寸和真實(shí)運(yùn)行工況加入了氣缸內(nèi)壓縮/膨脹的效果，并考慮了燃料的兩次噴射，其中第一次噴射形成了較均勻的混合氣，第二次燃料噴射增加了混合物分層。研究發(fā)現(xiàn)，PPC的燃燒過(guò)程非常復(fù)雜，是均質(zhì)壓燃、預(yù)混燃燒和擴(kuò)散燃燒三種主要燃燒模式的結(jié)合。在兩次燃料噴射之間的區(qū)域?yàn)榻瘜W(xué)計(jì)量比燃燒，是氮氧化物的生成區(qū)；而在化學(xué)計(jì)量比(φ)大于2的區(qū)域，混合不充分聚集了大量未燃碳?xì)浜虲O。文章使用Marching cube算法捕捉了三維火焰鋒面隨時(shí)間的變化。最后，使用反應(yīng)鋒面上高斯曲率(kg)與平均曲率(km)的聯(lián)合概率密度函數(shù)(PDF)以及平均曲率隨時(shí)間變化的概率密度函數(shù)，揭示了球形火焰鋒面和馬鞍形火焰鋒面的存在，前者占主要地位，并且隨著燃燒的進(jìn)行，負(fù)曲率增加，主要是因?yàn)橹行牡娜剂蠞鈪^(qū)在逐漸消耗。

直接數(shù)值模擬；部分預(yù)混燃燒；PRF70；燃燒特性；PPC發(fā)動(dòng)機(jī)

www.whxb.pku.edu.cn

1　引言

為了滿足越來(lái)越嚴(yán)格的排放法規(guī)，新一代的內(nèi)燃機(jī)必然向著高效、清潔的方向發(fā)展。傳統(tǒng)的汽油機(jī)效率較低而柴油機(jī)又由于碳煙排放過(guò)高使得研究人員不得不探索新的內(nèi)燃機(jī)燃燒模式。其中，部分預(yù)混燃燒(PPC)是新型的低溫燃燒模式的代表，它在柴油機(jī)中采用高辛烷值燃料，通過(guò)兩次或多次噴射，實(shí)現(xiàn)噴油與燃燒過(guò)程的分離，同時(shí)實(shí)現(xiàn)較高的熱效率、較低NOx、碳煙(soot)的排放，是一種介于傳統(tǒng)柴油機(jī)燃燒與均質(zhì)壓燃的一種燃燒模式1－7。

在PPC發(fā)動(dòng)機(jī)中，進(jìn)氣道首先噴入一定燃油形成較均勻的混合氣，在壓縮沖程缸內(nèi)直噴剩余燃油，不僅可以形成一定的濃度分層，也達(dá)到了引燃的目的。為了避免傳統(tǒng)柴油機(jī)邊噴油邊著火，PPC發(fā)動(dòng)機(jī)通常采用較高辛烷值的燃料，有利于油空混合和延長(zhǎng)滯燃期，并且由于多次噴射也避免了出現(xiàn)局部過(guò)濃區(qū)域，較好地避開soot生成2,5。

Manente等8－10使用不同辛烷值燃料(乙醇、汽油、汽柴油摻混等)研究了噴油時(shí)刻、噴油比例對(duì)PPC燃燒模式的影響，并指出最優(yōu)的混合氣濃度可以通過(guò)兩次噴油實(shí)現(xiàn)，而且高辛烷值燃料優(yōu)于低辛烷值燃料，因?yàn)榍罢呔哂懈L(zhǎng)的滯燃期、使得NOx很低但是未燃碳?xì)?UHC)和CO有所增加。作者進(jìn)一步指出當(dāng)兩次燃料噴射質(zhì)量相等時(shí)NOx和soot可以同時(shí)達(dá)到最優(yōu)。Hanson等11在一臺(tái)重型柴油機(jī)上使用汽油實(shí)現(xiàn)PPC燃燒，采用3:7的燃油質(zhì)量噴射比例(第一次比第二次)也發(fā)現(xiàn)了NOx和soot此消彼長(zhǎng)的關(guān)系：隨著第二次燃料噴射質(zhì)量的增加，NOx會(huì)降低，而HC、soot會(huì)增加。這與Manente等的結(jié)論一致?；趯?duì)于燃料的探索，Lv 等12采用進(jìn)氣道噴射PRF燃料，缸內(nèi)直噴正庚烷的模式詳細(xì)研究了噴油比例、時(shí)間、燃料特性、化學(xué)計(jì)量比等參數(shù)對(duì)燃燒和排放影響，并發(fā)現(xiàn)使用進(jìn)氣道噴射PFR50燃料并且適當(dāng)提高預(yù)混比例有助于提高燃油經(jīng)濟(jì)性和降低NOx。然而，最近的實(shí)驗(yàn)也發(fā)現(xiàn)部分預(yù)混燃燒模式向大負(fù)荷的拓展受到NOx和CO不能同時(shí)降低的制約13,14。另一方面，激光誘導(dǎo)熒光法通過(guò)測(cè)量甲醛、羥基等關(guān)鍵中間組分15，定量化化學(xué)計(jì)量比和燃料活性分布16，使得燃燒過(guò)程進(jìn)一步可視化，但是更具體的燃燒過(guò)程和更多組分的測(cè)量是激光測(cè)試手段和臺(tái)架實(shí)驗(yàn)所不能完成的。為了更進(jìn)一步地認(rèn)識(shí)PPC燃燒中多尺度的物理、化學(xué)過(guò)程，優(yōu)化放熱率，降低污染物排放，需要借助數(shù)值模擬的手段。

Jia和Xie17發(fā)展了一個(gè)包含內(nèi)核心區(qū)、邊界層區(qū)等在內(nèi)的六區(qū)模型研究了庚烷在HCCI發(fā)動(dòng)機(jī)中的燃燒和排放特性，使用了詳細(xì)的庚烷化學(xué)動(dòng)力學(xué)機(jī)理。研究表明該多區(qū)模型能夠較準(zhǔn)確地預(yù)測(cè)HC、CO和NO的排放。由于是分區(qū)計(jì)算，簡(jiǎn)化了每個(gè)區(qū)域中不同網(wǎng)格間的傳熱傳質(zhì)，沒有描述其詳細(xì)的燃燒過(guò)程。Wang等18使用三維數(shù)值模擬軟件FIRE耦合簡(jiǎn)化的化學(xué)動(dòng)力學(xué)模型(89組分的異辛烷機(jī)理)研究了缸內(nèi)直噴高辛烷值燃料的PPC燃燒模式，發(fā)現(xiàn)第二次噴射的燃料首先形成富燃區(qū)，隨后引燃周圍較均勻的第一次噴射的燃料，兩種不同燃燒模式的疊加使得放熱更加平穩(wěn)，有利于燃燒相位的控制。不同的數(shù)值方法都可以在一定程度上揭示發(fā)動(dòng)機(jī)燃燒性能，但是，直接數(shù)值模擬(DNS)是研究發(fā)動(dòng)機(jī)內(nèi)湍流/化學(xué)相互作用最準(zhǔn)確的方法。目前，發(fā)動(dòng)機(jī)內(nèi)的DNS研究主要采用二維(2D)計(jì)算，Chen19，Bansal和Im20，Yoo21，El-Asrag 和Ju22，Luong23等研究了帶有溫度和濃度不均勻性的氫氣或者庚烷等燃料的著火過(guò)程，發(fā)現(xiàn)對(duì)于單階段放熱的燃料氫氣，其溫度不均勻性能夠有效地改善HCCI發(fā)動(dòng)機(jī)放熱率過(guò)快的問(wèn)題，濃度不均勻性也可以在一定程度上緩解過(guò)高的壓升率，但是作用要小于溫度分層，原因是滯燃期主要受化學(xué)反應(yīng)速率影響，后者受溫度影響很大，而化學(xué)計(jì)量比對(duì)滯燃期的敏感性主要與燃料特性有關(guān)。對(duì)于具有兩階段著火過(guò)程的庚烷，其放熱率的改善與初始平均溫度和溫度分層的大小有關(guān)，需要具體分析23。Zhang等24通過(guò)一系列2D DNS計(jì)算定性地研究了PPC發(fā)動(dòng)機(jī)中兩次燃料噴射質(zhì)量之比對(duì)后續(xù)著火過(guò)程和污染物生成的影響規(guī)律；發(fā)現(xiàn)當(dāng)兩次噴射質(zhì)量相等時(shí)可同時(shí)實(shí)現(xiàn)NOx和UHC最少18，第二次噴油質(zhì)量越多，滯燃期和CA50推遲、壓升率降低，CO生成越多，反之則NOx生成量增加，形成了CO和NOx此消彼長(zhǎng)的關(guān)系。然而，由于采用多次噴射，PPC發(fā)動(dòng)機(jī)內(nèi)的濃度分層與HCCI發(fā)動(dòng)機(jī)中的分層完全不同，特別是在近上止點(diǎn)附近，有大范圍燃料濃區(qū)存在，因此燃燒過(guò)程也會(huì)與HCCI發(fā)動(dòng)機(jī)中的不同，而且燃料與湍流的相互作用對(duì)混合氣形成以及PPC燃燒過(guò)程至關(guān)重要，采用2D DNS不能準(zhǔn)確地描述湍流對(duì)分層著火的影響，因此要認(rèn)識(shí)PPC燃燒的本質(zhì)，應(yīng)該使用三維(3D)DNS解析小尺度湍流渦結(jié)構(gòu)以及三維火焰鋒面的結(jié)構(gòu)。

本研究基于作者前期使用2D DNS對(duì)PPC發(fā)動(dòng)機(jī)中兩次燃料噴射質(zhì)量之比對(duì)著火滯燃期和污染物趨勢(shì)的預(yù)測(cè)24，進(jìn)一步計(jì)算了當(dāng)兩次噴射質(zhì)量相等時(shí)的3D DNS算例。詳細(xì)描述了三維反應(yīng)鋒面的結(jié)構(gòu)以及存在的不同燃燒模式；更深入地分析了PPC燃燒的物理過(guò)程，湍流/化學(xué)相互作用，定量化火焰鋒面的曲率特性，包括平均曲率和高斯曲率的統(tǒng)計(jì)值以及它們隨時(shí)間的變化。

2　數(shù)值方法和算例設(shè)置

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

描述PPC發(fā)動(dòng)機(jī)內(nèi)湍流燃燒的過(guò)程可使用連續(xù)方程、N-S方程以及組分和能量的輸運(yùn)方程。本文采用in-house DNS程序，使用笛卡爾坐標(biāo)系下正交均勻網(wǎng)格，其數(shù)值格式的時(shí)間/空間精度通過(guò)了網(wǎng)格/時(shí)間步的精度測(cè)試，并且已經(jīng)在多種燃燒器中使用25－27。計(jì)算采用低馬赫數(shù)假設(shè)，忽略聲波的壓縮效應(yīng)，計(jì)算的時(shí)間步長(zhǎng)可以遠(yuǎn)遠(yuǎn)大于一般的可壓流求解器使用的時(shí)間步長(zhǎng)，同時(shí)保證系統(tǒng)的穩(wěn)定性。計(jì)算采用理想氣體假設(shè)，忽略體積力、熱/質(zhì)擴(kuò)散效應(yīng)和熱傳導(dǎo)輻射，組分的質(zhì)量擴(kuò)散使用Fickian定律；時(shí)間積分項(xiàng)采用二階對(duì)稱分裂算子方法，即在兩個(gè)半步的擴(kuò)散和對(duì)流時(shí)間步之間，積分一個(gè)完整的化學(xué)反應(yīng)速率時(shí)間步；時(shí)間步的推進(jìn)采用二階龍哥庫(kù)塔結(jié)合二階Adam-Bashforth(AB2)格式；所有的空間差分采用六階中心格式，除了組分和能量方程中的對(duì)流項(xiàng)采用五階加權(quán)非震蕩WENO格式。數(shù)值方法的具體描述參考文獻(xiàn)25。

為了模擬真實(shí)PPC發(fā)動(dòng)機(jī)中的壓縮/膨脹過(guò)程，質(zhì)量方程、動(dòng)量方程和能量方程分別加入了源項(xiàng)。經(jīng)過(guò)理論推導(dǎo)，如果連續(xù)方程以守恒形式求解，此壓縮源項(xiàng)將會(huì)以一種空間均勻分布的方式加入到連續(xù)方程；如果將其它的輸運(yùn)方程以非守恒形式(物質(zhì)導(dǎo)數(shù)形式)求解，意味著速度、組分質(zhì)量分?jǐn)?shù)和溫度以跟隨拉格朗日粒子的形式運(yùn)動(dòng)而不受影響，那么輸運(yùn)方程中的壓縮源項(xiàng)對(duì)于方程是沒有貢獻(xiàn)的。最終，連續(xù)方程中的壓縮源項(xiàng)只與平均密度的變化有關(guān)，可以直接通過(guò)發(fā)動(dòng)機(jī)曲柄連桿關(guān)系式得到，其幾何參數(shù)將在下一節(jié)介紹。更詳細(xì)的理論推導(dǎo)可參見文獻(xiàn)28。

2.2初始設(shè)置

本研究的計(jì)算條件為0.614 mm×0.614 mm× 0.614 mm的正方形容器，四面采用周期性邊界條件，均勻網(wǎng)格，網(wǎng)格數(shù)512×512×512，網(wǎng)格精度是1.2 μm。計(jì)算的時(shí)間步長(zhǎng)根據(jù)CFL準(zhǔn)則(Umax× Δt/Δx

其中，k0=2π/l0，為積分長(zhǎng)度尺度；u?為湍流強(qiáng)度，u?=(2K/3)1/2；K為湍動(dòng)能。

初始隨機(jī)的濃度分層也通過(guò)能量譜給定并疊加在平均濃度場(chǎng)之上，如圖1(a)所示，中心區(qū)域?yàn)榈诙稳剂蠂娚鋮^(qū)，濃度較高。根據(jù)實(shí)驗(yàn)中的已知參數(shù)：EGR與空氣的質(zhì)量比為1:1，第二次與第一次燃料噴射質(zhì)量之比為1.02，混合物總的化學(xué)計(jì)量比為0.695，則中心處二次燃料形成的球形半徑、化學(xué)計(jì)量比以及周圍區(qū)域化學(xué)計(jì)量比分布由關(guān)系式得到27：φ(x,y,z)=φ1+φ0e－br(x,y,z)/R，其中，r(x, y,z)是任一點(diǎn)距離計(jì)算區(qū)域中心點(diǎn)的長(zhǎng)度；R是計(jì)算區(qū)域長(zhǎng)度的一半；b是中心區(qū)域高濃度燃油球形體積的半徑。最終得到分布參數(shù)b=6.12，中心區(qū)域φ0=72φ1，其中φ1=0.455。第二次噴油之后最大的化學(xué)計(jì)量比為73φ1=33.2。從圖1(a)中看出，二次燃料噴射形成的中心濃區(qū)的球形半徑約為正方體邊長(zhǎng)的六分之一。關(guān)于化學(xué)計(jì)量比分布的具體推導(dǎo)可參考文獻(xiàn)27。

溫度分布與燃料濃度分布成一定的比例關(guān)系，PPC發(fā)動(dòng)機(jī)中的大渦模擬計(jì)算30顯示噴油結(jié)束后燃料蒸汽與溫度滿足二次多項(xiàng)式：T=1400Yf2－1600Yf+1200。其中，T為溫度，Yf為正庚烷、異辛烷質(zhì)量分?jǐn)?shù)之和。根據(jù)這一關(guān)系式得到初始溫度場(chǎng)中，第一次燃料噴射區(qū)溫度范圍在1075至1175 K之間，第二次燃料噴射區(qū)的溫度范圍在850 K至1075 K之間。雖然比實(shí)際PPC發(fā)動(dòng)機(jī)缸內(nèi)對(duì)應(yīng)的燃料噴射溫度略高，但是計(jì)算中仍然可以撲捉到一些低溫反應(yīng)過(guò)程，并且使得三維DNS計(jì)算量保持在可以承受的范圍內(nèi)。圖1(b)顯示了初始渦量場(chǎng)(λ2)的三維結(jié)構(gòu)。λ2定義為S2+Ω2=SijSji+ΩijΩji中對(duì)應(yīng)的三個(gè)特征值的中間值，其中，Sij和Ωij分別為變形率張量和旋轉(zhuǎn)率張量。λ2方法通常用來(lái)表示旋轉(zhuǎn)渦量的中心31。圖1中顯示的Z軸為壓縮方向：即發(fā)動(dòng)機(jī)壓縮沖程方向。由數(shù)值方法中的討論可知，此壓縮過(guò)程不引起網(wǎng)格的變化，僅僅改變物質(zhì)密度。

圖1　(a)三維(3D)空間中初始異辛烷質(zhì)量分?jǐn)?shù)的分布；(b)λ2渦量場(chǎng)等值面Fig.1　(a)Initial field of iso-octane mass fraction in three-dimensional(3D)space;(b)λ2-vortex iso-surface

發(fā)動(dòng)機(jī)中壓縮/膨脹過(guò)程是通過(guò)標(biāo)準(zhǔn)的曲柄連桿關(guān)系式得到。發(fā)動(dòng)機(jī)的幾何尺寸和設(shè)計(jì)依照Sandia D12 PPC發(fā)動(dòng)機(jī)的參數(shù)。發(fā)動(dòng)機(jī)幾何參數(shù)為：曲柄桿長(zhǎng)615 mm，缸徑127.5 mm，行程154 mm，壓縮比17.1，轉(zhuǎn)速1100 r·min－1。EGR中各組分的質(zhì)量分?jǐn)?shù)分別為：水0.0426，二氧化碳0.0923，氧氣0.115，氮?dú)?.727。計(jì)算中采用了簡(jiǎn)化的PRF化學(xué)動(dòng)力學(xué)機(jī)理，包含33個(gè)組分、38步基元反應(yīng)，也包含了PRF70燃燒過(guò)程中的低溫和高溫反應(yīng)機(jī)理32。本研究中，一個(gè)三維DNS的PPC算例使用4096核酷睿CPU(主頻為2.3 GHz)并行運(yùn)轉(zhuǎn)30天，共消耗350萬(wàn)計(jì)算機(jī)小時(shí)。

3　結(jié)果與討論

3.1PPC燃燒過(guò)程

圖2顯示了三維PPC燃燒過(guò)程中體積平均的總放熱率(HR)、壓力和CO質(zhì)量分?jǐn)?shù)隨時(shí)間的變化。可以看出，在初始時(shí)刻0.05 ms之前放熱率為負(fù)，表示反應(yīng)一開始有吸熱現(xiàn)象，主要是因?yàn)檎榈牡蜏鼗瘜W(xué)反應(yīng)32－34。之后，放熱率增加很快，最大值在0.188 ms，對(duì)應(yīng)了計(jì)算區(qū)域里最大的CO質(zhì)量分?jǐn)?shù)。隨后0.27 ms放熱率變得很低，主要的燃燒過(guò)程已經(jīng)結(jié)束，相應(yīng)的CO的消耗速率也降低。圖2中也顯示最大壓升率對(duì)應(yīng)了放熱率最大的時(shí)刻。

圖2　體積平均的放熱率(HR)、壓力和CO質(zhì)量分?jǐn)?shù)(w(CO))隨時(shí)間的變化Fig.2　Temporal evolution of total volume averaged heat release rate(HR),pressure,and mass fraction of CO(w(CO))

圖3　二維(x－y)中間平面上放熱率和庚烷質(zhì)量分?jǐn)?shù)的空間分布Fig.3　Contour plots of n-heptane mass fraction and HR in x－y space

為了進(jìn)一步分析PPC發(fā)動(dòng)機(jī)內(nèi)詳細(xì)的燃燒過(guò)程，圖3展示了三維計(jì)算結(jié)果中二維中心平面上庚烷的質(zhì)量分?jǐn)?shù)和放熱率在不同時(shí)刻的空間分布。在0.0137 ms，由于兩次噴油和湍流的共同作用，庚烷的分布有明顯的分層，初始時(shí)刻按照?qǐng)A形布置的第二次燃料噴射的區(qū)域已經(jīng)與周圍溫度較高、濃度較低的區(qū)域發(fā)生質(zhì)量和熱量的擴(kuò)散，并開始一系列低溫化學(xué)反應(yīng)。其中最先發(fā)生的是烷烴的脫氫加氧反應(yīng)，比如庚烷基C7H15O2最先出現(xiàn)在第二次噴射燃料濃區(qū)的邊緣和第一次燃料噴射的區(qū)域，而甲醛CH2O主要出現(xiàn)在第一次燃料噴射區(qū)域，主要是因?yàn)榈谝淮稳剂蠂娚鋮^(qū)溫度較高，庚烷基這樣的低溫反應(yīng)產(chǎn)物一旦生成很快會(huì)被進(jìn)一步分解成小分子烷基，比如CH2O；而第二次噴射的燃料溫度略低，需要等到溫度上升到一定范圍，KETO和C7H14OOH才開始分解為小分子烷烴和烯烴。由此也看到在0.074 ms，第一次燃料噴射區(qū)首先有部分熱量釋放。需要指出的是，燃料PRF70中庚烷摩爾質(zhì)量只有30%，因此在圖2中不能看到明顯的低溫放熱峰值。到了0.171 ms，周圍濃度分布較均勻的區(qū)域進(jìn)行著HCCI模式的燃燒，在第一次與第二次燃料噴射之間的部分，化學(xué)計(jì)量比在1附近，放熱最強(qiáng)。到了0.2 ms，燃料很濃的區(qū)域面積縮小至很小一部分，對(duì)應(yīng)著放熱率圖中被很薄的反應(yīng)層包圍的區(qū)域，這里一邊在進(jìn)行著自燃(大量的氧氣，C7H15O2，CH2O未顯示)一邊形成了預(yù)混反應(yīng)鋒面(在下一節(jié)會(huì)具體分析)。在兩次燃料噴射之間的部分由于是接近化學(xué)計(jì)量比燃燒，燃燒效率很高，很快就燃燒完全，放熱大大減少。在0.237 ms，庚烷幾乎消耗殆盡，在放熱率的圖像上出現(xiàn)了兩個(gè)明顯的薄反應(yīng)層，中間一個(gè)對(duì)應(yīng)了之前的預(yù)混反應(yīng)鋒面，外圍大的反應(yīng)層是由于CO和OH的相互擴(kuò)散形成的擴(kuò)散反應(yīng)層。具體過(guò)程是：燃料濃區(qū)形成大量CO得不到氧化，而周圍化學(xué)計(jì)量比較低的區(qū)域燃燒之后生成大量中間產(chǎn)物OH，并且有殘余氧氣，這些組分相互擴(kuò)散、氧化，形成了非預(yù)混的火焰鋒面。由此也看出，PPC燃燒過(guò)程非常復(fù)雜，存在三種不同的燃燒模式，即：在第一次燃料噴射之后形成的濃度較低的均質(zhì)壓燃模式的燃燒，緊接著是第二次燃料噴射后在中心區(qū)域形成的高濃度預(yù)混燃燒，最后是由于CO擴(kuò)散在兩次燃料噴射之間的區(qū)域形成的擴(kuò)散燃燒。目前，還沒有一種燃燒模型能夠同時(shí)準(zhǔn)確描述這三個(gè)過(guò)程，這也是今后工作的另一個(gè)方面，發(fā)展一個(gè)能夠耦合不同燃燒模式的PPC燃燒模型。

為了進(jìn)一步理解反應(yīng)區(qū)詳細(xì)的燃燒過(guò)程，圖4a和4b分別顯示了庚烷和CO的擴(kuò)散和化學(xué)反應(yīng)速率的一維空間分布，對(duì)應(yīng)了圖3放熱率中虛線A 和C處。擴(kuò)散速率和化學(xué)反應(yīng)速率分別定義為，其中，k為任意組分，這里取庚烷和CO，ρ為氣體密度。Dk和Rk的定義與Chen等19定義的budget term非常相像，只是不考慮未燃?xì)怏w密度和組分濃度梯度的變化；同樣，與Bansal和Im20定義的更加類似，其選定HO2為組分k的值，Da代表Dk和Rk二者的比值。這幾種方法目的都是為了描述反應(yīng)鋒面上是火焰?zhèn)鞑ミ€是自燃著火過(guò)程，根據(jù)燃料的不同可選取不同的組分。通常選取的是一維層流火焰結(jié)構(gòu)中與放熱率空間分布最吻合的那一組分，來(lái)代表反應(yīng)層。在本算例中采用兩種組分是因?yàn)閷?duì)于具有兩階段放熱的燃料庚烷，其主放熱階段的峰值對(duì)應(yīng)CO的生成，而在低溫放熱階段主要為庚烷的脫氫加氧過(guò)程，因此在這一時(shí)刻選用庚烷來(lái)確定反應(yīng)鋒面。例如在0.171 ms沿著虛線A，庚烷的擴(kuò)散速率和化學(xué)反應(yīng)速率幾乎在同一個(gè)量級(jí)，并且后者只在圖3中放熱率對(duì)應(yīng)的薄反應(yīng)層區(qū)為非零值，包圍著燃料的濃區(qū)，也說(shuō)明這里有火焰鋒面存在，這一現(xiàn)象與虛線B處相同，可以看到中心區(qū)域更明顯的薄反應(yīng)層區(qū)，所以圖4略去了B處對(duì)應(yīng)的庚烷擴(kuò)散與化學(xué)反應(yīng)速率分布。到了0.237 ms虛線C處，對(duì)應(yīng)著薄反應(yīng)層之間(0.3－0.4 mm)處，CO的化學(xué)反應(yīng)速率為正，說(shuō)明CO會(huì)不斷生成，在周圍區(qū)域CO的反應(yīng)速率為負(fù)，說(shuō)明CO在不斷消耗。由于第一次和第二次燃料噴射之間的區(qū)域?yàn)榻瘜W(xué)計(jì)量比燃燒，氧氣幾乎消耗完，大量的CO需要依賴后期周圍氧氣和OH的擴(kuò)散才足以燃燒完，由此可見，如果要減少最終CO的生成，需要增大湍流強(qiáng)度促進(jìn)混合。

圖4　庚烷和一氧化碳組分的擴(kuò)散項(xiàng)與化學(xué)反應(yīng)項(xiàng)分別在0.171 ms,0.237 ms的一維空間的分布，對(duì)應(yīng)于圖3放熱率二維分布中的虛線A(a)和C(b)Fig.4　Diffusion and reaction in terms of n-heptane and CO at 0.171 ms and 0.237 ms in one dimensional distribution,corresponding to dashed linesA(a)and C(b)in Fig.3,respectively

圖5　溫度與化學(xué)計(jì)量比(φ)在三個(gè)不同時(shí)刻的散點(diǎn)圖Fig.5　Scatter plot of temperature and equivalence ratio (φ)at three different time instances

圖5定量化溫度和化學(xué)計(jì)量比在不同時(shí)刻的散點(diǎn)分布。在初始時(shí)刻，化學(xué)計(jì)量比的分布范圍很寬(φ=0.6－6)，對(duì)應(yīng)的溫度范圍在1200至820 K之間。到了0.179 ms，化學(xué)計(jì)量比在2以下的區(qū)域溫度都有所增加，從圖3中可以看到，此時(shí)大部分放熱來(lái)自于第一次和第二次燃料噴射之間的區(qū)域(化學(xué)計(jì)量比在1附近)，同時(shí)，周圍HCCI模式的燃燒也在進(jìn)行；而高濃度區(qū)(φ>2)僅僅有一些低溫化學(xué)反應(yīng)，溫度升高非常緩慢。隨著反應(yīng)的進(jìn)一步進(jìn)行到0.237 ms時(shí)，φ<2的區(qū)域溫度都在1600 K以上，特別是化學(xué)計(jì)量比燃燒后溫度最高達(dá)到2100 K；在HCCI燃燒模式的區(qū)域(φ<0.9)，溫度在1600－2000 K之間；在φ>2的區(qū)域，自燃著火伴隨著火焰?zhèn)鞑?，溫度增加，但是低?600 K。由此推斷，大量的NOx會(huì)在兩次燃料噴射之間的區(qū)域φ=1附近出現(xiàn)，而CO會(huì)在燃料濃區(qū)φ>2的區(qū)域聚集，這與文獻(xiàn)13的結(jié)論類似，即UHC和CO出現(xiàn)在第二次燃料噴射的濃區(qū)。

從之前的分析可以看出，PPC燃燒模式非常復(fù)雜，其三維火焰鋒面的結(jié)構(gòu)就更為復(fù)雜，圖6展示了使用Marching cube算法35得到的三維空間內(nèi)溫度等值面(T=1600 K)隨時(shí)間的變化過(guò)程。可以看到，三維空間中反應(yīng)鋒面上有嚴(yán)重的扭曲和褶皺。從0.171 ms，計(jì)算區(qū)域中存在一個(gè)較小的著火區(qū)，對(duì)應(yīng)著圖3中第二次燃料噴射的濃區(qū)，此時(shí)周圍也在發(fā)生較均勻的自燃著火，但是溫度未達(dá)到1600 K，沒有捕捉到相應(yīng)的反應(yīng)鋒面。從不同的角度甚至可以看到散布著一些孤立的火核(該圖未顯示)。隨后在0.181 ms，出現(xiàn)了更大的的反應(yīng)鋒面即周圍第一次燃料噴射區(qū)域HCCI模式的著火，如果仔細(xì)看能發(fā)現(xiàn)在中心區(qū)域存在上一時(shí)刻出現(xiàn)的較小的火核。到了0.189 ms兩個(gè)溫度等值面共存的現(xiàn)象更加明顯：外圍大的反應(yīng)鋒面不斷向外傳播，中間小的反應(yīng)鋒面繼續(xù)向內(nèi)收縮。0.2 ms時(shí)，第一次噴射的燃料基本消耗完，只剩下中心處燃料濃區(qū)的一小部分庚烷，對(duì)應(yīng)了圖3中0.2 ms時(shí)庚烷和放熱率的分布。可以推測(cè)，0.23 ms以后，1600 K的溫度等值面會(huì)越來(lái)越小，直到燃料完全消耗。

圖6　溫度等值面隨時(shí)間變化的三維圖像Fig.6　Temporal evolution of 3D temperature iso-surface

3.2PPC反應(yīng)鋒面的曲率特性

為了更進(jìn)一步理解湍流和化學(xué)反應(yīng)對(duì)反應(yīng)鋒面變形產(chǎn)生的影響，圖7顯示了反應(yīng)鋒面上(T= 1600 K的等值面)高斯曲率(kg)和平均曲率(km)的散點(diǎn)分布，顏色代表聯(lián)合概率密度(PDF)(請(qǐng)參考網(wǎng)絡(luò)版)，虛線代表kg=0，表示二維空間上的圓柱形曲面。三維曲面上，平均曲率定義為兩個(gè)主曲率的算術(shù)平均：km=(k1+k2)/2，代表曲面是凸還是凹；高斯曲率定義為兩個(gè)主曲率的乘積：kg=k1×k2，正值(k1，k2同號(hào))表示鋒面為球形，負(fù)值(k1，k2異號(hào))表示鋒面為馬鞍形。0.171 ms時(shí)，中心區(qū)域|km|,|kg|趨于0，說(shuō)明大部分反應(yīng)鋒面比較平緩；km傾向于正值，意味著反應(yīng)前鋒凸向未燃區(qū)域，此時(shí)，火核正在向外傳播。圖中右上方的散點(diǎn)(kg>> 0,km>>0)表現(xiàn)出較大的曲率并形成球形的反應(yīng)鋒面，代表周圍燃料較均勻的區(qū)域中正在生成新的火核。圖中也可以看出，所有反應(yīng)鋒面上的散點(diǎn)均落在kg=km2這條拋物線上和其下方，這也是主曲率非實(shí)根的數(shù)學(xué)邊界；對(duì)于落在kg=km2拋物線上的點(diǎn)代表了此反應(yīng)鋒面為球形，并且兩個(gè)主曲率相等(km1=km2)；而kg<0的點(diǎn)代表馬鞍形的反應(yīng)鋒面，從圖中看出其概率很小。隨著時(shí)間的推移，馬鞍形反應(yīng)鋒面增多，代表火核間的交匯和融合(對(duì)比0.2 ms)，但是球形鋒面仍占主要部分。這一時(shí)刻中kg和km的分布與文獻(xiàn)中對(duì)于HCCI點(diǎn)火過(guò)程的描述一致36。到了0.2 ms，此時(shí)周圍HCCI模式的燃燒基本結(jié)束，留下中間混合氣較濃、代表第二次噴射的燃料。從圖3、圖6得知，這里形成了一個(gè)預(yù)混火焰，反應(yīng)鋒面一邊自燃一邊以預(yù)混形式傳播，燃料的不斷消耗導(dǎo)致火焰鋒面收縮、正曲率減小。可以推斷，隨著燃料的耗盡，中間濃區(qū)的反應(yīng)鋒面越來(lái)越小，最終火核越來(lái)越小，其負(fù)曲率會(huì)進(jìn)一步增大。

圖7　反應(yīng)鋒面上平均曲率(km)與高斯曲率(kg)在0.171和0.2 ms時(shí)的關(guān)系Fig.7　Correlation of mean curvature and Gaussian curvature on the reaction surface at 0.171 and 0.2 ms

圖8　不同時(shí)刻下反應(yīng)鋒面上平均曲率(km)的PDF分布Fig.8　PDF distribution in terms of kmat four time instances

圖8定量化平均曲率在不同時(shí)刻的概率密度函數(shù)分布。這里選取的反應(yīng)鋒面仍然是T=1600 K時(shí)的溫度等值面。在0.135 ms，兩次燃料噴射之間的區(qū)域剛剛形成火核，km的峰值向右偏移，說(shuō)明此時(shí)火核向外傳播。如果是層流火焰鋒面?zhèn)鞑?，km完全為正并且有唯一值，等于其火焰鋒面對(duì)應(yīng)的半徑的倒數(shù)37。相比之下，圖8中不同時(shí)刻對(duì)應(yīng)的km都存在一定范圍的正值和負(fù)值，也說(shuō)明湍流的作用使火焰面發(fā)生扭曲，形狀復(fù)雜，曲率變化很大。到了0.171 ms，火核繼續(xù)向四周傳播(見圖6)，引燃周圍HCCI模式的燃燒以及中心區(qū)的預(yù)混火焰，導(dǎo)致最大平均曲率減小和負(fù)曲率的增加。0.2 ms時(shí)刻，負(fù)曲率繼續(xù)增大，表明中心濃區(qū)預(yù)混火焰不斷縮小，見圖6。圖8中km的概率分布形狀對(duì)稱并且隨時(shí)間變化其峰值移向負(fù)值，這一現(xiàn)象與文獻(xiàn)中對(duì)球形火焰?zhèn)鞑?duì)應(yīng)的km的PDF描述一致37,38。

4　結(jié)論

使用三維直接數(shù)值模擬方法研究了部分預(yù)混發(fā)動(dòng)機(jī)條件下PRF70/空氣/EGR混合物的燃燒過(guò)程。文章通過(guò)對(duì)二維中心平面上庚烷質(zhì)量分?jǐn)?shù)和放熱率的空間分布，以及三維反應(yīng)鋒面隨時(shí)間變化的分析，詳細(xì)展示了PPC燃燒的物理過(guò)程：包含第一次燃料噴射區(qū)域形成的均質(zhì)燃燒，第二次燃料噴射區(qū)域的預(yù)混燃燒以及燃燒后期兩次燃料噴射之間的擴(kuò)散燃燒。文中使用庚烷輸運(yùn)方程中的budget term，發(fā)現(xiàn)了火焰鋒面和自燃著火鋒面的位置；最后，不同時(shí)刻下高斯曲率和平均曲率的聯(lián)合概率密度函數(shù)說(shuō)明了PPC燃燒過(guò)程中主要存在球形的火焰鋒面，在燃燒后期km的峰值由正變?yōu)樨?fù)說(shuō)明HCCI模式的燃燒結(jié)束，反映了中心濃區(qū)預(yù)混火焰鋒面在變化。本研究加深了人們對(duì)PPC燃燒現(xiàn)象的理解，并對(duì)于進(jìn)一步發(fā)展PPC發(fā)動(dòng)機(jī)內(nèi)的湍流燃燒模型提供了必要的數(shù)據(jù)支持。

致謝：本工作的計(jì)算資源使用了瑞典SNIC中心的Lindgren服務(wù)器、Lund大學(xué)的Lunarc服務(wù)器以及中國(guó)天津天河一號(hào)服務(wù)器。感謝瑞典隆德大學(xué)白雪松教授和喻日新副教授對(duì)本工作的指導(dǎo)、建議和支持。

References

(1)Yao,M.F.;Zheng,Z.L.;Liu,H.F.Prog.Energy Combust.Sci. 2009,35(5),398.doi:10.1016/j.pecs.2009.05.001

(2)Noehre,C.;Andersson,M.;Johansson,B.;Hultqvist,A.SAE 2006,2006-01-3412.doi:10.4271/2006-01-3412

(3)Jia,M.;Li,Y.;Xie,M.Z.Energy Procedia 2015,66,3. doi:10.1016/j.egypro.2015.02.018

(4)Kim,K.;Kim,D.;Jung,Y.;Bae,C.Fuel 2013,109(7),616. doi:10.1016/j.fuel.2013.02.060

(5)Kalghatgi,G.;Risberg,P.;Angstrom,H.SAE 2007,2007-01-0006.doi:10.4271/2007-01-0006

(6)Borgqvist,P.;Tuner,M.;Mello,A.;Tunestal,P.;Johansson,B. SAE 2012,2012-01-1578.doi:10.4271/2012-01-1578

(7)Lv,X.C.;Chen,W.;Huang,Z.J.Combust.Sci.Technol.2005, 11(6).[呂興才,陳偉,黃震.燃燒科學(xué)與技術(shù),2005,11 (6),493.]doi:10.3321/j.issn:1006-8740.2005.06.002.

(8)Manente,V.;Zander,C.G.;Johansson,B.;Tunestal,P.SAE 2010,2010-01-2198.doi:10.4271/2010-01-2198

(9)Manente,V.;Tunestal,P.;Johansson,B.SAE 2010,2010-01-0871.doi:10.4271/2010-01-0871

(10)Manente,V.;Johansson,B.;Tunestal,P.J.Eng.Gas Turb. Power 2010,132(1),082802.doi:10.1115/ICES2009-76165.

(11)Hanson,R.;Splitter,D.;Reitz,C.R.SAE 2009,2009-01-1442. doi:10.4271/2009-01-1442

(12)Lv,X.C.;Shen,Y.;Zhang,Y.;Zhou,X.;Ji,L.;Yang,Z.; Huang,Z.Fuel 2011,90(5),2026.doi:10.1016/j. fuel.2011.01.026

(13)Yang,D.;Wang,Z.;Wang,J.X.;Shuai,S.J.Appl.Energy 2011, 88(9),2949.doi:10.1016/j.apenergy.2011.03.004

(14)Zhang,F.;Xu,H.;Rezaei,S.Z.;Kalghatgi,G.;Shuai,S.J.SAE 2012,2012-01-1138.doi:10.4271/2012-01-1138

(15)Ma,X.;Wang,Z.;Jiang,C.;Jiang,Y.;Xu,H.M.;Wang,J.X. Fuel 2014,134(9),603.doi:10.1016/j.fuel.2014.06.002

(16)Tang,Q.L.;Geng,C.;Li,M.K.;Liu,H.F.;Yao,M.F.Acta Phys.-Chim.Sin.2015,31(12),2269.[唐青龍,耿超,李明坤,劉海峰,堯命發(fā).物理化學(xué)學(xué)報(bào),2015,31(12),2269.] doi:10.3866/PKU.WHXB20151008

(17)Jia,M.;Xie,M.Z.J.Combust.Sci.Technol.2005,11(3),261.[賈明,解茂昭.燃燒科學(xué)與技術(shù),2005,11(3),261.] doi:10.3321/j.issn:1006-8740.2005.03.013

(18)Wang,Z.;Shuai,S.J.;Wang,J.X.;Tian,G.H.Fuel 2006,85 (12－13),1831.doi:10.1016/j.fuel.2006.02.013

(19)Chen,J.H.;Hawkes,E.R.;Sankaran,R.;Mason,S.D.;Im,H. G.Combust.Flame 2006,145(1),128.doi:10.1016/j. combustflame.2005.09.017

(20)Bansal,G.;Im,H.G.Combust.Flame 2011,158(11),2105. doi:10.1016/j.combustflame.2011.03.019

(21)Yoo,C.S.;Lu,T.;Chen,J.H.;Law,C.K.Combust.Flame 2011,158(9),1727.doi:10.1016/j.combustflame.2011.01.025

(22)El-Asrag,H.A.;Ju,Y.Combust.Flame 2014,161(1),256. doi:10.1016/j.combustflame.2013.07.012

(23)Luong,M.B.;Yua,G.H.;Lu,T.;Chung,S.H.;Yoo,C.S. Combust.Flame 2015,162(12),4566.doi:10.1016/j. combustflame.2015.09.015

(24)Zhang,F.;Yu,R.;Bai,X.S.Appl.Energy 2015,149,283.doi: 10.1016/j.apenergy.2015.03.058

(25)Yu,R.;Yu,J.F.;Bai,X.S.J.Comput.Phys.2012,231(16), 5504.doi:10.1016/j.jcp.2012.05.006

(26)Yu,J.F.;Yu,R.;Fan,X.Q.;Christensen,M.;Konnov,A.A.; Bai,X.S.Combust.Flame 2013,160(7),1276.doi:10.1016/j. combustflame.2013.02.011

(27)Zhang,F.;Yu,R.;Bai,X.S.Proc.Combust.Inst.2014,35(3), 2975.doi:10.1016/j.proci.2014.09.004

(28)Zhang,F.Direct Numerical Simulations of Low Temperature Combustion in IC Engine Related Conditions.Ph.D. Dissertation,Lund University,Sweden,2013

(29)Passot,T.;Pouquet,A.J.Fluid Mech.1987,118(1),441. doi:10.1017/S0022112087002167

(30)Solsjo,R.;Jangi,M.;Chartier,C.;Andersson,O.;Bai,X.S. Proc.Combust.Inst.2012,34(1),3031.doi:10.1016/j. proci.2012.06.111

(31)Day,M.S.;Bell,J.B.Combust.Theor.Model 2000,4(4),535. doi:10.1088/1364-7830/4/4/309

(32)Tsurushima,T.Proc.Combust.Inst.2009,32(2),2835. doi:10.1016/j.proci.2008.06.018

(33)Curran,H.J.;Gaffuri,P.;Pitz,W.J.;Westbrook,C.K.A Combust.Flame 1998,114(1－2),149.doi:10.1016/S0010-2180(97)

(34)Zheng,Z.L.;Yao,M.F.Transactions of CSICE 2004,22(3), 227.[鄭朝蕾,堯命發(fā).內(nèi)燃機(jī)學(xué)報(bào),2004,22(3),227.] doi:10.3321/j.issn:1000-0909.2004.03.006.

(35)Lorensen,W.E.Comput.Graph.1987,21(4),163. doi:10.1145/37401.37422.

(36)Yu,R.;Bai,X.S.Combust.Flame 2013,160(9),1716. doi:10.1016/j.combustflame.2013.03.025.

(37)Gashi,S.;Hult,J.;Jenkins,K.W.;Chakraborty,N.;Cant,S.; Kaminski,C.F.Proc.Combust.Inst.2005,30(1),809. doi:10.1016/j.proci.2004.08.003.

(38)Haq,M.H.;Sheppard,C.G.W.;Woolley,R.;Greenhalgh,D. A.;Lockett,R.D.Combust.Flame 2002,131(1－2),1. doi:10.1016/S0010-2180(0200383-8

Three-Dimensional Direct Numerical Simulation of Partially Premixed Combustion in Engine-Related Conditions

ZHANG FanYAO Ming-Fa*
(State Key Laboratory of Engines,Tianjin University,Tianjin 300072,P.R.China)

Three-dimensional direct numerical simulation is conducted to simulate the auto-ignition of the highoctane fuel PRF70 under partially premixed combustion(PPC)engine conditions.Askeletal primary reference fuel(PRF)chemical kinetic mechanism is adopted,including 33 species and 38 elementary reactions. Compression/expansion effects caused by piston motion,the real engine geometry,and the working conditions are considered.The simulation includes two injections,the first being used to form a relatively uniform base mixture and the second to form a stratified mixture and trigger the ignition.It is found that the combustion process in PPC engines is a rather complex combination of homogeneous combustion,rich premixed and diffusioncontrolled combustion.The region between the two injections is near stoichiometry,resulting in the formation of NOx,while abundant CO is retained in the region with equivalence ratio(φ)>2,which needs to diffuse to meet the oxidizer and burn in a diffusion flame.The marching cube method is used to extract the 3D flame surface and show the temporal evolution of the reaction front.Finally,the joint PDF of the Gaussian curvature(kg)and principle mean curvature(km)and temporal evolution of the probability density function(PDF)in terms of kmshow that kmplays a more important role and becomes negative as time evolves because of the consumption of rich premixed flame in the center.

Direct numerical simulation;Partially premixed combustion;PRF70;Combustion characteristics;PPC engine

January 14,2016;RevisedApril 19,2016;Published on Web:April 22,2016.

O641；TK421.2

10.3866/PKU.WHXB201604223

*Corresponding author.Email:y_mingfa@tju.edu.cn;Tel:+86-22-23736295.

The project was supported by the Funds for International Cooperation and Exchange of National Natural Science Foundation of China (51320105008)and National Natural Science Foundation of China(51506146).

國(guó)家自然科學(xué)基金國(guó)際(地區(qū))合作與交流項(xiàng)目(51320105008)以及國(guó)家自然科學(xué)基金青年基金(51506146)資助

?Editorial office ofActa Physico-Chimica Sinica

［Article］