吳錦濤， 董剛， 栗保明

(南京理工大學(xué) 瞬態(tài)物理國家重點(diǎn)實(shí)驗(yàn)室， 江蘇 南京 210094)

0 引言

沖壓加速推進(jìn)技術(shù)作為一種新型發(fā)射技術(shù)，可以有效地提高彈丸的發(fā)射初速，為超遠(yuǎn)程火炮和超高速穿甲效應(yīng)提供技術(shù)支持。沖壓推進(jìn)過程可以以亞爆轟速度、跨爆轟速度和超爆轟速度為模式工作，然而，其中涉及的可燃預(yù)混氣爆轟過程十分復(fù)雜[1]，因此開展相關(guān)研究具有重要的理論意義和應(yīng)用前景。

從數(shù)值模擬的角度看，要精確模擬氣相爆轟過程，往往要使用包含多組分和多基元反應(yīng)的詳細(xì)化學(xué)反應(yīng)機(jī)理。然而，在實(shí)際問題中由于組分之間的化學(xué)反應(yīng)時(shí)間尺度相差懸殊[2-3]，反應(yīng)源項(xiàng)的求解通常存在嚴(yán)重的剛性現(xiàn)象，導(dǎo)致計(jì)算量很大。為提高化學(xué)反應(yīng)的計(jì)算效率，傳統(tǒng)的處理方法常常涉及到對(duì)化學(xué)反應(yīng)機(jī)理的簡(jiǎn)化[4-7]，但簡(jiǎn)化機(jī)理一般依賴于所研究的具體物理問題，缺乏通用性，而且對(duì)于某些重要的痕量物質(zhì)往往無法預(yù)測(cè)或是嚴(yán)重失真，因而受到不同程度的限制。

1997年，Pope[8]提出了動(dòng)態(tài)自適應(yīng)建表(ISAT)方法，以減少反應(yīng)流數(shù)值模擬中化學(xué)反應(yīng)源項(xiàng)的計(jì)算時(shí)間。不同于傳統(tǒng)的機(jī)理簡(jiǎn)化方法，該方法采用在內(nèi)存中建立數(shù)據(jù)表的方式對(duì)化學(xué)反應(yīng)的數(shù)據(jù)進(jìn)行存儲(chǔ)、查找和取回，通過這些操作來替代原有反應(yīng)源項(xiàng)的直接積分(DI)運(yùn)算。由于數(shù)據(jù)操作花費(fèi)的時(shí)間遠(yuǎn)小于DI的時(shí)間，ISAT方法可以有效地提高化學(xué)反應(yīng)源項(xiàng)的計(jì)算效率。在低馬赫數(shù)定常燃燒問題中，流場(chǎng)內(nèi)一個(gè)步長(zhǎng)前后的化學(xué)反應(yīng)熱力學(xué)狀態(tài)變化不大，大部分?jǐn)?shù)據(jù)可以在表中查找成功并取回，因此ISAT方法在此類問題中具有很高的加速比。例如：在均質(zhì)壓燃(HCCI)發(fā)動(dòng)機(jī)[9]中，其化學(xué)反應(yīng)加速比相對(duì)于DI能夠達(dá)到100倍以上；而在混合攪拌發(fā)生器(PSRs)[10]中，建表方法可縮短反應(yīng)計(jì)算時(shí)間約1/42.

然而，國內(nèi)外針對(duì)ISAT方法在可壓縮流的非定常燃燒問題中的應(yīng)用卻報(bào)道較少。事實(shí)上在這類問題中，考慮到每個(gè)時(shí)刻流場(chǎng)內(nèi)熱力學(xué)狀態(tài)相對(duì)上一時(shí)刻變化較大，大量數(shù)據(jù)可能無法從數(shù)據(jù)表取回，因此就會(huì)有大量節(jié)點(diǎn)不斷插入表中為后續(xù)計(jì)算提供可查找的數(shù)據(jù)，這嚴(yán)重影響了查表效率并增加了表的尺寸，因此ISAT算法的加速效果必定會(huì)受到不同程度的限制。Dong等[11]考慮流場(chǎng)存在間斷的非定常燃燒特點(diǎn)，將ISAT算法應(yīng)用到二維氣相爆轟波引發(fā)和傳播的數(shù)值模擬中，在0階取回判據(jù)的條件下，其化學(xué)加速比最大可達(dá)到17.88，最大總加速比可達(dá)到4.75. 在Dong等[12]后續(xù)研究中，還進(jìn)一步提出了兩種節(jié)點(diǎn)刪除(全局節(jié)點(diǎn)刪除和局部節(jié)點(diǎn)刪除)機(jī)制來避免表的過度增長(zhǎng)，并針對(duì)爆轟波傳播的瞬態(tài)過程計(jì)算對(duì)比了兩種刪除機(jī)制下ISAT方法的效率。

