鈍體火焰穩定器后燃燒不穩定數值模擬

仝照旭，韓啟祥

（南京航空航天大學能源與動力學院，南京 210016）

0 引言

隨著航空發動機向高推重比、高可靠性的方向發展，其加力燃燒室內部的燃燒不穩定現象備受關注[1-3]。鈍體穩定器作為最通用的火焰穩定器之一，廣泛應用于沖壓發動機及渦噴發動機加力燃燒室中[4-5]。鈍體穩定器后燃燒不穩定現象極易受來流條件的影響，隨著燃燒室來流溫度升高，來流馬赫數加大，鈍體穩定器后發生燃燒不穩定的幾率大幅增加，給航空發動機高效、穩定運行帶來重大隱患。因此，燃燒不穩定現象自從被發現以來，廣泛受到工程設計者的關注。但由于燃燒不穩定現象受到多種因素相互作用[6-8]，如化學反應，湍流流動等，且涉及熱力學、傳熱學、聲學、流體力學等多個學科相互交叉，一直難以解決。

國內外學者們通常采用數值模擬對燃燒不穩定現象開展研究。Matthew 等[9]針對軸對稱火焰的縱向振蕩問題，采用2 維模型得到的計算結果與試驗結果吻合程度較好，相較于3 維模型極大簡化了計算資源；劉衛東等[10]應用改進的PISO 算法對液氫/液氧火箭發動機徑向燃燒不穩定現象進行數值模擬，結果表明對該型發動機聲腔應用諧振器阻尼后，壓力擾動波最終被系統抑制；朱旻明等[11]通過U-RANS/PDF 方法對鈍體穩定器后漩渦脫落特性進行分析，相對準確地預測了穩定器尾緣渦脫落頻率，但發生燃燒不穩定時穩定器后的局部雷諾數較高、特征時間較短（一般在μs量級），其流場具有很強的非定常特性及3維效應。這些因素的限制使得國內外關于加力燃燒室燃燒不穩定的研究工作主要圍繞時均特征展開，導致對燃燒不穩定的瞬態特征和作用機理缺乏系統及深入的認識。

本文基于商用軟件Fluent，采用大渦模擬（Large Eddy Simulation，LES）方法對帶鈍體穩定器的模型加力燃燒室進行3 維熱態數值模擬。通過改變來流溫度、來流馬赫數、當量比3 種影響燃燒的關鍵因素來分析其對鈍體穩定器后燃燒不穩定的影響。

1 研究對象與數值計算方法

1.1 幾何模型

選擇加力燃燒室中常見的鈍體穩定器，尺寸如圖1 所示。其主要參數為：槽寬D=45 mm、頂角α=30°、頂角半徑R=12 mm、壁厚b=3mm。大渦模擬對網格精度和計算資源要求極高，有必要對計算域進行簡化。本文采用Euqenio[12]對帶鈍體穩定器的模型燃燒室進行簡化：對鈍體高度方向取45 mm（1D），簡化后鈍體高度方向2 個面取平移周期邊界，3維計算域如圖2所示。

圖1 鈍體穩定器尺寸

圖2 3維計算域

計算域在笛卡爾坐標下表示，選取鈍體穩定器尾緣中心為原點，x方向為長度方向，y方向為寬度方向，z方向為高度方向。3 維計算域尺寸為750 mm×136 mm×45 mm。

1.2 網格劃分與獨立性驗證

采用ICEM 軟件劃分計算域結構化網格。為了精準捕捉鈍體穩定器壁面的流動分離現象，對其周圍劃分邊界層網格，鈍體穩定器壁面第1 層網格高度Δy≤0.02 mm，邊界層網格增長率為1.2，滿足LES 對于近壁網格的要求（y+～1）。整體網格劃分從鈍體穩定器出發，向周圍均勻過渡。網格總體質量≥0.85，網格最小正交性≥0.7。

對冷態流場選取3種數量的網格進行無關性驗證：mesh1（41.8萬）、mesh2（117.2萬）、mesh3（240.1萬）。3種網格穩定器尾緣x=15 mm處時均流向速度如圖3所示。從圖中可見，當網格量從117.2萬增大到240.1萬時，由于網格量改變造成的計算誤差基本可以忽略，綜合求解精度以及計算成本選取117.2萬（mesh2）進行數值計算。中心截面網格剖分如圖4所示。

