基于動(dòng)態(tài)分區(qū)概念的高超聲速燃燒大渦模擬1)

姚衛(wèi) 劉 杭 張 政 肖雅彬 岳連捷,?

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

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

引言

2004 年,X-43 A 在35 km 高度的臨近空間實(shí)現(xiàn)了接近馬赫數(shù)10 的高超聲速飛行[1],取得了大氣層內(nèi)吸氣式飛行的最快速度.近年來(lái),高馬赫數(shù)飛行器引起了世界各國(guó)的高度重視,成為各空天科技強(qiáng)國(guó)競(jìng)相爭(zhēng)奪的科技制高點(diǎn).飛行馬赫數(shù)上限達(dá)到12 甚至更高的超燃沖壓發(fā)動(dòng)機(jī)是實(shí)現(xiàn)可重復(fù)使用天地往返飛行器的關(guān)鍵[2].目前反射激波風(fēng)洞和激波膨脹風(fēng)洞是能夠產(chǎn)生馬赫數(shù)大于8 高焓氣流的兩種主要途徑.由于高焓氣流地面模擬的困難,目前世界上僅有為數(shù)不多的高馬赫數(shù)風(fēng)洞,例如美國(guó)的NASA HYPULSE 激波風(fēng)洞[3]、日本的JAXA HIEST 自由活塞驅(qū)動(dòng)激波風(fēng)洞[4]和澳大利亞昆士蘭大學(xué)T4 反射激波風(fēng)洞[5].其中在T4 風(fēng)洞中針對(duì)飛行馬赫數(shù)7.5～ 12[6-12]工況開(kāi)展了一系列的矩形轉(zhuǎn)橢圓(REST)高超聲速燃燒室縮比模型測(cè)試.由于高超聲速實(shí)驗(yàn)測(cè)量手段和可獲得數(shù)據(jù)均十分有限,計(jì)算流體力學(xué)(CFD)成為燃燒內(nèi)部流場(chǎng)分析和發(fā)動(dòng)機(jī)性能優(yōu)化的必需手段.目前文獻(xiàn)中已經(jīng)開(kāi)展了馬赫數(shù)8,12,13 和15[13-15]等對(duì)應(yīng)飛行條件下的高超聲速超燃沖壓發(fā)動(dòng)機(jī)模擬,極大彌補(bǔ)了實(shí)驗(yàn)測(cè)量信息的不足.然而,隨著飛行馬赫數(shù)的提高,湍流強(qiáng)度增加,流動(dòng)和化學(xué)反應(yīng)特征時(shí)間相互重疊區(qū)域增加.同時(shí),為了增加流動(dòng)滯留時(shí)間以提高混合和燃燒效率,相比于馬赫數(shù)4～ 8的常規(guī)超燃沖壓發(fā)動(dòng)機(jī),馬赫數(shù)大于8 的高馬赫數(shù)超燃沖壓發(fā)動(dòng)機(jī)流道尺寸大幅增加,意味著巨大的全尺寸模擬計(jì)算代價(jià).高馬赫數(shù)條件下內(nèi)流與外流的耦合效應(yīng)也顯著增強(qiáng),為了提高模擬置信度一般計(jì)算區(qū)域需要至少包含部分飛行器前體.上述獨(dú)特性意味著高馬赫數(shù)超燃沖壓發(fā)動(dòng)機(jī)數(shù)值模擬在兼顧模型保真度(復(fù)雜多重湍流-化學(xué)反應(yīng)交互作用模式的表征)與計(jì)算代價(jià)(耦合詳細(xì)化學(xué)反應(yīng)動(dòng)力學(xué)的高解析度大渦模擬)兩方面均面臨著極大的挑戰(zhàn).

空間上復(fù)雜而時(shí)間上多變的湍流化學(xué)反應(yīng)交互作用關(guān)系(turbulence-chemistry interaction,TCI)是超聲速燃燒室內(nèi)部湍流燃燒的主要特點(diǎn).在超燃沖壓發(fā)動(dòng)機(jī)中,其內(nèi)部流場(chǎng)不僅存在混合控制為主的快速反應(yīng)區(qū)和化學(xué)反應(yīng)控制為主的慢速反應(yīng)區(qū),甚至還存在預(yù)混、非預(yù)混和部分預(yù)混等復(fù)雜的湍流燃燒模式[16-20].跨越多個(gè)時(shí)間和空間尺度的湍流和化學(xué)反應(yīng)不僅自身的模擬極為困難,而且其交互作用所引起的復(fù)雜耦合效應(yīng)進(jìn)一步加劇了其模擬挑戰(zhàn)性.湍流燃燒模擬的核心在于對(duì)復(fù)雜的湍流-化學(xué)反應(yīng)交互作用關(guān)系準(zhǔn)確、高效建模.目前的湍流燃燒模型大致可分為兩類(lèi):有限速率類(lèi)和守恒標(biāo)量類(lèi).有限速率類(lèi)模型包括近似層流模型(QL)[21]、部分?jǐn)嚢璺磻?yīng)器模型(PaSR)[22]、渦耗散/破碎模型(EDC/EBU)[23]、增厚火焰面模型[24]和線(xiàn)性渦模型(LEM)[25]等.守恒標(biāo)量類(lèi)模型包括各類(lèi)火焰面模型[26-28]、條件矩封閉模型(CMC)[29]、假定[30-31]/輸運(yùn)[32]概率密度函數(shù)模型等.有限速率類(lèi)模型適用范圍較廣但一般需要在每個(gè)網(wǎng)格節(jié)點(diǎn)上求解所有的組分輸運(yùn)方程,因而計(jì)算量較大.守恒標(biāo)量類(lèi)模型通過(guò)對(duì)流動(dòng)和化學(xué)反應(yīng)的解耦操作實(shí)現(xiàn)了湍流燃燒的降維求解,一般可以顯著降低計(jì)算量.然而,流動(dòng)-化學(xué)反應(yīng)可以解耦的前提是化學(xué)反應(yīng)受湍流的影響模式單一,以便選擇適合的守恒變量-熱力學(xué)參數(shù)空間分布去表征當(dāng)前化學(xué)反應(yīng)狀態(tài).以火焰面模型為例,當(dāng)湍流化學(xué)反應(yīng)交互處于火焰面或薄反應(yīng)區(qū)模式(戴姆克拉數(shù)Da＞10,卡爾洛維奇數(shù)Ka＜100)時(shí),適宜采用計(jì)算簡(jiǎn)便的穩(wěn)態(tài)火焰面模型;而當(dāng)湍流化學(xué)反應(yīng)模式處于慢反應(yīng)區(qū)模式(Ka＞100,Da＜10)時(shí),此時(shí)化學(xué)反應(yīng)速率小于湍流混合速率,需要采用非穩(wěn)態(tài)火焰面(UFM)或火焰面/進(jìn)度變量(FPV)模型[33];在發(fā)動(dòng)機(jī)燃燒室中還經(jīng)常存在局部預(yù)混或部分預(yù)混模式,此時(shí)需要采用火焰面模型的其他變種如G 方程模型[34]或火焰面生成流形模型(FGM)[35]等.

與亞聲速流動(dòng)相比,超聲速流動(dòng)有如下特性:(1)滯留時(shí)間極短(毫秒量級(jí)),和碳?xì)淙剂系幕瘜W(xué)反應(yīng)特征時(shí)間相當(dāng);(2)雷諾數(shù)一般高于5 × 105,最小湍流尺度遠(yuǎn)小于火焰面厚度.這意味著超聲速燃燒中湍流和化學(xué)反應(yīng)因特征時(shí)間和特征尺度相當(dāng),存在強(qiáng)烈的耦合效應(yīng).超聲速燃燒場(chǎng)在不同局部區(qū)域存在差異性極大的湍流-化學(xué)反應(yīng)交互作用模式,反應(yīng)狀態(tài)參數(shù)對(duì)守恒標(biāo)量的依賴(lài)存在極大的空間不均勻性,因此采用單一火焰面去描述整場(chǎng)湍流反應(yīng)狀態(tài)將會(huì)帶來(lái)嚴(yán)重的誤差,甚至一般認(rèn)為火焰面模型不適用于超聲速燃燒模擬.盡管目前文獻(xiàn)中已經(jīng)發(fā)展了基于多重火焰面的表征火焰面模型[36],但該模型仍然把整場(chǎng)作為一個(gè)火焰面分區(qū),并不能很好地區(qū)分受當(dāng)?shù)赝牧饔绊懙幕瘜W(xué)反應(yīng)狀態(tài).而如果能夠根據(jù)局部流場(chǎng)狀態(tài)自適應(yīng)地采用單獨(dú)的守恒標(biāo)量反應(yīng)狀態(tài)參數(shù)空間去表征,即局部表征火焰面模型的概念,則可極大地提高整場(chǎng)計(jì)算的準(zhǔn)確性.局部火焰面模型的理論依據(jù)是,盡管反應(yīng)狀態(tài)參數(shù)對(duì)守恒標(biāo)量的依賴(lài)存在空間不均勻性,但針對(duì)局部流場(chǎng)可以假設(shè)該依賴(lài)關(guān)系具有統(tǒng)計(jì)均一性.分區(qū)模擬從物理本質(zhì)上提高了火焰面模型在超聲速燃燒模擬中的適用性,從而為高效高保真超聲速燃燒大渦模擬提供了一個(gè)獨(dú)特的解決思路.從針對(duì)各局部分區(qū)分別建立守恒標(biāo)量與反應(yīng)狀態(tài)參數(shù)之間關(guān)聯(lián)的角度看,分區(qū)湍流燃燒模型與條件矩封閉模型類(lèi)似[29].但是分區(qū)湍流燃燒模型的分區(qū)是基于局部流場(chǎng)的湍流-化學(xué)反應(yīng)交互作用模式并且隨流場(chǎng)演化而動(dòng)態(tài)調(diào)整,而條件矩封閉模型一般是單純基于幾何坐標(biāo)位置(多為笛卡爾正交坐標(biāo))的固定分區(qū).對(duì)于瞬態(tài)流場(chǎng)而言,局部湍流化學(xué)反應(yīng)交互作用模式是動(dòng)態(tài)變化的.在分區(qū)火焰面模型中,各分區(qū)的守恒標(biāo)量-狀態(tài)參數(shù)空間分布由一個(gè)表征火焰面(RIF)[37]描述.表征火焰面描述的并非實(shí)際物理上的火焰面結(jié)構(gòu),而是基于統(tǒng)計(jì)規(guī)律建立的近似映射關(guān)系.因此該分區(qū)表征火焰面模型并不局限于傳統(tǒng)的薄火焰面假設(shè),而是普遍適用于非平衡、瞬態(tài)和非均勻反應(yīng)過(guò)程的模擬.分區(qū)湍流燃燒模擬的前提是分區(qū)內(nèi)湍流化學(xué)反應(yīng)交互作用模式局部近似均勻,不合適的分區(qū)會(huì)導(dǎo)致表征火焰面與實(shí)際的狀態(tài)參數(shù)在守恒標(biāo)量空間內(nèi)分布差別較大.為了準(zhǔn)確描述局部化學(xué)反應(yīng)動(dòng)力學(xué)狀態(tài),需要自適應(yīng)地進(jìn)行分區(qū)動(dòng)態(tài)調(diào)整,以便使得各分區(qū)內(nèi)化學(xué)反應(yīng)受湍流的影響模式局部單一.

