分區(qū)參數(shù)對超聲速湍流燃燒動(dòng)態(tài)分區(qū)火焰面模擬的影響1)

南江浪 張 政 姚 衛(wèi) ,3) 劉鳳君

* (北京動(dòng)力機(jī)械研究所,北京 100074)

? (中國科學(xué)院力學(xué)研究所高溫氣體動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室,北京 100190)

** (中國科學(xué)院大學(xué)工程科學(xué)學(xué)院,北京 100049)

引言

近幾十年來,超燃沖壓發(fā)動(dòng)機(jī)作為高超聲速飛行器最有前景的發(fā)動(dòng)機(jī)之一,是世界各國競相發(fā)展的熱點(diǎn)領(lǐng)域之一[1].超燃沖壓發(fā)動(dòng)機(jī)具有結(jié)構(gòu)簡單、性能突出等特點(diǎn),且發(fā)動(dòng)機(jī)能夠提供較高比沖,滿足高超聲速飛行的性能要求,因此超燃沖壓發(fā)動(dòng)機(jī)在遠(yuǎn)程快速的民用運(yùn)輸和可重復(fù)的太空探索往返飛行方面表現(xiàn)出了巨大的潛力[2].一般來說,超燃沖壓發(fā)動(dòng)機(jī)地面實(shí)驗(yàn)具有極大挑戰(zhàn),特別在高馬赫飛行條件下更加難以實(shí)現(xiàn)復(fù)雜且長時(shí)間的測試.此外,受限于目前的測量技術(shù)發(fā)展水平和高昂的地面實(shí)驗(yàn)費(fèi)用,超燃沖壓發(fā)動(dòng)機(jī)的設(shè)計(jì)階段進(jìn)行飛行實(shí)驗(yàn)是難以實(shí)現(xiàn)的.因此,高保真計(jì)算流體力學(xué)(CFD)數(shù)值模擬已經(jīng)成為超燃沖壓發(fā)動(dòng)機(jī)研究和設(shè)計(jì)領(lǐng)域的必要方法[3].

由于超聲速燃燒流動(dòng)過程的復(fù)雜性和高昂的計(jì)算成本,模擬超聲速燃燒通常具有較大難度[4].在超聲速燃燒室內(nèi)強(qiáng)沖擊波與其他流動(dòng)模式,如邊界層[5]、剪切層[6]和渦流[7]進(jìn)行復(fù)雜相互作用,這些過程與燃燒化學(xué)反應(yīng)耦合形成了一個(gè)復(fù)雜的可壓縮湍流燃燒過程,并在燃燒室內(nèi)隨時(shí)空迅速演變[8].因此,超燃沖壓發(fā)動(dòng)機(jī)中的燃燒是極不穩(wěn)定的,并且包含豐富的特征時(shí)間尺度和空間尺度,這給數(shù)值模擬帶來了巨大的困難[4].為準(zhǔn)確捕捉發(fā)動(dòng)機(jī)內(nèi)復(fù)雜的極限燃燒行為,如點(diǎn)熄火,超亞燃模態(tài)轉(zhuǎn)換等,采用詳細(xì)化學(xué)反應(yīng)與大渦模擬(LES)成為超燃沖壓發(fā)動(dòng)機(jī)數(shù)值仿真的必然趨勢[9].

然而,在實(shí)際工程應(yīng)用中采用包含千萬或億級網(wǎng)格量的高分辨率LES 需花費(fèi)巨大的計(jì)算時(shí)間成本,尤其對于包含幾十甚至幾百步基元反應(yīng)的燃燒模擬,其化學(xué)反應(yīng)求解將占據(jù)50%以上的CPU 時(shí)間[10].為了減少化學(xué)反應(yīng)求解的計(jì)算時(shí)間,目前主要通過,加速化學(xué)基元反應(yīng)直接求解,如化學(xué)反應(yīng)機(jī)理簡化,化學(xué)反應(yīng)建表,及基于人工神經(jīng)網(wǎng)絡(luò)或GPU加速等方法,或采用守恒標(biāo)量類湍流燃燒模型,如各類火焰面模型[11]、條件矩封閉模型(CMC)[12]等.上述模型將流動(dòng)與燃燒解耦,將反應(yīng)標(biāo)量通過分布函數(shù)映射于守恒標(biāo)量(如混合分?jǐn)?shù)、化學(xué)反應(yīng)進(jìn)度變量)的狀態(tài)空間,實(shí)現(xiàn)湍流燃燒的降維求解,無需采用有限速率類模型(如PaSR[13]和 EDC[14]等)中的方法對每個(gè)計(jì)算網(wǎng)格節(jié)點(diǎn)進(jìn)行耗時(shí)的化學(xué)反應(yīng)求解,從而實(shí)現(xiàn)包含詳細(xì)化學(xué)反應(yīng)機(jī)理的湍流燃燒高效模擬.

