基于動力學關系的非壓縮激波的Godunov型格式

王志剛

(阜陽師范學院 數學與統計學院，安徽 阜陽 236037)

基于動力學關系的非壓縮激波的Godunov型格式

王志剛

(阜陽師范學院 數學與統計學院，安徽 阜陽 236037)

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

單個守恒律方程；動力學關系；壓縮激波；非壓縮激波；Godunov格式

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

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

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

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

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

其中?是動力學函數，?-1為?的反函數。

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

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

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

1 非經典Riemann解

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

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

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

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

Ⅰ.當 ul＞0，

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

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

Ⅱ.當ul≤0，

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

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

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

2 數值格式

本部分主要研究非壓縮激波的一種新的Godunov型格式。為了捕捉非壓縮激波，此格式在重構函數之前，先預估非經典激波的位置，僅在非經典激波處使用動力學關系進行函數重構，并在相鄰兩個網格采用分片常數重構，而在其它網格采用文獻[15]的重構方法，以此提高對經典波的近似精度。

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

數值格式為Godunov格式，包括函數重構、發展和求網格平均三個步驟。

2.1 函數重構

為了捕捉非壓縮激波，先定義經典網格和非經典網格。

定義1如果

則稱網格[xj-1/2,xj+1/2)為非經典網格，其它網格都稱為經典網格。

注1一個非經典壓縮波離散化后，網格將會有兩種情形：一個非經典網格或兩個相鄰的非經典網格。對于兩個非經典網格的情形，當非經典壓縮波向右傳播時，僅認定左邊的網格為非經典的；當非經典壓縮波向左傳播時，僅認定右邊的網格為非經典的。

▲非經典網格及其相鄰網格

若[xj-1/2,xj+1/2)為非經典網格，則[xj-3/2,xj-1/2)、[xj-1/2,xj+1/2)和[xj+1/2,xj+3/2)上的重構函數為

▲其它網格

[xj-1/2,xj+1/2)為經典網格，且相鄰的網格也是經典的。采用文獻[15]中的重構方法，即

2.2 發展

求Cauchy問題

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

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

2.3 求網格平均

計算t=tn+1時刻的數值解

其中，vl和vr為激波的左右狀態，s為激波速度。

3 數值實驗

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

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

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

計算t=0.03時刻的數值解，見圖1。

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

計算t=0.03時刻的數值解，見圖2。

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

計算t=0.01時刻的數值解，見圖3。初值Ⅳ：正弦函數初值

計算t=0.5時刻的數值解，見圖4。

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

圖1 單個非經典激波

圖2 非經典激波+經典激波

圖3 非經典激波+疏散波

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

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

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

計算t=0.03時刻的數值解，見圖5。

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

計算t=0.002,0.0028和0.1時的數值解，見圖6。

算例2中的流函數是一個非單調函數，故波的傳播方向既有向左的，也有向右的，有時也可能相互作用。初值Ⅰ產生的波一個向左傳播，一個向右傳播，而初值Ⅱ產生的波相互作用后又產生一個向左的非壓縮波和一個向右的疏散波。從數值結果可看出：本文設計的格式對非單調流函數的情形仍然有效，彌補了文獻[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

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

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