圖3 3種網格穩定器尾緣x=15 mm處流向時均速度

圖4 中心截面網格剖分

1.3 計算模型與算例驗證

1.3.1 計算模型

湍流模型選擇LES 模型，其中亞格子模型選擇Wale 模型。因為Wale 模型是一種類似Smagorinsky模型的代數模型，但克服了Smagorinsky 模型耗散偏大的缺陷且能夠相對準確地預測層流到湍流的過渡。

離散相模型選擇非穩態顆粒追蹤，離散相顆粒選擇適合液體燃料的液滴（Droplet），阻力系數Cd選擇動力學曳力模型（Dynamic Drag Model）。離散相模型普適性好，幾乎適用所有的流動情形，二次霧化選擇波動破碎模型（Wave Break Up），噴嘴類型選擇實心噴嘴（Solid Cone）。

在鈍體穩定器前150 mm 處垂直來流方向設置6 個噴嘴來提供燃燒所需的燃料，供油位置（黑色圓圈）如圖5所示。

圖5 供油位置

燃燒模型選取有限速率/渦耗散模型（Finite-Rate/Eddy-Dissipation，FR/ED），該模型結合了Arrhenius 公式和渦耗散方程，避免了Eddy-Dissipation 模型出現的提前燃燒問題，并且結合了動力學因素和湍流因素，廣泛用于航空發動機燃燒室的燃燒模擬。采用航空煤油2步總包反應[13]。

本文采用有限體積法求解納維-斯托克斯方程（Navier-Stokes，N-S）方程，壓力和速度耦合采用SIMPLE 算法，壓力項采用2 階迎風格式，動量項采用中心差分格式，時間離散采用中心差分格式。為保證計算精度，時間步長選取1×10-6s。

1.3.2 標準算例驗證

本文通過標準算例來驗證選擇的計算模型能否準確預測穩定器后的流場、氣流的脈動頻率以及火焰的分布形態。其中Giacomazzi 等[12]和Fureby[14]的研究對象為閉口鈍體穩定器，Barry 等[15]的研究對象為開口鈍體穩定器。

1.3.2.1 穩定器后流場結構對比

本文與Giacomazzi 的模擬結果的對比如圖6 所示。從圖中可見，二者在穩定器后流域的漩渦數量均為6個。從左向右前3個漩渦在流域中的位置幾乎一致，后3 個漩渦的空間分布稍有差異，但總體分布位置相似。故本文的計算方法可以相對準確地預測鈍體穩定器后的流場結構。

圖6 本文與Euqenio模擬結果的對比

1.3.2.2 穩定器后監測點的脈動頻率

穩定器后監測點脈動頻率見表1。從表中可見，本文的計算結果與Fureby的試驗結果的誤差為5.7%，且 與Giacomazzi 的模擬結果的誤差相差較小。說明本文的計算方法可以相對準確地預測鈍體穩定器后氣流的脈動頻率。

表1 穩定器后監測點脈動頻率

1.3.2.3 穩定器后火焰形態

鈍體穩定器尾跡火焰形態如圖7 所示。從圖中可見，本文的模擬結果與Barry 等[15]通過高速攝影得到的穩定器后火焰形態的一致性較好。

圖7 鈍體穩定器尾跡火焰形態

1.4 物性參數與邊界條件

計算域中流體介質為空氣或者航空煤油多組分混合物。介質的物性參數設置包括密度、定壓比熱容cp、動力粘度μ和導熱系數λ。對于空氣分別采用分段多項式（Piecewise-polynomial）、Sutherland 公式、動力學理論（Kinetic-theory）進行計算；對于混合物采用混合定律計算；對于密度項ρ均選擇理想不可壓縮氣體。

邊界條件：采用速度進口、壓力出口，壁面采用固壁無滑移邊界，展向2個表面采用平移周期邊界。

2 穩定器后燃燒不穩定的影響因素

給出本文的基準參數，然后在此基礎上對部分參數進行修改來研究鈍體穩定器后燃燒不穩定的影響因素。

