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

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

0 引言

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

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

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

1 控制方程和解耦算法

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

(1)

其中

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

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

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

(2)

進而可將能量方程改寫為

(3)

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

(4)

最終將控制方程改寫為

(5)

其中

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

(6)

(7)

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

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

(8)

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

2 并行計算

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

圖1給出了網(wǎng)格分區(qū)的一個例子。方形區(qū)域內(nèi)的三角形網(wǎng)格被分為A、B兩個分區(qū)，如圖1(a)所示。初始分區(qū)A、B中沒有重合的單元。由于使用空間2階格式進行計算，在計算一個單元的通量時，需要用到直接相鄰單元的物理量及其梯度；而計算直接相鄰單元的物理量梯度時需要第二層相鄰單元的物理量。這里稱兩層相鄰單元為計算分區(qū)的“外部單元”，稱分區(qū)原有單元為“內(nèi)部單元”，分別如圖1(b)、(c)所示。其中，標記為1的是直接相鄰單元，標記為2的是第二層相鄰單元。為了減少計算過程中的通信次數(shù)，需要存儲外部單元的物理量。因此，各個進程存儲的子域網(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中容易看出，對于多個計算分區(qū)，一個分區(qū)的外部單元必定是另一個分區(qū)的內(nèi)部單元，而一個分區(qū)的內(nèi)部單元可能對應(yīng)著多個分區(qū)的外部單元。計算中，在每一步時間推進后，外部單元上的物理量需要從所對應(yīng)的內(nèi)部單元獲得，以保證不同子網(wǎng)格間流場物理量的一致性。因此，在進行網(wǎng)格分區(qū)時，需要建立所有分區(qū)間的內(nèi)部單元和外部單元的映射關(guān)系。

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

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

3 數(shù)值算例

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

3.1 計算模型

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

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

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

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

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

圖3 周期性振蕩燃燒實驗紋影圖Fig.3 Experiment schlieren of periodic oscillating combustion

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

圖4 駐點線上密度隨時間變化曲線Fig.4 Temporal variation of density along stagnation line

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

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

圖6 駐點壓力隨時間變化曲線Fig.6 Temporal variation of the pressure at stagnation point

3.3 計算效率分析

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

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

(a)總計算時間

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

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

表2 CPU計算時間對比

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

表3 墻上時間對比

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

算例墻上時間 串行14核16核32核t/s516282704536234533635

4 結(jié)論

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