自適應湍流耦合建表燃燒模型的振蕩燃燒數值模擬

陳濤，徐興平，張宏達，韓省思，*

1．南京航空航天大學 能源與動力學院，南京 210016

2．中國航發沈陽發動機研究所，沈陽 110015

由于環境保護法規的逐漸嚴格，現代燃氣輪機及航空發動機等動力設備的燃燒室設計越來越多采用預混、部分預混的貧燃燃燒方式［1］。但是，貧燃預混燃燒的燃燒不穩定性某些情況下較差，易于激發熱聲耦合振蕩燃燒模態，導致振蕩燃燒問題。振蕩燃燒具有很強的破壞性，甚至可能伴隨有回火、固體形變等多種問題［2］。對于振蕩燃燒物理機理的研究已經成為國際熱點和難點問題之一。

針對復雜燃燒系統的振蕩燃燒試驗技術復雜，存在成本較高、數據不完備等不足。隨著數值模擬技術的快速發展，其逐漸成為燃氣輪機及航空發動機等燃燒設計研發的重要輔助工具之一。對于非穩定燃燒，特別是熱聲耦合振蕩燃燒問題的高精度、高效數值預測顯得尤為重要。程林［3］采用雷諾應力和渦耗散模型，對軸向旋流器外加激勵條件下的振蕩燃燒特性開展了多參數的數值研究，表明振蕩幅值對熱聲耦合具有重要作用，頻率的逐漸增大對釋熱脈動作用并不顯著。Han 等［4-5］采用燃燒和聲波二者解耦的方法，基于低馬赫數大渦模擬，以火焰描述函數耦合低階聲學網絡模型的方式，數值預測了劍橋大學鈍體駐定燃燒室的極限環振蕩模態，計算結果與試驗測量吻合較好。Laera 等［6］通過有限元方法構建了有熱源的非線性火焰模型，通過亥姆霍茲求解器對縱向模態，及環形燃燒室分別進行了極限環參數數值研究，計算結果良好，且可用于復雜外形的工業燃燒器。Shahi 等［7-8］在通用可壓縮求解器上分別測試了SAS（Scale-Adaptive Simulation）湍流模型和渦耗散、PDF（Probability Density Function）建表和BVM（Burning Velocity Model）模型等多種燃燒模型耦合對LIMOUSINE 自激振蕩的數值預測效果，其中BVM 燃燒模型對振蕩燃燒的頻率預測精度較好。該方法對于長火焰的單一縱向模態燃燒室具有較好的應用前景。Chen 等［9］采用大渦模擬和修正的FGM（Flamelet Generated Manifolds）燃燒模型對旋流燃燒室自激工況進行了數值計算，振蕩火焰形態和試驗比較接近，但預測的振蕩頻率低于試驗值。Fredrich等［10］將Pant 等［11］構建的一種新的可壓縮LES（Large Eddy Simulation）-PDF 數值計算方法應用于經典的PRECCINSTA 單旋流燃燒室，并成功捕獲了完全自激的極限環振蕩燃燒模態，預測的熱聲耦合振蕩特征頻率和試驗數據吻合良好。Xia 等［12］則在CFD 和有限元波動方程耦合的方法 中 ，采 用 URANS（Unsteady Reynolds-Averaged Navier-Stokes）湍流模型探究了LIMOUSINE 燃燒室自激振蕩的臨界工況模態，成功預測了特征聲學主頻，但振蕩幅值的預測相比試驗值偏高較多。

前期研究中發現，振蕩燃燒中振蕩幅值的數值預測結果受到多種因素的顯著影響，如湍流模型、湍流燃燒模型、數值耗散、壁面邊界、出口邊界等。通常情況下準確預測振蕩幅值難度較大。除了上述湍流模型和燃燒模型的影響外，文獻中對其他影響因素也做了深入的分析和討論。高階數值格式［13-14］可以更精確地捕捉燃燒室中火焰位置和反應區域，對于徑向聲學模態具有更好的預測，也就是對于高振蕩頻率的幅值預測具有更好的適用性。壁面邊界［15-16］通常和聲學邊界［17-18］共同作用，主要反映為壁面的幾何結構和傳熱特性對聲學邊界阻抗值和燃燒溫度及組分壁面反應的影響?？傊?，對于振蕩燃燒，頻率的預測和多種耦合機制相關，但不論哪種振蕩機制，壓力脈動振幅的預測具有更大的難度。實際研究中，特別是試驗研究中，通常需要精細地設計試驗裝置，以削弱或控制壁面邊界、出口邊界等因素對振蕩燃燒的影響。針對本文研究的經典的LIMOUSINE 部分預混燃燒室振蕩燃燒問題，前期數值研究結果中對數值差分格式、壁面邊界、出口邊界等因素的影響已經提供了一定的參考。但是對于如湍流模型和湍流燃燒模型的影響還有待深入研究。