基準參數為：速度進口：Ma=0.1，900 K；壓力出口：0.1 MPa；當量比：0.76。

在穩定器后中心截面沿流向方向設置5 個監測點用于分析燃燒室下游火焰的脈動情況。其流向坐標分別為：point-1（x=0.01 m）、point-2（x=0.05 m）、point-3（x=0.10 m）、point-4（x=0.20 m）、point-5（x=0.30 m）。監測點位置如圖8所示。

圖8 監測點位置

2.1 來流溫度的影響

通常情況下，加力燃燒室的進口即是主燃燒室的出口，進入加力燃燒室的氣體為經過主燃燒室燃燒并經過渦輪膨脹作功后的高溫混氣。相較于主燃燒室來說，加力燃燒室來流混氣含氧量較低，燃燒條件相對較差，但從另一角度看，隨著來流溫度的升高，混氣含氧量雖然降低，但溫度提高對燃燒有利。選擇3 種來流溫度600、900、1200 K，分析來流溫度變化對鈍體穩定器后燃燒不穩定的影響。

通過簡單的熱力學分析[16]估算主燃燒室燃燒后的混氣溫度與混氣成分之間的關系（假設燃油完全燃燒），計算得到的結果見表2。

表2 不同溫度時來流氣體成分的質量分數

不同來流溫度時鈍體穩定器后溫度分布如圖9所示。從圖中可見，不同來流溫度對應鈍體穩定器后火焰形態差異較大。當來流溫度為600 K 時，穩定器后的火焰近似呈直線狀，燃燒不穩定性主要體現在內部火焰與外部混氣的交界處，整個交界處呈波浪形火焰，但此時火焰面還未出現大尺度旋渦，穩定器尾部火焰幾乎完全呈現K-H 不穩定；當來流溫度為900 K時，火焰形態由2 部分組成：穩定器尾緣向后延伸一小段長度的剪切層火焰，在剪切層火焰后原本的直線狀火焰逐漸變為由BVK 不穩定主導的旋渦火焰；當來流氣體溫度達到1200 K 時，穩定器尾部火焰幾乎完全呈現由BVK 不穩定主導的旋渦火焰，并且燃燒室下游火焰的擺動幅度很大。

圖9 不同來流溫度時鈍體穩定器后溫度分布

通過上述分析可知，來流溫度是影響穩定器后火焰形態的1 個重要因素。其深層原因之一是由于未氣化的液態煤油不易改變方向，因此在來流溫度較低時，低溫油霧“包裹”在鈍體穩定器后火焰表面抑制其向旋渦火焰轉變；而當來流溫度較高時，液態煤油會被高溫來流氣體加熱而迅速氣化，僅有少部分未來得及氣化的低溫油霧順著穩定器兩側向下游流去，但很快就完全氣化。因此在來流溫度較高時，剪切層火焰向旋渦火焰轉變的流向位置會提前。Gokulakrishnan等[17]同樣以相似結構的鈍體穩定器為研究對象，選擇丙烷為燃料在來流溫度為600 K 時穩定器后火焰形態與本文以液體煤油為燃料、來流溫度為900 K 時的類似。

不同來流溫度時鈍體穩定器后速度分布如圖10所示。從圖中可見，穩定器后火焰由K-H 不穩定主導時，穩定器后流場速度分布分層較為明顯。高速區主要集中在火焰與壁面之間，低速區主要集中在穩定器尾緣；而當火焰由K-H 不穩定向BVK 不穩定轉變時，穩定器后非對稱漩渦將極大影響流場中速度分布，從中心截面的瞬時速度分布可見，當BVK 不穩定出現后，在穩定器下游產生速度很高的渦團，使得原本相對規則的速度分布變得不均勻。

圖10 不同來流溫度時鈍體穩定器后速度分布

2.2 來流馬赫數的影響

來流馬赫數通過影響湍流流動、油氣混合、油氣駐留時間等因素，進而影響鈍體穩定器后火焰的燃燒不穩定現象。分別選取來流馬赫數Ma=0.1、0.2、0.3，分析來流馬赫數變化對穩定器后燃燒不穩定的影響。