上述研究表明，盡管ISAT方法本身通用性好，但其計(jì)算效率仍依賴于物理問題。隨著近年來大規(guī)模并行計(jì)算的日益普及，ISAT方法也逐漸應(yīng)用于并行計(jì)算中[13]，但針對(duì)可壓縮非定常問題，目前尚未提出有效的并行算法。而以往滿足可壓縮非定常問題的算法則是基于串行計(jì)算的[12]，因此，發(fā)展基于ISAT的適合于可壓縮非定常燃燒問題的有效并行算法十分必要?？蓧嚎s氣相爆轟波精細(xì)結(jié)構(gòu)的捕捉是非定?？蓧嚎s反應(yīng)流計(jì)算的重要問題之一。氣相爆轟波是化學(xué)反應(yīng)放熱支持下的激波，其化學(xué)反應(yīng)特性決定了波陣面反應(yīng)區(qū)的大小并因而影響了爆轟波的傳播過程。另一方面，由于爆轟波陣面主要由各類強(qiáng)間斷和接觸間斷所構(gòu)成，其發(fā)展是不穩(wěn)定的，要想精確捕捉精細(xì)結(jié)構(gòu)和不穩(wěn)定發(fā)展過程，必須對(duì)計(jì)算網(wǎng)格尺寸有十分苛刻的要求，因而計(jì)算量巨大。因此，針對(duì)這類問題來發(fā)展基于ISAT的并行算法并分析其性能(計(jì)算精度和計(jì)算效率)是非常適合的。

本文以二維氣相爆轟波傳播的過程為研究對(duì)象，在ISAT方法的基礎(chǔ)上，建立了基于多線程信息傳遞接口(MPI)并行架構(gòu)的化學(xué)反應(yīng)數(shù)據(jù)動(dòng)態(tài)存儲(chǔ)、刪除加速算法，考慮了不同的建表方法、建表容差判據(jù)和流場(chǎng)計(jì)算格式等因素的影響，通過數(shù)值模擬分析了并行加速算法的計(jì)算精度和計(jì)算效率。

1 計(jì)算方法

1.1 動(dòng)態(tài)存儲(chǔ)、刪除的ISAT算法

在含有K種組分的反應(yīng)體系中，其熱力學(xué)狀態(tài)函數(shù)所構(gòu)成的狀態(tài)空間可表達(dá)為

Φ=(φ1,…，φk, …，φK,φK+1,φK+2)=
(Y1,…，Yk,…，YK,T,ρ),

(1)

式中：Yk為組分k的質(zhì)量分?jǐn)?shù)；T為體系溫度；ρ為體系密度。對(duì)燃燒過程，流場(chǎng)內(nèi)每個(gè)網(wǎng)格(顆粒)內(nèi)的熱力學(xué)狀態(tài)在規(guī)定步長(zhǎng)Δt內(nèi)發(fā)生的化學(xué)變化可以表達(dá)為

(2)

(2)式給出了狀態(tài)函數(shù)積分的求解過程，稱為DI，其中S代表化學(xué)反應(yīng)速率源項(xiàng)，由于其具有強(qiáng)剛性和強(qiáng)非線性，因此求解將花費(fèi)大量CPU時(shí)間。

在ISAT算法中，DI可被數(shù)據(jù)表的操作所取代。計(jì)算過程中首先在內(nèi)存中建立一個(gè)空白數(shù)據(jù)表，并且將初次計(jì)算得到的Φ及其對(duì)應(yīng)的變化量ΔΦ存入數(shù)據(jù)表的一個(gè)節(jié)點(diǎn)中，對(duì)后續(xù)計(jì)算中的任一變量Φ′，在數(shù)據(jù)表中先查找與之近似的Φ，并將表中滿足容差范圍的節(jié)點(diǎn)值Φ+ΔΦ取回作為Φ′的計(jì)算結(jié)果，以避免DI計(jì)算。否則，進(jìn)行(2)式DI并將得到的結(jié)果存入數(shù)據(jù)表的新節(jié)點(diǎn)中。由于對(duì)表中節(jié)點(diǎn)的存儲(chǔ)和取回操作所花費(fèi)的時(shí)間遠(yuǎn)小于DI所花費(fèi)的時(shí)間，這種方法可明顯節(jié)約計(jì)算的CPU時(shí)間[8]。對(duì)容差范圍的判定，本文采用了常數(shù)近似[11-12]以替代原始ISAT中的線性近似[8]方法，以進(jìn)一步節(jié)約對(duì)數(shù)據(jù)表中節(jié)點(diǎn)的操作時(shí)間。同時(shí)，設(shè)定容差判據(jù)εtol以保證待查數(shù)據(jù)能夠落入誤差域內(nèi)，否則將DI結(jié)果及其誤差域作為一個(gè)數(shù)據(jù)簇存入表中新節(jié)點(diǎn)上。事實(shí)上，容差判據(jù)的大小決定了ISAT對(duì)節(jié)點(diǎn)數(shù)據(jù)的容忍范圍，直接影響數(shù)據(jù)查找和取回的效率。

氣相爆轟波計(jì)算中，不同節(jié)點(diǎn)代表的熱力學(xué)狀態(tài)不同，因而其被取回利用的次數(shù)也不同。為避免數(shù)據(jù)表占用內(nèi)存超限，設(shè)定最大容量M，計(jì)算過程中，一旦表內(nèi)節(jié)點(diǎn)總量超過M，則對(duì)數(shù)據(jù)表進(jìn)行清理。本文采用平衡二叉樹(BBT)結(jié)構(gòu)來構(gòu)造數(shù)據(jù)表，以確保對(duì)表中節(jié)點(diǎn)數(shù)據(jù)查找、取回、存儲(chǔ)以及刪除等操作能夠高效完成。為了保留表中利用率高的節(jié)點(diǎn)，通過非遞歸中序遍歷方法對(duì)表內(nèi)節(jié)點(diǎn)進(jìn)行查找，將表中從未取回過的節(jié)點(diǎn)進(jìn)行刪除[12]，以此來控制數(shù)據(jù)表的大小。

