低速修正的可壓縮求解器對湍流模擬精度的影響

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

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

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

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

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

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

1 控制方程及計算方法

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

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

(1)

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

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

(2)

式中：下標i、j、k分別代表x、y和z三方向的離散節點，i±1/2、j±1/2、k±1/2分別表示對應節點左右半節點；上標“～”表示半節點處x、y和z三方向的無黏通量，通常使用重構格式和通量格式聯合計算獲得；上標v表示半節點處x、y和z三方向的黏性通量，使用中心差分格式求解。

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

通過對上述流場計算過程的分析可知，方程中空間項的模擬精度主要受通量格式和重構格式的雙重影響。在低速問題中，2類格式共同影響計算結果。

1.2 通量格式

Roe格式[11]是Godunov類方法，是目前經典的通量格式，通過求解線化Riemann問題獲得全場的數值近似解,其具有間斷分辨率高、穩定性強、計算效率高等優點，廣泛使用于可壓縮流動求解中。以x方向為例，Roe格式可以寫成

(3)

低速情況下，原始Roe格式壓強與馬赫數的變化關系和控制方程不同，這將導致求解精度降低、收斂速度變慢等問題。針對這一問題，可引入當地馬赫數[6]，即

(4)

式中：Uloc為當地流速；aloc為當地聲速。

利用當地馬赫數修正原始Roe格式左右界面的速度差ΔU，將其替換成min(Maloc,1)ΔU。通過這一修正，在馬赫數趨近于0時，格式的壓強與馬赫數的關系和原控制方程相同。修正后的格式即為LMRoe格式[6]。

1.3 重構格式

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

2m-1階WENO格式的表達式可寫為

(5)

(6)

非線性權系數ω的表達式為

(7)

式中：Ck為理想權系數，對于五階WENO格式，其值為C0=0.1，C1=0.6，C2=0.3；ISk為衡量當地模板光滑性的光滑因子，詳細表達式見文獻[12]。

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

五階線性緊致格式的表達式為

(8)

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

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

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

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

2 泰勒-格林渦算例

2.1 算例設置

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

算例的初始流場為

(9)

下文的分析中主要用到了體平均動能Ek及其耗散率εEk這2個平均量。體平均動能的計算式為

(10)

式中：u為速度矢量；Ω表示積分域。其對應的動能耗散率的計算式為

(11)

Ek和εEk能夠反映流動發展過程中流場動能變化的情況，從而衡量模擬結果的精度。

2.2 計算結果與討論

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

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

圖4和圖5為中等網格下WENO5和compact5兩種重構格式、Roe和LMRoe兩種通量格式兩兩組合的計算結果。與粗網格結果相比，密網格下的計算結果精度更高，且2種通量格式計算結果間的差別變小。

為了定量比較不同格式的計算結果差異，以網格間距2π/256的密網格compact5計算獲得的參考解為基準，得到不同計算結果與基準解間的相對誤差，公式為

(12)

式中：g(t)表示某格式計算結果中某個變量隨時間的變化；gc(t)表示密網格compact5計算結果中對應變量隨時間的變化。

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

(13)

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

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

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

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

利用如下公式：

(14)

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

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

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

圖8和圖9為2套網格下原始Roe格式和LMRoe格式的計算差別χ′柱狀圖?？梢钥闯?，隨著網格加密，2種通量格式的計算結果間差別變小，表現出向精確解收斂的特征。

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

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

重構格式體平均動能動能耗散率WENO30.360.57 WENO50.630.61WENO70.690.75compact51.041.00

表2 網格間距2π/128時不同通量格式結果的誤差比值Table 2 Ratio of two flux schemes’ result errors with grid space being 2π/128

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

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

3 結 論

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

1) 當流場中存在低速流動時，使用低速修正的通量格式能夠一定程度上改進計算結果。

2) 通量格式對模擬結果的影響是與重構格式耦合的。當重構格式的精度較低時，通量格式對結果的影響顯著，但當重構格式的精度足夠高后，通量格式對結果的影響不明顯。

3) 隨著網格加密，有無低速修正的通量格式計算結果都呈現出向精確解收斂的特征。使用較粗網格時，低速修正的通量格式對計算結果的改進更明顯。

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