高東碩，金 捷，沈 碩，李 敏，曾 家，王 方

(北京航空航天大學能源與動力工程學院仿真研究中心，北京 100191)

在航空發(fā)動機和燃氣輪機燃燒器中，旋流火焰以良好的點火性能和穩(wěn)定性得到廣泛應(yīng)用，且旋流火焰能在非常有限的小體積內(nèi)進行高效率的能量轉(zhuǎn)化[1].由于雷諾數(shù)較高，旋流一般處于湍流狀態(tài)，流動發(fā)展過程中常出現(xiàn)多種不穩(wěn)定現(xiàn)象如漩渦破碎和進動渦核(PVC)[2].進動渦核由剪切層的不穩(wěn)定性產(chǎn)生，是旋流流場中的一種不穩(wěn)定的大尺度擬序結(jié)構(gòu)，PVC會影響附近流體使之出現(xiàn)進動特征.Widenhorn 等[3]對一種雙旋流模型燃燒室無反應(yīng)冷態(tài)流進行了數(shù)值研究與激光多普勒測速(LDA)結(jié)果作對比分析，驗證了位于內(nèi)剪切層的頻率為1 514 Hz 的進動渦核的存在.而PVC 引起的速度脈動會導致燃燒場的不穩(wěn)定性，Canepa 等[4]利用3-D 多普勒速度儀對一個雙同向旋轉(zhuǎn)的徑向旋流器引起的不穩(wěn)定流場進行了測量，表明回流區(qū)附近存在大尺度的渦脫落，渦脫落過程會影響燃燒的穩(wěn)定性.

在湍流燃燒模擬中，精度高適用性強的湍流燃燒模型對于湍流燃燒的強烈非線性相互作用的解釋更為合理.LES[5-6](large eddy simulation)作為現(xiàn)階段描述細節(jié)最多而計算量可行的湍流模擬方法[7]，能夠捕捉流場中的低頻振蕩運動，適合用于模擬有非定常特征的旋流流動.PDF[8](probability density function transport equation)是可以耦合詳細化學機理的火焰模型，LES-TPDF 相結(jié)合能夠得到比較詳細的湍流與化學反應(yīng)時間尺度，可以描述詳細的瞬時流場和溫度場，因此在有旋流場中得到廣泛應(yīng)用.

以往文獻表明當量比的擾動與氣動力學引起的渦旋形成和破碎過程會造成燃燒室的不穩(wěn)定，對當量比的擾動的影響研究已經(jīng)有很多，而對PVC 引起的旋渦與燃燒的不穩(wěn)定關(guān)系的研究較少.本文預(yù)期通過對湍流火焰的數(shù)值模擬，分析燃燒場的數(shù)據(jù)，捕捉流場中的局部高溫區(qū)，討論進動渦核對燃燒不穩(wěn)定性的影響，進而分析局部高溫區(qū)演化機理.

1 計算程序與模型

1.1 程序介紹

AECSC(Aero Engine Combustor Simulation Code)是基于LES-TPDF 方法的三維亞聲速兩相湍流燃燒數(shù)值模擬程序，來源于帝國理工William Jones教授研究組，北京航空航天大學數(shù)值仿真中心發(fā)展了氣相、兩相和大規(guī)模并行版本.

AECSC 程序是FORTRAN 程序，采用模塊化編程，拓展性好，計算精度高，原作者Jones 利用該程序檢驗并發(fā)展了許多模型，取得了良好成果.程序采用大渦模擬方法求解連續(xù)方程和動量方程，采用概率密度函數(shù)方法求解能量方程和組分方程.大渦模擬亞格子應(yīng)力模型使用動力Smagorinsky-Lilly 模型，概率密度函數(shù)方程使用歐拉隨機場解法.已有文獻表明其能夠模擬旋流流場，如Jones 等[9]在2013 年使用LES-TPDF 方法對旋流模型燃燒室進行兩相模擬，與實驗對比良好，火焰結(jié)構(gòu)得到準確捕捉.

1.2 計算模型

