李鈺航，李文龍，吳寶元

(1.液體火箭發(fā)動機技術重點實驗室，西安航天動力研究所，西安 710100；2.航天推進技術研究院，西安 710100)

0 引言

沖壓發(fā)動機一般工作在飛行馬赫數(shù)大于2的條件下，來流空氣經過進氣道壓縮之后，速度仍然很高，遠大于燃料的火焰?zhèn)鞑ニ俣?。因此，沖壓發(fā)動機燃燒室需要采用特殊的火焰穩(wěn)定結構，保證燃料在燃燒室內有足夠的停留時間，使得霧化、摻混、燃燒等復雜過程得以穩(wěn)定地進行[1-2]。

鈍體穩(wěn)焰器具有結構簡單、工作可靠等優(yōu)勢，在工程中得到了廣泛應用[3]。其原理是利用鈍體后方的氣流低速區(qū)駐定火焰，并以此作為點火源點燃剩余燃料，從而在高速氣流中實現(xiàn)穩(wěn)定燃燒。然而，在實際工程研制過程中，由于設計的不合理，常會出現(xiàn)火焰吹脫、穩(wěn)焰器燒蝕以及燃燒不穩(wěn)定等問題。為了解決這些問題，有必要對穩(wěn)焰器尾流的流動與燃燒過程進行深入研究[4]。鈍體穩(wěn)焰器尾流的流動結構較復雜，其中低速回流區(qū)與高速主流之間存在剪切層，具有流動本質不穩(wěn)定特性[5-6]。早期的實驗研究觀察到振蕩燃燒時，穩(wěn)焰器尾流渦結構與穩(wěn)定燃燒時存在明顯區(qū)別，并提出了將非定常渦脫落作為燃燒不穩(wěn)定的重要誘因[7-8]。隨后，針對鈍體預混燃燒的一些研究，又對比了冷態(tài)與熱態(tài)[9-10]、不同鈍體形狀[11]等條件下流動的差異，指出化學反應放熱能夠改變尾流渦結構。近期對于鈍體尾緣噴注燃料的非預混火焰的研究[12-13]又表明，反應區(qū)與剪切層的相對位置對于尾流結構具有顯著影響，燃料和反應放熱區(qū)的分布是火焰動力學的決定因素。

由于實驗研究周期長、費用高且獲得的數(shù)據(jù)有限，計算流體力學方法已成為鈍體穩(wěn)焰器研究的重要手段。已有的研究表明，傳統(tǒng)的雷諾平均湍流模擬(RANS)方法在流動分離區(qū)域嚴重高估了渦粘性；大渦模擬方法(LES)則對壁面附近網(wǎng)格尺寸的要求很高[14]；分離渦方法(DES)能夠在對近壁區(qū)流動?；耐瑫r，準確處理流動分離區(qū)的渦粘性，適合于對鈍體穩(wěn)焰器的數(shù)值模擬[15]。湍流燃燒方面，基于火焰面概念的燃燒模型能夠在考慮詳細化學反應機理的同時，通過預先建表過程減小計算量，適用于預混燃燒、非預混燃燒及部分預混燃燒的模擬。

本文以Volvo鈍體穩(wěn)定火焰的試驗構型為對象進行數(shù)值模擬，采用改進延遲分離渦模型(IDDES)與火焰面生成流形方法(FGM)處理湍流流動與化學反應，研究鈍體穩(wěn)焰器的冷態(tài)與熱態(tài)流動特性，分析燃燒放熱對時均流場及鈍體尾流中非定常渦結構的影響。

1 數(shù)學物理模型

1.1 流動控制方程

針對湍流預混燃燒問題，本文采用的控制方程包含連續(xù)性方程、動量方程、混合分數(shù)以及進度變量輸運方程：

(1)

(2)

(3)

(4)

上式均經過濾波操作。其中，DZ與DC分別為混合分數(shù)及進度變量的擴散系數(shù)；Dt為相應的湍流擴散系數(shù)；ωC為進度變量源項；τij為亞格子應力項，其表達式為

(5)

(6)

式中Δ為濾波尺度；CZ為常量系數(shù)，CZ=0.07。

1.2 湍流模型