1.2 MPI并行計(jì)算中的建表策略

對(duì)于包含詳細(xì)化學(xué)反應(yīng)機(jī)理的氣相爆轟數(shù)值模擬問題，由于所涉及的流場(chǎng)結(jié)構(gòu)比較復(fù)雜且對(duì)網(wǎng)格分辨率要求很高，普通的串行執(zhí)行方式已經(jīng)無法滿足計(jì)算需求，故本文采用MPI方法對(duì)整個(gè)計(jì)算過程進(jìn)行多線程并行處理。

MPI并行環(huán)境下，ISAT建立熱力學(xué)數(shù)據(jù)表的思路可以分為兩種：第1種是各個(gè)線程分別獨(dú)自建立數(shù)據(jù)表，在表所對(duì)應(yīng)的分區(qū)內(nèi)各線程對(duì)表中數(shù)據(jù)同時(shí)進(jìn)行插入、查找和取回等操作，整個(gè)計(jì)算過程中各數(shù)據(jù)表彼此互不關(guān)聯(lián)；第2種則是在整個(gè)計(jì)算區(qū)域內(nèi)建立一個(gè)數(shù)據(jù)表，所有線程在計(jì)算過程中共享該表中的數(shù)據(jù)并依次對(duì)其進(jìn)行操作。

第1種建表方法的優(yōu)勢(shì)在于各線程可同時(shí)分別對(duì)其分區(qū)內(nèi)的數(shù)據(jù)進(jìn)行操作，避免了數(shù)據(jù)操作的順序等待時(shí)間，因此具有更高的計(jì)算效率；第2種建表方法通過表中數(shù)據(jù)在所有分區(qū)共享的方式避免了第一種方法中不同數(shù)據(jù)表內(nèi)數(shù)據(jù)重復(fù)的不足，但共享數(shù)據(jù)表的劣勢(shì)主要表現(xiàn)在對(duì)數(shù)據(jù)進(jìn)行操作時(shí)各線程需要排隊(duì)等待。

由于不同的建表策略各有優(yōu)缺點(diǎn)，本文將分析兩種建表策略在氣相爆轟波傳播數(shù)值計(jì)算中的性能。

2 二維氣相爆轟的計(jì)算模型

2.1 控制方程和數(shù)值格式

針對(duì)二維氣相爆轟波傳播問題，本文采用理想氣體模型，控制方程為二維帶化學(xué)反應(yīng)的多組分Euler方程：

(3)

計(jì)算時(shí)采用了分裂算法，將控制方程(3)式中的流動(dòng)過程和化學(xué)反應(yīng)過程解耦，針對(duì)每個(gè)步長(zhǎng)，先不考慮化學(xué)反應(yīng)，只進(jìn)行流動(dòng)過程的計(jì)算，再利用經(jīng)過流動(dòng)過程的新參數(shù)計(jì)算化學(xué)反應(yīng)源項(xiàng)。分別采用基于ISAT的并行加速算法和DI方法來求解反應(yīng)源項(xiàng)。流動(dòng)過程采用了高精度加權(quán)本質(zhì)非振蕩(WENO)格式進(jìn)行求解，為分析不同數(shù)值格式對(duì)并行加速算法的計(jì)算精度和計(jì)算效率的影響，分別采用了5階WENO格式和9階WENO格式兩種數(shù)值格式，它們的具體形式可參見文獻(xiàn)[16]. 最后，采用3階顯式Runge-Kutta方法進(jìn)行時(shí)間推進(jìn)，計(jì)算的定步長(zhǎng)為1.0×10-9s.

2.2 計(jì)算模型

圖1為氣相爆轟波傳播的初始二維計(jì)算區(qū)域示意圖，計(jì)算區(qū)域Ly×Lx為9 mm×3 mm的長(zhǎng)方形。采用均勻正交網(wǎng)格覆蓋整個(gè)計(jì)算域，網(wǎng)格數(shù)量為900×300(y×x)，網(wǎng)格尺寸Δx=Δy=0.01 mm. 沿著y方向?qū)⒘鲌?chǎng)分成15個(gè)相同大小的物理分區(qū)，各進(jìn)程(進(jìn)程1～進(jìn)程15，見圖1)在對(duì)應(yīng)的分區(qū)并行計(jì)算。計(jì)算區(qū)域左邊界、右邊界采用周期性邊界條件，上邊界、下邊界采用零梯度邊界條件，為保證爆轟波穩(wěn)定傳播的距離足夠長(zhǎng)，令計(jì)算區(qū)域隨著爆轟波傳播一起向上方移動(dòng)。計(jì)算結(jié)果表明，本文所使用的計(jì)算區(qū)域大小能夠完整地捕捉爆轟波陣面的精細(xì)結(jié)構(gòu)，從而達(dá)到無限長(zhǎng)管道中二維爆轟波穩(wěn)定傳播的數(shù)值模擬效果。