近年來(lái),分區(qū)類(lèi)燃燒模擬方法在發(fā)動(dòng)機(jī)燃燒模擬中得到了較多關(guān)注[38-45].目前的分區(qū)燃燒模擬主要是依據(jù)化學(xué)反應(yīng)熱力學(xué)狀態(tài)對(duì)燃燒場(chǎng)進(jìn)行分區(qū)或聚類(lèi),以便采用簡(jiǎn)化化學(xué)反應(yīng)計(jì)算模型(如凍結(jié)反應(yīng)、均勻混合和平衡態(tài)假設(shè)等)或基于分區(qū)對(duì)化學(xué)反應(yīng)進(jìn)行整體求解.簡(jiǎn)單的分區(qū)一般基于幾何空間分為未燃區(qū)、反應(yīng)區(qū)/火焰鋒面區(qū)和已燃區(qū)等[38-42,45],復(fù)雜的分區(qū)對(duì)具有相似反應(yīng)狀態(tài)的單元(可以幾何不相鄰)進(jìn)行聚類(lèi)操作[43-44,46-48].然而現(xiàn)有分區(qū)燃燒模擬并未考慮湍流效應(yīng),確切地說(shuō)應(yīng)該稱(chēng)之為“化學(xué)分區(qū)”.已經(jīng)有部分研究[49-50]在多燃燒模型耦合計(jì)算方面開(kāi)展了初步的探索,但是其所采用的方法需要提前初步評(píng)估各模型性能以指定分區(qū),難以在后續(xù)計(jì)算中動(dòng)態(tài)調(diào)整分區(qū),因此不適用于瞬態(tài)性較強(qiáng)的超聲速燃燒流場(chǎng).對(duì)各模型結(jié)合當(dāng)前算例進(jìn)行預(yù)測(cè)評(píng)估也需要額外的計(jì)算代價(jià),一定程度上降低了多模型耦合計(jì)算的效率優(yōu)勢(shì).此外,有限速率類(lèi)模型因計(jì)算代價(jià)過(guò)大,難以耦合詳細(xì)的化學(xué)反應(yīng)機(jī)理,與火焰面模型耦合使用時(shí)會(huì)產(chǎn)生“短板效應(yīng)”,嚴(yán)重降低整體計(jì)算效率,因此具有不同理論基礎(chǔ)的湍流燃燒模型分區(qū)耦合應(yīng)用價(jià)值有限.

本文針對(duì)超聲速燃燒模擬提出了一種基于當(dāng)?shù)赝牧骰瘜W(xué)反應(yīng)交互作用關(guān)系進(jìn)行自適應(yīng)流場(chǎng)分區(qū)和基于局部表征火焰面的分區(qū)湍流燃燒模擬方法.該方法不僅可以極大提高存在復(fù)雜多重湍流燃燒模式的超聲速燃燒模擬的保真度,而且通過(guò)對(duì)局部湍流化學(xué)反應(yīng)的解耦可取得極高的計(jì)算效率,從而突破長(zhǎng)期以來(lái)制約超聲速燃燒模擬的流動(dòng)-化學(xué)反應(yīng)多尺度作用耦合建模這一關(guān)鍵瓶頸.研究進(jìn)一步結(jié)合該方法對(duì)不同構(gòu)型的高超聲速超燃沖壓發(fā)動(dòng)機(jī)進(jìn)行了全尺寸內(nèi)外流耦合一體化大渦模擬,并對(duì)其中的流動(dòng)、混合和燃燒機(jī)理進(jìn)行了深入分析,對(duì)比了兩種構(gòu)型高超聲速燃燒室的性能,展示了動(dòng)態(tài)分區(qū)湍流燃燒模擬方法在超聲速燃燒模擬基礎(chǔ)研究和超聲速燃燒室優(yōu)化設(shè)計(jì)方面的應(yīng)用.

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

1.1 動(dòng)態(tài)分區(qū)火焰面模型(DZFM)

定義當(dāng)前局部分區(qū)內(nèi)的表征火焰面結(jié)構(gòu)函數(shù)為Qα=〈Yα|ξ(x,t)=η,x∈zone〉,其中 η 為混合分?jǐn)?shù)空間的取樣變量,x代表空間物理坐標(biāo),x∈zone 表示條件平均統(tǒng)計(jì)局限于當(dāng)前區(qū)域.瞬態(tài)質(zhì)量分?jǐn)?shù)Yα與表征質(zhì)量分?jǐn)?shù)Qα的關(guān)聯(lián)為

考慮相變效應(yīng)和組分間差異擴(kuò)散的描述質(zhì)量守恒、瞬態(tài)混合分?jǐn)?shù) ξ 和組分質(zhì)量分?jǐn)?shù)Yα的完整控制方程如下

其中 ρ 為密度,U為速度向量,Dξ和Dα分別代表混合物、單獨(dú)組分的擴(kuò)散速率,Wα為化學(xué)反應(yīng)速率,為相變/凝固等引起的質(zhì)量增添,Yl,α為其他相(如液滴)中的組分 α 的質(zhì)量分?jǐn)?shù),ξl為其他相的整體混合分?jǐn)?shù),對(duì)于純?nèi)剂弦旱慰扇?ξl=1.

將式(1)差分后代入式(4),并使用式(2)和式(3),可得

對(duì)式(5)取條件平均 ξ (x,t)=η,和位于當(dāng)前分區(qū)內(nèi)x∈zone,可得表征火焰面條件組分Qα的最終控制方程為

Evap表示相變效應(yīng)貢獻(xiàn)項(xiàng),在不涉及兩相相變的情況下為零.χ是標(biāo)量耗散率定義為(?ξ)2,ρη=〈ρ|η〉,分區(qū)條件擴(kuò)散分區(qū)條件平均的引入意味著在每個(gè)分區(qū)內(nèi)統(tǒng)計(jì)上遵循式(6)定義的守恒性.式(6)中左端第二項(xiàng)表示局部小火焰在相鄰區(qū)域之間流動(dòng)的對(duì)流輸運(yùn),其物理意義為使得下游小火焰繼承上游小火焰的反應(yīng)狀態(tài),這是準(zhǔn)確模擬點(diǎn)火過(guò)程和火焰抬升現(xiàn)象的關(guān)鍵.式(6)右端第1 項(xiàng)和第2 項(xiàng)分別模擬了因標(biāo)量耗散和組分間差異擴(kuò)散引起的混合分?jǐn)?shù)空間內(nèi)組分?jǐn)U散,第3 項(xiàng)表示分區(qū)內(nèi)組分因化學(xué)反應(yīng)而發(fā)生的演化.

分區(qū)內(nèi)的化學(xué)反應(yīng)受分區(qū)條件平均溫度控制,在本模型中采用基于當(dāng)前分區(qū)內(nèi)瞬態(tài)焓統(tǒng)計(jì)的方法獲得[51].通過(guò)忽略焓脈動(dòng) (H′～O(0))[52-54],焓的概率密度分布 〈H|η〉 可以假設(shè)為以當(dāng)前Favre 平均值為中心的狄拉克 δ 函數(shù).因此條件焓的 〈H|η〉 分布可以基于歷史統(tǒng)計(jì)方法獲得

其中n為符合=η 的取樣空間數(shù)據(jù)點(diǎn).Thornber等[55]開(kāi)展的高分辨率大渦模擬計(jì)算中采用了類(lèi)似的統(tǒng)計(jì)方法以獲取燃燒模型所需要的條件平均量.條件溫度通過(guò)條件統(tǒng)計(jì)平均焓與混合分?jǐn)?shù)空間組分信息獲得,QT=f(〈H|η〉,Qα) .該統(tǒng)計(jì)模型極大節(jié)省了直接求解條件能量方程的計(jì)算成本,以及對(duì)超聲速可壓縮條件下條件能量源項(xiàng)建模的復(fù)雜度.該統(tǒng)計(jì)條件溫度模型的適用前提是分區(qū)內(nèi)包含足夠多的時(shí)間或空間采樣樣本,以及數(shù)值模擬能夠捕捉脈動(dòng)瞬態(tài)值,因此統(tǒng)計(jì)模型方法僅適用于大渦模擬(LES)和直接數(shù)值模擬(DNS)等高解析度模擬方法[55].

上式表明式 (8) 右端第1 項(xiàng)eY的作用為在每個(gè)分區(qū)對(duì)應(yīng)的采樣空間內(nèi)重新分布條件平均量,并且統(tǒng)計(jì)上重布效果相互抵消,分區(qū)內(nèi)整體質(zhì)量分?jǐn)?shù)脈動(dòng)守恒.通過(guò)動(dòng)態(tài)更新火焰面分區(qū)使得每個(gè)火焰面分區(qū)對(duì)應(yīng)一個(gè)較窄范圍的混合分?jǐn)?shù)子空間η ∈.其中為當(dāng)前分區(qū)平均混合分?jǐn)?shù),Δ ξ 為瞬態(tài)混合分?jǐn)?shù)的脈動(dòng)范圍.通過(guò)進(jìn)一步細(xì)化分區(qū),當(dāng)前分區(qū)混合分?jǐn)?shù)變化幅度趨近于零 Δ ξ →0,分區(qū)內(nèi)標(biāo)量概率密度分布坍縮為一個(gè)以為中心的狄拉克 δ 函數(shù).據(jù)此有,這意味著式 (8) 中的條件脈動(dòng)輸運(yùn)項(xiàng)通過(guò)動(dòng)態(tài)自適應(yīng)分區(qū)可以被忽略.

DZFM 模型在四維空間(三維物理空間加一維混合分?jǐn)?shù)空間)內(nèi)求解火焰面的演化.傳統(tǒng)的火焰面模型可以看作是一個(gè)把整個(gè)計(jì)算區(qū)域當(dāng)作一個(gè)單一分區(qū)的特殊分區(qū)火焰面模型,屬于低維流型模型(LDMM) 的一種.全組分輸運(yùn)湍流燃燒模型(如PaSR[22])需在三維物理空間離散網(wǎng)格單元逐一求解組分輸運(yùn)方程和化學(xué)反應(yīng),計(jì)算代價(jià)一般遠(yuǎn)高于火焰面等流動(dòng)-化學(xué)反應(yīng)解耦模型.DZFM 模型實(shí)現(xiàn)了流動(dòng)和化學(xué)反應(yīng)在當(dāng)前分區(qū)內(nèi)的局部解耦,因而只需在火焰面分區(qū)單元上求解火焰面結(jié)構(gòu)函數(shù)隨空間的演化.由于火焰面分區(qū)的空間分辨率通常遠(yuǎn)小于CFD 網(wǎng)格數(shù)目,所需求解的化學(xué)反應(yīng)體系大幅減少,因此整體上提高了計(jì)算效率.注意DZFM 模型的火焰面分區(qū)需隨著流場(chǎng)而動(dòng)態(tài)演化,以滿(mǎn)足湍流化學(xué)反應(yīng)交互模式局部單一的前提假設(shè).

1.2 分區(qū)自適應(yīng)化學(xué)(Z-DAC)