相較于亞聲速燃燒,高超聲速燃燒多為高雷諾數(shù)流動(dòng),化學(xué)反應(yīng)與湍流脈動(dòng)時(shí)間尺度相當(dāng),湍流小尺度渦與反應(yīng)區(qū)相互作用.因此,基于薄火焰面系綜的火焰面類模型在高超聲速湍流燃燒的適應(yīng)性受到挑戰(zhàn).而CMC 類模型基于統(tǒng)計(jì)上的條件平均概念構(gòu)建映射于條件參數(shù)空間內(nèi)的火焰面(條件矩),不受物理上薄火焰面假設(shè)的限制,具有更廣泛的適用性[12].為準(zhǔn)確封閉CMC 類模型反應(yīng)源項(xiàng),需滿足CMC 網(wǎng)格反應(yīng)參數(shù)在其條件均值附近小脈動(dòng)的要求.對于低速不可壓縮擴(kuò)散火焰,簡單的單條件一階CMC 模型即可獲得與實(shí)驗(yàn)較好吻合的結(jié)果[15].對于包含復(fù)雜的點(diǎn)火、熄火等臨界燃燒過程,如超聲速湍流燃燒問題,CMC 模型的準(zhǔn)確封閉可采用雙條件矩方法,以降低反應(yīng)參數(shù)偏離條件均值的脈動(dòng),或通過二階CMC 模型對脈動(dòng)項(xiàng)進(jìn)行封閉[16].雙條件矩及二階矩方法雖然提升了模型的準(zhǔn)確性,但其一方面增加了求解方程的數(shù)量,另一方面對于高階矩的?；忾]仍未成熟.此外,對于包含詳細(xì)反應(yīng)機(jī)理的LES 模擬,采用雙條件或二階矩將顯著提高計(jì)算花銷和模型復(fù)雜程度,不利于模型的工程化應(yīng)用.近年來,基于動(dòng)態(tài)分區(qū)概念,Yao 等[17]采用與傳統(tǒng)CMC 模型中靜態(tài)幾何坐標(biāo)分區(qū)方法不同的基于熱力學(xué)參數(shù)動(dòng)態(tài)分區(qū)的CMC 網(wǎng)格聚類劃分方法,提出動(dòng)態(tài)分區(qū)火焰面模型(DZFM).通過混合分?jǐn)?shù)以及其他與反應(yīng)狀態(tài)強(qiáng)相關(guān)的參數(shù)對數(shù)值計(jì)算網(wǎng)格進(jìn)行聚類分區(qū),使各分區(qū)內(nèi)反應(yīng)狀態(tài)參數(shù)的條件平均脈動(dòng)降低,在局部分區(qū)內(nèi)實(shí)現(xiàn)條件參數(shù)(混合分?jǐn)?shù))與反應(yīng)參數(shù)的強(qiáng)相關(guān),從而實(shí)現(xiàn)分區(qū)火焰面對局部反應(yīng)狀態(tài)的準(zhǔn)確表征.當(dāng)前基于DZFM 方法的超聲速湍流燃燒模擬普遍采用混合分?jǐn)?shù)與流向坐標(biāo)雙參數(shù)分區(qū)的方法,分區(qū)參數(shù)空間的選取對DZFM 模型準(zhǔn)確性的影響尚未進(jìn)行研究.本文采用德國航空航天中心(DLR)研發(fā)的超燃沖壓發(fā)動(dòng)機(jī)模型進(jìn)行數(shù)值模型測試研究[18].其結(jié)構(gòu)簡單,實(shí)驗(yàn)數(shù)據(jù)豐富,國內(nèi)外學(xué)者已采用DLR燃燒室開展了大量的超聲速湍流燃燒模型開發(fā)與驗(yàn)證研究[19].本文基于DZFM,對DLR 支板穩(wěn)焰的吸氣式超聲速氫氣燃燒室進(jìn)行改進(jìn)延遲分離渦模擬(IDDES),研究不同分區(qū)參數(shù)的選取對DZFM 模型準(zhǔn)確性的影響.

1 物理模型與數(shù)值方法

1.1 控制方程與湍流模型

對包含反應(yīng)的三維可壓縮N-S (Navier-Stokes)方程進(jìn)行求解.所有的變量(ρ,ui,Ht,ξ) 通過Favre平均分解為大渦模擬空間過濾后直接求解的平均量和未求解的脈動(dòng)量f′,由此得到的質(zhì)量、動(dòng)量、能量、混合分?jǐn)?shù)及其脈動(dòng)的控制方程如下

式中t為時(shí)間,xi和ui分別表示i方向的坐標(biāo)和速度分量,表示平均速度,表示平均壓力,表示總焓即絕對焓、動(dòng)能與湍動(dòng)能之和.為混合分?jǐn)?shù),為混合分?jǐn)?shù)脈動(dòng)均方值,DT表示混合物的熱擴(kuò)散系數(shù),Dξ代表混合平均質(zhì)量擴(kuò)散系數(shù),Dα表示組分 α 在混合物中的質(zhì)量擴(kuò)散系數(shù).R=Ru/W表示基于氣體混合物摩爾質(zhì)量,W為氣體常數(shù),其中通用氣體常數(shù),Ru=8.314 J/(mol·K) .注意,這里并不求解通常意義上的Favre 平均質(zhì)量分?jǐn)?shù)的輸運(yùn)方程,而是求解條件組分的輸運(yùn)方程Qα,二者的關(guān)系為

式中,η 為混合分?jǐn)?shù)空間內(nèi)的取樣變量,P(η) 為依賴于的 β 函數(shù)形式的概率密度函數(shù).

