基于OpenFOAM的正十二烷噴霧燃燒數值模擬

李文剛, 肖隱利, 路易聘, 曹志博, 吳娟

(西北工業大學 動力與能源學院, 陜西 西安 710072)

隨著世界范圍內能源危機與環境污染的加重，高效能低排放已成為航空發動機設計需要考慮的重要因素。開發先進的燃燒控制方法是實現低排放的重要手段，為此，需要全面掌握燃燒特性，如點火延遲時間、火焰浮升高度和火焰穩定機理等。為了實現這一目標，ECN[1]為國際試驗和數值模擬的合作提供了一個開放可存取的數據庫和論壇。

ECN將正十二烷的噴霧燃燒稱為Spray A，由于常用柴油的化學成分復雜，其熱力學性質和化學性質很難獲得。因此，ECN剛開始數值模擬時采用易于操作且化學反應動力學已被廣泛驗證的正庚烷替代柴油。然而典型柴油碳鏈是由10～25個碳原子組成，因此正庚烷并不能較好地反映柴油的熱物理性質，尤其是柴油的沸騰和蒸發等性質，而這些性質對準確預測柴油噴霧燃燒的點火延遲時間和火焰浮升高度等參數至關重要。而正十二烷有更長的碳鏈，熱物理性質和輸運性質明確，更適合替代柴油燃料。近年來，許多學者和研究人員意識到開發正十二烷合適的化學反應機理的重要性，Luo等[2]開發并驗證了一個包含106組分420基元反應的正十二烷骨架機理，采用非定常雷諾時均法(URANS)驗證了射流噴霧燃燒；斯坦福大學的Narayanaswamy等[3]開發了一個包含257組分1 521基元反應的正十二烷反應機理；Ranzi等[4]開發了一個130組分2 395基元反應的正十二烷反應機理。清華大學的Yao等[5]開發并驗證了一個54組分269基元反應的正十二烷骨架機制。

Spray A有大量可用的試驗數據，包括環境溫度、環境密度、燃油注入壓力等參數變化相對應的試驗數據，ECN團隊利用高速攝像技術得到了正十二烷噴霧的瞬態過程，利用Mie散射和Rayleigh散射技術對無反應情況下氣體的混合情況進行定量測量[6]。

針對Spray A，不同的研究人員采用不同的湍流模型和燃燒模型進行了廣泛的數值研究。早些年，由于計算機資源的限制，大部分的數值研究是基于非定常雷諾時均法(URANS)下的求解器進行的，Errico等[7]采用與參考文獻[2]相同的正十二烷化學機理，比較了良好混合(well-mixed)燃燒模型和多代表交互小火焰(multiple Representative Interactive Flamelet，mRIF)燃燒模型分別在Spray A不同環境溫度和氧氣濃度條件下的情況，結果表明，采用well-mixed模型能夠準確地預測點火延遲，但高估了熱釋放率，預測的燃燒室壓力偏高，進而影響對污染物排放的準確預測，相比之下mRIF燃燒模型結果更好，與試驗結果有較高的一致性；Bhattacharjee和Haworth[8]研究了不同環境參數對Spray A的影響，此研究使用二維RANS結合輸運概率密度函數(TPDF)的方法來模擬反應擴散過程。

大渦模擬(LES)近些年被廣泛地應用于Spray A的數值研究中，Salehi等[9]首次將LES與多映射調節(multiple mapping conditioning,MMC)模型相結合應用于Spary A，發現Spray A火焰穩定發生在具有顯著湍流化學相互作用的區域，在靠近噴嘴出口液相會影響火焰；Wehrfritz等[10]基于OpenFOAM平臺，采用LES和火焰面生成流(flamelet generated manifold,FGM)方法對不同的環境氧氣濃度條件下的Spray A進行模擬，使用拉格朗日粒子追蹤(Lagrangian particle tracking)方法描述液滴顆粒，比較了參考文獻[3]和參考文獻[4]中2種正十二烷詳細機理的表現，發現2種化學機制在點火特性存在顯著差異，但對于穩定火焰來說，差別不大，該研究首次量化了燃油噴嘴出口與火焰浮升高度之間的冷焰結構(cool flame)，還發現在較低環境氣量濃度下，火焰浮升高度預測很好，21%的氧氣濃度下，對火焰浮升高度預測過高；Cheng等[11]采用LES方法對900 K和1 000 K 2個環境溫度的Spray A進行數值模擬，預測了Spray A兩級點火行為，結果表明，第一級著火發生在貧燃混氣中，第二級著火在富油混合氣中發生，貧燃混氣中的第一級點火促進了富燃料混合氣中的第一級和第二級點火，用自燃和火焰傳播與低溫點火耦合2種機理來解釋燃燒過程的火焰浮升高度和穩定性，2種機理相互競爭，關鍵取決于環境溫度。