高精度數值計算得益于更加完備的計算數據，對于振蕩燃燒物理機理的探究具有較明顯的優勢。目前，高精度的自激振蕩燃燒的數值預測較多采用大渦模擬，而大渦模擬對于復雜的燃燒室計算網格要求較高，從而導致計算量很大。基于此，發展高效率且具有高精度的計算方法尤為重要。

目前，振蕩燃燒的數值計算可分為解耦方法和直接耦合方法兩大類。解耦的振蕩燃燒預測方法將非定常燃燒和聲學求解進行解耦。非定常燃燒需要大量的CFD 參數計算以構建足夠樣本和可靠性的火焰描述函數。若采用大渦模擬進行計算，工作量較大，對于工業應用尚有不足。直接耦合計算方法是在可壓縮求解器上，直接通過高精度的非定常計算獲得振蕩燃燒的結果，即在同一個計算框架內一次完成某個工況下的非定常燃燒和聲學的耦合計算過程，對于工業應用具有較大的吸引力［19］。

基于此，本文擬進一步發展耦合的振蕩燃燒數值模擬方法，以在保證計算精度的前提下進一步提升計算效率。在自適應超大渦模擬方法框架下，耦合可壓縮的建表燃燒模型FGM 方法，評估不同的燃燒源項封閉方法。針對經典的LIMOUSINE 部分預混燃燒室開展自激熱聲耦合振蕩燃燒數值研究，驗證其對燃燒不穩定的預測精度，探究在部分預混燃燒下熱聲振蕩的模態特性和內在機制。

1 數學物理模型

1.1 自適應湍流模型SATES

本文采用的自適應湍流模型（Self-Adaptive Turbulence Eddy Simulation，SATES）框架是由Han 和Krajnovi?［20-22］長期發展并廣泛驗證的VLES（Very-Large Eddy Simulation）模型方法。嚴格的理論推導表明，Han 和Krajnovi?［20-22］發展的VLES 方法在網格足夠精細的情況下和LES完全等價，并且可以根據當地計算網格尺度、湍流未解的大尺度、最小的Kolmogorov 尺度三者之間的關系自動調節湍流?；蠼獾谋壤?。由于湍流?；^程自適應調節，因此本文將該方法稱為自適應湍流模擬方法SATES。

自適應湍流模型SATES 方法基于低馬赫數牛頓流體假設求解Favre 過濾的控制方程，其中殘余湍流應力τij的?；贐SL（Baseline）k-ω模型進行求解。通過分辨率控制函數Fr對湍流黏性系數μt進行重新模化，使其達到網格分辨率相關的殘余尺度。分辨率控制函數Fr是SATES方法的核心，具體形式為

式中：Lc、Li和Lk分別為截斷網格長度尺度、積分長度尺度和Kolmogorov 長度尺度，3 個湍流尺度的表達式為

其中：Cx為重要的模型參數，是連接自適應湍流模型SATES 和傳統的大渦模擬LES 之間的橋梁。在最初的SATES 模型中，通過引入LES 和RANS 模型常數進行模化，得到的模型常數為Cx=0.61［20］。

在本研究中，為進一步優化SATES 模型對于復雜物理過程的求解，在前期理論推導基礎上［22］，通過理論分析，得到了更加統一的模型參數模化形式，即采用依據流場參數動態確定模型參數：

式（5）即為本研究新發展的優化的SATES模型參數求解模型。可以看出，與原始模型中Cx=0.61 不同，新模型的參數隨流場變化存在時空演化，從而可以更加準確地描述湍流的演化過程。

在BSLk-ω湍流模型的基礎上，通過分辨率控制函數Fr模化后的湍流黏性系數表達式為

湍動能k和比耗散率ω的輸運方程與經典的BSLk-ω模型完全一致，具體為

如式（1）所示，分辨率控制函數Fr的取值根據網格分辨率的不同在0～1.0 之間平滑取值，以此定量確定湍流的?；潭取．擣r→1 時，μt→，即SATES 將趨近于底層的非定常RANS 湍流模型；當Fr→0 時，μt→0，即SATES方法將趨近于直接數值模擬（DNS）方法。其余模型參數及變量定義可參考文獻［22-24］。

1.2 FGM-FSD 燃燒模型

