基于超大渦模擬的燃燒室氣動性能仿真研究進展

張宏達 ，韓省思 ，劉太秋 ，朱 健 ，馬宏宇 ，任祝寅

（1.中國航發(fā)沈陽發(fā)動機研究所，沈陽 110015；2.南京航空航天大學能源與動力學院，南京 210016；3.清華大學航空發(fā)動機研究院，北京 100084）

0 引言

燃燒室是航空發(fā)動機的核心部件之一，燃燒室氣動性能對整機性能起到至關(guān)重要的作用[1-2]。為了滿足發(fā)動機對燃燒室日益嚴苛的技術(shù)要求，燃燒室內(nèi)部燃燒組織模式以及流動特性都變得十分復雜。擴壓器減速增壓過程可能面臨強逆壓梯度下流動分離；氣流經(jīng)過多級旋流裝置形成大尺度渦旋結(jié)構(gòu)，一方面促進液體燃料霧化蒸發(fā)，與燃料形成強脈動、非定常混合氣，另一方面生成氣動回流區(qū)駐定火焰；主燃/摻混孔多股射流與火焰筒內(nèi)橫向來流相互作用，形成反向旋轉(zhuǎn)渦對，對湍流混合有重要影響。在流動基礎上，霧化蒸發(fā)、摻混、化學反應、湍流與火焰相互作用等多尺度物理化學過程強烈耦合，共同決定了燃燒室氣動特性，這些物理化學過程的高精度建模和計算一直是國內(nèi)外研究的熱點問題。

燃燒室中霧化蒸發(fā)、摻混和燃燒過程是在湍流流動環(huán)境下發(fā)展演化的，因此流動是燃燒室氣動性能仿真的基礎[3]。湍流基本特征是由于非線性對流過程導致流動參數(shù)呈現(xiàn)隨機脈動，湍流包含眾多渦系結(jié)構(gòu)，不同渦系在長度和時間尺度的跨度巨大，并且隨著雷諾數(shù)增大，各尺度之間的跨度急劇增大。根據(jù)湍流渦系結(jié)構(gòu)被直接求解的比例，湍流模擬方法劃分為直接數(shù)值模擬（Direct Numerical Simulation，DNS）、雷諾平均（Reynolds-Averaged Navier-Stokes，RANS）、大渦模擬（Large Eddy Simulation，LES）和混合湍流模擬方法。工程中廣泛應用的RANS 方法對湍流平均場進行求解，采用模型模化全部湍流脈動信息，計算量小，但精度較差，對于燃燒室中強旋流和非定常流動過程，RANS 不能滿足精細化設計要求。Pitsch[4]指出，LES 計算量介于RANS 和DNS 之間，目前應用于中低雷諾數(shù)非受限空間湍流燃燒計算。由于燃燒室近壁區(qū)湍流尺度小，流動雷諾數(shù)高，僅燃燒室單個頭部的LES計算所需網(wǎng)格量在數(shù)億至數(shù)十億量級，如此高的計算資源消耗限制了LES 在燃燒室仿真中廣泛使用[5]。基于超大渦模擬方法（Very Large Eddy Simulation，VLES）和混合雷諾平均-大渦模擬方法（Hybrid RANS-LES Method）框架建立高精度計算模型和方法，是數(shù)值仿真的重要趨勢[5-6]。Han 等[7-9]發(fā)展的VLES方法解決了傳統(tǒng)LES中過濾網(wǎng)格尺度和求解湍流尺度匹配限制導致的計算效率低問題，實現(xiàn)了湍流多尺度特性、瞬態(tài)演化特性以及網(wǎng)格分辨率之間的耦合建模，VLES依據(jù)渦系結(jié)構(gòu)演化的實時特征，自適應調(diào)整湍流求解和模型模化之間的比例，在保證計算精度的條件下，顯著降低計算成本。

盡管如此，與傳統(tǒng)LES相比，VLES的理論及特性尚未得到廣泛研究和使用。本文系統(tǒng)介紹VLES 的模化理論及其在燃燒室相關(guān)多種物理場景下的使用效果，促進VLES 在航空發(fā)動機燃燒室仿真領(lǐng)域的規(guī)模化應用。

1 超大渦模擬方法

