低階煤低溫干餾爐中低溫燃燒室溫度場(chǎng)分布模擬及結(jié)構(gòu)優(yōu)化

劉 俊，張永發(fā)，王 影，陳 磊，徐 英，趙海濱

（太原理工大學(xué)煤科學(xué)與技術(shù)教育部和山西省重點(diǎn)實(shí)驗(yàn)室，山西 太原 030024）

褐煤水分高、熱值低、易風(fēng)化和自燃，單位能量運(yùn)輸成本高，不利于長(zhǎng)距離輸送和貯存。褐煤熱解提質(zhì)可獲得寶貴的熱解煤氣、焦油和潔凈半焦[1-4]。因此國(guó)內(nèi)外研究者對(duì)褐煤低溫?zé)峤饧夹g(shù)進(jìn)行了大量的研究開發(fā)[5-9]。作者課題組開發(fā)了一種高效采油低溫干餾爐[10]，其結(jié)構(gòu)與現(xiàn)代焦?fàn)t結(jié)構(gòu)相似，高效采油低溫炭化爐中的傳熱方式為內(nèi)外混熱式，熱解溫度范圍為600～700℃。為適應(yīng)600～700℃熱解的要求，炭化爐中的燃燒室采用自主研發(fā)的旋流式低溫燃燒室。該低溫燃燒室的燃燒器采用擴(kuò)散鼓風(fēng)式燃燒器[11]，燃燒所需空氣由鼓風(fēng)機(jī)供給，燃燒過(guò)程屬于擴(kuò)散燃燒，根據(jù)對(duì)火焰形狀以及燃燒強(qiáng)度的要求采用套管式燃燒器；燃燒室以焦?fàn)t煤氣為燃料采用旋流式燃燒室，燃?xì)饧叭紵a(chǎn)物形成旋流使其內(nèi)部溫度分布均勻。

計(jì)算流體力學(xué)CFD是通過(guò)計(jì)算機(jī)模擬流體流動(dòng)、傳熱和燃燒等物理現(xiàn)象的技術(shù)，在研究燃燒器內(nèi)部流動(dòng)、多組分?jǐn)U散、化學(xué)反應(yīng)等復(fù)雜現(xiàn)象以及燃燒器優(yōu)化設(shè)計(jì)中起著重要的作用。人們已對(duì)大型工業(yè)燃燒器進(jìn)行了三維全尺寸數(shù)值模擬與優(yōu)化調(diào)整，優(yōu)化后的燃燒器運(yùn)行工況良好，解決了原燃燒經(jīng)常出現(xiàn)的問(wèn)題[12-14]。文獻(xiàn)檢索表明，目前未見(jiàn)低溫燃燒室（低于850℃）結(jié)構(gòu)及其相關(guān)研究。本文采用Fluent 6.3 CFD軟件，對(duì)自主研發(fā)的褐煤熱解高效采油低溫炭化爐中旋流式低溫燃燒室內(nèi)部溫度場(chǎng)分布進(jìn)行模擬計(jì)算，并進(jìn)行了燃燒室結(jié)構(gòu)優(yōu)化。

1 實(shí)驗(yàn)部分

1.1 燃料煤氣的組成

采用焦?fàn)t煤氣作為燃料，其組成如表1所示。

表1 焦?fàn)t煤氣主要組分的體積分?jǐn)?shù)

1.2 燃燒室

燃燒室結(jié)構(gòu)為長(zhǎng)方體結(jié)構(gòu)（1500 mm×200 mm×1500 mm），煙氣出口（150 mm×150 mm）位于燃燒室中央，燃燒器采用套管式噴頭，噴頭分別位于燃燒室右側(cè)頂部和左側(cè)底部。燃燒室具體結(jié)構(gòu)如圖1。對(duì)其進(jìn)行網(wǎng)格劃分，目前常采用區(qū)域法或分區(qū)網(wǎng)格[15]，基本思想是，根據(jù)外形特點(diǎn)把復(fù)雜的物理域或復(fù)雜的拓?fù)浣Y(jié)構(gòu)劃分成簡(jiǎn)單的網(wǎng)格（結(jié)構(gòu)塊網(wǎng)格）。燃燒器焦?fàn)t煤氣和空氣進(jìn)口處采用T-grid方法進(jìn)行網(wǎng)格劃分，形成四面體網(wǎng)格，長(zhǎng)方體部分即爐膛采用Hex/Wedge copper方法[16]進(jìn)行網(wǎng)格劃分，形成六面體網(wǎng)格。網(wǎng)格劃分情況如圖1所示。

