求解雙曲守恒律方程的高分辨率熵相容格式

任 炯， 封建湖， 劉友瓊， 梁 楠

（1.長安大學理學院， 陜西西安 710064；2.西北工業大學航空學院流體力學系， 陜西西安 710072）

求解雙曲守恒律方程的高分辨率熵相容格式

任 炯1，2， 封建湖1，*， 劉友瓊1， 梁 楠1

（1.長安大學理學院， 陜西西安 710064；2.西北工業大學航空學院流體力學系， 陜西西安 710072）

為提高熵相容格式的精度，利用限制器機制構造高分辨率格式，將構造的通量限制器插入熵相容格式，得到一類高分辨率熵相容格式.構造Euler方程高分辨率熵相容格式時，對熵相容格式中的幾個參數做簡單調整，提高了接觸間斷處的分辨率.將所得格式的數值結果與熵相容格式的數值結果比較表明，構造的高分辨率熵相容格式具有穩健和基本無振蕩等特性.

雙曲守恒律；熵相容格式；限制器；高分辨率

0 引言

在一維情況下，雙曲守恒律方程的一般形式為

其中（x，t）∈R×R＋，u：R×R＋→RN（N≥1）是守恒變量，f（u）是通量函數.該類方程作為流體力學中重要的物理模型，在航空航天、氣象、海洋等科學工程領域都有廣泛的應用，其數值求解的方法一直是國內外眾多學者致力于研究的方向.對于非線性雙曲守恒律方程，即使所給的初始條件光滑，其解在時間推進的某個時刻也可能會出現間斷，產生激波.間斷解的出現使雙曲守恒律方程的解違背了古典解的理論，于是Lax在1954年提出了弱解［1］的概念，允許間斷解存在，但弱解不唯一，這就需要從方程的實際物理背景出發給出限制條件來確定物理相關的解.這個工作也是Lax完成的，他從物理學中的熱力學第二定律出發，提出：如果弱解u滿足熵穩定條件［2］

則u就是唯一的具有物理意義的解，其中E（u）是u的凸函數，滿足E″（u）＞0，被稱為熵函數，F（u）稱為熵通量函數．滿足式（2）的E（u）和F（u）被稱為熵對.為了滿足熵穩定條件，Tadmor首先研究了熵穩定和數值粘性之間的關系，在文獻［3］和［4］中通過引入熵變量和熵勢的概念，構造了一類二階熵守恒格式，該格式的數值通量保持總熵不變，適用于光滑解的計算，但是當有間斷解出現時，由于熵守恒格式缺乏熵的耗散機制，使其數值解表現了嚴重的不穩定性，這個特點可由第3節中的數值算例直觀說明.在此基礎上Tadmor又提出：一個三點格式只需含有比熵守恒格式更多的粘性則是熵穩定的.所以要達到熵穩定，只需在該熵守恒格式的基礎上適當地加入一些粘性.但究竟應該添加多少數值粘性，一直以來都是相關研究人員思索的問題，其中Tadmor，Zhong［5－6］和Roe都做了相應的工作，但以Roe的方法最為理想和實用，Roe是在熵守恒格式的基礎上加上他提出的Roe格式的數值粘性，獲得了一階精度的熵穩定ERoe格式［7］，隨后在2009年，他和Ismail進一步通過分析得到：解在跨過激波時產生了激波強度立方倍的熵增，在此基礎上發展了對數值粘性項更精確量化的熵相容格式［8］，但由于其數值粘性項只有一階精度，從而導致該格式比熵守恒格式的精度有所降低.鑒于此，本文采用傳統的構造二階總變差減少（total variation diminishing，簡稱TVD）格式的方法［9－10］，即利用通量限制器機制提高熵相容格式的精度，達到高分辨率的要求.本文通過推廣常用的Superbee限制器構造了新的通量限制器，并利用該限制器構造得到高分辨率熵相容格式.由于構造的通量限制器的曲線落在二階TVD區域內，且過（1，1）點（見圖1（a）），所以能夠保證新構造的高分辨率熵相容格式滿足二階精度，并且在限制器的開關調節作用下，該格式能夠達到自適應性：在解的光滑區域使用二階精度熵守恒格式，同時避免非物理解的產生；而在間斷區域采用熵相容格式，保證有足夠的數值粘性抑制振蕩.另外本文在構造Euler方程高分辨率熵相容格式時，對其相應的熵相容格式中的幾個參數做了簡單調整，以改善接觸間斷處的捕捉效果.最后在第3節通過幾個數值算例直觀形象地說明了新格式具有穩健性、高精度性和無振蕩性等特點.

