李 虎 羅 勇 劉旭亮 武從海 韓帥斌 王益民

(中國空氣動力研究與發(fā)展中心,空氣動力學國家重點實驗室,四川綿陽 621000)

(中國空氣動力研究與發(fā)展中心,計算空氣動力研究所,四川綿陽 621000)

引言

氣動噪聲是流體非定常運動引起的,主要有以下特征: 一是微尺度特性,近場聲擾動的幅值比流體壓力脈動的幅值小幾個數量級[1-2];二是多尺度特性,湍流結構產生的噪聲具有很寬的頻率范圍;三是聲波近似等熵,無耗散和無色散,傳播距離很遠.氣動噪聲的微尺度、多尺度和近似等熵特性對數值模擬中的離散格式要求特別高,需要格式具有高階精度以及寬頻域的低耗散和低色散特性.

激波噪聲是由流場中的湍流結構和激波的相互作用產生.因此,對于激波噪聲問題,還要求數值格式具有魯棒的激波捕捉能力.然而,同時捕捉激波和聲波本質上是不相容的,因為前者的模擬需要消除激波周圍的偽振蕩波,而后者的模擬又需要準確分辨所有聲波.到目前為止,還不存在適合于計算激波噪聲的理想數值格式.為了實現激波噪聲的精準模擬和預測,需要發(fā)展具有高階精度、低色散、低耗散和魯棒激波捕捉能力的空間離散格式.

加權本質無振蕩(weighted essentially nonoscillatory,WENO)格式[3]是一種經典的高階精度激波捕捉格式,以WENO-JS 格式[4]為代表,廣泛應用于包含激波和復雜光滑結構可壓縮流動問題的數值模擬[5-9].雖然WENO-JS 格式能夠捕捉強激波,但對于氣動聲學問題耗散過大.針對WENO-JS 格式耗散過大、極值點附近精度降階和計算定常激波不收斂等問題,學者們發(fā)展了一系列改進格式,如WENO-M[10],WENO-Z[11],WENO-S[12],WENOSYMBO[13],WENO-SYMCU[14],WENO-ZS[15]和TENO[16]等.

緊致格式也是一種典型的高精度高分辨率數值方法.Lele[17]系統研究了兩類線性中心緊致格式: 網格結點型緊致格式(cell-noded compact schemes,CNCS)和網格中心型緊致格式(cell-centered compact schemes,CCCS).CNCS 格式具有零耗散和低色散特性,廣泛應用于不含激波的湍流和氣動聲學問題數值模擬.然而,在實際的計算中,零耗散會導致計算不穩(wěn)定,需引入人工耗散或者濾波.

為了能夠捕捉激波,基于Lele[17]提出的網格中心型緊致格式,結合WENO 格式的激波捕捉思想,Deng 等[18]設計了五階精度的加權緊致非線性格式(weighted compact nonlinear scheme,WCNS).WCNS 格式具有與WENO-JS 格式相當的激波捕捉能力和較小的耗散,以及更靈活的通量選擇.Zhang等[19]將WCNS 格式拓展到了更高階精度,并且在構造網格中心處的數值通量時,是直接對通量進行插值,而不是對守恒變量或原始變量進行插值后再構造通量.

CNCS 格式和CCCS 格式在構造過程中,都只用了設計模板上的部分信息,格式精度和分辨率不能達到最優(yōu).Liu 等[20]設計了一類新型線性中心緊致格式(cemtral compact schemes with spectral-like resolution,CCSSR),該格式在求解網格結點上的函數導數時,將網格結點上的函數值和網格中心處的函數值全部引入,精度最高可達14 階,并且接近譜方法的分辨率.隨后,劉旭亮等[21]以六階精度的線性中心緊致格式(CCSSR-L-6) 為基礎,通過迎風/對稱混合型加權非線性插值(WENO 插值)計算網格中心處的數值通量,構造了能夠捕捉強激波的迎風/對稱混合型加權非線性緊致格式,記為CCSSRHW-6 格式.

