李 維，王 穩(wěn)，黃勇成，汪 映

(1.西京學(xué)院 汽車學(xué)院，陜西 西安710123;2.西安交通大學(xué) 能動(dòng)學(xué)院，陜西 西安710049)

采用均質(zhì)充量壓縮燃燒(HCCI)方式可使發(fā)動(dòng)機(jī)在中低負(fù)荷具有高的熱效率和低的PM和NOx排放，但HCCI還未從根本上解決著火時(shí)刻和燃燒速度的控制問(wèn)題，因此HCCI燃燒技術(shù)目前只局限在內(nèi)燃機(jī)的中低負(fù)荷范圍內(nèi)應(yīng)用［1］。采用在同一循環(huán)中使用預(yù)混進(jìn)氣和缸內(nèi)直噴的復(fù)合燃燒(PCCIDI)方式可縮短滯燃期、減少擴(kuò)散燃燒部分的油量、降低最高燃燒溫度與壓力［2］。采用PCCI-DI燃燒方式不但可進(jìn)一步提高二甲醚發(fā)動(dòng)機(jī)的熱效率和降低NOx排放，而且發(fā)動(dòng)機(jī)仍可在較寬廣的轉(zhuǎn)速和負(fù)荷下運(yùn)行［3］。

目前，對(duì)二甲醚噴霧特性已經(jīng)有了一些較深入的研究［4-5］。但利用二甲醚發(fā)動(dòng)機(jī)缸內(nèi)燃燒過(guò)程的數(shù)值模擬來(lái)闡述二甲醚發(fā)動(dòng)機(jī)性能的研究還不是很多。由二甲醚與柴油的理化性質(zhì)差異較大，因此直接使用原來(lái)適用于柴油燃料的壓燃式發(fā)動(dòng)機(jī)過(guò)程的數(shù)值模擬并不完全適合用來(lái)研究二甲醚在發(fā)動(dòng)機(jī)缸內(nèi)的燃燒，但二者作為壓縮燃燒的燃燒過(guò)程還有著許多相似之處，可以在原有方法的基礎(chǔ)上進(jìn)行一些改造來(lái)研究二甲醚的燃燒過(guò)程。為此，筆者采用數(shù)值模擬的方法開(kāi)展PCCI-DI二甲醚發(fā)動(dòng)機(jī)燃燒過(guò)程的研究。

為了對(duì)空氣運(yùn)動(dòng)、噴油特性和燃燒室?guī)缀涡螤畹纫蛩刂g的配合進(jìn)行數(shù)學(xué)分析和模型，以達(dá)到其最佳匹配實(shí)現(xiàn)高效低污染燃燒，要應(yīng)用多維模型［4］。KIVA-Ⅲ是專門為模擬內(nèi)燃機(jī)工作過(guò)程而設(shè)計(jì)的軟件，筆者以KIVA-Ⅲ程序基礎(chǔ)上，通過(guò)擴(kuò)展和修改其中的噴霧、燃燒等模型來(lái)進(jìn)行二甲醚發(fā)動(dòng)機(jī)PCCI-DI燃燒過(guò)程的三維數(shù)值模擬。

數(shù)值模擬過(guò)程中，氣相模型由缸內(nèi)氣體流動(dòng)的守恒偏微分方程組、湍流模型方程和狀態(tài)方程等控制方程構(gòu)成。湍流模型采用標(biāo)準(zhǔn)的κ-ε雙方程模型。本文僅就噴霧、燃燒模型做一介紹。

1 噴霧模型

1.1 噴霧方程

一般對(duì)噴霧的描述是采用概率統(tǒng)計(jì)的方法，液滴概率分布函數(shù)f定義如下:

表示在t時(shí)刻位置x上，單位體積內(nèi)，具有速度處于(v，v+dv)范圍內(nèi)，半徑處于(r，r+dr)范圍內(nèi)，溫度處于(Td，Td+dTd)范圍內(nèi)，球形偏移參數(shù)處于(y，y+dy)和范圍內(nèi)的液滴的概率意義上的個(gè)數(shù)。

液滴概率分布函數(shù)f的控制方程為:

式中:t是時(shí)間;v為速度;r為液滴尺寸;F=dv/dt為液滴加速度;R=dr/dt為液滴尺寸r隨時(shí)間變化的速率;Td表示液滴溫度;y表示液滴圓球度和分別為由于液滴的碰撞和破碎而產(chǎn)生的源項(xiàng)。

但液滴概率分布函數(shù)的變量非常多，直接求解其控制方程幾乎不可能。筆者采用隨機(jī)質(zhì)點(diǎn)方法，用離散的液滴來(lái)描述噴霧，并借助于蒙特卡洛方法，用一定數(shù)量的顆粒代表整群顆粒的采樣，每一個(gè)用于計(jì)算的顆粒被認(rèn)為代表一群具有相同特征的顆粒。于是，連續(xù)分布的液滴概率分布函數(shù)可離散為:

每一個(gè)顆粒p由Np個(gè)具有相同的位置xp、速度vp、直徑rp、溫度Td、振蕩參數(shù)yp和y·p的液滴組成。這些液滴運(yùn)動(dòng)軌跡相同，它們?cè)诳臻g中的運(yùn)動(dòng)受以下動(dòng)力方程控制:

液滴的加速度F由氣相阻力和重力決定:

式中:u為氣相速度;u'為氣相湍流速度;v為噴霧粒子速度;ρd為液滴密度;阻力系數(shù)CD是根據(jù)小球在空氣中運(yùn)動(dòng)的Stokes關(guān)系，并結(jié)合實(shí)驗(yàn)的修正而得到:

μair()可由式(6))求得:

式中:A1和A2為常數(shù)。

1.2 液滴破碎模型

KIVA程序中采用的TAB(Taylor Analogy Breakup)液滴破碎模型是Taylor將液滴振動(dòng)及變形與彈性質(zhì)量系統(tǒng)比擬的基礎(chǔ)上得到的。對(duì)黏度甚小的二甲醚而言，其噴霧貫穿度小于柴油，噴霧錐角比柴油大［4-5］。因此，TAB模型對(duì)黏度甚小的二甲醚不易得到準(zhǔn)確的結(jié)果。筆者采用Kelvin-Helmholtz(KH)Wave模型來(lái)進(jìn)行計(jì)算［6］。

假設(shè)在液滴表面施加一個(gè)極小的偏移量，即

式中:η'0是擾動(dòng)的初始振幅;?是波動(dòng)的增長(zhǎng)速率;k為波數(shù);R表示為括號(hào)內(nèi)復(fù)數(shù)函數(shù)的實(shí)數(shù)部分。列出連續(xù)方程和運(yùn)動(dòng)方程，并進(jìn)行求解，可得到波動(dòng)增長(zhǎng)速率和波長(zhǎng)的方程。其中波動(dòng)的最大增長(zhǎng)率Ω和相對(duì)應(yīng)的波長(zhǎng)Λ與氣液兩相物理特性的關(guān)系如下:

式中:We2=ρ2U2a/σ，是氣體的Weber數(shù);We1=ρ1U2a/σ ，是液體的 Weber 數(shù);Z=We0.51/Re1，是Ohnesorge數(shù)［Re1=Ua/v，是液體的Reynold數(shù)(U為液氣相對(duì)速度;a為液滴的半徑)］;T=ZWe0.52。

新生成的液滴半徑為:

式中:B0=0.61。

原液滴半徑的變化遵從以下關(guān)系:

式中:τ為破碎時(shí)間，由式(12)計(jì)算:

式中:B1為可調(diào)模型常數(shù)。

當(dāng)液滴每次破碎后，子液滴具有與父液滴相同的溫度和位置，速度矢量在父液滴速度方向上的分量u為父液滴的速度U，在其垂直平面上2個(gè)速度分量分別為cosφ，θ為噴霧錐角，φ 按統(tǒng)計(jì)規(guī)律在(0，2π)隨機(jī)選取。

1.3 液滴蒸發(fā)模型

假定液滴是球?qū)ΨQ的，液滴內(nèi)部溫度均勻，則根據(jù)液滴與周圍工質(zhì)之間的能量平衡可得液滴溫度Td隨時(shí)間t的變化率方程:

式中:L(Td)為液滴汽化潛熱;cl為液體的比熱;Qd為液滴表面單位面積上的導(dǎo)熱率。

液滴半徑變化率R由Frossling關(guān)系式給出:

液滴表面燃油蒸氣的質(zhì)量百分比Y*1的表達(dá)式為:

式中:W1為燃油蒸氣的分子量;W0為除燃油蒸氣外所有其它物質(zhì)的局部平均分子量;p為工質(zhì)總壓力;p0(Td)為在溫度Td下燃油蒸氣的分壓。

液滴表面單位面積上導(dǎo)熱率Qd表達(dá)式可由Ranz-Marshall關(guān)系式給出:

通過(guò)隱式迭代求解式(14)可以得到液滴溫度Td和尺寸r。

通過(guò)將t時(shí)刻處于位置x的所有液滴的質(zhì)量、動(dòng)量、和能量的變化率分別相加，就可以得到氣相控制方程中由于噴霧引起的的源項(xiàng)

由于DME的沸點(diǎn)甚低，密度又小，試驗(yàn)中噴嘴啟噴壓力在16.0 MPa，同時(shí)噴霧是噴在高溫的燃?xì)猱?dāng)中［2］，二甲醚的噴霧貫穿度減小，因此，不存在燃油碰壁現(xiàn)象，所以噴霧模型中不考慮噴霧碰壁的影響。

2 燃燒模型

燃燒過(guò)程的模擬是建立在包括45個(gè)組分的351個(gè)反應(yīng)的反應(yīng)機(jī)理之上［5-6］，可將這一機(jī)理編入KIVA程序代碼中用于DME模擬計(jì)算。但KIVA程序中所采用的湍流混合控制燃燒模型難以使用詳細(xì)的化學(xué)反應(yīng)機(jī)理，因此筆者采用部分?jǐn)嚢璺磻?yīng)器(Partially Stirred Reactor，PaSR)燃燒模型。

2.1 質(zhì)量守恒方程

系統(tǒng)內(nèi)發(fā)生的化學(xué)反應(yīng)可表達(dá)為:

式中:xs代表1 摩爾的組分s;ν'sr，ν″sr分別為正、逆反應(yīng)的化學(xué)計(jì)量系數(shù);Ns為參與反應(yīng)組分的數(shù)量;Nr為發(fā)生化學(xué)反應(yīng)的數(shù)量。

在連續(xù)方程中化學(xué)反應(yīng)的源項(xiàng)可以表示為:

化學(xué)反應(yīng)速率由式(18)給出:

式中:fm為組分m的摩爾濃度變化率;cs=ρs/Ws為第s種組分的摩爾濃度(其中ρs，和Ws分別為組分s的密度和分子量);分別為正向和逆向的反應(yīng)

式中:nr，Ar和Era分別為與空間因子、碰撞頻率和活化能相關(guān)的經(jīng)驗(yàn)參數(shù)。

2.2 參考組分方法

由于采用詳細(xì)的化學(xué)反應(yīng)機(jī)理模擬燃燒，而不同的化學(xué)反應(yīng)的特征時(shí)間相差很大，因此化學(xué)反應(yīng)源項(xiàng)的計(jì)算是剛性的(即強(qiáng)烈依賴于時(shí)間)。為此必須解決計(jì)算時(shí)間步長(zhǎng)選取這一問(wèn)題，于是引入了參考組分方法。

考慮下面基元反應(yīng)

式中:νi，ci分別為反應(yīng)計(jì)量系數(shù)和摩爾濃度。由于大多基元反應(yīng)的計(jì)量系數(shù)不大于2，為了簡(jiǎn)化，令νi=1。于是上述基元反應(yīng)的速率為:

速度常數(shù)，其表達(dá)式為修正的阿累尼烏斯形式:

式中:α=kf(c2+c1)+kb(c4+c3)。采用半隱式格式離散式(23)得:

將式(24)代入式(22)得到:

為了解決這個(gè)問(wèn)題，引入?yún)⒖冀M分方法。其中，參考組分定義為在化學(xué)反應(yīng)中其濃度變?yōu)樨?fù)數(shù)的危險(xiǎn)性最大的組分，參考組分被化學(xué)反應(yīng)所消耗，且具有最低的濃度。為了便于討論，假定參考組分cr為c1，則c2＞＞c1，于是有:

將式(26)代入式(25)得:

采用參考組分方法后，對(duì)于每個(gè)化學(xué)反應(yīng)自動(dòng)滿足化學(xué)平衡條件，不必將反應(yīng)速度很快的平衡反應(yīng)和反應(yīng)速度較慢的動(dòng)力反應(yīng)分開(kāi)來(lái)處理，消除了化學(xué)反應(yīng)中組分濃度變?yōu)樨?fù)數(shù)的危險(xiǎn)性，而且還給出了化學(xué)反應(yīng)的特征時(shí)間尺度τchem為:

雖然根據(jù)式(24)、式(28)可以求得化學(xué)反應(yīng)產(chǎn)生的源項(xiàng)，但沒(méi)有考慮湍流的影響。為了同時(shí)考慮化學(xué)反應(yīng)和湍流對(duì)燃燒的影響，可引入部分?jǐn)嚢璺磻?yīng)器燃燒模型。

部分?jǐn)嚢璺磻?yīng)器(Partially Stirred Reactor，PaSR)燃燒模型是Golovitchev等人在渦耗散概念(Eddy-Dissipation Concept，EDC)燃燒模型的基礎(chǔ)上發(fā)展而來(lái)，而且在EDC燃燒模型中結(jié)合使用詳細(xì)的化學(xué)反應(yīng)機(jī)理，自燃和燃燒的建模便沒(méi)有了本質(zhì)的區(qū)別，因此沒(méi)有必要對(duì)自燃過(guò)程再建立模型［7］。

在PaSR模型中，每個(gè)計(jì)算單元被分成反應(yīng)區(qū)和非反應(yīng)區(qū)，如圖1。

圖1 PaSR反應(yīng)器示意Fig.1 Schematic diagram of PaSR Reactor

反應(yīng)區(qū)視為完全攪拌反應(yīng)器(Perfectly Stirred Reactor，PSR)，占整個(gè)計(jì)算單元的體積分?jǐn)?shù)為k*。為了便于分析，反應(yīng)器中定義了以下3個(gè)平均濃度:進(jìn)口處組分平均濃度c0，反應(yīng)區(qū)組分濃度c和反應(yīng)器出口處組分濃度c1。在一個(gè)時(shí)間步長(zhǎng)τ內(nèi)，計(jì)算過(guò)程可按以下2步進(jìn)行:第1步，化學(xué)反應(yīng)使反應(yīng)區(qū)中組分濃度從初始的c0變化到c;第2步，在湍流作用下反應(yīng)區(qū)(組分濃度為c)與未反應(yīng)區(qū)(組分濃度為c0)混合使組分平均濃度變?yōu)閏1，第2步經(jīng)歷的時(shí)間為τmix。于是經(jīng)歷一個(gè)時(shí)間步長(zhǎng)τ后，計(jì)算單元中組分平均濃度為:

式(29)中次網(wǎng)格中反應(yīng)區(qū)濃度c為未知量，必須用單元中已知信息來(lái)表達(dá)。為此假定在計(jì)算步長(zhǎng)τ內(nèi)計(jì)算單元中組分濃度從c0變化到c1的平均速率與在混合時(shí)間τmix內(nèi)反應(yīng)區(qū)組分濃度從c變化到c1的平均速率相等，并與式(17)表達(dá)的化學(xué)反應(yīng)的源項(xiàng)相對(duì)應(yīng)，于是有:

由式(29)和式(30)可得:

式(30)表明計(jì)算時(shí)間步長(zhǎng)大于混合時(shí)間，當(dāng)k*趨近于1時(shí)整個(gè)計(jì)算單元為反應(yīng)區(qū)，而當(dāng)k*趨近于0時(shí)，整個(gè)計(jì)算單元沒(méi)有發(fā)生化學(xué)反應(yīng)。確定了反應(yīng)區(qū)體積分?jǐn)?shù)k*后，需消除未知量c以求得c1。

