基于解耦算法的化學(xué)非平衡流并行數(shù)值模擬

田壯壯，王光學(xué),2，張懷寶，張 帆

(1. 中山大學(xué)物理學(xué)院，廣州 510006；2. 國防科技大學(xué)空天科學(xué)學(xué)院，長沙 410073；3. 中山大學(xué)航空航天學(xué)院，廣州 510006)

0 引言

解耦算法又稱為時(shí)間分裂算法或算子分裂法[1]。在此類方法的計(jì)算中，半離散化得到的常微分方程組被分解為剛性和非剛性兩部分，而后進(jìn)行獨(dú)立求解。劉君等[2-3]構(gòu)造了一種化學(xué)非平衡流解耦算法，這種算法將化學(xué)非平衡流的內(nèi)能中與溫度無關(guān)的能量，即有效零點(diǎn)能或化學(xué)焓分離出來，添加到化學(xué)反應(yīng)源項(xiàng)中。進(jìn)而引入等效比熱比，構(gòu)造了在形式上與量熱完全氣體完全相同的能量與壓力的關(guān)系式，并用等效比熱比計(jì)算等效聲速。進(jìn)行時(shí)間分裂后，物理問題被分解成流動(dòng)(非剛性)和化學(xué)反應(yīng)(剛性)兩部分，流動(dòng)控制方程的求解與量熱完全氣體的求解方法完全相同，而且流動(dòng)方程與化學(xué)反應(yīng)方程都具有明確的物理意義。

利用這一解耦算法，劉君[3]對(duì)沖壓加速器中的化學(xué)非平衡流進(jìn)行了數(shù)值模擬。其中，流場(chǎng)空間離散采用ENO格式[4]，化學(xué)反應(yīng)采用梯形公式計(jì)算；并分別采用7組分8反應(yīng)的氫氣/空氣反應(yīng)機(jī)理和13組分32反應(yīng)的甲烷/空氣反應(yīng)機(jī)理描述化學(xué)反應(yīng)過程；隨后又采用19組分65反應(yīng)甲烷/空氣反應(yīng)機(jī)理進(jìn)行了模擬[5]?；谕瑯拥慕怦罘椒?，在國內(nèi)首先對(duì)鈍頭體激波誘導(dǎo)振蕩燃燒現(xiàn)象進(jìn)行了數(shù)值模擬，得到了準(zhǔn)確的非定常振蕩燃燒現(xiàn)象[6]；在化學(xué)動(dòng)力學(xué)模型對(duì)氫氣/空氣超聲速燃燒模擬的影響方面進(jìn)行了研究[7]。應(yīng)用上述解耦算法，孫明波等[8]采用5階WENO格式[9]進(jìn)行空間離散，對(duì)鈍頭體激波誘導(dǎo)振蕩燃燒進(jìn)行了模擬。

通過采用上述的解耦算法，化學(xué)非平衡流求解器對(duì)流動(dòng)控制方程(非剛性)和化學(xué)反應(yīng)源項(xiàng)(剛性)分別采用合適的方法求解，流動(dòng)方程求解不受化學(xué)反應(yīng)時(shí)間尺度的影響，從而在保證計(jì)算精度的前提下提高了計(jì)算效率；但是在復(fù)雜問題的數(shù)值計(jì)算中，化學(xué)反應(yīng)源項(xiàng)的存在導(dǎo)致整體計(jì)算量顯著大于量熱完全氣體數(shù)值模擬的計(jì)算量。本文對(duì)化學(xué)非平衡流解耦算法的非結(jié)構(gòu)網(wǎng)格有限體積法求解器進(jìn)行MPI并行化，以便為更復(fù)雜的化學(xué)反應(yīng)流動(dòng)數(shù)值模擬提供工具。

1 控制方程和解耦算法

本文計(jì)算中求解的微分形式二維軸對(duì)稱歐拉方程可以寫為

(1)

其中

式中，Q表示守恒變量，S和Qf表示反應(yīng)源項(xiàng)和軸對(duì)稱方程源項(xiàng)；E和F分別表示x和y兩個(gè)方向的對(duì)流通量；σi表示第i個(gè)組分的化學(xué)反應(yīng)速率。此外，u和v分別為x和y方向的速度，p為壓力，ρ為混合氣體的總密度，ρi為第i個(gè)組分的密度。單位質(zhì)量總能是

