利用高精度迎風(fēng)緊致群速度控制法模擬激波與柱形界面相互作用*

姚清河，陳耀欽，薛 莉，朱慶勇

(中山大學(xué)應(yīng)用力學(xué)與工程系，廣東 廣州 510275)

眾所周知，對于氣動方程，不管初始值如何光滑，其解都可能出現(xiàn)間斷,采用高于一階精度格式求解激波問題時，其數(shù)值解在激波附近都會產(chǎn)生非物理振蕩。

在流場激波捕捉技術(shù)中，已經(jīng)發(fā)展了許多行之有效的方法，如TVD 格式、ENO 格式、WENO 格式、CSCM格式、NND 格式、耗散比擬方法等[1-7],這些格式都有較高的精度和激波分辨率。群速度控制法從物理角度出發(fā)分析非物理振蕩產(chǎn)生原因[8-9]，并提出改進激波解的辦法。在這個思想基礎(chǔ)上，本文從五階迎風(fēng)緊致格式出發(fā)，引進一個函數(shù)用以控制數(shù)值解的群速度，構(gòu)造帶有群速度控制的五階迎風(fēng)緊致格式。該格式在激波前后表現(xiàn)為不同的性質(zhì)，使得激波前后的振蕩向激波靠攏，從而達到消除激波振蕩的目的。同時，文中進一步對所構(gòu)造的格式的精度及數(shù)值解的行為進行了分析。

對于接觸間斷兩邊的流體，一方面以上格式始終存在數(shù)值耗散而導(dǎo)致接觸間斷被抹平，另一方面Euler 方法在模擬時需要進行特別的處理，如MAC，VOF和Level Set法等。本文在前人基礎(chǔ)上，采用五階迎風(fēng)緊致格式求解描述界面的Level Set 方程，采用Ghost法處理界面附近的邊界條件[10]。本方法在激波與柱形界面相互作用的計算中取得良好的結(jié)果，數(shù)值結(jié)果與實驗是吻合的。

1 控制方程

1.1 流場控制方程

針對激波與流體界面相互作用的問題，由歐拉坐標(biāo)系中二維非定?？蓧嚎s守恒性Euler 方程來描述流場。具體地

(1)

其中

,,

流體的狀態(tài)方程為

p=(γ-1)(E-ρ/2(u2+v2))

其中，ρ是密度，u,v分別是x,y方向上的速度分量，p是流體壓力，E是單位體積流體的總能量，e是比內(nèi)能。

1.2 Level Set 方程

二維界面Level Set方程

φt+uφx+vφy=0

(2)

其中，φ(x,t)是x到界面Γ(t)的符號距離，u=(u,v)是流體速度，t是時間。但是，由于數(shù)值方法的內(nèi)在效應(yīng)，即使只進行了幾個時間步求解，φ(x,t)將不再滿足符號距離函數(shù)的定義了。為了保持這一良好性質(zhì)，我們采用一種所謂重新初始化(Reinitialization)的手段[11]，也就是要改造φ(x,t)使重新成為x到界面Γ(t)的符號距離。重新初始化方程如下

(3)

2 數(shù)值方法

2.1 標(biāo)量方程的空間離散方法

首先考慮如下的模型方程及對應(yīng)的半離散化方程

,f=cu,c=const

(4)

(5)

這里Fj/Δx是一階導(dǎo)數(shù)?f/?x的差分逼近?？紤]取如下五個網(wǎng)格點相聯(lián)系的差分逼近[11]

ω1Fj+1+ω0Fj+ω-1Fj-1=β-2fj-2+β-1fj-1+

β0fj+β1fj+1+β2fj+2

(6)

,,，

這樣(6)式可寫為

(7)

Δx5)

(8)

即逼近精度的量級為Δx5。

簡單采用(7)式,該方法經(jīng)過激波時存在非物理振蕩。本文采用如下群速度控制方法以克服激波附近非物理振蕩，提高激波分辨率，同時提高了光滑區(qū)的計算精度。本文在(5)式右端加入修正項以控制其群速度。

