化學(xué)非平衡解耦算法精度分析及其改進

劉 君, 董海波, 張文昊

(大連理工大學(xué) 航空航天學(xué)院, 遼寧 大連 116024)

0 引 言

工程熱物理領(lǐng)域認為低溫燃燒技術(shù)可有效解決NOX排放問題，理論上實現(xiàn)這一技術(shù)的關(guān)鍵是控制中間鏈式反應(yīng)，為此國內(nèi)外學(xué)者采用質(zhì)譜儀等先進測試手段對燃燒過程進行細致研究，發(fā)現(xiàn)碳氫燃料的化學(xué)反應(yīng)機理非常復(fù)雜，例如生物柴油包含1034組元4236反應(yīng)[1]。在進行飛行器常規(guī)外流場模擬時如果不考慮湍流模型方程，基于量熱完全氣體假設(shè)的流動控制方程組由5個偏微分方程構(gòu)成，如果對上述生物柴油的化學(xué)非平衡流動進行模擬需要增加1033個組元方程，為解決化學(xué)非平衡流的剛性問題常采用涉及到矩陣求逆的隱式算法，國外文獻比較了計算量與求解變量后得出的結(jié)論是呈3次方關(guān)系[2-3]，由此推算求解變量和方程數(shù)目增加到1038個所引起的計算量增加接近107倍，有理由認為采用精細化學(xué)動力學(xué)模型開展三維化學(xué)非平衡流模擬是比湍流DNS模擬具有更大的挑戰(zhàn)。

目前解決這一難題的主要途徑是在原始測量數(shù)據(jù)基礎(chǔ)上建立簡化反應(yīng)機理。國家科技部2017年度“高性能計算”重點專項指南中關(guān)于“發(fā)動機數(shù)值裝置原型系統(tǒng)—數(shù)值發(fā)動機”的指標是“構(gòu)建基于國產(chǎn)航空燃料化學(xué)反應(yīng)簡化機理的燃燒模型，簡化機理物種數(shù)在35個左右，模擬點火延遲預(yù)測誤差不大于30%，燃燒模型能夠比較準確模擬流動與化學(xué)反應(yīng)相互作用”。已有的研究表明采用不同簡化準則得到的模型存在著一定的迥異，其明顯特征表現(xiàn)為預(yù)測點火延遲特性有較大差異[1]。以反應(yīng)機理相對簡單的H2和Air燃燒為例，通過對國內(nèi)外文獻調(diào)研發(fā)現(xiàn)，采用忽略中間鏈式反應(yīng)的3組元2反應(yīng)的單向總包反應(yīng)模型[4]，即使采取網(wǎng)格加密、時間小步長等措施提高計算精度，也無法模擬出與試驗相符合的點火延遲、周期性振蕩燃燒等物理現(xiàn)象，得到點火延遲大致結(jié)構(gòu)需要采用7組元8反應(yīng)[5]機理，模擬出周期性振蕩燃燒需要9組元19反應(yīng)[2]或13組元33反應(yīng)[6-7]機理。由此看來，采用35個左右組元的簡化機理實現(xiàn)誤差不大于30%的點火延遲特性預(yù)測確實是很高的指標要求。反應(yīng)簡化機理缺乏普適性，如果燃料的產(chǎn)地和生產(chǎn)商改變需要重新研究。如果能夠提高化學(xué)非平衡流問題的計算效率，從而可以考慮更為復(fù)雜的化學(xué)反應(yīng)機理，以此來解決研制數(shù)值發(fā)動機遇到的預(yù)測精度和計算量的難題，也是一個很好的技術(shù)途徑。本文主要介紹化學(xué)非平衡流問題計算方法的研究進展。

1 化學(xué)非平衡流解耦算法理論基礎(chǔ)

20世紀60 年代，計算流體力學(xué)主要以顯式計算格式為主，為了解決化學(xué)反應(yīng)源項引起的剛性問題，采用解耦算法在空間點上把控制方程組分解為流動和化學(xué)反應(yīng)兩部分，為滿足反應(yīng)算子的計算穩(wěn)定性條件需要迭代上千步才能匹配流動算子時間迭代一步，由于其計算效率偏低，難以滿足工程需要。20世紀70 年代，在解決航天器再入過程中遇到的熱障問題推動下，化學(xué)非平衡流計算方法取得較大進展，構(gòu)造出反應(yīng)源項線性化的隱式算法，包括流動算子也采用隱式方法計算的全隱算法和流動算子采用顯式方法計算的點隱算法。模擬點火延遲特性對不同算法存在較高的時間精度要求，全隱算法只能結(jié)合雙時間步法來模擬非定常流動，由此增加的計算量難以接受。通過下面的分析可以看出點隱算法實現(xiàn)時間二階精度也較為困難。