財(cái)務(wù)會(huì)計(jì)與管理會(huì)計(jì)有著緊密的聯(lián)系，工作中也存在相似性，推動(dòng)轉(zhuǎn)型并不是對(duì)財(cái)務(wù)會(huì)計(jì)進(jìn)行否定，而是將財(cái)務(wù)會(huì)計(jì)以及管理會(huì)計(jì)的優(yōu)點(diǎn)進(jìn)行融合。在保證企業(yè)會(huì)計(jì)信息公開、透明的同時(shí)，為企業(yè)戰(zhàn)略的制定提供強(qiáng)有力的科學(xué)依據(jù)。因此應(yīng)當(dāng)對(duì)使用的財(cái)務(wù)資料格式進(jìn)行適當(dāng)?shù)脑O(shè)計(jì)和調(diào)整。

引言中介紹的解耦算法定義等效內(nèi)能為

(2)

進(jìn)而可將能量方程改寫為

(3)

基于同樣的關(guān)系，還可以定義等效比熱比

(4)

最終將控制方程改寫為

(5)

其中

基于時(shí)間分裂法，可以將上述基于等效內(nèi)能的控制方程分裂為流動(dòng)方程和化學(xué)反應(yīng)方程

(6)

(7)

對(duì)上述兩個(gè)方程分別進(jìn)行求解，并且依據(jù)Strang分裂格式

Q′,n+1=Lc(Δt/2)Lf(Δt)Lc(Δt/2)Q′,n

(8)

的順序進(jìn)行計(jì)算。其中，Lc和Lf分別為化學(xué)反應(yīng)算子和流動(dòng)算子，n表示第n個(gè)時(shí)間步。如引言所述，流動(dòng)算子可以是量熱完全氣體流動(dòng)求解中采用的各類2階時(shí)間精度離散方法(顯式[4]或隱式[8])，而化學(xué)反應(yīng)算子可以采用各類剛性常微分方程組求解算法。兩類算子通過式(8)所示的順序進(jìn)行計(jì)算，可以確?；瘜W(xué)非平衡流求解的整體二階時(shí)間精度。

2 并行計(jì)算

在數(shù)值模擬中，區(qū)域分解法是一種常用的并行算法。其基本思想是：將整個(gè)計(jì)算區(qū)域分成N個(gè)子域并相應(yīng)地分配給N個(gè)計(jì)算進(jìn)程，每個(gè)進(jìn)程計(jì)算各自子域的流場(chǎng)；計(jì)算時(shí)各個(gè)進(jìn)程之間進(jìn)行必要的數(shù)據(jù)通信。為提高區(qū)域分解法的效率，需要盡量使各個(gè)進(jìn)程的負(fù)載均衡，并減少進(jìn)程間的通信。本文使用常用的METIS分區(qū)軟件來實(shí)現(xiàn)區(qū)域分解，獲得大小均衡的計(jì)算子域，從而保證分區(qū)間的負(fù)載均衡；同時(shí)，要使分區(qū)間的連接盡可能簡(jiǎn)單，從而減少分區(qū)間的通信。

圖1給出了網(wǎng)格分區(qū)的一個(gè)例子。方形區(qū)域內(nèi)的三角形網(wǎng)格被分為A、B兩個(gè)分區(qū)，如圖1(a)所示。初始分區(qū)A、B中沒有重合的單元。由于使用空間2階格式進(jìn)行計(jì)算，在計(jì)算一個(gè)單元的通量時(shí)，需要用到直接相鄰單元的物理量及其梯度；而計(jì)算直接相鄰單元的物理量梯度時(shí)需要第二層相鄰單元的物理量。這里稱兩層相鄰單元為計(jì)算分區(qū)的“外部單元”，稱分區(qū)原有單元為“內(nèi)部單元”，分別如圖1(b)、(c)所示。其中，標(biāo)記為1的是直接相鄰單元，標(biāo)記為2的是第二層相鄰單元。為了減少計(jì)算過程中的通信次數(shù)，需要存儲(chǔ)外部單元的物理量。因此，各個(gè)進(jìn)程存儲(chǔ)的子域網(wǎng)格是分區(qū)內(nèi)部單元和外部單元的總和，分別如圖1(d)、(e)所示。