(9)

(10)

1)情形一，ssj+1=ssj-1=1：

σ(-cos3α+6cos2α-15cosα+10)，

σ(3cos3α-8cos2α+5cosα)

我們有kr≥0,?α∈[0,π]。D(α)≥1,?α∈[0,α*]，這里0<α*<π。

2)情形二，ssj+1=1,ssj-1=-1：

σ(2cos2α-8cosα+6)，

我們有kr≥0,?α∈[0,π]。

3)情形三，ssj+1=-1,ssj-1=1：

σ(-2cos3α+10cos2α-22cosα+14)，

我們有kr≥0,?α∈[0,π]。

4)情形四，ssj+1=-1,ssj-1=-1：

σ(-cos3α+6cos2α-15cosα+10)，

σ(-3cos3α+8cos2α-5cosα)

我們有kr≥0,?α∈[0,π]。D(α)≤1,?α∈[0,α*]，這里0<α*<π。

在上述四種情況中，情形一和四的范數(shù)ki-αL2大于情形二、三的范數(shù)ki-αL2。情形二、三的范數(shù)ki-αL2相等。

根據(jù)以上分析，我們可以得到以下結(jié)論，在情形一中，Lg能使得振蕩波向激波靠攏；在情形二、情形三中，F(xiàn)j+σLg都是正耗散，有利于振蕩的消除，F(xiàn)j的群速度不受影響；情形四中，Lg也能使得振蕩波向激波靠攏。

帶群速度控制的五階迎風(fēng)緊致格式Fj+σLg的色散性質(zhì)、耗散性質(zhì)和群速度，分別如圖1，圖2，圖3所示(σ=0.5)。

圖1 帶群速度控制的五階迎風(fēng)緊致格式的色散性質(zhì)

圖2 帶群速度控制的五階迎風(fēng)緊致格式的耗散性質(zhì)

圖3 帶群速度控制的五階迎風(fēng)緊致格式的群速度

2.2 二維空間離散方法

在對實際流動問題進行數(shù)值模擬時，往往會碰到區(qū)域邊界形狀復(fù)雜和物理變量在區(qū)域內(nèi)變化大的情況。這時，為了較好的進行模擬，需要先建立曲線坐標(biāo)網(wǎng)格，即把不規(guī)則物理域先映射到規(guī)則計算域，然后才能在計算域上進行差分計算。當(dāng)然，相應(yīng)的流體力學(xué)方程組也需要變換到計算域。

假設(shè)物理域上的(t,x,y)到矩形計算域上的(τ,ξ,η)之間的坐標(biāo)變換為

(11)

則可以把方程(1)變換成如下方程

(12)

其中

ξxηy-ξyηx，/J，

(13)

(14)

特征值矩陣為

采用Steger-Warming 通量分裂技術(shù)

±ε2)1/2)/2],l=1,2,3,4；

逼近(12)式且具有群速度控制項的半離散差分格式為

(15)

2.3 Level Set方程的離散方法

本文采用五階迎風(fēng)緊致格式對二維界面Level Set方程(2)式進行空間離散。重新初始化方程是一個Hamilton-Jacobi方程，方程(3)可以寫為如下形式：

φt+H(φx,φy)=0

(16)

方程(16)的半離散形式如下

(17)

(18)

其中

,,

I(a,b)=[min(a,b),max(a,b)],

(19)

本文利用Ghost方法將Level Set 函數(shù)與物理量的控制方程進行耦合，詳見文[10]。本文在時間方向采用三階R-K方法。

3 數(shù)值結(jié)果

一個Mach數(shù)為1.22的激波通過一個氣泡相遇，無量綱計算區(qū)域(0,325)×(-44.5,44.5)。網(wǎng)格為325×81，無量綱化初值條件是

x>225:ρ=1.374 6,u=-0.394,