為了能夠捕捉強激波,增強格式在計算時的穩(wěn)定性,CCSSR-HW-6 格式在其構造網格中心處數值通量的過程中包含兩級加權: 一是五點迎風模板插值和六點對稱模板插值的加權;二是內部子模板插值的加權.兩級加權的權系數都不是固定值,而是非線性的函數.兩級加權雖然增強了格式的耗散自適應能力和激波捕捉能力,但也會使格式形式變得復雜,非線性效應增強.數值格式的非線性效應會導致流場中產生非物理的寄生波[22-24],嚴重影響氣動噪聲問題的高保真計算.在大多數激波噪聲問題中,激波的強度并不是很高,CCSSR-HW-6 格式的超強激波捕捉能力并不總是必需的,這就為一定程度上放寬激波捕捉能力的要求,進而降低格式的非線性效應和提高格式的分辨能力留出了空間.

本文通過優(yōu)化CCSSR-HW-6 格式中第一級加權的自由權系數,使其取能夠讓格式空間分辨能力達到最優(yōu)的固定值,從而簡化格式的構造,減小格式的耗散誤差和非線性效應,使其更適合于激波噪聲問題的數值模擬.

1 加權優(yōu)化緊致格式設計

1.1 迎風/對稱混合型加權非線性緊致格式

現有的CCSSR-HW-6 格式是在具有類譜分辨率的六階精度線性中心緊致格式(CCSSR-L-6)的基礎上,結合迎風/對稱混合型WENO 插值技術發(fā)展而來.圖1 是格式構造模板.具體構造過程[21]如下.

圖1 迎風/對稱混合型加權非線性緊致格式的構造模板Fig.1 Candidate stencils of hybrid upwind/symmetric weighted nonlinear compact scheme

六階精度的線性中心緊致格式(CCSSR-L-6)的表達式為

在內部點上,采用隱式格式,系數取值為

在邊界點上,采用顯式格式,系數取值為

在CCSSR-L-6 格式中,網格中心上的函數值或數值通量是未知的.為了能夠捕捉激波,基于網格結點上的函數值,采用迎風/對稱混合型WENO 插值得到網格中心上的函數值.迎風/對稱混合型WENO插值技術包含兩級加權: 一是五點迎風模板插值和六點對稱模板插值的加權;二是內部四個子模板插值的加權.

(1) 第一級加權是五點迎風模板插值和六點對稱模板插值的加權,具體如下

(2) 第二級加權是四個內部子模板上插值的加權,具體如下

其中,dr是子模板上的第二級權系數.

子模板上的3 階精度線性插值多項式為

聯立式(4)和式(7)可建立兩級加權的權系數之間的關系

以上內容即為CCSSR-HW-6 格式的線性部分.

(3) 非線性權的構造,具體如下:

第一級權系數 σ 定義如下[21]

式中,rc是常數閾值,其值越大,數值格式的耗散越大.CCSSR-HW-6 格式將其取值為rc=0.5.

相應的非線性權[21]為

式中,ε是小值正數,作用是避免分母為零,具體取值為ε=1.0×10-40.引入參數C的目的是增大最優(yōu)線性權的貢獻,其值的增大或減小僅會改變數值格式的耗散誤差.參數C的值越大,數值格式的耗散誤差越小.在處理含有激波的問題時,過大的C值可能會導致離散格式的數值不穩(wěn)定,CCSSR-HW-6 格式將參數C取值為C=5.

內部四個子模板上的光滑指示因子[21]為

參考光滑指示因子[21]定義為

此外,非線性權會導致格式在特定的極值點處損失精度.為解決CCSSR-HW-6 格式在極值點處可能存在的精度降階問題,CCSSR-HW-6 格式還引入了映射函數[10]

使得非線性權 ωr在臨界點處盡量接近線性權dr,從而還原格式的近似精度.

CCSSR-HW-6 格式的兩級加權確實增強了格式的耗散自適應能力和激波捕捉能力,但也使格式的形式變得復雜,非線性效應增強.對于大多數激波噪聲問題,激波的強度并不是很高.因此,可以通過固定CCSSR-HW-6 格式中第一級加權的權系數取值,簡化格式構造,降低格式的耗散誤差,減弱格式的非線性效應,提高格式的分辨能力.

