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

張　帆　堯命發(fā)

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

?

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

張帆堯命發(fā)*

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

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

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

www.whxb.pku.edu.cn

1　引言

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

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

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

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

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

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

2.1數(shù)值計算方法

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

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

2.2初始設(shè)置

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

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

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

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

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

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

3　結(jié)果與討論

3.1PPC燃燒過程

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

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

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

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

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

圖4　庚烷和一氧化碳組分的擴散項與化學反應(yīng)項分別在0.171 ms,0.237 ms的一維空間的分布，對應(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　溫度與化學計量比(φ)在三個不同時刻的散點圖Fig.5　Scatter plot of temperature and equivalence ratio (φ)at three different time instances

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

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

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

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

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

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

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

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

4　結(jié)論

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

致謝：本工作的計算資源使用了瑞典SNIC中心的Lindgren服務(wù)器、Lund大學的Lunarc服務(wù)器以及中國天津天河一號服務(wù)器。感謝瑞典隆德大學白雪松教授和喻日新副教授對本工作的指導(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).[呂興才,陳偉,黃震.燃燒科學與技術(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ā).物理化學學報,2015,31(12),2269.] doi:10.3866/PKU.WHXB20151008

(17)Jia,M.;Xie,M.Z.J.Combust.Sci.Technol.2005,11(3),261.[賈明,解茂昭.燃燒科學與技術(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)燃機學報,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).

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

?Editorial office ofActa Physico-Chimica Sinica

［Article］