秦文瑾，周孟林，邵 宇

(上海理工大學(xué)機(jī)械工程學(xué)院，上海 200093)

航空發(fā)動(dòng)機(jī)科技的迅速發(fā)展影響著我國武器裝備、民航運(yùn)輸?shù)阮I(lǐng)域的發(fā)展趨勢.航空煤油是航空發(fā)動(dòng)機(jī)的主要燃料，它是一類組成復(fù)雜的碳?xì)浠衔铮渫轭惢衔锏奶荚哟蠖鄶?shù)為十至十六，主要含有單烷烴、環(huán)烷烴、芳香烴以及多環(huán)烷烴等不同類型的元素，涉及到的燃燒基本類型機(jī)理幾百種成分以及上千個(gè)基元反應(yīng)[1].航空煤油多組分替代燃料JW(Jet-A world average)[2]是由美國5 個(gè)不同廠家生產(chǎn)的Jet-A 航空煤油樣品混合而成，可以很好地表征Jet-A 航空煤油的物理化學(xué)性質(zhì)[3]，航空煤油噴霧霧化情況和燃燒效果有著非常重要的聯(lián)系.因此，了解多組分替代燃料JW 噴霧特性有利于后續(xù)提高燃油經(jīng)濟(jì)性和改善排放性能.一些學(xué)者針對Jet-A 航空煤油和航空煤油替代燃料開展了相應(yīng)研究.王維[4]針對航空煤油，在定容彈實(shí)驗(yàn)的基礎(chǔ)上，采用直接攝影法和陰影法研究了燃料的理化性質(zhì)對噴霧特性的影響.Yang等[5]通過對多組分Jet-A 綜合替代燃油的模擬結(jié)果與真實(shí)航空煤油的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行比較，討論了環(huán)境溫度和壓力，化學(xué)反應(yīng)機(jī)理以及液滴蒸發(fā)模型對航空煤油射流火焰結(jié)構(gòu)的影響.

近年來，高效科學(xué)的數(shù)學(xué)分析方法逐漸在數(shù)據(jù)分析中推廣應(yīng)用.由Rowley 等[6]及Schmid 等[7]提出的動(dòng)態(tài)模態(tài)分解方法(dynamic mode decomposition，DMD)主要應(yīng)用于湍流分析領(lǐng)域.該方法基于Koopman 分析理論，可以從固定時(shí)間間隔的瞬時(shí)場數(shù)據(jù)中提取重要?jiǎng)討B(tài)結(jié)構(gòu)，每個(gè)DMD 模態(tài)都代表特定的頻率信息[8].例如，王廣[9]基于DMD 方法，對具有特定幾何參數(shù)的擴(kuò)壓器內(nèi)部流場進(jìn)行解構(gòu)分析，捕捉到重要的渦結(jié)構(gòu).Mezi? 等[10]對一臺6 缸柴油機(jī)燃燒進(jìn)行了數(shù)值模擬，采用動(dòng)態(tài)模態(tài)分析的方法對其進(jìn)氣歧管流場進(jìn)行研究，提取了其特定頻率的流場結(jié)構(gòu)用于流場特性的分析.謝海潤等[11]將DMD 應(yīng)用于圓柱/葉片附近壓力場，以研究尾渦激振現(xiàn)象.Peng等[12]基于DMD 研究了氫燃料超燃沖壓發(fā)動(dòng)機(jī)燃燒室展向位置的頻率特性，通過表征主頻和增長因子研究燃燒不穩(wěn)定性.

本文采用大渦模擬方法對航空煤油多組分替代燃料JW 的射流噴霧現(xiàn)象進(jìn)行數(shù)值計(jì)算，研究其特性.并將DMD 方法應(yīng)用于該燃料射流噴霧濃度場特性分析，通過對瞬態(tài)濃度場數(shù)據(jù)庫進(jìn)行加工處理，獲得了不同頻率下的燃料濃度場結(jié)構(gòu)信息，豐富了人們對航空煤油多組分替代燃料JW 濃度場特性的認(rèn)識.

1 研究對象