雖然Chomiak和Karlsson[12]采用了部分攪拌反應器(PaSR)模型研究正庚烷的噴霧射流燃燒，但還未發現基于PaSR燃燒模型的正十二烷噴霧射流燃燒數值模擬報道，基于此，本文采用PaSR燃燒模型與一個包含54組分269基元反應的正十二烷骨架機制相結合，對Spray A進行數值研究。

1 物理模型和數值模擬方法

1.1 物理模型與網格劃分

Spray A噴霧燃燒器和截面結構示意圖如圖1所示，該燃燒器是由Sandia國家實驗室(SNL)構造的一個特征長度為108 mm的定容立方體，燃油噴嘴安裝在立方體上表面正中心處。數值模擬工況嚴格按照試驗條件設置，試驗具體條件如表1所示，有反應和無反應2種情況燃燒室氣體組分如表2所示。

圖1 SNL噴霧燃燒器及截面結構示意圖

表1 Spray A試驗條件H

表2 燃燒室氣體組分(體積分數) %

對計算網格的劃分,為了提高數值模擬計算效率,定容燃燒室簡化為一個長方體,其底面為邊長70 mm的正方形,長方體高108 mm。非均勻結構網格使用OpenFOAM網格生成字典blockMeshDict生成,利用網格加密字典refineMeshDict進行一次加密,加密區域為一個底面與燃油注入平面重合的長方體,底面為邊長16 mm的正方形,正方形中心與燃油噴嘴重合,高為78 mm,其二維截面如圖3所示。以便更好預測噴霧射流性能參數,總體六面網格數量約27.7萬,加密后最小網格尺寸為0.25 mm,縱橫比約2.3。

圖2 三維結構網格 圖3 三維網格加密區域

1.2 PaSR燃燒模型

PaSR模型通過求解組分輸運方程的方法來封閉湍流模型中的反應源項。該模型假設火焰結構比計算網格單元小得多,一個網格單元可分為反應區和非反應區,在非反應區,只發生反應物之間的混合,在反應區是一個完全攪拌反應區,在反應區中,化學時間尺度要遠遠大于湍流時間尺度。

假設整個計算網格單元是完全攪拌反應PSR,某反應物的初始平均濃度為c0,k,在時間t內反應物平均濃度從c0,k變化為ck,那么在時間t內,PSR網格單元中的平均反應速率如下表示

(1)

PaSR計算網格反應物平均濃度想要達到ck,需要額外的時間τ,那么平均反應濃度ck可以表示為

(2)

則PaSR計算網格反應物平均反應濃度c1,可表示為

c1=k*c+(1-k*)c0

(3)

式中：

(4)

k*為在一個PaSR計算網格中PSR反應區所占的體積分數,由于τ′與反應區與和非反應區的組分之間的混合過程相關,所以將τ′通常被稱為混合時間尺度τmix。在k-ε中可以通過公式(5)得到混合時間尺度

(5)

式中，Cmix是經驗確定的模型常數,其值為0.03。

1.3 Reitz-KHRT霧化模型