fm(c)在c1附近的Taylor展開(kāi)式為:

由式(29)得到c的表達(dá)式后代入式(32)并結(jié)合式(30)得:

由式(33)可得到PaSR模型中計(jì)算化學(xué)反應(yīng)源項(xiàng)的方程

2.2.1 湍流混合時(shí)間 τmix的確定

τmix表征的是PaSR反應(yīng)器中新鮮、未發(fā)生反應(yīng)的混合物與已燃?xì)怏w間混合過(guò)程的特征時(shí)間。對(duì)于考慮湍流與化學(xué)反應(yīng)相互作用的燃燒模型，正確選擇微混合時(shí)間至關(guān)重要。從較大的漩渦到分子水平，表征湍流時(shí)間尺度的參數(shù)有許多。Kj?ldman對(duì)Kolmogorov尺度，Taylor尺度以及 Kolmogorov和Taylor尺度的幾何平均等3種時(shí)間尺度進(jìn)行了研究，發(fā)現(xiàn)采用Taylor時(shí)間尺度結(jié)果最佳。本文也采用Taylor時(shí)間尺度計(jì)算微混合時(shí)間。對(duì)于κ-ε湍流模型有:

式中:Cmix為模型常數(shù)。

2.2.2 化學(xué)反應(yīng)特征時(shí)間尺度τchem的確定

化學(xué)反應(yīng)時(shí)間尺度有多種定義方式。本文中化學(xué)反應(yīng)時(shí)間尺度的定義基于參考組分cr參與反應(yīng)的時(shí)間(參考組分具有最短的壽命)，即:

由式(16)可得:

對(duì)于第r個(gè)反應(yīng)，其反應(yīng)速率為:

由于反應(yīng)速率是一個(gè)與所有組分濃度相關(guān)的高度非線性函數(shù)，為了簡(jiǎn)化計(jì)算，利用前面介紹的計(jì)算反應(yīng)速率的參考組分法，在參考組分cr附近對(duì)式(35)線性化后有:

式中:上標(biāo)0表示計(jì)算初始值，為已知量。對(duì)式(40)求導(dǎo)得:

為了獲得下一步長(zhǎng)的反應(yīng)速率，采用半隱式格式離散式(41)可得［8-9］:

其中:上標(biāo)1表示下一步長(zhǎng)的值。于是可以得到下一步長(zhǎng)的反應(yīng)速率和反應(yīng)時(shí)間尺度:

將第r個(gè)反應(yīng)的反應(yīng)速率即式(43)，代入PaSR方程［即式(34)］可得:

3 邊界條件

一般固體壁面附近的邊界層很薄，用實(shí)際計(jì)算網(wǎng)格難以求解，在湍流流場(chǎng)中，邊界可采用湍流壁面函數(shù)代替。

油膜的蒸發(fā)使氣相產(chǎn)生了垂直于壁面的速度，因而改變了壁面油膜上方的湍流邊界層。定性地講，有蒸發(fā)的湍流邊界層相對(duì)于無(wú)蒸發(fā)的邊界層對(duì)質(zhì)量、動(dòng)量、能量的輸運(yùn)有抑制作用。因此，必需建立一個(gè)適用于油膜蒸發(fā)的壁面函數(shù)才能較精確地描述邊界層上質(zhì)量、動(dòng)量、能量的輸運(yùn)。

推導(dǎo)壁面函數(shù)的過(guò)程中，作了以下2點(diǎn)假設(shè):

①整個(gè)輸運(yùn)是湍流擴(kuò)散和蒸發(fā)所引起的對(duì)流輸運(yùn)的總和，且與壁面法向坐標(biāo)無(wú)關(guān);②象無(wú)蒸發(fā)邊界層那樣，假設(shè)湍流擴(kuò)散率與壁面距離成正比。

［y+=(νl為層流動(dòng)力黏性);K為湍動(dòng)能;cμ=0.09;Scl，ScT分別為層流和湍流的Schmidt數(shù);κ =0.433，為 Karmann 常數(shù)］;Yv為在y+處燃油蒸汽質(zhì)量百分比;Yv0為在溫度為Ts下飽和燃油蒸汽質(zhì)量百分比。