湍流模擬方法對計算資源消耗及模型影響如圖1所示，RANS、LES和VLES方法均通過湍流建模實現(xiàn)流動仿真。需要說明的是，VLES 最早的明確定義由Pope[6]給出，指的是“計算網(wǎng)格尺度太粗糙從而使直接求解的湍流動能小于總湍流動能的80%"。同時Pope給出的LES 含義[6]是“計算網(wǎng)格很精細從而使直接求解的湍流動能大于總湍流動能的80%”。盡管如此，需要特別說明的是，本文介紹的VLES 是在前期基礎上重新模化發(fā)展的新型計算方法。雖然命名一樣，但新的VLES 方法與Pope定義的VLES 方法具有本質(zhì)上的不同。從圖中可見，傳統(tǒng)的湍流模式按計算精度依次為RANS、URANS、hybrid RANS/LES、LES、DNS，在新的模型框架下，湍流模式劃分按計算精度依次為RANS、VLES、DNS，即VLES 方法實現(xiàn)了多種傳統(tǒng)湍流模式的統(tǒng)一，不同模型在實際計算中依據(jù)局部特征自適應光滑過渡和轉(zhuǎn)換。

圖1 湍流模擬方法對計算資源消耗及模型的影響

與RANS 和LES 方法相比，VLES 建模中考慮了更多的特征尺度，不同湍流模擬方法建模中的特征尺度見表1。從表中可見，積分尺度、湍流截斷尺度和Kolmogorov尺度（最小湍流尺度）等3個具有代表性的湍流尺度，能夠整體考慮流場的多尺度特征。RANS方法僅考慮了流場中最大尺度（即積分尺度），未描述湍流小尺度結(jié)構(gòu)；LES 方法考慮了與網(wǎng)格尺度直接相關(guān)的截斷尺度，但模型中沒有體現(xiàn)湍流尺度的影響，因此LES結(jié)果嚴重依賴于計算網(wǎng)格，為得到高精度的計算結(jié)果，需要劃分很細密的網(wǎng)格。

表1 不同湍流模擬方法建模中的特征尺度

VLES 方法的本質(zhì)是通過引入分辨率控制函數(shù)Fr，對傳統(tǒng)的RANS湍流應力進行重新模化，通過雷諾應力衰減可以得到亞格子尺度湍流應力張量[10]。Speziale[11]最早提出了分辨率控制函數(shù)Fr的模型，但原始模型存在以下缺陷：（1）雷諾應力衰減過度，無法依據(jù)局部網(wǎng)格分辨率合理演化為RANS 模式，難以產(chǎn)生有效的壁面模化；（2）無法保證流場計算達到正確的LES模式；（3）當雷諾數(shù)很高時，網(wǎng)格尺度對于分辨率控制函數(shù)的影響將大幅降低，流場計算基本為RANS 模式，無法達到LES 模式。為了克服上述缺陷，Han 等[7-9]提出了一種全新的超大渦模擬理論架構(gòu)，基于湍流能譜（如圖2所示），直接采用3種特征尺度構(gòu)建分辨率控制函數(shù)

圖2 超大渦模擬建模湍流能譜分區(qū)

式中：E（L）為長度為L的湍流具有的湍流動能能量密度。

Han等[7-9]通過理論分析和推導，由模化和總體湍動能譜的積分，得到了分辨率控制函數(shù)的通用函數(shù)形式

式中：β為常數(shù)，取值為0.002；Lc、Li和Lk分別為截斷網(wǎng)格尺度、積分尺度和Kolmogorov尺度，3個湍流尺度的表達式分別為

式中：Δx、Δy、Δz分別為x、y、z方向的網(wǎng)格尺度；k、ω和β*分別為模化的湍流動能、比耗散率和模型常數(shù)，與RANS 中BSLk-ω模型的相應參數(shù)一致[7]；ν為流體運動黏性；Cx為VLES 中重要的模型參數(shù)，是連接VLES和傳統(tǒng)大渦模擬LES之間的橋梁。在最初的VLES模型中[9-10]，通過引入LES 和RANS 模型常數(shù)進行模化，得到模型參數(shù)Cx= 0.61。

考慮燃燒室內(nèi)“受限空間、大分離、強旋流”的流動特點，優(yōu)化VLES 模型對于復雜物理過程的求解，推導出統(tǒng)一的模型參數(shù)模化形式[12]，即依據(jù)流場參數(shù)動態(tài)確定模型參數(shù)

式中：Cμ為標準k-ε模型參數(shù)；Cs為Smagorinsky LES模型參數(shù)|為應變率。

可以看出，與原始模型中Cx= 0.61 不同，新模型中Cx隨流場變化存在時空演化，可以更加準確地描述湍流的演化過程。

VLES方法實現(xiàn)了對復雜湍流流動的高效高精度數(shù)值模擬，在此基礎上，進一步耦合燃燒、霧化、污染物、聲學等過程的建模，可以實現(xiàn)對燃燒室多物理過程的高精度數(shù)值模擬。

