全隱LU-SGS算法在高超聲速熱化學非平衡流剛性問題中的應用*

蔣 浩，柳 軍，王君媛，黃 偉，杜 洋

(國防科技大學 空天科學學院， 湖南 長沙 410073)

下一代天地往返空天飛行器具有高馬赫數、寬速域的特點，穿越空域從對流層、平流層變化到臨近空間，更為復雜的飛行任務剖面給飛行器氣動、結構和熱防護系統設計等帶來巨大挑戰，在設計中必須考慮高溫氣體效應的影響[1]。當飛行器高速再入時，在頭部脫體激波及邊界層的強黏性干擾作用下，空氣被加熱至數千度高溫，其中的氧氣和氮氣組分發生振動激發、離解甚至電離，成為由分子、原子、離子和電子組成的混合氣體，且當流動特征時間與能量松弛、組分化學反應特征時間相互比擬時，流動被稱作熱化學非平衡流。美國在早期航天飛機設計中依據完全氣體假設，在試飛實驗中出現了配平攻角高出設計值一倍多的氣動異常現象[2]。

因此，在依靠計算流體力學(Computational Fluid Dynamics, CFD)工具進行高超聲速飛行器氣動預測時，常規的完全氣體假設不再適用，必須采用考慮熱化學非平衡的高溫氣體假設，但熱化學非平衡計算存在嚴重的數值剛性問題，即數值計算失穩或收斂困難。其中，第一類是由源項帶來的剛性，這是由于局部流場的流動特征時間可能與能量松弛、化學反應特征時間量級尺度存在極大差異，使得控制方程中相應源項的量級相差太大，進而在時間推進求解時出現收斂困難甚至發散的問題。第二類是由網格加密帶來的剛性，熱化學非平衡流動計算特別是熱流計算對壁面網格加密要求極高，通常要求當地聲速網格雷諾數小于10[3]；另外，激波/邊界層干擾等復雜流動要求在局部流場參數梯度較大的區域進行網格加密，這進一步加劇了網格剛性的問題。基于上述兩種因素分析，熱化學非平衡流程序相比完全氣體程序，穩定性和魯棒性較差，為滿足穩定性限制，熱化學非平衡流動計算庫朗數 (Courant-Friedrichs-Lewy, CFL)一般取值較小，特別在采用顯式時間推進格式時，該問題極為突出[4]。為此，在熱化學非平衡流實際計算中一般選用穩定性較好的隱式格式，通過預處理方法以放寬線性方程組迭代求解對時間步長的限制。

美國從20世紀80年代開始熱化學非平衡CFD研究，為了解決化學反應源項的剛性問題，Bussing 和Murman[5]提出將化學反應源項進行隱式處理，由于隱式項不需進行空間差分離散處理，只考慮本單元數據對本單元殘差的貢獻，故該算法被稱作點隱式格式。在此基礎上， Eberhardt和Imlay[6]提出了化學反應源項Jacobian矩陣的對角化形式，簡化了化學反應源項的隱式處理。為進一步提高CFL數，有必要增加對流項的隱式處理，Yoon和Jameson[7]提出了LU-SGS(lower-upper symmetric Gauss-Seidel)算法，該算法通過正負分裂對流項Jacobian矩陣，將線性方程組左手項分解為三個子矩陣，因而避免了復雜矩陣求逆過程。LU-SGS算法最初被用于跨聲速的流動求解，之后被推廣到高超聲速流動計算中，成為熱化學非平衡流定常計算中最為通用的隱式時間推進算法。Chen和Wang[8]在完全氣體LU-SGS算法基礎上發展了BLU-SGS (block LU-SGS)算法，該算法在保證隱式系統對角特性的基礎上，在每次更新中嵌入內迭代過程，通過內迭代過程中的上掃和下掃引入非對角塊的貢獻。另外，相比LU-SGS這類線性迭代求解算法，學者還研究了GMRES (generalized minimum residual) 等非線性迭代求解算法[9-10]，進一步提高了高超聲速流動計算效率。