在流場(chǎng)求解過(guò)程中動(dòng)態(tài)地依據(jù)當(dāng)?shù)責(zé)峄瘜W(xué)反應(yīng)狀態(tài)進(jìn)行實(shí)時(shí)機(jī)理簡(jiǎn)化的方法稱(chēng)之為動(dòng)態(tài)自適應(yīng)化學(xué)(DAC)[56-57].自適應(yīng)化學(xué)方法通過(guò)發(fā)展并即時(shí)應(yīng)用適用于當(dāng)前局部和瞬態(tài)流場(chǎng)狀態(tài)的簡(jiǎn)化機(jī)理可以在保證化學(xué)反應(yīng)模擬準(zhǔn)確性的前提下提高化學(xué)反應(yīng)求解效率.超聲速燃燒中為了準(zhǔn)確捕捉自點(diǎn)火和穩(wěn)焰模態(tài)等流動(dòng)-化學(xué)反應(yīng)強(qiáng)耦合瞬態(tài)特征,一般需要采用盡可能詳細(xì)的化學(xué)反應(yīng)機(jī)理.對(duì)于需要直接積分(DI)的化學(xué)反應(yīng)流求解,其計(jì)算代價(jià)近似與組分?jǐn)?shù)目的平方成正比[58].限于當(dāng)前的計(jì)算技術(shù)水平,目前超聲速燃燒大渦模擬或直接數(shù)值模擬可接受的反應(yīng)組分?jǐn)?shù)目一般最多不超過(guò)50 個(gè)組分[59].過(guò)度簡(jiǎn)化的機(jī)理往往降低了機(jī)理的適用范圍,而適用范圍較廣的機(jī)理尺寸又因尺寸過(guò)大無(wú)法應(yīng)用于高解析度的三維反應(yīng)流模擬.此外,流場(chǎng)中存在復(fù)雜而多變的熱化學(xué)環(huán)境,難以提前針對(duì)所有環(huán)境發(fā)展準(zhǔn)確的簡(jiǎn)化反應(yīng)機(jī)理[56].自適應(yīng)化學(xué)方法通過(guò)凍結(jié)對(duì)目標(biāo)組分和特性關(guān)聯(lián)較小的反應(yīng)路徑和組分實(shí)現(xiàn)了在不同時(shí)間和空間流場(chǎng)中應(yīng)用不同的反應(yīng)機(jī)理,其局部所應(yīng)用的反應(yīng)機(jī)理通常為更詳細(xì)母反應(yīng)機(jī)理的一個(gè)子集.由于局部子機(jī)理的應(yīng)用目標(biāo)熱化學(xué)環(huán)境變窄,因此更容易在較少組分和反應(yīng)步的情況下實(shí)現(xiàn)高化學(xué)反應(yīng)保真度.這里局部子機(jī)理既可以采用當(dāng)前局部分區(qū)狀態(tài)進(jìn)行實(shí)時(shí)簡(jiǎn)化的方式,也可以采用預(yù)估狀態(tài)進(jìn)行預(yù)先簡(jiǎn)化建庫(kù)供后續(xù)查詢(xún)的方式.

自適應(yīng)化學(xué)準(zhǔn)確應(yīng)用的關(guān)鍵是選取合適的分區(qū).以往的研究中一般采用根據(jù)空間坐標(biāo)分區(qū)的方法甚至直接采用并行分區(qū)的方式.這種方式無(wú)法保證當(dāng)前分區(qū)的熱化學(xué)反應(yīng)狀態(tài)是均一的,所簡(jiǎn)化的局部子機(jī)理甚至無(wú)法適用于當(dāng)前分區(qū)內(nèi)所有CFD節(jié)點(diǎn)狀態(tài).對(duì)于采用并行分區(qū)的方法,雖然建模框架搭建較為簡(jiǎn)便,但是同樣無(wú)法保證當(dāng)前分區(qū)內(nèi)熱化學(xué)狀態(tài)的均一性,此外分區(qū)方式(如剖分算法、剖分?jǐn)?shù)目)的改變也會(huì)影響計(jì)算結(jié)果.

本文提出了一種基于當(dāng)?shù)刈赃m應(yīng)分區(qū)的動(dòng)態(tài)自適應(yīng)化學(xué)方法(zonal dynamic adaptive chemsitry,ZDAC).該方法依據(jù)流場(chǎng)的反應(yīng)狀態(tài)變量(如混合分?jǐn)?shù)和溫度等)將計(jì)算區(qū)域劃分為若干分區(qū),每個(gè)分區(qū)內(nèi)的熱化學(xué)反應(yīng)狀態(tài)相對(duì)均一.在本文中選取了混合分?jǐn)?shù)來(lái)表征當(dāng)前組分的變化,用溫度來(lái)表征當(dāng)前分區(qū)的化學(xué)反應(yīng)活躍性.如果包含對(duì)壓力敏感的化學(xué)反應(yīng),可進(jìn)一步拓展加入壓力指標(biāo)體現(xiàn)壓力不均勻性對(duì)反應(yīng)路徑的影響.每個(gè)分區(qū)均對(duì)應(yīng)某一特定的混合分?jǐn)?shù)子空間、較窄的溫度子空間以及壓力子空間,因此可認(rèn)為分區(qū)內(nèi)的熱化學(xué)反應(yīng)狀態(tài)均勻,從而為自適應(yīng)子機(jī)理的匹配性與適用性提供了保障.由于混合分?jǐn)?shù)、溫度和壓力等反應(yīng)狀態(tài)參數(shù)是隨時(shí)間動(dòng)態(tài)變化的,因此Z-DAC 的分區(qū)同DZFM一樣,也是隨流場(chǎng)更新而不斷動(dòng)態(tài)更新的.通過(guò)判斷當(dāng)前CFD 網(wǎng)格單元熱化學(xué)反應(yīng)狀態(tài)參數(shù),自動(dòng)匹配對(duì)應(yīng)Z-DAC 分區(qū).注意,Z-DAC 分區(qū)同樣可以是不規(guī)則形狀,甚至包含了離散的“孤島”形CFD 網(wǎng)格單元.在應(yīng)用中,反應(yīng)狀態(tài)參數(shù)空間的網(wǎng)格劃分在等當(dāng)量比、臨界點(diǎn)/熄火溫度等反應(yīng)活躍性劇變的區(qū)間需要局部加密以提高關(guān)鍵燃燒反應(yīng)區(qū)域流動(dòng)-化學(xué)耦合求解的準(zhǔn)確性.

Z-DAC 中的實(shí)時(shí)機(jī)理簡(jiǎn)化方法采用Lu 和Law[60]提出的直接關(guān)系圖法(DRG),該方法基于組分之間的相互作用關(guān)系圖譜,評(píng)估其他組分對(duì)于某幾類(lèi)關(guān)鍵目標(biāo)組分生成或者消耗速率的影響大小,從而基于給定的誤差閾值去除影響較弱的組分和相關(guān)反應(yīng)步.DRG 法尤其對(duì)于規(guī)模較大的機(jī)理具有高效率和高可靠性.本文采用了基于誤差傳遞的直接關(guān)系圖法(DRGEP)[61],其考慮不同反應(yīng)路徑的誤差傳遞過(guò)程,即考慮包含若干基元反應(yīng)路徑的整體重要性而非單個(gè)基元反應(yīng)的重要性.根據(jù)反應(yīng)過(guò)程選擇合適的目標(biāo)組分,如燃料、氧氣、氮?dú)夂筒糠株P(guān)鍵自由基.根據(jù)基元反應(yīng)關(guān)系圖確定目標(biāo)組分與所有其他相關(guān)組分的反應(yīng)路徑,建立所有目標(biāo)組分的依賴(lài)路徑圖譜.在每條路徑的相鄰組分之間計(jì)算直接相關(guān)系數(shù)(DIC)

其中Nr為當(dāng)前反應(yīng)機(jī)理中所有反應(yīng)的步數(shù),ωi為第i步反應(yīng)的凈反應(yīng)速率,νAi為反應(yīng)i中組分A 的計(jì)量系數(shù),PA和CA分別為組分A 的生成和消耗速率.進(jìn)而累計(jì)乘積路徑上所有相鄰組分之間的DIC得出路徑相關(guān)系數(shù)(PIC).定義總相關(guān)系數(shù)(OIC)為所有反應(yīng)路徑PIC 的最大值.由于OIC 區(qū)分了受檢組分與目標(biāo)物組分的相對(duì)重要性,OIC 低于指定OIC 閾值的組分被視為不重要組分,涉及該組分的反應(yīng)路徑被移除.較之原始DRG 方法,DRGEP 方法中OIC 考慮了包含多步基元反應(yīng)的反應(yīng)路徑的相對(duì)重要性,反應(yīng)路徑上距離目標(biāo)組分越遠(yuǎn)的組分重要性越低.調(diào)節(jié)OIC 閾值可以確定原始機(jī)理中需要去除的不重要組分,直至使給定工況條件下簡(jiǎn)化機(jī)理與原始機(jī)理計(jì)算的點(diǎn)火延遲誤差達(dá)到最大允許偏差.DRGEP 簡(jiǎn)化中用于起始路徑搜索的目標(biāo)組分選擇對(duì)最終反應(yīng)機(jī)理保留路徑影響較大.目標(biāo)組分一般選取為燃料、氧化劑和主要產(chǎn)物(一般為H2O 和CO2).如果有需要與實(shí)驗(yàn)對(duì)比的關(guān)鍵中間自由基組分如OH 和CH2O 等,也需要加入目標(biāo)組分以在機(jī)理簡(jiǎn)化中保留這些中間組分及其關(guān)聯(lián)反應(yīng)路徑.一般而言,當(dāng)非最終產(chǎn)物(H2O 和CO2)中C/H-O 元素當(dāng)量比高于某一特定值(如0.1)時(shí)則認(rèn)為處于較活躍的燃燒反應(yīng)區(qū),此時(shí)需加入CO 和HO2等鏈?zhǔn)疥P(guān)鍵中間組分;當(dāng)反應(yīng)區(qū)溫度高于1000 K 時(shí)則需加入NOx等氮氧化物組分.

1.3 分區(qū)并行自適應(yīng)建表(Z-ISAT)

碳?xì)淙剂系脑敿?xì)反應(yīng)動(dòng)力學(xué)模型往往包含成百上千個(gè)基元反應(yīng),其在實(shí)際的三維燃燒計(jì)算中的應(yīng)用需要巨大的計(jì)算資源.即便采用骨架或簡(jiǎn)化機(jī)理,現(xiàn)有的計(jì)算機(jī)技術(shù)水平仍難以為繼.為了進(jìn)一步提高化學(xué)反應(yīng)源項(xiàng)的計(jì)算速度,局部自適應(yīng)建表(ISAT)方法[62-63]通過(guò)對(duì)局部熱力學(xué)初始和最終狀態(tài)建立映射存儲(chǔ)從而大大節(jié)約了計(jì)算時(shí)間.對(duì)于固定的反應(yīng)機(jī)理,特定時(shí)間步長(zhǎng)之后反應(yīng)熱力學(xué)最終狀態(tài)是初始狀態(tài)的唯一解,因此初始和最終狀態(tài)之間存在一一對(duì)應(yīng)的映射關(guān)系.該方法首先按照常規(guī)方法通過(guò)積分計(jì)算某一化學(xué)熱力學(xué)初始狀態(tài)經(jīng)過(guò)特定反應(yīng)時(shí)間后的最終熱力學(xué)狀態(tài),然后存儲(chǔ)起始狀態(tài)和最終狀態(tài).在后續(xù)的計(jì)算中,通過(guò)在存儲(chǔ)鏈表中對(duì)比當(dāng)前熱力學(xué)狀態(tài)參數(shù)與存儲(chǔ)參數(shù),將相似狀態(tài)對(duì)應(yīng)最終狀態(tài)的插值近似作為當(dāng)?shù)靥囟ǚ磻?yīng)時(shí)間之后的最終狀態(tài),以替代冗長(zhǎng)費(fèi)時(shí)的積分過(guò)程.理論上熱力學(xué)狀態(tài)的映射存儲(chǔ)可以一次完整建立.然而映射表的尺寸隨反應(yīng)組分?jǐn)?shù)成冪次增長(zhǎng),如果完整構(gòu)建熱力學(xué)狀態(tài)參數(shù)空間,即便是較小尺寸的反應(yīng)系統(tǒng)存儲(chǔ)都需要巨大的存儲(chǔ)空間,且后續(xù)的查找效率也較低.事實(shí)上,大多數(shù)反應(yīng)系統(tǒng)僅涉及完整熱力學(xué)狀態(tài)參數(shù)空間的一個(gè)較小子空間.該子空間的范圍取決于實(shí)際燃燒條件,如燃燒反應(yīng)機(jī)理、熱物性參數(shù)和湍流特性等.由于事先無(wú)法確定,必須在模擬過(guò)程中動(dòng)態(tài)建立.隨著映射表的逐步建立,化學(xué)反應(yīng)的計(jì)算速率也逐步提高.對(duì)于準(zhǔn)穩(wěn)態(tài)反應(yīng)流,在存儲(chǔ)空間允許的范圍內(nèi)大部分化學(xué)反應(yīng)都可以通過(guò)查表和插值近似得出,從而極大地加速了復(fù)雜化學(xué)反應(yīng)計(jì)算.