為了對鈍體后方的大范圍分離流動進行準確的建模，同時兼顧計算量，本文選取基于k-ωSST的改進延遲分離渦湍流模型：

(7)

-ρβω2-ρ(F1-1)CDkω+Sω

(8)

(9)

其余參數(shù)可參考文獻[16]。

1.3 燃燒模型

本文研究鈍體穩(wěn)焰器的預混火焰，為了考慮詳細化學反應機理，采用火焰面生成流形方法(FGM)處理湍流預混燃燒。該方法在進度變量空間上求解一維預混層流火焰面方程：

(10)

(11)

其中，χC為進度變量C的標量耗散率，定義為

(12)

(13)

式中A為反應速率常數(shù)；ρu為未燃混合氣密度；St為湍流火焰速度，定義可參考文獻[17]。

湍流化學交互作用(TCI)通過概率密度函數(shù)方法處理。

1.4 計算域及網(wǎng)格劃分

針對Volvo鈍體冷流與燃燒試驗[18]進行數(shù)值模擬研究，計算域及其邊界條件設置如圖1所示。丙烷燃料在鈍體上游噴入并與空氣充分混合，保證氣流經過鈍體時為完全預混狀態(tài)。從鈍體上游0.32 m開始選取1.0 m長度的長方體計算域，其橫截面為0.12 m ×0.10 m的矩形。鈍體為寬度0.04 m的實心正三棱錐。流場域上下實體面為絕熱無滑移壁面，左右兩側采用周期性邊界條件以模擬展向無限長的結構。入口條件給定流量，出口采用壓力出口邊界。采用SIMPLE算法進行非定常計算，時間步內迭代至殘差小于10-4時，認為當前時間步收斂，計算時間步長為2×10-6s，時均統(tǒng)計時間為1 s。

設置三種計算工況，分別對應于冷態(tài)流動以及兩種不同當量比的預混火焰，具體如表1所示。

表1 算例工況設置

對計算域的網(wǎng)格劃分采用全局六面體網(wǎng)格，具體如圖2所示，網(wǎng)格總量約為1 600 000。為保證計算的準確性，對壁面附近以及鈍體尾緣剪切層位置進行網(wǎng)格加密，所有壁面法向第一層網(wǎng)格特征尺寸y+<10。

2 計算結果及分析

2.1 數(shù)值模擬結果與實驗數(shù)據(jù)的校驗

為了驗證本文計算結果的正確性，將冷態(tài)與熱態(tài)流動模擬結果中多個位置的軸向速度及溫度分布與實驗測量數(shù)據(jù)進行了對比。

圖2 Volvo算例的網(wǎng)格劃分

圖3展示了冷態(tài)流場中心軸線上軸向時均速度與實驗的對比結果，數(shù)值計算分別采用了IDDES方法與k-ωSST模型的非定常雷諾平均(URANS)方法。結果顯示，IDDES與實驗吻合良好，能夠準確預測流動分離區(qū)的速度分布。URANS則大幅高估了分離區(qū)的回流速度以及分離區(qū)下游的主流速度。其主要原因是鈍體尾流剪切層中大尺度脈動的各向異性特征較強，而URANS方法基于線性渦粘假設，無法解析大尺度渦向小尺度渦的能量級串過程。IDDES方法在自由流動區(qū)直接求解大尺度渦，在壁面附近對邊界層進行模化，模擬此類存在局部分離的流動時具有很好的效果。

圖3 流向平均速度分布與實驗結果對比

圖4給出了冷態(tài)工況C1與燃燒工況H2中三個不同軸向位置處的速度分布與實驗測量值的對比結果，兩工況計算結果均與實驗結果吻合較好。由于工況H2存在燃燒放熱，氣體升溫膨脹增加了主流速度，熱物性變化使得剪切層位置的速度梯度更大。

圖5給出了燃燒工況H2流場中兩個軸向位置處溫度分布與實驗結果對比。對于x/D=0.95位置，計算結果對中心回流區(qū)溫度計算準確，而對鈍體尾緣后方溫度梯度的預示過大。主要原因是鈍體近場位置的火焰鋒面由剪切層主導，湍流與燃燒的交互作用較強，概率密度函數(shù)方法對剪切層位置高溫區(qū)域的橫向擴散存在低估。隨著流動向下游發(fā)展，剪切層不斷增厚，溫度徑向變化趨勢更加平緩。因此，x/D=9.4位置的計算結果更準確。

