李彥蘇， 張坤， 何承軍， 閻超

(1. 北京機(jī)電工程研究所, 北京 100074； 2. 北京航空航天大學(xué)航空科學(xué)與工程學(xué)院, 北京 100083)

準(zhǔn)確捕捉特點(diǎn)各異的流動(dòng)結(jié)構(gòu)、獲得高精度流場(chǎng)數(shù)據(jù)是利用數(shù)值仿真開展湍流問題研究的前提。在高速可壓縮湍流流動(dòng)中，間斷(如激波)和多尺度渦結(jié)構(gòu)共存[1]，要求數(shù)值格式既能夠穩(wěn)定捕捉激波，又具有較高的分辨率。而邊界層等局部低速區(qū)的存在對(duì)數(shù)值方法的精度、頻譜特性及穩(wěn)定性也提出了嚴(yán)苛的要求。最近研究發(fā)現(xiàn)[2]，當(dāng)馬赫數(shù)趨近0時(shí)，可壓縮求解器的空間離散方法(主要是通量格式)往往無法正確逼近原始方程。因此，在求解高速流動(dòng)中的局部低速區(qū)域時(shí)，可能出現(xiàn)耗散大、收斂慢等問題。

針對(duì)這一問題，學(xué)者提出了預(yù)處理矩陣技術(shù)[3]，以改變?cè)刂品匠痰臄?shù)學(xué)性質(zhì)，再用現(xiàn)有通量格式進(jìn)行計(jì)算。數(shù)值實(shí)驗(yàn)表明，該方法一定程度上改善了低速流動(dòng)的求解精度。但該方法通常需要設(shè)置經(jīng)驗(yàn)參數(shù)，參數(shù)的選擇對(duì)計(jì)算效率及精度影響較大，尤其在流場(chǎng)中流速跨度較大的情況下難以保證魯棒性。另一部分學(xué)者則提出了對(duì)現(xiàn)有通量格式(如Roe格式、AUSM系列格式)的修正技術(shù)，較為常見的有：修正了AUSM+格式中壓力通量的AUSM+up格式[4]，無需調(diào)節(jié)人工參數(shù)的SLAU、SLAU2格式[5]，以及通過對(duì)Roe格式理論分析獲得的LMRoe格式[6]、F-Roe格式[7]、Roe-m格式[8]等。這些研究對(duì)不同馬赫數(shù)下普通通量格式和全速域格式的計(jì)算結(jié)果進(jìn)行了深入對(duì)比和分析，馬赫數(shù)的模擬范圍可低至10-3[6]。

值得注意的是，在研究低速問題計(jì)算時(shí)，為了簡(jiǎn)化模型，學(xué)者們大多圍繞Euler方程開展研究，并使用階次較低的計(jì)算精度以突顯低速修正對(duì)計(jì)算精度的影響。但在開展湍流模擬時(shí)，使用的計(jì)算方法精度較高，網(wǎng)格量較大，這些因素將與求解器低速修正耦合，對(duì)湍流模擬產(chǎn)生綜合影響。

因此，本文重點(diǎn)研究有無低速修正的可壓縮求解器對(duì)湍流模擬的影響，旨在定量分析不同精度、網(wǎng)格量下，低速修正對(duì)模擬精度的改進(jìn)效果，評(píng)估低速修正的可壓縮求解器的使用范圍，為今后計(jì)算格式的構(gòu)造和應(yīng)用提供參考。

1 控制方程及計(jì)算方法

1.1 三維可壓縮Navier-Stokes方程及離散方法

本文的控制方程為三維可壓縮Navier-Stokes方程組[9]，即

(1)

式中：Q為守恒變量；F、G、H分別代表x、y和z三方向的通量，下標(biāo)c表示對(duì)流通量，下標(biāo)v表示黏性通量。

將式(1)寫成半離散形式為

(2)