需要說明的是，混合RANS/LES 方法存在的一個經(jīng)典問題是RANS 和LES 模型在模式轉(zhuǎn)換中存在“灰區(qū)”現(xiàn)象，在該區(qū)域中計算的湍流存在“對數(shù)率不匹配”問題[5]。該問題導致從RANS 向LES 的過渡過慢，從而使得對混合過程的計算偏慢。VLES方法極大改進了“灰區(qū)”現(xiàn)象，顯著優(yōu)化了“對數(shù)率不匹配”問題。盡管如此，在網(wǎng)格非常粗糙的情況下，VLES方法還存在一定的“對數(shù)率不匹配”問題，對計算結(jié)果的精度有一定影響。

2 燃燒室典型物理過程仿真

2.1 強旋流動超大渦模擬

航空發(fā)動機燃燒室通常采用多級旋流、強旋流等流場組織形式[13]，旋流流動是燃燒室中最基本的流動形態(tài)。由于旋流在流向和切向的流動均占主導，相比傳統(tǒng)的管流、槽道流和射流流動等，旋流的湍流脈動具有更強的各項異性，因此旋流的數(shù)值仿真對湍流模擬方法提出了很大的挑戰(zhàn)。Xia 等[8]采用VLES 方法計算了經(jīng)典的管內(nèi)強旋流動算例；Dellenback 等[14]針對該算例開展了流場試驗，具有詳細的試驗數(shù)據(jù)。計算算例的流動雷諾數(shù)為1.0×105（基于圓管直徑），旋流數(shù)為1.23。計算采用2 套結(jié)構(gòu)化網(wǎng)格，稀疏網(wǎng)格（M1）總數(shù)約為90 萬，加密網(wǎng)格（M2）總數(shù)約為510萬，進一步將計算獲得的統(tǒng)計矩結(jié)果與試驗結(jié)果進行比較，以驗證VLES方法的計算精度。

為了清晰地顯示強旋流動下渦旋生成和破碎的演化過程，采用Q 準則等值面來識別旋流場中3 維湍流大尺度結(jié)構(gòu)。由RANS BSLk-ω模型和VLES 模型計算得到的Q 準則等值面識別的湍流大尺度結(jié)構(gòu)的壓力分布如圖3 所示。從圖中可見，傳統(tǒng)的RANS 模型只能得到宏觀的大尺度結(jié)構(gòu)，無法捕捉到湍流小尺度結(jié)構(gòu)；而VLES 方法能夠捕捉到湍流相干結(jié)構(gòu)的瞬態(tài)演化過程，即在出口處存在1 個典型的旋進渦核（Precessing Vortex Core，PVC），渦核在進入突擴段后由于剪切K-H 不穩(wěn)性和流場膨脹導致的壓力突變產(chǎn)生劇烈的破碎，由大尺度湍流轉(zhuǎn)變?yōu)樾〕叨韧牧鳌Ｍ瑫r，可以發(fā)現(xiàn)中心PVC 結(jié)構(gòu)經(jīng)過破碎后，在下游處的渦旋尺度變小，即渦管直徑明顯變小。此外，疏密2套網(wǎng)格下VLES 方法計算得到的湍流大尺度結(jié)構(gòu)分布基本一致，區(qū)別在于加密網(wǎng)格可以明顯地捕捉到更多的渦旋破碎小尺度結(jié)構(gòu)。

圖3 Q準則等值面識別的湍流大尺度結(jié)構(gòu)的壓力分布[8]

強旋流動下不同下游位置處周向平均速度和脈動速度沿徑向分布不同方法計算結(jié)果與試驗結(jié)果的比較如圖4 所示。圖中橫、縱坐標分別為無量綱距離和無量綱速度，其中D1為進口圓管直徑，Uin為進口平均速度。從圖中可見，流場表現(xiàn)為典型的類Rankin復合渦逐漸向單一剛體渦過渡。比較計算和試驗結(jié)果，可以發(fā)現(xiàn)VLES 方法對強旋流的周向速度預測具有很高的計算精度，與試驗測量的分布吻合良好。傳統(tǒng)的RANS 方法對旋流計算存在非常大的偏差，無法正確地預測出旋流流場以及湍流脈動的空間演化。相比而言，VLES方法對復雜強旋流動下平均速度場、脈動速度場預測以及空間演化具有非常高的精度，即使在相對稀疏的網(wǎng)格分辨率下依然可以保證較高的計算精度。對于周向平均速度的預測，在疏密2 套網(wǎng)格分辨率下VLES方法的計算結(jié)果基本一致。

圖4 強旋流動下不同下游位置處周向平均速度和脈動速度沿徑向分布不同方法計算結(jié)果與試驗結(jié)果的比較[8]

2.2 湍流燃燒超大渦模擬