ISAT 映射表的建立過(guò)程如下:

(1) 計(jì)算開(kāi)始時(shí),ISAT 鏈表為空.直接使用剛性反應(yīng)求解器基于反應(yīng)機(jī)理對(duì)初始組分Y0、溫度T0和壓力P0在特定時(shí)間步Δt上進(jìn)行求解積分,得到最終組分Y1和溫度T1(假設(shè)恒壓反應(yīng)).ISAT 樹(shù)狀鏈表中單個(gè)葉節(jié)點(diǎn)存儲(chǔ)信息為:初始熱力學(xué)參數(shù)?0=(Y0,T0,P0),最終熱力學(xué)參數(shù) ?1=(Y1,T1,P1),映射梯度矩陣A=??1/??0和在熱力學(xué)參數(shù)空間內(nèi)原點(diǎn)位于 ?0的精度控制超橢球面(EOA).

(2) 設(shè)有某初始熱力學(xué)參數(shù)空間向量 ?q,0,在ISAT 樹(shù)形鏈表中查詢(xún)與之最接近的節(jié)點(diǎn)設(shè)為 ?0,根據(jù)映射梯度矩陣線(xiàn)性插值計(jì)算近似最終結(jié)果?q,1=.如果 ?q,1位于EOA 范圍內(nèi),則上述線(xiàn)性插值的結(jié)果在制定誤差 εtot允許范圍之內(nèi),所計(jì)算結(jié)果將被取用.否則,則直接積分計(jì)算得到 ?DI,1,確定映射誤差,其中B為縮放矩陣.若 ε ＜εtot,則選用 ?q,1,同時(shí)擴(kuò)張EOA 范圍以包括狀態(tài)向量 ?q,0;否則 ?q,0生成新葉節(jié)點(diǎn).

(3) 計(jì)算中按上述步驟依次查詢(xún)或建立新的葉節(jié)點(diǎn),最終形成一個(gè)二叉樹(shù)形鏈表.當(dāng)新的葉節(jié)點(diǎn)插入時(shí),最相似原始葉節(jié)點(diǎn)的位置將變成一個(gè)二叉結(jié)點(diǎn),重新連接原始葉節(jié)點(diǎn)和新的葉節(jié)點(diǎn).在二叉結(jié)點(diǎn)上將計(jì)算并存儲(chǔ)其所屬兩片葉節(jié)點(diǎn)的超級(jí)分割平面.計(jì)算中將以這個(gè)超級(jí)分割平面作為判據(jù)從樹(shù)形鏈表的根部二叉結(jié)點(diǎn)依次往下層分支二叉結(jié)點(diǎn)上判斷與 ?q,0最相似的葉節(jié)點(diǎn).

初始計(jì)算時(shí),大部分樹(shù)形鏈表的操作都是增加葉節(jié)點(diǎn)(增操作)或擴(kuò)張?jiān)泄?jié)點(diǎn)EOA(擴(kuò)操作).當(dāng)大部分狀態(tài)參數(shù)子空間都被樹(shù)形鏈表的EOA 區(qū)域覆蓋時(shí),讀取操作將更加頻繁.盡管初期的增、擴(kuò)操作會(huì)耗費(fèi)額外的計(jì)算時(shí)間,但隨著鏈表覆蓋范圍的完善,計(jì)算效率將逐步提高.注意,EOA 擴(kuò)操作包含了一定程度的近似,即并不能保證EOA 范圍內(nèi)的所有點(diǎn)的誤差均在指定誤差范圍之內(nèi),這是由于映射關(guān)系的非線(xiàn)性特性.實(shí)際計(jì)算表明EOA 范圍內(nèi)絕大部分狀態(tài)參數(shù)子空間的誤差都在指定的控制誤差范圍之內(nèi)[64-65],并且局部極個(gè)別的誤差越界對(duì)整體計(jì)算精度影響可以忽略,因此ISAT 的EOA 誤差控制方法是符合計(jì)算需求的.

隨著計(jì)算中增操作的進(jìn)行,ISAT 鏈表的尺寸會(huì)越來(lái)越大,尤其是對(duì)于非穩(wěn)態(tài)流場(chǎng).過(guò)大的二叉樹(shù)鏈表尺寸一方面耗費(fèi)存儲(chǔ)空間另一方面增加了二叉樹(shù)搜索深度.由于先前歷史中添加的葉節(jié)點(diǎn)被提取的概率逐漸降低,可以定期對(duì)鏈表進(jìn)行清理維護(hù)以保證每個(gè)葉節(jié)點(diǎn)都有較高的提取率.本研究中實(shí)時(shí)記錄每個(gè)節(jié)點(diǎn)被提取的次數(shù)并動(dòng)態(tài)排名,根據(jù)排名確定最常用葉節(jié)點(diǎn)雙向鏈表.當(dāng)整個(gè)鏈表的尺寸超過(guò)設(shè)定上限后則清空整個(gè)二叉樹(shù)形鏈表,再以最常用雙向鏈表中的節(jié)點(diǎn)信息重新建立二叉樹(shù)形鏈表.計(jì)算中發(fā)現(xiàn),一般最常用鏈表的尺寸上限設(shè)為整個(gè)樹(shù)形鏈表尺寸上限20%～ 50%可保持較好的加速比.

在并行計(jì)算中,一般采取針對(duì)各并行分區(qū)建立獨(dú)立的ISAT 鏈表,相互之間沒(méi)有數(shù)據(jù)交換,即純粹本地建表(purely local processing,PLP) 策略[66].PLP 策略的缺陷在于在強(qiáng)瞬態(tài)、非均勻燃燒流場(chǎng)中無(wú)法保證反應(yīng)計(jì)算載荷均衡性.特別是在超聲速燃燒流場(chǎng)中,不同分區(qū)之間熱化學(xué)反應(yīng)狀態(tài)分布差異極大,并且隨時(shí)間演化而劇烈變化.此時(shí)不同分區(qū)ISAT 表的尺寸和成功取回率均差異較大.一般主要反應(yīng)區(qū)所需維護(hù)的ISAT 鏈表尺寸較大,而無(wú)反應(yīng)區(qū)的ISAT 鏈表尺寸甚至最小僅包含一個(gè)葉節(jié)點(diǎn).此外,當(dāng)并行分區(qū)方式或數(shù)目改變時(shí),原有的ISAT 表需重新建立,無(wú)法重復(fù)利用.

本文提出了基于熱化學(xué)反應(yīng)狀態(tài)進(jìn)行動(dòng)態(tài)分區(qū)的ISAT 建表策略,即Z-ISAT (zonal ISAT).該方法首先依據(jù)混合分?jǐn)?shù)、溫度、壓力、反應(yīng)進(jìn)度變量等當(dāng)?shù)胤磻?yīng)狀態(tài)參數(shù)對(duì)流場(chǎng)進(jìn)行動(dòng)態(tài)分區(qū),并在每個(gè)分區(qū)建立單獨(dú)的ISAT 表.由于每個(gè)分區(qū)的熱化學(xué)反應(yīng)狀態(tài)相對(duì)單一,因此ISAT 表一旦建成即可保持尺寸長(zhǎng)期穩(wěn)定,熱化學(xué)狀態(tài)點(diǎn)覆蓋率和成功取回率效率也可大幅提升,降低了ISAT 鏈表樹(shù)清理和重整頻率.不同于從初始起便固定的并行分區(qū),該基于熱化學(xué)狀態(tài)的動(dòng)態(tài)分區(qū)建表策略尤其適用于強(qiáng)瞬態(tài)和空間非均勻燃燒場(chǎng).計(jì)算中為保持各建表分區(qū)與熱化學(xué)狀態(tài)子范圍空間的嚴(yán)格對(duì)應(yīng)關(guān)系,需定期動(dòng)態(tài)更新分區(qū)包含的CFD 計(jì)算節(jié)點(diǎn).同樣地,ISAT 分區(qū)為不規(guī)則區(qū)域并且可以包含孤立的CFD 節(jié)點(diǎn).對(duì)應(yīng)每種熱力學(xué)狀態(tài)分區(qū)的ISAT 表可以存儲(chǔ)以便在計(jì)算重啟或相似工況計(jì)算時(shí)再利用.

無(wú)通訊PLP-ISAT 策略適用于反應(yīng)區(qū)較集中而穩(wěn)定的燃燒情形.超燃沖壓發(fā)動(dòng)機(jī)中,反應(yīng)區(qū)一般集中于射流尾跡或凹腔,但瞬態(tài)反應(yīng)區(qū)存在劇烈的時(shí)空變化,甚至在特殊的配置條件下還存在振蕩燃燒模態(tài)[67].主反應(yīng)區(qū)在各熱力學(xué)狀態(tài)分區(qū)之間頻繁移動(dòng),意味著相應(yīng)的ISAT 表存在頻繁的節(jié)點(diǎn)插入和增長(zhǎng)操作,一方法增加了ISAT 表的尺寸即內(nèi)存需求,另一方面降低了查表獲取率,加劇了直接積分需求.為此本研究中發(fā)展了基于云計(jì)算的ISAT 并行策略.該云計(jì)算策略將本地?zé)峄瘜W(xué)反應(yīng)狀態(tài)分為頻繁取用和臨時(shí)兩類(lèi)狀態(tài),只有頻繁取用的節(jié)點(diǎn)才存儲(chǔ)ISAT 鏈表,而臨時(shí)狀態(tài)點(diǎn)則通過(guò)打包分配到其余空閑節(jié)點(diǎn)進(jìn)行直接積分計(jì)算求解.這里區(qū)分頻繁取用的判據(jù)為ISAT 中被調(diào)用5 次以上的狀態(tài)點(diǎn).臨時(shí)狀態(tài)點(diǎn)的打包分發(fā)同樣會(huì)占據(jù)內(nèi)存帶寬與通訊時(shí)間,為此研究中發(fā)展了主節(jié)點(diǎn)統(tǒng)計(jì)各從節(jié)點(diǎn)待分配狀態(tài)點(diǎn)數(shù)目,繪制狀態(tài)點(diǎn)重布向量地圖,各從節(jié)點(diǎn)再依據(jù)此重布向量地圖直接采用點(diǎn)對(duì)點(diǎn)通訊交互數(shù)據(jù)的方式,從而避免了集中分發(fā)所帶來(lái)的額外通訊代價(jià).該云計(jì)算策略既保證了ISAT 表較小的尺寸和極高的查詢(xún)獲取率,又通過(guò)并行分發(fā)策略提高了強(qiáng)瞬態(tài)流場(chǎng)中化學(xué)反應(yīng)計(jì)算求解的載荷均衡性.

