鈍體穩(wěn)定火焰尾流特性的數(shù)值模擬研究*

李鈺航，李文龍，吳寶元

(1.液體火箭發(fā)動機(jī)技術(shù)重點(diǎn)實(shí)驗(yàn)室，西安航天動力研究所，西安 710100；2.航天推進(jìn)技術(shù)研究院，西安 710100)

0 引言

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

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

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

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

1 數(shù)學(xué)物理模型

1.1 流動控制方程

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

(1)

(2)

(3)

(4)

上式均經(jīng)過濾波操作。其中，DZ與DC分別為混合分?jǐn)?shù)及進(jìn)度變量的擴(kuò)散系數(shù)；Dt為相應(yīng)的湍流擴(kuò)散系數(shù)；ωC為進(jìn)度變量源項(xiàng)；τij為亞格子應(yīng)力項(xiàng)，其表達(dá)式為

(5)

(6)

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

1.2 湍流模型

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

(7)

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

(8)

(9)

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

1.3 燃燒模型

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

(10)

(11)

其中，χC為進(jìn)度變量C的標(biāo)量耗散率，定義為

(12)

(13)

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

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

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

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

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

表1 算例工況設(shè)置

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

2 計算結(jié)果及分析

2.1 數(shù)值模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)的校驗(yàn)

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

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

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

圖3 流向平均速度分布與實(shí)驗(yàn)結(jié)果對比

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

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

(a)Case C1

(b)Case H2

圖5 不同位置的徑向平均溫度分布與實(shí)驗(yàn)結(jié)果對比

2.2 時均流場分析

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

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

圖7給出了兩個當(dāng)量比不同的燃燒工況的平均溫度分布?？梢?，對于低當(dāng)量比工況H1，火焰最高溫度約1530 K，高溫產(chǎn)物集中于鈍體后方的回流區(qū)內(nèi)，在下游位置燃燒產(chǎn)物與兩側(cè)低溫來流相互摻混，使得流場的平均溫度降低至1200 K左右。對于當(dāng)量比更高的工況H2，火焰最高溫度約1700 K。由于高當(dāng)量比可燃?xì)獾娜紵艧岣鼮槌浞?，高溫區(qū)域一直延伸至下游。此外，H2工況時均溫度場的剪切層厚度較薄，且?guī)缀跗叫杏趤砹鞣较颉Ｏ啾戎?，H1工況時均溫度場的剪切層更厚且呈現(xiàn)較大擴(kuò)張角，這是由于火焰鋒面在橫向出現(xiàn)大幅度脈動導(dǎo)致的，與剪切層內(nèi)的渦動力學(xué)特性相關(guān)，將在下文進(jìn)一步討論。

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

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

2.3.1 流場渦結(jié)構(gòu)分析

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

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

圖8 速度梯度張量第二不變量等值面識別的渦結(jié)構(gòu)

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

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

(14)

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

圖9 渦量源項(xiàng)的橫向分布

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

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

(15)

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

(16)

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

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

(a)Case C1

(b)Case H1

(c)Case H2

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

3 結(jié)論

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

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

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

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