本文采用詳細化學熱力學建表的火焰面生成流型湍流燃燒模型（FGM）耦合火焰面密度方法（FSD）。FGM 假定湍流火焰的標量演化可以用層流火焰中的若干標量演化進行近似。通過求解簡化的層流火焰方程，采用有限個特征標量描述化學反應過程，結合假定的PDF 模型描述湍流燃燒的相互作用。預先建立包含溫度、組分等信息的化學熱力學表，在考慮了詳細化學反應機理的同時，減少了直接求解標量輸運方程的個數，極大地降低了計算量。

但是FGM 燃燒模型將湍流和燃燒解耦之后，對于強烈的湍流燃燒相互作用過程的計算需要進行驗證和考量?；鹧婷婷芏确椒‵SD 可以保證部分預混或者預混火焰中火焰傳播速度的正確性，從而可以宏觀上較好地處理湍流和燃燒的相互作用。但是通?；鹧婷婷芏饶Ｐ虵SD ?；问捷^多，不同的FSD 模型形式對結果的影響較大。振蕩燃燒屬于強烈的非定常燃燒過程，因此對燃燒模型的可靠性要求更高。Lecocq等［25-26］成功將FGM 和FSD 這2 種方法耦合，并數值探究了PRECCINSTA 燃燒室的預混燃燒工況，取得了較好的效果。

針對火焰面密度模型FSD，Boger 等［27］提出了一種通用的火焰面密度形式：

其中：Ξ為褶皺因子，表征湍流脈動對火焰面的影響，是火焰面密度模型構建的核心參數之一。

本文分別采用Fureby［28］與Zimont 和Lipatnikov［29］提出的模型方法：

式中：模型的其他具體參數及定義可參考文獻［27-30］，這里不再贅述。

在FGM 燃燒模型中，湍流燃燒源項的封閉有多種形式，通常為有限速率方法和湍流火焰速度方法。有限速率方法封閉為

式中：P是反應進度變量和混合物分數的聯合概率密度函數，有限速率采用積分形式封閉源項。

湍流火焰速度方法則具有和火焰面密度相似的定義，具體形式為

通過湍流火焰速度封閉的方法可以將建表燃燒模型FGM 和FSD 方法進行耦合。同時，相比有限速率方法，褶皺因子可以通過校準實現更加精確地預測火焰面位置。本研究中通過式（10）～式（12）對其進行封閉。

通常FGM 燃燒模型采用假定Beta 分布積分獲得湍流化學熱力學表，一般在非絕熱狀態下需要計算聯合混合物分數、未歸一化反應進度變量及其方差以及熱焓即可。但對于不穩定的振蕩燃燒直接計算，需要考慮聲波的求解，因此需要采用可壓縮的求解形式。而解耦方法中的密度、壓力和溫度將不會受到聲學的影響，因此，本文采用考慮可壓縮性的FGM 建表方法?？蓧嚎s的FGM 建表將修正后的流體密度和壓力利用理想氣體狀態方程重新耦合求解，從而可以保證聲學波動和燃燒過程二者之間的強烈耦合。

2 試驗工況與數值設置

本研究將第1 節所述數值計算方法應用于經典的LIMOUSINE 燃燒室［31］，并對其多個工況進行數值模擬研究。LIMOUSINE 燃燒室是荷蘭特溫特大學（University of Twente）為研究熱聲耦合振蕩燃燒發展的部分預混燃燒室，本文采用V3 燃燒室版本，其基本結構如圖1 所示。燃燒器上游和下游分別由2 個不同寬度的矩形管道組成，中間由一個正三角體分隔。上游矩形通道為25 mm×150 mm，軸向長275 mm，根據文獻［7-8］的研究，當發生不穩定燃燒時，聲學在鈍體后流場對燃燒產生較大影響，但對于上游的影響通常較小。因此，本文計算中簡化上游長度為80 mm；下游矩形通道為50 mm×150 mm，軸向長780 mm。正三角型鈍體邊長為22 mm?？諝馔ㄟ^上游注入燃燒室內，純甲烷燃料通過正三角體上的1 mm 孔注入后與空氣快速摻混，并在鈍體面后進行燃燒。本文計算中將甲烷噴孔整體簡化為1 mm 寬的狹縫。

圖1 LIMOUSINE 燃燒室幾何結構Fig.1 Geometric diagram of LIMOUSINE combustion chamber

試驗研究中探究了多種不同的振蕩燃燒工況，燃燒功率包括20、40、50、60 kW，覆蓋了完全發展的極限環振蕩燃燒到穩定燃燒狀態。不同燃燒功率下，燃料流量不同；相同燃燒功率下，空氣流量不同，即空氣系數λ不同（空氣系數的定義同余氣系數）。本文主要計算了3 個典型的振蕩燃燒工況，具體參數見表1。

表1 數值計算工況Table 1 Working condition of numerical calculation