1.3 工況和邊界條件

噴嘴中心圓通道（內(nèi)管）走焦?fàn)t煤氣，圓環(huán)通道（外管）走空氣，通常焦?fàn)t燃燒室空氣過(guò)量系數(shù)[17]控制在1.1～1.25范圍內(nèi)，本實(shí)驗(yàn)空氣過(guò)量系數(shù)為定值1.25。燃?xì)膺M(jìn)口、空氣進(jìn)口、煙氣出口和壁面的邊界條件見(jiàn)表2。

2 結(jié)果與討論

2.1 數(shù)學(xué)模型的建立

2.1.1 湍流模型

本模擬對(duì)象選取帶旋流修正的湍流模型（Realizablek-ε模型[18]）作為湍流模型?？蓪?shí)現(xiàn)k-ε模型的湍動(dòng)能及其耗散率輸運(yùn)方程如式（1）、式（2）。

表2 模擬工況的邊界條件及參數(shù)

式中，Gk表示由于平均速度梯度引起的湍動(dòng)能產(chǎn)生；Gb表示由于浮力影響引起的湍動(dòng)能產(chǎn)生；YM表示可壓縮湍流脈動(dòng)膨脹對(duì)總的耗散率的影響；C2和C1ε是常數(shù)；σk和σε分別是湍動(dòng)及其耗散率的湍流普朗特?cái)?shù)。

在Fluent中，作為默認(rèn)值常數(shù)，C1ε=1.44，C2=1.9，σk=1.0，σε=1.2。該模型適合的流動(dòng)類型比較廣泛，包括有旋均勻剪切流、自由流（射流和混合層）、腔道流動(dòng)和邊界層流動(dòng)，特別是可實(shí)現(xiàn)k-ε模型對(duì)圓口射流和平板射流模擬中，能給出較好的射流擴(kuò)張角，所以本模擬選取帶旋流修正的湍流模型。

2.1.2 輻射模型

焦?fàn)t煤氣和空氣的非預(yù)混燃燒最高溫度為2400 K，因此輻射傳熱不可忽視，在計(jì)算時(shí)必須考慮輻射傳熱[16]。本模擬采用P-1輻射模型，P1模型僅適用于光學(xué)厚度αL（在計(jì)算輻射傳輸時(shí)，兩個(gè)給定高度層之間的單位截面鉛直氣柱內(nèi)特定的吸收或發(fā)射物質(zhì)的質(zhì)量）較大的場(chǎng)合，即僅適用于計(jì)算域較大的場(chǎng)合。對(duì)于該燃燒室內(nèi)的流動(dòng)，L為燃燒室的直徑，αL＞1。P-1輻射模型[19]是P-N模型中最簡(jiǎn)單的類型，P-N模型的出發(fā)點(diǎn)是把輻射強(qiáng)度展開成為正交的球諧函數(shù)（正交序列級(jí)數(shù)）。如果只取正交球諧函數(shù)的前4項(xiàng)，對(duì)于輻射熱流qr，能得到式（3）。

式中，a為吸收系數(shù)，σs為散射系數(shù)，G為入射輻射，C為線性各相異性相位函數(shù)系數(shù)。

2.1.3 化學(xué)組分輸送和反應(yīng)流的模型

模擬燃燒的模型有通用有限速率模型、非預(yù)混模型、預(yù)混模型和部分預(yù)混模型4種。焦?fàn)t煤氣和空氣以相異流進(jìn)入燃燒室，因此本模擬采用非預(yù)混模型[20]。

在非預(yù)混燃燒模型里，并不是解每一個(gè)組分輸送方程，而是解一個(gè)或兩個(gè)守恒標(biāo)量（混合分?jǐn)?shù)）的輸送方程，然后從預(yù)測(cè)的混合分?jǐn)?shù)分布推導(dǎo)出每一個(gè)組分的濃度。在一定假設(shè)條件下，流體的瞬時(shí)熱化學(xué)狀態(tài)與一個(gè)守恒量有關(guān)，這個(gè)守恒量就是混合分?jǐn)?shù)，用ω表示，見(jiàn)式（4），是來(lái)自燃料流的元素質(zhì)量分?jǐn)?shù)。