計(jì)算初始時(shí)刻，在計(jì)算域y為0～3 mm范圍內(nèi)設(shè)定ZND形式的爆轟波結(jié)構(gòu)，波陣面上方的區(qū)域(y為3～9 mm)內(nèi)充滿了初始?jí)毫0=1 atm和初始溫度T0=300 K、當(dāng)量比為2∶1的H2、O2預(yù)混氣，以維持二維爆轟波的自持傳播。然而，由于ZND模型是理論上爆轟波傳播的穩(wěn)定解，它本身是一維結(jié)構(gòu)解，作為本文二維計(jì)算的初始條件，需要在數(shù)值誤差擾動(dòng)下經(jīng)過很長(zhǎng)時(shí)間才能發(fā)展為不穩(wěn)定的爆轟波結(jié)構(gòu)。因此，為盡快誘導(dǎo)爆轟波形成不穩(wěn)定的胞格結(jié)構(gòu)，本文在一維ZND波陣面上方的中心位置處設(shè)定一道細(xì)長(zhǎng)條形的擾動(dòng)區(qū)(初始?jí)毫?8p0，初始溫度T0不變)，以加速胞格結(jié)構(gòu)的出現(xiàn)。計(jì)算開始后，爆轟波陣面以穩(wěn)定爆速自下向上傳播，受擾動(dòng)作用之后發(fā)生強(qiáng)烈失穩(wěn)并在波面上出現(xiàn)若干三波點(diǎn)，經(jīng)過一段時(shí)間的衰減，爆轟波由過驅(qū)動(dòng)狀態(tài)逐漸轉(zhuǎn)化為穩(wěn)定傳播狀態(tài)。

2.3 計(jì)算驗(yàn)證

考慮到計(jì)算初期由于人為設(shè)定擾動(dòng)區(qū)所導(dǎo)致的過驅(qū)動(dòng)狀態(tài)下爆轟波傳播過程不具有物理上的真實(shí)性，本文從穩(wěn)定傳播后的某個(gè)時(shí)刻(步長(zhǎng)10 000)開始對(duì)計(jì)算結(jié)果進(jìn)行統(tǒng)計(jì)和處理。結(jié)果表明，達(dá)到穩(wěn)定傳播狀態(tài)后的爆轟波傳播速度vD=2 812.0 m/s，與Chapmam-Jouget(C-J)理論爆轟波傳播速度vC-J=2 843.2 m/s的結(jié)果非常接近，證明了本文計(jì)算爆轟波傳播狀態(tài)的準(zhǔn)確性。

圖2給出了穩(wěn)定傳播條件下典型步長(zhǎng)(步長(zhǎng)20 000)時(shí)的爆轟波陣面結(jié)構(gòu)和胞格寬度。從圖2(a)中可以明顯觀察到爆轟波是由入射激波(SW)、馬赫桿(MS)、橫波(TW)及其它們所匯聚形成的三波點(diǎn)(TP)所組成。經(jīng)統(tǒng)計(jì)，爆轟波在穩(wěn)定傳播過程中，陣面上約存在14個(gè)TP，相應(yīng)的胞格寬度尺寸約為0.428 mm. 圖2(b)給出了本文計(jì)算得到的爆轟波胞格寬度與文獻(xiàn)[17]中實(shí)驗(yàn)測(cè)得的等當(dāng)量比H2、O2(與本文相同)爆轟波胞格寬度隨著初始?jí)毫ψ兓谋容^，從中可以發(fā)現(xiàn)胞格尺寸隨著初始?jí)毫Φ脑黾佣鴾p小，本文在p0=1 atm下的計(jì)算結(jié)果與文獻(xiàn)[17]的變化規(guī)律是一致的，證明了本文計(jì)算的爆轟波結(jié)構(gòu)是合理的。

3 結(jié)果與討論

本文以等當(dāng)量比H2、O2氣相爆轟波傳播過程為研究對(duì)象，針對(duì)其中的化學(xué)反應(yīng)過程，分析了建表方法、容差判據(jù)和數(shù)值格式精度等因素對(duì)基于并行計(jì)算的動(dòng)態(tài)存儲(chǔ)、刪除ISAT算法性能影響，性能評(píng)價(jià)指標(biāo)包括計(jì)算精度和計(jì)算效率。基于該算法的不同算例參數(shù)見表1.

表1 算例及對(duì)應(yīng)參數(shù)

表1中，建表方法I代表每個(gè)線程各自獨(dú)立建表的方式，方法II代表所有線程共享一個(gè)數(shù)據(jù)表的建表方式，參見1.2節(jié)。其中，方法I中設(shè)定獨(dú)立建表的數(shù)據(jù)最大容量M=3×106，而方法II中設(shè)定表中數(shù)據(jù)最大容量M=4.5×107，從而保證方法I中所有數(shù)據(jù)表(15個(gè))的總?cè)萘亢头椒↖I中的單個(gè)數(shù)據(jù)表容量相等。

3.1 建表方法

圖3給出了步長(zhǎng)20 000時(shí)刻爆轟波陣面附近的流場(chǎng)壓力及溫度分布圖，以比較兩種不同建表方法(算例1和算例2)的計(jì)算精度，為了便于比較，圖3中還給出了DI對(duì)應(yīng)的計(jì)算結(jié)果。從圖3可以看出，兩種建表方法的計(jì)算結(jié)果和DI相比，陣面上的壓力和溫度分布以及TP個(gè)數(shù)均保持一致，證明建表方法在本文氣相爆轟波計(jì)算中的準(zhǔn)確度，而波系結(jié)構(gòu)顯示出的輕微不規(guī)則性是由9階WENO格式較低的耗散特性引起的，與加速算法本身沒有關(guān)系。