(a)Case C1

(b)Case H2

圖5 不同位置的徑向平均溫度分布與實驗結果對比

2.2 時均流場分析

圖6給出了三種工況的時均軸向速度場。對于冷態(tài)工況，兩側的高速氣流經過鈍體之后，速度衰減較明顯，速度在橫向的分布很快趨于均勻，鈍體后方所形成的回流區(qū)較短。對于燃燒工況，釋熱膨脹對氣流的加速作用使得流道出口的速度更高，同時鈍體兩側氣流也不再減速。此外，燃燒工況中存在的氣流膨脹并未使得回流區(qū)縮短，反而使其加長，這與火焰鋒面兩側物性的劇烈變化及其對于剪切層所產生的影響相關。

圖6 三種工況的時均軸向速度場

圖7給出了兩個當量比不同的燃燒工況的平均溫度分布。可見 ，對于低當量比工況H1，火焰最高溫度約1530 K，高溫產物集中于鈍體后方的回流區(qū)內，在下游位置燃燒產物與兩側低溫來流相互摻混，使得流場的平均溫度降低至1200 K左右。對于當量比更高的工況H2，火焰最高溫度約1700 K。由于高當量比可燃氣的燃燒放熱更為充分，高溫區(qū)域一直延伸至下游。此外，H2工況時均溫度場的剪切層厚度較薄，且?guī)缀跗叫杏趤砹鞣较?。相比之下，H1工況時均溫度場的剪切層更厚且呈現(xiàn)較大擴張角，這是由于火焰鋒面在橫向出現(xiàn)大幅度脈動導致的，與剪切層內的渦動力學特性相關，將在下文進一步討論。

圖7 燃燒工況的時均溫度場

2.3 瞬態(tài)流場分析

2.3.1 流場渦結構分析

鈍體尾流存在著復雜的渦結構，采用速度梯度張量第二不變量(Q-criterion)的等值面，可識別這些特征結構，進而分析鈍體尾流的流動特性以及燃燒放熱所產生的影響。

首先，對冷態(tài)流動進行分析。圖8中的C1工況展現(xiàn)了兩種典型的流動不穩(wěn)定結構。其一是鈍體兩側尾緣剪切層中的Kelvin-Helmholtz(K-H)不穩(wěn)定性，這種不穩(wěn)定的特征表現(xiàn)為高速主流與回流區(qū)交界面之間的二維小尺度渦。此小尺度渦作為界面上的初始擾動，在向下游發(fā)展的過程中不斷放大，導致了交界面的翻轉并產生了包含展向的三維結構。其二是鈍體遠場的Benard-Von Karman(BVK)不穩(wěn)定性，這種不穩(wěn)定的特征表現(xiàn)為尾跡區(qū)存在反對稱分布的大尺度渦結構。值得一提的是，BVK不穩(wěn)定是一種二次不穩(wěn)定性，由一次不穩(wěn)定發(fā)展而來[4]。剪切層中的渦結構在發(fā)展過程中放大，上下兩側的渦團在回流區(qū)截止位置發(fā)生融合，最終形成了BVK不穩(wěn)定的渦結構。

圖8 速度梯度張量第二不變量等值面識別的渦結構

對于熱態(tài)工況，火焰鋒面位于鈍體尾緣剪切層中，燃燒放熱導致的氣體膨脹和物性變化改變了剪切層的動力學特征。從宏觀角度看，熱態(tài)流場中最明顯的變化是尾流中反對稱結構的消失。由圖8的H1工況可見，上下兩側剪切層中卷起兩排并列的小尺度渦。這些渦結構的發(fā)展較為緩慢，不會迅速增長并形成反對稱的BVK渦結構，而是以近似對稱的形式間歇性地向下游對流。對于H2工況，由于當量比的增加，火焰溫度上升，且氣體膨脹比增加，渦量的間歇特征減弱?？傮w上看，燃燒放熱抑制了BVK不穩(wěn)定性，且當量比越大，抑制作用越強。

