李 虎 羅 勇 劉旭亮 武從海 韓帥斌 王益民
(中國(guó)空氣動(dòng)力研究與發(fā)展中心,空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,四川綿陽 621000)
(中國(guó)空氣動(dòng)力研究與發(fā)展中心,計(jì)算空氣動(dòng)力研究所,四川綿陽 621000)
氣動(dòng)噪聲是流體非定常運(yùn)動(dòng)引起的,主要有以下特征: 一是微尺度特性,近場(chǎng)聲擾動(dòng)的幅值比流體壓力脈動(dòng)的幅值小幾個(gè)數(shù)量級(jí)[1-2];二是多尺度特性,湍流結(jié)構(gòu)產(chǎn)生的噪聲具有很寬的頻率范圍;三是聲波近似等熵,無耗散和無色散,傳播距離很遠(yuǎn).氣動(dòng)噪聲的微尺度、多尺度和近似等熵特性對(duì)數(shù)值模擬中的離散格式要求特別高,需要格式具有高階精度以及寬頻域的低耗散和低色散特性.
激波噪聲是由流場(chǎng)中的湍流結(jié)構(gòu)和激波的相互作用產(chǎn)生.因此,對(duì)于激波噪聲問題,還要求數(shù)值格式具有魯棒的激波捕捉能力.然而,同時(shí)捕捉激波和聲波本質(zhì)上是不相容的,因?yàn)榍罢叩哪M需要消除激波周圍的偽振蕩波,而后者的模擬又需要準(zhǔn)確分辨所有聲波.到目前為止,還不存在適合于計(jì)算激波噪聲的理想數(shù)值格式.為了實(shí)現(xiàn)激波噪聲的精準(zhǔn)模擬和預(yù)測(cè),需要發(fā)展具有高階精度、低色散、低耗散和魯棒激波捕捉能力的空間離散格式.
加權(quán)本質(zhì)無振蕩(weighted essentially nonoscillatory,WENO)格式[3]是一種經(jīng)典的高階精度激波捕捉格式,以WENO-JS 格式[4]為代表,廣泛應(yīng)用于包含激波和復(fù)雜光滑結(jié)構(gòu)可壓縮流動(dòng)問題的數(shù)值模擬[5-9].雖然WENO-JS 格式能夠捕捉強(qiáng)激波,但對(duì)于氣動(dòng)聲學(xué)問題耗散過大.針對(duì)WENO-JS 格式耗散過大、極值點(diǎn)附近精度降階和計(jì)算定常激波不收斂等問題,學(xué)者們發(fā)展了一系列改進(jìn)格式,如WENO-M[10],WENO-Z[11],WENO-S[12],WENOSYMBO[13],WENO-SYMCU[14],WENO-ZS[15]和TENO[16]等.
緊致格式也是一種典型的高精度高分辨率數(shù)值方法.Lele[17]系統(tǒng)研究了兩類線性中心緊致格式: 網(wǎng)格結(jié)點(diǎn)型緊致格式(cell-noded compact schemes,CNCS)和網(wǎng)格中心型緊致格式(cell-centered compact schemes,CCCS).CNCS 格式具有零耗散和低色散特性,廣泛應(yīng)用于不含激波的湍流和氣動(dòng)聲學(xué)問題數(shù)值模擬.然而,在實(shí)際的計(jì)算中,零耗散會(huì)導(dǎo)致計(jì)算不穩(wěn)定,需引入人工耗散或者濾波.
為了能夠捕捉激波,基于Lele[17]提出的網(wǎng)格中心型緊致格式,結(jié)合WENO 格式的激波捕捉思想,Deng 等[18]設(shè)計(jì)了五階精度的加權(quán)緊致非線性格式(weighted compact nonlinear scheme,WCNS).WCNS 格式具有與WENO-JS 格式相當(dāng)?shù)募げú蹲侥芰洼^小的耗散,以及更靈活的通量選擇.Zhang等[19]將WCNS 格式拓展到了更高階精度,并且在構(gòu)造網(wǎng)格中心處的數(shù)值通量時(shí),是直接對(duì)通量進(jìn)行插值,而不是對(duì)守恒變量或原始變量進(jìn)行插值后再構(gòu)造通量.
CNCS 格式和CCCS 格式在構(gòu)造過程中,都只用了設(shè)計(jì)模板上的部分信息,格式精度和分辨率不能達(dá)到最優(yōu).Liu 等[20]設(shè)計(jì)了一類新型線性中心緊致格式(cemtral compact schemes with spectral-like resolution,CCSSR),該格式在求解網(wǎng)格結(jié)點(diǎn)上的函數(shù)導(dǎo)數(shù)時(shí),將網(wǎng)格結(jié)點(diǎn)上的函數(shù)值和網(wǎng)格中心處的函數(shù)值全部引入,精度最高可達(dá)14 階,并且接近譜方法的分辨率.隨后,劉旭亮等[21]以六階精度的線性中心緊致格式(CCSSR-L-6) 為基礎(chǔ),通過迎風(fēng)/對(duì)稱混合型加權(quán)非線性插值(WENO 插值)計(jì)算網(wǎng)格中心處的數(shù)值通量,構(gòu)造了能夠捕捉強(qiáng)激波的迎風(fēng)/對(duì)稱混合型加權(quán)非線性緊致格式,記為CCSSRHW-6 格式.
為了能夠捕捉強(qiáng)激波,增強(qiáng)格式在計(jì)算時(shí)的穩(wěn)定性,CCSSR-HW-6 格式在其構(gòu)造網(wǎng)格中心處數(shù)值通量的過程中包含兩級(jí)加權(quán): 一是五點(diǎn)迎風(fēng)模板插值和六點(diǎn)對(duì)稱模板插值的加權(quán);二是內(nèi)部子模板插值的加權(quán).兩級(jí)加權(quán)的權(quán)系數(shù)都不是固定值,而是非線性的函數(shù).兩級(jí)加權(quán)雖然增強(qiáng)了格式的耗散自適應(yīng)能力和激波捕捉能力,但也會(huì)使格式形式變得復(fù)雜,非線性效應(yīng)增強(qiáng).數(shù)值格式的非線性效應(yīng)會(huì)導(dǎo)致流場(chǎng)中產(chǎn)生非物理的寄生波[22-24],嚴(yán)重影響氣動(dòng)噪聲問題的高保真計(jì)算.在大多數(shù)激波噪聲問題中,激波的強(qiáng)度并不是很高,CCSSR-HW-6 格式的超強(qiáng)激波捕捉能力并不總是必需的,這就為一定程度上放寬激波捕捉能力的要求,進(jìn)而降低格式的非線性效應(yīng)和提高格式的分辨能力留出了空間.
本文通過優(yōu)化CCSSR-HW-6 格式中第一級(jí)加權(quán)的自由權(quán)系數(shù),使其取能夠讓格式空間分辨能力達(dá)到最優(yōu)的固定值,從而簡(jiǎn)化格式的構(gòu)造,減小格式的耗散誤差和非線性效應(yīng),使其更適合于激波噪聲問題的數(shù)值模擬.
現(xiàn)有的CCSSR-HW-6 格式是在具有類譜分辨率的六階精度線性中心緊致格式(CCSSR-L-6)的基礎(chǔ)上,結(jié)合迎風(fēng)/對(duì)稱混合型WENO 插值技術(shù)發(fā)展而來.圖1 是格式構(gòu)造模板.具體構(gòu)造過程[21]如下.
圖1 迎風(fēng)/對(duì)稱混合型加權(quán)非線性緊致格式的構(gòu)造模板Fig.1 Candidate stencils of hybrid upwind/symmetric weighted nonlinear compact scheme
六階精度的線性中心緊致格式(CCSSR-L-6)的表達(dá)式為
在內(nèi)部點(diǎn)上,采用隱式格式,系數(shù)取值為
在邊界點(diǎn)上,采用顯式格式,系數(shù)取值為
在CCSSR-L-6 格式中,網(wǎng)格中心上的函數(shù)值或數(shù)值通量是未知的.為了能夠捕捉激波,基于網(wǎng)格結(jié)點(diǎn)上的函數(shù)值,采用迎風(fēng)/對(duì)稱混合型WENO 插值得到網(wǎng)格中心上的函數(shù)值.迎風(fēng)/對(duì)稱混合型WENO插值技術(shù)包含兩級(jí)加權(quán): 一是五點(diǎn)迎風(fēng)模板插值和六點(diǎn)對(duì)稱模板插值的加權(quán);二是內(nèi)部四個(gè)子模板插值的加權(quán).
(1) 第一級(jí)加權(quán)是五點(diǎn)迎風(fēng)模板插值和六點(diǎn)對(duì)稱模板插值的加權(quán),具體如下
(2) 第二級(jí)加權(quán)是四個(gè)內(nèi)部子模板上插值的加權(quán),具體如下
其中,dr是子模板上的第二級(jí)權(quán)系數(shù).
子模板上的3 階精度線性插值多項(xiàng)式為
聯(lián)立式(4)和式(7)可建立兩級(jí)加權(quán)的權(quán)系數(shù)之間的關(guān)系
以上內(nèi)容即為CCSSR-HW-6 格式的線性部分.
(3) 非線性權(quán)的構(gòu)造,具體如下:
第一級(jí)權(quán)系數(shù) σ 定義如下[21]
式中,rc是常數(shù)閾值,其值越大,數(shù)值格式的耗散越大.CCSSR-HW-6 格式將其取值為rc=0.5.
相應(yīng)的非線性權(quán)[21]為
式中,ε是小值正數(shù),作用是避免分母為零,具體取值為ε=1.0×10-40.引入?yún)?shù)C的目的是增大最優(yōu)線性權(quán)的貢獻(xiàn),其值的增大或減小僅會(huì)改變數(shù)值格式的耗散誤差.參數(shù)C的值越大,數(shù)值格式的耗散誤差越小.在處理含有激波的問題時(shí),過大的C值可能會(huì)導(dǎo)致離散格式的數(shù)值不穩(wěn)定,CCSSR-HW-6 格式將參數(shù)C取值為C=5.
內(nèi)部四個(gè)子模板上的光滑指示因子[21]為
參考光滑指示因子[21]定義為
此外,非線性權(quán)會(huì)導(dǎo)致格式在特定的極值點(diǎn)處損失精度.為解決CCSSR-HW-6 格式在極值點(diǎn)處可能存在的精度降階問題,CCSSR-HW-6 格式還引入了映射函數(shù)[10]
使得非線性權(quán) ωr在臨界點(diǎn)處盡量接近線性權(quán)dr,從而還原格式的近似精度.
CCSSR-HW-6 格式的兩級(jí)加權(quán)確實(shí)增強(qiáng)了格式的耗散自適應(yīng)能力和激波捕捉能力,但也使格式的形式變得復(fù)雜,非線性效應(yīng)增強(qiáng).對(duì)于大多數(shù)激波噪聲問題,激波的強(qiáng)度并不是很高.因此,可以通過固定CCSSR-HW-6 格式中第一級(jí)加權(quán)的權(quán)系數(shù)取值,簡(jiǎn)化格式構(gòu)造,降低格式的耗散誤差,減弱格式的非線性效應(yīng),提高格式的分辨能力.
本文以CCSSR-HW-6 格式線性部分的修正波數(shù)為基礎(chǔ),采用修正波數(shù)的誤差積分函數(shù)[13,25-26]為空間分辨能力優(yōu)化目標(biāo)函數(shù),通過優(yōu)化,固定第一級(jí)加權(quán)權(quán)系數(shù)的取值,使格式的空間分辨能力達(dá)到最優(yōu).格式優(yōu)化過程具體如下.
第一步,推導(dǎo)CCSSR-HW-6 內(nèi)點(diǎn)格式和邊界格式線性部分的無量綱修正波數(shù).
CCSSR-HW-6 格式線性部分無量綱修正波數(shù)的計(jì)算方法如下:
考慮純粹的諧波函數(shù)
這里,x是位置,k是波數(shù).
諧波函數(shù)的解析導(dǎo)數(shù)為
定義n為任意整數(shù),Δx是網(wǎng)格間距,則有
一般有限差分方法的近似導(dǎo)數(shù)可表示為
其中,系數(shù)an是由差分格式中的系數(shù)決定.
定義無量綱波數(shù)為 κ=kΔx,結(jié)合式(16)~式(18)可得無量綱修正波數(shù)的一般表達(dá)式
利用上述方法可以得到CCSSR-HW-6 格式線性部分的無量綱修正波數(shù).
對(duì)于CCSSR-HW-6 內(nèi)點(diǎn)格式,無量綱修正波數(shù)表達(dá)式為
式(20)中的系數(shù)如下
對(duì)于CCSSR-HW-6 邊界格式,無量綱修正波數(shù)表達(dá)式為
式中的系數(shù)如下
第二步,確定空間分辨率優(yōu)化目標(biāo)函數(shù).
數(shù)值格式修正波數(shù)的實(shí)部表征了色散誤差(相位誤差),虛部則表征了耗散誤差(幅值誤差).為了對(duì)CCSSR-HW-6 格式進(jìn)行優(yōu)化,確定空間分辨率優(yōu)化目標(biāo)函數(shù)為修正波數(shù)的誤差積分函數(shù).修正波數(shù)的誤差積分函數(shù)綜合考慮了數(shù)值格式的色散誤差和耗散誤差.
無量綱修正波數(shù)的誤差積分函數(shù)[13,25-26]表達(dá)式為
其中,χ,ν ,ξ 和 η 是優(yōu)化控制參數(shù).
第三步,設(shè)置恰當(dāng)?shù)膬?yōu)化控制參數(shù).
優(yōu)化控制參數(shù)在誤差積分函數(shù)中分別具有不同的作用[13]: 參數(shù)χ控制相位誤差和幅值誤差的權(quán)值;參數(shù) ξ 控制了為確保格式穩(wěn)定而增加耗散的正弦項(xiàng);參數(shù) ν 影響數(shù)值格式在低波數(shù)區(qū)的誤差;參數(shù) η 影響數(shù)值格式在高波數(shù)區(qū)的誤差.
在本文中,優(yōu)化控制參數(shù)的取值與文獻(xiàn)[13]相同,分別為第四步,執(zhí)行優(yōu)化程序,得到優(yōu)化結(jié)果.
選擇優(yōu)化算法,執(zhí)行優(yōu)化程序,確定第一級(jí)加權(quán)的權(quán)系數(shù)取值,使得修正波數(shù)的誤差積分函數(shù)取最小值.
第一級(jí)加權(quán)的權(quán)系數(shù)優(yōu)化結(jié)果:
內(nèi)點(diǎn)格式優(yōu)化系數(shù)
邊界格式優(yōu)化系數(shù)
經(jīng)過空間分辨能力優(yōu)化后的迎風(fēng)/對(duì)稱混合型加權(quán)非線性緊致格式稱之為加權(quán)優(yōu)化緊致格式(weighted optimization compact scheme,WOCS).
通過一維標(biāo)量方程算例對(duì)加權(quán)優(yōu)化緊致格式(WOCS 格式)的收斂精度進(jìn)行數(shù)值驗(yàn)證.時(shí)間推進(jìn)采用三步三階TVD Runge-Kutta 方法[4].
一維線性標(biāo)量方程
初始條件為
計(jì)算采用周期邊界條件,總時(shí)長(zhǎng)為t=1;計(jì)算網(wǎng)格由10 逐步加密為160;計(jì)算網(wǎng)格為10 時(shí),初始CFL 數(shù)為0.5,隨著網(wǎng)格的每一次加密,對(duì)于r階格式,CFL 數(shù)以因子遞減.
表1 列出了WOCS 格式求解該問題時(shí)的數(shù)值誤差和精度階數(shù).結(jié)果顯示,WOCS 格式的收斂精度高于5 階.
表1 加權(quán)優(yōu)化緊致格式的數(shù)值誤差和精度階數(shù)Table 1 Numerical errors and accuracy order of the weighted optimization compact scheme
數(shù)值格式截?cái)嗾`差的精度階數(shù)只能提供數(shù)值解向精確解的漸進(jìn)收斂速率信息,不能給出有限計(jì)算網(wǎng)格上的實(shí)際誤差[27];而差分格式的波傳播特性(譜特性)能給出計(jì)算網(wǎng)格所支持的Fourier 模態(tài)的演化[27].
對(duì)于波數(shù)為k的單一Fourier 模態(tài),非線性格式不僅會(huì)在相同的波數(shù)k上產(chǎn)生線性響應(yīng),還會(huì)在其他波數(shù)上產(chǎn)生非線性響應(yīng)[24,28].對(duì)于非線性格式不存在解析的譜關(guān)系式,將廣泛應(yīng)用于線性格式的Fourier 分析直接應(yīng)用到非線性格式的空間分辨能力分析是行不通的.針對(duì)非線性格式發(fā)展的近似色散關(guān)系(approximate dispersion relation,ADR)技術(shù)[27]實(shí)際上得到的只是非線性格式修正波數(shù)的線性響應(yīng)部分.基于空間導(dǎo)數(shù)的修正波數(shù)計(jì)算方法[24,28]可以得到數(shù)值格式的線性響應(yīng)和非線性響應(yīng),具體如下.
非線性格式對(duì)于波數(shù)k的修正波數(shù)為
圖2 是WOCS 格式和原CCSSR-HW-6 格式修正波數(shù)線性響應(yīng)部分的對(duì)比.線性響應(yīng)的實(shí)部和虛部分別表征色散特性和耗散特性.圖2(a)顯示,WOCS格式和原CCSSR-HW-6 格式具有相同的色散誤差;圖2(b)顯示,WOCS 格式在中高波數(shù)區(qū)的耗散誤差明顯小于原CCSSR-HW-6 格式.總的來說,通過優(yōu)化使CCSSR-HW-6 格式第一級(jí)加權(quán)的權(quán)系數(shù) σ 取固定值,提高了格式在中高波數(shù)區(qū)的空間分辨能力.
圖2 修正波數(shù)的線性響應(yīng)部分Fig.2 The linear response part of the modified wavenumber
圖3 是WOCS 格式和原CCSSR-HW-6 格式的修正波數(shù)非線性響應(yīng)部分的Fourier 系數(shù)對(duì)比.計(jì)算網(wǎng)格點(diǎn)數(shù)為N=64,在該網(wǎng)格所支持所有 Fourier 模態(tài)中,模態(tài)5、模態(tài)10 和模態(tài)15 是屬于中低波數(shù)的模態(tài)(模態(tài)k表示該模態(tài)的波數(shù)是k).對(duì)于任意模態(tài)k,非線性格式除了在波數(shù)k上會(huì)產(chǎn)生線性響應(yīng),還會(huì)在其他波數(shù) (3+2n)kmodN(N是整數(shù))上產(chǎn)生非線性響應(yīng)(非物理的寄生波)[24].非線性響應(yīng)的實(shí)數(shù)部分代表寄生波的波數(shù)響應(yīng),虛數(shù)部分代表寄生波的幅值響應(yīng).圖3 顯示,在中低波數(shù)區(qū)間,WOCS格式非線性響應(yīng)的實(shí)數(shù)部分與原CCSSR-HW-6 格式相當(dāng),而非線性響應(yīng)的虛數(shù)部分則明顯弱于原CCSSR-HW-6 格式.綜合來看,與原CCSSR-HW-6 格式相比,WOCS 格式降低了中低波數(shù)區(qū)的非線性響應(yīng).
圖3 修正波數(shù)的非線性響應(yīng)部分Fig.3 The nonlinear response part of the modified wavenumber
本節(jié)針對(duì)加權(quán)優(yōu)化緊致格式(WOCS 格式)在激波噪聲問題計(jì)算中的適用性進(jìn)行測(cè)試,包括激波捕捉能力、分辨能力和非線性效應(yīng).數(shù)值實(shí)驗(yàn)包括一維Euler 方程算例、二維Euler 方程算例和二維Navier-Stokes 方程算例.無黏數(shù)值通量計(jì)算采用HLL 型Riemann 求解器[29]或Lax-Fredrichs 通量分裂方法[4],時(shí)間推進(jìn)采用三步三階TVD Runge-Kutta 方法[4].
(1)激波-聲波相互作用
激波-聲波相互作用[30]是氣動(dòng)聲學(xué)中的重要基準(zhǔn)問題之一.計(jì)算條件如下.
計(jì)算域?yàn)閇0,1],計(jì)算網(wǎng)格數(shù)為N=401,CFL 數(shù)為0.5,計(jì)算時(shí)長(zhǎng)為t=30Tλ;初始時(shí)刻,激波位置為xs=0.5,激波馬赫數(shù)為Ms=2,氣體從左向右流動(dòng),激波左右兩側(cè)的流動(dòng)變量滿足激波關(guān)系式.計(jì)算域右邊界設(shè)置為無反射邊界條件,在左邊界處添加聲擾動(dòng)
式中參數(shù)取值如下
圖4 是激波-聲波相互作用的數(shù)值解及其與參考解的比較.由于解析解是未知的,這里選取五階精度WCNS-E-5 格式[18]在計(jì)算網(wǎng)格數(shù)為N=1001 時(shí)求解該問題得到的數(shù)值解作為參考解.圖4(b) 和圖4(c) 是圖4(a) 的局部放大圖.圖中的縱坐標(biāo)δP(t)=P(x,t)-P(x,0)表示聲擾動(dòng).聲波穿過激波后其幅值會(huì)被放大.結(jié)果顯示: 在分辨穿過激波后的聲波時(shí),WOCS 格式計(jì)算得到的聲波幅值與參考解更為吻合,其分辨能力優(yōu)于原CCSSR-HW-6 格式.
圖4 激波-聲波相互作用問題的數(shù)值結(jié)果Fig.4 Numerical results of shock-sound interaction problem
(2)激波-密度波相互作用
激波-密度波相互作用,即Titarev-Toro 問題[31],也是驗(yàn)證數(shù)值格式激波捕捉能力和流場(chǎng)細(xì)微結(jié)構(gòu)分辨能力的標(biāo)準(zhǔn)算例,在Titarev-Toro 問題中,向右運(yùn)動(dòng)的激波(馬赫數(shù)為1.1)與正弦波形式的高頻密度擾動(dòng)相互作用,產(chǎn)生的流場(chǎng)既包含光滑解,也包含間斷.該問題可以看作是激波/湍流相互作用的簡(jiǎn)化模型.初始條件如下
計(jì)算條件: 計(jì)算網(wǎng)格數(shù)為N=1001,CFL 數(shù)為0.5,計(jì)算時(shí)間為t=5.0;邊界條件設(shè)置為計(jì)算域左、右邊界處的守恒變量不隨時(shí)間變化,這在本算例所感興趣的計(jì)算時(shí)長(zhǎng)內(nèi)是適宜的.
圖5 是Titarev-Toro 問題的數(shù)值解及其與參考解的比較.由于解析解是未知的,這里選取五階精度WCNS-E-5 格式[18]在計(jì)算網(wǎng)格數(shù)為N=10 001 時(shí)求解該問題得到的數(shù)值解作為參考解.圖5(b)和圖5(c)是圖5(a)的局部放大圖.結(jié)果顯示: WOCS 格式數(shù)值解的幅值與參考解更為接近,其對(duì)高頻波的分辨能力遠(yuǎn)優(yōu)于原CCSSR-HW-6 格式.這是因?yàn)橄噍^于原CCSSR-HW-6 格式,WOCS 格式顯著降低了高波數(shù)區(qū)的耗散誤差,如圖2(b)所示.
圖5 Titarev-Toro 問題的數(shù)值結(jié)果Fig.5 Numerical results of Titarev-Toro problem
選擇激波-渦量波相互作用[32-33]作為二維Euler方程測(cè)試算例,該問題是激波-湍流相互作用的二維簡(jiǎn)化模型.在此算例中,馬赫數(shù)為Ms=8 的運(yùn)動(dòng)激波與散度為零的傾斜渦量場(chǎng)發(fā)生相互作用.計(jì)算域范圍: (x,y)∈[-1.5,1.5]×[-1,1],均勻計(jì)算網(wǎng)格,網(wǎng)格間距為 Δx=Δy=0.01 ;計(jì)算時(shí)長(zhǎng)t=0.2.
初始時(shí)刻,激波位于x=-1 處,傳播速度為
這里的cs1是波前聲速(下標(biāo)1 表示波前狀態(tài),下標(biāo)2 表示波后狀態(tài)).
波前(x≥-1)初始狀態(tài)為
波后(x≤-1)初始狀態(tài)為
邊界條件設(shè)置: 在y方向,上、下兩端采用周期邊界條件;在x方向,左端為超聲速入流條件,右端為無反射邊界條件.
圖6 是采用WOCS 格式計(jì)算激波-渦量波相互作用得到的t=0.2 時(shí)刻的渦量空間分布.圖7 是t=0.2時(shí)刻直線y=0 上的渦量分布及其與參考解的比較.參考解是在計(jì)算網(wǎng)格數(shù)為N=961 × 641 條件下采用CCSSR-HW-6 格式計(jì)算該問題時(shí)得到的數(shù)值解.結(jié)果顯示: WOCS 格式的數(shù)值解與參考解要更為接近,分辨能力優(yōu)于原CCSSR-HW-6 格式.
圖6 激波-渦量波相互作用在t =0.2 時(shí)刻的渦量空間分布Fig.6 The spatial distribution of vorticity in the shock-vorticity wave interaction at t =0.2
圖7 激波-渦量波相互作用在t =0.2 時(shí)刻沿直線y=0 的渦量分布Fig.7 The vorticity distribution along the line y=0 in the shockvorticity wave interaction at t =0.2
二維Navier-Stokes 方程測(cè)試算例選取的是激波-強(qiáng)旋渦相互作用問題[6,34].激波-旋渦相互作用是重要的激波噪聲理論研究模型,文獻(xiàn)[6,34]對(duì)此問題進(jìn)行了深入細(xì)致的研究.此外,該問題也被廣泛用作高精度數(shù)值方法的驗(yàn)證算例[35].下面簡(jiǎn)要介紹該問題計(jì)算細(xì)節(jié).
計(jì)算域范圍: (x,y)∈[-20,8]×[-12,12],均勻計(jì)算網(wǎng)格,網(wǎng)格點(diǎn)數(shù)為N=701 × 601,計(jì)算總時(shí)長(zhǎng)t=10.普朗特?cái)?shù)Pr=0.75,參考溫度T∞=300 K .雷諾數(shù)Re=800,具體定義為Re=ρ∞a∞R/μ∞,其中,R是渦核半徑,ρ∞,a∞和 μ∞分別是激波前流動(dòng)的平均密度、聲速和黏性系數(shù).初始條件設(shè)置如下.
激波馬赫數(shù)為Ms=1.2,固定在xs=0 的位置.激波前(x≥0)的狀態(tài)參數(shù)為
激波后(x<0)的狀態(tài)參數(shù)為
初始時(shí)刻,旋渦位于激波波前,其從右向左穿越激波的過程中與之發(fā)生相互作用,產(chǎn)生并向外輻射聲波.初始渦心位置是 (xv,yv)=(2,0),旋渦強(qiáng)度為Mv=1.0.旋渦的速度、密度和壓力分布如下
邊界條件設(shè)置: 在y方向,上、下兩端采用周期邊界條件;在x方向,右端為超聲速入流條件,左端為無反射邊界條件.
對(duì)于氣動(dòng)聲學(xué)問題,脹量(速度的散度,?·V)和數(shù)值陰影(?2ρ)常用來表征聲場(chǎng)和流動(dòng)結(jié)構(gòu).圖8是采用WOCS 格式計(jì)算激波-強(qiáng)旋渦相互作用得到的t=10 時(shí)刻的脹量空間分布,從中可以看到向外輻射的聲波.圖9 是t=10 時(shí)刻脹量和數(shù)值陰影沿徑向(θ=π/4,θ 的定義見圖8)的分布及其與參考解的比較.參考解是在計(jì)算網(wǎng)格數(shù)為N=1401 × 1201條件下采用CCSSR-HW-6 格式計(jì)算該問題時(shí)得到的數(shù)值解.對(duì)于脹量,WOCS 和CCSSR-HW-6 兩種格式的計(jì)算結(jié)果一致;而對(duì)于數(shù)值陰影,WOCS 格式的計(jì)算結(jié)果優(yōu)于原CCSSR-HW-6 格式的計(jì)算結(jié)果,其與參考解要更為接近.
圖8 激波-強(qiáng)旋渦相互作用在t=10 時(shí)刻的脹量場(chǎng)Fig.8 The dilatation field of shock-strong vortex interaction at t=10
圖9 激波-強(qiáng)旋渦相互作用在t =10 時(shí)刻的脹量和數(shù)值陰影沿徑向的分布Fig.9 The radial distribution of dilatation and numerical shadowgraph in shock-strong vortex interaction at t=10
圖10 是分別采用WOCS 格式和CCSSR-HW-6 格式計(jì)算激波-強(qiáng)旋渦相互作用得到的t=10 時(shí)刻的數(shù)值陰影空間分布.結(jié)果顯示: 在CCSSR-HW-6格式計(jì)算得到的數(shù)值陰影中存在非常明顯的非物理振蕩,而WOCS 格式計(jì)算得到的數(shù)值陰影則更為“干凈”,上述非物理振蕩被顯著地減弱.文獻(xiàn)[22-24]指出激波捕捉格式在計(jì)算時(shí)的非線性執(zhí)行會(huì)導(dǎo)致非線性效應(yīng),使得數(shù)值解中產(chǎn)生此類非物理波動(dòng).前文3.2 節(jié)中的格式非線性響應(yīng)分析已經(jīng)表明,與原CCSSR-HW-6 格式相比,經(jīng)優(yōu)化后,WOCS 格式的非線性響應(yīng)顯著降低,由此可以解釋這兩種格式計(jì)算結(jié)果之間的差異.
圖10 激波-強(qiáng)旋渦相互作用在t=10 時(shí)刻的數(shù)值陰影Fig.10 Numerical shadowgraph of shock-strong vortex interaction at t=10
基于對(duì)稱模板的迎風(fēng)/對(duì)稱混合型加權(quán)非線性緊致格式在其構(gòu)造過程中包含兩級(jí)加權(quán),并且兩級(jí)加權(quán)的權(quán)系數(shù)都是非線性函數(shù),格式非線性效應(yīng)較強(qiáng).為了減弱格式的非線性效應(yīng),改進(jìn)格式的譜特性,本文基于六階精度的迎風(fēng)/對(duì)稱混合型加權(quán)非線性緊致(CCSSR-HW-6)格式,以修正波數(shù)的誤差積分函數(shù)為優(yōu)化目標(biāo)函數(shù),通過優(yōu)化第一級(jí)加權(quán)的權(quán)系數(shù),使其取分辨能力最優(yōu)的固定值,建立了加權(quán)優(yōu)化緊致格式,記為WOCS 格式.針對(duì)WOCS 格式的精度、譜特性和激波噪聲問題計(jì)算的適用性進(jìn)行了數(shù)值驗(yàn)證.
精度驗(yàn)證表明WOCS 格式的精度高于5 階.譜特性分析包括了線性響應(yīng)和非線性響應(yīng).線性響應(yīng)分析表明WOCS 格式在中高波數(shù)區(qū)的耗散特性明顯優(yōu)于CCSSR-HW-6 格式.非線性響應(yīng)分析表明:與CCSSR-HW-6 格式相比,WOCS 格式降低了中低波數(shù)區(qū)的非線性響應(yīng).
激波噪聲計(jì)算適用性的測(cè)試算例包括了激波-聲波相互作用、激波-密度波相互作用、激波-渦量波相互作用以及激波-強(qiáng)旋渦相互作用.數(shù)值實(shí)驗(yàn)結(jié)果表明: 與原CCSSR-HW-6 格式相比,WOCS 格式不僅降低了耗散誤差,提高了對(duì)高頻波的分辨能力,而且顯著地減弱了格式的非線性效應(yīng)和由此導(dǎo)致的數(shù)值解的非物理振蕩,較為適合于激波噪聲問題的數(shù)值模擬.