亞格子尺度(SGS) 能量和組分?jǐn)U散項(xiàng)和SGS黏性耗散項(xiàng)為小量可忽略[20].在式(2)～式(4)中剩余的與亞格子尺度相關(guān)的項(xiàng)有雷諾應(yīng)力 (τij)、SGS能量通量 (ΨT,j)以及混合分?jǐn)?shù)通量 (Ψξ,j),需要采用相應(yīng)的模型將上述SGS 項(xiàng)封閉.其中雷諾應(yīng)力項(xiàng)采用Boussinesq 渦黏性假設(shè)模擬,即假設(shè)其正比于可解尺度的應(yīng)變率張量

式中 νsgs為SGS 湍流黏性,ksgs為SGS 湍動(dòng)能.對于待封閉的亞格子混合分?jǐn)?shù)及能量通量,可采用濃度梯度假設(shè)模化,則有以 及Ψξ,j=式中湍流普朗特?cái)?shù)Prt=0.89 以及湍流施密特?cái)?shù)Sct=0.5[21].考慮到直接求解壁面邊界層所需網(wǎng)格尺度對于大尺度的發(fā)動(dòng)機(jī)模擬而言計(jì)算花銷巨大,采用RANS/LES 混合模型,即改進(jìn)延遲分離渦模擬(IDDES)[22]方法降低壁面邊界層求解計(jì)算量.其對壁面邊界層的黏性子層采用k-ω 剪切應(yīng)力輸運(yùn)雷諾時(shí)均模型[23]求解,而外流區(qū)域采用大渦模擬方法求解.

控制方程中的氣體混合物的熱力學(xué)及輸運(yùn)參數(shù)使用CHEMKIN-II 程序計(jì)算[24],其采用NISTJANAF 熱物性數(shù)據(jù)庫[25]以及CHEMKIN 格式的輸運(yùn)數(shù)據(jù)庫.黏性系數(shù)、比熱以及導(dǎo)熱率等均僅為溫度的函數(shù).混合物的平均黏性系數(shù)以及導(dǎo)熱率分別采用修正的Wilke 定律[26]以及摩爾平均計(jì)算得到.

1.2 湍流燃燒模型

采用動(dòng)態(tài)分區(qū)火焰面模型DZFM,通過基于混合分?jǐn)?shù)空間的分區(qū)條件平均組分輸運(yùn)方程及相應(yīng)的條件概率密度分布將化學(xué)反應(yīng)與湍流流動(dòng)進(jìn)行解耦求解,具體過程如下.

計(jì)算域網(wǎng)格通過特定的參數(shù)進(jìn)行動(dòng)態(tài)的聚合分區(qū),求解每個(gè)分區(qū)內(nèi)條件平均組分濃度(Qα=〈Yα|η=ξ(x,t),x∈zone〉)的輸運(yùn)而非基于計(jì)算網(wǎng)格求解式(4)中的平均組分濃度.條件平均組分濃度與瞬態(tài)組分濃度間的關(guān)系可以表示為Yα(x,t)=Qα(η=,其中表示瞬態(tài)值相對條件平均值的脈動(dòng).通過多參數(shù)分區(qū)和分區(qū)逐步細(xì)化,可以認(rèn)為在分區(qū)內(nèi)滿足,即分區(qū)內(nèi)條件反應(yīng)狀態(tài)的統(tǒng)計(jì)均一性假設(shè).由此可實(shí)現(xiàn)對組分輸運(yùn)非線性源項(xiàng)的準(zhǔn)確求解,以及對的一階或高階項(xiàng)進(jìn)行合理簡化.將Qα代入組分輸運(yùn)方程并結(jié)合相應(yīng)的高階脈動(dòng)項(xiàng)簡化及高雷諾數(shù)近似后得到Qα的輸運(yùn)方程如下[27]

條件平均溫度QT通過對分區(qū)內(nèi)的歷史統(tǒng)計(jì)方法計(jì)算得到[28-29],而非單獨(dú)求解一個(gè)包含一系列未封閉項(xiàng)的條件焓輸運(yùn)方程,從而可節(jié)約計(jì)算時(shí)間.此處的平均溫度QT僅控制條件平均反應(yīng)速率,流場的溫度則通過計(jì)算網(wǎng)格中的總焓輸運(yùn)方程(式(3))得到.由式(13)求解得到分區(qū)火焰面的Qα后,可通過對其在計(jì)算網(wǎng)格上的積分得到無條件的組分質(zhì)量分?jǐn)?shù).在此,假定分區(qū)內(nèi)的混合分?jǐn)?shù)的概率密度分布函數(shù)為 β 函數(shù),,其中由方程式(4)和式(5)求解得到.

湍流流動(dòng)和燃燒的數(shù)值求解,基于開源平臺OpenFOAM[30]開發(fā)的可壓縮湍流反應(yīng)求解器Amber進(jìn)行.對流項(xiàng)采用二階低馬赫數(shù)修正的混合KNP/中心格式[31-33],其可準(zhǔn)確捕捉遠(yuǎn)離激波處的湍流精細(xì)結(jié)構(gòu)同時(shí)保持激波間斷附近的穩(wěn)定性.面通量的插值采用總變差不增(TVD)型的限制器Minmod.時(shí)間推進(jìn)采用二階精度的Crank-Nicholson 格式.該求解器已在諸多超聲速湍流燃燒模擬中得到驗(yàn)證[10,34-35].

