隨機(jī)場(chǎng)TPDF在非結(jié)構(gòu)網(wǎng)格平臺(tái)的開發(fā)及檢驗(yàn)

王 方 ，楊興峰 ，方 存 ，許笑顏 ，莫 毅 ，高浩卜

（1.北京航空航天大學(xué)能源與動(dòng)力工程學(xué)院，北京 100191；2.北京航空航天大學(xué)江西研究院，南昌 330096；3.北京航空航天大學(xué)成都航空動(dòng)力創(chuàng)新研究院，成都 611930；4.中國航空發(fā)動(dòng)機(jī)研究院，北京 101300）

0 引言

航空發(fā)動(dòng)機(jī)燃燒室?guī)缀谓Y(jié)構(gòu)非常復(fù)雜，先進(jìn)航空發(fā)動(dòng)機(jī)設(shè)計(jì)對(duì)燃燒室提出更精確的燃燒組織要求，對(duì)燃燒室數(shù)值模擬提出了更高的精度要求[1]，對(duì)湍流燃燒模型也提出了更高的要求[2]。概率密度函數(shù)輸運(yùn)方程湍流燃燒模型（Transport Probability Density Function Turbulence Combustion Model，TPDF）可以耦合詳細(xì)化學(xué)反應(yīng)機(jī)理，高精度解析湍流燃燒過程[3]。

概率密度函數(shù)輸運(yùn)方程湍流燃燒模型[4]主要分為3種：速度-標(biāo)量聯(lián)合概率密度函數(shù)輸運(yùn)方程；標(biāo)量聯(lián)合概率密度函數(shù)輸運(yùn)方程；速度、標(biāo)量和耗散率三者聯(lián)合的概率密度函數(shù)輸運(yùn)方程。目前，3 種TPDF中應(yīng)用最廣泛的是標(biāo)量聯(lián)合的概率密度函數(shù)輸運(yùn)方程。在標(biāo)量聯(lián)合概率密度函數(shù)輸運(yùn)方程解法中，歐拉隨機(jī)場(chǎng)法是一種概率密度函數(shù)的求解方法[5]。Jones等[5-7]采用射流、鈍體、旋流、燃燒室等算例檢驗(yàn)了歐拉隨機(jī)場(chǎng)解法的TPDF 應(yīng)用；王方等[8-9]在BOFFIN 程序的基礎(chǔ)上，發(fā)展了適用于航空發(fā)動(dòng)機(jī)兩相湍流燃燒的AECSC 軟件，并針對(duì)化學(xué)反應(yīng)機(jī)理、兩相蒸發(fā)模型等進(jìn)行了研究；Gong 等[10]和Rihab 等[11]合作，在開源軟件OpenFOAM 平臺(tái)上開發(fā)了TPDF 燃燒模型代碼，求解聲速問題。TPDF可望成為燃燒室高精度模擬的基礎(chǔ)模型之一[12]。