本文選取的計算對象是一個雙徑向旋流模型燃燒室，由德國航空航天中心的Meier 等設(shè)計[10-11].實驗[11]將測量區(qū)域簡化為燃燒室主燃區(qū)，即在燃燒器的近場和受液體燃料主要影響的區(qū)域.為了方便實驗測量，燃燒室由4 塊石英玻璃構(gòu)成，因此實驗是在光學可達的高壓單扇型燃燒室中進行的.利用LDA測量冷態(tài)流場速度，相位多普勒分析儀(PDA)用于測量霧化液滴的尺寸和速度，使用平面激光誘導熒光(PLIF)方法使火焰區(qū)域溫度分布可視化.

GTMC 的幾何結(jié)構(gòu)如圖1 所示，燃燒室截面為102 mm×102 mm 的正方形，長264 mm，實驗使用的燃料為液態(tài)航空煤油Jet-A.空氣通過中心旋流器和外旋流器進入燃燒室，兩個旋流器的徑向通道均為12 個，空氣質(zhì)量流量比接近1∶1.65，其設(shè)計細節(jié)如圖2 所示.圖3 展示了空氣和航空煤油進口的設(shè)計細節(jié)，煤油由兩個位置相對的燃油管道先進入一個環(huán)形的燃油通道，然后通過圓周均布的 36 個面積0.2 mm2大小的孔進入垂直槽，實驗徑向通道中高溫空氣會在垂直槽中給燃料帶來顯著的溫升和預(yù)蒸發(fā).燃料穿過垂直槽后以31°的傾角進入燃燒室，這會使得燃油的徑向動量很小，避免了燃油從出口處脫離.

圖1 GTMC的幾何結(jié)構(gòu)Fig.1 Geometric structure of GTMC

圖2 GTMC中心旋流器和外旋流器軸向截面圖Fig.2 Axial sections of GTMC center cyclone and outer cyclone

圖3 空氣和航空煤油進口的設(shè)計細節(jié)Fig.3 Design details of air and aviation kerosene inlets

1.3 網(wǎng)格信息

本文使用SolidWorks 軟件創(chuàng)建幾何模型，使用ICEM CFD 軟件劃分三維結(jié)構(gòu)化網(wǎng)格，網(wǎng)格共劃分為144 個區(qū)塊，分別劃分了網(wǎng)格數(shù)量為220 萬、297 萬與440 萬的3 套網(wǎng)格，在Fluent 軟件中使用3 套網(wǎng)格對模型燃燒室的冷態(tài)流場進行模擬，發(fā)現(xiàn)220 萬網(wǎng)格的計算結(jié)果不是很好，而297 萬與440 萬網(wǎng)格的模擬結(jié)果基本一致，綜合考量計算量與精確度，選擇297 萬的網(wǎng)格進行后續(xù)計算[12]，見圖4.在選用的297 萬網(wǎng)格中，軸向(流向)方向是Y 方向，圖4(a)展示了X＝0 截面的網(wǎng)格，在燃燒室軸向方向上由于前段有著非常強烈的兩股氣流的摻混，這里的網(wǎng)格進行了加密，中后段的流動和反應(yīng)都不是很強烈，網(wǎng)格調(diào)得稀疏一些.圖4(b)展示了內(nèi)外旋流器部分的網(wǎng)格細節(jié)，其中藍色和綠色的方形部分為旋流器的空氣進口.燃燒室的內(nèi)外旋流器采用5 層O-block 進行加密.

圖4 GTMC整體和旋流器部分的計算網(wǎng)格Fig.4 Computational grid of overall GTMC and swirler section

1.4 邊界條件

實驗的模型燃燒室工況見表1.本算例的燃料為航空煤油Jet-A，模擬時真實組分使用分子式C12H23代替，程序采用烷烴通用的四步七組分反應(yīng)機理[13].

表1 實驗工況Tab.1 Experimental conditions

2 模擬結(jié)果及高溫區(qū)分析

2.1 燃燒場模擬結(jié)果