1.2 分辨能力優(yōu)化

本文以CCSSR-HW-6 格式線性部分的修正波數為基礎,采用修正波數的誤差積分函數[13,25-26]為空間分辨能力優(yōu)化目標函數,通過優(yōu)化,固定第一級加權權系數的取值,使格式的空間分辨能力達到最優(yōu).格式優(yōu)化過程具體如下.

第一步,推導CCSSR-HW-6 內點格式和邊界格式線性部分的無量綱修正波數.

CCSSR-HW-6 格式線性部分無量綱修正波數的計算方法如下:

考慮純粹的諧波函數

這里,x是位置,k是波數.

諧波函數的解析導數為

定義n為任意整數,Δx是網格間距,則有

一般有限差分方法的近似導數可表示為

其中,系數an是由差分格式中的系數決定.

定義無量綱波數為 κ=kΔx,結合式(16)～式(18)可得無量綱修正波數的一般表達式

利用上述方法可以得到CCSSR-HW-6 格式線性部分的無量綱修正波數.

對于CCSSR-HW-6 內點格式,無量綱修正波數表達式為

式(20)中的系數如下

對于CCSSR-HW-6 邊界格式,無量綱修正波數表達式為

式中的系數如下

第二步,確定空間分辨率優(yōu)化目標函數.

數值格式修正波數的實部表征了色散誤差(相位誤差),虛部則表征了耗散誤差(幅值誤差).為了對CCSSR-HW-6 格式進行優(yōu)化,確定空間分辨率優(yōu)化目標函數為修正波數的誤差積分函數.修正波數的誤差積分函數綜合考慮了數值格式的色散誤差和耗散誤差.

無量綱修正波數的誤差積分函數[13,25-26]表達式為

其中,χ,ν ,ξ 和 η 是優(yōu)化控制參數.

第三步,設置恰當的優(yōu)化控制參數.

優(yōu)化控制參數在誤差積分函數中分別具有不同的作用[13]: 參數χ控制相位誤差和幅值誤差的權值;參數 ξ 控制了為確保格式穩(wěn)定而增加耗散的正弦項;參數 ν 影響數值格式在低波數區(qū)的誤差;參數 η 影響數值格式在高波數區(qū)的誤差.

在本文中,優(yōu)化控制參數的取值與文獻[13]相同,分別為第四步,執(zhí)行優(yōu)化程序,得到優(yōu)化結果.

選擇優(yōu)化算法,執(zhí)行優(yōu)化程序,確定第一級加權的權系數取值,使得修正波數的誤差積分函數取最小值.

第一級加權的權系數優(yōu)化結果:

內點格式優(yōu)化系數

邊界格式優(yōu)化系數

經過空間分辨能力優(yōu)化后的迎風/對稱混合型加權非線性緊致格式稱之為加權優(yōu)化緊致格式(weighted optimization compact scheme,WOCS).

2 精度驗證

通過一維標量方程算例對加權優(yōu)化緊致格式(WOCS 格式)的收斂精度進行數值驗證.時間推進采用三步三階TVD Runge-Kutta 方法[4].

一維線性標量方程

初始條件為

計算采用周期邊界條件,總時長為t=1;計算網格由10 逐步加密為160;計算網格為10 時,初始CFL 數為0.5,隨著網格的每一次加密,對于r階格式,CFL 數以因子遞減.

表1 列出了WOCS 格式求解該問題時的數值誤差和精度階數.結果顯示,WOCS 格式的收斂精度高于5 階.

表1 加權優(yōu)化緊致格式的數值誤差和精度階數Table 1 Numerical errors and accuracy order of the weighted optimization compact scheme

3 譜特性分析

數值格式截斷誤差的精度階數只能提供數值解向精確解的漸進收斂速率信息,不能給出有限計算網格上的實際誤差[27];而差分格式的波傳播特性(譜特性)能給出計算網格所支持的Fourier 模態(tài)的演化[27].

