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

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

* (中國科學(xué)院力學(xué)研究所高溫氣體動力學(xué)國家重點實驗室,北京 100190)

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

引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

ISAT 映射表的建立過程如下:

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

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

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

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

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

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

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

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

1.4 控制方程

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

其中Sct=1 為湍流施密特數(shù).

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

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

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

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

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

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

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

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

1.6 測試工況

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

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

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

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

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

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

2 結(jié)果與討論

2.1 REST 算例驗證

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖10 基于Borghi 圖統(tǒng)計分析了以氫氣為燃料的高超聲速燃燒室中各分區(qū)的湍流化學(xué)反應(yīng)交互作用關(guān)系(TCI).其中戴姆克拉數(shù)Da定義為湍流特征時間 τt與化學(xué)反應(yīng)特征時間 τc之比,Da=τt/τc.基于大渦模擬數(shù)據(jù)計算湍流特征時間時需要包括可解部分的湍流脈動kres[98],即 τt=(kres+kt)/ε .化學(xué)反應(yīng)特征時間 τc為CEMA 指數(shù)的導(dǎo)數(shù).卡爾洛維奇數(shù)Ka定義為化學(xué)反應(yīng)特征時間與最小湍流尺度之比Ka=τc/τk,其中Kolmogorov 特征時間 τ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)計分析中可以看出,超過90% 的絕大部分分區(qū)處于快速反應(yīng)假設(shè)成立的火焰面模式,5%左右的少部分分區(qū)處于薄反應(yīng)區(qū)模式,小于3%的分區(qū)處于慢速反應(yīng)區(qū).氫氣超聲速燃燒與碳氫燃料超聲速燃燒最大的區(qū)別在于,碳氫超聲速燃燒絕大部分處于薄反應(yīng)區(qū)(大于50%)和慢速反應(yīng)區(qū)(小于30%)[99].可見即便對于來流馬赫數(shù)極高的高超聲速燃燒室,氫氣反應(yīng)速率仍然遠快于滯留時間,并且燃燒室中極高的滯止溫度也進一步提高了氫氣的反應(yīng)速率,在混合即燃燒的情況下決定燃燒效率的瓶頸在于高效增混.高超聲速燃燒室中雷諾數(shù)高達108,最小湍流特征尺度為微米量級,這意味直接求解高超聲速燃燒的直接數(shù)值模擬計算代價仍然是目前計算機技術(shù)無法承受的,同時也預(yù)示高超聲速流場中跨越極廣尺度的湍流與速率較慢的碳氫燃燒化學(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é)論

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

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

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

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

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

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