為了研究VLES 方法預測湍流燃燒問題的可行性[15-16]，發(fā)展了基于VLES 方法耦合火焰面生成流型建表（Flamelet Generated Manifolds，F(xiàn)GM）[17-18]的湍流燃燒模型，基本思想是假設湍流火焰在局部具有1 維層流火焰結(jié)構(gòu)，湍流火焰面是一系列層流火焰面的系綜平均，因此可以將高維組分空間映射到由幾個特征變量（混合物分數(shù)、反應進度變量等）組成的低維流型，在考慮詳細反應機理的條件下，大幅減少了輸運方程的求解數(shù)量，進而明顯降低了計算成本。具體實施過程是基于混合物分數(shù)和反應進度變量構(gòu)建FGM層流數(shù)據(jù)表，通過假定概率密度函數(shù)方法對層流數(shù)據(jù)表的積分來考慮湍流燃燒間的相互作用，從而獲得湍流數(shù)據(jù)表，數(shù)值計算中求解混合物分數(shù)、反應進度變量以及相應方差的輸運方程，通過查詢湍流數(shù)據(jù)表獲得燃燒場信息。

采用VLES 結(jié)合FGM 建表的湍流燃燒模型對美國Sandia 試驗室試驗測量的甲烷/空氣湍流射流火焰（Flame D）開展了數(shù)值計算，并與試驗測量數(shù)據(jù)進行了定量比較。Sandia Flame D 算例（雷諾數(shù)為22400）的燃料物質(zhì)為體積比1∶3的甲烷與空氣完全混合氣體，燃料進口速度約為49.9 m/s，伴流速度約為11.4 m/s。值班火焰是已燃的甲烷與空氣混合氣體，伴流物質(zhì)是純空氣。計算采用結(jié)構(gòu)化網(wǎng)格，網(wǎng)格總數(shù)約為190萬。

不同組分質(zhì)量分數(shù)平均值沿軸線的分布如圖5所示。圖中橫、縱坐標分別為無量綱距離（D2為進口射流圓管直徑）和無量綱質(zhì)量分數(shù)。從圖中可見，VLES方法對燃燒過程主要組分的預測與試驗結(jié)果總體符合良好。不同下游位置處溫度在混合物分數(shù)空間的散點分布如圖6 所示，從圖中可見，VLES 方法預測的散點分布趨勢與試驗結(jié)果基本一致，僅計算的溫度極值比試驗值略高。VLES 計算獲得的瞬時渦量、溫度和分辨率控制函數(shù)分布如圖7 所示，其中實線取為Zst=0.351。從圖中可見，核心射流區(qū)域表現(xiàn)為強湍流脈動，隨著流場向下游發(fā)展，渦系結(jié)構(gòu)的尺度逐漸變大。從圖7（b）、（c）中可見，在絕大部分的化學反應區(qū)域，分辨率控制函數(shù)處于0和1之間，表明當?shù)鼐W(wǎng)格分辨率可以捕捉到大尺度湍流，僅通過模型模化小尺度湍流，此時VLES 表現(xiàn)為近似大渦模擬求解模式。在射流剪切層和下游的火焰外沿區(qū)域，分辨率控制函數(shù)接近于1，表明計算網(wǎng)格的截斷濾波尺度大于當?shù)氐耐牧鞒叨龋藭rVLES 表現(xiàn)為非定常雷諾平均求解模式。綜上所述，可知VLES 方法可以根據(jù)渦系結(jié)構(gòu)演化實時特征實現(xiàn)多種湍流求解模式的轉(zhuǎn)變，能夠較準確地預測湍流火焰中非定常燃燒過程。

圖5 不同組分質(zhì)量分數(shù)平均值沿軸線的分布[15]

圖6 不同下游位置處溫度在混合物分數(shù)空間的散點分布[15]

圖7 VLES計算的瞬時渦量、溫度、分辨率控制函數(shù)分布[15]

2.3 完整霧化過程超大渦模擬

航空發(fā)動機燃燒室中所用燃料多為液體燃料，液體燃料進入燃燒室經(jīng)歷了一次霧化和二次霧化過程。模擬液體燃料的完整霧化過程面臨諸多難點，包括氣液兩相拓撲界面構(gòu)型的捕捉、液柱變形和破裂、液帶和液絲向液滴的破碎演變、湍流流動和液滴的相互作用等。黃子威[19]發(fā)展了基于VLES 方法耦合VOFDPM 混合霧化計算方法的完整霧化過程仿真模型，實現(xiàn)了燃料霧化從連續(xù)液體到離散液滴的全流程數(shù)值模擬。

采用新發(fā)展的霧化過程仿真模型對經(jīng)典的橫向來流下液柱霧化過程開展了高精度數(shù)值計算，并與公開文獻中的試驗結(jié)果[20]和大渦模擬計算結(jié)果[21]進行了詳細比較。計算算例中氣相為空氣，速度分別為77.89、110.0 m/s，液相為液態(tài)水，速度為8.6 m/s，對應的韋伯數(shù)分別為100、200。為了更好地模擬二次破碎的過程，破碎模型采用了開爾文-亥姆霍茲和瑞利-泰勒（Kelvin-Helmholtz Rayleigh-Taylor，KHRT）模型。