對于波數為k的單一Fourier 模態(tài),非線性格式不僅會在相同的波數k上產生線性響應,還會在其他波數上產生非線性響應[24,28].對于非線性格式不存在解析的譜關系式,將廣泛應用于線性格式的Fourier 分析直接應用到非線性格式的空間分辨能力分析是行不通的.針對非線性格式發(fā)展的近似色散關系(approximate dispersion relation,ADR)技術[27]實際上得到的只是非線性格式修正波數的線性響應部分.基于空間導數的修正波數計算方法[24,28]可以得到數值格式的線性響應和非線性響應,具體如下.

非線性格式對于波數k的修正波數為

3.1 線性響應

圖2 是WOCS 格式和原CCSSR-HW-6 格式修正波數線性響應部分的對比.線性響應的實部和虛部分別表征色散特性和耗散特性.圖2(a)顯示,WOCS格式和原CCSSR-HW-6 格式具有相同的色散誤差;圖2(b)顯示,WOCS 格式在中高波數區(qū)的耗散誤差明顯小于原CCSSR-HW-6 格式.總的來說,通過優(yōu)化使CCSSR-HW-6 格式第一級加權的權系數 σ 取固定值,提高了格式在中高波數區(qū)的空間分辨能力.

圖2 修正波數的線性響應部分Fig.2 The linear response part of the modified wavenumber

3.2 非線性響應

圖3 是WOCS 格式和原CCSSR-HW-6 格式的修正波數非線性響應部分的Fourier 系數對比.計算網格點數為N=64,在該網格所支持所有 Fourier 模態(tài)中,模態(tài)5、模態(tài)10 和模態(tài)15 是屬于中低波數的模態(tài)(模態(tài)k表示該模態(tài)的波數是k).對于任意模態(tài)k,非線性格式除了在波數k上會產生線性響應,還會在其他波數 (3+2n)kmodN(N是整數)上產生非線性響應(非物理的寄生波)[24].非線性響應的實數部分代表寄生波的波數響應,虛數部分代表寄生波的幅值響應.圖3 顯示,在中低波數區(qū)間,WOCS格式非線性響應的實數部分與原CCSSR-HW-6 格式相當,而非線性響應的虛數部分則明顯弱于原CCSSR-HW-6 格式.綜合來看,與原CCSSR-HW-6 格式相比,WOCS 格式降低了中低波數區(qū)的非線性響應.

圖3 修正波數的非線性響應部分Fig.3 The nonlinear response part of the modified wavenumber

4 數值實驗

本節(jié)針對加權優(yōu)化緊致格式(WOCS 格式)在激波噪聲問題計算中的適用性進行測試,包括激波捕捉能力、分辨能力和非線性效應.數值實驗包括一維Euler 方程算例、二維Euler 方程算例和二維Navier-Stokes 方程算例.無黏數值通量計算采用HLL 型Riemann 求解器[29]或Lax-Fredrichs 通量分裂方法[4],時間推進采用三步三階TVD Runge-Kutta 方法[4].

4.1 一維Euler 方程算例

(1)激波-聲波相互作用

激波-聲波相互作用[30]是氣動聲學中的重要基準問題之一.計算條件如下.

計算域為[0,1],計算網格數為N=401,CFL 數為0.5,計算時長為t=30Tλ;初始時刻,激波位置為xs=0.5,激波馬赫數為Ms=2,氣體從左向右流動,激波左右兩側的流動變量滿足激波關系式.計算域右邊界設置為無反射邊界條件,在左邊界處添加聲擾動

式中參數取值如下

圖4 是激波-聲波相互作用的數值解及其與參考解的比較.由于解析解是未知的,這里選取五階精度WCNS-E-5 格式[18]在計算網格數為N=1001 時求解該問題得到的數值解作為參考解.圖4(b) 和圖4(c) 是圖4(a) 的局部放大圖.圖中的縱坐標δP(t)=P(x,t)-P(x,0)表示聲擾動.聲波穿過激波后其幅值會被放大.結果顯示: 在分辨穿過激波后的聲波時,WOCS 格式計算得到的聲波幅值與參考解更為吻合,其分辨能力優(yōu)于原CCSSR-HW-6 格式.

圖4 激波-聲波相互作用問題的數值結果Fig.4 Numerical results of shock-sound interaction problem

(2)激波-密度波相互作用

