一種改進的HLLEM格式及其激波穩定性分析*

胡立軍 袁海專 杜玉龍

1) (衡陽師范學院數學與統計學院,衡陽 421002)

2) (湘潭大學數學與計算科學學院,湘潭 411105)

3) (北京航空航天大學數學科學學院,北京 100191)

使用低耗散激波捕捉格式對高超聲速流動問題進行數值模擬時經常會遭受不同形式的激波不穩定性.本文基于二維無黏可壓縮Euler方程,對低耗散HLLEM格式進行激波穩定性分析.結果表明: 激波面橫向通量中切向速度的擾動增長誘發了格式的不穩定性.通過增加耗散來治愈HLLEM格式的激波不穩定性.為了避免引入過多的耗散進而影響剪切層的分辨率,定義激波探測函數和亞聲速區探測函數,使得只有在計算激波層亞聲速區的橫向數值通量時才增加耗散,其余地方的數值通量依然采用低耗散的HLLEM格式來計算.穩定性分析和數值模擬的結果表明,改進的HLLEM格式不僅保留了原格式高分辨率的優點,還大大提高了格式的魯棒性,在計算強激波問題時能夠有效地抑制不穩定現象的發生.

1 引 言

盡管流體力學數值方法的研究在最近幾十年取得了巨大的進步,然而依然存在一些挑戰性的問題,涉及激波、剪切層、激波/湍流邊界層干擾等現象的高超聲速流動問題的數值模擬就是其中之一[1,2].構造精確、高效、魯棒性好的激波捕捉格式是精確模擬高速流動問題的關鍵.近似黎曼解法器是一類流行的激波捕捉方法.一個好的黎曼解法器不僅要能夠精確地分辨各種間斷,還要具有很好的魯棒性,能夠抑制激波不穩定現象的發生.然而,能夠精確分辨接觸間斷和剪切波的低耗散黎曼解法器(例如Roe,HLLC和HLLEM)在計算強激波問題時都會產生嚴重的不穩定現象,從而影響數值模擬的可靠性.因此,分析低耗散數值格式激波不穩定性產生的根源并且在保留高分辨率優點的同時去治愈格式的激波不穩定性成為了計算流體力學數值方法研究的一個重要問題.

利用線性擾動分析方法,Quirk[3]指出壓力和密度場的相互作用誘發了某些格式的不穩定性,并且在強激波附近使用耗散的HLLE格式來消除Roe格式在數值模擬中的不穩定現象.Gressier和Moschetta[4]利用特殊的擾動分析得出結論: 能夠精確捕捉接觸間斷的數值格式通常都會遭受激波不穩定性的困擾.Dumbser等[5]設計了一種矩陣穩定性分析方法來分析格式激波不穩定性產生的根源并且發現在計算二維穩態激波時數值格式的不穩定性與數值激波的結構密切相關.謝文佳等[6]使用奇偶失聯線性方法對某些迎風格式進行穩定性分析,結果表明: 具有充足剪切黏性的數值格式可以避免一類激波不穩定現象(“紅玉”現象)的發生,通過分析,Shen 等[7]也得到了同樣的結論.通過增加數值耗散來治愈格式的不穩定性已經成為了一種主流方法.一般來說,有以下三種方式來增加數值耗散: 將耗散格式和低耗散格式進行混合,增加數值格式的人工黏性,控制數值格式的耗散機制.Kim等[8]構造了一種HLLC-HLLE混合格式,通過定義開關函數將強激波附近數值通量的計算由不穩定的HLLC格式切換成穩定的HLLE 格 式.Wu 等[9]和 Shen 等[10]構 造 了一 種Roe-HLLE混合格式來消除Roe格式的不穩定性.與全面混合不同的是,他們僅僅在質量方程和切向動量方程將這兩種格式進行混合.Ren[11]提出了一種旋轉Roe黎曼解法器,通過自動地引入人工黏性來抑制原始Roe格式不穩定現象的發生.Rodionov[12]用人工黏性系數代替分子黏性系數來治愈Godunov-型數值格式的“紅玉”現象并且在文獻[13,14]中分別討論了高階Godunov-型格式和三維情形下的激波不穩定性.Chen等[15]發現Roe格式和HLLE格式在亞聲速區剪切速度的差異可以解釋這兩種數值格式不同的激波行為,并且在激波附近亞聲速區增加Roe格式的剪切黏性來抑制不穩定現象的發生.Simon和Mandal[16,17]使用反擴散控制和選波修正兩種不同的策略來控制HLLC格式的耗散機制進而抑制格式的不穩定性,并且使用同樣的策略來構造激波穩定的HLLEM格式.Xie等[18]利用相鄰網格的壓力比定義開關函數來控制HLLEM格式的耗散水平,進而治愈其不穩定性.與增加數值耗散來抑制激波不穩定性的流行方法不同的是,Fleischmann等[19]通過修改Roe格式特征值的計算進而得到一種耗散更小、魯棒性更好的Roe-M格式.

