龍格庫(kù)塔間斷有限元方法在計(jì)算爆轟問(wèn)題中的應(yīng)用

?龍格庫(kù)塔間斷有限元方法在計(jì)算爆轟問(wèn)題中的應(yīng)用 張磊1,2[收稿日期] xxxx-xx-xxx；[修改日期] xxxx-xx-xx [基金項(xiàng)目] 國(guó)家973計(jì)劃（2005CB321703）和國(guó)家自然科學(xué)基金（10531080,10729101）資助項(xiàng)目 ， 袁禮1 1. 中國(guó)科學(xué)院數(shù)學(xué)與系統(tǒng)科學(xué)研究院計(jì)算數(shù)學(xué)所，科學(xué)與工程計(jì)算國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京100190； 2. 中國(guó)礦業(yè)大學(xué)(北京)理學(xué)院，北京100083 [摘要] 構(gòu)造了求解帶源項(xiàng)守恒律方程組的龍格庫(kù)塔間斷有限元（RKDG）方法，并分別結(jié)合源項(xiàng)的Strang分裂法和無(wú)分裂法數(shù)值求解模型守恒律方程和反應(yīng)歐拉方程。為了和有限體積型WENO方法進(jìn)行比較，設(shè)計(jì)了計(jì)算源項(xiàng)的WENO重構(gòu)格式。對(duì)一維帶源項(xiàng)守恒律的計(jì)算表明，對(duì)于非剛性問(wèn)題，RKDG方法比有限體積型FVWENO方法的誤差更小，而對(duì)于剛性問(wèn)題，RKDG方法對(duì)于間斷面位置的捕捉更為精確。對(duì)于一二維爆轟波問(wèn)題的計(jì)算結(jié)果表明，RKDG方法對(duì)于爆轟波結(jié)構(gòu)的分辨和爆轟波位置的捕捉能力更強(qiáng)。 [關(guān)鍵詞] 龍格庫(kù)塔間斷有限元方法；爆轟波；反應(yīng)Euler方程；剛性源項(xiàng) [中圖分類號(hào)] O241.82 [文獻(xiàn)標(biāo)識(shí)碼] A 引言 在非平衡氣體動(dòng)力學(xué)中，爆轟問(wèn)題通常被描述成無(wú)粘性的化學(xué)反應(yīng)流動(dòng)問(wèn)題，其物理模型是一個(gè)非齊次雙曲守恒律方程組，通常被稱作反應(yīng)歐拉方程組（reactive Euler equations），其中非齊次項(xiàng)（源項(xiàng)）通常被解釋為由于化學(xué)反應(yīng)引起的混合組分的質(zhì)量變化率[1]。 最簡(jiǎn)單的反應(yīng)歐拉方程組假定混合氣體僅僅由兩種組分構(gòu)成：已燃?xì)怏w（burnt gas）和未燃?xì)怏w（unburnt gas）[2]。當(dāng)混合氣體達(dá)到點(diǎn)火溫度時(shí)，未燃?xì)怏w通過(guò)一個(gè)不可逆的化學(xué)反應(yīng)轉(zhuǎn)化為已燃?xì)怏w，因此，混合氣體狀態(tài)可以用一個(gè)標(biāo)量，即未燃?xì)怏w的質(zhì)量分?jǐn)?shù)Y來(lái)表示，進(jìn)一步假設(shè)混合氣體各組分具有相同的比熱比和比氣體常數(shù)R，則二維反應(yīng)歐拉方程組可以寫成如下形式 , （1） 其中，各個(gè)變量和函數(shù)的具體形式為 式中是混合氣體的密度，分別是沿方向的速度分量，分別是壓力，總能量，溫度。稱為化學(xué)反應(yīng)率，通常用Arrhenius模型 （2） 其中是化學(xué)反應(yīng)率因子，是活化能。此外混合氣體的狀態(tài)方程和溫度分別為 （3） 其中為單位質(zhì)量未燃?xì)怏w發(fā)生化學(xué)反應(yīng)所釋放的熱量。 在實(shí)際爆轟波問(wèn)題中，化學(xué)反應(yīng)的時(shí)間尺度遠(yuǎn)小于流體流動(dòng)的時(shí)間尺度，所以非齊次方程組具有很強(qiáng)的剛性，給數(shù)值求解帶來(lái)很大困難。無(wú)論對(duì)源項(xiàng)用不用算子分裂方法，采用顯格式還是隱格式，在處理間斷解問(wèn)題時(shí)，都可能得到非物理解。這是因?yàn)楫?dāng)源項(xiàng)不能用足夠的空間和時(shí)間分辨率進(jìn)行求解時(shí)，激波捕捉方法會(huì)得到錯(cuò)誤的爆轟波傳播速度[3,4,5]。文[5]發(fā)現(xiàn)，當(dāng)空間分辨率不夠時(shí)，一個(gè)同時(shí)含有已燃?xì)怏w和未燃?xì)怏w的網(wǎng)格中會(huì)發(fā)生虛假反應(yīng)，數(shù)值爆轟波會(huì)以每個(gè)時(shí)間步一個(gè)空間步長(zhǎng)的速度傳播。文[6]也指出計(jì)算的數(shù)值誤差極可能促成溫度敏感化學(xué)反應(yīng)的提早發(fā)生。因此，要獲得正確的爆轟波位置，常用的解決方案是采用合適的點(diǎn)火溫度模型來(lái)抑制虛假反應(yīng)[3]，或者采用網(wǎng)格自適應(yīng)方法[7, 8, 9]來(lái)保證反應(yīng)區(qū)域內(nèi)有足夠多的網(wǎng)格點(diǎn)，或者采用高精度高分辨率方法來(lái)保證對(duì)反應(yīng)區(qū)有足夠高的分辨率。由于后兩種途徑有更好的普適性而受到更多的關(guān)注。 由Cockburn和Shu等人發(fā)展的龍格庫(kù)塔間斷有限元方法（Runge-Kutta Discontinuous Galerkin Method，RKDG）是一類具有高精度和高分辨率的數(shù)值方法[10]，在解決含有間斷現(xiàn)象的問(wèn)題中發(fā)揮著越來(lái)越重要的作用，它被廣泛地發(fā)展和應(yīng)用于水動(dòng)力學(xué)，氣動(dòng)力學(xué)，波傳播，半導(dǎo)體中的電荷傳輸?shù)葐?wèn)題。該方法既保持了一般有限元方法和有限體積方法的優(yōu)點(diǎn)，又克服了各自的不足。該方法可采用局部高階插值的方法構(gòu)造基函數(shù)，具有靈活處理間斷和邊界條件以及可顯式求解的能力，克服了一般有限元方法不適于間斷問(wèn)題的缺點(diǎn)，以及一般有限體積方法必須通過(guò)擴(kuò)大模板進(jìn)行重構(gòu)來(lái)提高精度的不足。RKDG方法所具有的單元上連續(xù)分布的高精度高分辨率逼近特性有望能更好地模擬爆轟波問(wèn)題。 本文應(yīng)用RKDG方法，結(jié)合化學(xué)源項(xiàng)的無(wú)分裂方法及算子分裂方法，對(duì)爆轟波問(wèn)題進(jìn)行數(shù)值模擬，并和有限體積型WENO（FVWENO）方法的計(jì)算結(jié)果進(jìn)行比較。由于FVWENO方法需要利用流動(dòng)變量單元平均值來(lái)計(jì)算源項(xiàng)的單元積分，本文采用Simpson積分公式，其積分點(diǎn)處的流動(dòng)變量采用WENO重構(gòu)以提高計(jì)算精度。對(duì)于帶源項(xiàng)的一維守恒律的數(shù)值測(cè)試結(jié)果表明，對(duì)于非剛性問(wèn)題，RKDG方法計(jì)算結(jié)果的誤差更??；而對(duì)于剛性問(wèn)題，RKDG方法對(duì)于間斷波位置的捕捉能力更強(qiáng)。對(duì)典型的一二維爆轟波算例，包括二維不穩(wěn)定爆轟波和爆轟波繞過(guò)90°拐角的衍射問(wèn)題的數(shù)值模擬結(jié)果顯示RKDG法仍有一定優(yōu)勢(shì)。 1 數(shù)值方法 1.1 空間離散 本文的空間離散采用RKDG方法[10]。設(shè)是區(qū)域的一個(gè)有限剖分。單元表示多邊形單元的一條邊界，表示單元邊界的外法向。是單元上的局部有限元空間，取作次多項(xiàng)式集合。，在間斷有限元空間 中尋找近似解，其中。 首先在單元上用試探函數(shù)乘以方程（1）的兩端，并用近似解代替方程（1）的精確解，用代替試探函數(shù)，將（1）寫成變分形式，并由分部積分和Green公式，可得 （4） 對(duì)于（4）式的后三項(xiàng)，采用數(shù)值積分計(jì)算 （5） （6） （7） 流通量采用與之相容的數(shù)值流通量代替，數(shù)值流通量定義為為某種形式的Riemann解算器。分別代表邊積分和單元積分的Gauss積分點(diǎn)個(gè)數(shù)，和分別代表邊積分和單元積分的Gauss積分點(diǎn)位置[10]。這樣得到半離散的數(shù)值格式 （8） 為了計(jì)算方便，在單元中取正交基函數(shù)（如勒讓德多項(xiàng)式），則質(zhì)量矩陣成為對(duì)角矩陣，故有限元解表示為 （9） 在式（8）中，取試探函數(shù)為所有基函數(shù)，得到半離散的常微分方程組 （10） 其中，為單元的質(zhì)量矩陣，和分別為 （11） 和 （12） 高階RKDG方法會(huì)產(chǎn)生數(shù)值振蕩，需要在每一個(gè)時(shí)間步后對(duì)數(shù)值解應(yīng)用斜率限制器，本文采用文獻(xiàn)[10]的TVB限制器。當(dāng)基函數(shù)的階數(shù)增加時(shí)，流通量的選擇對(duì)于數(shù)值結(jié)果的影響較小[10]，所以本文采用局部的Lax-Friedrichs流通量。 1.2 無(wú)分裂方法和分裂方法 對(duì)于空間離散后得到的常微分方程組（10），如果源項(xiàng)的剛性不太強(qiáng)，可以直接采用常微分方程組的求解方法進(jìn)行時(shí)間推進(jìn)。本文的無(wú)分裂方法是采用三階強(qiáng)穩(wěn)定保持的龍格庫(kù)塔（SSPRK）方法[11]直接求解方程組（10），該方法在空間離散取階基函數(shù)時(shí)，可以保證光滑解的三階時(shí)間和空間精度。但當(dāng)源項(xiàng)的剛性比較強(qiáng)時(shí)，通常要對(duì)流動(dòng)項(xiàng)和化學(xué)源項(xiàng)分裂求解，即對(duì)于方程組（10）分別求解流動(dòng)項(xiàng)的半離散方程 （13） 和化學(xué)源項(xiàng)的半離散方程 （14） 本文采用的分裂法是具有二階時(shí)間精度的Strang分裂方法[12] （15） 其中和分別表示以和為時(shí)間步長(zhǎng)求解（14）和（13）的算子。 在分裂方法中，守恒律的半離散方程組（13）仍采用顯式的三階的SSPRK方法求解，而源項(xiàng)的方程組（14）采用常微分方程求解器VODE求解[13]。 1.3 有限體積型WENO方法 本文將比較RKDG方法和高精度的FVWENO方法[14]。在§1.1節(jié)介紹的空間離散中，如果基函數(shù)只取為1，則可得到守恒律方程組（1）的有限體積方法。在二維矩形網(wǎng)格上，半離散有限體積法的方程為 （16） 上式中對(duì)流項(xiàng)的離散采用維數(shù)分裂的基于特征變量的五階WENO格式[14]。在應(yīng)用無(wú)分裂方法計(jì)算時(shí)，由于FVWENO方法的待求變量是流動(dòng)變量的網(wǎng)格單元平均值，如果直接用它計(jì)算（16）式中源項(xiàng)的積分，則只有二階空間精度。為盡量提高精度，本文對(duì)一維問(wèn)題采用四階的Simpson公式來(lái)計(jì)算積分 （17） 其中積分點(diǎn)的值和采用基于特征變量的五階WENO重構(gòu)得到。對(duì)于的重構(gòu)會(huì)產(chǎn)生負(fù)權(quán)，故需要對(duì)重構(gòu)結(jié)果進(jìn)行處理[15]。 對(duì)于二維問(wèn)題，可以采用（17）式的張量積形式，這時(shí)積分具有四階精度，但涉及重構(gòu)積分點(diǎn)處值，計(jì)算公式相對(duì)復(fù)雜。為計(jì)算方便，本文采用一種具有三階精度的二維Simpson積分 （18） 其中由方向的單元平均值重構(gòu)得到，由方向的單元平均值重構(gòu)得到，采用兩個(gè)方向重構(gòu)值的平均值。 2 數(shù)值算例 2.1 精度測(cè)試 第一個(gè)例子是一個(gè)帶有源項(xiàng)的Burgers方程 （19） 這個(gè)方程的解是連續(xù)的，其解析解為 （20） 式中控制方程的剛性，即為初值。為了測(cè)試算法的精度，取，此時(shí)問(wèn)題的剛性不強(qiáng)。計(jì)算區(qū)間取為[-15,25]，RKDG方法中TVB限制器的常數(shù)取為5，計(jì)算到0，CFL數(shù)取為0.18。表1是數(shù)值解的誤差及其數(shù)值精度階。 表1 帶有源項(xiàng)的Burgers方程 Table 1 Burgers equation with source term N Strang split method Unsplit method DG2 R WENO5 R DG2 R WENO5 R WENOM R 32 7.66´10-5 － 4.79´10-3 － 7.55´10-5 － 4.78´10-3 － 5.25´10-3 － 64 5.29´10-6 3.86 9.62´10-4 2.31 4.77´10-6 3.98 9.64´10-4 2.31 5.04´10-4 3.38 128 4.58´10-7 3.53 2.42´10-4 1.99 1.91´10-7 4.64 2.43´10-4 1.99 3.19´10-5 3.98 256 8.37´10-8 2.45 5.76´10-5 2.07 1.36´10-8 3.81 5.76´10-5 2.07 2.41´10-6 3.72 512 1.84´10-8 2.18 1.43´10-5 2.01 1.42´10-9 3.26 1.43´10-5 2.01 2.80´10-7 3.11 1024 4.33´10-9 2.09 3.56´10-6 2.00 1.69´10-10 3.07 3.56´10-6 2.00 3.53´10-8 2.99 2048 1.05´10-9 2.04 8.90´10-7 2.00 2.09´10-11 3.02 8.90´10-7 2.00 4.46´10-9 2.99 注：DG2表示P2基函數(shù)的RKDG方法，WENO5表示五階FVWENO方法， WENOM表示五階FVWENO方法+源項(xiàng)四階精度積分公式（17） 從表1中可以看出，用Strang分裂方法，由于在時(shí)間上只有二階精度，在時(shí)間步長(zhǎng)與空間步長(zhǎng)同階的情況下，即使采用高精度的基函數(shù)DG方法也只能達(dá)到二階精度。而無(wú)分裂的元DG方法和修正FVWENOM方法可以達(dá)到和SSPRK方法同階的三階精度，但沒(méi)有采用高精度源項(xiàng)積分的FVWENO只達(dá)到二階精度，無(wú)論是分裂方法還是無(wú)分裂方法，RKDG方法在網(wǎng)格規(guī)模相同的情況下，其數(shù)值解的誤差明顯要比FVWENO方法的小得多。 2.2 間斷計(jì)算能力測(cè)試 第二個(gè)例子是一個(gè)帶有源項(xiàng)的對(duì)流方程 （21） 取初值為 （22） 是間斷的初始位置，這個(gè)例子的解析解是一個(gè)以速度1向右傳播的間斷解。 文[5]指出，當(dāng)網(wǎng)格數(shù)不夠多時(shí)，計(jì)算的間斷波速度變慢甚至是不動(dòng)的。為了衡量方法的間斷捕捉能力，定義數(shù)值平均速度 其中為一維空間步長(zhǎng)，為計(jì)算中止時(shí)間，和分別是時(shí)刻初值和數(shù)值解的求和。 計(jì)算中取，CFL數(shù)取0.1，計(jì)算區(qū)間為[0,2]，初始間斷位置為。 表2 帶有剛性源項(xiàng)的對(duì)流方程 Table 2 Linear advection equation with stiff source term N Strang split method Unsplit method DG2 WENO5 DG2 WENO5 WENOM 256 5.3´10-4 4.7´10-4 0.7106 4.8´10-4 4.8´10-4 512 0.7249 4.0´10-2 0.9029 7.8´10-4 0.4032 1024 0.9332 0.8145 0.9689 0.8094 0.8327 2048 0.9877 0.9398 0.9956 0.9387 0.9421 從表2中可以看出，在網(wǎng)格規(guī)模相同的情況下，RKDG方法的數(shù)值平均速度比FVWENO的更接近于1，說(shuō)明前者對(duì)于間斷波的捕捉更為準(zhǔn)確。而對(duì)于同一數(shù)值格式，當(dāng)網(wǎng)格數(shù)較少時(shí)，計(jì)算的數(shù)值平均速度是不正確的。只有當(dāng)網(wǎng)格數(shù)足夠多時(shí)，各種格式得到的數(shù)值速度才趨近于1。無(wú)分裂的DG2和WENOM比對(duì)應(yīng)的分裂方法的結(jié)果稍微準(zhǔn)確些，這與無(wú)分裂方法比相應(yīng)的分裂方法的精度高一階有關(guān)。 2.3 一維ZND爆轟波 第三個(gè)算例是一維穩(wěn)定爆轟波問(wèn)題，初始時(shí)刻，假定爆轟波的von Neumann點(diǎn)位于。首先考慮一個(gè)CJ-ZND爆轟波情形。未燃?xì)怏w狀態(tài)為 狀態(tài)方程及Arrhenius模型中各參數(shù)分別為。這些參數(shù)下理論解半反應(yīng)長(zhǎng)度為1[3]。計(jì)算區(qū)間取為[-40,40]，網(wǎng)格數(shù)N=400，該網(wǎng)格數(shù)可以分辨反應(yīng)區(qū)，計(jì)算的CFL取為0.18。該問(wèn)題的理論爆轟波速度為5.4419。計(jì)算到時(shí)間，此時(shí)爆轟波傳播到27.2附近。圖1是計(jì)算的密度分布的局部放大圖，可以看出RKDG方法的解比FVWENO的更接近于精確解，且間斷面過(guò)渡網(wǎng)格點(diǎn)要少。但是由于RKDG方法使用了TVB型的限制器（限制器常數(shù)取為1.0），算出的von Neumann點(diǎn)比FVWENO的略低，而分裂方法對(duì)于解的影響幾乎看不出。 圖1 CJ-ZND爆轟波結(jié)果局部放大圖 Fig.1 Detail view of CJ-ZND detonation wave 圖2 強(qiáng)ZND爆轟波結(jié)果局部放大圖 Fig.2 Detail view of strong ZND detonation wave 其次，我們考慮一個(gè)強(qiáng)爆轟情形，其未燃?xì)怏w狀態(tài)和前面的CJ-ZND爆轟波波前的相同，狀態(tài)方程及Arrhenius模型中各參數(shù)分別為 。這時(shí)半反應(yīng)長(zhǎng)度為0.01[3]，計(jì)算區(qū)間為[-10,110]，網(wǎng)格數(shù)N=800。該網(wǎng)格數(shù)對(duì)這個(gè)問(wèn)題是“under-resolved”，即無(wú)法分辨反應(yīng)過(guò)程中組分的空間變化，只能捕捉爆轟波的位置。該情形下理論強(qiáng)爆轟波的速度為4.9093，計(jì)算到時(shí)間，此時(shí)爆轟波傳播到98.186附近。圖2是局部放大圖，圖中實(shí)線是用4000個(gè)網(wǎng)格求得的理論解?？梢钥闯鯮KDG方法的間斷面過(guò)渡網(wǎng)格點(diǎn)數(shù)比FVWENO的更少，而分裂方法對(duì)于兩種數(shù)值方法的影響不大。 2.4 二維不穩(wěn)定爆轟波 第四個(gè)算例考慮一個(gè)初始擾動(dòng)了的CJ-ZND爆轟波經(jīng)過(guò)一段時(shí)間發(fā)展之后，它的化學(xué)反應(yīng)區(qū)里具有胞格這樣一種規(guī)則的結(jié)構(gòu)。初始時(shí)刻，該ZND爆轟波的von Neumann點(diǎn)位于處，右側(cè)的未燃物狀態(tài)為 狀態(tài)方程及Arrhenius模型中各參數(shù)分。初值的周期擾動(dòng)在方向，初始流動(dòng)狀態(tài)其中為在上述參數(shù)下的精確CJ-ZND解。物理區(qū)域?yàn)?。網(wǎng)格數(shù)為，CFL數(shù)取0.18，上下邊界采用反射邊界條件，左右邊界分別采用入流和出流邊界條件。圖3和圖4分別是元RKDG方法和五階WENO方法結(jié)合Strang分裂法的結(jié)果。由于所用的網(wǎng)格足夠密，二者的結(jié)果差別看不出。在兩種方法得到的圖中，都可以清楚的看到三波點(diǎn)沿著爆轟波陣面的橫向移動(dòng)[17]。 圖3 P2-RKDG結(jié)合Strang分裂模擬二維不穩(wěn)定爆轟波結(jié)果，從左到右、從上到下依次為以Dt=1/30為時(shí)間間隔，從t=3/30到8/30的密度等值線圖 Fig.3 Numerical results of 2d unstable detonation wave by P2-RKDG with Strang split. Sequences from left to right and top to bottom correspond to those from t=3/30 to 8/30 with timestep Dt=1/30 圖4 五階FVWENO結(jié)合Strang分裂模擬二維不穩(wěn)定爆轟波結(jié)果， 從左到右、從上到下依次為以Dt=1/30為時(shí)間間隔，從t=3/30到8/30的密度等值線圖 Fig.4 Numerical results of 2d unstable detonation wave by fifth order FVWENO with Strang split. Sequences from left to right and top to bottom correspond to those from t=3/30 to 8/30 with timestep Dt=1/30 2.5 CJ-ZND爆轟波繞過(guò)拐角的衍射現(xiàn)象 第五個(gè)算例為一個(gè)CJ-ZND爆轟波繞過(guò)一個(gè)的拐角，計(jì)算區(qū)域?yàn)?，其中是一固壁，初始時(shí)刻一向右傳播的CJ-ZND爆轟波位于固壁上面處。右側(cè)的未燃?xì)怏w狀態(tài)為狀態(tài)方程及Arrhenius模型中的參數(shù)分別為 計(jì)算取到時(shí)間，CFL數(shù)取為0.18，網(wǎng)格數(shù)為。上下邊界和壁面采用反射邊界條件，左右邊界分別采用入流和出流邊界條件。 當(dāng)爆轟波傳播到拐角后，支撐爆轟波燃燒的激波在壁面附近的強(qiáng)度減弱，溫度降低，沿著繞射激波，激波后面的溫度將會(huì)明顯地降低到比相應(yīng)ZND波的溫度低，這時(shí)化學(xué)反應(yīng)將會(huì)停止。隨著未反應(yīng)激波垂直于豎直墻的表面繼續(xù)往前傳播，緊跟其后的流體溫度將會(huì)上升，這時(shí)依賴于活化能的強(qiáng)度，如果波后溫度相對(duì)于活化能足夠高，燃料可被重新點(diǎn)燃，爆燃波被重新激發(fā)，并在遠(yuǎn)處與水平向右傳播的爆轟波相連接。圖5(a)和5(b)分別是用階RKDG和五階WENO結(jié)合Strang分裂得到的結(jié)果，可以清楚的看到爆轟波繞過(guò)拐角后爆燃波和前導(dǎo)激波分離的現(xiàn)象。我們將FVWENO在網(wǎng)格上的計(jì)算結(jié)果作為更準(zhǔn)確的參考解，并取的截面進(jìn)行了比較，發(fā)現(xiàn)在網(wǎng)格上RKDG的結(jié)果和參考解的符合程度要略好于FVWENO的結(jié)果。 (a) RKDG方法的密度等值線圖 (a) Density contour of P2 RKDG method (b) FVWENO方法的密度等值線圖 (b) Density contour of 5th order FVWENO 圖5 CJ-ZND爆轟波繞過(guò)90°拐角 Fig. 5 CJ-ZND detonation wave passing a 90° corner 圖6 x=25截面的密度 Fig. 6 Density of slice plane x=25 3 結(jié)論 本文采用龍格庫(kù)塔間斷有限元方法對(duì)帶有源項(xiàng)的守恒律方程組進(jìn)行了數(shù)值求解并和有限體積型WENO方法進(jìn)行了比較。數(shù)值結(jié)果表明，間斷有限元方法不僅在光滑的區(qū)域具有較高精度，而且具有非常好的間斷捕捉能力，因而可以更好的模擬爆轟波的ZND結(jié)構(gòu)。即使對(duì)于強(qiáng)剛性問(wèn)題，由于RKDG法在單元內(nèi)通過(guò)選取適當(dāng)?shù)幕瘮?shù)提高數(shù)值精度，能夠在較少的網(wǎng)格數(shù)下得到正確的爆轟波位置，克服了一般有限體積方法在求解強(qiáng)剛性問(wèn)題時(shí)由于網(wǎng)格分辨率不夠而容易產(chǎn)生虛假反應(yīng)的現(xiàn)象，而且由于RKDG方法不要求近似解的全局連續(xù)性，因此非常易于自適應(yīng)和并行計(jì)算，有望成為模擬復(fù)雜化學(xué)反應(yīng)流動(dòng)的一種有效方法。 參考文獻(xiàn) [1] Tosatto L, Vigevano L. Numerical solution of under-resolved detonations [J]. J Comput Phys, 2008, 227: 2317-2343. [2] Helzel C, LeVeque R, Warneke G. A modified fractional step method for the accurate approximation of detonation waves [J], SIAM J Sci Comput, 1999, 22: 1489-1510. [3] Berkenbosch A C. Capturing detonation waves for reactive Euler equations: [Ph.D. thesis][D]. Edndhoven, The Netherlands. Technische Universiteit Eindhoven. 1995. [4] Colella P, Majda A, Roytburd V. Theoretical and numerical structure for reacting shock waves [J]. SIAM J Sci Stat Comput, 1986, 7: 1059-1080. [5] LeVeque R J, Yee H C. A study of numerical methods for hyperbolic conservation laws with stiff source term [J]. J Comput Phys, 1990, 86: 187-210. [6] Oran E S, Boris J P. Numerical simulation of reactive flow [M]. Elsevier, New York, 1987, Chap.13. [7] Yuan L, Tang T. Resolving the shock-induced combustion by an adaptive mesh redistribution method [J]. J Comput Phys, 2007, 224: 587-600. [8] Azarenok B, Tang T. Second-order Godunov scheme for reactive flow calculations on moving meshes [J]. J Comput Phys, 2005, 206: 48-80. [9] Deiterding R. Parallel adaptive simulation of multi-dimensional detonation structures. PhD thesis [D], Brandenburgische Technische Universität Cottbus, Sep. 2003. [10] Cockburn B. Discontinuous Galerkin methods for convection-dominated problem [A]. //Barth T J, Deconinck H. High-Order Methods for Computational Physics, Vol. 9, LNCSE[C], Berlin: Springer, 1999: 69-224. [11] Gottlieb S, Shu C W, Tadmor E. Strong stability-preserving high-order time discretization methods [J]. SIAM REVIEW, 2001, 43: 89-112. [12] LeVeque R J. Finite volume methods for hyperbolic problems [M], Cambridge: Cambridge University Press, 2002. [13] Brown P N, Byrne G D, Hindmarsh A C. VODE: A variable coefficient ODE solver [J]. SIAM J. Sci. Stat. Comput, 1989, 10: 1038-1051. [14] Shu C W. High order ENO and WENO schemes [A]. //Barth T J, Deconinck H. High-order Methods for Computational Physics, Vol. 9, LNCSE[C], Berlin: Springer, 1999: 439-582. [15] Shi J, Hu C Q, Shu C W. A technique of treating negative weights in WENO schemes [J]. J Comput Phys, 2002, 175: 108-127. [16] Bao W, Jin S. The random projection method for stiff detonation capturing [J]. SIAM J Sci Comput, 2001, 22: 1000-1026. [17] Lian Y S, Xu K. A gas-kinetic scheme for multimaterial flows and its application in chemical reactions [J]. J Comput Phys, 2000, 163: 349-375. 作者簡(jiǎn)介： 1. 張磊(1977-)，男，北京，講師，在讀博士生，主要從事間斷有限元方法及在流動(dòng)計(jì)算中應(yīng)用的研究。通信地址：中國(guó)礦業(yè)大學(xué)(北京)理學(xué)院，北京100083 第一作者照片： 其他： ------------------------------------------ 身份證號(hào)，130603197704170319 13