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

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

0 引言

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

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

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

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

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

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

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

1.2 TPDF簡介

詳細的歐拉隨機場的TPDF 介紹可以參考文獻[5]和[25]?？紤]密度加權(quán)平均的標量聯(lián)合概率密度函數(shù)輸運方程為

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

式中左側(cè)3 項分別為時間項、對流項和化學反應項，均是可以精確計算的。右側(cè)2 項則需要通過梯度擴散模型和小尺度混合模型進行?；玫阶罱K封閉形式的TPDF方程為

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

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

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

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

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

1.4 算例介紹和邊界條件

1.4.1 Flame D算例

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

表1 Flame D算例邊界條件

1.4.2 TECFLAM算例

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

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

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

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

1.5 計算網(wǎng)格和化學反應機理

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

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

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

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

2.1 Flame D算例

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

圖5 Flame D時均溫度

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

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

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

圖8 Saturne單步機理與多步機理時均溫度結(jié)果對比

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

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

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

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

圖10 URANS瞬時軸向速度和時均軸向速度

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

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

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

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

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

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

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

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

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

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

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

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

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

3 真實燃燒室算例測試

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

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

圖18 中心截面速度

圖19 中心截面溫度

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

4 結(jié)論

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

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

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