1.3 測試工況及計(jì)算設(shè)置

DLR 中心支板超燃沖壓發(fā)動(dòng)機(jī)幾何結(jié)構(gòu)如圖1所示[18].實(shí)驗(yàn)中空氣由高50 mm,寬40 mm 的進(jìn)口通入燃燒室.自燃燒室入口60 mm 后燃燒室呈3°擴(kuò)張角.燃燒室中心處一長為32 mm,底部高6 mm 的楔形支板置于入口后35 mm 處.支板底部均勻分布15 個(gè)直徑為1.0 mm 間距為2.4 mm 的圓形噴孔以注入燃料-氫氣.實(shí)驗(yàn)空氣來流進(jìn)入燃燒室前經(jīng)過加熱補(bǔ)氧以模擬飛行條件下的氣體狀態(tài),污染空氣中氧氣和水的質(zhì)量分?jǐn)?shù)分別為0.232 和0.032,其余為氮?dú)?經(jīng)過拉瓦爾噴管膨脹加速后以馬赫數(shù)Ma=2.0,靜溫T=340 K,靜壓P=0.1 MPa 通入燃燒室內(nèi).實(shí)驗(yàn)中壓縮氫氣從噴孔中以聲速噴出,其靜溫和靜壓分別為250 K 和0.1 MPa.

圖1 DLR 中心支板超燃沖壓發(fā)動(dòng)機(jī)幾何模型(單位:mm)Fig.1 Schematic of DLR scramjet combustor with center strut (unit:mm)

考慮到燃燒室的結(jié)構(gòu)的對稱性,燃燒室僅取中心3 個(gè)噴孔寬度的區(qū)域作計(jì)算域以減少計(jì)算開銷.計(jì)算域內(nèi)的網(wǎng)格結(jié)構(gòu)如圖2 所示,總網(wǎng)格數(shù)量2.78×107.大部分計(jì)算域使用非結(jié)構(gòu)多面體網(wǎng)格剖分,在支板與主流邊界層區(qū)域使用部分楔形和棱柱形網(wǎng)格填充,對支板、中心射流以及燃燒反應(yīng)所在燃燒室中心區(qū)域(燃燒反應(yīng)區(qū))采用平均0.15 mm 網(wǎng)格進(jìn)行加密,邊界層網(wǎng)格由15 層指數(shù)增長型膨脹層網(wǎng)格生成,無量綱近壁距離y* <1.采用3 套單元數(shù)量分別為1.5×107(縮寫為15.3 M)、2.95×107(縮寫為29.5 M)和4.82×107(縮寫為48.2 M)的網(wǎng)格進(jìn)行網(wǎng)格獨(dú)立性驗(yàn)證,計(jì)算得到的冷流狀態(tài)下壁面沿程壓力分布如圖3 所示,結(jié)果表明較粗的兩套網(wǎng)格相對最細(xì)網(wǎng)格的平均相對誤差分別為2.8%和0.7%,最大誤差分別為33.8%和5.7%.因此,采用2.95×108網(wǎng)格可在較小的計(jì)算量下滿足網(wǎng)格獨(dú)立性要求.

圖2 計(jì)算域整體及燃料噴口附近網(wǎng)格結(jié)構(gòu)Fig.2 Grid structures of overall computational domain and near the fuel nozzle

圖3 網(wǎng)格獨(dú)立性驗(yàn)證Fig.3 Grid independence verification

2 結(jié)果與討論

2.1 數(shù)值模型驗(yàn)證

在無反應(yīng)流條件下,對比實(shí)驗(yàn)紋影[18]與數(shù)值紋影如圖4(a).由圖可以看出,本文模擬結(jié)果較好地復(fù)現(xiàn)了實(shí)驗(yàn)紋影中呈現(xiàn)出的復(fù)雜波系結(jié)構(gòu),包括由自中心支板前緣形成的兩道入射斜激波,分別經(jīng)過上下壁面的反射以及透射通過中心射流區(qū)域后再次入射到上下壁面,其中在上壁面形成小分離泡,并由此分別形成分離激波和反射激波.在支板尾部上下轉(zhuǎn)折點(diǎn)形成膨脹扇,隨后在中心射流的作用下流向偏轉(zhuǎn)形成上下兩道壓縮波并分別向上下壁面入射進(jìn)而發(fā)生反射.此外,中心燃料噴孔噴出的高速氫氣與支板尾部的低速空氣流發(fā)生剪切作用并在K-H 不穩(wěn)定性和后續(xù)的激波干擾下發(fā)展成一系列復(fù)雜渦結(jié)構(gòu).為定量驗(yàn)證分析數(shù)值計(jì)算結(jié)果,針對不同流向位置中心截面處平均速度分布,與實(shí)驗(yàn)結(jié)果及Fureby等[19]的LES 計(jì)算結(jié)果進(jìn)行了對比,如圖4(b)所示.圖中數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好,尤其在中下游區(qū)域.而在靠近氫氣射流出口(x=11 mm)處,數(shù)值模擬結(jié)果相較于實(shí)驗(yàn)結(jié)果中心射流速度低,且兩側(cè)低速尾跡區(qū)速度更低,Fureby 等[19]的計(jì)算結(jié)果也觀察到相同的趨勢.