式中，ωi為元素i的質(zhì)量分?jǐn)?shù)；下角標(biāo)o為氧化劑流入口處的值；f為燃料流入口處的值?；旌戏?jǐn)?shù)是在所有組分（CO2，H2O，O2等）里，燃燒和未燃燒燃料流元素（C，H等）的局部質(zhì)量分?jǐn)?shù)。該方法主要用于模擬湍流擴(kuò)散火焰。在守恒標(biāo)量方法中，通過(guò)概率密度函數(shù)PDF[21]來(lái)考慮湍流的影響。層流火焰面模型是非預(yù)混模型的擴(kuò)展，它考慮了從化學(xué)平衡狀態(tài)形成的空氣動(dòng)力學(xué)的應(yīng)力誘導(dǎo)分離。該模擬非預(yù)混燃燒的概率密度函數(shù)PDF文件中燃燒平均溫度如圖2。

平均混合分?jǐn)?shù)為零時(shí)代表空氣流，平均混合分?jǐn)?shù)為1時(shí)代表燃料流。從圖2可以看出，在平均混合分?jǐn)?shù)在0.1左右時(shí)平均溫度為最高溫度2400 K。

2.2 燃燒模擬溫度分布及燃燒室結(jié)構(gòu)優(yōu)化

2.2.1 無(wú)擋板燃燒室的溫度分布模擬與實(shí)驗(yàn)驗(yàn)證

在圖1所示無(wú)擋板非預(yù)混燃燒室結(jié)構(gòu)中，23 m3/h的焦?fàn)t煤氣與126 m3/h空氣在燃燒室內(nèi)發(fā)生燃燒反應(yīng)。煙氣出口所在平面的溫度分布圖如圖3。

從圖3可以觀察到，高溫區(qū)主要集中于燃燒室煙氣出口平面外圍，燃燒室頂部和底部溫度偏高，燃燒室中央?yún)^(qū)域溫度達(dá)不到650～750℃，整體溫度分布不均勻，經(jīng)過(guò)傳熱后將導(dǎo)致炭化室中央?yún)^(qū)域的半焦不熟，達(dá)不到生產(chǎn)褐煤半焦的工業(yè)要求。為驗(yàn)證模擬結(jié)果的準(zhǔn)確性，所示結(jié)構(gòu)的燃燒室內(nèi)，同樣條件下進(jìn)行燃燒實(shí)驗(yàn)。燃燒工況和邊界條件如表3，溫度測(cè)量采用B型熱電偶。

表3 燃燒工況和邊界條件

在工況一條件下，實(shí)際測(cè)量值和模擬計(jì)算值的對(duì)比如圖4所示。從圖4可以看出，模擬計(jì)算溫度值比實(shí)際測(cè)量值高，誤差波動(dòng)幅度為50～70℃。其原因在于實(shí)驗(yàn)燃燒室保溫隔材料的熱耗散，所以模擬計(jì)算結(jié)果可以與實(shí)際測(cè)量值較好的吻合，能夠準(zhǔn)確反映燃燒室內(nèi)的溫度分布情況。實(shí)驗(yàn)結(jié)果與模擬計(jì)算結(jié)果能夠吻合，但燃燒室在不加擋板時(shí)，其中央?yún)^(qū)域溫度不能滿足熱解需要溫度，因此需要進(jìn)一步調(diào)整結(jié)構(gòu)使燃燒室頂部和底部的高溫區(qū)移動(dòng)至中間區(qū)域，使燃燒室整體溫度更加均勻。

2.2.2 加擋板燃燒室最優(yōu)結(jié)構(gòu)選取與實(shí)驗(yàn)驗(yàn)證

為了使燃燒室整體溫度均勻，在其內(nèi)部增加兩塊彎擋板，保持焦?fàn)t煤氣和空氣流量分別為23 m3/h和126 m3/h下進(jìn)行模擬燃燒，流體流動(dòng)情況改變。

（1）擋板豎直位置的確定 保持兩塊擋板的水平位置為x=375 mm和圓心角為100°的情況下，調(diào)整其豎直位置，調(diào)整后的燃燒溫度分布情況如圖5。