本文采用均勻網格上的半離散有限體積格式

用Δx表示空間步長，Δt表示時間步長，由于該格式與時間相關，所以時間方向的離散可以采用高精度的方法.數值算例在時間演化上采用三階Runge-Kutta方法［11］：

1 通量限制器的構造

對原始的Superbee限制器［12］

式（3）中的α，β，γ這三個參數在各自的范圍內可以連續變化，從而產生了一類限制器，我們將此類限制器稱為廣義Superbee限制器，簡記為GSbee類限制器.顯然，這類限制器都在Sweby［14］給出的限制器的二階TVD區域（見圖1（a））內，且都滿足φ（1）＝1（即過（1，1）點），能達到二階的要求，隨著α，β，γ的連續變化，GSbee類限制器覆蓋了整個二階TVD區域.將這類限制器應用于下節將研究的標量高分辨率熵相容格式在一維Burgers方程間斷初值問題上進行大量的數值試驗后，發現對于該問題，參數α，β，γ的變化使Gsbee類限制器的功效表現出某些規律：

1）當固定β和γ，α變化：α越大限制器功效越好；

2）當固定α和γ，β變化：β越大限制器功效越好；

3）當固定α和β，γ變化：γ越小限制器功效越好；

為簡便，這個試驗過程就不一一展示了，而在此，基于以上試驗結論的啟發，我們分別取α，β，γ的極限，即α＝1，β＝2，γ＝1，得到一個新的限制器：

這個限制器可以簡化為

此時，該限制器的形式類似于Minmod［15］限制器

且容易看出Superbee和Minmod限制器都是GSbee類限制器的特殊情況，即

當α＝1，β＝2，γ＝2時，是Superbee限制器；

當α＝1，β＝1，γ＝1時，是Minmod限制器.

在二階TVD區域內，當0＜θ≤1時，新的限制器是Superbee限制器的部分，當θ＞1時，是Minmod限制器的部分，所以將該新的限制器稱為S-M限制器（其中S是Superbee的首字母，M是Minmod的首字母）.Superbee限制器、Minmod限制器和S-M限制器各自在二階TVD區域中的曲線分別如圖1（b）－（d）中的粗黑線所示.

圖1 二階TVD區域和限制器圖示Fig.1 Second order TVD regions and limiters

從圖1（d）可以看到S-M限制器在二階TVD區域的邊界上，且通過點（1，1），能夠保持格式的二階精度.為方便于下節構造高分辨率熵相容格式和在編程實現時減少邏輯語句的輸入，我們將S-M限制器用下面的形式表示

2 高分辨率熵相容格式

2.1 一維標量守恒律方程

2.1.1 熵守恒格式

的積分形式［3－4］

Tadmor同時證明了熵守恒格式具有二階精度.

2.1.2 熵相容格式

熵守恒格式保持了總熵不變，而要滿足熵穩定，總熵必須有所耗散，Ismail［8］通過在熵守恒通量的基礎

Tadmor在文獻［3－4］中由該積分形式出發推導出數值通量滿足熵守恒的條件

為了更好地推廣到方程組的情況，（8）式也可以寫成如下形式［8］

其中，符號［·］＝（·）j＋1－（·）j，ˉa＝（aj＋aj＋1）／2，aj＝f′（uj）．但是這個迎風項的耗散量卻不足以抵消解在跨過激波時所產生的激波強度立方倍的熵增，于是Ismail在上述熵穩定格式（9）的迎風數值粘性項的平均特征速度中補充了一個與特征速度差分的絕對值成比例的量來使熵的耗散更精確，從而得到了熵相容通量

其中，上角標E為熵Entropy的縮寫，C為相容Consistent的縮寫，α＝1／6，熵相容格式是目前對熵的變化估計得最精確的一種熵穩定格式，但是由于其數值粘性項只有一階精度，所以導致熵相容格式的精度比熵守恒格式有所降低.

2.1.3 高分辨率熵相容格式