韋伯數(shù)100 工況下VLES 預測的完整霧化過程如圖8 所示。從圖中可見，在初始區(qū)域形成了薄片狀液柱，隨后液柱破裂成液帶和液絲，并在氣動力作用下破碎成液滴，液滴通過二次破碎進一步破碎成更小的液滴。韋伯數(shù)100工況下VLES 計算的流向速度和展向渦量分布如圖9 所示。從圖中可見，在液柱的背風側(cè)，存在著典型的低速回流區(qū)。通過瞬時渦量分布可以發(fā)現(xiàn)，液柱的背風側(cè)表現(xiàn)為強渦旋結(jié)構(gòu)，低速回流區(qū)內(nèi)強湍流運動有助于液柱薄片的破裂和液滴的形成。

圖8 韋伯數(shù)100工況下VLES預測的完整霧化過程[19]

圖9 韋伯數(shù)100工況下VLES計算的流向速度和展向渦量分布[19]

不同韋伯數(shù)下初始射流直徑與液柱開始破碎的液體射流最小流向尺寸之比如圖10 所示。圖中di為液柱開始破碎時液體射流的最小流向尺寸，D3為初始液體射流直徑。從圖中可見，VLES 計算結(jié)果與試驗結(jié)果吻合良好，優(yōu)于文獻[21]中大渦模擬計算結(jié)果。

圖10 不同韋伯數(shù)下初始射流直徑與液柱開始破碎的液體射流最小流向尺寸之比[19]

2.4 碳煙生成超大渦模擬

碳煙（soot）是碳氫燃料燃燒的中間產(chǎn)物，是通過一系列復雜物理化學作用團聚而成的顆粒污染物，碳煙顆粒的排放會形成嚴重的環(huán)境問題。在航空發(fā)動機燃燒室中，碳黑的形成會減少可燃燒的碳，從而降低燃燒效率，并且因為輻射換熱的增加導致火焰筒壁面溫度升高，縮短熱端部件的壽命，因此碳煙排放是燃燒室重要的性能評價指標[22]。此外，在典型的湍流火焰中盡管碳煙顆粒的濃度很低，其體積分數(shù)多在百萬分之一（×10-6）量級，但碳煙的發(fā)射輻射卻是CO2、H2O等氣體發(fā)射輻射的2 倍多，因此準確預測碳煙濃度分布對輻射換熱以及火焰溫度的準確預測極其重要[23]。

碳煙顆粒的形成主要由前驅(qū)體形成、顆粒成核、顆粒表面生長和顆粒氧化4 個步驟組成。對碳煙顆粒進行數(shù)值計算的模型雖然很多，但是計算精度均難以保證，是目前國際上研究的難點問題。采用基于VLES 方法耦合FGM 建表的湍流燃燒模型對湍流火焰進行數(shù)值模擬，結(jié)合半經(jīng)驗的碳煙生成模型，即Moss-Brooks 模型[24]，利用經(jīng)驗公式求解碳煙成核、表面生長和氧化的速率方程，結(jié)合氧化模型和羥基計算模型，初步構(gòu)建了對碳煙顆粒進行仿真預測的數(shù)值模擬方法。

以德國宇航中心試驗研究的乙烯/空氣湍流射流火焰[25]為對象開展了數(shù)值計算，雷諾數(shù)為10000，燃燒工況下全局當量比為0.48，火焰可見高度約為400 mm。VLES 計算和試驗測量的碳煙顆粒體積分數(shù)分布如圖11 所示。從圖中可見，VLES 計算的碳煙顆粒體積分數(shù)峰值約為0.53×10-6，與試驗測量值吻合良好，且體積分數(shù)的空間分布預測與試驗結(jié)果也吻合較好。驗證結(jié)果表明，VLES 方法初步具備了對碳煙顆粒生成的仿真預測能力。

圖11 VLES計算和試驗測量的碳煙顆粒體積分數(shù)分布

2.5 燃燒不穩(wěn)定超大渦模擬

為了滿足低排放要求，民機燃燒室設計通常采用預混或部分預混貧油燃燒組織方式。但是，貧燃預混的燃燒穩(wěn)定性較差，易于激發(fā)熱聲耦合振蕩燃燒模態(tài)，導致燃燒不穩(wěn)定現(xiàn)象的發(fā)生。燃燒不穩(wěn)定具有很強的破壞性，可能伴隨有回火、固體變形等問題，是燃燒室設計面臨的突出難題。