美國阿貢實(shí)驗(yàn)室構(gòu)建了邊長為108 mm 的立方體光學(xué)定容燃燒室，其設(shè)有4 個(gè)光學(xué)觀察通道，在其內(nèi)進(jìn)行了多種工況下JW 無氧環(huán)境噴霧實(shí)驗(yàn)[13]，定容燃燒室如圖1 所示.本文在該光學(xué)定容燃燒室的基礎(chǔ)上，構(gòu)建了一個(gè)定容燃燒室，其尺寸為上底直徑40 mm、下底直徑80 mm、高100 mm，如圖2 所示.基于CONVERGE 數(shù)值模擬軟件，對JW 噴霧現(xiàn)象進(jìn)行了數(shù)值模擬.為保證計(jì)算精度和計(jì)算效率，網(wǎng)格劃分采用自適應(yīng)加密和局部加密兩種加密策略，網(wǎng)格數(shù)量為400 萬左右，如圖3 所示.初始邊界條件參照實(shí)驗(yàn)參數(shù)進(jìn)行設(shè)置.采用0.09 mm 直徑的單孔噴嘴，噴射壓力為 150 MPa，模擬噴射脈寬時(shí)間為1.5 ms，燃油溫度為363 K.采用KH-RT[14]液滴破碎模型和Frossling[15]蒸發(fā)模型.基于JW 航空煤油中存在的主要烴基，采用Chen 等[16]開發(fā)的4 種組分替代燃料模型，該替代燃料的組成：正十二烷NC12H26、異十六烷IC16H34、反式十氫萘C10H18、甲苯C7H8，其摩爾分?jǐn)?shù)分別為0.3、0.36、0.246、0.094.最后采用PISO 算法求解壓力-速度耦合問題.

圖1 ECN光學(xué)定容燃燒室Fig.1 ECN optical constant volume combustion chamber

圖2 定容燃燒室Fig.2 Constant volume combustion chamber

圖3 網(wǎng)格示意Fig.3 Schematic diagram of grid

2 計(jì)算模型的建立

2.1 動(dòng)態(tài)模態(tài)分解

DMD 算法也是一個(gè)可以進(jìn)行Koopman 算譜分析的數(shù)值驅(qū)動(dòng)方式，它可以通過虛擬無限維的結(jié)構(gòu)線性系統(tǒng)描述為實(shí)際有限維的本構(gòu)非線性體系.經(jīng)DMD 分析后獲得的每一條模態(tài)都對應(yīng)一個(gè)特征頻率，因此可以直接獲取單頻的動(dòng)態(tài)模態(tài)及其所對應(yīng)的原始場結(jié)構(gòu)，而不需再通過其他頻率分析來獲得模態(tài)的特征頻率.

DMD 算法需要擁有恒定時(shí)間間隔的數(shù)據(jù)集輸入，首先將瞬時(shí)物理場表示為矩陣形式：

式中：v 為任意一個(gè)物理量，vi代表第i 個(gè)物理場的采樣數(shù)據(jù)構(gòu)成的列向量.為由第1 個(gè)到第N-1 個(gè)物理場數(shù)據(jù)構(gòu)成的向量集合；為由第2 個(gè)到第N個(gè)物理場數(shù)據(jù)構(gòu)成的向量集合.

為了使用非線性系統(tǒng)的線性逼近，需要通過一個(gè)假設(shè)的矩陣A 來實(shí)現(xiàn)連個(gè)物理場之間的投影關(guān)聯(lián)，即：

此時(shí)矩陣的特征值、特征向量和能量放大因子等屬性則直接反映了物理場序列集合的系統(tǒng)特性.對矩陣進(jìn)行奇異值分解：

式中：矩陣U、W 是酉矩陣；∑是奇異值構(gòu)成的對角矩陣.可獲得系統(tǒng)矩陣A 的低維近似系統(tǒng)矩陣S 的表達(dá)形式：

此時(shí)而DMD 模態(tài)則可以表示為

通過上述操作，原瞬態(tài)物理場可被分解為

式中：αi為DMD 模態(tài)的幅值，代表對應(yīng)物理場的能量或強(qiáng)度.

使用DMD 方法分解數(shù)據(jù)得到的特征值是復(fù)數(shù)，實(shí)部和虛部則分別反映了不同的信息，實(shí)部和虛部分別對應(yīng)著物理場的能量強(qiáng)度和頻率.同時(shí)把所有特性值按照數(shù)值大小進(jìn)行繪圖處理，可以得出，位于單位圓之內(nèi)的特征值，代表了其物理結(jié)構(gòu)對應(yīng)的能量隨時(shí)間愈來愈小，位于單位圓之上的特征值，代表了其物理結(jié)構(gòu)的能量不隨時(shí)間變化，是穩(wěn)定的，而位于單位圓之外的特征值，其能量則會隨時(shí)間愈來愈大.