Reitz-KHRT模型是一種二次霧化模型,該模型認為經過初次霧化的較大液滴在噴霧場中同時受到2種不穩定波的作用,一種是K-H(Kelvin-Helmholtz)不穩定波,另外一種為R-T(Rayleigh-Taylor)不穩定波。K-H模型認為液滴受到空氣動力作用,引起K-H不穩定波增長導致氣液界面不穩定性,最終導致液滴分裂破碎;R-T模型認為是由于受到R-T波的不穩定性增長導致液滴分裂破碎(R-T波是指在氣液界面上沿著流動方向的波)。在子液滴霧化時,K-H模型和R-T模型是競爭關系,對一個將要破碎的液滴,判斷其增長最快的擾動波長度與液滴直徑的大小關系,若增長最快的擾動波長度小于液滴直徑,液滴發生RT破碎;否則判斷液滴韋伯數We是否大于臨界韋伯數Wecr(通常取Wecr值12.5),若是,則發生K-H破碎。

1.4 數值模擬方法

本研究的求解器基于開源CFD框架OpenFOAM下的標準求解器SprayFoam。SprayFoam是一個瞬態求解器,適用于可壓縮的湍流氣液兩相噴霧燃燒;噴霧采用歐拉-拉格朗日方法進行模擬,氣相在歐拉框架下使用非定常雷諾時均法(URANS)進行描述,液相在拉格朗日框架下使用拉格朗日粒子追蹤的方法(LPT)進行處理;湍流模型采用Realizek-ε模型,燃燒模型采用PaSR模型,化學機理采用一個包含54組分269基元反應的正十二烷的骨架機制[5]。計算過程中動態調整時間步長,約5×10-7s,最大庫朗數不超過0.25,確保計算的穩定性,模擬結束時間設置為1.5 ms,此時,無反應情況下液體穿透長度穩定在10 mm左右,有反應情況下火焰形成準穩態火焰結構,對模擬結果不會產生影響。

液相霧化子模型中燃油注入采用Blob模型,液滴受到作用力為球形阻力模型,初次霧化采用LISA模型,二次霧化采用Reitz-KHRT模型,蒸發模型為liquidEvaporation模型;燃油注入持續1.5 ms,總共注入18 000 000個液體顆粒包(liquid parcels),液體顆粒包由一群具有相同直徑和動量的液體顆粒組成;液滴的分布采用Rosin-Rammler模型,其中,設置液滴直徑最大值與噴嘴內徑相同,分布指數n設置為3;噴霧半角設置為10°。

壓力速度方程由PIMPLE算法耦合,PIMPLE算法結合了SIMPLE和PISO算法,適用于非定常求解器,確保了穩定性和準確性。擴散項采用二階中心差分格式,對流項采用二階TVD中心差分格式,時間項采用一階隱式Euler格式。

2 計算結果與分析

2.1 網格獨立性驗證

歐拉-拉格朗日框架下的噴霧數值計算結果與網格大小有密切關系,因此需要驗證網格獨立性。為了評估噴霧性能對網格的依賴性,對無反應情況的Spray A采用4種不同的二維軸對稱非均勻結構網格進行數值模擬,網格信息如表3所示。Nordin等[13]發現網格依賴性主要受徑向(即垂直于噴射方向)上的單元網格尺寸的影響,如果徑向上的單元網格尺寸太大,在液體霧化區域內得到的氣體速度分布將是平坦的,較細的單元網格可以更準確地表示速度場,然而,噴射方向上的單元網格尺寸對液體蒸汽穿透長度和火焰浮升高度有很大影響,因此在計算域網格劃分時,縱橫比不能太大。

表3 二維網格信息

4種網格對液體穿透長度和蒸汽穿透長度數值模擬結果與SNL試驗結果進行對比,如圖4所示,發現隨著網格尺寸減小,液體穿透長度增大,但當網格數量達到一定數量時,再細化網格卻不增加液體顆粒包時,不會獲得更好的模擬結果,只會增加計算成本。另外,對于液體穿透長度,mesh3和mesh4與試驗結果有很好的一致性,說明網格最小尺寸為0.25 mm時,數值結果已經能夠對噴霧性能做出較為準確的預測,這也驗證了對三維網格的劃分精度足夠,不會對模擬結果造成影響。

圖4 不同網格下液體和蒸氣穿透長度

2.2 無反應情況射流噴霧性能