盡管眾多研究人員對于數值格式的激波不穩定性進行了大量的研究,但是關于激波不穩定性產生的根源仍然存在爭議,因此很多治愈方法都存在一定的局限性.本文以低耗散HLLEM格式為例,采用線性穩定性分析方法來推導兩種特殊擾動情形下的擾動發展方程,并且結合相應的數值實驗來探究激波不穩定性產生的根源.根據穩定性分析的結論提出一種改進的HLLEM格式來治愈原格式的不穩定性并且進行數值實驗來驗證改進格式的魯棒性和精度.

2 控制方程組和HLLEM黎曼解法器

2.1 二維歐拉方程組

考慮二維理想氣體流動的守恒律方程組

式中,守恒量U,矢通量F(U) 和G(U) 的定義分別為:

其中ρ為密度,u和v分別為全局坐標系中x方向和y方向的流體速度,p為靜壓力,E為總能.使用理想氣體狀態方程來使守恒律系統封閉:

其中比熱比γ=1.4.求解方程組(1)的守恒型數值方法可以表示成

其中 ?t為時間步長,?x和 ?y為空間步長,Fi+1/2,j和Gi,j+1/2為x方向和y方向的數值通量.

2.2 HLLEM黎曼解法器

Einfeldt等[20]在HLLE通量的基礎上添加與線性退化場相對應的反擴散項,從而得到了一種可以精確分辨接觸間斷和剪切波的HLLEM格式,其通量表達式為

其中,SL和SR表示左、右波速,R2,3表示通量雅克比矩陣與線性退化場相對應的右特征向量,α2,3為對應的波強,δ2,3表示控制耗散大小的系數,他們的計算公式分別為:

與其他可以捕捉接觸間斷和剪切波的低耗散格式一樣,使用HLLEM格式計算多維強激波問題時也會產生不同形式的激波不穩定現象,這大大限制了HLLEM格式在高超聲速流動問題數值模擬中的應用.接下來首先探究HLLEM格式激波不穩定性產生的根源,進而采用針對性的方法來治愈其不穩定性.

3 HLLEM 格式激波不穩定性的機理分析

3.1 線性穩定性分析

流體介質以速度u0=M0a0沿著x正方向(縱向) 運動,沿y方向 (橫向) 的速度v0=0 ,其中M0為馬赫數,為當地聲速.為了弄清楚是哪個方向的擾動增長造成了HLLEM格式的不穩定性,考慮在初始分布上增加兩類特殊的擾動.

首先,在x方向上增加奇偶對稱擾動,擾動后的流場狀態為:

在超聲速(M0>1 )的情況下,擾動發展方程為:

其中σx=?t/?x.(11)式的特征值為:

當σx<1/(u0+a0) 時 ,此時所有擾動量都能有效衰減.選取兩組不同的初始擾動量和初值:

經過時間迭代20步,每個擾動量的發展趨勢如圖1所示,其中σx取 0.2.

從圖1 可以看到,當x方向存在擾動時,在超聲速情形下,所有擾動量都會衰減到0,此時數值格式是穩定的.

圖1 超聲速情形下擾動量的發展曲線 (a)初始擾動量和初值(I);(b)初始擾動量和初值(II)Fig.1.Evolutionary curves of perturbation under supersonic condition: (a) Initial perturbation values and initial conditions (I);(b) initial perturbation values and initial conditions (II).