燃燒不穩(wěn)定的數(shù)值計算可分為解耦方法和直接耦合方法2 大類。解耦的燃燒不穩(wěn)定預測方法將非定常燃燒和聲學求解進行解耦分開。非定常燃燒需要大量的數(shù)值計算樣本，以構(gòu)建可靠的火焰描述函數(shù)。若采用大渦模擬計算方法，其計算資源消耗過大。直接耦合計算方法是在可壓縮求解方式基礎上，直接通過高精度的非定常計算獲得燃燒不穩(wěn)定的結(jié)果，即在同一個計算框架內(nèi)一次完成給定工況下非定常燃燒和聲學的耦合計算過程[26]。

在燃燒不穩(wěn)定解耦數(shù)值模擬研究方面，Huang等[27]發(fā)展了基于VLES方法耦合增厚火焰計算方法的燃燒不穩(wěn)定計算模型，實現(xiàn)了對聲波激勵下非定常燃燒過程的準確預測。計算算例為劍橋大學發(fā)展的鈍體駐定乙烯/空氣完全預混火焰，當量比為0.55，雷諾數(shù)約為17000。聲波激勵下非定常火焰動力學特性的VLES 計算結(jié)果和試驗結(jié)果對比如圖12 所示。從圖中可見，在進口激勵過程中，火焰在內(nèi)外剪切層處翻卷，演化形成反向旋轉(zhuǎn)渦對。在此過程中，蘑菇狀火焰輪廓的演變隨相位角變化而不斷發(fā)展。VLES計算結(jié)果很好地復現(xiàn)了試驗中觀測到的火焰演化特性。不同計算方法和試驗測量得到的160 Hz 聲波激勵下熱釋放率響應振幅和相位差對比如圖13 所示。圖中Q'、Q?分別為燃燒的脈動熱釋放量與平均熱釋放量，A為正弦聲波激勵的振幅，圖13（b）縱坐標為聲波激勵下燃燒的瞬態(tài)熱釋放信號與進口速度激勵信號之間的相位差。從圖中可見，VLES 方法的預測精度與大渦模擬精度[28]相當，均與試驗值吻合良好。而非定常RANS 方法雖然預測出了響應非線性的趨勢，但是計算得到的定量結(jié)果與試驗值偏差較大。對于相位差結(jié)果（圖13（b）），VLES 方法預測的相位差隨擾動振幅的變化趨勢與試驗結(jié)果基本一致，而大渦模擬結(jié)果[28]并未很好地預測到上述趨勢。

圖12 聲波激勵下非定常火焰動力學特性的VLES計算結(jié)果和試驗結(jié)果對比[27]（圖形左側(cè)為試驗結(jié)果，右側(cè)為VLES計算結(jié)果）

在燃燒不穩(wěn)定直接耦合數(shù)值模擬研究方面，陳濤等[12]發(fā)展了基于VLES 方法耦合FGM 建表的燃燒不穩(wěn)定直接數(shù)值模擬方法，通過湍流火焰速度封閉的方法將FGM 建表模型和火焰面密度模型（Flame Surface Density，F(xiàn)SD）進行耦合，可壓縮的FGM 建表模型將流體密度和壓力利用理想氣體狀態(tài)方程重新耦合求解，從而保證聲學波動和燃燒過程二者之間的強耦合。采用發(fā)展的燃燒不穩(wěn)定直接耦合計算方法對經(jīng)典的LIMOUSINE 燃燒室[29]開展了數(shù)值計算研究，該燃燒室由荷蘭特溫特大學為研究熱聲耦合燃燒不穩(wěn)定而設計，為部分預混燃燒室。在算例中，燃燒器上游和下游分別由2 個不同寬度的矩形管道組成，中間由1 個正三角體分隔。上游矩形通道為25 mm×150 mm，計算中上游長度設為80 mm。下游矩形通道為50 mm×150 mm，軸向長度為780 mm。算例中熱功率為40 kW，空氣系數(shù)為1.4。燃料甲烷的質(zhì)量流量為8×10-4kg/s，空氣的流量為0.019152 kg/s。研究表明，VLES方法準確預測了燃燒不穩(wěn)定的特征頻率，同時壓力脈動振幅的預測結(jié)果也與試驗值吻合較好。VLES計算和試驗測量的功率譜密度如圖14 所示，不穩(wěn)定周期內(nèi)溫度的演化特性如圖15 所示。從圖中可見，隨著聲學特性的振蕩，伴隨著甲烷在鈍體突擴噴口間隙處的反復阻塞流動過程，形成較弱的回火特征。同時，隨著聲學不穩(wěn)定周期的循環(huán)，新鮮混合物來回往復流動，火焰隨之發(fā)展或閃回，燃料在中心回流區(qū)的反復卷吸過程中，形成了明顯的高溫渦團結(jié)構(gòu)。