按照表1 中的冷態(tài)邊界條件，使用AECSC 程序模擬冷態(tài)流場，實驗測量了無反應(yīng)流的速度分布，取距離噴射器出口軸向高度h＝10 mm 處的軸向、徑向與切向速度與實驗數(shù)據(jù)進行對比，如圖5，3 個時均速度均與實驗測量速度分布趨勢相同，且曲線的極值點位置相近，驗證了AECSC 程序?qū)δＰ腿紵依鋺B(tài)模擬的準確性.

圖5 冷態(tài)工況下軸向、切向與徑向時均速度沿徑向分布Fig.5 Radial distribution of time-averaged velocity in axial，tangential and radial directions under cold condition

對于熱態(tài)工況，先按表1 中邊界條件模擬出沒有噴油的冷態(tài)流場，以此為基礎(chǔ)進行噴油并使用點火模型進行點火和燃燒.實驗沒有對進口燃油的溫度，速度和d32進行測量，考慮到燃油在垂直槽中的溫升，將初始液滴 d32設(shè)置為 20，液滴初始溫度設(shè)為340 K，液滴初始速度為50 m/s.使用一個隨機場求解PDF 模型.在點火成功，燃燒流場基本穩(wěn)定后統(tǒng)計了3 個周期，約30 萬步.在實驗中使用PDA 測量了霧化液滴的尺寸和速度，圖6 展示的是在距離噴射器出口軸向高度h＝10 mm 處液滴的軸向、徑向與切向速度與實驗數(shù)據(jù)的對比，可以看出程序模擬結(jié)果與實驗相差不大，而在內(nèi)回流區(qū)(IRZ)實驗測得的軸向速度為負值，而程序模擬的軸向速度為零，猜測是由于在模擬時采用的是圓錐空心霧化方式，直徑較小的液滴分布在中心，直徑較大的液滴分布在邊緣，當小液滴進入溫度較高的區(qū)域時會快速蒸發(fā)，從而導致以上現(xiàn)象.

圖6 熱態(tài)工況下液滴軸向、切向與徑向時均速度沿徑向分布Fig.6 Radial distribution of time-averaged velocity of droplets in axial，tangential and radial directions under combustion condition

溫度場可以直觀地表征燃燒狀態(tài)和火焰形態(tài)，圖7 展示了流場在X＝0 截面瞬時溫度云圖，表明了沿流向方向溫度的分布情況.圖8 分別展示了在X＝0截面上實驗測量的溫度云圖和AECSC 程序計算得到的時間平均溫度云圖，展示的均是圖7 中粉色方框區(qū)域.實驗中流場溫度由OH 的絕對密度計算得到，而OH 的絕對密度由PLIF 和吸收同時測量，由于OH 濃度隨溫度迅速下降，該方法可獲得的低溫極限在1 500 K 左右，因此實驗溫度中色條覆蓋范圍是1 500 K 至2 200 K.圖8 中程序模擬結(jié)果很好地再現(xiàn)了實驗的V 型火焰，且高溫區(qū)的溫度值與實驗結(jié)果很接近均在2 200 K 附近.模擬結(jié)果很好地還原了低溫的內(nèi)、外回流區(qū)，其中內(nèi)、外回流區(qū)位置見圖9，內(nèi)、外回流區(qū)的大部分區(qū)域溫度低于1 500 K 或在1 500 K 附近，與實驗結(jié)果基本吻合，說明在內(nèi)、外回流區(qū)沒有發(fā)生明顯的化學反應(yīng).

圖7 X＝0截面瞬時溫度云圖Fig.7 Instantaneous temperature map at section X＝0

圖8 X＝0 截面時均溫度云圖，包括AECSC 程序計算和實驗結(jié)果Fig.8 Time-averaged temperature map at section X＝0，including AECSC program calculation and experimental results