盡管在完全氣體計算中，Tysinger[11]、曹文斌[12]等研究發現，對黏性項隱式處理可提高計算效率，然而由于高溫熱化學非平衡流動控制方程包含多組分及多溫度，相比完全氣體其控制方程變量和方程數目增多，其中黏性項Jacobian矩陣推導過程復雜，國內外高溫程序對黏性項的隱式處理則較少見諸報道[4,13-15]。另外，國內趙慧勇[16]在用于燃燒計算的化學非平衡流程序中植入了考慮黏性項隱式處理的LU-SGS算法，但未考慮熱非平衡效應且未對其加速特性進行討論。

本文基于結構網格有限差分法，總結提出熱化學非平衡流動計算中強數值剛性帶來的收斂困難問題；給出熱化學非平衡流黏性項Jacobian矩陣的推導結果，并在時間推進算法中對該部分矩陣采取對角簡化形式，實現對化學反應源項、對流項、黏性項的全隱處理，將完全氣體條件下的FLU-SGS及BLU-SGS算法推廣到熱化學非平衡流動計算中；針對高超聲速熱化學非平衡二維圓柱實驗和軸對稱返回艙實驗算例，對比算法改進前后的加速收斂特性。

1 熱化學非平衡流控制方程

(1)

U=(ρ,ρi,ρu,ρv,ρet,ρev)T

(2)

其中：E和F是對流項；Ev和Fv是黏性項；W是源項；ρ是總密度；ρi分別為11個組元的密度；u和v是速度分量；et和ev分別為混合氣體的單位質量總能和振動能。將上述控制方程進行坐標變換，得到計算坐標系下的形式：

(3)

2 熱化學非平衡流全隱LU-SGS算法

2.1 黏性項隱式處理

熱化學非平衡流原始LU-SGS算法[13-15]在隱式處理方面僅考慮化學反應源項和對流項，即：

(4)

其中，Δt是時間步長。本文在式(4)基礎上，新增對黏性項的隱式處理：

(5)

由于對流項僅是守恒變量的函數，將其進行線化處理，可利用時間方向上的泰勒展開：

(6)

(7)

=Δt·RHS

(8)

(9)

其由空間離散得到，對流項離散格式采用AUSMPW，黏性項離散格式采用二階中心差分。

由于高超聲速定常流動計算對時間精度要求不高，且黏性項Jacobian矩陣形式復雜，求逆需耗費大量計算時間，因此通常采用相應的矩陣譜半徑進行對角近似處理，以保證左手項部分對角占優。經推導(推導過程見本文OSID拓展閱讀資料)，熱化學非平衡流計算坐標系k方向上黏性項Jacobian矩陣的譜半徑為：

(10)

其中，

Di是組元擴散系數，K、Kv分別是混合氣體平動和振動熱傳導系數，Cv是定容比熱比，Tv是振動溫度，ev,i是分子組元的振動能，s=mol代表所有的分子組元，詳細定義參見文獻[20]。

2.2 FLU-SGS算法

參考原始LU-SGS 算法，對式(8)進行近似LU分解，再使用一次對稱Gauss-Seidel迭代，則時間推進可化為三個步驟：

1)向前掃描：

(11)

2)向后掃描：

(12)

3)守恒量更新：

(13)

此外，本文采用當地時間步方法計算每個網格點的Δt。

2.3 BLU-SGS算法

與LU-SGS算法不同的是，BLU-SGS算法在每次守恒量更新步內還包含了smax次松弛迭代。若smax=1，則退化到LU-SGS算法；對于不同的算例，內迭代數smax存在著不同的最優值[11]。第s次內迭代表達式如下：

1)向前掃描：

(14)

2)向后掃描：

(15)

(16)

3 數值結果

3.1 二維圓柱繞流

3.1.1 來流條件及計算結果

本算例基于高焓激波風洞HEG中開展的準二維圓柱繞流實驗[21]，風洞來流焓值21.7 MJ/kg，計算采用的氣體模型是11組元空氣，壁面條件為等溫全催化壁面，圓柱半徑R為45 mm，詳細風洞來流條件和計算條件見表1。網格分布最終采用249×39，經過網敏感性分析后，最終壁面第一層距離采用3×10-7m，當地聲速網格雷諾數達到6.7，滿足熱流計算網格收斂性條件[3]。