本文采用全結構化網格離散計算域，燃燒室內網格分布比較均勻。在甲烷和空氣摻混的突擴結構處略有加密，燃燒室內最小網格尺寸為0.1 mm?？紤]到本研究中保留了整體下游的結構，燃燒室長度較長，因此燃燒室網格沿軸向以1.15 的比例逐漸增大。依據自適應湍流模型SATES 方法的優勢，本文經過多套計算網格測試后，最終網格總數約為280 萬。正三角型鈍體附近網格示意如圖2 所示。同時，為了驗證計算網格對計算結果的影響，本文將55 萬、140 萬、280 萬網格分別標記為網格M1、M2、M3。

圖2 計算域部分網格示意圖Fig.2 Schematic diagram of partial grid of computational domain

極限環熱聲耦合振蕩燃燒數值預測影響參數較多，聲學邊界是影響參數之一，對振蕩幅值具有明顯的影響［32］。本文研究中空氣進口和甲烷進口均采用質流量進口邊界，保證和試驗中質流量一致。根據Rochette 等［31］的試驗測量結果，出口邊界在壓力出口的基礎上可以設為聲學完全開放邊界；其余固體壁面均為無滑移壁面邊界條件。具體進口參數設置參見表2。

表2 數值計算進口邊界條件Table 2 Inlet boundary conditions for numerical calculation

熱聲耦合振蕩燃燒的重要信息隱含于壓力脈動信號中，本文數值預測中采用和試驗一樣的多點監測方法，分別在上游y=-10 mm，下游y=200，500，650 mm 設置4 個壓力監測點，相對位置如圖3 所示。其中，PMT200 mm 的流向位置和試驗研究中PMT4 的測點位置一致，也是振蕩壓力時序信號的核心測點位置。

圖3 壓力監測點位置示意圖Fig.3 Schematic diagram of pressure monitoring points

數值計算中為保證非定常計算具有較高的采樣頻率，設置非定常計算時間步長為8×10-5s。燃燒模型中詳細化學反應建表采用GRI 3.0 詳細化學反應機理。基于有限體積方法進行離散，時間項采用二階中心隱式格式推進；動量通量采用有界二階中心格式離散；標量通量采用二階迎風格式離散。壓力速度耦合算法采用SIMPLEC算法。

3 結果與討論

3.1 不同湍流燃燒模型預測結果

本研究在振蕩燃燒工況C1 上，對比研究了不同燃燒源項封閉方法下的燃燒模型對振蕩燃燒特征主頻和振幅的數值預測效果。分別將有限速率模型、Zimont 和Fureby 封閉的FSD 模型方法簡記為W1、W2 和W3。

在工況C1 的試驗測量結果中，燃燒室中聲學振蕩特征主頻為234.06 Hz，壓力振蕩振幅為801.09 Pa。3 種不同的源項封閉模型計算的特征主頻快速傅里葉變換（FFT）結果對比如圖4 所示。可以看出，有限速率（W1）、Zimont（W2）和Fureby（W3）對振蕩主頻的預測結果比較，分別為227.22、249.28、236.59 Hz，和試驗測量值236.06 Hz 的相比誤差均較小。但是，3 種源項封閉方法預測的壓力振蕩幅值預測差異較大，有限速率（W1）、Zimont（W2）和Fureby（W3）預測的壓力脈動振幅預測數值逐漸減小，分別為8 050.2、3 081.21、975.28 Pa，相比于試驗測量值801.09 Pa，Fureby（W3）模型方法預測的壓力脈動振幅最接近試驗結果，誤差約為17.9%，而其他2 個燃燒模型的結果誤差＞380%。

圖4 3 種燃燒模型（W1、W2 和W3）預測的特征主頻結果對比（工況C1）Fig.4 Comparison of results of characteristic principal frequency predicted by three combustion models（W1， W2 and W3） （Condition C1）

結合圖5 中壓力和體積熱釋放速率的信號時序圖來看，3 種燃燒模型方法的體積熱釋放速率和壓力變化趨勢比較接近，即二者相位存在強耦合。此時，燃燒火焰和聲波之間形成了正反饋回路，燃燒室處于振蕩自激模態。但是，有限速率（W1）和Zimont（W2）模型的體積熱釋放的振蕩幅值很大，是Fureby（W3）燃燒模型的若干倍。因此，在W1燃燒模型的數值預測過程中，呈現的燃燒振蕩更加激烈，W2 次之，W3 最低也最接近試驗值。同時發現在W1 和W2 的預測中，體積熱釋放率存在負值的情況，且隨著預測振幅偏差的加劇，體積熱釋放速率的負值表現越強，也就是說，在數值計算中，若計算的振蕩過于劇烈，超過了實際振蕩幅值太多，會使計算出現局部熄火的特性。

