王志剛

(阜陽師范學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，安徽 阜陽 236037)

基于動力學(xué)關(guān)系的非壓縮激波的Godunov型格式

王志剛

(阜陽師范學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，安徽 阜陽 236037)

非壓縮激波產(chǎn)生于大量的物理問題，如相變動力學(xué)、磁流體動力學(xué)、Camassa-Holm模型和燃燒系統(tǒng)等。與壓縮激波不同的是，它不僅滿足單個的熵不等式，還能滿足一些動力學(xué)關(guān)系。為了計(jì)算非壓縮激波的數(shù)值解，本文設(shè)計(jì)了一種Godunov型格式,包括函數(shù)重構(gòu)、發(fā)展和求網(wǎng)格平均三個步驟。在函數(shù)重構(gòu)時，先對非壓縮激波的位置進(jìn)行預(yù)估，在其相鄰網(wǎng)格上利用動力學(xué)關(guān)系進(jìn)行重構(gòu)，而在其它網(wǎng)格上采用數(shù)值解和數(shù)值熵進(jìn)行重構(gòu)。數(shù)值實(shí)驗(yàn)表明，此格式不僅對非壓縮激波有較好的分辨率，而且對經(jīng)典波也有較高的精度。

單個守恒律方程；動力學(xué)關(guān)系；壓縮激波；非壓縮激波；Godunov格式

單個守恒律方程是一類非常重要的發(fā)展方程，在物理學(xué)和力學(xué)等方面具有廣泛的應(yīng)用，一般可表示為

其中u=u(t,x)為未知函數(shù)，f(u)為流函數(shù)。無論初值多么光滑，方程(1)的解都可能在有限時間內(nèi)出現(xiàn)間斷（在物理上被稱為激波或接觸間斷）。因此，方程(1)的解應(yīng)被理解為分布意義上的弱解。然而，弱解又不能被初值和邊界條件唯一確定。為了確定唯一的物理解，Lax[1]，Oleinik[2]，Kruzkov[3]等引入了熵條件，即方程(1)的物理解還應(yīng)滿足如下的熵條件

對所有的熵-熵流對(η(u),F(u))成立。其中，η(u)為凸函數(shù)，F(xiàn)′(u)=η′(u)f′(u)。在條件(2)的限制下，方程(1)的一個間斷解

滿足條件(4)的間斷解被稱為壓縮激波。但是，當(dāng)流函數(shù)為非凸非凹函數(shù)時，在相變動力學(xué)、磁流體動力學(xué)、Camassa-Holm模型和燃燒系統(tǒng)等實(shí)際問題中，還存在一類重要的間斷解，它僅滿足如下的熵條件

對單個的熵-熵流對(η(u),F(u))成立，而不再滿足條件(4)，故稱該類間斷解為非壓縮激波。為了確保解的唯一性，非壓縮激波還必須滿足一種動力學(xué)關(guān)系[4]，即

其中?是動力學(xué)函數(shù)，?-1為?的反函數(shù)。

壓縮激波解的數(shù)值模擬已受到計(jì)算專家們的熱切關(guān)注，設(shè)計(jì)出了多種計(jì)算方法，例如：有限差分法、有限體積法和不連續(xù)的有限元法等。與壓縮激波不同，非壓縮激波不再滿足解的單調(diào)性和解的有限變差不增性[4]。因此，這些經(jīng)典的數(shù)值方法很難成功應(yīng)用到非壓縮激波。直到現(xiàn)在，非壓縮激波的數(shù)值模擬仍然具有較大的挑戰(zhàn)性。非壓縮激波解數(shù)值模擬的關(guān)鍵是動力學(xué)關(guān)系的離散，目前主要有兩種方法：擴(kuò)散界面法和自由界面法。擴(kuò)散界面法不直接離散動力學(xué)關(guān)系，而是針對方程(1)的原始模型——帶粘性項(xiàng)和散射項(xiàng)的高階方程，通過數(shù)值粘性，設(shè)計(jì)與之等價的數(shù)值格式。這類方法對小強(qiáng)度的激波是非常有效的，但也帶來了一些非物理性震蕩[5-9]。自由界面法是在激波跟蹤法或Glimm格式的框架下，直接使用動力學(xué)關(guān)系進(jìn)行離散。這類格式對滿足動力學(xué)關(guān)系的非經(jīng)典激波具有較高的分辨率，但對黎曼解的要求過高，不利于復(fù)雜問題的計(jì)算[10-13]。