這雖然不影響格式的自適應性，但可能影響限制器的作用范圍，幸運的是我們在第2節構造的S-M限制器（4）滿足φ（θj＋1／2）≤1，所以完全可以不用考慮加絕對值的問題，因此本文將采用格式（12）進行數值試驗，并將該格式稱為ECL格式.需要進一步說明的是該格式是在原始熵相容格式的基礎上通過限制器的開關作用合理地調節解在各個區域上的熵的耗散量，并沒有改變熵的耗散方向.目前還無法從理論上嚴格證明格式（12）和（13）的穩定性，但是數值算例中并沒有出現違反熵條件的現象.另外需要注意的是，文獻［10，14］中基于一階迎風格式和二階Lax-Wendroff格式構造高分辨率格式時，提到：為盡可能減少格式的數值耗散，提高對間斷的分辨率，需要最大化含有限制器的反擴散通量（即限制器曲線需盡可能靠近二階TVD區域的上邊界），而從（13）式看到，本文構造的高分辨率熵相容格式中，含有限制器的通量是格式的耗散通量，所以在滿足TVD條件的情況下，若要減少格式的耗散，就應該最小化耗散通量（即限制器曲線需盡可能靠近二階TVD區域的下邊界）.基于此分析，再結合圖1（d）中S-M限制器的曲線，我們認為：利用S-M限制器構造高分辨率熵相容格式是合理的.第3節的數值算例部分進一步對此估計做了數值上的驗證.

對于一維Burgers方程

2.2 一維Euler方程

考慮氣體動力學Euler方程

其中u＝［ρ，ρu，E］T是守恒型向量，f（u）＝［ρu，ρu2＋p，u（E＋p）］T是通量函數．ρ，u，p和E分別為密度、速度、壓強和總能，狀態方程為p＝（γ－1）（E－（ρu2）／2），γ＝1.4是比熱比，是聲速.選用Euler方程的熵對E（u）＝－ρS和F（u）＝－mS，其中m＝ρu是動量，S＝ln（pρ－γ）是Euler方程的熵，熵變量為

2.2.1 熵守恒格式

本文采用Roe的方法［8］得到Euler方程的熵守恒數值通量

關于對數平均的計算詳見文獻［8］.

2.2.2 熵相容格式

與標量情形類似，Ismail在Euler方程的熵穩定通量［7］本文將在下小節采用此通量進行高分辨率格式的構造，其中是Euler方程的右特征向量矩陣（詳見文獻［8］）.

2.2.3 高分辨率熵相容格式

3 數值算例

我們以熵守恒和熵相容格式的數值結果作為參照來說明新的高分辨率熵相容格式的特性；同時也說明我們在第2節構造的S-M限制器（4）可以用來構造高分辨率熵相容格式.

符號約定：

C：熵守恒格式；EC：標量熵相容格式（10）；EC2：Euler方程組熵相容格式（15）；ECL：標量高分辨率熵相容格式（12）；EC2L：Euler方程組高分辨率熵相容格式（17）；Exact：每個算例的精確解.

3.1 標量數值試驗

3.1.1 一維Burgers方程連續初值問題（精度測試）

考慮一維Burgers方程的初值問題，在區域［－2，2］上定義初始條件為

取u0＝0，u1＝0.5，CFL＝0.4，空間網格數為N＝40，采用周期邊界條件，圖2（a）和（b）分別給出了時間t＝0.32和t＝0.96的數值結果圖，且每個圖中都分別畫出了熵相容格式EC和高分辨率熵相容格式ECL的數值解以及精確解，從這兩幅圖我們可以看到，在解的光滑區域，兩種格式對解的捕捉效果相似，但當解的梯度變大時，ECL格式相對EC格式表現出了更銳利的效果，并且在表1中，通過兩種格式精度的比較，也說明了通過通量限制器構造的高分辨率熵相容格式在一定程度上能夠提高熵相容格式的精度.

圖2 Burgers方程連續初值問題的數值結果Fig.2 Numerical results of continuous IVP of Burgers equation（N＝40）

表1 初值問題3.1.1在時間t＝0.32和t＝0.96的誤差的L1－范數Table 1 Initial value problem 3.1.1，L1-norms of errors at t＝0.32 and t＝0.96

3.1.2 一維Burgers方程間斷初值問題