v=0,p=1.569 8；

x≤225:ρ=1,u=0,v=0,p=1；

(x-175)2+y2≤252:ρ=3.153 8,

u=0,v=0,p=1

圖4 初始時刻流場示意圖

Level Set 方程的初值條件為

邊界條件：左邊界取出流邊界條件，右邊界取入流邊界條件，即給定激波后的壓力、密度和速度值，上、下邊界取固壁邊界條件。由圖5給出t=18時刻的密度分布圖并與實驗結(jié)果相比較。圖6給出t=38時刻的密度分布圖并與實驗結(jié)果進行比較。

圖5 t=18時刻的等密度分布圖

圖6 t=38時刻的等密度分布圖

圖7、圖8分別給出不同時刻的密度分布圖。本文與Hass[12]等人的實驗結(jié)果相比較發(fā)現(xiàn)，模擬出來的氣泡變形與實驗圖形基本相符合，當(dāng)激波與氣泡相互作用時，產(chǎn)生了折射波和反射波，在激波掃后界面被壓扁,并隨時間的推移，氣泡的右邊基本保持原來的形狀，而左邊凹陷。本方法在激波與柱形界面相互作用的計算中取得良好的結(jié)果，數(shù)值結(jié)果與實驗是吻合的。

圖7 t=108時刻的等密度分布圖

圖8 t=320時刻的等密度分布圖

4 結(jié) 論

本文基于迎風(fēng)緊致群速度控制方法求解流場控制方程，采用群速度方法克服激波附近非物理振蕩，采用Level Set方法追蹤運動界面。本方法在激波與柱形界面相互作用的計算中取得良好的結(jié)果，數(shù)值結(jié)果與實驗是吻合的。

參考文獻：

[1]YEE H C,WARMING R F,HARTEN A.Implicit total variation diminishing TVD schemes for steady-state calculations [J].Journal of Computational Physics,1985,57(3): 327-360.

[2]HARTEN A,ENGQUIST B,OSHER S,et al.Uniformly high order accurate essentially non-oscillatory schemes III [J].Journal of Computational Physics,1987,71(2):231-303.

[3]SHU C W.Numerical experiments on the accuracy of ENO and modified ENO schemes [J].Journal of Scientific Computing,1990,5(2): 127-149.

[4]LIU X D,OSHER S,CHAN T.Weighted essentially non-oscillatory schemes [J].Journal of Computational Physics,1994,115 (1): 200-212.

[5]LOMBARD C K,BARDINA J,VENKATAPATHY E,et al.Multi-dimensional formulation of CSCM—an upwind flux difference eigenvector split method for the compressible Navier-Stokes equations [C]// New York: American Institute of Aeronautics and Astronautics,1983: 649-664.

[6]張涵信.無波動無自由參數(shù)的耗散差分格式 [J].空氣動力學(xué)報,1988,2: 143-165.

[7]馬延文,傅德薰.計算空氣動力學(xué)中一個新的激波捕捉—耗散比擬法 [J].中國科學(xué)：A輯(數(shù)學(xué)),1992,35(3): 263-271.

[8]FU D X,MA Y W.A high order accurate difference scheme for complex flow fields [J].Journal of Computational Physics,1997,134: 1-15.

[9]ZHU Q Y,LI Y.An upwind compact approach with group velocity control for compressible flow fields [J].International Journal for Numerical Methods in Fluids,2004,44: 463-482.

[10]FEDKIW R,ASLAM T,MERRIAN B,et al.A non-oscillatory Eulerian approach to interfaces in multimaterial flow (the ghost fluid method ) [J].Journal of Computational Physics,1999,152(2): 457-492.

[11]傅德薰.流體力學(xué)數(shù)值模擬 [M].北京: 國防工業(yè)出版社,1994.

[12]HASS J F,STURTEVANT B.Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities [J].Journal of Fluid Mechanics,1987,181: 41-76.