為了評估當前的湍流和噴霧模型設置,本文針對無反應情況Spray A,對數值模擬的液體和氣體穿透長度、燃料混合分數分布結果與SNL試驗數據進行對比驗證。圖5分別顯示了液體和蒸氣穿透長度隨注入時間變化的試驗與模擬結果,在本文中,液體穿透長度定義為從噴嘴出口到累積注入燃油總質量95%的最遠液滴位置的軸向距離。與試驗數據相比,該噴霧模型很好地預測了液體和蒸氣穿透長度,只是在燃油注入早期液體和蒸氣穿透長度與試驗有很小的偏差,這是由于該噴霧模型不能很好反映燃油初次霧化造成的。

圖5 無反應情況下液體和蒸氣穿透長度 圖6 混合分數徑向分布

燃料混合分數的準確分布對噴霧燃燒的點火延遲時間和火焰浮升高度的預測至關重要。在2個軸向位置(25 mm和45 mm)提取燃料混合分數Z,模擬結果與試驗結果進行比較,如圖6所示?？梢钥吹侥M結果很好地預測了混合分數的徑向分布,與試驗結果有很好的一致性。只是在Y=25 mm處,預測結果略低于試驗結果,但偏差值在試驗誤差范圍內。

2.3 噴霧點火過程

通過分析有反應情況下的Spray A數值模擬結果研究燃油霧化后點火過程,圖7顯示了有反應情況下在燃燒室內壓力升高隨時間變化的曲線,其中藍色虛線表示本文模擬得到的燃燒室壓力升高隨時間的變化曲線,綠色實線表示將藍色曲線數據乘修正系數0.42得到的修正壓力升高曲線,這是由于在試驗中燃燒室是邊長為108 mm的立方體,而模擬中為了降低計算成本,只關注中心射流及其霧化的計算域,從而將計算域簡化為底面邊長為70 mm的正方形,高為108 mm的長方體,體積的變化導致密度改變,所以得到的壓力升高曲線斜率要遠遠大于試驗數據曲線斜率;因此,假設燃油蒸發后迅速擴散,即燃燒室內氣體密度均勻,根據連續方程和理想氣體狀態方程得到

(6)

式中:ρ1,A1分別為假設模擬燃燒室氣體密度和底面積;ρ,A為假設試驗燃燒室氣體密度和底面積。

從圖7中看到,燃油通過壓力旋流噴嘴不斷被注入燃燒室,由于高壓(150 MPa)和較高的環境溫度(900 K),燃油迅速霧化蒸發,在t<0.43 ms這一過程,只是發生燃油霧化蒸發和低溫反應,釋放熱量小,壓力幾乎沒有發生變化;在0.43 ms時燃油自燃,在隨后的時間中,由于燃燒放出大量熱量,燃燒室內壓力隨時間單調增加;在主燃燒階段,燃油蒸汽與空氣混合控制火焰的傳播;修正壓力升高曲線與試驗數據有很好的一致性,點火延遲時間也得到了很好的預測。

圖7 壓力升高曲線

圖8顯示了無反應情況和有反應情況下不同時刻混合分數分布云圖,其中黃色實線為當量混合分數等值線(Zst=0.045),對于有反應情況,0.4 ms是燃油高溫點火之前時刻,1.5 ms為燃油燃燒后一時刻。發現高溫點火前2種情況下的噴霧徑向擴散基本相同,隨后,在軸向和徑向上,反應噴霧情況下的混合物分數比無反應情況下(圖中為1.5 ms)擴散到更大的區域。這是由于反應情況下的徑向速度高于無反應情況下的徑向速度,這使得燃料蒸氣在徑向上的傳輸更快;反應情況下的,高溫熱氣的膨脹加速了噴霧的徑向擴散速度,這也使得反應情況下蒸氣穿透長度更長。

圖8 混合分數分布云圖

圖9顯示了反應情況下不同時刻氣體溫度云圖,黃色實線表示當量混合分數等值線,從中可以看到在反應前即0.4 ms時刻,氣體最高溫度約為1 200 K,初始燃燒室內環境溫度為900 K,這表明在燃燒前燃燒室內已經進行了低溫放熱反應;在0.5 ms時刻,已經發生燃燒,可以看到在距離燃油噴嘴下游約為20 mm處形成一個高溫火核,該火核以當量混合分數等值線為中心溫度逐漸向四周擴散降低,火核隨后長大,火焰沿著化學當量混合分數線傳播,在下游化學當量混合分數周圍形成高溫區,表明混合控制擴散火焰建立。