1.4 控制方程

其中 ν 為依賴(lài)于溫度的運(yùn)動(dòng)黏性系數(shù),可解尺度的應(yīng)變率張量為

其中Sct=1 為湍流施密特?cái)?shù).

動(dòng)態(tài)分區(qū)火焰面模型[70]中求解了四維空間(三維物理空間加一維混合分?jǐn)?shù)空間)中的條件組分輸運(yùn)方程 (17) 而不是傳統(tǒng)的Favre 平均組分輸運(yùn)方程.平均組分質(zhì)量分?jǐn)?shù)通過(guò)對(duì)混合分?jǐn)?shù)空間內(nèi)火焰面變量Qα進(jìn)行概率密度函數(shù)加權(quán)平均得到

其中模型參數(shù)Cvar在亞聲速燃燒中校準(zhǔn)為0.1[73],但超聲速燃燒中該值需進(jìn)一步校準(zhǔn).

1.5 物理模型與求解器細(xì)節(jié)

大部分發(fā)動(dòng)機(jī)燃燒室模擬都必須要考慮壁面的摩擦和冷卻效應(yīng),相應(yīng)地近壁區(qū)域的湍流模擬對(duì)發(fā)動(dòng)機(jī)性能預(yù)測(cè)尤其重要.目前近壁邊界層模擬主要有3 種方法:(1)直接數(shù)值模擬求解所有湍流尺度、(2)近壁湍流模型、(3)壁面函數(shù)方法.在存在慣性子區(qū)的湍流中,直接求解最小湍流尺度的單維度方向網(wǎng)格量需求為雷諾數(shù)的冪次方Re3/4,求解最小耗散時(shí)間步需求為Re1/2[74].然而壁面附近不存在明顯的慣性子區(qū),因此近壁邊界層的模擬對(duì)網(wǎng)格有更加嚴(yán)苛的需求.在保持相同無(wú)量綱近壁距離(=Δynuτ/ν)的情況下,近壁邊界層直接數(shù)值模擬求解所需要的網(wǎng)格數(shù)量為

高馬赫數(shù)發(fā)動(dòng)機(jī)中存在劇烈變化的溫度和壓力,理想氣體模型僅在一定范圍內(nèi)適用.為此物性計(jì)算采用在熱力學(xué)狀態(tài)空間內(nèi)分區(qū)計(jì)算,具體為低于300 K,低于1 kPa 或高于臨界壓力時(shí)密度計(jì)算采用真實(shí)氣體模型,如廣義對(duì)應(yīng)狀態(tài)法則[79];高溫和中等壓力條件下則采用計(jì)算較為簡(jiǎn)便的理想氣體模型.熱擴(kuò)散率和組分?jǐn)U散率計(jì)算分別采用JANAF 格式熱物性數(shù)據(jù)庫(kù)[80]和CHEMKIN 格式輸運(yùn)屬性庫(kù)[81].混合物平均黏性和熱傳導(dǎo)率的計(jì)算分別采用修正版Wilke 定律和摩爾加權(quán)平均[82].

計(jì)算平臺(tái)采用基于OpenFOAM 開(kāi)源框架[83]開(kāi)發(fā)的可壓縮密度基超聲速燃燒求解器Amber (曾用版本名AstroFoam).該求解器采用低耗散修正[84-85]二階Kurganov-Tadmor 中心迎風(fēng)求解格式.求解器已經(jīng)在超聲速燃燒工況模擬中得到廣泛驗(yàn)證[16-18,70,86-89].

圖1 顯示了計(jì)算中流動(dòng)求解和湍流燃燒求解的耦合關(guān)系.其中流動(dòng)模塊Amber 求解式(12)～ 式(16) 基礎(chǔ)N-S 方程以及基于Spalart-Allmaras 的IDDES 湍流模型[90].與PaSR[22]等直接求解三維物理空間內(nèi)的平均組分輸運(yùn)方程不同,DZFM 湍流燃燒模型在四維空間內(nèi)求解分區(qū)火焰面組分輸運(yùn)方程,平均組分則通過(guò)對(duì)條件組分進(jìn)行概率密度函數(shù)加權(quán)平均得到.分區(qū)化學(xué)反應(yīng)速率由條件溫度控制,其通過(guò)統(tǒng)計(jì)方法估算得到,從而避免了復(fù)雜的可壓縮非絕熱條件溫度建模與求解.為體現(xiàn)可壓縮效應(yīng),平均溫度由平均焓和平均組分反推得出.化學(xué)反應(yīng)采用Jachimowski[91]針對(duì)超聲速燃燒開(kāi)發(fā)的13 組分33 步反應(yīng)詳細(xì)氫氣機(jī)理.本研究中整個(gè)計(jì)算域被分割成18 000 個(gè)火焰面分區(qū),其中在主流動(dòng)方向分割為200 個(gè)薄片狀區(qū)域,進(jìn)而對(duì)每個(gè)薄片區(qū)域依據(jù)當(dāng)?shù)鼗旌戏謹(jǐn)?shù)分割為90 個(gè)“年輪”狀環(huán)帶形子區(qū)域.火焰面分區(qū)可以是不規(guī)則、非連續(xù)區(qū)域,其隨混合分?jǐn)?shù)流場(chǎng)而動(dòng)態(tài)更新以保證每個(gè)分區(qū)對(duì)應(yīng)一個(gè)窄域的混合分?jǐn)?shù)子空間.

圖1 流動(dòng)模型與湍流燃燒模型耦合示意圖Fig.1 Coupling diagram of flow model and turbulent combustion model

1.6 測(cè)試工況

本文中選取了高超聲速飛行中兩種典型的發(fā)動(dòng)機(jī)燃燒室構(gòu)型進(jìn)行對(duì)比測(cè)試.該發(fā)動(dòng)機(jī)目標(biāo)飛行馬赫數(shù)10,飛行高度40 km.如圖2(a) 所示,模型發(fā)動(dòng)機(jī)包括進(jìn)氣道、隔離段、燃燒室和尾噴管4 部分.兩種構(gòu)型發(fā)動(dòng)機(jī)的主要區(qū)別在于燃燒室中采用不同的被動(dòng)增混方式,分別為中心支板方式(如圖2(b) 所示)和壁面撐擋方式(如圖2(c) 所示).發(fā)動(dòng)機(jī)的整體長(zhǎng)度為5.61 m,計(jì)算中包括了一部分長(zhǎng)度為0.67 m的外流部分以模擬自由流條件下的進(jìn)氣道特性.進(jìn)氣道長(zhǎng)度2.49 m,采用基于Busemann 基準(zhǔn)流場(chǎng)的流線(xiàn)追蹤法[92]設(shè)計(jì),其中邊界層黏性修正采用參考溫度法.進(jìn)氣道的幾何收縮比為10,靜壓壓縮比約為62.隔離段平滑過(guò)渡進(jìn)氣道出口到圓形超聲速燃燒室,燃燒室包括一段0.985 m 長(zhǎng)的微擴(kuò)張圓筒以過(guò)渡到單邊擴(kuò)張的尾噴管.中心支板燃燒室共計(jì)有24 個(gè)燃料噴孔,其中6 個(gè)燃料噴孔位于連接中心支板的支架側(cè)面上,其余18 個(gè)燃料噴孔以3 個(gè)為一組位于中心支板側(cè)壁.噴孔直徑從0.8～ 1.5 mm 直接變化以適應(yīng)于當(dāng)?shù)貦M向流動(dòng)速度增強(qiáng)射流穿透深度.壁面撐擋超聲速燃燒室同樣采用了24 個(gè)噴孔,2 個(gè)為一組均勻環(huán)繞周向分別從撐擋部件頂部(3 組)和燃燒室壁面(9 組)直接噴注.

圖2 燃燒室中心支板和壁面撐擋結(jié)構(gòu)幾何尺寸Fig.2 Geometric dimensions of the strut and pylon structures

本文中數(shù)值測(cè)試了上述兩種典型構(gòu)型高馬赫數(shù)超燃沖壓發(fā)動(dòng)機(jī)在飛行馬赫數(shù)10 和海拔40 km高度的性能.所模擬飛行工況參數(shù)如表1 所示,其中來(lái)流溫度和壓力均依據(jù)對(duì)應(yīng)高度處大氣層參數(shù).燃料流均為聲速?lài)娮?總溫為室溫,流量保持總包當(dāng)量比為1.0.

表1 測(cè)試飛行工況參數(shù)Table 1 Flight test conditions

計(jì)算域包括外流和發(fā)動(dòng)機(jī)內(nèi)部區(qū)域,其中中心支板發(fā)動(dòng)機(jī)總網(wǎng)格數(shù)量為1.25 億,壁面撐擋發(fā)動(dòng)機(jī)總網(wǎng)格數(shù)量為1.4 億.考慮到復(fù)雜的壁面突起結(jié)構(gòu)和中心支板部件,大部分區(qū)域(96.6%體積)采用無(wú)結(jié)構(gòu)四面體網(wǎng)格剖分,僅在拐角處采用少量楔形(3.39%體積分?jǐn)?shù))、棱柱形和金字塔形(共計(jì)0.01%體積分?jǐn)?shù))填充.平均網(wǎng)格尺度0.8 mm,在燃料噴孔附近局部加密至0.1 mm,近壁邊界層采用15 層向內(nèi)逐漸收縮的棱柱膨脹層,無(wú)量綱近壁距離y+＜1.網(wǎng)格質(zhì)量分析表明99.4% 的網(wǎng)格正交性超過(guò)0.5,99%的網(wǎng)格偏斜度小于0.5.研究針對(duì)中心支板燃燒室開(kāi)展了5413 萬(wàn)、7176 萬(wàn)、1.047 7 億和1.251 億共4 套網(wǎng)格的獨(dú)立性驗(yàn)證,其余3 套較粗網(wǎng)格與最精細(xì)網(wǎng)格的平均壁面壓力誤差分別為1.7%,1.2%和0.8%,隨網(wǎng)格尺度呈指數(shù)遞減,因此可認(rèn)為最精細(xì)網(wǎng)格的計(jì)算結(jié)果具有較好的網(wǎng)格獨(dú)立性,如圖3 所示.

圖3 中心支板燃燒室壁面平均壓力的網(wǎng)格獨(dú)立性分析Fig.3 Grid independence analysis on average wall pressure

2 結(jié)果與討論

2.1 REST 算例驗(yàn)證