圖6 和圖7 為3 種燃燒模型在同一時刻下的溫度標量場和軸向速度場對比分布。云圖分布同信號分析結果相一致。W1 燃燒模型在高強度的壓力振幅下，回火嚴重，高溫燃燒向貼著三角型鈍體壁面兩側向甲烷燃料進口蔓延。W2 和W3 模型的振幅預測相對W1 模型較低，受聲波振幅的往復影響，但未形成嚴重的回火現象，火焰高溫燃氣均向下游區域伸展的同時表現為強烈的褶皺結構。但W2 模型預測的火焰穩定比W3模型結果較差，這與W2 模型的振幅較強相關；類似地，W1 和W2 模型預測的火焰渦團卷吸相比W3 模型更加明顯，體現在速度場的中心回流區變化很大。速度向下游發展中，W3 模型受聲學影響主要體現在下游剪切層的不穩定，而W1 和W2 模型的湍流失穩均向上游狹縫靠近，特別是W1 模型，已經在三角型鈍體的上游位置形成了新的回流區，這會對下游形成阻塞，進一步降低燃料和新鮮氧化劑的摻混效果，更加增強火焰的不穩定性。從試驗結果來看，W3 模型預測結果更加接近于真實物理過程，但W1 模型和W2 模型的預測結果也體現了振蕩加劇后的強烈非定?；鹧嫜莼?。

圖6 同一時刻下W1、W2 和W3 模型預測的溫度分布瞬時圖對比結果（工況C1）Fig.6 Comparison results of temperature distribution instantaneous maps predicted by W1， W2 and W3 models at the same time （C1）

圖7 同一時刻下W1、W2 和W3 模型預測的軸向速度分布瞬時圖對比結果（工況C1）Fig.7 Comparison results of instantaneous maps of axial velocity distribution predicted by W1， W2 and W3 models at the same time（Condition C1）

3 種燃燒模型的不同預測結果表明，部分預混的不穩定振蕩燃燒，特別是極限振蕩環狀態下，數值預測結果對湍流及湍流燃燒模型具有很強的敏感性。針對源項的封閉方法，直接體現在相同條件下，不同的源項方法預測的火焰燃燒速度差異較大，燃燒速度過快/過慢，聲學某一模態可能和燃燒過程、高溫已燃氣體的發展相互耦合，熱聲正反饋回路的形成會逐漸增強，如有限速率（W1）和Zimont（W2）2 種源項模型的結果。但這與本文中LIMOUSINE 燃燒室的試驗結果是不一致的。Fureby（W3）封閉源項燃燒模型正確預測了本文中的熱聲耦合振蕩燃燒模式，即燃燒速度相對較慢，部分燃料在受聲學振幅影響卷入卷出速度回流區，向下游發展的同時相對緩慢地燃燒完全，聲學模態和燃燒過程剪切層相耦合，因此，Fureby（W3）模型預測的振蕩特征幅值和頻率均與試驗結果比較相近。

進一步地在C1 振蕩燃燒工況的W3-Fureby湍流燃燒模型上，分析了計算網格對計算結果的影響。計算網格M1、M2 和M3（分別包含55 萬、140 萬和280 萬網格）的計算對比如圖8 所示。網格M1 預測的結果為頻率217.69 Hz，壓力脈動振蕩幅值為670.35 Pa；網格M2 預測的結果為頻率211.86 Hz，壓力脈動振蕩幅值為582.37 Pa。粗糙網格M1 和M2 預測值比較接近，相比試驗值均偏小。網格較為粗糙時，對振蕩燃燒中強烈非定常性的渦系捕捉變弱，無法捕捉小尺度的渦系演化。因此壓力脈動的振幅精度變差。得益于SATES-Fureby 湍流燃燒模型，粗糙網格M1、M2 的預測結果也比較合理，但M3 的預測精度更高，特別是振蕩頻率的預測?？紤]到不同工況的渦系特征差異，以較細的M3 計算網格開展后續多工況的振蕩燃燒計算具有更高的計算精度。

圖8 不同網格在工況C1 下數值預測結果對比Fig.8 Comparison of numerical prediction results of different grids under Condition C1

3.2 不同工況數值結果

根據3 種不同燃燒模型源項封閉方法的計算結果對比，發現Fureby 模型在振蕩燃燒特征主頻和振蕩幅值的預測上均優于其他2 種方法。因此，本文下述章節研究采用了Fureby 褶皺因子模型，對LIMOUSINE 燃燒室繼續開展多工況的熱聲耦合振蕩燃燒研究。根據試驗結果，工況C1 和C3 為極限環振蕩燃燒工況，C3 的振蕩比C1 更加強烈，C2 為火焰穩定和不穩定的過渡狀態［31］。