2.2 大渦模擬湍流模型

大渦模擬方法(LES)主要采用了空間過濾方法，對于較大尺寸的渦團(tuán)則采用了直接數(shù)值模擬的方法，而較小尺寸的渦團(tuán)則主要是嵌入模擬的方法，所以大渦模擬可以對湍流流場進(jìn)行更加有效的模擬，可以更加簡單準(zhǔn)確地捕捉大渦團(tuán)、小渦團(tuán).另外，大渦模擬方法在計(jì)算效率和計(jì)算準(zhǔn)確性兩方面也是較其他方式有更好的效果.

本文基于大渦模擬方法，采用歐拉-拉格朗日方法進(jìn)行燃料噴霧的數(shù)值模擬，采用基于動(dòng)態(tài)結(jié)構(gòu)的湍流模型[17]，該模擬將亞網(wǎng)格應(yīng)力直接描述為亞網(wǎng)格湍動(dòng)能函數(shù)，并通過計(jì)算亞網(wǎng)格湍動(dòng)能輸運(yùn)方程對模型進(jìn)行封閉處理.并利用對N-S 方程的濾波處理，獲得了大渦模擬的控制方程.

式中：τij亞網(wǎng)格應(yīng)力，亞網(wǎng)格應(yīng)力是一個(gè)未知量，基于流場各向同性渦黏假設(shè)，亞網(wǎng)格應(yīng)力可根據(jù)大尺度流場的應(yīng)變率確定，即

式中：Sij是大尺度變形率張量；νt是亞網(wǎng)格渦黏系數(shù).亞網(wǎng)格渦黏系數(shù)個(gè)湍動(dòng)能輸運(yùn)方程分別為：

其中：k 是亞網(wǎng)格湍動(dòng)能；常數(shù) Ck和Cε需要根據(jù)實(shí)際的瞬態(tài)流動(dòng)問題來設(shè)定.

2.3 液滴破碎模型

液滴破碎模型選用KH-RT 模型，該基于氣液界面擾動(dòng)波的不穩(wěn)定理論，包括KH 和RT 兩個(gè)子模型，其中KH 模型主導(dǎo)一次破碎過程，而RT 模型則與KH 模型結(jié)合共同主導(dǎo)二次破碎過程.當(dāng)液體射流從噴嘴噴出，氣液交界面的剪切作用形成的KH 擾動(dòng)波不斷增長，使得燃油液柱開始剝落形成液滴.隨后，燃油液滴與環(huán)境氣體間的大密度梯度形成RT 擾動(dòng)波，在KH 和RT 擾動(dòng)波的作用下，燃油液滴開始進(jìn)一步的破碎.

KH 破碎模型中，液滴破碎由油束表面產(chǎn)生的波的最大生產(chǎn)率ΩΚΗ及其最大波長ΛKH控制.液滴破碎時(shí)間τKH表示為

式中：1B 表示KH 破碎特征時(shí)間系數(shù)；r0表示為破碎前的液滴半徑.破碎開始后的液滴半徑為

式中：0B 表示尺寸系數(shù)，該值直接影響破碎后的液滴半徑.RT 破碎模型建立在氣液密度梯度導(dǎo)致的慣性力基礎(chǔ)上，液滴破碎由慣性力產(chǎn)生的波的最大波動(dòng)頻率ΩRT及其對應(yīng)波長ΛRT控制.液滴破碎時(shí)間尺度τRT表示為

式中：Cτ表示RT 破碎特征時(shí)間系數(shù).當(dāng)波動(dòng)時(shí)間大于τRT時(shí)，二次破碎再次發(fā)生，液滴變得更小，直至到達(dá)極限.此時(shí)，液滴極限半徑rc表示為

式中，cRT表示尺寸系數(shù)，該值同樣直接影響破碎后的液滴半徑.除了上述的特征時(shí)間系數(shù)和尺寸系數(shù)外，CONVERGE 中，KH 和RT 模型中還分別增加了速度系數(shù)1C 和長度系數(shù) Cbl，用于調(diào)整液滴破碎，從而優(yōu)化噴霧效果.具體參數(shù)的選擇如表1 所示.

