歐拉方程數值求解的高精度通量分裂方法

鄭秋亞，蘇寧亞，梁益華

(1.長安大學 理學院，陜西 西安 710064；2.中國航空計算技術研究所 航空氣動力數值模擬重點實驗室，陜西 西安 710068)

0 引言

自20世紀50年代至今，計算流體力學發展出多種數值方法。優良的數值方法具備計算效率高、穩健性強、間斷處分辨率高等特點(如迎風格式)。迎風格式根據網格上的信息傳播方向對方程進行離散化處理，使方程所模擬的物理現象特征就與離散化方法結合起來。該格式結合物理背景得來，現今仍被廣泛使用與改進。迎風格式又分為矢量通量分裂(FVS)格式[1]和通量差分分裂(FDS)格式[2]。兩種格式各有優缺點：FDS格式有數值耗散小、間斷分辨率和黏性分辨率高的優點，缺點是計算效率較低、穩健性較差，在高速流動情況下容易出現紅斑現象[2-4]；FVS格式簡單易用、計算量小，激波附近不容易發生非物理振蕩，但是由于復雜度降低，導致其在間斷處分辨率低，特別是遇到接觸間斷和剪切波時。

隨后大量學者將兩種格式結合起來，得到一類新的混合迎風格式，其即有FVS格式高計算效率和穩健性的優點，又保留了FDS格式的高分辨率特點。混合迎風格式有兩種基本思想，第1種是對流迎風分裂(AUSM)格式[5-9]。AUSM格式根據物理現象將無黏守恒通量分裂成對流項和壓力項兩部分。AUSM格式中數值耗散項的系數以標量形式而不是以矩陣來計算的，這種形式在提高計算效率的同時又保留了FVS格式的簡單特性。

另一種混合迎風格式是對流迎風和分壓(CUSP)格式[10-11]。其中對流項的總能量為總焓的稱為總焓對流迎風和分壓(H-CUSP)格式[12]。另外一種主流的CUSP格式是總能對流迎風和分壓(E-CUSP)格式[13-17]，其對流項的總能量是總能。

E-CUSP格式是從特征分析出發，將雅可比矩陣中的特征值u±a分裂成對流速度u和聲速±a兩部分，得到相應的對流項和壓力項。在構造界面通量時，是將無黏通量分解為守恒型通量和壓力通量，避免了矩陣運算，提高了計算效率。但是該格式對于間斷處的抹平計算還是不太理想。

對于具有沖擊波或接觸波的問題，用本質無振蕩(ENO)格式或WENO來處理是非常有成效的。其中WENO通過使用所有候選模板的凸組合來代替ENO中最平滑的模板，WENO格式比ENO格式具有更多優勢。例如，模板中使用平滑的數值通量，具有更好的收斂速度，在平滑區域中接近高階最優精度。至今為止，WENO[18-20]在對流擴散方程中特別是在雙曲守恒律方程中，WENO被廣泛應用。

本文的目的是在空間方向上將Zha的E-CUSP格式與WENO相耦合，來捕獲間斷并在平滑區域實現高階精度，在時間方向上采用4階總變差遞減(TVD)Runge-Kutta方法[21]。本文首先推導E-CUSP格式，將其分別與3階WENO和5階WENO耦合(分別用WENO-3-E-CUSP、WENO-5-E-CUSP表示)。其次進行數值實驗，將E-CUSP格式、WENO-3-E-CUSP格式和WENO-5-E-CUSP格式進行對比，通過數值結果分析發現，新的格式更加準確、穩健和有效，從而表明采用WENO建立的高分辨率模型可以更好地模擬激波管問題。新的格式能夠保證在解的光滑區域精度更高，在解的間斷區域保持陡峭的間斷過渡和本質無振蕩性質。最后對全文進行總結和展望。

1 數值方法

1.1 控制方程

控制方程為一維Euler方程，其守恒形式為

(1)

1.2 E-CUSP格式

利用有限體積法，對(1)式在區間單元上進行積分，可得

(2)

E-CUSP格式的主要思想就是將通量分解為守恒型通量和壓力型通量兩項。這樣處理可以有效避免矩陣運算帶來的時間消耗，大大減少計算量，提高計算效率。(1)式中雅克比矩陣可表示為

F=TΛTTU.

(3)

由(3)式可得

其中Fc表示守恒型通量，Fc=u[ρ,ρu,ρe]T;Fp表示壓力型通量，Fp=[0,p,pu]T.

守恒型通量Fc可表示為

(4)

在壓力型通量Fp中，動量方程中的壓強p可表示為

p=P+pl+P-pr,

(5)

(6)

最終E-CUSP格式中的通量F可表示為

(7)

1.3 WENO格式

WENO格式用于評估變量Ul和Ur.其中變量Ul的WENO格式可以寫成：