為了更深入解釋火焰放熱對剪切層的影響機理，借助如下的渦量輸運方程[4,20]來進行分析：

(14)

上式右端各項表示不同作用對渦量的影響。第一項表示由渦結構的拉伸或收縮引起的渦量空間輸運；第二項表示氣體膨脹產生的影響；第三項為斜壓作用生成項，由壓力梯度與密度梯度不共線而引起；第四項為粘性引起的渦量生成。其中，斜壓及氣體膨脹項與粘性項符號相反，互成競爭作用。圖9給出了x/D=1位置各渦量源項沿橫向的分布?？梢娫诩羟袑游恢?，與冷流工況相比，熱態(tài)工況中火焰鋒面的釋熱作用引起了負的斜壓渦量-(p×ρ)/ρ2和膨脹渦量-ω·u，其作用與粘性渦量生成項×[1/(ρ·τ)]相反，所以導致了鈍體剪切流動渦量被抵消或削弱。這便解釋了上述熱態(tài)流場中剪切層渦卷起受到抑制以及BVK不穩(wěn)定消失現(xiàn)象的成因。此外，高當量比工況H2的燃燒釋熱更劇烈。因此，其膨脹渦量相比H1工況更為明顯，導致了前文所述對BVK不穩(wěn)定性更強的抑制作用。

圖9 渦量源項的橫向分布

2.3.2 流動不穩(wěn)定性的功率譜分析

對于上述的兩種流動不穩(wěn)定性，已有的一些研究給出了其相應頻率的計算方法。其中，BVK不穩(wěn)定性的特征頻率為

(15)

根據(jù)本文算例進行估算，相應頻率約在 110～153 Hz范圍內。對于K-H不穩(wěn)定性，特征頻率計算公式為

(16)

式中Srδ為采用當?shù)剡吔鐚雍穸榷x的斯特勞哈爾數(shù)，其值對于湍流邊界層約為0.044，邊界層厚度δ的計算方法可參考文獻[20]。根據(jù)本文算例設置進行估算，相應頻率約在971.3 Hz附近。

為了分析鈍體尾流中的流動不穩(wěn)定性，在鈍體尾緣(y=0.02 m)后方，設置多個測點，以監(jiān)測局部位置的瞬時速度。圖10給出了三種工況中鈍體近場(x/D=1)及遠場(x/D=4)位置的橫向無量綱脈動速度v/U0的頻譜分布。

(a)Case C1

(b)Case H1

(c)Case H2

在冷態(tài)流動工況C1中，兩個測點處的主導頻率均為f=120.9 Hz，該頻率對應于BVK不穩(wěn)定。對于低當量比燃燒工況H1，近場測點主頻837.9 Hz對應于K-H不穩(wěn)定，但同時也存在127.8Hz的低頻突峰；遠場測點則僅剩下121.1 Hz的低頻突峰。這說明H1工況中BVK不穩(wěn)定雖受到一定程度的削弱，但并未完全消失，而是在更靠下游的位置起主導作用。對于高當量比燃燒工況H2，近場測點主頻為952.7 Hz，與上文估算的K-H不穩(wěn)定頻率十分接近；遠場測點主頻降低為640.5 Hz。兩測點均未見BVK不穩(wěn)定對應的低頻特征頻率，表明其已被完全抑制。

3 結論

(1)采用IDDES湍流模型并結合FGM燃燒模型，能較準確地捕捉鈍體穩(wěn)焰器尾流的流動特征。

(2)鈍體冷態(tài)尾流中存在BVK與K-H兩種流動不穩(wěn)定性。其中，BVK不穩(wěn)定是由鈍體兩側尾緣剪切層相互融合而發(fā)展成的二次不穩(wěn)定。

(3)燃燒釋熱產生的體積膨脹與斜壓效應能夠抵消鈍體尾緣剪切流動引起的渦量生成，從而抑制尾流中的BVK不穩(wěn)定結構；來流中燃料當量比越高，抑制作用越明顯。

(4)可以通過改變穩(wěn)焰器幾何結構或者燃料局部當量比來調節(jié)穩(wěn)焰器后方尾流的渦動力學特性，進而消除由非定常渦脫落誘發(fā)的相應頻率范圍內的燃燒不穩(wěn)定，改進燃燒室性能。