由圖5可以看出，在（?375，450，0）和（375，?450，0）處安裝兩塊半徑、厚度和圓心角依次為100 mm、10 mm、100°的擋板，觀察在該結(jié)構(gòu)下的燃燒溫度分布，較不加擋板時(shí)，該燃燒室中央?yún)^(qū)域溫度增加了100℃，燃燒室豎拋平面外圍溫度高，但是燃燒室中央?yún)^(qū)域溫度仍然達(dá)不到650～750℃，需繼續(xù)增加擋板的豎直高度。

由圖6可見(jiàn)，在（?375，500，0）和（375，?500，0）處安裝兩塊半徑、厚度和圓心角依次為100 mm、10 mm、100°的擋板，觀察在該結(jié)構(gòu)下的燃燒溫度分布，燃燒室中央溫度較不加擋板時(shí)增加了100℃，高溫區(qū)略微向燃燒室中央移動(dòng)，但燃燒室中央?yún)^(qū)域溫度仍達(dá)不到600～650℃，需繼續(xù)增加擋板的豎直高度。

由圖7可見(jiàn)，在（?375，550，0）和（375，?550，0）處安裝兩塊半徑、厚度和圓心角依次為100 mm、10 mm、100°的擋板，觀察在該結(jié)構(gòu)下的燃燒溫度分布，燃燒室中央?yún)^(qū)域溫度均為526.85℃，雖然燃燒室中央?yún)^(qū)域溫度達(dá)不到650～750℃，但該結(jié)構(gòu)下的燃燒室中央低溫區(qū)較圖6縮小，所以選取y=550 mm處設(shè)置擋板。

銀行理財(cái)業(yè)務(wù)受到投資者的青睞，產(chǎn)品銷售數(shù)量突飛猛進(jìn)，不同銀行及業(yè)務(wù)的產(chǎn)品性質(zhì)相同，理財(cái)產(chǎn)品千篇一律，各銀行無(wú)特色及特長(zhǎng)，競(jìng)爭(zhēng)力較為薄弱，同時(shí)也造成理財(cái)市場(chǎng)產(chǎn)品過(guò)于單一，銀行理財(cái)產(chǎn)品同質(zhì)化嚴(yán)重，不能對(duì)多層次客戶提供需求。

（2）擋板的水平位置的確定 保持兩塊擋板的豎直位置y=550 mm和圓心角角度為100°的情況下，調(diào)整其水平位置，調(diào)整后的燃燒溫度分布情況如圖8、圖9。

由圖8可見(jiàn)，在（?175，550，0）和（175，?550，0）處安裝兩塊半徑、厚度和圓心角依次為100 mm、10 mm、100°的擋板，觀察在該結(jié)構(gòu)下的燃燒溫度分布，燃燒室中央?yún)^(qū)域溫度為526.85℃，該結(jié)構(gòu)下的燃燒室中央低溫區(qū)較圖7無(wú)顯著變化，需要繼續(xù)調(diào)整擋板的水平位置。

由圖9可見(jiàn)，在（?25，550，0）和（25，?550，0）處安裝兩塊半徑、厚度和圓心角依次為100 mm、10 mm、100°的擋板，觀察在該結(jié)構(gòu)下的燃燒溫度分布，燃燒室中央?yún)^(qū)域溫度較圖8升高100℃，能夠達(dá)到熱解要求溫度650～750℃，但煙氣出口溫度為526.85℃，達(dá)不到內(nèi)熱式傳熱要求溫度。下一步調(diào)整擋板圓心角角度進(jìn)行結(jié)構(gòu)優(yōu)化。

（3）擋板圓心角角度的確定 保持兩塊擋板的水平位置x=±25 mm和豎直位置y=±550 mm的情況下，調(diào)整擋板的圓心角角度，調(diào)整后的燃燒溫度分布情況如圖8。

由圖10可見(jiàn)，在（?25，550，0）和（25，?550，0）處安裝兩塊半徑、厚度和圓心角依次為100 mm、10 mm、60°的擋板，觀察在該結(jié)構(gòu)下的燃燒溫度分布，燃燒室中央?yún)^(qū)域溫度為626.85℃， 該結(jié)構(gòu)下的燃燒室中央低溫區(qū)較無(wú)擋板時(shí)溫度降低，既不能滿足熱解要求溫度，又達(dá)不到內(nèi)熱式傳熱要求溫度。因此需繼續(xù)調(diào)整彎擋板角度。