式中：下標(biāo)i、j、k分別代表x、y和z三方向的離散節(jié)點(diǎn)，i±1/2、j±1/2、k±1/2分別表示對(duì)應(yīng)節(jié)點(diǎn)左右半節(jié)點(diǎn)；上標(biāo)“～”表示半節(jié)點(diǎn)處x、y和z三方向的無黏通量，通常使用重構(gòu)格式和通量格式聯(lián)合計(jì)算獲得；上標(biāo)v表示半節(jié)點(diǎn)處x、y和z三方向的黏性通量，使用中心差分格式求解。

空間離散變量求解完成后，方程[9]變成常微分方程，使用三階TVD型Runge-Kutta方法[10]進(jìn)行時(shí)間推進(jìn)，獲得下一時(shí)刻流場(chǎng)變量Q的分布。

通過對(duì)上述流場(chǎng)計(jì)算過程的分析可知，方程中空間項(xiàng)的模擬精度主要受通量格式和重構(gòu)格式的雙重影響。在低速問題中，2類格式共同影響計(jì)算結(jié)果。

1.2 通量格式

Roe格式[11]是Godunov類方法，是目前經(jīng)典的通量格式，通過求解線化Riemann問題獲得全場(chǎng)的數(shù)值近似解,其具有間斷分辨率高、穩(wěn)定性強(qiáng)、計(jì)算效率高等優(yōu)點(diǎn)，廣泛使用于可壓縮流動(dòng)求解中。以x方向?yàn)槔?，Roe格式可以寫成

(3)

低速情況下，原始Roe格式壓強(qiáng)與馬赫數(shù)的變化關(guān)系和控制方程不同，這將導(dǎo)致求解精度降低、收斂速度變慢等問題。針對(duì)這一問題，可引入當(dāng)?shù)伛R赫數(shù)[6]，即

(4)

式中：Uloc為當(dāng)?shù)亓魉?；aloc為當(dāng)?shù)芈曀佟?/p>

利用當(dāng)?shù)伛R赫數(shù)修正原始Roe格式左右界面的速度差ΔU，將其替換成min(Maloc,1)ΔU。通過這一修正，在馬赫數(shù)趨近于0時(shí)，格式的壓強(qiáng)與馬赫數(shù)的關(guān)系和原控制方程相同。修正后的格式即為L(zhǎng)MRoe格式[6]。

1.3 重構(gòu)格式

為了全面分析通量格式與不同階次、分辨率的重構(gòu)格式組合后的模擬效果，本文使用的重構(gòu)格式包括三/五/七階WENO格式和五階線性緊致格式。其中，WENO格式通過低階模板的加權(quán)組合達(dá)到高階精度，同時(shí)在間斷區(qū)域降低相應(yīng)模板的權(quán)重，從而穩(wěn)定捕捉激波，常用于對(duì)計(jì)算精度要求較高的高速流動(dòng)模擬；五階線性緊致格式能夠在相同階數(shù)的情況下獲得更高的分辨率。

2m-1階WENO格式的表達(dá)式可寫為

(5)

(6)

非線性權(quán)系數(shù)ω的表達(dá)式為

(7)

式中：Ck為理想權(quán)系數(shù)，對(duì)于五階WENO格式，其值為C0=0.1，C1=0.6，C2=0.3；ISk為衡量當(dāng)?shù)啬０骞饣缘墓饣蜃?，詳?xì)表達(dá)式見文獻(xiàn)[12]。

其他階數(shù)的WENO格式可以通過類似方法獲得。

五階線性緊致格式的表達(dá)式為

(8)

從式(8)可見，五階線性緊致格式不需要計(jì)算非線性系數(shù)，計(jì)算量較小。但線性格式對(duì)數(shù)據(jù)光滑性的要求非常高，在流場(chǎng)具有間斷的可壓縮流動(dòng)數(shù)值模擬中容易產(chǎn)生非物理振蕩，甚至發(fā)散。因此，單純的線性格式難以在可壓縮復(fù)雜流動(dòng)數(shù)值模擬中直接應(yīng)用。

為了表述方便，稱三/五/七階WENO格式為WENO3、WENO5和WENO7，稱五階線性緊致格式為compact5。