射流火焰和旋流火焰是燃燒室湍流燃燒情況的典型代表。Flame D 是非預(yù)混射流火焰，測(cè)量數(shù)據(jù)多且可靠，是目前湍流火焰預(yù)測(cè)方法檢驗(yàn)的必經(jīng)之路。比如Nik 等[13]檢驗(yàn)速度標(biāo)量聯(lián)合濾波密度函數(shù)的方法；Renzo 等[14]驗(yàn)證FPV 湍流燃燒模型；Lysonko 等[15]檢驗(yàn)渦耗散模型；Ge 等[16]檢驗(yàn)稀疏拉格朗日MMC 模型。此外，還有學(xué)者用Flame D 研究輻射模型[17]影響、RCCE 應(yīng) 用[18]等；TECFLAM 是Dreizler 等[19]開 發(fā) 的 甲烷/空氣非受限強(qiáng)旋流貧油預(yù)混燃燒器；Freitag 等[20-22]利用直接數(shù)值模擬和大渦模擬（Large Eddy Simulation，LES）對(duì)TECFLAM 冷態(tài)和熱態(tài)工況進(jìn)行研究，并利用關(guān)聯(lián)函數(shù)的方法對(duì)冷態(tài)流動(dòng)工況中渦團(tuán)尺度進(jìn)行了分析；Wegner等[23]利用非穩(wěn)態(tài)雷諾應(yīng)力模擬（Unsteady Reynolds Averaged Navier-Stokes，URANS）和LES 對(duì)TECFLAM 冷態(tài)工況進(jìn)行了模擬，分析了URANS 在模擬旋進(jìn)渦核結(jié)構(gòu)的適用性，并給出了URANS和LES 處理旋進(jìn)頻率的不同方法；隋春杰[24]用LES 和改進(jìn)代數(shù)2 階矩湍流燃燒模型，結(jié)合甲烷2 步化學(xué)反應(yīng)機(jī)理，成功模擬了TECFLAM 冷態(tài)和熱態(tài)工況，指出入口當(dāng)量比擾動(dòng)對(duì)預(yù)混燃燒模擬的影響；在Rolls-Royce 公司對(duì)PRECISE-UNS 軟件發(fā)展測(cè)試工作中，對(duì)非結(jié)構(gòu)網(wǎng)格下的火焰面類、TPDF 等模型進(jìn)行了射流、旋流和燃燒室測(cè)試研究[12]。

Saturne 軟件是開源軟件，能夠應(yīng)用非結(jié)構(gòu)網(wǎng)格。在Saturne 平臺(tái)構(gòu)建TPDF 并對(duì)射流、旋流、燃燒室進(jìn)行模擬測(cè)試，可以檢驗(yàn)TPDF 對(duì)航空發(fā)動(dòng)機(jī)燃燒室的適用性，判斷能否高精度求解復(fù)雜結(jié)構(gòu)燃燒室內(nèi)部燃燒過程。本文在非結(jié)構(gòu)網(wǎng)格體系中研發(fā)隨機(jī)場(chǎng)TPDF方法及配合的計(jì)算單元，測(cè)試模擬效果，以促進(jìn)航空發(fā)動(dòng)機(jī)燃燒室模擬技術(shù)難題的解決。

1 非結(jié)構(gòu)網(wǎng)格平臺(tái)、TPDF開發(fā)和算例介紹

1.1 非結(jié)構(gòu)網(wǎng)格平臺(tái)

Saturne 是由法國電力公司研究院（EDF R&D）設(shè)計(jì)開發(fā)的通用CFD 軟件，可求解3 維Navier-Stokes 方程。Saturne軟件基于有限體積離散，使用Fortran和C語言作為主要編程語言，支持非結(jié)構(gòu)化網(wǎng)格并行計(jì)算。其涵蓋包括湍流模型在內(nèi)的多個(gè)物理模塊，并且可以兼容目前大多數(shù)類型的網(wǎng)格文件，如cgns、gambit等格式，具有積極的工程應(yīng)用價(jià)值。

1.2 TPDF簡(jiǎn)介

詳細(xì)的歐拉隨機(jī)場(chǎng)的TPDF 介紹可以參考文獻(xiàn)[5]和[25]。考慮密度加權(quán)平均的標(biāo)量聯(lián)合概率密度函數(shù)輸運(yùn)方程為

式中：?P為Favre 濾波后的概率密度函數(shù)；ρˉ為空間濾波后的密度；?uj為Favre 濾波后的j方向速度；ψ為標(biāo)量；α為組分；F為聯(lián)合概率密度函數(shù)；ω?α為組分α的化學(xué)反應(yīng)源項(xiàng)；Ns為組分?jǐn)?shù)量；?α為與α有關(guān)的標(biāo)量；ψα為?α構(gòu)成的樣本空間；ψ為ψα構(gòu)成的點(diǎn)的集合；μ為動(dòng)力粘度；σ為施密特?cái)?shù)。

式中左側(cè)3 項(xiàng)分別為時(shí)間項(xiàng)、對(duì)流項(xiàng)和化學(xué)反應(yīng)項(xiàng)，均是可以精確計(jì)算的。右側(cè)2 項(xiàng)則需要通過梯度擴(kuò)散模型和小尺度混合模型進(jìn)行模化，得到最終封閉形式的TPDF方程為