在區域x∈［－1，1］上定義初始條件

采用周期邊界條件，取CFL＝0.4，空間網格數為N＝50，計算到時間t＝0.3．這個問題的精確解主要包括兩部分：左邊由－1到1產生了一個稀疏波，右邊由1到－1產生了一個定常激波.對本問題，我們采用了熵守恒C格式、熵相容EC格式和高分辨率熵相容ECL格式進行數值試驗，其結果見圖3，圖3（a）為C的結果，圖3（b）為EC的結果，圖3（c）為ECL的結果，可以看到在x＝1／3的激波位置，熵守恒格式表現了較強的色散效應（即產生了強烈的振蕩），這是由于熵守恒格式沒有任何耗散機制的緣故；而EC格式對此現象進行了完美的修正，這正是熵相容格式在熵守恒格式基礎上添加合適數量的數值粘性所起的作用，只是由于數值粘性項只有一階精度，所以使得EC格式的精度相對C格式有所降低，也使得稀疏波段抹平較為嚴重；基于前面兩種格式的特點，本文構造的ECL格式在稀疏波段和波頭、波尾處都比EC有很大的改善，且同時保持了EC格式對激波的良好捕捉效果，無振蕩的產生，滿足高分辨率的特點.通過觀察圖3 （d）中三種格式的總熵ΔxΣiu2i／2隨時間的變化及其和準確解總熵變化的比較，可以看到C格式的總熵不變，保持熵守恒，而EC和ECL的總熵都是耗散的，保證了熵穩定，且其耗散都多于精確解的耗散，所以都有抹平現象，但相對ECL，EC格式的耗散更多，抹平得也較為嚴重.

3.2 一維Euler方程組數值試驗

3.2.1 一維Euler方程Lax激波管問題

取計算區域為［－0.5，0.5］，初始條件為

圖3 Burgers方程間斷初值問題的數值結果Fig.3 Numerical results of discontinuous IVP of Burgers equation（t＝0.3，N＝50）

采用齊次Neumann邊界條件，取CFL＝0.3，空間網格數為N＝200，計算到時間t＝0.16．其精確解從左到右依次包括稀疏波、接觸間斷和激波.采用C、EC2、EC2L格式對密度的計算結果進行比較，見圖4（a）－（c）.同樣由于C格式缺乏耗散機制而表現出強烈的振蕩現象，而其他兩種格式的解都沒有產生振蕩，其中由于EC2格式是一階精度，在解的間斷位置抹平得較為嚴重，而EC2L格式在解的激波和接觸間斷處都表現出比較銳利的捕捉效果.圖4（d）中三種格式和精確解的總熵－ΔxΣi（ρiSi）隨時間變化的比較也可體現，C格式的總熵不變，精確解、EC2和EC2L格式的總熵都有所耗散，當然精確解的熵耗散最少也最合適，EC2的多于EC2L，所以比EC2L抹平得嚴重.

3.2.2 一維Euler方程Sod激波管問題

取計算區域為［－0.5，0.5］，初始條件為

采用齊次Neumann邊界條件，取CFL＝0.3，空間網格數為N＝200，計算到時間t＝0.1．精確解包括三部分，分別為稀疏波、接觸間斷和激波.密度的計算結果如圖5（a）－（c）所示，與前述算例一樣，熵守恒C格式的解有明顯的振蕩，EC2格式對解的捕捉有很大的改進，但間斷處的抹平現象同樣比較嚴重，而EC2L格式對此進行了良好的改善.圖5（d）是各格式總熵變化的比較，與前算例的結論一樣：C的總熵不變，精確解的熵耗散最適中，EC2的熵耗散比EC2L的多，且這兩者的熵耗散都比精確解的多.

3.2.3 一維Euler方程低密度流問題

在區域［－0.5，0.5］上取初始條件

圖4 一維Euler方程Lax激波管問題的數值結果Fig.4 Numerical results of 1D Euler equation in Lax shock tube problem（t＝0.16，N＝200）

圖5 一維Euler方程Sod激波管問題的數值結果Fig.5 Numerical results of 1D Euler equation in Sod shock tube problem（t＝0.1，N＝200）

