二維移動網(wǎng)格矢通量分裂法

李彬彬，鄭素佩，王 令

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

0 引言

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

1 控制方程

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

(1)

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

(2)

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

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

2 網(wǎng)格剖分

(3)

邊界上函數(shù)積分為

(4)

(5)

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

(6)

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

(7)

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

3 基于變分原理的移動網(wǎng)格法

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

3.1 二維移動網(wǎng)格法

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

(8)

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

(9)

令G=ωI，對應(yīng)的歐拉-拉格朗日方程為

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

(10)

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

(11)

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

(12)

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

(13)

3.2 基于迎風(fēng)格式的物理量守恒映射

(14)

(15)

(16)

為了提高間斷處的分辨率，構(gòu)造二階精度的迎風(fēng)型物理量守恒映射，

(17)

(18)

(19)

(20)

3.3 基于非結(jié)構(gòu)網(wǎng)格的矢通量分裂法

(21)

式中：zci為Ik界面的中點(diǎn)坐標(biāo)。旋轉(zhuǎn)不變性的數(shù)值通量Hk為

(22)

4 數(shù)值算例

算例1雙馬赫反射問題

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

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

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

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

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

初值為

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

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

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

初值為

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

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

5 結(jié)語

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