為了定量比較兩種建表方法的計(jì)算精度，圖4給出了步長(zhǎng)20 000時(shí)刻加速算法及DI計(jì)算所得爆轟波陣面附近的組分OH和H2O2(痕量物質(zhì))分布曲線。圖4結(jié)果顯示，兩種建表方法(算例1、算例2)對(duì)應(yīng)的計(jì)算結(jié)果均展現(xiàn)出與DI相似的分布規(guī)律，峰值大小及所處的位置十分接近。從痕量物質(zhì)H2O2的分布曲線可以看出，建表加速算法對(duì)小濃度自由基(10-5量級(jí))的捕捉仍具有較高的準(zhǔn)確性，這是以簡(jiǎn)化機(jī)理為基礎(chǔ)的加速方法所不具備的。

MPI環(huán)境下ISAT加速算法的核心在于通過建立和查找數(shù)據(jù)表的途徑取代熱力學(xué)狀態(tài)的DI，從而能有效地縮短計(jì)算化學(xué)反應(yīng)過程所需的CPU時(shí)間，因此，計(jì)算效率是檢驗(yàn)算法性能的一個(gè)重要指標(biāo)。圖5給出了不同建表方法(算例1和算例2)和DI計(jì)算化學(xué)反應(yīng)所累計(jì)消耗的CPU時(shí)間tc曲線及其對(duì)應(yīng)的加速比Rt(定義為DI方法與ISAT方法計(jì)算化學(xué)反應(yīng)的CPU時(shí)間比)。從圖5(a)可以看出：各線程獨(dú)自建表的ISAT方法(算例1)耗時(shí)較DI明顯縮短，表明該方法在氣相爆轟數(shù)值模擬中的優(yōu)越性；而各線程共享數(shù)據(jù)表的ISAT方法(算例2)在計(jì)算前期(步長(zhǎng)16 230之前)整體效率低于第1種建表思路，但總體上節(jié)約源項(xiàng)計(jì)算的CPU時(shí)間依然比較可觀。之后，該曲線的斜率突然上升并不斷接近DI，甚至在步長(zhǎng)17 560之后超過DI的時(shí)間，這是因?yàn)?6 230時(shí)刻表中節(jié)點(diǎn)總數(shù)達(dá)到最大容量M后ISAT方法對(duì)數(shù)據(jù)表進(jìn)行了節(jié)點(diǎn)刪除操作，接著需要通過DI對(duì)表中節(jié)點(diǎn)進(jìn)行補(bǔ)充。該過程在共享數(shù)據(jù)表的方法中屬于順序操作，相比于算例1同步操作所耗費(fèi)的時(shí)間大幅上升，之后每經(jīng)歷一次表的節(jié)點(diǎn)刪除，時(shí)間曲線都會(huì)出現(xiàn)一個(gè)對(duì)應(yīng)的拐點(diǎn)(如步長(zhǎng)17 490和18 810時(shí)刻)，因此該方法在提高計(jì)算效率上存在一定的局限性。

由圖5(b)可知：加速比曲線定量表明了算例1的優(yōu)勢(shì)，經(jīng)過一小段時(shí)間的波動(dòng)之后最終穩(wěn)定下來的源項(xiàng)計(jì)算時(shí)間大約只占DI的40%；算例2由于數(shù)據(jù)表結(jié)構(gòu)相對(duì)于算例1中小容量數(shù)據(jù)表更復(fù)雜，因此其加速比一開始就低于算例1，同時(shí)，查找、取回和存儲(chǔ)的順序操作導(dǎo)致其加速比一直呈現(xiàn)不斷下降的趨勢(shì)并在計(jì)算后期低于算例1. 考慮到兩種并行ISAT建表方法的計(jì)算精度相當(dāng)，顯然，各線程獨(dú)自建表的ISAT方法(算例1)無疑是更優(yōu)的選擇。

3.2 容差判據(jù)

容差判據(jù)在整個(gè)并行ISAT計(jì)算方法中，直接決定了數(shù)據(jù)表中節(jié)點(diǎn)的取回精度，分析不同容差對(duì)ISAT方法的影響可以兼顧其在氣相爆轟波傳播問題中的計(jì)算精度和計(jì)算效率。從圖3不同容差加速方法(算例1和算例3)的壓力和溫度分布可以看出，相比于DI，在容差較大(算例3)的情況下，加速算法對(duì)數(shù)據(jù)的近似取回更多，這種近似作為數(shù)值擾動(dòng)不斷累加，導(dǎo)致其不同數(shù)據(jù)表之間的差異性變得更明顯，因而導(dǎo)致了波陣面下方出現(xiàn)大小與分區(qū)尺寸相同的低壓區(qū)(見圖3(a)中算例3)，這說明其相對(duì)于小容差加速算法(算例1)會(huì)損失一部分計(jì)算精度。