激波-密度波相互作用,即Titarev-Toro 問題[31],也是驗證數值格式激波捕捉能力和流場細微結構分辨能力的標準算例,在Titarev-Toro 問題中,向右運動的激波(馬赫數為1.1)與正弦波形式的高頻密度擾動相互作用,產生的流場既包含光滑解,也包含間斷.該問題可以看作是激波/湍流相互作用的簡化模型.初始條件如下

計算條件: 計算網格數為N=1001,CFL 數為0.5,計算時間為t=5.0;邊界條件設置為計算域左、右邊界處的守恒變量不隨時間變化,這在本算例所感興趣的計算時長內是適宜的.

圖5 是Titarev-Toro 問題的數值解及其與參考解的比較.由于解析解是未知的,這里選取五階精度WCNS-E-5 格式[18]在計算網格數為N=10 001 時求解該問題得到的數值解作為參考解.圖5(b)和圖5(c)是圖5(a)的局部放大圖.結果顯示: WOCS 格式數值解的幅值與參考解更為接近,其對高頻波的分辨能力遠優(yōu)于原CCSSR-HW-6 格式.這是因為相較于原CCSSR-HW-6 格式,WOCS 格式顯著降低了高波數區(qū)的耗散誤差,如圖2(b)所示.

圖5 Titarev-Toro 問題的數值結果Fig.5 Numerical results of Titarev-Toro problem

4.2 二維Euler 方程算例

選擇激波-渦量波相互作用[32-33]作為二維Euler方程測試算例,該問題是激波-湍流相互作用的二維簡化模型.在此算例中,馬赫數為Ms=8 的運動激波與散度為零的傾斜渦量場發(fā)生相互作用.計算域范圍: (x,y)∈[-1.5,1.5]×[-1,1],均勻計算網格,網格間距為 Δx=Δy=0.01 ;計算時長t=0.2.

初始時刻,激波位于x=-1 處,傳播速度為

這里的cs1是波前聲速(下標1 表示波前狀態(tài),下標2 表示波后狀態(tài)).

波前(x≥-1)初始狀態(tài)為

波后(x≤-1)初始狀態(tài)為

邊界條件設置: 在y方向,上、下兩端采用周期邊界條件;在x方向,左端為超聲速入流條件,右端為無反射邊界條件.

圖6 是采用WOCS 格式計算激波-渦量波相互作用得到的t=0.2 時刻的渦量空間分布.圖7 是t=0.2時刻直線y=0 上的渦量分布及其與參考解的比較.參考解是在計算網格數為N=961 × 641 條件下采用CCSSR-HW-6 格式計算該問題時得到的數值解.結果顯示: WOCS 格式的數值解與參考解要更為接近,分辨能力優(yōu)于原CCSSR-HW-6 格式.

圖6 激波-渦量波相互作用在t =0.2 時刻的渦量空間分布Fig.6 The spatial distribution of vorticity in the shock-vorticity wave interaction at t =0.2

圖7 激波-渦量波相互作用在t =0.2 時刻沿直線y=0 的渦量分布Fig.7 The vorticity distribution along the line y=0 in the shockvorticity wave interaction at t =0.2

4.3 二維Navier-Stokes 方程算例

二維Navier-Stokes 方程測試算例選取的是激波-強旋渦相互作用問題[6,34].激波-旋渦相互作用是重要的激波噪聲理論研究模型,文獻[6,34]對此問題進行了深入細致的研究.此外,該問題也被廣泛用作高精度數值方法的驗證算例[35].下面簡要介紹該問題計算細節(jié).

計算域范圍: (x,y)∈[-20,8]×[-12,12],均勻計算網格,網格點數為N=701 × 601,計算總時長t=10.普朗特數Pr=0.75,參考溫度T∞=300 K .雷諾數Re=800,具體定義為Re=ρ∞a∞R/μ∞,其中,R是渦核半徑,ρ∞,a∞和 μ∞分別是激波前流動的平均密度、聲速和黏性系數.初始條件設置如下.

激波馬赫數為Ms=1.2,固定在xs=0 的位置.激波前(x≥0)的狀態(tài)參數為