最近，Lagoutiere等基于單個非壓縮激波的傳播特征，提出了一種新的Godunov型格式[14]。為了分辨出非經(jīng)典激波，(i)當(dāng)和時，并不一定會產(chǎn)生非經(jīng)典激波，故格式也會計(jì)算出一些錯誤的數(shù)值解；(ii)由于是變化的，為了便于實(shí)現(xiàn)，此格式要求流函數(shù)必須為單調(diào)函數(shù)；(iii)此格式對經(jīng)典激波和疏散波僅有一階的精度。

綜上，本文將針對非壓縮激波設(shè)計(jì)一種新的Godunov型格式。為了捕捉非壓縮激波，此格式在重構(gòu)函數(shù)之前，先預(yù)估非經(jīng)典激波的位置，僅在非經(jīng)典激波處使用動力學(xué)關(guān)系進(jìn)行函數(shù)重構(gòu)，并在相鄰兩個網(wǎng)格采用分片常數(shù)重構(gòu)，而在其它網(wǎng)格采用文獻(xiàn)[15]的重構(gòu)方法，以此提高對經(jīng)典波的近似精度。數(shù)值實(shí)驗(yàn)表明：此格式對非壓縮激波具有較好的分辯率，對經(jīng)典波也具有較高的精度，且對非單調(diào)的凹凸流函數(shù)也有效。

1 非經(jīng)典Riemann解

本部分將介紹問題(1)-(5)-(6)的Riemann解。為了簡單，僅考慮流函數(shù)為凹凸函數(shù)，即

例如：f(u)=u3+au為凹凸函數(shù)。首先，定義一些函數(shù)。令函數(shù)?#∶R→R，滿足

記其反函數(shù)為 ?-#∶R→R 。令函數(shù) ?0∶R→R滿足

問題(1)-(5)-(6)的非經(jīng)典Riemann解如下：

Ⅰ.當(dāng) ul＞0，

ⅰ.如果ur≥ul，則解為連接ul和ur的一個疏散波；

ⅳ.如果 ur≥?(ul)，則解為連接 ul到 ?(ul)的一個非經(jīng)典激波與連接?(ul)到ur的一個疏散波。

Ⅱ.當(dāng)ul≤0，

ⅰ.如果ur≤ul，則解為連接ul和ur的一個疏散波；

ⅲ.如果 ur∈[?T(ul),?(ul))，則解為連接 ul到?(ul)的一個非經(jīng)典激波與連接?(ul)到ur的一個經(jīng)典的激波；

ⅳ.如果 ur≥?(ul)，則解為連接 ul到 ?(ul)的一個非經(jīng)典激波與連接?(ul)到ur的一個疏散波。

2 數(shù)值格式

本部分主要研究非壓縮激波的一種新的Godunov型格式。為了捕捉非壓縮激波，此格式在重構(gòu)函數(shù)之前，先預(yù)估非經(jīng)典激波的位置，僅在非經(jīng)典激波處使用動力學(xué)關(guān)系進(jìn)行函數(shù)重構(gòu)，并在相鄰兩個網(wǎng)格采用分片常數(shù)重構(gòu)，而在其它網(wǎng)格采用文獻(xiàn)[15]的重構(gòu)方法，以此提高對經(jīng)典波的近似精度。

令h和τ為空間步長和時間步長，且xj=jh,xj±1/2=(j±1/2)h,tn=nτ,其 中 j∈Z ，n∈Z+。數(shù)值解和數(shù)值熵是對解u(x,tn)和熵U(u(x,tn))在[xj-1/2,xj+1/2)上網(wǎng)格平均的近似，即

數(shù)值格式為Godunov格式，包括函數(shù)重構(gòu)、發(fā)展和求網(wǎng)格平均三個步驟。

2.1 函數(shù)重構(gòu)

為了捕捉非壓縮激波，先定義經(jīng)典網(wǎng)格和非經(jīng)典網(wǎng)格。

定義1如果

則稱網(wǎng)格[xj-1/2,xj+1/2)為非經(jīng)典網(wǎng)格，其它網(wǎng)格都稱為經(jīng)典網(wǎng)格。