基于本研究的數(shù)值模型和方法框架對(duì)國(guó)際上已公開(kāi)發(fā)表的澳大利亞昆士蘭大學(xué)(簡(jiǎn)稱(chēng)UQ)矩形變橢圓高超聲速燃燒室(REST)開(kāi)展了計(jì)算驗(yàn)證.昆士蘭大學(xué)基于T4 反射激波風(fēng)洞 (RST)開(kāi)展了一系列的自由和半自由高馬赫數(shù)超燃沖壓發(fā)動(dòng)機(jī)測(cè)試,在馬赫數(shù) 7.5[6],8[7],8.7[8],10[9],直到12[10-12]均觀(guān)察到了正凈推力.與本研究類(lèi)似,REST 超燃沖壓發(fā)動(dòng)機(jī)的自由射流模擬[15]均包含了部分飛行器前體以模擬內(nèi)外流一體化耦合過(guò)程.驗(yàn)證REST 算例對(duì)應(yīng)飛行馬赫數(shù)12 和海拔36 km,燃料為當(dāng)量比1.0 常溫聲速?lài)娮錃?模擬中采用了1.15 億單元的無(wú)結(jié)構(gòu)網(wǎng)格.如圖4(a)所示,基于當(dāng)前計(jì)算平臺(tái)和數(shù)值框架所預(yù)測(cè)的平均壓力與實(shí)驗(yàn)測(cè)量值吻合較好,同時(shí)捕捉到了與文獻(xiàn)中預(yù)測(cè)類(lèi)似的雙峰結(jié)構(gòu).壓力峰值位于燃燒室噴注點(diǎn)下游的微擴(kuò)張段,同時(shí)進(jìn)氣道噴注點(diǎn)附近也有初始燃燒但產(chǎn)生的壓升弱于弓形激波反射產(chǎn)生的第二峰值壓升.如圖4(b)～ 圖4(e) 所示,計(jì)算中包含了0.5 m 長(zhǎng)的飛行器前體部分以模擬邊界層等外流效應(yīng)對(duì)發(fā)動(dòng)機(jī)內(nèi)流的影響.上游噴注燃料的燃燒較為微弱,從溫度和產(chǎn)物分布看劇烈燃燒反應(yīng)主要發(fā)生于下游噴注點(diǎn)之后的微擴(kuò)張燃燒段.馬赫數(shù)分布表明燃燒室處于超燃沖壓模式,近壁面處特別是燃燒室上側(cè)壁面處的燃燒造成了連續(xù)而狹長(zhǎng)的局部亞聲速區(qū)域.數(shù)值紋影表明下游噴注點(diǎn)之前的進(jìn)氣道和隔離段存在多重反射激波并與剪切層相互作用,而在劇烈燃燒下游因體積膨脹效應(yīng)無(wú)明顯激波存在但可見(jiàn)豐富的湍流漩渦結(jié)構(gòu).

圖4 基于相同計(jì)算框架的REST 馬赫數(shù)12 高超聲速燃燒室模擬Fig.4 Simulation of REST Mach 12 hypersonic combustor based on the same computational framework

2.2 中心支板與壁面撐擋燃燒室性能對(duì)比

圖5 展示了兩種構(gòu)型高超聲速燃燒室的內(nèi)部流場(chǎng)結(jié)構(gòu).從馬赫數(shù)分布來(lái)看,兩種構(gòu)型燃燒室均在支板或撐擋上下游出現(xiàn)了較大的亞聲速區(qū)域,并在支板或撐擋對(duì)側(cè)壁面誘發(fā)了明顯的邊界層分離.由于中心支板結(jié)構(gòu)縮減,加大了流通面積,因此其對(duì)來(lái)流的減速效應(yīng)更加明顯,表現(xiàn)為支板處的超聲速射流核心縮窄至近乎壅塞.壁面撐擋燃燒室由于凹腔的動(dòng)量交換和穩(wěn)焰助燃均進(jìn)一步降低了局部馬赫數(shù),因此其下游的亞聲速區(qū)域延伸更長(zhǎng)距離.從溫度分布來(lái)看,兩種構(gòu)型均出現(xiàn)了較長(zhǎng)區(qū)域的噴注點(diǎn)前部燃燒,主要集中在燃料噴注撐擋一側(cè),表明部分燃料在回流裹挾下被帶到上游并在較高的邊界層滯止溫度(＞1500 K)下被點(diǎn)燃.中心支板燃燒室主要燃燒區(qū)域位于較靠下游區(qū)域,表明燃料混合距離較長(zhǎng).而壁面撐擋燃燒室中凹腔的穩(wěn)焰作用比較明顯,上半部分在凹腔剪切層中出現(xiàn)了火焰,下半部分在凹腔內(nèi)部回流區(qū)實(shí)現(xiàn)了穩(wěn)焰.可見(jiàn)即使對(duì)于當(dāng)?shù)伛R赫數(shù)整體大于2.5 的高超聲速燃燒室,凹腔的穩(wěn)焰作用依然較為顯著.壁面撐擋燃燒室的H2O 含量略高于中心支板回流區(qū),同時(shí)整體2400 K 以上的溫度也略高于整體2000 K 左右的中心支板回流區(qū),表明壁面撐擋回流區(qū)反應(yīng)強(qiáng)度更高.從兩個(gè)構(gòu)型流場(chǎng)對(duì)比來(lái)看,需要在燃燒室設(shè)計(jì)中一方面縮短兩個(gè)噴注點(diǎn)的橫向空間間距,另一方面盡可能讓射流噴注點(diǎn)遠(yuǎn)離外擴(kuò)形壁面,因其流動(dòng)轉(zhuǎn)向會(huì)增加下游尾跡的橫向距離.CEMA 指數(shù)λe是衡量燃燒化學(xué)反應(yīng)劇烈程度的重要指標(biāo)[15,93-94],為化學(xué)反應(yīng)速率雅可比矩陣的最大特征值,通常正值表示爆炸模式,負(fù)值表示趨于平衡的緩燃模式.正值發(fā)生于預(yù)混點(diǎn)火區(qū)域,而負(fù)值發(fā)生于擴(kuò)散控制和后點(diǎn)火區(qū)域.圖中所示為帶符號(hào)對(duì)數(shù)值signλe× lg|λe|,本研究中燃燒室中均為負(fù)值,表明燃燒主要由擴(kuò)散控制,高速流動(dòng)中極短的滯留時(shí)間通常難以維持高溫預(yù)混氣團(tuán)至點(diǎn)火狀態(tài).其中對(duì)應(yīng)溫度高于2000 K 的上游回流區(qū)、支板后主反應(yīng)區(qū)和凹腔等燃燒劇烈處的λe均為108s-1量級(jí)的特征反應(yīng)速率.尾噴管內(nèi)因溫度普遍低于1500 K,即便仍有未完全反應(yīng)的燃燒,反應(yīng)特征速率降為105～ 106s-1.由于反應(yīng)機(jī)理中包含了部分O2,N2和H2O 的離解反應(yīng),因此空氣來(lái)流區(qū)域也存在特征時(shí)間104s-1量級(jí)的反應(yīng).

圖5 中心支板與壁面撐擋高超聲速發(fā)動(dòng)機(jī)內(nèi)部流場(chǎng)Fig.5 The internal flow fields of strut and pylon high-Ma scramjets

圖6 展示了沿程流向各截面提取參數(shù)對(duì)比.壁面撐擋燃燒室靜壓峰值85 kPa,為中心支板燃燒室壓力峰值60 kPa 的1.5 倍.繼進(jìn)氣道壓縮引起的壓升之后,回流區(qū)燃燒進(jìn)一步提高了壓力;其中壁面撐擋燃燒室對(duì)應(yīng)回流區(qū)反應(yīng)更加劇烈因此壓升效應(yīng)更加明顯并產(chǎn)生了一個(gè)較小的峰值.壁面撐擋燃燒室的壓力峰值位于凹腔處,由劇烈燃燒誘導(dǎo);并在其后緩慢下降至尾噴管入口.而中心支板燃燒室的壓力峰值由支板或撐擋誘導(dǎo)的斜激波及其在壁面的反射激波引起,支板下游未見(jiàn)明顯的燃燒誘導(dǎo)峰值;其后燃燒室壓力出現(xiàn)了一個(gè)40 kPa 的平臺(tái)區(qū),燃燒近乎等壓燃燒.與冷態(tài)燃燒室不同,燃燒條件下混合效率的計(jì)算應(yīng)采用元素法[15],即統(tǒng)計(jì)氫元素與氧元素的比例,以考慮已反應(yīng)部分的燃料以及組分間差異擴(kuò)散效應(yīng).回流區(qū)兩種構(gòu)型燃燒室的混合效率接近.說(shuō)明支板或撐擋之后,壁面撐擋燃燒室的混合效率迅速升高遠(yuǎn)超中心支板燃燒室,這是由于壁面撐擋燃燒室的射流穿透深度較高,并且從圖6(c) 凹腔前后燃料濃度分布可知凹腔頂部剪切層和內(nèi)部回流區(qū)也起到了一定程度的增混作用.混合效率在燃燒室下游靠近尾噴管入口處(x=4.2 m)趨于穩(wěn)定值,表明在尾噴管內(nèi)平均馬赫數(shù)5.0 以上的高速流動(dòng)發(fā)生的混合極少.壁面撐擋燃燒室的最終混合效率為85.3%,大幅高于中心支板燃燒室的61.2%.燃燒效率的曲線(xiàn)與混合效率高度相似,最終燃燒效率分別為83.9%和60.7%,僅略低于混合效率,表明高超聲速燃燒室因滯止溫度較高且氫氣的特征反應(yīng)時(shí)間極短(微秒級(jí))存在混合即燃燒的特點(diǎn).這也與CEMA 分析中指出的燃燒多為擴(kuò)散混合控制的觀(guān)察結(jié)果相符.從這個(gè)角度上說(shuō),高超聲速燃燒室提高燃燒效率的關(guān)鍵在于增強(qiáng)混合.總壓計(jì)算中考慮了真實(shí)氣體比熱隨溫度變化的特性,而理想情況下的非等熵關(guān)系式.根據(jù)等熵關(guān)系式積分可得

圖6 沿流向的準(zhǔn)一維性能指標(biāo)Fig.6 Quasi-one-dimensional performance indices along the flow path

其中p和T為當(dāng)前靜壓、靜溫,Tt為總溫,R為氣體常數(shù).高超聲速燃燒室的總壓損失極高,中心支板燃燒室的總壓損失為99.3%,壁面撐擋燃燒室的總壓損失為99%,意味著產(chǎn)生了接近2 倍的熵增.兩種構(gòu)型燃燒室的總壓損失曲線(xiàn)基本一致,其中壁面撐擋附近因局部動(dòng)量交換更為劇烈導(dǎo)致總壓損失更加迅速.在設(shè)計(jì)中,進(jìn)氣道型面微調(diào)為平滑過(guò)渡到燃燒室入口,因此中心支板燃燒室的進(jìn)氣道略短導(dǎo)致總壓損失起始延后,但在保持相同幾何壓縮的條件下進(jìn)氣道末端的總壓損失一致.定義軸向沖量為

其中A為流向面積,ρ為密度,U為流向速度.根據(jù)動(dòng)量守恒,軸向推力F=ΔI,圖6(e)中下降曲線(xiàn)(負(fù)梯度)表示阻力,上升曲線(xiàn)段表示正推力.一個(gè)很明顯的結(jié)果是進(jìn)氣道因壓縮來(lái)流是產(chǎn)生阻力的主要部件;對(duì)應(yīng)x=3.48 m 處,中心支板和壁面撐擋均產(chǎn)生了較大的阻力;壁面撐擋燃燒室x=3.64 m 處的軸向沖量突躍位于凹腔,同時(shí)產(chǎn)生的推力和阻力,總體上表現(xiàn)為阻力;近等直燃燒段產(chǎn)生了微弱阻力;尾噴管是獲得推力的主要部件.從以入口沖量無(wú)量綱化的沖量曲線(xiàn)對(duì)比可見(jiàn),兩種構(gòu)型燃燒室進(jìn)氣道階段產(chǎn)生的阻力和尾噴管產(chǎn)生的推力均比較接近,主要差別在于燃燒段的阻力差別較大:特別是中心支板產(chǎn)生的阻力顯著高于壁面撐擋,凹腔雖然產(chǎn)生了一定程度的阻力但也有效增加了混合和燃燒,其利弊需要進(jìn)一步探索.

