李彬彬，鄭素佩，王 令

(長安大學 理學院 陜西 西安 710064)

0 引言

在對非線性雙曲守恒律方程進行數值求解時，即使在初始條件充分光滑的條件下，在某一時刻解也可能存在激波或稀疏波等間斷，尋找滿足物理意義的唯一弱解成為研究雙曲守恒律方程的關鍵。文獻[1-3]對矢通量分裂法進行了研究，并給出了通量如何選擇的建議。數值結果的好壞不僅與格式有關，還與網格的分布相關。根據解的性質對網格節(jié)點進行合理的分布，對高效、精確地計算非線性雙曲守恒律方程具有極其重要的作用，為此文獻[4]首次提出了自適應技術。目前，引起學者廣泛關注的是拉格朗日法細化后的移動網格法。移動網格法將方程的求解和網格的移動作為兩個完全獨立的部分，該算法具有易于推廣的特點。自20世紀80年代以來，文獻[5]開始利用自適應移動網格法提高分辨率。已有的大多數研究都是采用移動網格法來求解一維的動力學方程[6]。2003年，文獻[7]提出了求解多維雙曲守恒律方程的移動網格算法，能有效地解決激波間斷問題。到目前為止，將移動網格法與矢通量格式進行耦合的研究還很少。因此，本文將移動網格與矢通量分裂格式進行耦合來求解二維歐拉方程組，通過構造新的監(jiān)控函數以及物理量守恒映射提高移動網格的通用性，既保證求解過程中網格合理分布和提高間斷處的分辨率，也保證原格式的高精度特點。此外，二維歐拉方程組的數值算例模擬結果表明了新算法具有良好的間斷捕捉能力和高分辨率。

1 控制方程

考慮如下二維守恒型歐拉方程[8]：

(1)

記H=[F,G]·n，則式(1)可簡寫為

(2)

時間方向上采用文獻[9]提出的TVD型龍格-庫塔方法進行離散，

式中：L是空間離散算子。

2 網格剖分

(3)

邊界上函數積分為

(4)

(5)

將式(4)和式(5)代入式(3)，可得

(6)

由中值定理，可得式(6)的近似值為

(7)

式中：Δlk是Ik(k=1,…,4)的長度；Hk是Ik(k=1,…,4)中點處的函數值。

3 基于變分原理的移動網格法

對文獻[7]中的監(jiān)測函數和物理量守恒映射方法進行改進，以提高移動網格法的通用性。

3.1 二維移動網格法

z=z(ζ),ζ∈Ωc。

(8)

在變分法中，網格映射使得在計算區(qū)域如下函數達到最小值，具體函數表示式為

(9)

令G=ωI，對應的歐拉-拉格朗日方程為

(ωxξ)ξ+(ωxη)η=0,ζ∈Ωc,

(10)

用中心差分格式離散式(10)，可得

(11)

采用高斯賽德爾迭代法求方程(11)的近似解，可得

(12)

為避免數值解產生較大的誤差，可用低通濾波器對監(jiān)控函數進行光滑化[11]，

(13)

3.2 基于迎風格式的物理量守恒映射

(14)

(15)

(16)

為了提高間斷處的分辨率，構造二階精度的迎風型物理量守恒映射，

(17)

(18)

(19)

(20)

3.3 基于非結構網格的矢通量分裂法

(21)

式中：zci為Ik界面的中點坐標。旋轉不變性的數值通量Hk為

(22)

4 數值算例

算例1雙馬赫反射問題

計算區(qū)域為[0,4]×[0,1]，在其底部有一個始于x=1/6的墻。在計算開始時刻，一個馬赫數為10的右激波通過點(1/6,0)，并與x軸成60°夾角。則初始解為

其中左狀態(tài)UL、右狀態(tài)UR和激波的位置h分別為

左邊界和右邊界分別設為入流邊界和出流邊界，需要根據激波運動對上邊界進行分類討論，在x=1/6時對下邊界采用無穿透絕熱條件。計算網格為160×40，計算時間推進到t=0.2。在數值模擬中，畫圖區(qū)域僅截取[0,3.1]×[0,1.0]。圖1為雙馬赫反射問題的移動網格演化圖和密度圖?？梢钥闯觯苿泳W格集中分布在點(2.7,0.41)和(2.6,0)附近；網格演化圖的間斷區(qū)域分布較多網格節(jié)點，過渡帶明顯變窄，具有較高的分辨率。

圖1 雙馬赫反射問題的移動網格演化圖和密度圖Figure 1 Moving grid evolution map and density map of double Mach reflection problem

算例2二維歐拉方程的黎曼問題1

初值為

其解包含兩個雙馬赫反射，計算區(qū)域為[0,1]×[0,1]，計算網格為100×100，兩側的邊界條件使用出流邊界條件，計算時間為t=0.25。圖2為黎曼問題1的移動網格演化圖和密度圖?？梢钥闯?，激波兩交點(0.92,0.3)和(0.3,0.92)附近網格進行了加密處理；隨著網格節(jié)點自適應加密，弧形激波過渡區(qū)明顯變窄。

圖2 黎曼問題1的移動網格演化圖和密度圖Figure 2 Moving grid evolution map and density map of Riemann problem 1

算例3二維歐拉方程的黎曼問題2

初值為

計算區(qū)域為[0,1]×[0,1]，計算網格為100×100，兩側的邊界條件使用出流邊界條件，計算時間為t=0.25。圖3為黎曼問題2的移動網格演化圖和密度圖?？梢钥闯?，弧形激波和馬赫反射等間斷區(qū)域分布較多的網格節(jié)點；移動網格間斷處進行了自適應加密，間斷處過渡帶相對于固定網格更窄，能更準確地捕捉間斷。

圖3 黎曼問題2的移動網格演化圖和密度圖Figure 3 Moving grid evolution map and density map of Riemann problem 2

5 結語

本文將移動網格算法與高分辨率矢通量分裂格式耦合，得到了一種求解二維雙曲守恒律方程的移動網格矢通量分裂格式法。三個數值算例的模擬結果表明，新算法能有效地捕捉間斷和提高間斷處的分辨率，且易于編程和向高維問題推廣。