(a)初始分區(qū)

(b)A分區(qū)的外部單元

(c)B分區(qū)的外部單元

(d)最終的A分區(qū)

(e)最終的B分區(qū)圖1 網(wǎng)格分區(qū)示例Fig.1 Schematic of grid partitioning

從圖1中容易看出，對(duì)于多個(gè)計(jì)算分區(qū)，一個(gè)分區(qū)的外部單元必定是另一個(gè)分區(qū)的內(nèi)部單元，而一個(gè)分區(qū)的內(nèi)部單元可能對(duì)應(yīng)著多個(gè)分區(qū)的外部單元。計(jì)算中，在每一步時(shí)間推進(jìn)后，外部單元上的物理量需要從所對(duì)應(yīng)的內(nèi)部單元獲得，以保證不同子網(wǎng)格間流場(chǎng)物理量的一致性。因此，在進(jìn)行網(wǎng)格分區(qū)時(shí)，需要建立所有分區(qū)間的內(nèi)部單元和外部單元的映射關(guān)系。

各個(gè)進(jìn)程都額外存儲(chǔ)了外部單元網(wǎng)格，因此，不同子域網(wǎng)格間是相互交錯(cuò)重疊的，各個(gè)子網(wǎng)格的網(wǎng)格數(shù)之和大于原始網(wǎng)格數(shù)。由此易知，即使排除負(fù)載不均衡和進(jìn)程間通信的影響，也難以達(dá)到理想的并行加速比。

在化學(xué)非平衡流的計(jì)算中，化學(xué)反應(yīng)源項(xiàng)求解的計(jì)算量往往顯著大于流動(dòng)求解的計(jì)算量。尤其是當(dāng)化學(xué)反應(yīng)機(jī)理趨于復(fù)雜時(shí)，化學(xué)反應(yīng)源項(xiàng)計(jì)算量在總計(jì)算量中的比例也隨之進(jìn)一步增大。化學(xué)反應(yīng)源項(xiàng)的求解僅依賴于離散網(wǎng)格單元內(nèi)的物理量，特別是在2階精度有限體積法的框架下，源項(xiàng)僅依賴于單元內(nèi)變量的平均值/格心值。因此，一個(gè)時(shí)間步內(nèi)化學(xué)反應(yīng)源項(xiàng)的求解不需要進(jìn)行進(jìn)程/分區(qū)間通信，從而能夠提高并行計(jì)算效率。

3 數(shù)值算例

為了驗(yàn)證化學(xué)非平衡流并行計(jì)算程序的準(zhǔn)確性，本文以Lehr[10]于1972年進(jìn)行的激波誘導(dǎo)氫氣/空氣周期性振蕩燃燒實(shí)驗(yàn)為依據(jù)，對(duì)計(jì)算程序進(jìn)行考核。在彈道靶實(shí)驗(yàn)中觀察到：隨著彈丸飛行速度的增加，激波誘導(dǎo)燃燒逐漸失穩(wěn)，出現(xiàn)周期性振蕩燃燒現(xiàn)象；進(jìn)一步增大彈丸飛行速度，激波誘導(dǎo)燃燒的流場(chǎng)又趨于穩(wěn)定。實(shí)驗(yàn)中得到了不同飛行條件下的紋影圖像和不穩(wěn)定燃燒的振蕩頻率。本文選取Ma=4.48條件下振蕩燃燒現(xiàn)象進(jìn)行測(cè)試。

3.1 計(jì)算模型

圖2 計(jì)算網(wǎng)格示意圖Fig.2 Schematic of computation grid

實(shí)驗(yàn)中彈丸頭部直徑D=15mm，本文計(jì)算中采用軸對(duì)稱模型，邊界條件設(shè)置如圖2所示，計(jì)算網(wǎng)格為500×400，網(wǎng)格最小尺度Δ=4.5×10-6m。超聲速入口邊界面元給定來流參數(shù)；由于該算例中黏性影響微弱，故彈丸表面采用滑移壁面邊界條件；下游出口采用超聲速出流邊界，來流參數(shù)如表1所示。計(jì)算選取修正的Jachimowski[11]9組分19反應(yīng)機(jī)理來描述H2/Air反應(yīng)過程，其中N2作為第三體影響化學(xué)反應(yīng)進(jìn)程。計(jì)算采用4步Rung-Kutta方法進(jìn)行時(shí)間離散，CFL數(shù)取0.6。