圖9 不同時刻溫度云圖

另外,發現擴散火焰并沒有附著在燃油噴嘴附近,而是在化學當量混合某處下游被點燃,表現為火焰浮升一定高度,火焰浮升高度隨時間有輕微波動,不同時刻火焰浮升位置在圖中用白色實線表示出來;火焰浮升高度是指燃油噴嘴與最靠近噴嘴的火焰前緣之間距離,試驗中通常使用激發態OH基達到最大值50%的位置定義火焰浮升高度,采用化學發光法測得,本文采用溫度達到1 600 K作為判斷標準。從圖9看到,火焰建立是在約距離噴嘴下游20 mm處,隨后,火焰穩定在距離噴嘴下游17 mm處,將不同時刻火焰浮升高度進行平均,得到平均火焰浮升高度為17 mm,這與SNL試驗得到16.4 mm幾乎一致,得到很好預測。這表明該噴霧燃燒模型很好地模擬了燃油點火過程中的瞬態行為和火焰浮升現象。

為了進一步說明噴霧燃燒點火過程,圖10給出了幾個關鍵組分質量分數分布云圖,其中黃色實線表示化學當量混合分數等值線,上文通過不同時刻溫度云圖得知,在霧化蒸發燃油蒸氣未被點燃前,燃燒室已經發生了反應,氣體溫度升高;確實如此,有反應情況的噴霧點火經歷2個階段,第一階段為低溫點火階段,如10a)圖和圖10b)所示,從0.3 ms到0.4 ms這一時間段內,酮過氧化氫(ketohydroperox-ide,KET,在這指OC12H23OOH)質量分數迅速降低,而H2O2含量迅速升高表示第一階段低溫點火開始[15];從0.4 ms到0.5 ms這一時間段內,如圖13中b)和c)所示,H2O2質量分數降低,OH含量上升表示第二階段高溫點火已經開始;在1.5 ms時刻,從OH分布可知,火焰沿著化學當量混合區域傳播,另外,值得注意的是,四副圖的右半為H2O2質量分數分布云圖,即使在高溫點火后,在靠近燃油噴嘴的化學當量混合分數區域內,H2O2仍然存在且含量變化不大,說明在高溫點火后,低溫點火會持續穩定進行。

圖10 不同組分質量分數分布云圖

圖11顯示了在有反應和無反應2種情況下流場中溫度隨局部混合分數分布散點圖。對于無反應噴霧情況,混合分數空間內的溫度分為兩部分,第一部分沿著紅色虛線密集分布,第二部分在紅色虛線下方,與第一部分相交于流場中混合分數的最大值處,此處的溫度約為585 K,此溫度為當前環境壓力下(約6 MPa)燃油的沸騰溫度;第一部分中的燃油蒸氣來自于燃油噴嘴出口附近油滴在沸騰溫度下的蒸發,由于蒸發的吸熱作用,當燃油與高溫空氣混合時,局部溫度會隨著混合分數的增加而降低。第二部分中的燃油蒸氣來自于溫度低于沸點的油滴,這些油滴的蒸發速率隨著溫度的升高而增加,因此在第二部分的燃油蒸氣局部混合分數也隨著溫度升高而變大。

圖11 流場溫度隨混合分數分布散點圖

在無反應噴霧情況下,不同時刻流場溫度分布變化不大,對于無反應噴霧情況,Mastorakos等[15]使用一個二次函數來擬合流場溫度與混合分數之間的關系

T=1 500Z2-1 500Z+Tam

(7)

式中:Z為混合分數,Tam為環境初始溫度,此算例中為900 K。

圖11中所有為紅色虛線為溫度與混合分數的二次函數解析擬合,該溫度與混合分數的關系與氣態燃料中不同,一般情況下,如果氣態燃料路易斯數(Lewis number)與該燃油蒸氣的路易斯數相同時,氣態燃料射流溫度與混合分數為線性關系,這種差異是由于燃油噴霧蒸發吸收了大量熱量,造成熱量損失造成的。