空氣氣流流經(jīng)內(nèi)外旋流器，從旋流器出口高速噴出，形成了具有較大的軸向、徑向和切向速度的高速旋流流場.兩級旋流旋向相同，旋流產(chǎn)生的速度梯度可以使兩股旋流得到充分混合，隨著流向距離的增加，旋流沿徑向逐漸向外擴張.圖9 是計算得到的X＝0 截面的軸向速度云圖，展示的是圖7 中黑色方框區(qū)域.在軸向速度云圖中可以清晰地看到回流區(qū)的邊界(vmean＝0)，其中內(nèi)回流區(qū)是由雙旋流器的存在而產(chǎn)生的，因為它位于“V 形”噴射流之間，預(yù)期會影響燃料蒸氣與預(yù)熱空氣的混合，因而具有重要的意義.外回流區(qū)靠近燃燒室壁邊，由于存在突擴，在拐角處與主流形成一個漩渦區(qū).混合氣體的燃燒主要受外回流區(qū)和中心回流區(qū)的控制，在二者之間的剪切區(qū)進行燃燒，因此回流區(qū)的位置可以預(yù)測火焰產(chǎn)生的位置.

圖9 X＝0截面軸向時均速度云圖和等值線Fig.9 Time-averaged velocity map in axial direction and isolines at section X＝0

2.2 局部高溫區(qū)分析

在30 萬步的基礎(chǔ)上繼續(xù)模擬，以200 步為一間隔輸出模擬結(jié)果，總共輸出30 000 步，用于分析瞬時流場的局部高溫區(qū).計算時間步長為1.18×10-7s，統(tǒng)計非定常流動時間為0.003 54 s，即統(tǒng)計非定常流動時間開始于0.034 s 終止于0.037 54 s.

圖10 是在X＝0 截面不同時間的瞬時溫度云圖，展示的也是圖7 中粉色方框區(qū)域，圖中的輪廓線均為燃料蒸氣(C12H23)質(zhì)量分數(shù)在10%的等值線.結(jié)合圖9、10 中可以看出瞬時溫度場一直在變化，燃燒主要發(fā)生在速度剪切層.

圖10 X＝0 截面瞬時溫度云圖和燃料(C12H23)質(zhì)量分數(shù)等值線Fig.10 Instantaneous temperature map and isolines of mass fraction of fuel (C12H23) at section X＝0

本文在分析時，將溫度高于1 700 K 的區(qū)域定義為高溫區(qū)，在軸向Y＝0 至Y＝0.13 m(二分之一軸向長度)的流場中溫度高于1 700 K 區(qū)域占全部區(qū)域的15%左右.在圖11 中X＝0 的中心截面高溫區(qū)主要存在燃燒室的前段 0～0.075 m(燃燒室軸向長度為0.264 m)即Y/D＝0～5 部分，其中D 為中心旋流器出口直徑，D＝0.015 m，Y 是軸向方向.在Y/D＝0～5區(qū)域，分析垂直軸向截面上的高溫區(qū)變化規(guī)律.

圖11 X＝0截面瞬時高溫分布和流線Fig.11 Instantaneous high temperature map and streamlines at section X＝0

2.2.1 進動渦核對高溫區(qū)的影響

從圖11 中可以看出篩選出的瞬時高溫區(qū)區(qū)域與瞬時流線中的渦旋結(jié)構(gòu)位置接近，說明高溫區(qū)的產(chǎn)生與剪切層中的渦團有關(guān)，瞬時流線圖揭示了剪切層中大尺度旋渦的存在.

在圖12 中發(fā)現(xiàn)剪切層左右兩側(cè)中流線所突出的漩渦中心并不位于同一高度，用黑色折線將旋渦中心連接，漩渦中心排列成鋸齒狀，這與旋轉(zhuǎn)對稱有關(guān)，此現(xiàn)象是旋流中PVC 形成的強烈標志.進動渦核是旋流流場中在高旋流數(shù)下典型的不穩(wěn)定性現(xiàn)象，位于IRZ 并繞中心軸進動，由剪切層的不穩(wěn)定性產(chǎn)生.PVC 是不穩(wěn)定的大尺度旋轉(zhuǎn)流動結(jié)構(gòu)，可以促進流場的混合，周圍大部分流體繞其作旋轉(zhuǎn)運動，使之出現(xiàn)進動特征[14].

圖12 X＝0截面瞬時流線(軸向速度染色)Fig.12 Instantaneous streamlines at section X＝0(colored by axial velocity)