式中：μT為湍流粘度；σT為湍流施密特?cái)?shù)；τ為湍流混合時(shí)間尺度；??α為Favre濾波后的?α。

設(shè)有N個(gè)隨機(jī)場(chǎng)的集合ξm(x,t)=[ξm1,…,ξmN]，該隨機(jī)場(chǎng)對(duì)于所有標(biāo)量都是歐拉形式的。采用隨機(jī)場(chǎng)解法的Ito形式求解TPDF方程，得到最終可以應(yīng)用到程序中的TPDF方程形式為

1.3 基于Saturne軟件的TPDF研發(fā)

TPDF 構(gòu)建框架如圖1 所示。從圖中可見，僅用紅框標(biāo)出的文件為修改的文件，主要目的是實(shí)現(xiàn)TPDF的獨(dú)立；全部標(biāo)紅的文件為新增的文件，這是湍流燃燒模型的主要求解文件，用于增加組分輸運(yùn)方程和能量方程以及數(shù)據(jù)的傳輸和其他功能構(gòu)建。本文在Saturne 開發(fā)的TPDF 構(gòu)建了質(zhì)量分?jǐn)?shù)輸運(yùn)場(chǎng)和總焓輸運(yùn)場(chǎng)，另外創(chuàng)建了一些不參與輸運(yùn)方程求解的包括物種生成率、溫度等在內(nèi)的屬性場(chǎng)。主要求解文件是pdfini、pdftss、pdfphy、pdftcl，分別用來實(shí)現(xiàn)模型初始化、源項(xiàng)求解、物理屬性更新、邊條處理。

圖1 TPDF構(gòu)建框架

1.4 算例介紹和邊界條件

1.4.1 Flame D算例

Flame D 的幾何結(jié)構(gòu)見文獻(xiàn)[23、26]。FUEL 區(qū)為主噴射口，直徑d=7.2 mm，噴射工質(zhì)為體積分?jǐn)?shù)25%的甲烷和75%的空氣組成的混氣。PILOT 區(qū)為高溫值班火焰，主要成分是已燃?xì)猓唧w的各物質(zhì)質(zhì)量分?jǐn)?shù)以及不同軸向位置截面的試驗(yàn)數(shù)據(jù)可以在Sandia實(shí)驗(yàn)室官方網(wǎng)站[27]查到。Flame D 算例邊界條件見表1。

表1 Flame D算例邊界條件

1.4.2 TECFLAM算例

TECFLAM 燃燒器幾何結(jié)構(gòu)如圖2 所示。甲烷和空氣先在燃燒器進(jìn)口位置充分摻混，然后經(jīng)旋流葉片進(jìn)入主要的燃燒區(qū)域。本文選定的熱態(tài)試驗(yàn)工況為30 kW的熱負(fù)荷下的試驗(yàn)工況，旋流數(shù)為0.75，甲烷空氣當(dāng)量比為0.83，預(yù)混氣總體積流量為37.92 m3/s[28]。冷態(tài)和熱態(tài)工況的旋流和伴流進(jìn)口邊界條件見表2、3，其余所有壁面均為無滑移絕熱壁面條件，出口為壓力出口條件。

圖2 TECFLAM 燃燒器幾何結(jié)構(gòu)[6]

表2 TECFLAM 冷態(tài)算例進(jìn)口邊界條件

表3 TECFLAM 熱態(tài)算例進(jìn)口邊界條件

1.5 計(jì)算網(wǎng)格和化學(xué)反應(yīng)機(jī)理