圖4 冷態(tài)流場實(shí)驗(yàn)和數(shù)值(a)紋影與(b)不同位置流向平均速度(U x)對比Fig.4 Comparison of (a) schlieren shadow and (b) average flow velocity (Ux) at different positions between experiment and numerical results

圖5 所示為燃燒狀態(tài)下實(shí)驗(yàn)與數(shù)值紋影所呈現(xiàn)的流場結(jié)構(gòu)以及OH-PLIF 與數(shù)值模擬所得OH 分布的對比.由圖可知數(shù)值與實(shí)驗(yàn)所得流場結(jié)構(gòu)基本吻合.與冷流狀態(tài)相比,燃燒狀態(tài)下的波系結(jié)構(gòu)明顯不同.其中中心支板前緣所形成的斜激波經(jīng)過壁面反射與中心區(qū)域火焰面接觸后發(fā)生反射,與冷流場情況下透射穿過中心射流后與壁面進(jìn)一步反射從而形成更復(fù)雜的波系結(jié)構(gòu)情況不同.此外,由于燃燒釋熱導(dǎo)致的膨脹作用,支板尾部的回流區(qū)增寬,導(dǎo)致支板端部角點(diǎn)附近的膨脹扇結(jié)構(gòu)顯著減弱,同時(shí)其與中心射流所夾斜激波也消失.此外,對比OH 分布可以發(fā)現(xiàn),數(shù)值模擬準(zhǔn)確捕捉了支板后部火焰的推舉與局部熄火現(xiàn)象,且較準(zhǔn)確地預(yù)測了剪切層兩側(cè)火焰面出現(xiàn)交互的位置與火焰結(jié)構(gòu).

圖5 燃燒狀態(tài)實(shí)驗(yàn)與數(shù)值紋影及OH 分布對比Fig.5 Comparison of schlieren and OH distribution between the experiment and simulation

圖6 通過對比實(shí)驗(yàn)與數(shù)值計(jì)算所得不同位置中心截面的平均溫度分布,從定量上分析對比DZFM模型對反應(yīng)流模擬的準(zhǔn)確性.由圖可知,在支板尾部附近(x=11 mm),兩側(cè)剪切層區(qū)域數(shù)值模擬所得火焰平均溫度與實(shí)驗(yàn)吻合較好,而在中心區(qū)域未能捕捉到中心射流明顯的升溫,在基于DLR 燃燒室的其他LES[19,36]結(jié)果中亦出現(xiàn)相似的結(jié)果,上述誤差可能與上述數(shù)值模擬中均采用中心部位3 個(gè)噴注孔作為周期性對稱計(jì)算域的簡化有關(guān).在x=58 mm 位置,數(shù)值模擬所得溫度分布及極值與實(shí)驗(yàn)結(jié)果相當(dāng),火焰寬度略大于實(shí)驗(yàn)結(jié)果.在x=166 mm 位置,數(shù)值模擬與實(shí)驗(yàn)結(jié)果吻合較好.總體而言,本文模型(DZFM) 所得結(jié)果與實(shí)驗(yàn)結(jié)果吻合良好,與其他LES 模擬的結(jié)果相當(dāng).

圖6 燃燒條件下實(shí)驗(yàn)與數(shù)值計(jì)算所得不同流向位置平均溫度分布Fig.6 Average temperature distribution at different flow directions obtained from experiments and numerical calculations under the combustion condition

2.2 分區(qū)參量選取分析

圖7 所示為OH 質(zhì)量分?jǐn)?shù)(YOH)及其基于混合分?jǐn)?shù)的條件均值(〈YOH|η〉)在混合分?jǐn)?shù)空間(ξ)內(nèi)的分布情況.由圖可知,OH 的分布在混合分?jǐn)?shù)空間內(nèi)具有明顯的聚集性,也即火焰面內(nèi)反應(yīng)狀態(tài)與混合分?jǐn)?shù)具有強(qiáng)相關(guān)性,但其相對條件均值仍存在顯著的脈動(dòng),尤其在化學(xué)當(dāng)量混合分?jǐn)?shù)(ξst)附近.這是因?yàn)镈LR 中心支板超聲速燃燒室內(nèi)湍流與化學(xué)反應(yīng)交互作用(TCI)關(guān)系復(fù)雜,受多種因素影響,在相同混合分?jǐn)?shù)條件下可能存在多種不同的火焰模式.理論上火焰面依賴于多個(gè)變量,即多條件矩.然而多條件矩不僅建模尤為復(fù)雜而且計(jì)算代價(jià)巨大.通過將影響火焰面的因素納入分區(qū)指標(biāo),對整個(gè)流場而言相當(dāng)于將多條件矩具體化,對局部流場而言相當(dāng)于將多條件矩簡化為僅依賴于混合分?jǐn)?shù)的單一條件矩.多重指標(biāo)分區(qū)的引入極大簡化了建模和計(jì)算,多重分區(qū)指標(biāo)的選取對保證計(jì)算結(jié)果的準(zhǔn)確性尤為重要.為限制反應(yīng)狀態(tài)參數(shù)在混合分?jǐn)?shù)空間內(nèi)的強(qiáng)脈動(dòng),可選取與反應(yīng)狀態(tài)以及與流動(dòng)狀態(tài)相關(guān)的參數(shù)作為分區(qū)的指標(biāo),將物理空間進(jìn)行動(dòng)態(tài)分區(qū).非預(yù)混火焰因邊混合邊燃燒而具有顯著的流向分區(qū)特征,因此流向坐標(biāo)X通常作為基礎(chǔ)分區(qū)變量.在同一流向坐標(biāo)(X)區(qū)域,可能同時(shí)存在未燃和已燃兩種反應(yīng)區(qū)域,這時(shí)需要引入表征化學(xué)反應(yīng)進(jìn)度的反應(yīng)進(jìn)度變量(Cz)來區(qū)別不同的反應(yīng)狀態(tài),針對每種狀態(tài)采用單獨(dú)的火焰面去描述,這對于準(zhǔn)確表征預(yù)混區(qū)域的反應(yīng)狀態(tài)尤為重要.這里采用最終生成產(chǎn)物(H2O)中H 元素的含量計(jì)算反應(yīng)進(jìn)度變量

