高階WENO格式數(shù)值粘性對(duì)模擬R-T不穩(wěn)定性的影響

王革，謝昌坦，張斌

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

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

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

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

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

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

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

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

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

2 數(shù)值算例及結(jié)果分析

本算例的計(jì)算區(qū)域?yàn)椋?，0.25］×［0，1.0］，計(jì)算網(wǎng)格為方形網(wǎng)格，其寬度用h表示.初始時(shí)刻界面位于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è)為反射邊界條件，上下邊界設(shè)為無反射邊界條件，上邊界參數(shù)為 ρ=1.0，u=0，p=2.5，v=0，下邊界參數(shù)為 ρ=2.0，u=0，p=1，v=0.計(jì)算至?xí)r刻 t=1.95.

2.1 無粘情況

無粘情況下的數(shù)值流場是通過求解Euler方程得到的.圖1為不同精度WENO格式和不同網(wǎng)格數(shù)量條件下區(qū)域［0，0.25］×［0.2，0.9］的 15 條密度等值線分布圖.為簡便，引入符號(hào)WENO-N-X，其中，N代表WENO格式的精度，N=3、5、7、9;X 代表計(jì)算網(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ù)量對(duì)其數(shù)值解有很大影響.WENO-3-C、WENO-3-m和 WENO-5-C的數(shù)值結(jié)果基本相同.在重力作用下，下層重流體以尖釘?shù)男问角秩肷蠈虞p流體，而上層輕流體以氣泡的形式侵入下層重流體，兩層流體交界面形成了清晰的蘑菇狀結(jié)構(gòu)，這種結(jié)構(gòu)也是R-T不穩(wěn)定性發(fā)展到后期的典型結(jié)構(gòu).由于此時(shí)數(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)定性增強(qiáng)，蘑菇的傘狀結(jié)構(gòu)在輕流體的作用下翻轉(zhuǎn)，并在K-H不穩(wěn)定性的作用下斷裂、碎化，輕重流體混合加劇，導(dǎo)致傘狀結(jié)構(gòu)內(nèi)部的流場極其復(fù)雜，兩層流體的交界面變得模糊.蘑菇的桿狀結(jié)構(gòu)處也出現(xiàn)了渦狀結(jié)構(gòu)，但因渦的形狀比較大、數(shù)量少且排列有序，流體交界面仍然清晰可辨.高精度密網(wǎng)格情況下，數(shù)值流場的對(duì)稱性結(jié)構(gòu)被破壞，如WENO-7-F和WENO-9-F所示，這完全是由Euler方程本身造成的，此時(shí)所得的數(shù)值解是非物理的［8-9］.

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

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

2.2 有粘情況

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

圖2 Re=30 000時(shí)密度等值線分布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í)數(shù)值解收斂，那么當(dāng)Re=30 000時(shí)，WENO-5-M、WENO-7-M、WENO-9-m及 WENO-N-F的數(shù)值流場已基本穩(wěn)定;Re=60 000時(shí)，WENO-7-M、WENO-9-m和 WENO-5-F、WENO-7-F、WENO-9-F的數(shù)值解已基本收斂，這說明此時(shí)數(shù)值格式的數(shù)值粘性已經(jīng)遠(yuǎn)遠(yuǎn)小于流體的物理粘性了.同一雷諾數(shù)條件下，不同精度WENO格式所獲得的穩(wěn)定流場的結(jié)構(gòu)基本相同，且格式的精度越高，獲得穩(wěn)定解所需的網(wǎng)格越少.同一精度條件下，雷諾數(shù)越高，流體的物理粘性越小，則掩蓋數(shù)值粘性的影響所需的網(wǎng)格也就越多，穩(wěn)定流場的結(jié)構(gòu)越復(fù)雜.

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

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

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

圖5 Re=60 000時(shí)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時(shí)WENO-N-X的L1誤差(×102)Table 2 The L1 error of WENO-N-X when Re=60 000

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

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

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

3 結(jié)論

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

［1］王繼海.二維非定常流和激波［M］.北京:科學(xué)出版社，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.