傅里葉分析可以衡量格式的分辨率[13-14]，即表征在網(wǎng)格量不足情況下格式的計(jì)算精度(見圖1)。圖中:k為波數(shù)，Re(k′)和Im(k′)分別為色散和耗散特性，橫、縱軸均使用π歸一化?？梢姡S著階數(shù)升高，WENO格式的分辨率提高。而五階緊致格式的分辨率遠(yuǎn)優(yōu)于七階WENO格式，表明了緊致格式在分辨率具有優(yōu)勢(shì)。

圖1 不同重構(gòu)格式的分辨率特性曲線Fig.1 Resolution properties of different reconstruction schemes

2 泰勒-格林渦算例

2.1 算例設(shè)置

為了定量分析不同格式的計(jì)算性能，本文選擇了經(jīng)典的泰勒-格林渦算例。算例的流動(dòng)演化包括了“層流—轉(zhuǎn)捩—湍流—衰減”這幾個(gè)湍流發(fā)展的典型階段[15]，能夠反映數(shù)值方法對(duì)湍流不同發(fā)展階段的模擬能力。同時(shí)，該算例的初始流場(chǎng)有解析表達(dá)式，可以避免“啟動(dòng)問題”(set-up problem)。此外，該算例本質(zhì)上是不可壓縮流動(dòng)，實(shí)際模擬時(shí)馬赫數(shù)非常小，低速通量格式在模擬中能夠體現(xiàn)其作用。

算例的初始流場(chǎng)為

(9)

下文的分析中主要用到了體平均動(dòng)能Ek及其耗散率εEk這2個(gè)平均量。體平均動(dòng)能的計(jì)算式為

(10)

式中：u為速度矢量；Ω表示積分域。其對(duì)應(yīng)的動(dòng)能耗散率的計(jì)算式為

(11)

Ek和εEk能夠反映流動(dòng)發(fā)展過程中流場(chǎng)動(dòng)能變化的情況，從而衡量模擬結(jié)果的精度。

2.2 計(jì)算結(jié)果與討論

在網(wǎng)格間距為2π/64、2π/128和2π/256的3套網(wǎng)格下開展算例計(jì)算，以分析不同網(wǎng)格量下不同格式對(duì)計(jì)算結(jié)果的影響。在間距為2π/256的密網(wǎng)格上使用compact5獲得的結(jié)果與文獻(xiàn)[16]使用DRP方法在5123自由度上的DNS結(jié)果精度相當(dāng)，作為參考解進(jìn)行計(jì)算結(jié)果誤差的定性和定量分析。在粗網(wǎng)格(間距2π/64)上使用3種階數(shù)WENO格式和compact5與Roe和LMRoe 2種通量格式分別組合，重點(diǎn)考察粗網(wǎng)格量下重構(gòu)格式精度不同時(shí)，低速修正通量格式的使用對(duì)計(jì)算精度的影響程度。在中等網(wǎng)格(間距2π/128)上使用WENO5、compact5與Roe和LMRoe 2種通量格式兩兩組合，進(jìn)一步研究網(wǎng)格量變化后通量格式的低速修正對(duì)計(jì)算結(jié)果的影響。

圖2和圖3為粗網(wǎng)格上不同重構(gòu)格式下獲得的體平均動(dòng)能及其耗散率隨時(shí)間變化曲線，其中文獻(xiàn)[16]的DNS結(jié)果(圖中“reference-DRP5123”)和compact5在密網(wǎng)格的結(jié)果(圖中“256-Roe-compact5”)為參考解。總體上看，在網(wǎng)格量相同的情況下，重構(gòu)格式的精度越高，分辨率越高，計(jì)算結(jié)果越接近參考解。相同重構(gòu)格式下，LMRoe格式的計(jì)算結(jié)果更接近參考解，可見在粗網(wǎng)格下，低速修正起到了提高計(jì)算精度的作用。

圖4和圖5為中等網(wǎng)格下WENO5和compact5兩種重構(gòu)格式、Roe和LMRoe兩種通量格式兩兩組合的計(jì)算結(jié)果。與粗網(wǎng)格結(jié)果相比，密網(wǎng)格下的計(jì)算結(jié)果精度更高，且2種通量格式計(jì)算結(jié)果間的差別變小。