由圖11可見(jiàn)，在（?25，550，0）和（25，?550，0）處安裝兩塊半徑、厚度和圓心角依次為100 mm、10 mm、80°的擋板，觀察在該結(jié)構(gòu)下的燃燒溫度分布，燃燒室中央?yún)^(qū)域溫度為526.85℃，該結(jié)構(gòu)下的燃燒室溫度分布既不能滿足熱解要求溫度，又達(dá)不到內(nèi)熱式傳熱要求溫度。因此需繼續(xù)調(diào)整彎擋板角度。

由圖13可見(jiàn)，在（?25，550，0）和（25，?550，0）處安裝兩塊半徑、厚度和圓心角依次為100 mm、10 mm、130°的擋板，觀察在該結(jié)構(gòu)下的燃燒溫度分布，燃燒室中央?yún)^(qū)域溫度為626.85～726.85℃，煙氣出口溫度為626.85℃，該結(jié)構(gòu)下的燃燒室中溫度分布650～750℃既能滿足熱解要求溫度500～600℃，又能達(dá)到到內(nèi)熱式傳熱要求溫度600～650℃，燃燒室整體溫度分布較為均勻，因此在該結(jié)構(gòu)下，燃燒室的溫度分布能夠達(dá)到低溫炭化爐配套燃燒室的傳熱要求溫度。

（4）燃燒室最優(yōu)結(jié)構(gòu)的確定與實(shí)驗(yàn)驗(yàn)證 通過(guò)對(duì)比不同水平、豎直位置以及擋板角度下燃燒室的溫度分布圖，最終確定圖13所示燃燒室結(jié)構(gòu)為燃燒室最優(yōu)結(jié)構(gòu)，在該結(jié)構(gòu)下燃燒室中溫度分布為650～750℃既能滿足熱解要求溫度500～600℃，又能達(dá)到到內(nèi)熱式傳熱要求溫度600～650℃，燃燒室整體溫度分布較為均勻，其不同豎直高度及水平位置（以燃燒室左下角為水平位置零點(diǎn)）下的平均溫度分布如圖14。燃燒室的最優(yōu)結(jié)構(gòu)如圖15。

為驗(yàn)證模擬結(jié)果的準(zhǔn)確性，同樣條件下進(jìn)行實(shí)驗(yàn)，23 m3/h焦?fàn)t煤氣和126 m3/h空氣在最優(yōu)結(jié)構(gòu)下的燃燒室內(nèi)進(jìn)行燃燒。燃燒工況和邊界條件如表4，溫度測(cè)量采用B型熱電偶。

在工況二條件下，實(shí)際測(cè)量值和模擬計(jì)算值的對(duì)比如圖16所示。從圖可以看出，模擬計(jì)算溫度值比實(shí)際測(cè)量值高，誤差波動(dòng)幅度為50～70℃。其原因在于實(shí)驗(yàn)燃燒室保溫隔材料的熱耗散，所以模擬計(jì)算結(jié)果可以與實(shí)際測(cè)量值較好的吻合，能夠準(zhǔn)確反映燃燒室內(nèi)的溫度分布情況。

表4 燃燒工況和邊界條件

3 結(jié)論

（1）旋流式燃燒室內(nèi)無(wú)擋板時(shí)，燃燒室內(nèi)溫度分布不均勻，中央?yún)^(qū)域溫度偏低，外圍溫度偏高。在燃燒室內(nèi)添加不同形狀的擋板，燃燒室外圍的高溫區(qū)熱量移動(dòng)至中央?yún)^(qū)域，在燃燒室內(nèi)坐標(biāo)（?25，550，0）和（25，?550，0）處增加2個(gè)關(guān)于燃燒室中央軸對(duì)稱且半徑為100 mm的130°圓弧形擋板，燃燒室內(nèi)部溫度分布比較均勻且能達(dá)到褐煤熱解高效采油低溫干餾爐的熱解溫度且溫度分布均勻。褐煤熱解旋流式低溫燃燒室平均溫度可控制在750～850℃，燃燒室結(jié)構(gòu)滿足褐煤低溫?zé)峤馓炕癄t操作溫度的要求。

