王革，謝昌坦，張斌

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

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

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

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

1 控制方程及其數(shù)值解法

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

式中:ρ、p分別表示流體的密度和壓力，u、v分別表示流體沿 x、y方向的速度分量，總能)，γ為比熱比，Pr為普朗特數(shù)，Re為雷諾數(shù).當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 數(shù)值算例及結果分析

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

2.1 無粘情況

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

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

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

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

2.2 有粘情況

圖2、3分別為雷諾數(shù)Re=30 000，60 000時的15條密度等值線分布圖.為了能更好的了解R-T不穩(wěn)定性的數(shù)值解隨網(wǎng)格數(shù)量的變化情況，取圖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不穩(wěn)定性，隨著網(wǎng)格數(shù)量的增加，WENO格式的數(shù)值粘性逐步減小，流體的物理粘性逐步占優(yōu)，L1誤差減小.若規(guī)定L1≤1.2×10-3時數(shù)值解收斂，那么當Re=30 000時，WENO-5-M、WENO-7-M、WENO-9-m及 WENO-N-F的數(shù)值流場已基本穩(wěn)定;Re=60 000時，WENO-7-M、WENO-9-m和 WENO-5-F、WENO-7-F、WENO-9-F的數(shù)值解已基本收斂，這說明此時數(shù)值格式的數(shù)值粘性已經(jīng)遠遠小于流體的物理粘性了.同一雷諾數(shù)條件下，不同精度WENO格式所獲得的穩(wěn)定流場的結構基本相同，且格式的精度越高，獲得穩(wěn)定解所需的網(wǎng)格越少.同一精度條件下，雷諾數(shù)越高，流體的物理粘性越小，則掩蓋數(shù)值粘性的影響所需的網(wǎng)格也就越多，穩(wěn)定流場的結構越復雜.

圖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的密度在不同網(wǎng)格數(shù)量條件下的分布圖，y=0.6處的密度沿x方向的分布如圖7所示.可以看出，圖6中的前3幅圖的密度分布隨網(wǎng)格的加密變化較大，如圖7(a)所示.這是因為在雷諾數(shù)很高的情況下，流體物理粘性很小，網(wǎng)格較粗時的數(shù)值粘性占絕對優(yōu)勢，初步加密網(wǎng)格后其值迅速減小但仍遠大于物理粘性.隨著網(wǎng)格的繼續(xù)加密，物理粘性的影響開始體現(xiàn)出來并且逐步占優(yōu)，流場結構隨著網(wǎng)格數(shù)量的加密變化緩慢，當網(wǎng)格加密到一定程度后，數(shù)值粘性便遠小于物理粘性了，流場結構達到穩(wěn)定狀態(tài)，如圖7(b)所示.穩(wěn)定流場的流體交界面比較平滑，蘑菇的傘狀結構內(nèi)部形成了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不穩(wěn)定性問題.數(shù)值結果表明，R-T不穩(wěn)定性的復雜流場結構的細節(jié)取決于其后期激發(fā)出的K-H不穩(wěn)定性的強弱.無粘情況下，K-H不穩(wěn)定性的強弱取決于所采用的數(shù)值格式的數(shù)值粘性，也就是格式的精度和網(wǎng)格的數(shù)量;有粘情況下，隨著精度的提高和網(wǎng)格的細化，當數(shù)值粘性遠小于物理粘性時，流場的細微結構唯一取決于流體的物理粘性，所得到的數(shù)值解是可信的.

［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.