為了定量比較不同格式的計(jì)算結(jié)果差異，以網(wǎng)格間距2π/256的密網(wǎng)格compact5計(jì)算獲得的參考解為基準(zhǔn)，得到不同計(jì)算結(jié)果與基準(zhǔn)解間的相對(duì)誤差，公式為

(12)

式中：g(t)表示某格式計(jì)算結(jié)果中某個(gè)變量隨時(shí)間的變化；gc(t)表示密網(wǎng)格compact5計(jì)算結(jié)果中對(duì)應(yīng)變量隨時(shí)間的變化。

由此獲得2套網(wǎng)格下格式的誤差并制成柱狀圖，如圖6和圖7所示。表1和表2則給出2種通量格式誤差的相對(duì)量，即

(13)

可以看出，對(duì)于WENO系列格式，隨著格式階數(shù)的提高，計(jì)算精度提高明顯；使用LMRoe格式的精度略高于相同網(wǎng)格量時(shí)的Roe格式；加密網(wǎng)格后計(jì)算精度顯著提高。而當(dāng)使用compact5作為重構(gòu)格式時(shí)，相同網(wǎng)格量下2種通量格式的計(jì)算差異不明顯，尤其是體平均動(dòng)能，LMRoe格式的計(jì)算誤差甚至大于Roe格式。Compact5和WENO格式的計(jì)算結(jié)果比較則可看出，compact5格式的計(jì)算精度高于所有WENO格式；特別是相同精度階數(shù)情況下，compact5在粗網(wǎng)格(2π/64)上的計(jì)算結(jié)果與WENO5在密網(wǎng)格(2π/128)上的最佳結(jié)果(LMRoe的計(jì)算結(jié)果)精度相當(dāng)。

圖2 體平均動(dòng)能隨時(shí)間變化曲線(網(wǎng)格間距2π/64)Fig.2 Variation of volume-averaged kinetic energy with time under grid space 2π/64

圖3 動(dòng)能耗散率隨時(shí)間變化曲線(網(wǎng)格間距2π/64)Fig.3 Variation of energy dissipation rate with time under grid space 2π/64

圖4 體平均動(dòng)能隨時(shí)間變化曲線(網(wǎng)格間距2π/64和2π/128)Fig.4 Variation of volume-averaged kinetic energy with time (grid space 2π/64 and 2π/128)

利用如下公式：

(14)

可以更加清晰地衡量2個(gè)通量格式的計(jì)算差別(下標(biāo)Roe和LMRoe表示計(jì)算使用的是原始Roe格式或LMRoe格式)。

圖5 動(dòng)能耗散率隨時(shí)間變化曲線(網(wǎng)格間距2π/64和2π/128)Fig.5 Variation of energy dissipation rate with time (grid space 2π/64 and 2π/128)

圖6 不同格式和網(wǎng)格量下體平均動(dòng)能誤差柱狀圖Fig.6 Histogram of volume-averaged kinetic energy error for different schemes with different grids

圖8和圖9為2套網(wǎng)格下原始Roe格式和LMRoe格式的計(jì)算差別χ′柱狀圖。可以看出，隨著網(wǎng)格加密，2種通量格式的計(jì)算結(jié)果間差別變小，表現(xiàn)出向精確解收斂的特征。