壁面剪切應(yīng)力方程為:

邊界層熱流量方程為:

式中:Prl，PrT分別為層流和湍流Prandtl數(shù)。邊界層速度方程為:

其中:U為距壁面y處的流速切向分量;υ為分子運(yùn)動(dòng)黏度;B為壁面粗糙度常量;u*為剪切速度。

在計(jì)算湍流流動(dòng)時(shí)，湍動(dòng)能κ和湍流耗散率ε按距壁面距離為y處計(jì)算，邊界條件為:

4 結(jié)語(yǔ)

1)針對(duì)二甲醚噴霧貫穿度小，噴霧錐角大的特點(diǎn)，可采用WAVE破碎模型替代KIVA程序中的TAB模型;

2)采用部分?jǐn)嚢璺磻?yīng)器(Partially Stirred Reactor，PaSR)燃燒模型替代程序中所用的湍流混合控制燃燒模型可以使用詳細(xì)的化學(xué)反應(yīng)機(jī)理并將其編入KIVA程序中用于DME模擬計(jì)算。

［1］朱馳，劉圣華.二甲醚均質(zhì)充量壓燃發(fā)動(dòng)機(jī)排放特性的試驗(yàn)研究［J］.內(nèi)燃機(jī)學(xué)報(bào)，2004，22(1):51-54.ZHU Chi，LIU Shen-hua.Experiment investigation of DME HCCI engineemissions［J］. Transactions ofCSICE，2004，22(1):51-54.

［2］李維，周龍保，汪映，等.DME部分預(yù)混壓燃與缸內(nèi)直噴復(fù)合燃燒的研究［J］.重慶交通大學(xué)學(xué)報(bào):自然科學(xué)版，2010，29(3):480-483.LI Wei，ZHOU Long-bao，WANG Ying，et al.Experimental investigation on combustion characteristics of partial premixed charge compression ignition–direct injection engine fueled with dimethyl ether［J］.Journal of Chongqing Jiaotong University:Natural Science，2010，29(3):480-483.

［3］李維，汪映.DME預(yù)混合引導(dǎo)進(jìn)氣實(shí)現(xiàn)PCCI-DI燃燒的試驗(yàn)研究［J］.車用發(fā)動(dòng)機(jī)，2008，177(4):15-18.LI Wei，WANG Ying.Experimental investigation on the realization of PCCI-DI combustion by DME pilot from the intake pipe［J］.Vehicle Engine，2008，177(4):15-18.

［4］王賀武，周龍保.二甲醚噴霧特性的研究［J］.西安交通大學(xué)學(xué)報(bào):自然科學(xué)版，2001，35(9):1480-1483.WANG He-wu，ZHOU Long-bao.Study on spray characteristics of dimethyl ether［J］.Journal of Xian Jiaotong University:Natural Science，2001，35(9):1480-1483.

［5］馬愛(ài)兵，劉建江.對(duì)二甲醚的噴霧特性研究［J］.新技術(shù)新工藝，2008，11(4):99-100.MA Ai-bing，LIU Jian-jiang. Study on the characteristics of Dimethyl Ether(DME)general spary［J］.New Technology ＆New Process，2008，11(4):99-100.

［6］蔣德明.內(nèi)燃機(jī)燃燒與排放學(xué)［M］.西安:西安交通大學(xué)出版社，2001:479-500.

［7］Golovitchev V I，Nordin N，Jarnicki R ，et al.3-D Diesel Spray Simulations Using a New Detailed Chemistry Turbulent Combustion Model［EB/OL］.［2000-01-1891］.http://papers.sae.org/2000-01-1891.

［8］Fischer S L，Dryer F L，Curran H J.The reaction kinetics of dimethyl ether I:High-temperature pyrolysis and oxidation in flow reactors［J］.International Journal of Chemical Kinetics，2000，32(10):714-740.

［9］Curran H J，F(xiàn)ischer S L，Dryer F L.The reaction kinetics of dimethyl etherII:Low-temperature oxida-tion in flow reactors［J］.International Journal of Chemical Kinetics，2000，32(10):741-759.