不同來流馬赫數時鈍體穩定器后溫度分布如圖11 所示。從圖中可見，從整體上看，3 種來流馬赫數時穩定器后火焰形態非常接近，穩定器后幾乎完全由BVK 不穩定主導的旋渦火焰占據，表明來流馬赫數對穩定器后的火焰形態影響較小。但從細節上看，當來流馬赫數較低時，穩定器后火焰面較為光滑，沒有出現大量褶皺；而隨著來流馬赫數逐漸提高，火焰面逐漸凹凸不平，越來越多的褶皺出現在火焰表面。表明隨著來流馬赫數提高，穩定器后的火焰表面趨于破碎。由于燃油從霧化、蒸發需要一定時間，隨著來流馬赫數提高，低溫燃油在燃燒室的駐留時間縮短。當馬赫數為0.1 時，穩定器尾緣低溫燃油（圖中藍色區域）基本消耗殆盡；而當來流馬赫數提高到0.2 時，少量低溫燃油出現在燃燒室中部；當來流馬赫數達到0.3時，甚至在燃燒室下游也有低溫燃油的存在。

圖11 不同來流馬赫數時鈍體穩定器后溫度分布

由于來流馬赫數不同，不能把三者的速度標尺設為一致。將3 種工況下速度標尺的上限設置為各自來流馬赫數的2 倍，比較3 種來流馬赫數時鈍體穩定器后速度分布的變化。

不同來流馬赫數時鈍體穩定器后的速度分布如圖12所示。從圖中可見，從整體上看，三者的速度分布相似，這與溫度分布的結果趨于一致。表明在由相同火焰不穩定性主導的區域具有相近的溫度分布以及速度分布。但從細節上看，在低來流馬赫數時，穩定器后方的渦團較為完整，尺寸較大且各渦團之間尺度相近；隨著來流馬赫數提高，穩定器后大尺度的渦團在高速氣流的沖擊下逐漸破碎成次級尺度渦團。

圖12 不同來流馬赫數時鈍體穩定器后速度分布

為了分析鈍體穩定器后火焰的脈動特征，設置監測點（圖7）。以來流馬赫數為0.3、監測點5 為例。將該點壓力脈動與放熱脈動的時域圖和頻域圖進行整理，如圖13所示。

圖13 監測點5壓力脈動與放熱脈動時域與頻域

從圖13（a）、（c）中可見，無論是壓力脈動還是放熱脈動，時域曲線均充滿毛刺，表明流場中除了大尺度渦以外還存在小尺度渦。對應到圖13（b）、（d）上可得，頻域曲線除了主頻外，還存在很多與主頻振幅可比的次頻。

將3 種工況下鈍體穩定器后監測點的時域曲線進行FFT 變換得到火焰壓力脈動和放熱脈動的頻譜信息，分別如圖14、15所示。

圖14 不同來流馬赫數時監測點壓力脈動頻域

從圖14 中可見，隨著來流馬赫數提高，燃燒不穩定的壓力脈動頻率呈提高趨勢。這是由于提升來流馬赫數會加快穩定器尾緣非對稱旋渦的生成，而這些非對稱漩渦會對穩定器下游流場產生顯著影響。壓力脈動的振幅同樣隨著來流馬赫數的提高而增大，這可以解釋為提高來流馬赫數會導致流場中各處的平均流速提高，進而導致氣流壓力脈動更劇烈。

對比圖14（a）和圖15（a）可見，發生燃燒不穩定現象時，火焰的放熱脈動和壓力脈動的頻率較為接近，這是因為放熱脈動和壓力脈動均會受流場中大尺度旋渦的影響。從圖15（b）中可見，一方面，當來流馬赫數從0.1提高到0.2時，放熱脈動的幅值呈增大趨勢，這是因為適當提高來流馬赫數使得穩定器下游油滴破碎程度更大，霧化效果更好。另一方面，當馬赫數從0.2 提高到0.3 時，放熱脈動幅值反而在3 種來流馬赫數下最小。這可以解釋為當來流馬赫數過高時，盡管對穩定器后油霧破碎和混合有利，但同時減少了油氣混合的時間和油霧在燃燒室的駐留時間，且過快的氣流會將燃燒不穩定產生的熱量過早帶出燃燒室，不利于熱量的累積。

圖15 不同來流馬赫數時監測點放熱脈動頻域