圖7 OH 質(zhì)量分?jǐn)?shù)(YOH)及其基于混合分?jǐn)?shù)空間的條件均值(〈 YOH|η〉)在混合分?jǐn)?shù)空間(ξ)的分布散點(diǎn)圖Fig.7 Scatter plot of the distribution of OH mass fraction (YOH) and its conditional mean based on mixed fraction space (〈 YOH|η〉) in mixed fraction space (ξ)

其中Yα表示組分 α 的質(zhì)量分?jǐn)?shù),Hα表示H 元素在組分 α 中的質(zhì)量分?jǐn)?shù).在同一X坐標(biāo)區(qū)域,也可能存在純?nèi)剂?、純空氣和部分混合物等不同的混合狀態(tài),不同的混合狀態(tài)對應(yīng)截然不同的反應(yīng)狀態(tài),因此需要引入混合分?jǐn)?shù)(ξ)以區(qū)分不同的混合狀態(tài).考慮到反應(yīng)釋熱與反應(yīng)和湍流渦量的強(qiáng)相關(guān)性,選擇熱釋放速率(dQ)作為參數(shù)之一.在超聲速流動(dòng)中,因氣動(dòng)加熱導(dǎo)致流體本身動(dòng)能的變化率與燃燒的熱釋放速率相當(dāng),例如激波前后會出現(xiàn)溫度以及壓力的劇烈變化,邊界層因減速也會出現(xiàn)高焓層.此外,對于超聲速湍流燃燒,流場的可壓縮性也對化學(xué)反應(yīng)有顯著的影響[37].因此有必要引入馬赫數(shù)(Ma)這一分區(qū)指標(biāo)以區(qū)分亞聲速和超聲速區(qū)域.增加某一分區(qū)指標(biāo)的數(shù)目,可以降低分區(qū)內(nèi)條件脈動(dòng)幅值從而提高計(jì)算的準(zhǔn)確性.然而僅僅通過增加單一的數(shù)目而不引入其他分區(qū)指標(biāo),難以保證影響當(dāng)前分區(qū)上反應(yīng)狀態(tài)的條件統(tǒng)計(jì)均一性.

分區(qū)參數(shù)在混合分?jǐn)?shù)與OH 質(zhì)量分?jǐn)?shù)空間內(nèi)的分布如圖8 所示.為實(shí)現(xiàn)在每個(gè)動(dòng)態(tài)分區(qū)內(nèi)條件脈動(dòng)足夠小,即混合分?jǐn)?shù)與化學(xué)反應(yīng)參量關(guān)聯(lián)關(guān)系(函數(shù)或映射)的統(tǒng)計(jì)均一性,分區(qū)參量在YOH與 ξ 所張成的二維空間(YOH-ξ)內(nèi)的分布應(yīng)盡量滿足其梯度最大方向同YOH與 ξ 間所構(gòu)成的映射關(guān)系曲線的切線方向垂直.由圖8 可知,熱釋放速率、馬赫數(shù)以及混合分?jǐn)?shù)在YOH-ξ 空間中的分布,在大部分區(qū)域梯度最大方向同YOH與 ξ 間所構(gòu)成的映射關(guān)系曲線的切線方向接近垂直,而反應(yīng)進(jìn)度變量在大部分區(qū)域接近于平行.因此,dQ,Ma以及X可作為分區(qū)參量,對計(jì)算網(wǎng)格進(jìn)行分區(qū)聚類.

圖8 熱釋放速率(dQ),馬赫數(shù)(Ma),流向坐標(biāo)(X)以及反應(yīng)進(jìn)度變量(Cz)在混合分?jǐn)?shù)與OH 質(zhì)量分?jǐn)?shù)空間內(nèi)的分布云圖Fig.8 Scatter plot of OH distributed in mixture fraction space rendered by heat release rate (dQ),Mach number (Ma),flow direction coordinate (X),and reaction progress variable (Cz)