表1 KH-RT模型參數(shù)Tab.1 The partmeters of KH-RT model

3 結(jié)果與討論

3.1 網(wǎng)格無關(guān)性驗(yàn)證

圖4 是JW 燃料在氣體環(huán)境溫度為1 000 K、壓力為 6.03 MPa 工況下，分別選取了最小尺寸為0.2 mm、0.1 mm、0.062 5 mm 的3 套網(wǎng)格尺寸進(jìn)行數(shù)值模擬得到的噴霧貫穿距與實(shí)驗(yàn)測量的對比.液相貫穿距和氣相貫穿距分別定義為噴嘴軸線上燃油液體體積分?jǐn)?shù)、氣相質(zhì)量分?jǐn)?shù)為0.1%位置到達(dá)噴嘴的軸向距離.從圖中可知，計(jì)算得到貫穿距總體趨勢與實(shí)驗(yàn)符合較好.噴霧前期，燃油從噴嘴噴出，其初始速度很大.由于受到環(huán)境氣體的阻力，導(dǎo)致燃油液滴速度不斷減小.液相貫穿距在噴射開始后，存在約0.05 ms 的快速增長期，隨即達(dá)到燃料噴射與蒸發(fā)的平衡狀態(tài)，液相貫穿距在極小的范圍內(nèi)波動(dòng).氣相貫穿距則先快速增大，然后增速放緩.與實(shí)驗(yàn)結(jié)果相比，3 種網(wǎng)格尺寸下液體貫穿距較為相似，但0.2 mm的網(wǎng)格下，氣相貫穿距誤差較大，0.1 mm 和0.062 5 mm 的計(jì)算結(jié)果較為接近，且與實(shí)驗(yàn)結(jié)果誤差相對較小.為了節(jié)省計(jì)算資源，本研究采用0.1 mm的網(wǎng)格進(jìn)行后文數(shù)值模擬工作.

圖4 JW實(shí)驗(yàn)與不同網(wǎng)格尺寸模擬的貫穿距Fig.4 Penetration of JW experiment and simulation with different mesh sizes

3.2 燃油液滴與蒸氣質(zhì)量分?jǐn)?shù)分布

圖5 進(jìn)一步給出了計(jì)算得到的JW 燃油液滴和蒸氣質(zhì)量分?jǐn)?shù)的空間分布，與實(shí)驗(yàn)中關(guān)于JW 燃料主要組分正十二烷燃油液滴和蒸氣質(zhì)量分?jǐn)?shù)的空間分布對比.因?yàn)镴W 射流的噴射與空氣間產(chǎn)生了強(qiáng)烈的剪切運(yùn)動(dòng)，導(dǎo)致JW 氣相濃度場產(chǎn)生了非常顯著的渦團(tuán)結(jié)構(gòu)，數(shù)值模擬結(jié)果很好地呈現(xiàn)了JW 射流形態(tài)的高瞬態(tài)非對稱性特點(diǎn).圖6 是射流基本趨于穩(wěn)定情況下，JW 射流軸向不同位置處燃油蒸氣混合分?jǐn)?shù)分布情況.從整體來看，在軸向距離噴嘴20 mm 處的混合分?jǐn)?shù)峰值大于30 mm 處的混合分?jǐn)?shù)，且30 mm 處JW 的的徑向分布寬度要比20 mm 處大，即燃料在噴射區(qū)間下游擴(kuò)散更廣，總體看來燃油蒸氣沿軸向呈傘狀向四周擴(kuò)散.另外由于湍流擾動(dòng)，混合分?jǐn)?shù)和溫度曲線都呈現(xiàn)空間波動(dòng)特點(diǎn).

圖5 正十二烷與JW的燃油蒸氣質(zhì)量分?jǐn)?shù)分布和燃油液滴粒徑分布Fig.5 The distribution of fuel mass fraction and droplet radius of n-dodecane and JW

圖6 平均混合分?jǐn)?shù)分布Fig.6 Mean mixing fraction distribution

3.3 濃度場的DMD分析