表1 激波誘導(dǎo)燃燒來流參數(shù)

3.2 數(shù)值模擬結(jié)果

圖3給出了激波誘導(dǎo)周期性振蕩燃燒的實(shí)驗(yàn)紋影圖。可以看出，彈丸前端存在光滑的脫體激波。由于彈丸飛行馬赫數(shù)略低于來流參數(shù)條件下CJ爆速，混合氣體經(jīng)過激波壓縮后增溫增壓，在下游發(fā)生周期性振蕩燃燒，從圖3中看出，燃燒陣面呈現(xiàn)出規(guī)則的波紋狀，在燃燒陣面和激波間存在著未燃混合氣體。

圖3 周期性振蕩燃燒實(shí)驗(yàn)紋影圖Fig.3 Experiment schlieren of periodic oscillating combustion

采用3.1節(jié)參數(shù)進(jìn)行數(shù)值模擬，得到了駐點(diǎn)流線上(即對(duì)稱邊界)密度隨時(shí)間變化曲線，如圖4所示。可以看出，激波和燃燒陣面呈現(xiàn)出顯著的周期性振蕩，來流沿駐點(diǎn)流線通過激波后增溫增壓，化學(xué)反應(yīng)速率增大，在誘導(dǎo)區(qū)內(nèi)經(jīng)過誘導(dǎo)時(shí)間后釋放能量，進(jìn)而發(fā)生化學(xué)反應(yīng)，形成燃燒陣面。彈丸前激波到燃燒陣面間的距離即為誘導(dǎo)區(qū)寬度，誘導(dǎo)區(qū)范圍由激波后的誘導(dǎo)時(shí)間決定，駐點(diǎn)流線上彈丸頭激波強(qiáng)度接近正激波，誘導(dǎo)時(shí)間最短，因此駐點(diǎn)線上誘導(dǎo)區(qū)寬度最小。振蕩燃燒過程中，隨機(jī)擾動(dòng)產(chǎn)生一道壓縮波，由于正激波后為亞聲速區(qū)域，因此該壓縮波向上游傳播至頭激波，增強(qiáng)了激波強(qiáng)度，同時(shí)形成向下游傳播的弱反射聲波和一道接觸間斷。隨著激波強(qiáng)度的增強(qiáng)，波后溫度升高，誘導(dǎo)距離縮短，燃燒陣面靠近頭激波，此時(shí)化學(xué)反應(yīng)速率增大，能量釋放速率的短暫增大又會(huì)產(chǎn)生新的向上游傳播的壓縮波，流場(chǎng)進(jìn)入新一周期的振蕩燃燒。此外，由于燃燒陣面的位置前移，由頭激波向下游傳播的接觸間斷到達(dá)原始燃燒陣面時(shí)，無可燃?xì)怏w支持化學(xué)反應(yīng)放熱，因此形成了向上游傳播的膨脹波，膨脹波到達(dá)激波陣面時(shí)與激波相互作用降低了激波強(qiáng)度，減緩了波后放熱速率，誘導(dǎo)區(qū)寬度重新增加，燃燒陣面向下游運(yùn)動(dòng)，恢復(fù)至初始位置。綜上所述，振蕩燃燒的根本原因在于壓縮波和膨脹波對(duì)激波及燃燒陣面的周期性影響。

圖4 駐點(diǎn)線上密度隨時(shí)間變化曲線Fig.4 Temporal variation of density along stagnation line

圖5顯示了流場(chǎng)及H2O2組分的密度分布?？梢钥闯觯紵嚸娲嬖谟诩げㄏ掠?，并未與激波耦合，因此無法形成穩(wěn)定爆轟結(jié)構(gòu)。在化學(xué)反應(yīng)過程中，燃燒陣面不再保持光滑，呈現(xiàn)周期性振蕩。圖6記錄了駐點(diǎn)處壓力隨時(shí)間變化曲線，提取30μs～45μs時(shí)間段內(nèi)數(shù)據(jù)進(jìn)行傅里葉變換，得到了周期性振蕩頻率為417kHz，與實(shí)驗(yàn)數(shù)據(jù)[10]425kHz和文獻(xiàn)計(jì)算結(jié)果[12]431kHz符合較好，證明在并行計(jì)算條件下求解器具有良好的計(jì)算精度。