注1一個非經(jīng)典壓縮波離散化后，網(wǎng)格將會有兩種情形：一個非經(jīng)典網(wǎng)格或兩個相鄰的非經(jīng)典網(wǎng)格。對于兩個非經(jīng)典網(wǎng)格的情形，當(dāng)非經(jīng)典壓縮波向右傳播時，僅認(rèn)定左邊的網(wǎng)格為非經(jīng)典的；當(dāng)非經(jīng)典壓縮波向左傳播時，僅認(rèn)定右邊的網(wǎng)格為非經(jīng)典的。

▲非經(jīng)典網(wǎng)格及其相鄰網(wǎng)格

若[xj-1/2,xj+1/2)為非經(jīng)典網(wǎng)格，則[xj-3/2,xj-1/2)、[xj-1/2,xj+1/2)和[xj+1/2,xj+3/2)上的重構(gòu)函數(shù)為

▲其它網(wǎng)格

[xj-1/2,xj+1/2)為經(jīng)典網(wǎng)格，且相鄰的網(wǎng)格也是經(jīng)典的。采用文獻(xiàn)[15]中的重構(gòu)方法，即

2.2 發(fā)展

求Cauchy問題

該問題為包含一族中心在xj，xˉj和xj+1/2的黎曼問題。為了避免波的相互作用，選取網(wǎng)格步長比滿足

對于非經(jīng)典網(wǎng)格，以xˉj為中心的波還會影響到邊界上。但是，由于重構(gòu)函數(shù)(13)的特點(diǎn)，邊界上的數(shù)值流函數(shù)是可計(jì)算的，確保了格式的可操作性。

2.3 求網(wǎng)格平均

計(jì)算t=tn+1時刻的數(shù)值解

其中，vl和vr為激波的左右狀態(tài)，s為激波速度。

3 數(shù)值實(shí)驗(yàn)

本部分將通過兩個數(shù)值例子驗(yàn)證格式的有效性。為了比較，也使用文獻(xiàn)[14]的格式進(jìn)行計(jì)算。本文數(shù)值格式計(jì)算的數(shù)值解用“new”表示，而文獻(xiàn)[14]的格式計(jì)算的數(shù)值解用“Lagoutiere”表示。

算例1考慮凹凸流 f(u)=u3+u，熵函數(shù)和熵流函數(shù)為：

選取動力學(xué)函數(shù)：?(u)=-βu,(β=0.75)，則?T(u)=-u-?(u)。選取λ=0.005,h=0.02。初值Ⅰ：生成單個非壓縮波的初值

計(jì)算t=0.03時刻的數(shù)值解，見圖1。

初值Ⅱ：生成一個非壓縮波和一個壓縮波的初值

計(jì)算t=0.03時刻的數(shù)值解，見圖2。

初值Ⅲ：生成一個壓縮波和一個疏散波的初值

計(jì)算t=0.01時刻的數(shù)值解，見圖3。初值Ⅳ：正弦函數(shù)初值

計(jì)算t=0.5時刻的數(shù)值解，見圖4。

從算例1的計(jì)算結(jié)果可看出：對非壓縮激波，本文設(shè)計(jì)的Godunov型格式和文獻(xiàn)[14]的格式具有幾乎相同的分辨率；而對經(jīng)典波，本文設(shè)計(jì)的格式的近似精度比文獻(xiàn)[14]的格式要好一些。

圖1 單個非經(jīng)典激波

圖2 非經(jīng)典激波+經(jīng)典激波

圖3 非經(jīng)典激波+疏散波

圖4 正弦初值問題在t=0.5時刻的數(shù)值解

算例2考慮凹凸流 f(u)=u3-15u，熵函數(shù)和熵流函數(shù)為：

初值Ⅰ：生成一個非壓縮波和一個疏散波的初值

計(jì)算t=0.03時刻的數(shù)值解，見圖5。

初值Ⅱ：一個壓縮和一個非壓縮波的相互作用

計(jì)算t=0.002,0.0028和0.1時的數(shù)值解，見圖6。

算例2中的流函數(shù)是一個非單調(diào)函數(shù)，故波的傳播方向既有向左的，也有向右的，有時也可能相互作用。初值Ⅰ產(chǎn)生的波一個向左傳播，一個向右傳播，而初值Ⅱ產(chǎn)生的波相互作用后又產(chǎn)生一個向左的非壓縮波和一個向右的疏散波。從數(shù)值結(jié)果可看出：本文設(shè)計(jì)的格式對非單調(diào)流函數(shù)的情形仍然有效，彌補(bǔ)了文獻(xiàn)[14]中格式的不足。

