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

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

(國(guó)防科技大學(xué) 空天科學(xué)學(xué)院， 湖南 長(zhǎng)沙 410073)

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

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

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

盡管在完全氣體計(jì)算中，Tysinger[11]、曹文斌[12]等研究發(fā)現(xiàn)，對(duì)黏性項(xiàng)隱式處理可提高計(jì)算效率，然而由于高溫?zé)峄瘜W(xué)非平衡流動(dòng)控制方程包含多組分及多溫度，相比完全氣體其控制方程變量和方程數(shù)目增多，其中黏性項(xiàng)Jacobian矩陣推導(dǎo)過程復(fù)雜，國(guó)內(nèi)外高溫程序?qū)︷ば皂?xiàng)的隱式處理則較少見諸報(bào)道[4,13-15]。另外，國(guó)內(nèi)趙慧勇[16]在用于燃燒計(jì)算的化學(xué)非平衡流程序中植入了考慮黏性項(xiàng)隱式處理的LU-SGS算法，但未考慮熱非平衡效應(yīng)且未對(duì)其加速特性進(jìn)行討論。

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

1 熱化學(xué)非平衡流控制方程

(1)

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

(2)

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

(3)

2 熱化學(xué)非平衡流全隱LU-SGS算法

2.1 黏性項(xiàng)隱式處理

熱化學(xué)非平衡流原始LU-SGS算法[13-15]在隱式處理方面僅考慮化學(xué)反應(yīng)源項(xiàng)和對(duì)流項(xiàng)，即：

(4)

其中，Δt是時(shí)間步長(zhǎng)。本文在式(4)基礎(chǔ)上，新增對(duì)黏性項(xiàng)的隱式處理：

(5)

由于對(duì)流項(xiàng)僅是守恒變量的函數(shù)，將其進(jìn)行線化處理，可利用時(shí)間方向上的泰勒展開：

(6)

(7)

=Δt·RHS

(8)

(9)

其由空間離散得到，對(duì)流項(xiàng)離散格式采用AUSMPW，黏性項(xiàng)離散格式采用二階中心差分。

由于高超聲速定常流動(dòng)計(jì)算對(duì)時(shí)間精度要求不高，且黏性項(xiàng)Jacobian矩陣形式復(fù)雜，求逆需耗費(fèi)大量計(jì)算時(shí)間，因此通常采用相應(yīng)的矩陣譜半徑進(jìn)行對(duì)角近似處理，以保證左手項(xiàng)部分對(duì)角占優(yōu)。經(jīng)推導(dǎo)(推導(dǎo)過程見本文OSID拓展閱讀資料)，熱化學(xué)非平衡流計(jì)算坐標(biāo)系k方向上黏性項(xiàng)Jacobian矩陣的譜半徑為：

(10)

其中，

Di是組元擴(kuò)散系數(shù)，K、Kv分別是混合氣體平動(dòng)和振動(dòng)熱傳導(dǎo)系數(shù)，Cv是定容比熱比，Tv是振動(dòng)溫度，ev,i是分子組元的振動(dòng)能，s=mol代表所有的分子組元，詳細(xì)定義參見文獻(xiàn)[20]。

2.2 FLU-SGS算法

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

1)向前掃描：

(11)

2)向后掃描：

(12)

3)守恒量更新：

(13)

此外，本文采用當(dāng)?shù)貢r(shí)間步方法計(jì)算每個(gè)網(wǎng)格點(diǎn)的Δt。

2.3 BLU-SGS算法

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

1)向前掃描：

(14)

2)向后掃描：

(15)

(16)

3 數(shù)值結(jié)果

3.1 二維圓柱繞流

3.1.1 來流條件及計(jì)算結(jié)果

本算例基于高焓激波風(fēng)洞HEG中開展的準(zhǔn)二維圓柱繞流實(shí)驗(yàn)[21]，風(fēng)洞來流焓值21.7 MJ/kg，計(jì)算采用的氣體模型是11組元空氣，壁面條件為等溫全催化壁面，圓柱半徑R為45 mm，詳細(xì)風(fēng)洞來流條件和計(jì)算條件見表1。網(wǎng)格分布最終采用249×39，經(jīng)過網(wǎng)敏感性分析后，最終壁面第一層距離采用3×10-7m，當(dāng)?shù)芈曀倬W(wǎng)格雷諾數(shù)達(dá)到6.7，滿足熱流計(jì)算網(wǎng)格收斂性條件[3]。

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

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

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

圖2 圓柱壁面熱流對(duì)比結(jié)果Fig.2 Comparison of the heat flux of cylinder wall

3.1.2 三類LU-SGS算法收斂速率對(duì)比

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

