高階WENO格式數值粘性對模擬R-T不穩定性的影響

王革，謝昌坦，張斌

(1.哈爾濱工程大學航天與建筑工程學院，黑龍江哈爾濱150001;2.北京航空航天大學宇航學院，北京100083)

受擾動的輕重流體的交界面，當受到方向由重流體指向輕流體的重力或慣性力作用時，擾動將發展，界面將失穩，2種物質將發生湍流混合，這種不穩定現象稱為Rayleigh-Taylor(R-T)不穩定性.當相互接觸的2層流體間存在切向速度差時，界面也會擾動發展，這種不穩定現象稱為 Kelvin-Helmholtz(K-H)不穩定性.一般在R-T不穩定性發展的后期，流體界面處會產生強剪切流，K-H不穩定性被激發并導致流體的小尺度混合.因此，R-T不穩定性常常伴隨著K-H不穩定性，而后者的出現加重了前者后期的非線性發展，加劇了界面附近流體的混合程度.上述2種不穩定性在天體物理、激光聚變、高速碰撞和加速等高技術領域中都是客觀存在的，甚至對水力機械以及各種發動機也都是非常重要的［1］.

由于R-T不穩定性自身的復雜性，一般只有采用高階精度格式才能獲得較好的模擬效果.1987年，Harten等提出了一種完全考慮TVD性質，通過節點模板的單調選擇、擴展來達到高分辨率的數值方法——ENO 格式［2］，而后 Osher、Shu 等避開節點模板選擇，利用加權思想構造了WENO格式［3-4］.大量的理論研究和實驗表明，ENO和WENO格式特別適合用于求解既包含激波、間斷又包含光滑區域的復雜流場.

任何數值格式都具有數值粘性，而且數值粘性的大小有時會對數值解產生很大影響.格式的精度越高，計算網格越密，則數值粘性越小.然而數值粘性的定量求解是一項極其繁雜的工作，這點對于高分辨率格式尤為突出.為此本文主要通過一個經典的二維R-T不穩定性算例定性研究了高階WENO格式的數值粘性，分別應用3、5、7、9階精度WENO格式求解二維無量綱Euler和Navier-Stokes方程獲得了不同網格尺度下無粘和有粘時的數值解.通過研究R-T不穩定性后期激發出的K-H不穩定性的強弱，分析了數值粘性和物理粘性對數值流場細微結構的影響趨勢.

1 控制方程及其數值解法

含有重力源項(假設重力方向豎直向上，重力加速度取)的二維無量綱Navier-Stokes方程為

式中:ρ、p分別表示流體的密度和壓力，u、v分別表示流體沿 x、y方向的速度分量，總能)，γ為比熱比，Pr為普朗特數，Re為雷諾數.當S(U)=0即Re→∞時，式(1)便退化為含有重力源項的Eluer方程.

采用文獻［5］介紹的時間分裂法求解式(1)，即將式(1)分解為Euler方程和源項2部分單獨求解，如下:

式(2)的空間項分別采用基于特征分解法的3、5、7、9 階有限差分型 WENO-LF 格式求解［6-7］，時間項采用三階TVD-Runge-Kutta法求解［3］.為保證與式(2)的精度一致，式(3)的空間項分別采用2、4、6、8 階中心差分格式，時間項采用三階Runge-Kutta法求解.

2 數值算例及結果分析

本算例的計算區域為［0，0.25］×［0，1.0］，計算網格為方形網格，其寬度用h表示.初始時刻界面位于y=0.5處，重氣體位于界面以下，參數為ρ=2.0，u=0.0，p=2y+1，v=-0.025c cos(8πx);輕氣體位于界面以上，參數為 ρ=1.0，u=0.01，p=y+1.5，v=-0.025c cos(8πx)，c 為音速，c=，γ =1.666 7，Pr=0.7.左右邊界設為反射邊界條件，上下邊界設為無反射邊界條件，上邊界參數為 ρ=1.0，u=0，p=2.5，v=0，下邊界參數為 ρ=2.0，u=0，p=1，v=0.計算至時刻 t=1.95.

2.1 無粘情況