對有反應噴霧情況,在反應初期(0.2 ms),2種情況的流場溫度分布差別很大,在有反應情況下此刻進行低溫反應,將有反應情況的散點圖與無反應情況的擬合函數曲線進行比較,發現有反應情況在低溫反應開始前混合物的溫度分布與無反應情況下的相似;在0.4 ms時,混合物分數高于0.05的大多數混合物溫度在1 000～1 500 K之間,這說明低溫反應在此已經充分發展;1.2 ms 時的散點圖表明,此時流場可分為兩部分,從火焰浮升高度位置分開,火焰浮升高度上方(溫度高于1 200 K)流場中的燃燒呈現典型的擴散火焰結構,其溫度最大值對應的化學計量混合分數為0.045,火焰浮升高下方流場中的燃燒顯示出充分發展的冷焰(cool flame)結構。

2.4 注入壓力對點火延遲時間和火焰浮升高度的影響

圖12顯示了在不同燃油注入壓力下,點火延遲時間和火焰浮升高度的變化曲線,研究了注入壓力分別為50,100,150 MPa的3種情況。除了燃油注入壓力不同外,其他試驗條件均相同,環境溫度900 K,環境氣體組分與表2有反應情況相同?？梢钥吹?模擬較好地預測了點火延遲時間和火焰浮升高度趨勢,但在3種情況下,模擬的點火延遲時間預測均稍高于試驗數據,這可能與計算網格有關,本文所使用的網格預測到的混合分數要稍低于試驗數據,混合分數對點火延遲時間和火焰浮升高度的精準預測非常重要,這可能是導致在不同試驗條件下,模擬點火延遲時間均高于試驗數據的原因。

圖12 注入壓力對點火時間和火焰浮升高度的影響

點火延遲時間和火焰浮升高度這兩個參數對注入壓力的變化并不敏感,總的趨勢是隨著注入壓力的升高,點火延遲時間降低,火焰浮升高度有所上升;但我們看到,不同的燃油注入壓力對應的點火延遲時間在0.5 ms附近,相鄰2個不同的注入壓力之間點火延遲時間相差小于0.05 ms, 這也導致火焰浮升高度變化不大;有趣的是,注入壓力升高,點火延遲時間降低,火焰浮升高度卻有所升高,通常來講,點火延遲時間降低,火焰浮升高度也會降低,但在這正好相反,這是由于注入壓力升高,燃油霧化蒸發速率變大,第一階段的低溫反應能夠更快開始,點火延遲時間也會降低,但注入壓力的升高會給燃油液滴更高的初始動量,這就使得低溫反應和高溫點火在距離燃油噴嘴更遠的地方發生,即表現為火焰浮升高度的升高。

3 結 論

本文采用URANS模型和有限速率PaSR模型對Spray A無反應和有反應2種情況進行數值計算,主要結論如下:

1) 噴霧性能對計算網格依賴性很高,對網格加密能夠提升液體穿透長度和蒸氣穿透長度,但網格數量達到一定值后,再進行加密網格不會對結果準確性有很大提升;發現最小網格尺寸在0.25 mm時獲得的結果與試驗有較好的一致性。

2) 無反應情況下,由于燃油蒸發造成熱量損失,流場中溫度分布與混合分數成二次函數關系。

3) 反應情況下,Spray A噴霧燃燒分為2個階段,低溫反應首先開始釋放熱量,使得化學當量混合下游區域溫度升高,在這階段壓力不會升高,溫度到達一定值后,高溫點火開始,燃燒釋放大量熱量,燃燒室壓力線性升高,且在高溫點火后低溫反應仍持續穩定進行,起到穩定火焰的作用。

4) 無反應和有反應2種情況相比,在噴霧早期,2種情況下的蒸氣穿透長度近似,在噴霧后期,反應情況下,高溫熱氣的膨脹加速了噴霧的徑向擴散速度,使得反應情況下蒸氣穿透長度更長,蒸氣徑向分布范圍更大。

5) 反應情況下,提高燃油注入壓力,點火延遲時間減小,火焰浮升高度增加。