上述分析表明，通量格式、重構(gòu)格式對(duì)計(jì)算精度的影響是耦合的，通量格式、重構(gòu)格式、網(wǎng)格量等各部分?jǐn)?shù)值誤差的總和構(gòu)成了最終的計(jì)算誤差。在網(wǎng)格較粗且重構(gòu)格式精度較低的情況下， LMRoe格式能夠顯著提高計(jì)算精度。結(jié)合LMRoe格式修正原理，這是由于在較粗網(wǎng)格下，低馬赫數(shù)下LMRoe格式求解的壓強(qiáng)與馬赫數(shù)的對(duì)應(yīng)關(guān)系與Navier-Stokes方程一致，導(dǎo)致求解誤差明顯降低，整體計(jì)算結(jié)果顯著改善；而原始Roe格式低速時(shí)，壓強(qiáng)與馬赫數(shù)對(duì)應(yīng)關(guān)系與Navier-Stokes方程不一致的問題在粗網(wǎng)格下凸顯，對(duì)整體計(jì)算結(jié)果的影響較大。但網(wǎng)格加密、重構(gòu)格式精度提高后，這兩部分誤差的減小一定程度上彌補(bǔ)了原始Roe格式帶來的誤差，使得總體的計(jì)算精度有所提升。且對(duì)于網(wǎng)格收斂的格式，當(dāng)網(wǎng)格足夠密時(shí)，能夠給獲得收斂的結(jié)果，也就是當(dāng)網(wǎng)格足夠密時(shí)，原始Roe格式也能夠得到精度足夠高的結(jié)果。因此，當(dāng)使用密網(wǎng)格、高分辨率重構(gòu)格式時(shí)，原始Roe格式和LMRoe格式的計(jì)算差異很小。

圖7 不同格式和網(wǎng)格量下動(dòng)能耗散率誤差柱狀圖Fig.7 Histogram of energy dissipation rate error for different schemes with different grids

重構(gòu)格式體平均動(dòng)能動(dòng)能耗散率WENO30.360.57 WENO50.630.61WENO70.690.75compact51.041.00

表2 網(wǎng)格間距2π/128時(shí)不同通量格式結(jié)果的誤差比值Table 2 Ratio of two flux schemes’ result errors with grid space being 2π/128

圖8 不同網(wǎng)格量下Roe和LMRoe通量格式的計(jì)算差異(體平均動(dòng)能)Fig.8 Calculation difference of Roe and LMRoe flux schemes with different amounts of grid (volume-averaged kinetic energy)

圖9 不同網(wǎng)格量下Roe和LMRoe通量格式的計(jì)算差異(動(dòng)能耗散率)Fig.9 Calculation difference of Roe and LMRoe flux schemes with different amounts of grid (energy dissipation rate)

3 結(jié) 論

通量格式和重構(gòu)格式的精度均會(huì)影響湍流流動(dòng)的模擬結(jié)果。本文利用泰勒-格林渦算例，研究了有無低速修正的通量格式對(duì)模擬結(jié)果的影響，結(jié)論如下：

1) 當(dāng)流場(chǎng)中存在低速流動(dòng)時(shí)，使用低速修正的通量格式能夠一定程度上改進(jìn)計(jì)算結(jié)果。

2) 通量格式對(duì)模擬結(jié)果的影響是與重構(gòu)格式耦合的。當(dāng)重構(gòu)格式的精度較低時(shí)，通量格式對(duì)結(jié)果的影響顯著，但當(dāng)重構(gòu)格式的精度足夠高后，通量格式對(duì)結(jié)果的影響不明顯。

3) 隨著網(wǎng)格加密，有無低速修正的通量格式計(jì)算結(jié)果都呈現(xiàn)出向精確解收斂的特征。使用較粗網(wǎng)格時(shí)，低速修正的通量格式對(duì)計(jì)算結(jié)果的改進(jìn)更明顯。

考慮到在實(shí)際工程應(yīng)用時(shí)，由于研究對(duì)象外形通常較復(fù)雜，需要使用魯棒性較高、耗散較大的重構(gòu)格式以保證計(jì)算穩(wěn)定性，且受到研究周期和計(jì)算條件的限制，網(wǎng)格量通常相對(duì)較小，使用低速修正的通量格式能夠提高計(jì)算精度。在開展湍流流動(dòng)機(jī)理研究時(shí)，對(duì)計(jì)算精度要求較高，通常使用高精度高分辨率重構(gòu)格式及較密的網(wǎng)格開展模擬，此時(shí)低速修正的通量格式對(duì)計(jì)算結(jié)果改進(jìn)有限，未經(jīng)過低速修正的原始通量格式也能夠獲得較好的計(jì)算結(jié)果。