如表3 所示，給出了3 種不同工況下本文的數值計算和試驗結果對比。3 種工況的特征主頻預測值與試驗吻合均較好，壓力脈動振幅預測結果中工況C3 和C1 誤差相對較大，工況C2 的誤差較小，但最大誤差不大于17.9%。特別地，也給出了針對工況C2 文獻［12］中采用URANS-SST湍流模型結合有限速率封閉源項的FGM 方法預測的結果對比，如表4 所示（Xia 等［12］）?？梢钥闯觯琔RANS 方法預測的工況C2 特征主頻峰值約為300 Hz，振幅約為3 000 Pa，相比本文的數值結果，URANS 頻率和振幅的預測精度誤差均較大?？赡艿脑蚴潜疚牟捎昧薙ATES 自適應湍流模型耦合更加精確的Fureby 燃燒模型，對湍流燃燒的脈動過程預測得更加準確，因此，對于非定常效應強的熱聲耦合不穩定燃燒，大渦方法相比URANS 方法精度更高。同樣，本文也給出了工況C3 的大渦模擬結果［33］，如表5 所示。Hernández 等［33］采用AVBP-Smagorinsky 對工況C3 的相近條件下的燃燒室作了熱聲耦合數值預測。相比于本文的工況C3，僅空氣系數一項由1.2 增大至1.25。因此，文獻［33］預測的振蕩主頻略高于試驗值，本文預測結果和LES 相比，二者預測的差異在合理的范圍之內，均與試驗結果吻合比較好。

表3 多工況下本研究數值預測與試驗結果［31］對比Table 3 Comparison of numerical prediction and test results［31］ under multiple operating conditions

表4 工況C2 本研究數值預測與URANS［12］結果對比Table 4 Comparison of numerical prediction and URANS［12］ results under Condition C2

表5 工況C3 本研究數值預測與LES［33］結果對比Table 5 Comparison of numerical prediction and LES［33］ results under Condition C3

圖9 為工況C2 中熱釋放速率和平均混合物分數的散點分布圖。熱釋放速率隨著混合物分數的變化呈現較寬的分布范圍，在平均混合物分數0.045 附近較為集中，即火焰燃燒在相對廣泛的摻混范圍中發生，進一步證實本文發展的數值方法較好地預測了LIMOUSINE 燃燒室部分預混的特性，也側面反映了流動摻混過程的準確預測對計算結果的顯著影響。

圖9 工況C2中熱釋放速率隨平均混合物分數的散點分布Fig.9 Scatter plot of heat release rate vs average mixture fraction under Condition C2

圖10 為工況C2 中本文計算結果和文獻［7］的時均流場結果對比，其中圖10（a）為PIV 拍攝結果，圖10（b）為文獻［7］中SAS 湍流模型計算結果?？梢钥闯?，本文的數值方法可以有效地捕捉到振蕩燃燒過渡狀態工況C2 的基本流場結構，包括突擴后典型的中心回流區和側回流區。但根據圖10 中0 速度線分布可知，本文預測的中心回流區長度略短于試驗結果［7］。

圖10 工況C2 中平均流場分布對比結果Fig.10 Comparison results of average flow field distribution in Condition C2

圖11 為工況C1 一個振蕩周期內，燃料甲烷的分布示意圖。該序列顯示了火焰隨著聲學振蕩特性的周期振蕩特征。甲烷自噴嘴進入下游，起初的摻混過程較弱，隨后甲烷經過突擴狹縫后在湍流回流區的作用下，自剪切層開始抖動，強化摻混后進行燃燒；與聲學耦合的特性進一步增強了振蕩的幅值。當正弦運動到負相位時，甲烷被阻塞在鈍體狹縫中一段很小的距離，此時，振蕩特性提前強化了燃料的摻混。當下一個相位來臨時，甲烷流出順暢，又重新在鈍體回流區下游的剪切層位置摻混并燃燒，如此往復。同時，與非定場湍流的中心回流區相耦合，甲烷燃料在中心回流區附近一側會卷吸進入中間位置，另一側則向下游伸長延展，最終不穩定火焰在這樣2 個往復過程中形成極限環模態的振蕩燃燒。

圖11 工況C1 中一個周期內燃料甲烷振蕩特性分布Fig.11 Distribution of fuel methane oscillation characteristics in one cycle under Condition C1

結合圖12 來看，工況C1 在隨著聲學特性振蕩時，伴隨著甲烷在鈍體突擴噴口間隙處的反復阻塞流動過程，形成較為輕微的回火特征。同時，隨著聲學振蕩周期的循環，新鮮混合物來回往復流動，火焰隨之發展或閃回，燃料在中心回流區的反復卷吸過程中，形成了明顯的高溫渦團結構。回火特征的數值預測，也顯示出振蕩燃燒的復雜特性。