對于波動模型方程：

采用顯式差分格式離散，時間二階精度的Taylor展開式：

利用方程(1)把二階時間導(dǎo)數(shù)轉(zhuǎn)化為空間導(dǎo)數(shù)：

化學(xué)非平衡流模型增加了源項：

點隱算法把源項在空間點線化：

源項進行線化處理以后具有時間二階精度，將方程(6)帶入到方程(4)得到點隱算法的計算格式：

因此，計算格式(7)達不到時間二階精度，還沒有看到按照方程(8)構(gòu)造高階格式需要處理上式中源項f(u)的時間和空間導(dǎo)數(shù)的算法。

采用解耦算法計算方程(5)本質(zhì)是分步求解如下波動方程和常微分方程組：

常微分方程有多種成熟計算方法，采用算子形式表示為：

un+1=1+Lf(Δt)un(10)

不管算子Lf具體形式，時間二階精度離散基礎(chǔ)依然為：

波動方程的時間二階精度離散算子如前，總算子按照如下Strang形式：

上式展開以后會得到方程組(9)的時間和空間導(dǎo)數(shù)，可以證明對于線形算子具有時間二階精度。因此，解耦算法易于構(gòu)造時間二階格式。

2 劉君解耦算法的優(yōu)化

包含有n個組元的化學(xué)非平衡流守恒型Euler方程組：

式中守恒變量和源項表達式：

U=(ρ,ρu,ρv,ρw,E,ρ1,…,ρn-1)T(14)

S=(0,…,0,σ1,…,σn-1)T(15)

然后用E′代替E對原始方程(13)進行改造，得到新的守恒變量和源項表達式為：

U0=(ρ,ρu,ρv,ρw,E′,ρ1,…,ρn-1)T=(U1,U2) (18)

進一步引入符號γ′把E′整理成如下形式：

新型解耦算法早期只能用于有限差分法，先后采用過VanLeer、HLLC、Roe、AUSM等格式進行對流項通量分裂，空間離散采用過NND2MWAF、ENO、五階WENO等格式[12-14]，時間離散包括二階顯式方法和隱式雙時間法[15]。近期課題組結(jié)合解耦算法的特點和有限體積法的空間平均特性提出一種優(yōu)化解耦算法，如圖1所示，流場從n時刻更新到n+1時刻，控制體k和相鄰單元之間的質(zhì)量交換通過求解流動算子Lf體現(xiàn)，流體微團在對流或者擴散作用下進出控制體過程中沒有發(fā)生化學(xué)反應(yīng)，組元質(zhì)量分數(shù)ci保持不變，通過記錄Δt時間步內(nèi)進出邊界的質(zhì)量流量，以此計算出n+1時刻控制體中組元i的密度值。假設(shè)在Δt時間步內(nèi)，流體微團從控制體單元k流出到單元3，從相鄰單元1和2流入。由于U1和U2方程解耦，可以單獨求解U1的方程，記錄進出邊界的質(zhì)量流量，例如，單元1外法向速度un，流量為m1=(ρunSΔt)1，S邊界面積，那么控制體k在n+1時刻i組元的密度為：

通過以上簡化處理在求解過程中不需要計算U2的偏微分方程組就可以得到U2中組元的密度。即無論多少組元，流動算子只求解基本變量構(gòu)成的、形式上與量熱完全氣體完全相似的5個方程組，目前優(yōu)化算法主要完成無粘流動方程的驗證，理論上通過記錄相鄰控制體單元邊界面組元擴散引起的流量實現(xiàn)粘性流動計算模擬，下一步準備開展相關(guān)編程和驗證工作。

圖1 優(yōu)化算法機理示意圖Fig.1 Schematic of optimization algorithm mechanism

本文作者在文獻[11]中采用優(yōu)化解耦算法對1972年Lehr[16]完成的H2/Air預(yù)混氣體條件下激波誘導(dǎo)燃燒試驗進行模擬，使用修正的Jachimowski(9組元19反應(yīng))反應(yīng)機理[17]，采用二維模型進行模擬，計算的工況包括馬赫數(shù)M=6.46的超爆轟模態(tài)和M=4.48、4.79兩種跨爆轟模態(tài)，比較計算得到的流場結(jié)構(gòu)與試驗數(shù)據(jù)符合很好，沿駐點流線的參數(shù)分布和文獻基本一致。

3 反應(yīng)算子的改進算法

從方程(13)出發(fā)，描述化學(xué)反應(yīng)的方程組可以表示為：

結(jié)合定容假設(shè)條件?ρ/?t=0可將方程(24)轉(zhuǎn)化為：