式中：s為候選模板數；qk(k=0,…,s-1)為不同模板中變量的重建；ωk為權重，

(8)

(9)

ε是為了避免分母變為0，Jiang等[20]的數值實驗表明，只要ε的選擇在10-5～10-7范圍內，對結果的影響不大，通常ε取10-6.

變量Ur的WENO格式與Ul類似，向右移動1個空間節點即可。

對于3階(s=2)WENO，有

對于5階(s=3)WENO，有

1.4 Runge-Kutta方法

2 數值實驗

2.1 算例1：一維Euler方程激波管問題

一維激波管非定常流動問題中包含激波、接觸間斷等流動特征，在數值計算中網格容易生成，初始條件、邊界條件容易處理，常被用來考察數值格式性能。

本文在區域[0,1]上求解初值問題：

(10)

采用Neumann邊界條件，取空間網格點數200個，計算中條件數CFL=0.1.

圖1所示為速度、壓力、馬赫數和密度在t=0.25時用E-CUSP格式、WENO-3-E-CUSP格式和WENO-5-E-CUSP格式分別求解一維Euler方程激波管問題得到的無量綱化數值計算結果。從圖1中可以看出，E-CUSP格式在激波和接觸間斷處的過渡較寬，WENO-3-E-CUSP格式和WENO-5-E-CUSP格式在接觸間斷和激波處的計算結果明顯優于E-CUSP格式。尤其在激波處過渡帶明顯變窄。圖1表明WENO-3-E-CUSP格式和WENO-5-E-CUSP格式的數值模擬輪廓基本相同，在激波和接觸間斷處WENO-5-E-CUSP格式要比WENO-3-E-CUSP格式所獲結果更加接近理論解。

圖1 一維Euler方程激波管問題的數值解Fig.1 Numerical solution of shock tube problem in one-dimensional Euler equation

2.2 算例2：一維Euler方程Lax激波管問題

在區域[0,1]上求解初值問題：

(11)

采用Neumann邊界條件，取網格點數200個，計算中條件數CFL=0.03.圖2所示為速度、壓力、馬赫數和密度在t=0.16時用E-CUSP格式、WENO-3-E-CUSP格式和WENO-5-E-CUSP格式分別計算一維Euler方程Lax激波管問題的無量綱化后的數值模擬結果。

圖2 一維Euler方程Lax激波管問題的數值解Fig.2 Numerical solution of the Lax shock tube problem in one-dimensional Euler equation

該問題的理論解由稀疏波、接觸間斷和激波3個波構成。從圖2的局部放大圖中，可以觀察到WENO-3-E-CUSP格式和WENO-5-E-CUSP格式在接觸間斷和激波處的計算結果依舊明顯優于E-CUSP格式，并且WENO-5-E-CUSP格式比WENO-3-E-CUSP格式所獲結果更加接近理論解，與算例1觀察到的結果相同，表明新的格式對解的捕捉更加精準。比較圖1和圖2可以看出，算例2中稀疏波的波尾部和激波的頂部有過沖現象。算例2的稀疏波與接觸間斷之間過渡帶較寬且出現明顯的振蕩。圖3所示為算例2中壓力的局部放大圖。

圖3 算例2壓力局部放大圖Fig.3 Partial enlarged detail of pressure in Example 2

由于E-CUSP格式的守恒變量經過WENO的重構，使得其計算量大大增加。表1展示了3種格式在算例1和算例2所用的時間，對比可發現新格式的計算效率較低、時間代價大。WENO-3-E-CUSP格式和WENO-5-E-CUSP格式的計算時間基本相同。算例2中，WENO-5-E-CUSP格式所用的計算時間甚至比WENO-3-E-CUSP格式要少。綜合以上分析可知，采用WENO-5與E-CUSP進行耦合是最佳的選擇。

表1 不同數值格式的CPU計算時間Tab.1 CPU computing times of different schemes s

3 結論

本文通過數值模擬一維Euler方程組，以相應的Riemann問題解為基礎，考察了低耗散E-CUSP格式耦合高精度WENO格式后的性能。算例結果表明，由于WENO的重構使得計算量大大增加，新格式的計算效率相比E-CUSP格式的計算效率較低。在較大時間代價的情況下，新格式對接觸間斷和激波的捕捉能力較強，尤其是對激波的捕捉能力僅需要3～4個網格單元。采用WENO-5耦合E-CUSP格式是最佳的選擇，實驗結果表明新的格式有更高的準確性和穩健性。

本文主要針對歐拉方程進行數值模擬，來驗證新格式的優勢。對計算流體力學中的實際問題，如磁流體動力學[22]、機翼展開[23]和空泡流動[24]等問題，應用新格式進行數值模擬將會產生更好的預期結果。將精度更高的WENO[25-26]與E-CUSP格式耦合，并結合自適應算法來提高計算效率，是進一步研究的方向。