圖13 中展示的是X＝0 截面瞬時放熱率云圖，文中瞬時放熱率的定義是單位時間內(nèi)單位體積所釋放的熱量.公式為：

式中：i 代表化學反應(yīng)中涉及的組分；Hi是對應(yīng)組分的焓值；ki是對應(yīng)組分的化學反應(yīng)速率.

在圖13 中發(fā)現(xiàn)瞬時放熱率比較高的區(qū)域與篩選出的高溫區(qū)前半段區(qū)域形狀類似，此處也是漩渦比較密集的區(qū)域，結(jié)合圖14 的軸向渦量圖，表明在渦量比較高的區(qū)域，PVC 對周圍流體的作用顯著，加強了已燃的高溫氣體和未燃的高溫氣體摻混，促進燃燒和放熱.而在圖11 篩選出的后半段高溫區(qū)域放熱率很小，說明高溫區(qū)域的形成主要原因不是燃燒放熱，而是PVC 旋轉(zhuǎn)流動產(chǎn)生的流體摻混.有研究表明放熱增強會導致PVC 頻率的增加和PVC 振蕩的幅值減小[14].從圖14 中也可以看出隨著軸向距離的增加，渦量值減小，PVC 對周圍流體的作用減弱.

圖13 X＝0截面瞬時放熱率云圖Fig.13 Instantaneous heat release rate map at section X＝0

圖14 X＝0截面瞬時軸向渦量圖Fig.14 Instantaneous axial vorticity map at section X＝0

2.2.2 高溫區(qū)的演化

圖15 是相同時間不同方向的PVC 瞬時視圖，采用旋流器出口附近的瞬時壓力p＝0.426 MPa 等值面顯示并以瞬時軸向速度染色.中央噴嘴的旋流器所產(chǎn)生的旋流運動驅(qū)動著旋渦芯，像施加的渦流一樣，整個結(jié)構(gòu)圍繞自身渦軸順時針旋轉(zhuǎn).

圖15 進動渦核展示(p＝0.426 MPa 等值面，軸向速度染色)Fig.15 Precession vortex core visualized(iso-surface p ＝0.426 MPa，colored by axial velocity)

為了展示PVC 在時間上的演變，圖16 進一步展示了相同的瞬時壓力等值面(0.426 MPa)XZ 方向上在4 個不同時間的PVC，螺旋所包圍的流動在軸向速度的作用下進行再循環(huán)，圖16 中漩渦變化了一個周期，PVC 的頻率在1 000～2 500 Hz.

圖16 不同時間下進動渦核展示(p＝0.426 MPa 等值面，軸向速度染色)Fig.16 Precession vortex core visualized at different time(iso-surface p＝0.426 MPa，colored by axial velocity)

在分析垂直軸向截面上的局部高溫區(qū)變化時，首先篩選溫度高于1 700 K 的區(qū)域，并將局部高溫區(qū)域旋轉(zhuǎn)一圈的時間定義為一個周期，在垂直軸向截面Y/D＝0～5 上，高溫區(qū)的變化頻率均在 1 000～2 500 Hz，頻率近似.圖17 中展示了軸向Y/D＝1.4截面上t＝0.034～0.034 72 s 的局部高溫區(qū)變化圖，圖中高溫區(qū)從1 到6 旋轉(zhuǎn)變化了一個周期，頻率為1 388 Hz 左右.圖18 中展示了軸向Y/D＝3.3 截面上t＝0.034 26～0.035 086 s 的局部高溫區(qū)變化圖，黑色實線圈出了截面上尺寸比較大的高溫區(qū)，可以發(fā)現(xiàn)高溫區(qū)旋轉(zhuǎn)變化了一圈，頻率為1 210 Hz 左右.在本文的統(tǒng)計中，從t＝0.34 s 開始高溫區(qū)變化的第1 個周期如下：Y/D＝1 截面高溫區(qū)變化頻率在1 410 Hz 附近；Y/D＝2 截面高溫區(qū)變化頻率在1 280 Hz 附近；Y/D＝4 截面高溫區(qū)變化頻率在1 140 Hz 附近.截面上高溫區(qū)的變化頻率在統(tǒng)計的時間內(nèi)會有所差異，但總體范圍與PVC 的振蕩頻率相近，均在1 000～2 500 Hz，說明流場中的渦旋結(jié)構(gòu)對燃燒有影響，且PVC 會影響附近流體，使之出現(xiàn)進動特征，導致燃燒場出現(xiàn)不穩(wěn)定性，使垂直軸向截面的局部高溫區(qū)域呈現(xiàn)周期性旋轉(zhuǎn).