在空間中，液體燃料噴霧瞬態(tài)濃度場為在整體性較大尺寸的宏觀形態(tài)結(jié)構(gòu)的基礎(chǔ)上，疊加著數(shù)量較多的局部尺度較小的形態(tài)結(jié)構(gòu)，并且不同尺度的噴霧結(jié)構(gòu)分別攜帶了不同的強(qiáng)度和頻率.僅通過對瞬態(tài)濃度場觀察，難以對液體燃料噴霧特性進(jìn)行深入理解.由于噴霧屬于液相，其發(fā)生破碎蒸發(fā)后形成氣相的濃度場，即氣相濃度場結(jié)構(gòu)很大程度上受液相噴霧場的影響控制，故也可認(rèn)為氣相濃度場結(jié)構(gòu)在一定程度上是對液體燃料噴霧宏觀形態(tài)的拓展.利用DMD方法對不同頻率的濃度場結(jié)構(gòu)進(jìn)行提取分析，可以預(yù)測濃度場結(jié)構(gòu)和強(qiáng)度的動(dòng)態(tài)變化，豐富對燃料噴霧的認(rèn)識，為后續(xù)對噴霧濃度場的人工干預(yù)奠定基礎(chǔ).

圖7 是使用DMD 方法對濃度場數(shù)據(jù)庫進(jìn)行后處理，可以得出特征值是呈現(xiàn)共軛分布的趨勢，實(shí)部和虛部分別對應(yīng)著濃度場的能量強(qiáng)度和運(yùn)動(dòng)頻率.同時(shí)把所有特性值按照數(shù)值大小進(jìn)行繪圖處理，可以得出，位于單位圓之內(nèi)的特征值，代表了濃度場結(jié)構(gòu)對應(yīng)的強(qiáng)度隨時(shí)間愈來愈小，位于單位圓之上的特征值，代表了其濃度場結(jié)構(gòu)的強(qiáng)度隨時(shí)間保持不變，而位于單位圓之外的特征值，其強(qiáng)度則會隨時(shí)間愈來愈大.為了便于研究分析，根據(jù)相對于單位圓位置的不同，在圖7 中選取了A、B、C、D 4 個(gè)不同頻率的代表點(diǎn)，且其頻率逐漸增加，特征值大小如表2示，DMD 模態(tài)空間分布如圖8 所示.從圖7 中可以看出，點(diǎn)A 幾乎落在單位圓上，其對應(yīng)的濃度場強(qiáng)度隨著時(shí)間的變化保持不變.且A 點(diǎn)的特征值原本是一個(gè)復(fù)數(shù)，但是其虛部為0，所以該點(diǎn)的DMD 模態(tài)也只有實(shí)部，沒有虛部.同時(shí)由于B 點(diǎn)、C 點(diǎn)和D 點(diǎn)的特征值既有實(shí)部和虛部，其對應(yīng)的模態(tài)也是有實(shí)部和虛部的.如圖8 所示，B 點(diǎn)、C 點(diǎn)和D 點(diǎn)對應(yīng)的模態(tài)質(zhì)量分?jǐn)?shù)較小的區(qū)域逐漸增多，其對應(yīng)的特征值虛部也逐漸增大，即代表其攜帶的濃度場頻率逐漸變大.從圖中可知，B 點(diǎn)和C 點(diǎn)特征值位于單位圓之內(nèi)，且經(jīng)過對數(shù)化處理后，實(shí)部特征值為負(fù)數(shù)，表示其模態(tài)對應(yīng)的濃度場結(jié)構(gòu)強(qiáng)度在時(shí)間域上愈來愈弱.而D 點(diǎn)則位于單位圓之外，特征值經(jīng)過對數(shù)化后實(shí)部為正數(shù)，表示其模態(tài)對應(yīng)的濃度場結(jié)構(gòu)強(qiáng)度在時(shí)間域上愈來愈強(qiáng).

表2 代表點(diǎn)特征值Tab.2 Eigenvalues of selected points

圖7 特征值分布Fig.7 Eigenvalues distribution

圖8 代表點(diǎn)DMD模態(tài)Fig.8 Represents the point DMD mode