圖4給出了兩種不同容差對(duì)應(yīng)的加速算法中組分OH和H2O2(痕量物質(zhì))的分布曲線。由圖4可以看出，相對(duì)于小容差(算例1)算法，大容差(算例3)算法與DI計(jì)算結(jié)果相比，雖然仍具有相似的規(guī)律，但流場(chǎng)內(nèi)組分含量值特別是峰值大小存在一定的差異，尤其在痕量物質(zhì)H2O2的分布曲線上，大容差的結(jié)果呈現(xiàn)出比較明顯的震蕩，表明其對(duì)波系結(jié)構(gòu)的刻畫相對(duì)于小容差計(jì)算而言不夠準(zhǔn)確。

為了分析容差放大判據(jù)對(duì)計(jì)算效率的影響，圖6給出了算例1和算例3的節(jié)點(diǎn)取回率Rn，這里Rn定義為數(shù)據(jù)表中數(shù)據(jù)取回次數(shù)和數(shù)據(jù)查找總次數(shù)的比值，該值反映了并行ISAT方法的數(shù)據(jù)查找效率。從圖6可以看出：兩種誤差放大判據(jù)對(duì)應(yīng)的加速算法中，線程1～線程14的數(shù)據(jù)取回率都高于85%，部分線程的取回率甚至在95%以上，證明并行ISAT方法在源項(xiàng)計(jì)算中的數(shù)據(jù)查找效率令人滿意；線程15因其所在分區(qū)包含爆轟波陣面在內(nèi)，TP周期性碰撞使得該分區(qū)熱力學(xué)狀態(tài)變化非常頻繁，進(jìn)而導(dǎo)致二者取回率均有所下降，其中，算例1取回率下降到約60%，算例3取回率下降到約74%，總體上大容差(算例3)算法相對(duì)于小容差(算例1)算法的查找效率更高。

為了進(jìn)一步分析不同容差判據(jù)對(duì)線程15查找效率的影響，圖7給出了該線程一次完整建表過程中節(jié)點(diǎn)取回次數(shù)的分布曲線?？紤]到不同容差下并行ISAT建立熱力學(xué)數(shù)據(jù)表的規(guī)律不同，為了保證一致性，二者均選取時(shí)間步步長(zhǎng)在14 100～14 940之間的典型建表過程。結(jié)果顯示，在包含M=3×106個(gè)節(jié)點(diǎn)的數(shù)據(jù)表中，算例1中約有700個(gè)節(jié)點(diǎn)的取回次數(shù)不小于1，而算例3中約有7 000個(gè)節(jié)點(diǎn)的取回次數(shù)不小于1. 剩余節(jié)點(diǎn)沒有在后續(xù)計(jì)算過程中被取回，因而視為“無效”節(jié)點(diǎn)被刪除，這部分調(diào)用次數(shù)高的節(jié)點(diǎn)在經(jīng)過平衡操作之后仍保留在表中，并在下一次建表過程繼續(xù)使用，從而保證節(jié)點(diǎn)高利用率的同時(shí)有效地節(jié)約內(nèi)存空間。對(duì)比二者可以發(fā)現(xiàn)，節(jié)點(diǎn)取回次數(shù)的分布規(guī)律相同，但算例3由于放大了容差范圍使其能夠被取回重復(fù)利用的節(jié)點(diǎn)數(shù)量和節(jié)點(diǎn)調(diào)用總次數(shù)均有所增加，再次證明了對(duì)容差區(qū)間的放大可以有效地提高并行ISAT方法的查找效率。

根據(jù)圖6和圖7對(duì)數(shù)據(jù)查找效率的對(duì)比，直觀上看，算例3的計(jì)算效率應(yīng)當(dāng)高于算例1，然而，通過圖5 CPU時(shí)間和加速比曲線的對(duì)比可以看出，除計(jì)算初期(步長(zhǎng)12 500之前)算例3具有較高的加速比外，后期不同容差的加速算法具有基本相同的加速比。顯然，在建立、刪除數(shù)據(jù)表的過程中，部分操作“拖累”了算例3的計(jì)算效率。

為探討是何種操作“拖累”了算例3的計(jì)算效率，圖8給出了不同容差條件下幾個(gè)典型線程對(duì)應(yīng)的數(shù)據(jù)表建立、刪除過程。從圖8(a)結(jié)果可以看出，小容差(算例1)計(jì)算中，線程15表的建立、刪除過程十分頻繁，這是由于爆轟波陣面位于此線程，其化學(xué)反應(yīng)和熱力學(xué)狀態(tài)變化劇烈，表的取回率較低(見圖6(a))，因此插入表中的節(jié)點(diǎn)迅速增加，從而導(dǎo)致了頻繁的表建立、刪除過程。相比之下，其他線程所處理的流場(chǎng)熱力學(xué)狀態(tài)變化較小，表中節(jié)點(diǎn)取回率高，表的刪除頻次較小(線程11)甚至不刪除(線程12～線程14)。圖8(b)結(jié)果則表明，容差判據(jù)放大后，線程15的表中節(jié)點(diǎn)取回率提高，表的建立、刪除頻率下降，但由于大容差導(dǎo)致的差異性變大，從上游線程(線程15)計(jì)算得到的結(jié)果傳播到下游線程(線程11～線程14)后，會(huì)導(dǎo)致下游線程取回率下降(見圖6(b))，因此這些線程的建立、刪除頻率有所提高，進(jìn)而使得這些線程中表的建立、刪除同步性顯著下降。因此，從計(jì)算效率來看，容差的區(qū)別并不能補(bǔ)償同步性變差所帶來的影響，故算例1和算例3的加速比之間并沒有顯著差別。