圖5 振蕩燃燒密度分布云圖Fig.5 Density contour of periodic oscillating combustion

圖6 駐點(diǎn)壓力隨時(shí)間變化曲線Fig.6 Temporal variation of the pressure at stagnation point

3.3 計(jì)算效率分析

本文的計(jì)算是在44核、主頻2.2GHz、內(nèi)存128GB配置的惠普Z840工作站上完成，分別進(jìn)行了串行、14核、16核和32核并行計(jì)算。

圖7和表2為4種計(jì)算條件下CPU計(jì)算時(shí)間對(duì)比?？梢钥闯?，隨著核心個(gè)數(shù)的增加，所占用CPU時(shí)間也相應(yīng)增大，主要是因?yàn)楹诵臄?shù)的增加引起不同進(jìn)程間通信時(shí)間增大。進(jìn)一步對(duì)比了流動(dòng)算子和化學(xué)反應(yīng)算子所占CPU時(shí)間，如圖7中虛線所示。可知不同進(jìn)程間通信時(shí)間的增大主要體現(xiàn)在流動(dòng)算子的求解，而化學(xué)反應(yīng)的求解受進(jìn)程數(shù)量的影響較小。這一結(jié)果也支撐了前文中的觀點(diǎn)：由于不需要進(jìn)行進(jìn)程/分區(qū)間通信，化學(xué)反應(yīng)源項(xiàng)的并行求解易于實(shí)現(xiàn)較高的效率。

(a)總計(jì)算時(shí)間

(b)流動(dòng)與反應(yīng)計(jì)算時(shí)間圖7 周期性振蕩燃燒算例CPU計(jì)算時(shí)間對(duì)比Fig.7 Comparison of the CPU time for the periodic oscillating combustion

圖8和表3為3種并行計(jì)算條件下墻上時(shí)間對(duì)比。隨著并行核心個(gè)數(shù)的增加，墻上時(shí)間顯著減小，大大縮短了計(jì)算周期。特別是采用32核并行計(jì)算的墻上時(shí)間相對(duì)于16核并行計(jì)算的墻上時(shí)間降低了46%；而3個(gè)并行計(jì)算算例的加速比接近線性。并行計(jì)算程序與串行程序相比并未實(shí)現(xiàn)線性加速，這是由于并行化程序的總墻上時(shí)間中包括了(串行程序不需要的)網(wǎng)格分區(qū)在內(nèi)的一系列功能，因此這也是程序進(jìn)一步優(yōu)化的重要內(nèi)容。

表2 CPU計(jì)算時(shí)間對(duì)比

圖8 周期性振蕩燃燒算例墻上時(shí)間對(duì)比Fig.8 Comparison of the Wall time for the periodic oscillating combustion

表3 墻上時(shí)間對(duì)比

Tab.3 Comparison of the Wall time for the periodicoscillating combustion

算例墻上時(shí)間 串行14核16核32核t/s516282704536234533635

4 結(jié)論

本文介紹了作者所在研究團(tuán)隊(duì)在化學(xué)非平衡流數(shù)值模擬領(lǐng)域的研究工作。其中所采用的數(shù)值算法是目前較為先進(jìn)的化學(xué)非平衡流解耦算法，能夠有效地改善數(shù)值計(jì)算剛性的問題，并且已經(jīng)得到多種工程問題的驗(yàn)證；基于該化學(xué)非平衡流求解器實(shí)現(xiàn)了MPI并行計(jì)算，并且基于經(jīng)典算例進(jìn)行了測(cè)試。結(jié)果表明，該求解器具有良好的計(jì)算精度和計(jì)算效率。后續(xù)進(jìn)一步將該求解器應(yīng)用于高超聲速飛行器工程中，拓展其GPU并行計(jì)算的能力，實(shí)現(xiàn)超大規(guī)模的化學(xué)非平衡流模擬和研究。