隨機場TPDF在非結構網格平臺的開發及檢驗

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

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

0 引言

航空發動機燃燒室幾何結構非常復雜，先進航空發動機設計對燃燒室提出更精確的燃燒組織要求，對燃燒室數值模擬提出了更高的精度要求[1]，對湍流燃燒模型也提出了更高的要求[2]。概率密度函數輸運方程湍流燃燒模型（Transport Probability Density Function Turbulence Combustion Model，TPDF）可以耦合詳細化學反應機理，高精度解析湍流燃燒過程[3]。

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

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

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

1 非結構網格平臺、TPDF開發和算例介紹

1.1 非結構網格平臺

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

1.2 TPDF簡介

詳細的歐拉隨機場的TPDF 介紹可以參考文獻[5]和[25]。考慮密度加權平均的標量聯合概率密度函數輸運方程為

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

式中左側3 項分別為時間項、對流項和化學反應項，均是可以精確計算的。右側2 項則需要通過梯度擴散模型和小尺度混合模型進行模化，得到最終封閉形式的TPDF方程為

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

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

1.3 基于Saturne軟件的TPDF研發

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

圖1 TPDF構建框架

1.4 算例介紹和邊界條件

1.4.1 Flame D算例

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

表1 Flame D算例邊界條件

1.4.2 TECFLAM算例

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

圖2 TECFLAM 燃燒器幾何結構[6]

表2 TECFLAM 冷態算例進口邊界條件

表3 TECFLAM 熱態算例進口邊界條件

1.5 計算網格和化學反應機理

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

圖3 Flame D算例網格

圖4 TECFLAM 算例網格

2 模擬結果與討論

2.1 Flame D算例

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

圖5 Flame D時均溫度

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

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

圖7 Flame D不同軸向位置時均甲烷質量分數徑向分布

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

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

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

2.2 TECFLAM 冷態工況研究

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

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

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

圖11 TECFLAM 冷態p =1.5 Pa壓力等值面

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

圖12 點（30 mm，20 mm）瞬態軸向速度時域與頻域

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

圖13 TECFLAM 冷態4個軸向位置的時均速度徑向分布

2.3 TECFLAM 熱態工況研究

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

圖14 TECFLAM 熱態時均溫度

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

圖15 TECFLAM 熱態4個軸向位置的時均速度徑向分布

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

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

3 真實燃燒室算例測試

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

圖17 燃燒室網格

圖18 中心截面速度

圖19 中心截面溫度

圖20 出口截面溫度與出口溫度分布試驗數據對比

4 結論

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

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

（3）TPDF模型可以實現工程燃燒室性能的仿真，具有一定的發展潛力和工程應用價值。