激波后(x＜0)的狀態(tài)參數為

初始時刻,旋渦位于激波波前,其從右向左穿越激波的過程中與之發(fā)生相互作用,產生并向外輻射聲波.初始渦心位置是 (xv,yv)=(2,0),旋渦強度為Mv=1.0.旋渦的速度、密度和壓力分布如下

邊界條件設置: 在y方向,上、下兩端采用周期邊界條件;在x方向,右端為超聲速入流條件,左端為無反射邊界條件.

對于氣動聲學問題,脹量(速度的散度,?·V)和數值陰影(?2ρ)常用來表征聲場和流動結構.圖8是采用WOCS 格式計算激波-強旋渦相互作用得到的t=10 時刻的脹量空間分布,從中可以看到向外輻射的聲波.圖9 是t=10 時刻脹量和數值陰影沿徑向(θ=π/4,θ 的定義見圖8)的分布及其與參考解的比較.參考解是在計算網格數為N=1401 × 1201條件下采用CCSSR-HW-6 格式計算該問題時得到的數值解.對于脹量,WOCS 和CCSSR-HW-6 兩種格式的計算結果一致;而對于數值陰影,WOCS 格式的計算結果優(yōu)于原CCSSR-HW-6 格式的計算結果,其與參考解要更為接近.

圖8 激波-強旋渦相互作用在t=10 時刻的脹量場Fig.8 The dilatation field of shock-strong vortex interaction at t=10

圖9 激波-強旋渦相互作用在t =10 時刻的脹量和數值陰影沿徑向的分布Fig.9 The radial distribution of dilatation and numerical shadowgraph in shock-strong vortex interaction at t=10

圖10 是分別采用WOCS 格式和CCSSR-HW-6 格式計算激波-強旋渦相互作用得到的t=10 時刻的數值陰影空間分布.結果顯示: 在CCSSR-HW-6格式計算得到的數值陰影中存在非常明顯的非物理振蕩,而WOCS 格式計算得到的數值陰影則更為“干凈”,上述非物理振蕩被顯著地減弱.文獻[22-24]指出激波捕捉格式在計算時的非線性執(zhí)行會導致非線性效應,使得數值解中產生此類非物理波動.前文3.2 節(jié)中的格式非線性響應分析已經表明,與原CCSSR-HW-6 格式相比,經優(yōu)化后,WOCS 格式的非線性響應顯著降低,由此可以解釋這兩種格式計算結果之間的差異.

圖10 激波-強旋渦相互作用在t=10 時刻的數值陰影Fig.10 Numerical shadowgraph of shock-strong vortex interaction at t=10

5 結論

基于對稱模板的迎風/對稱混合型加權非線性緊致格式在其構造過程中包含兩級加權,并且兩級加權的權系數都是非線性函數,格式非線性效應較強.為了減弱格式的非線性效應,改進格式的譜特性,本文基于六階精度的迎風/對稱混合型加權非線性緊致(CCSSR-HW-6)格式,以修正波數的誤差積分函數為優(yōu)化目標函數,通過優(yōu)化第一級加權的權系數,使其取分辨能力最優(yōu)的固定值,建立了加權優(yōu)化緊致格式,記為WOCS 格式.針對WOCS 格式的精度、譜特性和激波噪聲問題計算的適用性進行了數值驗證.

精度驗證表明WOCS 格式的精度高于5 階.譜特性分析包括了線性響應和非線性響應.線性響應分析表明WOCS 格式在中高波數區(qū)的耗散特性明顯優(yōu)于CCSSR-HW-6 格式.非線性響應分析表明:與CCSSR-HW-6 格式相比,WOCS 格式降低了中低波數區(qū)的非線性響應.

激波噪聲計算適用性的測試算例包括了激波-聲波相互作用、激波-密度波相互作用、激波-渦量波相互作用以及激波-強旋渦相互作用.數值實驗結果表明: 與原CCSSR-HW-6 格式相比,WOCS 格式不僅降低了耗散誤差,提高了對高頻波的分辨能力,而且顯著地減弱了格式的非線性效應和由此導致的數值解的非物理振蕩,較為適合于激波噪聲問題的數值模擬.