在亞聲速( 0

(14)式的特征值為:

當σx<1/(a0+u0) 時 ,1(i=1,2,3,4) ,此時所有擾動量均能有效衰減.同樣選取兩組不同的初始擾動量和初值:

經過時間迭代20步,每個擾動量的發展趨勢如圖2所示,其中σx取 0.4.

從圖2 可以看到,當x方向存在擾動時,在亞聲速情形下,所有擾動量也會衰減到0.因此當流體沿著x方向運動時,該方向的擾動不會誘發HLLEM格式的不穩定性.

接下來,考慮在y方向增加奇偶對稱擾動,x方向上流場無擾動.此時的流場狀態為

由于在x方向沒有擾動,因此數值通量Fi?1/2,j=Fi+1/2,j=0.使用 HLLEM格式計算數值通量Gi,j?1/2和Gi,j+1/2,進而得到擾動量的發展方程為:

(18)式的特征值為:

圖2 亞聲速情形下擾動量的發展曲線 (a)初始擾動量和初值(III);(b)初始擾動量和初值(IV)Fig.2.Evolutionary curves of perturbation under subsonic condition: (a) Initial perturbation values and initial conditions (III);(b) initial perturbation values and initial conditions (IV).

其中σy=?t/?y.當σy<1/a0時,和是衰減的,而密度擾動和切向速度擾動不會衰減.設初值ρ0=1.4,u0=1 ,p0=1 ,選取以下四組不同的初始擾動量來觀察在該擾動下各個物理量的擾動發展趨勢

3.2 激波不穩定性的多維特性

關于激波不穩定性是一種一維激波的異常現象還是一種多維現象仍然是有爭議的.Chauvat等[21]和Kitamura等[22]發現激波不穩定性與內部激波結構有關并且主張它是一種一維激波的異?，F象.然而,純粹的一維算例不會遭受激波不穩定現象的困擾,并且3.1節的線性穩定性分析的結論表明:y方向的擾動會造成x方向的激波不穩定性,反之亦然.因此,本文認為它是一種多維現象.

圖3 橫向擾動下擾動量的發展曲線 (a)初始擾動量 (V);(b)初始擾動量 (VI);(c)初始擾動量 (VII);(d) 初始擾動量 (VIII)Fig.3.Evolutionary curves of perturbation in transverse direction: (a) Initial perturbation values (V);(b) initial perturbation values (VI);(c) initial perturbation values (VII);(d) initial perturbation values (VIII).

圖4 二維 Sedov 爆轟波問題的密度等值線 (a) HLLEM;(b) x-HLLE+y-HLLEM (c) x-HLLEM+y-HLLEFig.4.Density contours for 2D Sedov blast wave problem: (a) HLLEM;(b) x-HLLE+y-HLLEM (c) x-HLLEM+y-HLLE.

接下來,計算文獻[19]中的二維Sedov爆轟波問題來證明激波不穩定性的多維特性.初始時,區域 [ 0,2.4]×[0,2.4]中心的壓力pin=3.5×105,其余地方的壓力pout=10?10,整個區域的密度為1,速度u0和v0均為 0.計算中使用 4 80×480 的笛卡爾網格,上、下、左、右都使用反射邊界條件.圖4為計算得到的二維Sedov爆轟波問題的密度等值線,可以看到,采用 HLLEM 格式進行計算時,在激波面和坐標軸平行的位置出現了明顯的“紅玉”現象.當x方向的通量采用耗散的HLLE格式計算時,y方向的“紅玉”現象消失了,反之亦然.該數值實驗的結果與線性穩定性分析的結論相一致,從而進一步證明了激波不穩定現象的多維特性.

4 一種改進的 HLLEM 格式

激波面橫向擾動的發展方程(18)式表明: 密度和切向速度的擾動不會衰減.由于激波不穩定性是一種多維現象,并且一維問題存在熵波但是沒有剪切波.因此我們認為: 激波面橫向通量中切向速度的擾動增長造成了HLLEM格式的激波不穩定性.通過增加橫向通量的耗散來抑制不穩定現象的發生,其具體表達式為

在耗散項的作用下,橫向擾動的發展方程式為:

(22)式表明橫向通量中的切向速度的擾動也會衰減.

4.1 激波探測函數

在實際的數值模擬中,只需要在激波附近增加耗散就可以抑制不穩定現象的發生,為了避免在不適當的位置增加耗散進而影響剪切層的分辨率,利用網格界面的壓力比來探測激波的位置:

顯然,h∈(0,1) ,且在激波附近壓力差較大,從而h→0.由前面的分析可知: 在x方向增加耗散可以抑制y方向的不穩定性,反之亦然.因此,計算x方向的數值通量Fi+1/2,j時,只需檢測與之相鄰的y方向的網格界面.計算y方向的數值通量Gi,j+1/2時,只需檢測與之相鄰的x方向的網格界面.即

由于不穩定性的嚴重程度會隨著激波強度的增加而增加,我們采用文獻[15]中的余弦函數來定義激波探測函數:

此外還給出另外兩種激波探測函數:

圖5畫出了這三種激波探測函數的曲線圖.從中可以看到三個函數的取值都會隨著激波強度的減弱而減小,且對于固定的h,滿足g1>g>g2.

圖5 三種激波探測函數的函數曲線Fig.5.Curves for three different shock-detecting functions.

4.2 M-HLLEM格式

Xu和Li[23]在分析定常激波問題時發現激波不穩定性是由激波層亞聲速區的擾動增長導致的,這一結論得到了眾多研究人員的支持[6,15,19].因此定義亞聲速區探測函數

式中M表示馬赫數.在亞聲速區 0<φ<1 ,在超聲速區φ=0.

因此,在一般的結構化四邊形網格下,最終的耗散項可以表示為

其中n=(nx,ny) 為網格界面單位外法向量,g和φ分別為激波探測函數和亞聲速區探測函數.因此,一種改進的HLLEM格式(M-HLLEM,Modified HLLEM)數值通量的表達式為

5 數值結果和分析

本節計算一些典型的算例來檢驗本文構造的M-HLLEM格式的魯棒性和精度.

5.1 隨機數值噪聲問題

馬赫數為6的平面激波沿著x方向運動,初始時激波位于x=5 ,右側波前的初始條件為(ρ0,u0,v0,p0)=(1.4,0,0,1),左側的波后狀態可以通過激波和馬赫數之間的關系式來得到.計算區域 [ 0,1000]×[0,20]被劃分成 1 000×20 的笛卡爾網格.這是一個純粹的一維激波,為了誘發不穩定性,在波前的初始分布上增加取值從?0.5×10?5到0.5×10?5的隨機數值擾動

其中αk(k=1,2,3,4) 為?0.5 到 0.5 之間的隨機數.

在該算例中,y方向速度的最大量值max(|v|)可以用來衡量不穩定性的嚴重程度.圖6展示了t=120時,原始的HLLEM格式以及與不同的激波探測函數結合的改進格式的密度等值線圖.圖7展示了y方向速度的最大量值隨著時間變化的曲線圖.可以看到,HLLEM格式的激波結構完全被破壞,且y方向速度的最大量值也從 1 0?5量級增長到100量級.而三種改進格式都得到了清晰的激波面,且y方向速度的最大量值也始終保持在初始時的10?5量級附近.三種激波探測函數的數值結果之間幾乎沒有差異,因此后面的數值算例均選用流行的余弦函數(25)式來計算其函數值.

圖6 隨機數值噪聲問題的密度等值線 (a) HLLEM;(b)M-HLLEM-g;(c) M-HLLEM-g1;(d) M-HLLEM-g2Fig.6.Density contours of random numerical noise problem: (a) HLLEM;(b) M-HLLEM-g;(c) M-HLLEM-g1;(d)M-HLLEM-g2.

圖7 隨機數值噪聲問題 y 方向速度的最大量值Fig.7.Maximum magnitude of velocity in y-direction of random numerical noise problem.

5.2 高超聲速繞柱流問題

高超聲速繞柱流問題經常用來檢驗數值格式是否會遭受一種嚴重的不穩定現象—“紅玉”現象(carbuncle).馬赫數為20的自由來流流經半徑為1的圓柱體,整個區域流場的初始條件(ρ0,u0,v0,p0)=(1.4,20,0,1),具體的計算區域、網格劃分和邊界條件可以參考文獻[24].本文采用320×40的結構化貼體四邊形網格來計算.圖8畫出了時間t=4 時,使用HLLEM和M-HLLEM格式計算得到的密度等值線.可以清晰地看到,HLLEM格式出現了嚴重的“紅玉”現象,而M-HLLEM格式消除了不穩定現象,得到了具有清晰激波面的穩態弓形激波.

圖8 高超聲速繞柱流問題的密度等值線 (a) HLLEM;(b) M-HLLEMFig.8.Density contours of hypersonic flow over a cylinder:(a) HLLEM;(b) M-HLLEM.

5.3 雙馬赫反射問題

Woodward和Colella[25]提出的雙馬赫反射問題是檢驗數值格式魯棒性的一個著名算例.馬赫數為10的斜激波向底部壁面運動,激波面與壁面所成 角 度 為 60°,計 算 區 域 [ 0,4]×[0,1]被 劃 分 成480×120的笛卡爾網格.詳細的初始條件、邊界條件可參考文獻 [25].圖9畫出了時間t=0.2 時,HLLEM和M-HLLEM格式的計算結果.可以看到,HLLEM格式出現了彎曲的馬赫桿和一個三角點,這是一種明顯的非物理現象.而M-HLLEM格式消除了不穩定現象,得到了清晰的馬赫桿.

圖9 雙馬赫反射問題的密度等值線 (a) HLLEM;(b) M-HLLEMFig.9.Density contours of double Mach reflection problem:(a) HLLEM;(b) M-HLLEM.

5.4 二維Sedov爆轟波問題

計算3.2節描述的二維Sedov爆轟波問題,由于該算例涉及高壓力比的強激波問題,因此可以用來作為檢驗數值格式性能的一個挑戰性測試.圖10畫出了時間t=0.1 時,HLLEM和 M-HLLEM格式的壓力等值線.可以清晰地看到,HLLEM格式在激波面和網格線平行的位置出現了四個明顯的“紅玉”現象.而M-HLLEM格式有效地抑制了不穩定現象的發生.

5.5 二維無黏接觸間斷問題

數值格式可以用來計算黏性流的前提條件是能夠精確捕捉接觸間斷.該算例用來檢驗改進的HLLEM格式捕捉接觸間斷的能力.區域[0,1]×[0,1]被均勻劃分成 1 0×10 的正方形網格,上、下兩部分不同速度的兩種流體的初始條件為:

圖10 二維 Sedov 爆轟波問題的壓力等值線 (a) HLLEM;(b) M-HLLEMFig.10.Pressure contours of 2D Sedov blast wave problem:(a) HLLEM;(b) M-HLLEM.

圖11展示了計算迭代1000步后,x=0.5 處的密度分布.可以看到,HLLE格式得到了具有較大耗散的解,而HLLEM和M-HLLEM格式都能準確捕捉該接觸間斷.這說明M-HLLEM格式保留了原HLLEM格式精確分辨接觸間斷的優點.

圖11 二維無黏接觸間斷問題的密度分布Fig.11.Density distribution of 2D inviscid contact discontinuity problem.

6 結 論

本文構造了一種激波穩定的HLLEM格式.線性穩定性分析和相關的數值實驗表明: 橫向通量中切向速度的擾動增長誘發了HLLEM格式的激波不穩定性.通過增加耗散來治愈格式的不穩定性,為了防止過多的耗散影響剪切層的分辨率,定義激波探測函數和亞聲速區探測函數,使得耗散項僅僅添加在激波層亞聲速區域的橫向通量上,其余的數值通量依然采用低耗散的HLLEM格式來計算,從而在不犧牲精度的前提下,提高格式的魯棒性.數值模擬的結果表明: 改進的HLLEM格式消除了激波不穩定現象并且保留了原始HLLEM格式精確分辨接觸間斷的優點,因此該格式可以廣泛應用于高超聲速可壓縮流的數值模擬中.采用類似的方法來改進其他的低耗散格式,結合高階精度重構方法并且將其推廣到三維情形來計算不同的守恒律系統可以作為未來的研究工作.

感謝中國科學院計算數學研究所袁禮研究員的討論.