2.3 當量比的影響

鈍體穩定器后的燃燒不穩定現象與當量比密切相關。選取當量比為0.5、0.6、0.7，分別當量比變化對穩定器后燃燒不穩定的影響。

不同當量比時鈍體穩定器后溫度分布如圖16 所示。從圖中可見，在不同當量比時鈍體穩定器后火焰形態趨于一致，均為剪切層火焰和旋渦火焰共存，2種火焰形態的過渡位置在流向方向位置也幾乎一致。表明當量比從0.5增大到0.7時，不會對穩定器尾部火焰形態造成顯著影響。

圖16 不同當量比時鈍體穩定器后溫度分布

不同當量比時鈍體穩定器后速度分布如圖17 所示。從圖中可見，從整體上看，3 種當量比時穩定器后速度分布基本一致，均在穩定器外側形成速度梯度，在穩定器尾緣形成一定范圍的低速區，且在低速區尾部出現K-H 不穩定與BVK 不穩定的過渡；但從細節來看，當量比為0.5 時Φ= 0.5，在BVK 主導的旋渦火焰區域出現尺寸較大的圓形或橢圓形的渦團；而當量比為0.6Φ= 0.6 時，這些渦團被拉伸，形狀呈條狀；隨著當量比增大到為0.7Φ=0.7時，被拉伸的渦團破碎成更小尺度的渦團。表明隨著當量比的增大，燃燒室下游的燃燒情況更為復雜，其中的火焰尺度更廣。

圖17 不同當量比時鈍體穩定器后速度分布

不同當量比時監測點壓力脈動和放熱脈動的頻域如圖18、19 所示。從圖中可見，隨著當量比從0.5增大到0.7，火焰的脈動頻率變化較小，表明在一定范圍內火焰的脈動頻率對當量比并不敏感。隨著當量比從0.5 增大到0.6，壓力脈動的幅值對應增大，而當量比從0.6 增大到0.7，壓力脈動幅值卻不增反減；放熱脈動的幅值變化趨勢與之基本相同。這可以解釋為穩定器后存在1 個局部最佳當量比，在當量比小于局部最佳當量比時，增大供油量會加快化學反應速率，提高燃燒室內平均溫度；而當量比大于局部最佳當量比時，再增大供油量將會使穩定器后局部富油，過量的低溫油霧不僅會抑制火焰燃燒而且使得燃燒室下游平均溫度降低。

圖18 不同當量比時監測點壓力脈動頻域

圖19 不同當量比時監測點放熱脈動頻域

本文選擇的FR/ED 模型只能用于單步或雙步反應，對于湍流燃燒的細節刻畫較為粗糙。根據Porumbel 等[18]的觀點：FR/ED 模型容易高估反應區的厚度，采用該模型計算得到剪切層火焰的長度與厚度可能偏高，這可能是導致當量比改變時穩定器后火焰形態變化較小的一種原因。但同時Porumbel也指出，在穩定器下游由于燃燒過程化學能的添加使FR/ED 模型在遠離穩定器的區域的預測能力有所提升。

3 結論

（1）來流溫度對穩定器后火焰形態有顯著影響。當來流溫度為600 K 時，穩定器后火焰形態由K-H 不穩定主導；而當來流溫度升到高900 K 時，K-H 不穩定與BVK 不穩定共存；當來流溫度達到1200 K 時，穩定器后火焰形態由BVK不穩定主導。

（2）來流馬赫數對穩定器后火焰形態影響較小；提升來流馬赫數會顯著提高火焰的脈動頻率，并且壓力脈動的幅值也有同樣的變化趨勢。但提升來流馬赫數會對燃燒產生2 種相反的影響：一方面，提高來流馬赫數會使得燃油破碎更好，油霧混合更充分，有利于燃燒；另一方面過高的來流馬赫數使得油霧的駐留時間縮短，并且將燃燒產生的熱量迅速帶走，對火焰的放熱脈動產生不利影響。

（3）燃油當量比從0.5增大到0.7時穩定器后火焰形態變化較小。穩定器后存在局部最佳當量比。當超過局部最佳當量比時繼續增大供油量，局部富油會抑制火焰的放熱脈動。