為進(jìn)一步驗(yàn)證上述分區(qū)參量對條件脈動(dòng)的抑制作用,分別選取上述單一參量以及組合參量進(jìn)行分區(qū),得到分區(qū)后單一分區(qū)內(nèi)化學(xué)反應(yīng)狀態(tài)參量(YOH)及其分區(qū)條件平均(QOH=〈YOH|η〉zone)在混合分?jǐn)?shù)空間內(nèi)的分布情況,如圖9 所示.由圖中YOH的分布可知,在中心支板尾部區(qū)域附近(x∈[0,20]),也即火焰抬升區(qū)域附近,由于火焰發(fā)生局部熄火和再點(diǎn)火,在小的流向空間分區(qū)內(nèi)仍然具有明顯的條件脈動(dòng).通過引入熱釋放速率(dQ)和馬赫數(shù)(Ma)作為分區(qū)指標(biāo)可有效降低條件脈動(dòng)幅值,其中熱釋放速率(dQ)分區(qū)參量對降低條件脈動(dòng)的作用更為顯著,尤其在貧燃區(qū)域(ξ <ξst=0.028).在流場中下游區(qū)域,僅以流場坐標(biāo)(X)作為分區(qū)指標(biāo)即可較好限制反應(yīng)參量在混合分?jǐn)?shù)空間內(nèi)的條件脈動(dòng)值,實(shí)現(xiàn)分區(qū)內(nèi)反應(yīng)狀態(tài)的均一性.同時(shí)在此基礎(chǔ)上引入馬赫數(shù)和熱釋放速率作分區(qū)指標(biāo)可進(jìn)一步抑制中下游區(qū)域的條件脈動(dòng).綜上可知,以流向坐標(biāo)與熱釋放速率作為分區(qū)指標(biāo)可顯著抑制全場的分區(qū)條件脈動(dòng).下文將進(jìn)一步分析不同分區(qū)策略對燃燒室模擬結(jié)果的影響.

圖9 各分區(qū)內(nèi)OH 基團(tuán)質(zhì)量分?jǐn)?shù)(YOH)與其分區(qū)條件均值(QOH=〈YOH|η〉zone)在混合分?jǐn)?shù)空間內(nèi)的離散分布情況Fig.9 Discrete distribution of OH mass fraction (YOH) and its condition mean (QOH=〈YOH|η〉zone) in mixture fraction space within each partition

2.3 多參量分區(qū)的影響分析

圖10 所示為選取不同分區(qū)指標(biāo)進(jìn)行多參量分區(qū)后,DZFM 模擬所得燃燒室內(nèi)不同流向位置火焰溫度分布.根據(jù)上一節(jié)中不同分區(qū)參數(shù)對條件脈動(dòng)的影響分析,分別對比了3 種由流向坐標(biāo)X,混合分?jǐn)?shù) ξ,馬赫數(shù)Ma以及熱釋放速率dQ構(gòu)成的組合分區(qū)方式:流向坐標(biāo)與混合分?jǐn)?shù)組合分區(qū)(X+ξ);流向坐標(biāo)、混合分?jǐn)?shù)以及馬赫數(shù)組合分區(qū)(X+ξ +Ma);流向坐標(biāo)、混合分?jǐn)?shù)以及熱釋放速率組合分區(qū)(X+ξ +dQ).由圖可知,不同分區(qū)指標(biāo)組合成的多參量分區(qū)方法對不同流向位置下火焰溫度的分布的影響不同,在支板尾部區(qū)域(X=11 mm)以X+ξ和X+ξ +Ma多參數(shù)分區(qū)的數(shù)值模擬所得火焰溫度在剪切層區(qū)域明顯高于實(shí)驗(yàn)結(jié)果,而X+ξ +dQ多分區(qū)方法所得剪切層火焰溫度更接近于實(shí)驗(yàn)結(jié)果.在X=58 mm 位置,3 種不同分區(qū)方式對火焰溫度分布無明顯影響.在下游區(qū)域(X=166 mm),以X+ξ +dQ和X+ξ +Ma進(jìn)行多參數(shù)分區(qū)的計(jì)算結(jié)果與實(shí)驗(yàn)更相符,而以X+ξ 組合進(jìn)行多參數(shù)分區(qū)所得火焰溫度較上述兩者更低且與實(shí)驗(yàn)偏差更大.對比采用X+ξ (51×91)雙分區(qū)參量方法與多分區(qū)參量方法X+ξ +Ma(51×91×19)和X+ξ +dQ(51×91×19),計(jì)算單流通時(shí)間(TFT=0.41 ms)后兩者的CPU 計(jì)算時(shí)間分別增加5.65%和5.28%.計(jì)算中針對空分區(qū)即沒有對應(yīng)CFD 網(wǎng)格的分區(qū)空間采取了“凍結(jié)”計(jì)算的處理,因此計(jì)算時(shí)間并沒有隨著分區(qū)總數(shù)目的增加而線性增加,從而顯著提高了計(jì)算效率.

圖10 不同多參量分區(qū)對沿流向位置中心處火焰平均溫度分布的影響Fig.10 The effect of different multiparameter zoning on the average temperature distribution of flame at the center of the flow direction position