采用齊次Neumann邊界條件，取CFL＝0.3，空間網格數N＝200，計算到時間t＝0.05．其精確解包括兩個強稀疏波和一個平凡的接觸間斷.解的中間部分壓力幾乎為0，接近于真空狀態，很多數值方法在該位置會出現壓力為負的情況，使計算無法進行，如熵守恒格式C計算到第11步時，壓力就為負了，圖6（a）只給出了第11步的數值結果.而其他兩種格式都能使計算順利完成，且都對解的結構有準確的捕捉，見圖6 （b）－（c）.對本算例，同樣有EC2L格式的結果優于EC2格式的結果，進一步體現了高分辨率熵相容格式的優點，圖6（d）是其總熵關于時間變化的對比，可見C格式前10步的總熵依然沒有變化，而EC2L的熵耗散幾乎與精確解的完全重合，但如果放大圖的比例還是能看出，精確解的熵耗散稍小些，EC格式比起這兩者，其耗散就比較多了.

圖6 一維Euler方程低密度流問題的數值結果Fig.6 Numerical results of 1D Euler equation in Low density problem（t＝0.05，N＝200）

3.2.4 一維Euler方程強稀疏波問題

在區域［－10，15］上，初始條件為

采用齊次Neumann邊界條件，取CFL＝0.3，空間網格數為N＝200，計算到時間t＝0.01．其精確解包括稀疏波，接觸間斷和激波.依然采用上述三種格式，密度的計算結果如圖7所示.熵守恒C格式有明顯的振蕩，EC2和EC2L格式在稀疏波段均有一個平滑的過渡，但在間斷處，同3.2.1和3.2.2的算例一樣，EC2格式的抹平較嚴重，而EC2L有很大的改進，充分體現了EC2L格式的優良性.圖7（d）顯示了各格式對應的總熵關于時間的變化，它們各自耗散的多少說明了與前面幾個算例相同的結論：在精確解熵耗散的對比下，熵守恒格式的熵不變；熵相容EC2格式的耗散最多，導致其數值結果在間斷處抹平嚴重；高分辨率熵相容EC2L格式的耗散少且較為適中，表現了銳利的捕捉效果.

3.2.5 一維Euler方程爆炸波問題

在計算區域［－0.5，0.5］上，初始條件為

圖7 一維Euler方程強稀疏波問題的數值結果Fig.7 Numerical results of 1D Euler equation in strong expansion problem（t＝0.01，N＝200）

采用反射邊界條件，取CFL＝0.4，空間網格數為N＝400，計算到時間t＝0.038．由于該問題的初始數據有兩次間斷，已經無法用經典的精確Riemann求解器計算其精確解，本文采用文獻［16］中的WENO-RF－3格式在800個空間網格上的數值解作為精確解，對比EC2和EC2L格式的數值結果，分別見圖8（a）和（b），很顯然EC2L比EC2有很大的改進，這也進一步說明了EC2L格式具有一定的實用性.

圖8 一維Euler方程爆炸波問題的數值結果Fig.8 Numerical results of 1D Euler equation in blast waves problem（t＝0.038，N＝400）

4 結論

通過推廣原始的Superbee限制器得到一類廣泛的GSbee類限制器，隨著α，β，γ的連續變化，這類限制器覆蓋了整個二階TVD區域；取α＝1，β＝2，γ＝1構造出S-M限制器；然后根據傳統的構造二階TVD格式的思想，利用此限制器構造得到高分辨率熵相容格式，同時對Euler方程的熵相容格式的幾個參數做了簡單的調整，使其對接觸間斷的捕捉更精確；在本文的第3節通過幾個數值算例將所得高分辨率熵相容格式的數值結果與熵相容格式和熵守恒格式的數值結果比較分析，從直觀上說明了S-M限制器的高分辨率熵相容格式的合理性和優良性.另外需要說明的是

1）本文只是在熵守恒和熵相容格式的基礎上，采用S-M限制器構造高分辨率格式，至于S-M限制器是否適合于其他經典格式的高分辨率格式構造，還有待進一步研究；

2）由于熵相容EC2格式對不同的問題需要調整一些參數的大小［8］，如αmin、αmax、β1和β2，所以得到的相應的高分辨率熵相容格式如果要對所有的問題（包括定常［8］和非定常問題）都適用，同樣需要調整各個參數的大小.為了進一步提高該格式的實用性，有待下一步從Euler方程熵相容格式本身的構造上進行改進，如修正其格式中的穩定項和熵增項；