圖14 VLES計算和試驗測量的功率譜密度[12]

圖15 不穩(wěn)定周期內(nèi)溫度特性分布[12]

3 燃燒室氣動性能集成仿真

3.1 GTMC雙旋流燃燒室超大渦模擬

選取具有高湍流度的雙旋流燃氣輪機模型燃燒室（Gas Turbine Model Combustor，GTMC）進行數(shù)值仿真研究。國際上基于該模型燃燒室已經(jīng)開展了大量的試驗[30]，具有詳細的流場和溫度場的試驗數(shù)據(jù)。該模型燃燒室主要適用于分析雙旋流流動和燃燒的特性以及仿真模型的驗證與發(fā)展，并已被用作數(shù)值仿真的經(jīng)典驗證算例[31-32]。

針對GTMC 雙旋流燃燒室，計算選取了試驗測量中Flame A 火焰，其中燃料為甲烷氣體，空氣進口和燃料進口的質(zhì)量流量分別為19.74、1.256 g/s。當前外旋流器的旋流數(shù)為0.9，雷諾數(shù)約為58000。計算采用疏密2 套結(jié)構(gòu)化網(wǎng)格，計算網(wǎng)格主要在近壁面和旋流剪切層處進行了加密，稀疏網(wǎng)格（M1）總數(shù)約為320萬，加密網(wǎng)格（M2）總數(shù)約為675萬。為了驗證不同湍流模擬方法對于受限空間高湍流度旋流燃燒室的計算精度，采用VLES 和傳統(tǒng)大渦模擬LES 方法開展了仿真計算，以便相互對比。

GTMC 雙旋流燃燒室中不同位置處軸向速度的平均值沿徑向的分布如圖16 所示。從圖中可見，試驗測量結(jié)果的軸向速度呈現(xiàn)旋流燃燒室中典型的雙峰結(jié)構(gòu)，并逐漸向下游擴張衰減。VLES 方法計算得到的軸向速度與試驗結(jié)果吻合得非常好，可以準確地預測旋流剪切層處的峰值位置及大小。相比而言，LES 方法的計算結(jié)果與試驗值偏差明顯，靠近旋流器出口處預測的軸向速度峰值較試驗結(jié)果明顯偏低，并在下游處雙峰快速衰減。進一步分析發(fā)現(xiàn)，VLES 和LES 計算得到的內(nèi)外旋流器旋流數(shù)差異較小，而分流比具有明顯的差異。VLES計算得到的分流比和試驗測量結(jié)果較為接近，而LES得到的分流比要明顯高于試驗值。表明LES 計算得到的外旋流器質(zhì)量流量相對于實際流量偏高，而內(nèi)旋流器質(zhì)量流量則相對偏低。

圖16 GTMC雙旋流燃燒室不同位置處的軸向速度平均值沿徑向的分布

GTMC 雙旋流燃燒室中不同位置處溫度的平均值沿徑向的分布如圖17 所示。從圖中可見，試驗測量結(jié)果在靠近旋流器出口處溫度較低，而在中心處由于存在旋流產(chǎn)生的中心回流區(qū)駐定火焰，進而存在溫度峰值。此外，兩側(cè)由于燃燒室突擴結(jié)構(gòu)形成的外回流區(qū)因卷吸一部分高溫燃氣而溫度較高。通過比較平均溫度的計算值和試驗值可以發(fā)現(xiàn)，VLES 方法可以很好地預測不同位置處的溫度分布，并較為準確地捕捉溫度峰值的大小與位置，與試驗結(jié)果符合較好。而LES對于溫度預測的偏差非常大，由于LES對于集氣腔中沖擊平板流動以及旋流器中流動預測的累計偏差，導致對于整個旋流場的預測偏差較大，無法正確地捕捉旋流火焰結(jié)構(gòu)。進一步表明LES 方法對于網(wǎng)格分辨率極其敏感，針對復雜燃燒室仿真計算需要極高的計算網(wǎng)格數(shù)和計算資源。LES和VLES計算得到的燃燒室瞬時和平均溫度分布如圖18 所示。GTMC 雙旋流燃燒室是典型的旋流氣動駐定火焰，從圖中可見，火焰懸定在旋流器出口上方，懸定的火焰尖點被稱為火焰駐點（Stagnation Point）。通過比較可以發(fā)現(xiàn)，靠近旋流出口處，LES 預測的火焰出現(xiàn)明顯的壁面附著現(xiàn)象，而VLES 預測的火焰則以一定的角度向下游擴張，呈現(xiàn)典型的V 型火焰拓撲形態(tài)。LES 對火焰拓撲形態(tài)的錯誤預測主要是由于網(wǎng)格分辨率不足，引起旋流場的計算偏差。