用ICEM 軟件繪制網(wǎng)格劃分O-BLOCK，F(xiàn)lame D網(wǎng)格軸向長(zhǎng)度0.54 m，徑向長(zhǎng)度0.3 m，總網(wǎng)格單元數(shù)目約為30 萬，對(duì)進(jìn)口和中心射流區(qū)進(jìn)行加密，F(xiàn)lame D 算例網(wǎng)格如圖3 所示。TECFLAM 算例采用文獻(xiàn)[24]的網(wǎng)格劃分策略，TECFLAM 算例網(wǎng)格如圖4 所示。文獻(xiàn)[22]采用長(zhǎng)度2D的環(huán)形通道使流動(dòng)發(fā)展，計(jì)算區(qū)域總軸向長(zhǎng)度為10D。D代表TECFLAM 中心鈍體的直徑30 mm。最終的網(wǎng)格單元總數(shù)目約為80萬。本文選用Zhou 等[29]的甲烷單步化學(xué)反應(yīng)機(jī)理，在mol-cm-g-s 單位制標(biāo)準(zhǔn)Arrhenius 表達(dá)式中，A=3.2×1014cm3/（mol·s），b=0.0，E=36560 kcal/mol，甲烷和氧氣的反應(yīng)級(jí)數(shù)均為1.0。

圖3 Flame D算例網(wǎng)格

圖4 TECFLAM 算例網(wǎng)格

2 模擬結(jié)果與討論

2.1 Flame D算例

RANS-TPDF 模擬的Flame D 時(shí)均溫度如圖5 所示。Flame D 火焰核心射流進(jìn)口流速較高，以噴嘴出口速度為特征速度的雷諾數(shù)超過20000，屬于典型高雷諾數(shù)部分預(yù)混射流火焰。其PILOT 區(qū)提供了1880 K 的高溫，甲烷被高溫點(diǎn)燃發(fā)生化學(xué)反應(yīng)釋放出大量熱，瞬時(shí)高溫可達(dá)2200～2300 K。隨著射流發(fā)展甲烷不斷被消耗，到中部截面處高溫區(qū)匯合，形成“蠟燭形”火焰。之后，高溫氣體繼續(xù)擴(kuò)散。

圖5 Flame D時(shí)均溫度

分別將模擬平均溫度和甲烷濃度與試驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，F(xiàn)lame D 不同軸向位置時(shí)均溫度徑向分布如圖6 所示，F(xiàn)lame D 不同軸向位置時(shí)均甲烷質(zhì)量分?jǐn)?shù)徑向分布如圖7 所示。從圖6、7 中可見，模擬結(jié)果總體與試驗(yàn)值有較好一致性。前3 個(gè)截面是燃燒發(fā)展區(qū)，流動(dòng)主導(dǎo)而反應(yīng)剛開始，吻合較好；中間3 個(gè)截面化學(xué)反應(yīng)起主要作用，放熱較為集中且有大量熱量累積。在測(cè)試中不同化學(xué)反應(yīng)機(jī)理有不同的溫度分布和最高溫度，Saturne 單步機(jī)理與多步機(jī)理時(shí)均溫度結(jié)果對(duì)比如圖8 所示，多步化學(xué)反應(yīng)機(jī)理預(yù)測(cè)結(jié)果明顯優(yōu)于單步化學(xué)反應(yīng)機(jī)理的。RANS 模擬對(duì)流動(dòng)預(yù)測(cè)存在較大誤差，單步甲烷機(jī)理放熱快速，二者在中間截面有相當(dāng)大的誤差。下游3 個(gè)截面主要受對(duì)流擴(kuò)散影響，吻合較好。

圖6 Flame D不同軸向位置時(shí)均溫度徑向分布

圖7 Flame D不同軸向位置時(shí)均甲烷質(zhì)量分?jǐn)?shù)徑向分布

圖8 Saturne單步機(jī)理與多步機(jī)理時(shí)均溫度結(jié)果對(duì)比

進(jìn)一步與相同條件下Fluent 中渦破碎（Eddy Break-Up，EBU）模型、E-A（EBU-Arrehnius）模型和有限速率（Finite-Rate，F(xiàn)-R）模型計(jì)算結(jié)果對(duì)比，不同燃燒模型在不同軸向位置的時(shí)均溫度徑向分布如圖9 所示，可見新軟件精度與E-A 和F-R 湍流燃燒模型相近，優(yōu)于EBU湍流燃燒模型的。