圖9 為DMD 模態(tài)的幅值與其對應(yīng)的特征值的關(guān)系，縱坐標(biāo)為幅值，表示每一個(gè)模態(tài)濃度場結(jié)構(gòu)的強(qiáng)度大小，橫坐標(biāo)分別為特征值實(shí)部和虛部的對數(shù)，分別代表了模態(tài)的耗散率和頻率，其中A、B、C、D 4個(gè)點(diǎn)如圖中標(biāo)記所示.其中A 點(diǎn)模態(tài)的幅值遠(yuǎn)大于其他模態(tài)的幅值，表示該部分濃度場結(jié)構(gòu)強(qiáng)度最大，B 點(diǎn)、C 點(diǎn)和D 點(diǎn)對應(yīng)模態(tài)的幅值比較小，表示這3個(gè)模態(tài)對應(yīng)濃度場結(jié)構(gòu)強(qiáng)度較弱.同時(shí)A 點(diǎn)模態(tài)的對應(yīng)的頻率為零，說明其對不同時(shí)刻的濃度場結(jié)構(gòu)的貢獻(xiàn)是恒定不變的，而其他3 個(gè)模態(tài)對應(yīng)的頻率都不為零，表示它們對不同時(shí)刻的濃度場結(jié)構(gòu)的貢獻(xiàn)是不同的，即構(gòu)成了濃度場的差異性.圖9(b)顯示A 點(diǎn)的特征值實(shí)部為零，即耗散率為零，表示該模態(tài)對應(yīng)的濃度場結(jié)構(gòu)的強(qiáng)度隨時(shí)間保持不變.由于B 點(diǎn)和C 點(diǎn)特征值實(shí)部為負(fù)數(shù)，表示其對應(yīng)的濃度場結(jié)構(gòu)的強(qiáng)度隨時(shí)間愈來愈弱，而D 點(diǎn)特征值實(shí)部為正數(shù)，表示該濃度場結(jié)構(gòu)的強(qiáng)度隨時(shí)間愈來愈強(qiáng).

圖9 DMD模態(tài)幅值與頻率和耗散率對應(yīng)的關(guān)系Fig.9 The relation between DMD mode amplitude and frequency and dissipation rate

4 結(jié)論

本文針對航空煤油多組分替代燃料JW 在定容燃燒室內(nèi)的噴霧過程，使用大渦模擬方法進(jìn)行了冷態(tài)噴霧數(shù)值模擬，探討了多組分替代燃料JW 的噴霧特性.利用動(dòng)態(tài)模態(tài)分解方法對定容燃燒室內(nèi)濃度場數(shù)據(jù)進(jìn)行了深度加工，預(yù)測了不同頻率濃度場結(jié)構(gòu)強(qiáng)度在噴霧過程中的變化情況，挖掘了濃度場噴霧結(jié)構(gòu)更深層次的信息.主要結(jié)論如下：

(1)使用大渦模擬方法數(shù)值模擬得到貫穿距總體趨勢與實(shí)驗(yàn)結(jié)果符合較好，較準(zhǔn)確地捕捉氣相和液相噴霧發(fā)展過程，并很好地預(yù)估了JW 噴霧的形態(tài)特征，反映了JW 燃料混合分?jǐn)?shù)分布高瞬態(tài)非對稱性.在射流基本趨于穩(wěn)定情況下，JW 燃料在噴射區(qū)間下游擴(kuò)散更廣.另外由于湍流擾動(dòng)，混合分?jǐn)?shù)呈現(xiàn)空間波動(dòng)的特點(diǎn).

(2)由于瞬態(tài)噴霧濃度場結(jié)構(gòu)是由多個(gè)不同頻率的結(jié)構(gòu)疊加構(gòu)成，但是不能直接在瞬態(tài)濃度場中體現(xiàn)，利用DMD 方法可對不同頻率的JW 噴霧蒸氣濃度場結(jié)構(gòu)進(jìn)行提取.對提取的不同頻率濃度場結(jié)構(gòu)進(jìn)行動(dòng)態(tài)發(fā)展分析，可以預(yù)測濃度場結(jié)構(gòu)的強(qiáng)度的變化，進(jìn)而豐富人們對JW 燃料噴霧濃度場結(jié)構(gòu)的認(rèn)識，為下一步對濃度場結(jié)構(gòu)的人工干預(yù)奠定基礎(chǔ)，進(jìn)而便于形成更加合理的油氣混合物，實(shí)現(xiàn)高效燃燒.