圖8計(jì)算結(jié)果也進(jìn)一步說明，各線程數(shù)據(jù)操作的同步性也是影響并行算法計(jì)算效率的重要原因之一，發(fā)展具有同步性的高效算法也是今后需要進(jìn)一步開展的工作。

3.3 數(shù)值格式精度

計(jì)算的數(shù)值格式可以影響流場(chǎng)分布，進(jìn)而會(huì)改變反應(yīng)體系的熱力學(xué)狀態(tài)分布，這也會(huì)影響加速算法的計(jì)算性能。本節(jié)采用了5階WENO計(jì)算格式替代前面的9階WENO計(jì)算格式(算例4)，以此分析數(shù)值格式的精度對(duì)加速算法影響。

圖9給出了步長(zhǎng)20 000時(shí)刻爆轟波陣面附近流場(chǎng)的壓力和溫度分布，其中，為了便于比較，DI的計(jì)算結(jié)果也是采用5階WENO格式計(jì)算的。結(jié)果顯示，采用5階WENO格式，所有計(jì)算結(jié)果中三波結(jié)構(gòu)分布得更整齊，爆轟波穩(wěn)定傳播之后陣面上的TP結(jié)構(gòu)為12個(gè)，相比于9階WENO格式計(jì)算的TP個(gè)數(shù)(圖3中為14個(gè))有所減少，這說明數(shù)值格式對(duì)計(jì)算結(jié)果有明顯影響。從圖9的5階WENO格式下算例4與DI計(jì)算結(jié)果的比較來看，并行ISAT加速算法的結(jié)果仍能與DI的結(jié)果保持一致，尤其是波陣面及波系結(jié)構(gòu)的形態(tài)與DI的結(jié)果吻合得均很好。

圖10給出了5階WENO格式下ISAT和DI計(jì)算的組分OH和H2O2分布曲線。從圖10可以看出，由于爆轟波陣面分布和波系結(jié)構(gòu)更加規(guī)則，并行ISAT算法得到的結(jié)果與DI吻合很好，即便對(duì)濃度很小的痕量物質(zhì)H2O2，二者在分布規(guī)律和組分含量的峰值大小、位置上均吻合得很好(見圖10(b))，再次表明了并行ISAT加速算法在計(jì)算精度上的優(yōu)勢(shì)。

圖11為5階WENO格式下算例4的加速比曲線。圖11結(jié)果顯示，隨著爆轟波的傳播，加速比曲線經(jīng)過一段時(shí)間的波動(dòng)后，最終穩(wěn)定在約2.2倍左右，這與9階WENO格式下算例1穩(wěn)定后的加速比非常接近，說明并行ISAT加速算法的計(jì)算效率對(duì)數(shù)值格式的依賴性不強(qiáng)。

綜合以上討論可知，采用9階WENO格式的算例1和5階WENO格式的算例4在計(jì)算精度和計(jì)算效率上均表現(xiàn)出較好的性能。表2總結(jié)了計(jì)算截止時(shí)刻(步長(zhǎng)20 000)兩個(gè)算例各項(xiàng)計(jì)算時(shí)間以及加速比的情況，并與相應(yīng)計(jì)算格式下的DI計(jì)算結(jié)果進(jìn)行了比較，其中，總時(shí)間ttot定義為整個(gè)計(jì)算花費(fèi)的總CPU時(shí)間，總加速比Rtot定義為DI計(jì)算總CPU時(shí)間與加速算法計(jì)算總CPU時(shí)間的比值，φc表示化學(xué)反應(yīng)計(jì)算花費(fèi)的CPU時(shí)間在總CPU時(shí)間中的占比。

表2 不同數(shù)值格式下的計(jì)算效率統(tǒng)計(jì)

從表2中的數(shù)據(jù)可以看出，兩種數(shù)值格式下，加速算法均能有效地提高整體計(jì)算效率，其總加速比都維持在1.5倍以上。相對(duì)于5階WENO格式而言，9階WENO下的并行加速算法可以將整個(gè)計(jì)算過程中化學(xué)反應(yīng)計(jì)算時(shí)間的占比從DI方法中的58.5%降低至36.7%，表明本文提出的并行算法在高精度格式下的優(yōu)勢(shì)更明顯。

4 結(jié)論

本文在并行計(jì)算的基礎(chǔ)上，針對(duì)非定?？蓧嚎s反應(yīng)流中化學(xué)反應(yīng)源項(xiàng)計(jì)算量巨大的特點(diǎn)，提出了各線程獨(dú)自建表和所有線程共有單表的兩種基于MPI環(huán)境的多線程并行策略，并成功地將其與帶節(jié)點(diǎn)刪除機(jī)制的ISAT方法相結(jié)合，取代原始的DI運(yùn)算，以在不損失計(jì)算精度的同時(shí)提高計(jì)算效率。為分析不同因素對(duì)計(jì)算結(jié)果的影響，本文從建表方式、容差范圍和數(shù)值計(jì)算格式3個(gè)方面給出了并行加速算法在反應(yīng)性瞬態(tài)流場(chǎng)中的計(jì)算精度和計(jì)算效率。所得結(jié)論如下：