圖12 工況C1 中一個周期內溫度振蕩特性分布Fig.12 Distribution of temperature oscillation characteristics in one cycle under Condition C1

在工況C1 中，發生極限環振蕩燃燒模態后，整個狹長的燃燒室呈現軸向聲學模態耦合，如圖13 工況C1 中4 個壓力測點的FFT 分布結果所示，各個測點的特征主頻一致，振幅隨著測點沿流向分布而逐漸降低。PMT1 測點的振幅1 121.1 Pa，作為上游水動力學的特征值，要高于聲學振蕩主頻對應的975.25 Pa。此后，越靠近下游，壓力脈動越弱。

圖13 工況C1 中4 個壓力測點的FFT 分布結果Fig.13 FFT distribution results of 4 pressure measuring points under Condition C1

工況C3 和C1 宏觀特性比較接近。但隨著燃燒功率的增加，振蕩燃燒時壓力脈動更加劇烈，如圖14 所示。工況C3 中試驗測定的聲學振蕩主頻為271.88 Hz、振幅為2 846.97 Pa，本文數值方法預測的結果分別為266.37 Hz、2 489.2 Pa，與試驗值較接近；對于上游的PMT1 水動力學測點位置，數值預測的壓力脈動幅值為2 828.3 Pa，大于下游的PMT200 mm 測點位置數值，但差異不大。在工況C3 中，上下游測點的壓力脈動幅值數值預測結果較為接近，可能為在高功率工況C3 下，數值計算聲學的影響比較強烈，已經傳導至上游鈍體臨近噴口處，實際試驗中PMT1 測點的位置要比數值計算中更加靠近上游空氣進口，對湍流水動力的影響反饋也更加集中。這與工況C1 中的結果是類似的，但由于工況C3 的燃燒功率更高，振蕩環模態更加充分，因此在FFT 的結果中也比工況C1 更加顯著。

圖14 工況C3 中2 個主要壓力測點的FFT 分布結果Fig.14 FFT distribution results of two main pressure measuring points under Condition C3

此外，工況C3 中OH 組分和溫度2 種標量場的周期性分布如圖15 和圖16 所示?？梢钥闯?，LIMOUSINE 燃燒室的自激振蕩燃燒模態激發時，火焰會向后推移一段距離（大約為1～2倍鈍體特征長度）后燃燒充分發展起來，這與試驗中觀察的結果比較相似［31］。OH 的分布也呈現明顯的周期性特征，伴隨著上述分析中的燃料摻混效果的反復，OH 組分在1～2 倍鈍體特征長度的區域里隨著剪切層和中心回流區而由中心向兩側壁面振蕩，隨后在3 倍鈍體特征長的區域內形成與上一段距離相反的擴張或收縮特征，這與湍流大尺度渦系結構在振蕩過程中對火焰的作用效果相吻合。

圖15 工況C3中一個周期內組份OH 標量場的分布特性Fig.15 Distribution characteristics of component OH scalar field in one cycle under Condition C3

圖16 工況C3 中一個周期內溫度標量場的分布特性Fig.16 Distribution characteristics of temperature scalar field in one cycle under Condition C3

回火特征在工況C3 中只有輕微表現，可能是隨著燃燒功率的增強，流量增大后湍流雷諾數升高，湍流對聲學振蕩特性的抵抗能力也隨之增強，對可能造成回火的已燃高溫氣體向上游的阻斷作用更明顯。

圖17 所示為工況C3 一個周期內湍流渦系的分布演化特性云圖。在波峰位置時，渦管特征明顯，湍流向中心回流區卷吸，并沿壁面向后發展，向上游形成肩部回流區，呈現典型的突擴后回流特征。當相位逐漸向波谷改變時，渦管被明顯抑制。

圖17 工況C3 中一個周期內湍流渦系的分布特性Fig.17 Distribution characteristics of turbulent vortices in one cycle under Condition C3

3.3 過渡振蕩模式工況雙頻特性結果

工況C1 和C3 均為極限環振蕩燃燒模態工況，工況C2 較為不同，是穩定和不穩定的臨界狀態。本文對工況C2 利用相同的數值模擬方法作了數值研究。