為進(jìn)一步分析不同多參量分區(qū)對火焰結(jié)構(gòu)的影響,研究了3 種分區(qū)方法下的中心截面的瞬態(tài)溫度分布及基于改進(jìn)型TFI (Takeno flame index)指數(shù)[38]的火焰模式分布.改進(jìn)型TFI具體定義如下

當(dāng)TFI≥ 0.5,為富燃預(yù)混火焰;當(dāng)TFI≤-0.5 時(shí),為貧燃預(yù)混火焰;其他情況(-0.5 ≤TFI≤0.5)則為擴(kuò)散火焰.由圖11 可以看出,在支板尾部氫氣射流區(qū)域以及中心射流兩側(cè)的剪切層區(qū)域溫度較低且處于預(yù)混火焰狀態(tài).富燃預(yù)混火焰區(qū)同時(shí)也是火焰推舉區(qū)域,在X+ξ 分區(qū)條件下這一區(qū)域火焰推舉高度更高且兩側(cè)剪切層區(qū)域火焰面連續(xù)無明顯的局部熄火行為.而X+ξ +Ma和X+ξ +dQ分區(qū)條件下,均能觀察到剪切層火焰的局部熄火現(xiàn)象,且中心預(yù)混火焰區(qū)火焰前緣向上游移動(dòng),尤其在X+ξ +dQ分區(qū)條件下局部熄火位置更接近支板尾部,這與圖10 中X=11 mm 處較低的平均溫度分布相符.同時(shí),上述分區(qū)的作用與圖9 中X∈[0,20]區(qū)間內(nèi)的OH 分布相符,即通過引入dQ和Ma作為分區(qū)指標(biāo)可以降低條件脈動(dòng),從而更準(zhǔn)確地捕捉火焰在支板尾部預(yù)混模式下的局部熄火行為.在燃燒室中部區(qū)域,燃燒由擴(kuò)散火焰模式主導(dǎo),因此X+ξ 分區(qū)可以準(zhǔn)確模擬當(dāng)?shù)厝紵隣顟B(tài)而無須引入Ma和dQ分區(qū),這與圖10 中X=58 mm 處的溫度曲線重疊相符,同時(shí)也驗(yàn)證了圖9 中X∈[50,60] 區(qū)間內(nèi)OH 與 ξ 的強(qiáng)相關(guān).而在燃燒室下游高雷諾數(shù)區(qū)域,火焰在局部高馬赫來流所致的強(qiáng)拉伸、剪切作用下出現(xiàn)局部熄火同時(shí)呈現(xiàn)零散分布的貧燃預(yù)混火焰區(qū).由圖10 和圖11 可以看出,在這一區(qū)域X+ξ 分區(qū)所得結(jié)果與X+ξ +Ma和X+ξ +dQ分區(qū)條件下所得結(jié)果有較大差異,后者具有更寬的燃燒高溫區(qū)且更符合實(shí)驗(yàn)結(jié)果.

圖11 不同多參量分區(qū)條件下中心截面處火焰溫度及火焰預(yù)混模式分布Fig.11 Distribution of flame temperature and flame premixing mode at the central cross-section under different multiparameter zoning conditions

3 結(jié)論

本文基于動(dòng)態(tài)分區(qū)火焰面模型(DZFM)采用詳細(xì)氫氣化學(xué)反應(yīng)機(jī)理對DLR 超燃沖壓發(fā)動(dòng)機(jī)內(nèi)湍流燃燒進(jìn)行改進(jìn)延遲分離渦(IDDES)模擬.模擬結(jié)果準(zhǔn)確捕捉了燃燒室冷流狀態(tài)下實(shí)驗(yàn)紋影中觀察到的復(fù)雜波系結(jié)構(gòu)以及中心支板穩(wěn)焰狀態(tài)下的火焰推舉行為,對火焰結(jié)構(gòu)的模擬結(jié)果與實(shí)驗(yàn)符合較好.通過分析化學(xué)反應(yīng)標(biāo)量(YOH)在混合分?jǐn)?shù)空間(ξ)的分布與不同分區(qū)參量在(YOH-ξ)空間內(nèi)梯度方向,發(fā)現(xiàn)流向坐標(biāo)X,馬赫數(shù)Ma以及燃燒熱釋放速率dQ可抑制分區(qū)內(nèi)反應(yīng)標(biāo)量的條件脈動(dòng).通過進(jìn)一步分析不同多參數(shù)分區(qū)對DZFM 模擬所得火焰結(jié)構(gòu)的影響,發(fā)現(xiàn)在出現(xiàn)局部熄火的低雷諾數(shù)富燃預(yù)混火焰區(qū)以及高雷諾數(shù)貧燃預(yù)混火焰區(qū),條件脈動(dòng)顯著增強(qiáng).在上述區(qū)域,僅以X和 ξ 進(jìn)行網(wǎng)格分區(qū)存在較大的誤差,而額外引入dQ進(jìn)行三維參數(shù)空間分區(qū)可通過抑制條件脈動(dòng)提升DZFM 模型在低雷諾數(shù)富燃預(yù)混火焰以及高雷諾數(shù)貧燃預(yù)混火焰區(qū)域的準(zhǔn)確性.