圖17 垂直軸向Y/D＝1.4截面上高溫區(qū)旋轉(zhuǎn)變化Fig.17 Rotation change of high temperature zone at section Y/D＝1.4 in vertical axis

圖18 垂直軸向Y/D＝3.3截面上高溫區(qū)旋轉(zhuǎn)變化Fig.18 Rotation change of high temperature zone at section Y/D＝3.3 in vertical axis

高溫區(qū)的旋轉(zhuǎn)變化現(xiàn)象隨著軸向距離的增加而越來越不明顯，這是由于隨著軸向距離的增加，渦量減小，附近流場的進動特征沿流向逐漸衰減，但變化頻率近似，說明這些位置的局部高溫區(qū)在空間上處在同一個大尺度流動結(jié)構(gòu)中.氣動力學引起的渦旋形成及破碎過程造成了燃燒室中燃燒的不穩(wěn)定，并且PVC 會通過多種方式影響旋流火焰的穩(wěn)定性，如PVC 會加強燃料和空氣的混合，強化未燃和已燃氣體的混合、卷曲、拉伸以及局部熄滅的化學反應(yīng)區(qū)，也有可能與火焰的熱聲耦合振蕩相互作用[16-17].

取剪切層內(nèi)點X＝0，Y＝0.021 m，Z＝0.021 m，即X＝0，Y/D＝1.4，Z/D＝1.4.圖19 展示了在統(tǒng)計的時間內(nèi)瞬時壓力脈動和放熱率脈動變化情況，在時域圖中可以看到其壓力脈動和放熱率脈動的波動具有一定的周期性，在頻域圖中，壓力脈動和放熱脈動對幅值貢獻比較大的頻率均在1 000～2 500 Hz，且在1 000～2 500 Hz 內(nèi)頻譜近似.燃燒場中的高溫區(qū)呈現(xiàn)徑向旋轉(zhuǎn)變化，存在變化頻率，并且壓力脈動和放熱率脈動的波動具有一定的周期性，說明在該情況下，燃燒出現(xiàn)不穩(wěn)定性現(xiàn)象.

圖19 壓力及放熱率脈動數(shù)據(jù)Fig.19 Fluctuation data of pressure and heat release rate

3 結(jié)論

本文使用LES-TPDF 方法對雙旋流模型燃燒室流動和燃燒狀態(tài)進行數(shù)值模擬，將溫度高于1 700 K的區(qū)域定義為高溫區(qū)，研究在軸向位置Y/D＝0～5的范圍內(nèi)燃燒室流場內(nèi)的局部高溫區(qū)演化現(xiàn)象，發(fā)現(xiàn)局部高溫區(qū)域與進動渦核位置接近.在本文研究情況中，主要結(jié)論如下：

(1) 在燃燒室頭部，燃燒主要發(fā)生在流場中的剪切層，在旋流流場的剪切層中存在進動渦核，這個不穩(wěn)定的大尺度旋轉(zhuǎn)流動結(jié)構(gòu)使周圍流體出現(xiàn)進動特征，加強已燃的高溫氣體和未燃的高溫氣體摻混，促進燃燒和放熱，影響高溫區(qū)的形成.

(2) 將局部高溫區(qū)域繞軸旋轉(zhuǎn)一圈的時間定義為一個周期，在垂直軸向Y/D＝0～5 不同截面上，高溫區(qū)的變化頻率均在1 000～2 500 Hz，與進動渦核頻率近似.說明雙旋流模型燃燒室流場中的渦旋結(jié)構(gòu)是高溫區(qū)早期形成的重要影響因素，導致燃燒場高溫區(qū)的不穩(wěn)定演化現(xiàn)象.