表1 二維圓柱算例來流條件

圖1為圓柱繞流無量綱壓力云圖，結果顯示本文程序能較好地捕捉到激波位置，流場中各物理量計算結果比較合理。圖2是圓柱壁面熱流Qw計算結果與實驗結果對比，可以看出壁面熱流計算值與實驗值符合良好，驗證了本文數值方法的正確性。

圖1 圓柱繞流無量綱壓力云圖Fig.1 Dimensionless pressure contour for cylinder case

圖2 圓柱壁面熱流對比結果Fig.2 Comparison of the heat flux of cylinder wall

3.1.2 三類LU-SGS算法收斂速率對比

高超聲速計算普遍采用來流條件作為初場，計算流場普遍經歷以下過程：從流動在壁面附近滯止、形成脫體激波或分流區，直至收斂呈現為穩定的流場結構，其中流場基本結構的建立是整個計算過程中數值剛性最強的步驟。因此，在計算初期只能采用較低的CFL數以避免發散，待激波、邊界層、分離等流動特征穩定建立后，才能逐步增大CFL數。在此，首先給出最大CFL數的定義：在該類格式下密度殘差Δρ能收斂到機器誤差(10-12)的CFL數的最大值。為比較三類算法在計算初期所能采用的最大CFL數，在數值實驗中采用恒定CFL數進行計算。

圖3為針對圓柱算例開展的關于三種LU-SGS算法收斂速率的數值實驗結果。結果發現，三類LU-SGS方法的最大CFL數不同，采用原始LU-SGS算法，最大的CFL數取值約為0.01，而FLU-SGS算法的最大CFL數為10量級，BLU-SGS的最大CFL數可取到100量級，達到機器誤差的迭代步約為FLU-SGS的73.4%。

圖3 圓柱算例三種LU-SGS算法收斂速率對比Fig.3 Comparison of the convergence rate between three LU-SGS algorithms for cylinder case

因此，增加黏性項隱式處理的FLU-SGS和BLU-SGS算法由于在計算中考慮了黏性干擾的流場信息，故能快速建立邊界層。兩類改進算法在20 000步前的收斂速率較大，而在此之后由于流動結構基本不變，收斂速率出現拐點，密度殘差的收斂曲線近似呈線性，直至收斂到機器誤差。另外，在近期文獻[14]報道中，也對本節同樣算例采用LU-SGS算法進行了熱化學平衡流場計算，其CFL數取值為5，可見本文改進算法對提升收斂速率有較大優勢。

3.1.3 CFL數對FLU-SGS算法收斂的影響

圖4 對比了CFL數為10、100和1 000下采用FLU-SGS算法的收斂速率。結果表明，隨著CFL數取大，前期收斂速率加快。盡管采用全隱算法后，初始CFL數能夠取到1 000，但當殘差降低到10-8后呈現周期振蕩現象；當CFL數進一步取到10后，該種振蕩現象能夠消除，并能降低到機器誤差。另外，雖然在大CFL數條件下的收斂后期，密度殘差會出現周期振蕩，但當其降低到10-6以下后，流動結構和壁面參數的實際變化很小，可近似認為收斂。

圖4 圓柱算例中CFL數對FLU-SGS收斂速率的影響Fig.4 Effect of CFL number on the convergence rate of FLU-SGS for cylinder case

3.1.4 內迭代數對BLU-SGS算法收斂的影響

圖5是內迭代數smax對BLU-SGS收斂速率的影響，可以看出在大CFL數條件下，增加smax抑制了計算后期的殘差振蕩，增加了收斂的穩定性。雖然smax越大，密度殘差隨迭代步的收斂速度加快，但是從計算機時角度看，增加smax耗時更多。

(a) 密度殘差隨迭代步的變化(a) Variation of density residual with iteration number

3.2 軸對稱返回艙繞流

3.2.1 來流條件及計算結果