（2）Realizablek-ε湍流模型、P-1輻射模型和非預(yù)混燃燒模型能夠準(zhǔn)確的用于模擬褐煤熱解高效采油低溫干餾爐中旋流式低溫燃燒室的燃燒狀況。在23 m3/h焦?fàn)t煤氣126 m3/h空氣工業(yè)條件下燃燒實(shí)驗(yàn)，結(jié)果表明模擬計(jì)算溫度值和實(shí)際測(cè)量值基本一致，誤差波動(dòng)幅度為50～70℃，滿足工業(yè)要求。

[1]尹立群.我國(guó)褐煤資源及其利用前景[J].煤炭科學(xué)技術(shù)，2004，32（8）：12-14.

[2]石文秀，張玉財(cái)，金管會(huì).淺談褐煤研究的必要性及褐煤的性質(zhì)[J].化工進(jìn)展，2012，31（s1）：88-94.

[3]Atwood M T，Schulman B L.Toscoal process emdash pyrolysis of western coals and lignites for char and oil production[J].Preprints ofPapers American Chemical Society Division of Fuel Chemical，1997，22（2）：233-252.

[4]李香蘭，崔新濤，張永發(fā).GC-MS在內(nèi)蒙褐煤型煤塊低溫煤焦油成分分析中的應(yīng)用[J].分析儀器，2012（3）：17-24.

[5]游偉，趙濤，章衛(wèi)星，等.美國(guó)低階煤提質(zhì)技術(shù)發(fā)展概述[J].化肥設(shè)計(jì)，2009，47（7）：5-9.

[6]張慧榮，張永發(fā).美國(guó)LFC和ACCP褐煤提質(zhì)工藝概述[J].能源與節(jié)能，2011，2（65）：71-74.

[7]Sass A.Garrett’s coal pyrolysis process[J].Chemical Engineering Progress，1974，70（1）：72-73.

[8]Carlson F B，Yardumian L H，Atwood M T.Toscoal process for low temperature pyrolysis of coal[J].Transaction s of the Society of Mining Engineers of AIME，1974， 256（2）：128-131.

[9]Atwood M T，Schulman B L.Toscoal process emdash pyrolysis of western coals and lignites for char and oil production[J].Preprints of Papers American Chemical Society Division of Fuel Chemical，1997，22（2）：233-252.

[10]張永發(fā)，郭紅生，陳磊，等.導(dǎo)氣磚及其導(dǎo)氣墻及應(yīng)用該磚墻的粉煤干餾炭化室：中國(guó)，201210119488.0[P].2012-4-23.

[11]劉蓉，劉文斌.燃?xì)馊紵c燃燒裝置[M].北京：機(jī)械工業(yè)出版社，2009.

[12]王家楣，彭峰.燃燒器三維流動(dòng)和燃燒的數(shù)值模擬及優(yōu)化結(jié)果[J].武漢理工大學(xué)學(xué)報(bào)，2004，26（3）：79-82.

[13]陸曉峰，李元青.基于CFD計(jì)算的燃燒器結(jié)構(gòu)改進(jìn)研究[J].石油學(xué)報(bào)，2011，10，27（5）：787-791.

[14]郭英峰，李保有，劉海軍.賽科裂解爐輻射段燃燒的數(shù)值模擬及優(yōu)化[J].化工進(jìn)展，2009，28（6）：1090-1094.

[15]黃志超，包忠詡，周天瑞，等.有限元網(wǎng)格劃分技術(shù)[J].南昌大學(xué)學(xué)報(bào)，2001，23（4）：25-31.

[16]Fluent Inc.Fluent users guide[M].USA ：Fluent Inc.，2003.

[17]《焦化設(shè)計(jì)參考資料》編寫組.焦化設(shè)計(jì)參考資料[M].北京：冶金工業(yè)出版社，1978.

[18]王福軍.計(jì)算流體動(dòng)力學(xué)分析[M].北京：清華大學(xué)出版社，2004.

[19]Bidi M，Hosseini R.Numerical analysis of methane–air combustion considering radiation effect[J].Energy Conversion and Management，2008，49：3634-3647.

[20]Carbonell D，Perez-Segarra C D.Flamelet mathematical models for non-premixed laminar combustion[J].Combustion and Flame，2009，156：334-347.

[21]Ziani Lotfi，Chaker Abla，Chetehouna Khaled，et al.Numerical simulationsofnon-premixed turbulentcombustionofCH4-H2mixtures using the PDF approach[J].International Journal of Hydrogen Energy，2013，38（20）：8597-8603.