表2 對(duì)比了兩種構(gòu)型燃燒室的總體性能.其中壁面撐擋燃燒室的未裝機(jī)凈推力為344 N,比沖比較理想為1234 s;而中心支板燃燒室凈推力較弱僅為99 N,比沖437 s 與常規(guī)火箭發(fā)動(dòng)機(jī)性能接近,體現(xiàn)不出吸氣式發(fā)動(dòng)機(jī)的優(yōu)勢(shì).壁面撐擋燃燒室的燃燒效率也較為理想,高于通常認(rèn)為的可產(chǎn)生推力的80% 準(zhǔn)則[95].壁面撐擋燃燒室的壓比16 也顯著高于中心支板燃燒室的11,由此相比中心支板產(chǎn)生了增幅36% 的無(wú)黏(壓力)推力.同時(shí),中心支板燃燒室的阻力比壁面撐擋燃燒室高23%.中心支板燃燒室進(jìn)氣道外形在保持幾何壓縮比的條件下做了一定的調(diào)整以適應(yīng)微擴(kuò)的隔離段,因此其進(jìn)氣道捕獲流量略低,但計(jì)算中保持總當(dāng)量比恒定.

表2 壁面撐擋和中心支板燃燒室總體性能對(duì)比Table 2 Overall performance comparison between pylon and strut combustors

2.3 動(dòng)態(tài)分區(qū)火焰面特性

如圖7 所示,動(dòng)態(tài)分區(qū)火焰面模型中火焰面結(jié)構(gòu)函數(shù)隨時(shí)間和空間而不斷演化以準(zhǔn)確表征當(dāng)?shù)厝紵瘜W(xué)反應(yīng)狀態(tài).分區(qū)條件平均溫度取決于當(dāng)前分區(qū)的焓增/減程度和反應(yīng)進(jìn)度,用于控制當(dāng)前分區(qū)反應(yīng)速率.隨著流向距離增加,分區(qū)條件平均溫度隨著反應(yīng)趨向平衡遞增,并在燃燒室末端x=4.478 m 處達(dá)到峰值,后續(xù)在尾噴管內(nèi)隨著部分總焓轉(zhuǎn)換為動(dòng)能又逐漸降低.上游的火焰面狀態(tài)會(huì)被下游分區(qū)繼承并繼續(xù)反應(yīng),因此隨著反應(yīng)趨于平衡產(chǎn)物H2O 的濃度持續(xù)增加.自由基O 由O2在高于2500 K 溫度的離解反應(yīng)生成,并為在H2-O2的鏈?zhǔn)椒磻?yīng)中起到重要作用的中間組分 d [O]/dt=k1[ H][O2]-k2[ H2][O],(k1和k2為反應(yīng)速率常數(shù)).一般認(rèn)為燃燒反應(yīng)過(guò)程中O 的濃度在寬溫度范圍內(nèi)相對(duì)恒定[96],離解產(chǎn)生的過(guò)量O 在燃燒段被逐漸消耗,在溫度低于離解溫度(2500 K)僅有燃燒反應(yīng)的尾噴管中保持相對(duì)恒定.火焰面的演化除了在當(dāng)前分區(qū)控制溫度作用下演化以外,還不斷與相鄰分區(qū)通過(guò)對(duì)流交互信息.以O(shè)2濃度為例,其在燃燒段被不斷消耗的同時(shí)也不斷從相鄰未燃來(lái)流分區(qū)卷吸高濃度的O2,因此在燃燒段其濃度甚至?xí)虝荷仙?與H2O 的不斷生成對(duì)應(yīng),O2的濃度在自上游至下游的流動(dòng)過(guò)程中被逐漸消耗,而在尾噴管中因卷吸作用減弱,O2在反應(yīng)中被逐漸消耗.

圖7 不同流向距離處火焰面函數(shù)變化Fig.7 Representative flamelet variations changes at different flow distance

動(dòng)態(tài)分區(qū)火焰面的動(dòng)態(tài)過(guò)程一方面是指火焰面的時(shí)空演化,另一方面火焰面分區(qū)也不斷自適應(yīng)于流場(chǎng)而變化,如圖8(a)所示.首先根據(jù)流向進(jìn)行空間分區(qū)以區(qū)別不同流動(dòng)階段的火焰面.精細(xì)的流向分區(qū)是準(zhǔn)確捕捉點(diǎn)火距離的關(guān)鍵,因其準(zhǔn)確區(qū)分了從純混合、點(diǎn)火到充分燃燒反應(yīng)甚至熄火等不同階段的反應(yīng)狀態(tài).分區(qū)獨(dú)立性分析[70]表明,當(dāng)分區(qū)空間間隔小于一定程度時(shí),流場(chǎng)計(jì)算結(jié)果不再依賴(lài)于分區(qū)數(shù)目.一般在超聲速燃燒室模擬中,一個(gè)合適的分區(qū)準(zhǔn)則為分區(qū)間隔應(yīng)不大于其點(diǎn)火距離,并在燃燒段保證至少50 個(gè)空間分區(qū).即空間上對(duì)流場(chǎng)進(jìn)行分區(qū)后,進(jìn)一步在各“樹(shù)墩”狀剖分段依據(jù)當(dāng)?shù)鼗旌戏謹(jǐn)?shù)進(jìn)行分區(qū).圖8(b) 展示了支板前、中、后典型位置處的混合分?jǐn)?shù)分區(qū)情形.其中在中段距離燃料噴注點(diǎn)附近,當(dāng)?shù)鼗旌戏謹(jǐn)?shù)覆蓋從0 到1 的分布,因此分區(qū)較為密集,而前后段由于擴(kuò)散效應(yīng)僅覆蓋低混合分?jǐn)?shù)段,因此分區(qū)數(shù)目減少.再往上游或下游的更遠(yuǎn)端,當(dāng)?shù)鼗旌衔餇顟B(tài)為完全未預(yù)混的純來(lái)流或已經(jīng)混合較均勻的狀態(tài),此時(shí)混合分?jǐn)?shù)分區(qū)數(shù)目縮減為單個(gè),即單一火焰面即可充分表征當(dāng)?shù)胤磻?yīng)狀態(tài).在依據(jù)混合分?jǐn)?shù)分區(qū)時(shí),為提高主反應(yīng)區(qū)的表征保真度,一般需要在當(dāng)量比混合分?jǐn)?shù)附近加密分區(qū),如圖8(b) 所示.從圖中還可以看到,混合分?jǐn)?shù)分區(qū)形狀不一定是規(guī)則的,甚至可以是分散的孤島集合,在計(jì)算空間對(duì)流時(shí)依據(jù)格林公式沿交界面進(jìn)行積分.當(dāng)流場(chǎng)不均勻性增大時(shí),可進(jìn)一步增加分區(qū)指標(biāo).例如根據(jù)當(dāng)?shù)販囟取ⅠR赫數(shù)等進(jìn)一步進(jìn)行分區(qū)剖分以提高分區(qū)條件均一性.

圖8 火焰面分區(qū)示意圖Fig.8 Diagram of flamelet division

2.4 分區(qū)自適應(yīng)化學(xué)特性統(tǒng)計(jì)

圖9 展示了分區(qū)自適應(yīng)化學(xué)(Z-DAC)中各分區(qū)反應(yīng)機(jī)理的統(tǒng)計(jì).從圖9(a) 中可見(jiàn)前13 步反應(yīng)為在Z-DAC 簡(jiǎn)化中幾乎被始終保留的反應(yīng)步,而14～ 19步反應(yīng)被使用的概率略低,約為90% 左右.20～ 33步反應(yīng)則在63.5% 的分區(qū)機(jī)理簡(jiǎn)化中被忽略.前13 步反應(yīng)為H2的氧化反應(yīng),在燃燒反應(yīng)區(qū)被完整保留,僅在純空氣來(lái)流區(qū)域被忽略.第14～ 19 步涉及過(guò)氧化氫(H2O2)的反應(yīng)僅在OH 濃度較高的活躍反應(yīng)區(qū)才被激活.而后14 步涉及N 元素的反應(yīng)則在1800 K 以上的高溫區(qū)域才會(huì)被激活.從各分區(qū)的反應(yīng)機(jī)理尺寸統(tǒng)計(jì)來(lái)看,64%的區(qū)域采用了31～ 33步反應(yīng),27%的區(qū)域采用了19 步反應(yīng),8%的區(qū)域采用了11～ 13 步反應(yīng).另有極少的區(qū)域采用了7 步反應(yīng),為涉及N2和O2離解反應(yīng)高溫且無(wú)燃料區(qū)域.可見(jiàn)分區(qū)自適應(yīng)化學(xué)方法在將近一半的計(jì)算域上降低了反應(yīng)求解計(jì)算代價(jià),特別是在無(wú)燃料區(qū)反應(yīng)機(jī)理的簡(jiǎn)化幅度更加明顯.

圖9 分區(qū)自適應(yīng)化學(xué)(Z-DAC)統(tǒng)計(jì)數(shù)據(jù)Fig.9 Zonal dynamic adaptive chemistry (Z-DAC) statistics

2.5 與傳統(tǒng)燃燒模型效率對(duì)比

表3 對(duì)比了DZFM 模型與以PaSR 為代表的傳統(tǒng)有限速率模型的計(jì)算效率.PaSR 模型除了需要對(duì)流場(chǎng)中所有CFD 網(wǎng)格單元逐一積分求解當(dāng)?shù)鼗瘜W(xué)反應(yīng)外,還需要求解全組分輸運(yùn)方程,因此計(jì)算代價(jià)一般遠(yuǎn)高于流動(dòng)-化學(xué)反應(yīng)解耦模型.超燃沖壓發(fā)動(dòng)機(jī)中的特征流通時(shí)間一般以毫秒為量級(jí),這里選取1 ms 流動(dòng)物理時(shí)間考察不同燃燒模型及加速方法的計(jì)算效率.在相同并行規(guī)模504 核(CPU 型號(hào)Intel Xeon CPU E5-2640,單核主頻2.40 GHz)的計(jì)算測(cè)試中,純PaSR 模型的單毫秒的計(jì)算周期需要24 d,采用Z-DAC 加速后計(jì)算周期縮短為約20 d,采用ZISAT 加速縮短為14 d.由于Z-DAC 機(jī)理簡(jiǎn)化與ISAT 建表查詢(xún)方法近乎相互獨(dú)立,因此二者耦合使用的加速效率近似于二者單獨(dú)加速比的乘積:332.92 h × (468.16 h/582.61 h)=267.52 h=11 d,僅略小于其耦合使用的12 d.ISAT 與DAC 的耦合使用的組合在部分文獻(xiàn)[97]中被稱(chēng)為建表自適應(yīng)化學(xué)(TDAC),本研究中的方法為基于分區(qū)的Z-TDAC.純DZFM的單毫秒計(jì)算周期縮短為52 h,對(duì)比純PaSR 的加速比高達(dá)11 倍.采用DZFM 模型,每個(gè)火焰面分區(qū)上僅有一維混合分?jǐn)?shù)空間的若干狀態(tài)點(diǎn)需要更新化學(xué)反應(yīng)狀態(tài).Z-DAC 盡管對(duì)當(dāng)前分區(qū)反應(yīng)機(jī)理進(jìn)行了最大程度的簡(jiǎn)化,然而因反應(yīng)狀態(tài)點(diǎn)數(shù)量大幅降低,因此其加速效率不如與傳統(tǒng)有限速率模型耦合時(shí)明顯,僅實(shí)現(xiàn)了5%的加速.同樣地,Z-ISAT 的加速率也為5%.二者耦合使用的Z-TDAC 加速率略高,也僅為8%.表3 列出了不同模型計(jì)算1 ms 物理時(shí)間對(duì)應(yīng)的CPU 時(shí)間.可見(jiàn)分區(qū)解耦湍流燃燒模擬將三維物理空間千萬(wàn)量級(jí)的化學(xué)反應(yīng)狀態(tài)點(diǎn)求解近似為較少的火焰面分區(qū)數(shù)目(本文中為1.8 萬(wàn)個(gè)分區(qū))乘以一維混合分?jǐn)?shù)空間內(nèi)的離散狀態(tài)點(diǎn)上的反應(yīng)求解,從而極大地提升了加速比,同時(shí)使得傳統(tǒng)的反應(yīng)計(jì)算加速方法加速效率相對(duì)不明顯.由于機(jī)理簡(jiǎn)化和建表查詢(xún)均會(huì)產(chǎn)生一定的近似誤差,因此在加速比較低的情況下可以考慮優(yōu)先提高精度.