無粘情況下的數值流場是通過求解Euler方程得到的.圖1為不同精度WENO格式和不同網格數量條件下區域［0，0.25］×［0.2，0.9］的 15 條密度等值線分布圖.為簡便，引入符號WENO-N-X，其中，N代表WENO格式的精度，N=3、5、7、9;X 代表計算網格的寬度 h，X=C(h=1/400)、M(h=1/800)、F(h=1/1 600).

從圖1可以看出，采用Euler方程模擬R-T不穩定性，數值格式的精度和網格的數量對其數值解有很大影響.WENO-3-C、WENO-3-m和 WENO-5-C的數值結果基本相同.在重力作用下，下層重流體以尖釘的形式侵入上層輕流體，而上層輕流體以氣泡的形式侵入下層重流體，兩層流體交界面形成了清晰的蘑菇狀結構，這種結構也是R-T不穩定性發展到后期的典型結構.由于此時數值格式的數值粘性較大，K-H不穩定性較弱，流體的交界面比較平滑.比較WENO-3-F、WENO-5-M，WENO-5-F，WENO-7-C，WENO-7-M，WENO-9-C，WENO-9-M可知，隨著格式精度的提高和網格數量的增加，數值格式的數值粘性減弱，K-H不穩定性增強，蘑菇的傘狀結構在輕流體的作用下翻轉，并在K-H不穩定性的作用下斷裂、碎化，輕重流體混合加劇，導致傘狀結構內部的流場極其復雜，兩層流體的交界面變得模糊.蘑菇的桿狀結構處也出現了渦狀結構，但因渦的形狀比較大、數量少且排列有序，流體交界面仍然清晰可辨.高精度密網格情況下，數值流場的對稱性結構被破壞，如WENO-7-F和WENO-9-F所示，這完全是由Euler方程本身造成的，此時所得的數值解是非物理的［8-9］.

圖1 無粘時的密度等值線分布Fig.1 The contour lines for density for non-physical viscosity

所以，當采用Euler方程模擬R-T不穩定性時，由K-H不穩定性導致的物質界面處復雜結構的細節強烈依賴于所采用數值格式的數值粘性，網格細化所得的數值解一般不具有強收斂性，而且當網格過度細密時，數值流場的復雜結構往往是非物理的.正確計算R-T不穩定性的方法應當采用帶有真實物理粘性的 Navier-Stokes方程［8，10］，如3.2 節所示，當數值格式的數值粘性遠遠小于流體的物理粘性時，所得的數值解一般是可信的.

2.2 有粘情況

圖2、3分別為雷諾數Re=30 000，60 000時的15條密度等值線分布圖.為了能更好的了解R-T不穩定性的數值解隨網格數量的變化情況，取圖2和圖3中y=0.6處即流場結構最復雜地方的密度沿x方向的分布為研究對象，如圖4、5所示，相應的L1誤差［10］列于表1、2.

圖2 Re=30 000時密度等值線分布Fig.2 The contour lines for density with Re=30 000

可以看出，采用Navier-Stokes方程模擬R-T不穩定性，隨著網格數量的增加，WENO格式的數值粘性逐步減小，流體的物理粘性逐步占優，L1誤差減小.若規定L1≤1.2×10-3時數值解收斂，那么當Re=30 000時，WENO-5-M、WENO-7-M、WENO-9-m及 WENO-N-F的數值流場已基本穩定;Re=60 000時，WENO-7-M、WENO-9-m和 WENO-5-F、WENO-7-F、WENO-9-F的數值解已基本收斂，這說明此時數值格式的數值粘性已經遠遠小于流體的物理粘性了.同一雷諾數條件下，不同精度WENO格式所獲得的穩定流場的結構基本相同，且格式的精度越高，獲得穩定解所需的網格越少.同一精度條件下，雷諾數越高，流體的物理粘性越小，則掩蓋數值粘性的影響所需的網格也就越多，穩定流場的結構越復雜.

圖3 Re=60 000時密度等值線分布Fig.3 The contour lines for density with Re=60 000