圖5 非壓縮波+疏散波

圖6 壓縮波和非壓縮波的相互作用

[1]Lax P D.Hyperbolic systems of conservation laws and the mathematical theory of shock waves[C].Regional Confer.,1973.

[2]Oleinik O.Discontinuous solutions of nonlinear differential equations[J].Transl.Amer.Math.Soc.,1963,26:95-172.

[3]Kruzkov S N,First order quasilinear equations with several independent variables[J].Mat.Sb.N.S.,1970,81:228-235.

[4]Lefloch P G.Hyperbolic systems of conservation laws:the theory of classical and nonclassical shock waves[M].LectureNotesinMathematics,ETH Zurich,Birkhauser,2002.

[5]Hayes B T,LeFloch P G.Nonclassical shocks and kinetic relations: fi nite difference schemes[J].SIAM J.Numer.Anal.,1998,35:2169-2194.

[6]LeFloch P G,Rohde C.High-order schemes,entropy inequalities,and nonclassical shocks[J].SIAM J.Numer.Anal.,2000,37:2023-2060.

[7]Chalons C,LeFloch P G.A fully discrete scheme for diffusive-dispersive conservation laws[J].Numerisch Math.,2001,89:493-509.

[8]Chalons C,LeFloch P G.High-order entropy conservative schemes and kinetic relations for van der waals fl uids[J].J Comput.Phys.,2001,167:1-23.

[9]Ernest J,LeFloch P G,Mishra S.Schemes with Well-Controlled Dissipation[J].SIAM J.Numer.Anal.,2015,53(1):674-699.

[10]Hou T Y,LeFloch P G,Zhong X.Converging methods for the computation of propagating solid-solid phase boundaries[J].J.Comput.Phys.,1996,124:192-216.

[11]Amadori D,Baiti P,LeFloch P G,Piccoli B.Nonclassical shocks and the Cauchy problem for non-convex conservation laws[J].J.Differential Equations,1999,151:345-372

[12]Hou T Y,LeFloch P G,Rosakis P.A level-set approach to the computation of twinning and phase transition dynamics[J].J Comput.Phys.,1999,150:302-331.

[13]Chalons C,LeFloch P G.Computing undercompressive waves with the random choice scheme[J].Nonclassical shock waves,Interfaces Free Boundaries,2003,5:129-158.

[14]Boutin B,Chalons C,Lagoutiere F A.Convergent and conservative schemes for nonclassical solutions based on kinetic relations[J].Interfaces and Free Boundaries,2008,10(3):399-421.

[15]Chen R S,Mao D K.Entropy-TVD scheme for nonlinear scalar conservation Laws[J].Journal of Scientific Computing,2011,47(2):150-169.

A Godunov-type scheme computing undercompressive shock waves based on kinetic relations

WANG Zhi-gang
(School of Mathematics and Statistics,Fuyang Normal University,Fuyang Anhui 236037,China)

Undercompressive shock waves arise from many physical problems:phase transitions,magnetohydrodynamics,Camassa-Holm model,combustion theory,etc.Different from compressive shock waves,it satisfies single entropy inequality and some kinetic relations.In order to computing undercompressive shock waves,we design a Godunov-type scheme,including reconstruction,evolution,cell-averaging.We forecast the position of undercompressive shocks and use kinetic relations to reconstruct functions near undercompressive shocks.In other cells,we use numerical solutions and numerical entropy to reconstruct functions.Numerical experiments show that our scheme has not only good resolution for undercompressive shocks,but also better accuracy for classical waves.

single conservation law;kinetic relations;compressive shocks;undercompressive shocks;Godunov scheme

O241.82

A

1004-4329（2017）02-001-05

10.14096/j.cnki.cn34-1069/n/1004-4329（2017）02-001-05

2017-03-20

國家自然科學(xué)基金（11401104）；中國博士后基金（2015M581579）；阜陽師范學(xué)院博士科研基金（FSB201201010）資助。

王志剛（1979－ ），男，博士，副教授，研究方向：守恒律方程的數(shù)值模擬。