圖9 不同燃燒模型在不同軸向位置的時(shí)均溫度徑向分布

2.2 TECFLAM 冷態(tài)工況研究

URANS 瞬時(shí)軸向速度和時(shí)均軸向速度如圖10所示。瞬時(shí)軸向速度在不停地波動(dòng)，時(shí)均軸向速度存在因強(qiáng)旋產(chǎn)生的中心回流區(qū)以及在噴嘴出口壁面附近因臺(tái)階效應(yīng)產(chǎn)生的小回流區(qū)。

圖10 URANS瞬時(shí)軸向速度和時(shí)均軸向速度

TECFLAM 冷態(tài)p=1.5 Pa 壓力等值面如圖11 所示，顯示為旋進(jìn)渦核結(jié)構(gòu)（Precessing Vortex Core，PVC）。對(duì)于PVC 旋進(jìn)頻率，可在試驗(yàn)和LES 中由湍流能譜[20]通過時(shí)間自相關(guān)得到，也可在URANS 中通過對(duì)固定點(diǎn)的速度時(shí)間序列進(jìn)行傅里葉變換得到[19]。根據(jù)文獻(xiàn)和測(cè)試結(jié)果，選定距噴嘴出口中心軸向30 mm 和徑向20 mm 為固定點(diǎn)。

圖11 TECFLAM 冷態(tài)p =1.5 Pa壓力等值面

點(diǎn)（30 mm，20 mm）瞬態(tài)軸向速度時(shí)域與頻域如圖12 所示。從圖中可見，PVC 結(jié)構(gòu)穩(wěn)定，旋進(jìn)周期約為0.0225 s，旋進(jìn)頻率約為45 Hz。Freitag 等[20]用LES預(yù)測(cè)的旋進(jìn)頻率約為42 Hz；Wegne 等[23]用LES 預(yù)測(cè)的旋進(jìn)頻率為50 Hz，URANS 為35 Hz；隋春杰等[24]用LES預(yù)測(cè)約為40 Hz。本文結(jié)果與前人研究一致。

圖12 點(diǎn)（30 mm，20 mm）瞬態(tài)軸向速度時(shí)域與頻域

TECFLAM 冷態(tài)4 個(gè)軸向位置的時(shí)均速度徑向分布如圖13所示。從圖中可見，對(duì)比距噴嘴出口1、30、60、90 mm 截面的時(shí)均軸向速度，模擬與試驗(yàn)結(jié)果在峰值大小和整體趨勢(shì)上接近。說明新軟件對(duì)復(fù)雜旋流的適用性好，為熱態(tài)模擬奠定了基礎(chǔ)。

圖13 TECFLAM 冷態(tài)4個(gè)軸向位置的時(shí)均速度徑向分布

2.3 TECFLAM 熱態(tài)工況研究

TECFLAM 熱態(tài)時(shí)均溫度如圖14 所示。從圖中可見，高溫區(qū)主要在強(qiáng)旋流產(chǎn)生的中心回流區(qū)內(nèi)，穩(wěn)定在中心鈍體下游。中心回流區(qū)將高溫燃?xì)鈴南掠尉砘兀Y(jié)合中心鈍體，給甲烷和空氣提供了1 個(gè)穩(wěn)定、低速、高溫的穩(wěn)定燃燒區(qū)域。從上游向下游，高溫區(qū)變寬后溫度不斷降低，寬度漸漸變窄。從中心高溫區(qū)到常溫空氣伴流區(qū)，有非常明顯的溫度梯度下降。

圖14 TECFLAM 熱態(tài)時(shí)均溫度