主要時序信號的監測結果如圖18 所示（圖中HRR 表示體積熱釋放速率），PMT200 mm 的燃燒室壓力時序信號的模態與工況C1 和C3 存在明顯差異，壓力在一個周期內上升時，分成2 段，第1 段斜率較高，第2 段斜率較低，之后以同一斜率下降至波谷。燃燒室更加靠后的2 個壓力測點PMT500 mm 和PMT650 mm 的信號演化則趨勢相同，且基本相位差較小。PMT200 mm 和燃燒室鈍體上游的PMT1 以及熱釋放速率的監測信號間存在明顯的相位差。熱釋放速率曲線存在二次波峰，并且間替分布。在PMT200 mm 壓力脈動從波谷向波峰的振幅過程中，熱釋放曲線既存在波峰，也存在波谷。相比完全極限環振蕩工況C1 和C3，工況C2 作為中間模態涉及的物理機制更為復雜。

圖18 工況C2中壓力和熱釋放速率的信號時序監測變化Fig.18 Signal time sequence monitor change of pressure and heat release rate under Condition C2

試驗中測量到工況C2 發生振蕩的特征頻率為212.81 Hz，振幅為676.54 Pa（在PMT200 mm 位置處測量）。如圖19 所示，工況C2 的數值預測結果與試驗相比，誤差較小。不同于前述工況C1 和C3 的模式，工況C2 存在明顯的雙峰預測結果，且只在PMT200 mm 位置處的聲學特征主頻同試驗保持一致。之后，聲波的能量則隨著向下游發展，逐漸同湍流的水動力學相耦合，在PMT500 mm 位置的高頻壓力振幅升高，隨后在PMT650 mm 處因臨近出口又略有降低。

圖19 工況C2 中4 個壓力測點的FFT 分布結果Fig.19 FFT distribution results of 4 pressure measuring points under Condition C2

結合圖20 中工況C2 的熱釋放速率周期分布云圖來看，燃燒室中最根部位置發生典型的渦團周期振蕩后，在一倍的三角型鈍體特征長度下游，熱釋放速率趨于穩定的貼壁以波浪式發展，更下游位置的燃燒脈動較弱，和聲學的振蕩耦合逐漸減弱，聲學模式轉而和下游的已燃氣體通過湍流的強烈耦合而影響整個燃燒室中低頻和高頻聲學頻率的振幅分布。

圖20 工況C2 的瞬態熱釋放速率周期性分布Fig.20 Periodic distribution of transient heat release rate under Condition C2

工況C2 作為過渡狀態，呈現復雜的熱聲耦合物理機制。如圖21 為工況C2 一個周期內的湍流渦系分布云圖。在一個周期內的渦管變化和圖17 中工況C3 的趨勢相同，均隨著相位波峰波谷的變化呈現增強和抑制狀態。盡管如此，在工況C2 中，存在明顯的雙層渦管發展特征，如圖21中紅圈部分所示。第1 層渦管在靠近上游位置，第2 層渦管自中心回流區中間，連接著第1 層渦管之后向下游發展，2 層渦管在一個周期內的某些時刻中同時存在。這也可能是工況C2 呈現雙峰特征的原因。這與Han 和Morgans［5］的LES 模擬結果觀測到的結果比較相似。

圖21 工況C2 一個周期內的湍流渦系分布Fig.21 Turbulent vortices distribution in one cycle under Condition C2

4 結 論

本文在自適應湍流模擬方法SATES 框架下，發展動態模型參數的優化模型，并結合3 種不同燃燒源項封閉方法耦合可壓縮詳細化學反應建表FGM 燃燒模型，對LIMOUSINE 部分預混燃燒室開展了振蕩燃燒數值預測和振蕩模態研究，得到主要結論如下：

1） 熱聲振蕩燃燒的數值預測結果對湍流燃燒封閉源項方法具有較高的敏感性，不同的源項封閉方法由于湍流火焰燃燒速度預測結果不同，直接影響聲學模態的耦合方式及振蕩燃燒的結果。

2） 在自適應湍流模型SATES ?；蚣芟拢邢匏俾誓Ｐ停╓1）、Zimont（W2）和Fureby（W3）3 種燃燒模型對振蕩特征主頻預測結果均較好，但有限速率和Zimont 方法預測的壓力振幅結果遠高于試驗值，Fureby 模型的結果則與試驗值吻合良好。

3） 部分預混的振蕩燃燒數值預測精度一方面受燃燒模型影響，另一方面受湍流摻混效應的影響，特別是對振蕩幅值的影響較大。燃燒過快或過慢可能會由于主體對象不同（如已燃氣體或燃燒釋熱），引起不同類型的熱聲耦合過程。

4） 不同工況下的振蕩燃燒模態預測存在差異。完全振蕩模式下數值預測趨于單一的軸向模態特征主頻分布；過渡振蕩模式下，數值預測的模態存在明顯的雙峰分布特性。

5） 振蕩燃燒的單峰和雙峰特征呈現不同的渦系耦合機理，雙峰特征中存在著雙層渦旋管結構。