最早采用梯形公式對反應(yīng)源項常微分方程組進行求解，應(yīng)用中發(fā)現(xiàn)存在初場敏感性和收斂不穩(wěn)定的缺陷，需要采取在穩(wěn)定無反應(yīng)流場基礎(chǔ)再耦合反應(yīng)算子Lr模塊、反應(yīng)劇烈需要小步長追趕等經(jīng)驗措施，后來改用VODE(常微分變系數(shù)方法)[18]和α-QSS(擬穩(wěn)態(tài)逼近方法)[19]兩種求解器解決了穩(wěn)定性問題，但是計算效率明顯降低。精細積分方法是由鐘萬勰院士[20-21]在1994年提出的，在常微分方程的求解上具有突出的優(yōu)勢：(1) 具有極高的計算精度，(2) 具有較好的計算魯棒性。本文探索采用精細積分方法求解反應(yīng)算子。

可以把方程(23)和(25)寫成同一矩陣形式：

其中X=ρ1,…,ρn,TT，H是系數(shù)矩陣，已知X(tk)時間推進一步tk+1=tk+Δt后：

Xk+1=exp(HΔt)·Xk(27)

為了提高Y=exp(HΔt)的計算精度，根據(jù)指數(shù)矩陣函數(shù)特性exp(HΔt)≡expHΔt/mm，將Δt分成非常小的時間區(qū)段τ=Δt/m，其中m為任意正整數(shù)，可選用m=2N。由于下式的Ya與單位矩陣I相比非常小，容易被計算機舍入操作：

exp(Hτ) ≈I+Hτ+(Hτ)2·

≈I+Ya(28)

精細積分算法先將Y做分解：

Y=(I+Ya)2N=(I+Ya)2N-1×(I+Ya)2N-1=… (29)

當(dāng)N次循環(huán)結(jié)束后，再進行Y=I+Ya運算，這樣可以保證計算精度。精細積分早期應(yīng)用于結(jié)構(gòu)動力學(xué)分析[22]中，由于其數(shù)值結(jié)果的高度精確，已經(jīng)在優(yōu)化控制[23]、偏微分方程的精細求解[24]、非線性動力系統(tǒng)剛性方程求解[25]等領(lǐng)域得到了廣泛應(yīng)用。

采用Jachimowski的9組元19反應(yīng)機理，計算初始條件為壓強P=101.325 kPa、溫度T=1200 K的預(yù)混燃氣爆轟過程。廣泛使用的一種剛性常微分方程求解器是VODE，采用BDF(Backward Differentiation Formula)方法求解，具有高階精度，常被用來驗證其它方法的精度。采用統(tǒng)一時間步長Δt=1×10-8s進行計算，兩種方法得到的溫度和組元質(zhì)量分數(shù)隨時間變化如圖2和圖3所示，幾乎一致，說明采用精細積分方法求解化學(xué)動力學(xué)方程組具有較高的精度。

圖2 溫度分布時程曲線Fig.2 Time history curve of temperature distribution

VODE是一種隱式方法，求解過程涉及到矩陣運算，在化學(xué)組元較多、矩陣維數(shù)較大時，會導(dǎo)致計算量較大，難以用于復(fù)雜反應(yīng)機理的化學(xué)非平衡流動模擬[26-27]。α-QSS是一種擬穩(wěn)態(tài)逼近方法，對于線形問題具有二階精度，對于非線性問題因源項局部線性化處理精度達不到二階，理論上通過基于精度自動選擇時間步長可以保持計算過程穩(wěn)定，但是在應(yīng)用中發(fā)現(xiàn)在高壓環(huán)境或者沒有稀釋氣體的情況下，α-QSS對于部分燃燒問題的求解存在不穩(wěn)定現(xiàn)象[28]。精細積分作為一種顯示方法，計算精度高，穩(wěn)定性好，在選取較大的時間步長時也可以得到較高的計算精度。為此，選取較大步長Δt=5×10-7s，比較三種方法的計算效率和精度，CPU計算時間和相對誤差如表1所示，相對誤差是指計算的穩(wěn)定溫度與參考值之間的誤差，取VODE在時間步長Δt=1×10-9s下計算的穩(wěn)定溫度作為參考值。與VODE和α-QSS兩種方法相比，精細積分在大時間步長情況下計算用時最短。此外，由圖4可看出精細積分計算結(jié)果與前兩種方法幾乎一致，說明精細積分方法在大步長的情況下，具有較高的計算效率。

圖3 H原子和O原子質(zhì)量分數(shù)時程曲線Fig.3 Time history curves of mass fraction of H and O

表1 不同計算方法在相同步長條件下的結(jié)果對比Table 1 Comparision of the results of different calculation methods under the condition of the same step size