本算例基于Hollis在HYPULSE膨脹管風洞對MP-1返回艙實驗模型的測量結果[22]。采用本文程序對該實驗進行復現，實驗模型是返回艙-支桿的構型，其中，支桿為固定裝置，返回艙半徑R=0.025 4 m。由于模型是軸對稱構型，圖6僅展示xy截面，其網格分布為208×119。來流條件為：馬赫數Ma=7.93，速度U∞=5.162 km/s；壓強P∞=1 824 Pa；平動溫度T∞=1 113 K；振動溫度Tv∞=1 113 K；來流雷諾數為668 000；來流組元為空氣(N2和O2質量分數分別為79%和21%)，壁面采用等溫、無滑移、非催化壁條件，壁面溫度為300 K。

圖6 返回艙算例無量綱x方向速度云圖 Fig.6 Dimensionless x-velocity contour for re-entry capsule case

圖6顯示了返回艙繞流的無量綱速度云圖，高速氣流撞擊返回艙前端產生脫體激波，波后形成強烈的壓縮，前體部分是附著流動，流動在返回艙肩點之后發生了分離，并在尾部區域和支桿之間形成大型分離區，分離區內部有兩對渦對，流動進而在支桿上再附，形成再附激波，支桿下游保持邊界層附著流動。

3.2.2 三類LU-SGS算法收斂速率對比

圖7對比了三類LU-SGS方法的收斂速率。由于返回艙繞流比圓柱繞流的流場結構更為復雜，在原始LU-SGS算法下的初始計算中，最大CFL數僅能取到10-4量級，因而收斂過程極為緩慢。采用FLU-SGS算法能將最大CFL數大幅提升到1的量級，但是后期收斂速率有所降低。BLU-SGS算法最大CFL數量級能達到10，能實現快速收斂，達到機器誤差的迭代步僅為FLU-SGS的8.5%。

圖7 返回艙算例三種LU-SGS算法收斂速率對比Fig.7 Comparison of the convergence rate between three LU-SGS algorithms for re-entry capsule case

3.2.3 CFL數對FLU-SGS算法收斂的影響

圖8對比了CFL數為1、 10和100下采用FLU-SGS算法的收斂速率。結果表明，隨著CFL數取大，前期收斂速率加快。同圓柱算例類似，在大CFL數下，殘差降低到10-3后呈現周期振蕩現象，當CFL數進一步取到1后，該種振蕩現象能夠消除，并能降低到機器誤差。由于采用二階格式，該類殘差振蕩現象與過激波參數在兩個網格間跳動等因素相關。

圖8 返回艙算例CFL數對FLU-SGS收斂速率的影響Fig.8 Effect of CFL number on the convergence rate of FLU-SGS for re-entry capsule case

3.2.4 內迭代數對BLU-SGS算法收斂的影響

圖9是內迭代數smax對BLU-SGS算法收斂速率的影響，可以看出在CFL數為10、smax=1時，BLU-SGS算法退化為FLU-SGS算法，內迭代數smax增加到3時才能抑制計算后期的殘差振蕩。

(a) 密度殘差隨迭代步的變化(a) Variation of density residual with iteration number

當內迭代數smax繼續增加到4后，總密度殘差隨迭代步的收斂速度加快，計算機時有所增加。因此在實際復雜工程外形計算中，為兼顧計算穩定性和計算效率，選用BLU-SGS算法最佳，且在滿足收斂到機器誤差的要求下，內迭代數smax視情況可選取2至4，因為該參數取值過大對收斂性沒有明顯的增益，且浪費計算機時。

4 結論

1)在高溫熱化學非平衡流計算中，考慮黏性項隱式處理的兩種全隱LU-SGS算法引入了跨越邊界層的黏性影響，能夠快速建立強黏性干擾初場，實現初始最大CFL數3～5個數量級的提升。

2)BLU-SGS算法通過增加內迭代數能提升計算穩定性，在高超聲速熱化學非平衡計算中準確性高、可靠性強、收斂速率快，適合復雜流場計算。

3)隨著CFL數升高，計算穩定性有所降低，殘差曲線表現小范圍周期振蕩，該現象與時間推進格式、對流格式、網格分布、流場非定常特性等因素相關，要進一步抑制該現象，需從構造適用于熱化學非平衡流動的數值格式、提高數值格式精度和提升網格質量入手。

致謝

中國空氣動力研究與發展中心趙慧勇研究員在公式推導方面提供了幫助和指導，謹致謝意！