圖17 GTMC雙旋流燃燒室不同位置處的溫度平均值沿徑向的分布

3.2 高溫升燃燒室超大渦模擬

張宏達[33]等采用VLES方法對高溫升燃燒室開展了非定常高精度數(shù)值模擬研究，該燃燒室由3 級旋流裝置和雙油路燃油噴射系統(tǒng)組成，較常規(guī)溫升燃燒室的流動和燃燒特性更加復雜。計算采用四面體非結(jié)構(gòu)網(wǎng)格，計算網(wǎng)格主要在頭部裝置和冷卻孔處進行了加密，網(wǎng)格總數(shù)約為1500 萬。選取高溫升、高熱負荷工況作為燃燒室數(shù)值模擬的工況點，關(guān)于數(shù)值方法的詳細介紹見文獻[33-34]。

高溫升燃燒室中心截面瞬時和統(tǒng)計平均速度分布VLES 計算結(jié)果如圖19 所示。為了表征回流區(qū)的形狀，用零速度等值線表示（圖中黑線）。從圖中可見，氣流在燃燒室擴壓器內(nèi)減速后，經(jīng)頭部旋流裝置形成旋流，進而生成中心氣動回流區(qū)a，中心回流區(qū)用于駐定火焰和小工況穩(wěn)定燃燒；同時火焰筒壁面和第3 級旋流氣體之間、以及第2 級旋流和第3 級旋流氣體之間分別形成了突擴回流區(qū)b 和c。主燃孔和摻混孔的多股射流與火焰筒內(nèi)橫向來流相互作用，形成了反向旋轉(zhuǎn)渦對，進一步強化了湍流輸運和混合。

圖19 高溫升燃燒室中心截面瞬時和統(tǒng)計平均速度分布VLES計算結(jié)果[33]

高溫升燃燒室中心截面的瞬時溫度場和統(tǒng)計平均溫度場VLES 計算結(jié)果如圖20所示。從圖中可見，氣流經(jīng)過頭部旋流裝置形成的大尺寸中心氣動回流區(qū)，一方面由于強湍流運動促進了液體燃料霧化蒸發(fā)，另一方面燃料氣體與來流空氣在中心回流區(qū)內(nèi)混合燃燒，導致溫度快速升高，是火焰筒內(nèi)溫度最高的區(qū)域。大孔射流用于強化高溫燃氣與射流空氣間的摻混，調(diào)控獲得滿足設計要求的出口溫度場。高溫升燃燒室出口徑向溫度分布曲線如圖21 所示。通過比較計算結(jié)果和試驗結(jié)果可以發(fā)現(xiàn)，由于傳統(tǒng)RANS 方法不能正確捕捉燃燒室中強旋流和非定常火焰結(jié)構(gòu)，導致出口溫度分布的預測偏差非常大。而VLES 方法的預測結(jié)果與試驗結(jié)果較為接近，較準確地捕捉了出口徑向溫度分布的峰值位置。VLES計算結(jié)果與試驗結(jié)果的偏差體現(xiàn)在火焰筒出口截面徑向相對位置60%以下，分析認為物理建模中未考慮實際的火焰筒冷卻結(jié)構(gòu)，可能導致火焰筒出口近壁氣膜與高溫燃氣間的混合程度不足，從而引起計算偏差。

圖20 高溫升燃燒室中心截面瞬時和統(tǒng)計平均溫度場VLES計算結(jié)果[33]

圖21 高溫升燃燒室出口徑向溫度分布[33]

4 總結(jié)

（1）針對燃燒室內(nèi)多尺度強耦合物理化學過程，介紹了高精度建模和數(shù)值模擬面臨的技術(shù)挑戰(zhàn)，對超大渦模擬的基本思想和模型理論進行了闡述，實現(xiàn)了從實驗室尺度到復雜工程應用場景下湍流流動與燃燒的高精度、高效率的數(shù)值模擬。

（2）根據(jù)超大渦模擬理論架構(gòu)，發(fā)展了VLES 耦合FGM 建表的湍流燃燒模型、VLES 耦合VOF-DPM霧化計算的完整霧化過程仿真模型、VLES 結(jié)合碳煙生成預測模型、燃燒不穩(wěn)定解耦和直接耦合數(shù)值計算模型，實現(xiàn)了燃燒室中旋流流動、燃燒、霧化、碳煙生成和燃燒不穩(wěn)定等典型多物理過程的高精度仿真預測。

（3）選取雙旋流模型燃燒室和高溫升燃燒室，開展了基于超大渦模擬的氣動性能集成仿真研究，為方案篩選、設計優(yōu)化和性能評估提供了工程可用的3 維非定常仿真手段，有望在實際航空發(fā)動機燃燒室設計中規(guī)模化應用。