圖4 不同計算方法在Δt=5×10-7s下的溫度分布時程曲線Fig.4 Time history curve of temperature distribution with different calculation methods at Δt=5×10-7s

結(jié)合流動算子模擬Lehr試驗，對于馬赫數(shù)Ma=6.46的超爆轟模態(tài)和Ma=4.79跨爆轟模態(tài)，精細積分方法的計算結(jié)果與文獻[11]中采用α-QSS的計算結(jié)果完全一致，在此不再展示。下面給出馬赫數(shù)Ma=3.55亞爆轟模態(tài)的計算結(jié)果。彈丸以1892 m/s飛行速度進入H2/O2(2H2+O2)的預(yù)混氣體，當(dāng)?shù)丨h(huán)境溫度和壓力分別為T=293K，p=24793Pa，相應(yīng)的爆轟速度為2550m/s。采用250×200的四邊形網(wǎng)格、Evan(7 組元8反應(yīng))機理[29]模擬。密度云圖中激波陣面和燃燒陣面位置與試驗結(jié)果比較如圖5所示，過激波后形成的誘導(dǎo)區(qū)和燃燒陣面等結(jié)構(gòu)符合較好，在二維情況下，控制方程組不包含w，因此圖6所示為求解4個偏微分方程的優(yōu)化算法和求解10個偏微分方程的原算法得到的密度云圖對比，從以上結(jié)果的對比中可以驗證優(yōu)化算法對于穩(wěn)定爆轟的流場具有較好的空間分辨率。

圖5 計算結(jié)果與試驗結(jié)果對比Fig.5 Comparison between the simulation and experimental results

圖6 亞爆轟密度云圖比較Fig.6 Comparison of density contour of the subdetonation case

誘導(dǎo)區(qū)長度對應(yīng)于激波波后的點火延遲時間。經(jīng)過誘導(dǎo)區(qū)以后，駐點流線上的壓力幾乎保持恒定，而溫度開始上升，密度下降。這種現(xiàn)象也可以清晰地從主要組元在駐點流線上的分布曲線中得到，如圖7所示，相應(yīng)的計算結(jié)果在文獻[30]也有所提及。從圖中可以看出，在誘導(dǎo)區(qū)相對應(yīng)的范圍內(nèi)，H2O組元質(zhì)量分數(shù)保持為0，說明沒有明顯的H2O產(chǎn)生；經(jīng)過燃燒陣面以后，H2O組元的質(zhì)量分數(shù)逐漸增加伴隨著O2和H2組元質(zhì)量分數(shù)迅速降低，這種現(xiàn)象說明經(jīng)過燃燒陣面以后急劇上升的溫度促使化學(xué)反應(yīng)更加劇烈，最終產(chǎn)生大量的H2O組元。

采用精細積分方法與α-QSS方法的CPU計算時間對比如圖8所示。在相同迭代步數(shù)條件下，流動算子的CPU計算時間幾乎保持一致；當(dāng)N=20時精細積分方法所用的CPU計算時間最長，要高于α-QSS方法；隨著N取值的逐漸降低，計算效率逐漸提高。表2列出了反應(yīng)算子在不同迭代步數(shù)條件下的CPU時間對比，當(dāng)N=10時精細積分方法與α-QSS方法的CPU計算用時非常接近，當(dāng)N=5時精細積分方法的計算用時只相當(dāng)于α-QSS方法的79%左右，說明精細積分方法可以進一步提高計算效率，同時針對剛性問題的求解，由精細積分方法構(gòu)造的反應(yīng)算子求解器具有更好的穩(wěn)定性。

圖7 組元質(zhì)量分數(shù)沿駐點流線分布Fig.7 Distribution of mass fraction along the stagnation streamline

圖8 計算步數(shù)與CPU時間對比曲線Fig.8 Calculated step number and CPU time curve

表2 反應(yīng)算子CPU時間對比Table 2 Comparison of CPU time

4 結(jié) 論

本文從求解化學(xué)非平衡流問題的難點入手，結(jié)合劉君解耦算法和有限體積法的平均特性，提出一種新型的優(yōu)化算法，針對經(jīng)典的激波誘導(dǎo)燃燒算例進行數(shù)值模擬，通過與不同文獻和試驗的結(jié)果相對比，驗證了該優(yōu)化算法的時空精度。同時，基于精細積分方法構(gòu)造了一種剛性問題求解器，其模擬結(jié)果與傳統(tǒng)的VODE方法和α-QSS方法的計算結(jié)果對比分析，適當(dāng)?shù)倪x取時間步長，可以充分發(fā)揮精細積分方法的計算優(yōu)勢。