圖3為針對(duì)圓柱算例開展的關(guān)于三種LU-SGS算法收斂速率的數(shù)值實(shí)驗(yàn)結(jié)果。結(jié)果發(fā)現(xiàn)，三類LU-SGS方法的最大CFL數(shù)不同，采用原始LU-SGS算法，最大的CFL數(shù)取值約為0.01，而FLU-SGS算法的最大CFL數(shù)為10量級(jí)，BLU-SGS的最大CFL數(shù)可取到100量級(jí)，達(dá)到機(jī)器誤差的迭代步約為FLU-SGS的73.4%。

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

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

3.1.3 CFL數(shù)對(duì)FLU-SGS算法收斂的影響

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

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

3.1.4 內(nèi)迭代數(shù)對(duì)BLU-SGS算法收斂的影響

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

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

3.2 軸對(duì)稱返回艙繞流

3.2.1 來流條件及計(jì)算結(jié)果

本算例基于Hollis在HYPULSE膨脹管風(fēng)洞對(duì)MP-1返回艙實(shí)驗(yàn)?zāi)Ｐ偷臏y(cè)量結(jié)果[22]。采用本文程序?qū)υ搶?shí)驗(yàn)進(jìn)行復(fù)現(xiàn)，實(shí)驗(yàn)?zāi)Ｐ褪欠祷嘏?支桿的構(gòu)型，其中，支桿為固定裝置，返回艙半徑R=0.025 4 m。由于模型是軸對(duì)稱構(gòu)型，圖6僅展示xy截面，其網(wǎng)格分布為208×119。來流條件為：馬赫數(shù)Ma=7.93，速度U∞=5.162 km/s；壓強(qiáng)P∞=1 824 Pa；平動(dòng)溫度T∞=1 113 K；振動(dòng)溫度Tv∞=1 113 K；來流雷諾數(shù)為668 000；來流組元為空氣(N2和O2質(zhì)量分?jǐn)?shù)分別為79%和21%)，壁面采用等溫、無滑移、非催化壁條件，壁面溫度為300 K。

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

圖6顯示了返回艙繞流的無量綱速度云圖，高速氣流撞擊返回艙前端產(chǎn)生脫體激波，波后形成強(qiáng)烈的壓縮，前體部分是附著流動(dòng)，流動(dòng)在返回艙肩點(diǎn)之后發(fā)生了分離，并在尾部區(qū)域和支桿之間形成大型分離區(qū)，分離區(qū)內(nèi)部有兩對(duì)渦對(duì)，流動(dòng)進(jìn)而在支桿上再附，形成再附激波，支桿下游保持邊界層附著流動(dòng)。

3.2.2 三類LU-SGS算法收斂速率對(duì)比

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

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

3.2.3 CFL數(shù)對(duì)FLU-SGS算法收斂的影響

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

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

3.2.4 內(nèi)迭代數(shù)對(duì)BLU-SGS算法收斂的影響

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

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

當(dāng)內(nèi)迭代數(shù)smax繼續(xù)增加到4后，總密度殘差隨迭代步的收斂速度加快，計(jì)算機(jī)時(shí)有所增加。因此在實(shí)際復(fù)雜工程外形計(jì)算中，為兼顧計(jì)算穩(wěn)定性和計(jì)算效率，選用BLU-SGS算法最佳，且在滿足收斂到機(jī)器誤差的要求下，內(nèi)迭代數(shù)smax視情況可選取2至4，因?yàn)樵搮?shù)取值過大對(duì)收斂性沒有明顯的增益，且浪費(fèi)計(jì)算機(jī)時(shí)。

4 結(jié)論

1)在高溫?zé)峄瘜W(xué)非平衡流計(jì)算中，考慮黏性項(xiàng)隱式處理的兩種全隱LU-SGS算法引入了跨越邊界層的黏性影響，能夠快速建立強(qiáng)黏性干擾初場(chǎng)，實(shí)現(xiàn)初始最大CFL數(shù)3～5個(gè)數(shù)量級(jí)的提升。

2)BLU-SGS算法通過增加內(nèi)迭代數(shù)能提升計(jì)算穩(wěn)定性，在高超聲速熱化學(xué)非平衡計(jì)算中準(zhǔn)確性高、可靠性強(qiáng)、收斂速率快，適合復(fù)雜流場(chǎng)計(jì)算。

3)隨著CFL數(shù)升高，計(jì)算穩(wěn)定性有所降低，殘差曲線表現(xiàn)小范圍周期振蕩，該現(xiàn)象與時(shí)間推進(jìn)格式、對(duì)流格式、網(wǎng)格分布、流場(chǎng)非定常特性等因素相關(guān)，要進(jìn)一步抑制該現(xiàn)象，需從構(gòu)造適用于熱化學(xué)非平衡流動(dòng)的數(shù)值格式、提高數(shù)值格式精度和提升網(wǎng)格質(zhì)量入手。

致謝

中國(guó)空氣動(dòng)力研究與發(fā)展中心趙慧勇研究員在公式推導(dǎo)方面提供了幫助和指導(dǎo)，謹(jǐn)致謝意！