圖4 Re=30 000時y=0.6處密度沿x方向的分布Fig.4 Distribution of density along x at y=0.6 with Re=30 000

表1 Re=30 000時WENO-N-X的L1誤差(×102)Table 1 The L1 error of WENO-N-X when Re=30 000

圖5 Re=60 000時y=0.6處的密度沿x方向的分布Fig.5 Distribution of density along x at y=0.6 when R e=60 000

表2 R e=60 000時WENO-N-X的L1誤差(×102)Table 2 The L1 error of WENO-N-X when Re=60 000

圖6給出了Re=140 000時WENO-9的密度在不同網格數量條件下的分布圖，y=0.6處的密度沿x方向的分布如圖7所示.可以看出，圖6中的前3幅圖的密度分布隨網格的加密變化較大，如圖7(a)所示.這是因為在雷諾數很高的情況下，流體物理粘性很小，網格較粗時的數值粘性占絕對優勢，初步加密網格后其值迅速減小但仍遠大于物理粘性.隨著網格的繼續加密，物理粘性的影響開始體現出來并且逐步占優，流場結構隨著網格數量的加密變化緩慢，當網格加密到一定程度后，數值粘性便遠小于物理粘性了，流場結構達到穩定狀態，如圖7(b)所示.穩定流場的流體交界面比較平滑，蘑菇的傘狀結構內部形成了2個較大且清晰的渦狀結構，桿狀結構的下部也形成了2個較小的渦狀結構.

圖6 Re=140 000時WENO-9的密度等值線分布Fig.6 Distribution of numerical density for WENO-9 with Re=140 000

圖7 Re=140 000時y=0.6處的密度沿x方向的分布Fig.7 Distribution of density along x at y=0.6 with Re=140 000

3 結論

本文通過利用高階WENO格式求解無粘Euler和Navier-Stokes方程模擬了二維可壓縮流體中的R-T不穩定性問題.數值結果表明，R-T不穩定性的復雜流場結構的細節取決于其后期激發出的K-H不穩定性的強弱.無粘情況下，K-H不穩定性的強弱取決于所采用的數值格式的數值粘性，也就是格式的精度和網格的數量;有粘情況下，隨著精度的提高和網格的細化，當數值粘性遠小于物理粘性時，流場的細微結構唯一取決于流體的物理粘性，所得到的數值解是可信的.

［1］王繼海.二維非定常流和激波［M］.北京:科學出版社，1994:348-471.WANG Jihai.Two dimensional unsteady flows and shock waves［M］.Beijing:Science Publishing Press，1994:348-471.

［2］HARTEN A，ENGQUIST B，OSHER S，et al.Uniformly high order essentially non-oscillatory schemes［J］.JComput Phys，1987，71:231-303.

［3］JIANG G S，SHU CW.Efficient implementation of weighted ENO schemes［J］.J Comput Phys，1996，126:202-228.

［4］LIU X D，OSHER S，CHAN T.Weighted essentially nonoscillatory schemes［J］.JComput Phys，1994，115:200-232.

［5］TORO E F.Riemann solvers and numericalmethods for fluid dynamics［M］.3rd Ed.Berlin:Springer，2009:531-542.

［6］KAMANISN A，DOUGALIS V A，EKATERINARIS JA.Effective computationalmethods for wave propagation［M］.NewYork:Chapman＆ Hall/CRC，2008:533-540.

［7］BALSARA D S，SHU CW.Monotonicity preserving weighted essentially nonoscillatory schemes with increasingly high order of accuracy［J］.J Comput Phys，2000，160:405-452.

［8］SHI J，ZHANG Y T，SHU CW.Resolution of high order WENO schemes for complicated flowstructures［J］.JComput Phys，2003，186:690-696.

［9］SAMTANEY R，PULLIN D I.On initial-value and self-similar solutions of the compressible Euler equations［J］.Physics of Fluids，1996，8:2650-2655.

［10］ZHANG Y T，SHIJ，SHU CW，etal.Numerical viscosity and resolution of high order weighted essentially non-oscillatory schemes for compressible flowswith high Reynolds numbers［J］.Physical ReviewE，2003，68:046709.