3）最后是關于熵相容格式和高分辨率熵相容格式向高維問題的推廣使用，由于高維情形下每個點的信息具有無限多的傳播方向，為了更準確地反映多維效應的影響，我們將在未來的工作中將本文的格式與旋轉Riemann求解器結合，實現對高維問題的計算研究.

［1］Lax P D.Weak solutions of non-linear hyperbolic equations and their numerical computations［J］.Comm Pure Appl Math，1954，7（1）：159－193.

［2］Lax P D.Hyperbolic systems of conservation laws and the mathematical theory of shock waves［C］∥11th of SIAM Regional Conferences Lectures in Applied Mathematics.Philadelphia：Society for Industrial and Applied Mathematics，1973：1－48.

［3］Tadmor E.The numerical viscosity of entropy stable schemes for systems of conservation laws［J］.Math Comp，1987，49 （179）：91－103.

［4］Tadmor E.Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems［J］.Acta Numer，2003，12（1）：451－512.

［5］Tadmor E，Zhong W.Entropy stable approximations of Navier-Stokes equations with no artificial numerical viscosity［J］.J Hyperbolic Differ Equ，2006，3（3）：529－559.

［6］Tadmor E，Zhong W.Energy preserving and stable approximations for the two-dimensional shallow water equations［M］∥Mathematics and computation，a contemporary view，Proc of the third Abel symposium.Berlin Heidelberg：Springer，2008：67 －94.

［7］Roe P L.Entropy conservative schemes for Euler equations［R］.Talk at HYP 2006，Lyon，France.2006.

［8］Ismail F，Roe P L.Affordable，entropy-consistent Euler flux functions II：Entropy production at shocks［J］.J Comput Phys，2009，228（15）：5410－5436.

［9］Toro E F.Riemann solvers and numerical methods for fluid dynamics：A practical introduction［M］.Springer，1997.

［10］Thomas J W.Numerical partial differential equations：Finite difference methods［M］.Springer，1995.

［11］Gottlieb S，Shu C W，Tadmor E.High order time discretizations with strong stability properties［J］.SIAM Review，2001，43 （1）：89－112.

［12］Roe P L.Some contributions to the modeling of discontinuous flows［C］∥Large-scale computations in fluid mechanics. Providence，RI：American Mathematical Society，1985：163－193.

［13］Yee H C.Construction of explicit and implicit symmetric TVD schemes and their applications［J］.J Comput Phys，1987，68 （1）：151－179.

［14］Sweby P K.High resolution schemes using flux limiters for hyperbolic conservation Laws［J］.SIAM J Num Analy，1984，21 （5）：995－1011.

［15］Sweby P K，Baines M J.On convergence of Roe's scheme for the general non-linear scalar wave equation［J］.J Comput Phys，1984，56（1）：135－148.

［16］Jiang G S，Shu C W.Efficient implementation of weighted ENO schemes［J］.J Comput Phys，1996，126（1）：202－228.

High Resolution Entropy Consistent Schemes for Hyperbolic Conservation Laws

REN Jiong1，2，FENG Jianhu1，LIU Youqiong1，LIANG Nan1

（1．College of Science，Chang'an University，Xi'an 710064，China；2．Department of Fluid Dynamics，College of Aeronautics，Northwestern Polytechnical University，Xi'an 710072，China）

To improve accuracy of entropy consistent schemes，we proposed high resolution entropy consistent schemes by inserting a new flux limiter into entropy consistent schemes.It uses limiter mechanism to construct high resolution schemes.In constructing high resolution entropy consistent schemes of Euler equations，we improve resolution of contact discontinuity by adjusting parameters of corresponding entropy consistent schemes.Several numerical experiments illustrate robustness and essentially non-oscillations of the schemes.

hyperbolic conservation laws；entropy consistent scheme；limiter；high resolution

date：2013－07－17；Revised date：2013－12－31

O354；O241.82

A

2013－07－17；

2013－12－31

國家自然科學基金（11171043）和長安大學中央高校基本科研業務費（CHD2102TD015）資助項目

任炯（1985－），女，山西呂梁，碩士生，從事科學與工程中的高性能計算技術研究，E-mail：rensj6962＠163.com

*通訊作者

1001-246X（2014）05-0539-13