表3 DZFM 模型與PaSR 模型計(jì)算1 ms 物理時(shí)間對(duì)應(yīng)的CPU 時(shí)間Table 3 CPU time of DZFM and PASR model calculation of 1 millisecond physical time

圖10 基于Borghi 圖統(tǒng)計(jì)分析了以氫氣為燃料的高超聲速燃燒室中各分區(qū)的湍流化學(xué)反應(yīng)交互作用關(guān)系(TCI).其中戴姆克拉數(shù)Da定義為湍流特征時(shí)間 τt與化學(xué)反應(yīng)特征時(shí)間 τc之比,Da=τt/τc.基于大渦模擬數(shù)據(jù)計(jì)算湍流特征時(shí)間時(shí)需要包括可解部分的湍流脈動(dòng)kres[98],即 τt=(kres+kt)/ε .化學(xué)反應(yīng)特征時(shí)間 τc為CEMA 指數(shù)的導(dǎo)數(shù).卡爾洛維奇數(shù)Ka定義為化學(xué)反應(yīng)特征時(shí)間與最小湍流尺度之比Ka=τc/τk,其中Kolmogorov 特征時(shí)間 τk=(ν/ε)1/2.據(jù)此有Re～(τt/τk)2=Da2·Ka2.根據(jù)Ka和Da可將湍流化學(xué)反應(yīng)交互作用關(guān)系分為3 種模式:(1) 火焰面模式Ka＜1 且Da＞10,(2) 薄反應(yīng)區(qū)模式:1 ＜Ka＜100且Da＞10,(3) 慢反應(yīng)模式:Ka＞100 且Da＜10.從兩種燃燒室的統(tǒng)計(jì)分析中可以看出,超過(guò)90% 的絕大部分分區(qū)處于快速反應(yīng)假設(shè)成立的火焰面模式,5%左右的少部分分區(qū)處于薄反應(yīng)區(qū)模式,小于3%的分區(qū)處于慢速反應(yīng)區(qū).氫氣超聲速燃燒與碳?xì)淙剂铣曀偃紵畲蟮膮^(qū)別在于,碳?xì)涑曀偃紵^大部分處于薄反應(yīng)區(qū)(大于50%)和慢速反應(yīng)區(qū)(小于30%)[99].可見(jiàn)即便對(duì)于來(lái)流馬赫數(shù)極高的高超聲速燃燒室,氫氣反應(yīng)速率仍然遠(yuǎn)快于滯留時(shí)間,并且燃燒室中極高的滯止溫度也進(jìn)一步提高了氫氣的反應(yīng)速率,在混合即燃燒的情況下決定燃燒效率的瓶頸在于高效增混.高超聲速燃燒室中雷諾數(shù)高達(dá)108,最小湍流特征尺度為微米量級(jí),這意味直接求解高超聲速燃燒的直接數(shù)值模擬計(jì)算代價(jià)仍然是目前計(jì)算機(jī)技術(shù)無(wú)法承受的,同時(shí)也預(yù)示高超聲速流場(chǎng)中跨越極廣尺度的湍流與速率較慢的碳?xì)淙紵瘜W(xué)反應(yīng)存在極其復(fù)雜的交互作用.

圖10 分區(qū)湍流化學(xué)反應(yīng)交互作用關(guān)系Borghi圖Fig.10 Borghi diagram of zonal turbulent chemical reaction interaction relationship

圖10 分區(qū)湍流化學(xué)反應(yīng)交互作用關(guān)系Borghi 圖(續(xù))Fig.10 Borghi diagram of zonal turbulent chemical reaction interaction relationship (continued)

3 結(jié)論

本文展示了基于動(dòng)態(tài)分區(qū)概念的高馬赫數(shù)全尺寸超燃沖壓發(fā)動(dòng)機(jī)的內(nèi)外流耦合一體化改進(jìn)延遲分離渦(IDDES)模擬研究.研究建立了完整的動(dòng)態(tài)分區(qū)燃燒模擬框架,包括動(dòng)態(tài)分區(qū)火焰面湍流燃燒模型(DZFM),以及分區(qū)自適應(yīng)化學(xué)(Z-DAC)和分區(qū)并行自適應(yīng)建表(Z-ISAT)等化學(xué)反應(yīng)加速模型.前者通過(guò)分區(qū)解耦的思想既準(zhǔn)確表征了當(dāng)?shù)赝牧骰瘜W(xué)交互作用關(guān)系,又有效提升了整場(chǎng)湍流燃燒的計(jì)算效率.后兩者通過(guò)在分區(qū)框架內(nèi)對(duì)化學(xué)反應(yīng)機(jī)理進(jìn)行動(dòng)態(tài)實(shí)時(shí)簡(jiǎn)化和建表查詢(xún),可進(jìn)一步提升當(dāng)前分區(qū)內(nèi)化學(xué)反應(yīng)的求解效率.沿不同空間位置處的火焰面結(jié)構(gòu)函數(shù)的變化表明動(dòng)態(tài)分區(qū)火焰面模型可以捕捉反應(yīng)混合物微團(tuán)內(nèi)溫度和組分等參數(shù)隨反應(yīng)進(jìn)程逐漸演化的物理過(guò)程.精細(xì)的流向分區(qū)準(zhǔn)確區(qū)分了從純混合、點(diǎn)火到充分燃燒反應(yīng)甚至熄火等不同階段的反應(yīng)狀態(tài),因而是準(zhǔn)確捕捉點(diǎn)火距離的關(guān)鍵.混合分?jǐn)?shù)分區(qū)形狀不必是規(guī)則的,甚至可以是分散的孤島集合,其在不同流向距離處的分布密集程度也自適應(yīng)于流場(chǎng)而變化.

首先通過(guò)對(duì)比有完備實(shí)驗(yàn)數(shù)據(jù)的馬赫12 REST標(biāo)準(zhǔn)高超聲速燃燒室模型,驗(yàn)證了分區(qū)模擬框架的保真性.分別基于1.25 和1.4 億網(wǎng)格動(dòng)態(tài)分區(qū)框架對(duì)馬赫數(shù)10 和海拔40 km 條件下全尺寸的中心支板(strut)和壁面撐擋型(pylon)兩類(lèi)構(gòu)型氫氣高超聲速燃燒室進(jìn)行了深入的數(shù)值分析,揭示了高超聲速燃燒中一系列獨(dú)特的規(guī)律和機(jī)理.

壁面撐擋燃燒室的未裝機(jī)凈推力為344 N,比沖比較理想為1234 s;而中心支板燃燒室凈推力較弱僅為99 N,比沖437 s.壁面撐擋燃燒室的燃燒效率也較為理想,高于通常認(rèn)為的可產(chǎn)生推力的80% 準(zhǔn)則.壁面撐擋燃燒室的壓比為16,也顯著高于中心支板燃燒室的壓比11,由此相比中心支板燃燒室產(chǎn)生了增幅為36%的無(wú)黏(壓力)推力.同時(shí),中心支板燃燒室的阻力也比壁面撐擋燃燒室高了23%.兩種構(gòu)型燃燒室的總壓損失均高達(dá)99%,其中壁面撐擋附近因局部動(dòng)量交換更為劇烈導(dǎo)致總壓損失更加迅速.

基于Borghi 圖的數(shù)值分析表明,上述兩類(lèi)燃燒室內(nèi)燃燒的模式均為擴(kuò)散主要控制的火焰面模式.超過(guò)90% 的絕大部分分區(qū)處于快速反應(yīng)假設(shè)成立的火焰面模式,遠(yuǎn)大于碳?xì)涑曀偃紵抑械募s10% 的比例.在混合即燃燒的情況下,提高燃燒效率的關(guān)鍵在于增強(qiáng)混合.壁面撐擋燃燒室中撐擋結(jié)構(gòu)的燃料噴注由于處于亞聲速區(qū)因而其穿透深度較中心支板鈍體表面噴注更深,從而取得了更好的混合效率,進(jìn)而提高了燃燒效率.

分區(qū)自適應(yīng)化學(xué)Z-DAC 在準(zhǔn)確地區(qū)分了各分區(qū)的化學(xué)反應(yīng)狀態(tài),其中H2的氧化反應(yīng)在燃燒反應(yīng)區(qū)被完整保留,涉及過(guò)氧化氫(H2O2) 的反應(yīng)僅在OH 濃度較高的活躍反應(yīng)區(qū)才被激活,而涉及N 元素的反應(yīng)則在1800 K 以上的高溫區(qū)域才會(huì)被激活.64% 的區(qū)域采用了近乎完整的33 步反應(yīng)機(jī)理,27%的區(qū)域采用了19 步反應(yīng),8%的區(qū)域采用了11～13 步反應(yīng),另有極少的高溫純空氣離解區(qū)域采用了7 步反應(yīng)機(jī)理.可見(jiàn)分區(qū)自適應(yīng)化學(xué)方法在將近一半的計(jì)算域上降低了反應(yīng)求解計(jì)算代價(jià),特別是在無(wú)燃料區(qū)反應(yīng)機(jī)理的簡(jiǎn)化幅度更加明顯.

相比于傳統(tǒng)的有限速率PaSR 模型,DZFM 模型將三維物理空間千萬(wàn)量級(jí)的化學(xué)反應(yīng)狀態(tài)點(diǎn)求解近似為較少的火焰面分區(qū)數(shù)目(本研究中為1.8 萬(wàn)個(gè)分區(qū))乘以一維混合分?jǐn)?shù)空間內(nèi)的離散狀態(tài)點(diǎn)上的反應(yīng)求解,實(shí)現(xiàn)了高達(dá)11 倍的加速比,同時(shí)使得傳統(tǒng)的反應(yīng)計(jì)算加速方法加速效率相對(duì)不明顯.由于機(jī)理簡(jiǎn)化和建表查詢(xún)均會(huì)產(chǎn)生一定的近似誤差,因此在加速比較低的情況下可以考慮優(yōu)先提高精度.