1) 相比于DI方法，各線程獨(dú)自建表的并行加速方法在保證計(jì)算精度的同時(shí)，化學(xué)反應(yīng)的計(jì)算加速比可達(dá)到2.17～2.43.

2) 本文提出的兩種加速方法計(jì)算精度與DI方法相比幾乎一致，但各線程獨(dú)自建表的算法具有更高的計(jì)算效率。

3) 不同容差范圍下并行ISAT算法的數(shù)據(jù)查找效率不同，但小容差下的加速算法在保證更高精度的同時(shí)，仍具有和大容差條件下相當(dāng)?shù)挠?jì)算效率。

4) 不同數(shù)值格式下，并行加速算法在保證計(jì)算精度的基礎(chǔ)上，具有大致相當(dāng)?shù)募铀俦?，表明本文算法的?jì)算效率對(duì)數(shù)值格式并不敏感。

基于并行ISAT方法的加速算法不依賴于具體的物理問題，有很好的可移植性，為數(shù)值研究沖壓加速推進(jìn)技術(shù)中的氣相爆轟過程提供了一種高效的計(jì)算手段。需要指出的是，針對(duì)其他不同類型的非定?？蓧嚎s反應(yīng)流問題，本文所提出的并行加速算法性能仍有待于進(jìn)一步考察。同時(shí)，后續(xù)工作將拓展現(xiàn)有的建表策略，以進(jìn)一步提高不同線程表操作之間的平衡性和同步性。此外，還要進(jìn)一步研究并行加速算法在非均勻網(wǎng)格或非均勻物理分區(qū)下的計(jì)算精度和計(jì)算效率。

參考文獻(xiàn)(References)

[1] 金志明. 高速推進(jìn)內(nèi)彈道學(xué)[M]. 北京: 國防工業(yè)出版社, 2001.

JIN Zhi-ming. High speed internal ballistics [M]. Beijing: National Defense Industry Press, 2001. (in Chinese)

[2] Blasi J M, Kee R J. In situ adaptive tabulation (ISAT) to accelerate transient computational fluid dynamics with complex heterogeneous chemical kinetics[J]. Computers & Chemical Engineering, 2016, 84(1): 36-42.

[3] 胡宗民, 牟乾輝, 張德良,等. 爆轟波在彎管內(nèi)傳播過程數(shù)值分析[J]. 計(jì)算物理, 2004, 21(5):408-414.

HU Zong-min, MU Qian-hui, ZHANG De-liang, et al. Numerical simulation of gaseous detonation wave propagation through bends with a detailed chemical reaction model [J]. Chinese Journal of Computational Physics, 2004, 21(5): 408-414. (in Chinese)

[4] Oran E S, Gamezo V N. Origin of the deflagration-to-detonation transition in gas-phase combustion[J]. Combustion and Flame, 2007, 148(1): 4-47.

[5] Law C K. Comprehensive description of chemistry in combustion modeling [J]. Combustion Science and Technology, 2005, 177(5): 845-870.

[6] Massias A, Diamantis D, Mastorakos E, et al. An algorithm for the construction of global reduced mechanisms with CSP data[J]. Combustion and Flame, 1999, 117(4): 685-708.

[7] Nafe J, Maas U. Hierarchical generation of ILDMs of higher hydrocarbons[J]. Combustion and Flame, 2006, 135(12): 17-26.

[8] Pope S B. Computationally efficient implementation of combustion chemistry using in situ adaptive tabulation[J]. Combustion Theory and Modelling, 1997, 1(1): 41-63.

[9] Zhou L, Wei H. Chemistry acceleration with tabulated dynamic adaptive chemistry in a realistic engine with a primary reference fuel [J]. Fuel, 2016, 171(1): 186-194.

[10] Adhikari S, Sayre A, Chandy A J. In situ adaptive tabulation (ISAT) for combustion chemistry in a network of perfectly stirred reactors (PSRs) [J]. Computers & Chemical Engineering, 2017, 97(1): 124-134.

[11] Dong G, Fan B C, Chen Y L. Acceleration of chemistry computations in two-dimensional detonation induced by shock focusing using reduced ISAT [J]. Combustion Theory and Modeling, 2007, 11(5): 823-837.

[12] Dong G, Fan B C. Chemistry acceleration modeling of detonation based on the dynamical storage/deletion algorithm [J]. Combustion Science and Technology, 2009, 181(9): 1207-1216.

[13] Lu L, Lantz S R, Ren Z, et al, Computationally efficient implementation of combustion chemistry in parallel PDF calculations[J]. Journal of Computational Physics, 2009, 228(15): 5490-5525.

[14] Mcbride B J, Gordon S. Computer program for calculating and fitting thermodynamic functions [M]. Washington, US: NASA， 1992.

[15] Kee R J, Grcar J F, Smooke M D, et al. PREMIX:a Fortran program for modeling steady laminar one-dimensional premixed flames, SAND85-8240 [M].Albuquerque, NM, US: Sandia National Laboratories, 1985.

[16] Shi J, Zhang Y T, Shu C W. Resolution of high order WENO schemes for complicated flow structures[J]. Journal of Computational Physics, 2003, 186(2): 690-696.

[17] Voitsekhovskii B V, Mitrofanov V V, Topchiyan M E. Structure of the detonation front in gases (survey) [J]. Combustion, Explosion and Shock Waves, 1969, 5(3): 267-273.