TECFLAM 熱態(tài)4 個(gè)軸向位置的時(shí)均速度徑向分布如圖15 所示。從圖中可見，模擬的最大軸向速度值和整體趨勢(shì)與試驗(yàn)值較為相符，特別是在上游截面，這與射流模擬結(jié)果一致。熱態(tài)工況和冷態(tài)相比，時(shí)均軸向速度正向最大值和負(fù)向最大值都有了提高，說明燃燒的存在使局部流場(chǎng)溫度提高密度減小，對(duì)局部氣流產(chǎn)生了加速效果。下游2 個(gè)截面的速度誤差一部分來源于與射流類似的情況，即RANS 模擬方法和簡(jiǎn)單反應(yīng)機(jī)理；另一部分來源于熱態(tài)旋流各向異性的輸運(yùn)特征。熱態(tài)工況遠(yuǎn)離中心區(qū)域，在徑向大于30 mm 的位置，模擬值與試驗(yàn)值有比較好的貼合。在化學(xué)反應(yīng)劇烈和旋流強(qiáng)的區(qū)域，湍流模型、湍流燃燒模型和化學(xué)反應(yīng)機(jī)理都受到挑戰(zhàn)。

圖15 TECFLAM 熱態(tài)4個(gè)軸向位置的時(shí)均速度徑向分布

TECFLAM 熱態(tài)4 個(gè)軸向位置的時(shí)均速度徑向分布與文獻(xiàn)[24]中改進(jìn)代數(shù)2 階矩模型（Modified Algebraic Second Order Moment，MASOM）和部分?jǐn)嚢璺磻?yīng)器模型（Partially Stirred Reactor，PaSR），文獻(xiàn)[30]火焰面生成流型模型（Flamelet Generated Manifolds，F(xiàn)GM）的結(jié)果對(duì)比如圖16 所示。本文模擬的時(shí)均溫度與試驗(yàn)值在最高溫和整體趨勢(shì)上都保持了較好的一致性。Saturne 模擬結(jié)果從中心火焰區(qū)到低溫伴流區(qū)溫度降低滯后且降低速度過快，說明模擬的火焰面結(jié)構(gòu)過薄，與試驗(yàn)存在差異。本文基于RANS 的模擬結(jié)果優(yōu)于文獻(xiàn)結(jié)果，特別是一些基于LES 的結(jié)果，說明TPDF的優(yōu)越性。

圖16 TECFLAM 熱態(tài)4個(gè)軸向位置的時(shí)均速度徑向分布與文獻(xiàn)[24]、[30]對(duì)比

3 真實(shí)燃燒室算例測(cè)試

基于一致性等研究基礎(chǔ)[31]，應(yīng)用新軟件模擬某型兩相燃燒室，采用377萬非結(jié)構(gòu)化網(wǎng)格，如圖17所示。空氣速度進(jìn)口為121.07 m/s、溫度為696.01 K，出口為充分發(fā)展條件，壁面絕熱無滑移。初始?jí)毫?90216.58 Pa，溫度為696.01 K。中心截面速度、溫度和出口溫度分布如圖18～20 所示。從圖中可見，出口溫度分布的預(yù)測(cè)數(shù)值和試驗(yàn)測(cè)量得到的數(shù)值比較接近，相對(duì)誤差不超過3.8%。

圖17 燃燒室網(wǎng)格

圖18 中心截面速度

圖19 中心截面溫度

圖20 出口截面溫度與出口溫度分布試驗(yàn)數(shù)據(jù)對(duì)比

4 結(jié)論

（1）根據(jù)文獻(xiàn)和模擬結(jié)果分析，典型射流火焰與旋流火焰中間截面誤差主要來源于湍流模型和反應(yīng)機(jī)理，新構(gòu)建的TPDF模型具有較高的準(zhǔn)確性。

（2）對(duì)某型工程燃燒室的模擬，流場(chǎng)與文獻(xiàn)計(jì)算結(jié)果分布一致，旋流器下游燃燒劇烈得到較高的溫度場(chǎng)，下游摻混區(qū)之后溫度逐漸降低，出口溫度與實(shí)驗(yàn)值接近，誤差小于3.8%。在兩相湍流燃燒情況下，模擬結(jié)果驗(yàn)證了新構(gòu)建的TPDF 模型在真實(shí)結(jié)構(gòu)燃燒室典型工況下的準(zhǔn)確性。

（3）TPDF模型可以實(shí)現(xiàn)工程燃燒室性能的仿真，具有一定的發(fā)展?jié)摿凸こ虘?yīng)用價(jià)值。