激波內(nèi)部結(jié)構(gòu)的數(shù)值求解方法1)

朱清波 周文元 楊慶春 徐 旭

* (北京航空航天大學(xué)宇航學(xué)院,北京 102206)

? (上?？臻g推進研究所,上海 201112)

引言

在氣體動力學(xué)中,激波常被簡化為零厚度的間斷面,經(jīng)過激波,流動參數(shù)發(fā)生階躍,熵不可逆地增加,這對大部分情況是很好的近似.然而真實的宏觀流動中不會出現(xiàn)數(shù)學(xué)意義上的不連續(xù),激波也不例外,從波前向波后狀態(tài)的過渡必然是連續(xù)的,而激波厚度則是一個有限值.圖1 展示了一個典型駐定正激波的內(nèi)部流動參數(shù)變化曲線,其中Ma1和Ma2分別為波前和波后馬赫數(shù),δ 為基于最大速度梯度定義的激波厚度.

圖1 激波內(nèi)部結(jié)構(gòu)示意圖Fig.1 Schematic diagram of the internal structure of shock waves

研究激波內(nèi)部流動參數(shù)變化過程的問題被稱為激波結(jié)構(gòu)問題,它不僅是一個重要的科學(xué)問題,還具有工程實用價值.在稀薄氣體動力學(xué)[1]、分子動理論[2]、非平衡態(tài)熱力學(xué)[3]和天體物理[4]等基礎(chǔ)科學(xué)領(lǐng)域中,激波結(jié)構(gòu)問題促進了對激波、輸運現(xiàn)象以及流體力學(xué)理論基礎(chǔ)的理解.在航空航天[5-6]、光學(xué)[7]以及核[8]等技術(shù)領(lǐng)域中,激波結(jié)構(gòu)的正確解析對解決一些工程問題至關(guān)重要.例如航天器再入[5-6]時,高層大氣極其稀薄導(dǎo)致激波厚度較大,準(zhǔn)確模擬激波內(nèi)部結(jié)構(gòu)是外流場仿真的關(guān)鍵之一;又比如,激波會使光發(fā)生偏折,正確描述激波結(jié)構(gòu)有助于修正它對光學(xué)系統(tǒng) (風(fēng)洞試驗的光學(xué)測量裝置、超聲速飛行器的光學(xué)觀瞄設(shè)備與機載激光武器等) 產(chǎn)生的不利影響[7].此外,駐定正激波是一個一維定常且邊界條件極其簡單的強非平衡流動,這使其成為驗證氣體動力學(xué)模型及相關(guān)數(shù)值方法準(zhǔn)確性與可靠性的理想場景,再加上豐富的實驗數(shù)據(jù)[9-10],激波結(jié)構(gòu)長期被用作嚴(yán)苛而敏感的標(biāo)準(zhǔn)測試算例[11-14].

激波結(jié)構(gòu)的求解有解析和數(shù)值兩種方式.一般來說,包括激波結(jié)構(gòu)問題在內(nèi)的大部分流動問題非常復(fù)雜,難以求得解析解.但Becker[15]發(fā)現(xiàn)當(dāng)氣體的Prandtl 數(shù)等于3/4 這個特殊值時,激波內(nèi)部流動所滿足的能量方程可以被簡化,總焓沿流向保持不變,進而獲得了第一個對任意強度激波有效的解析解.Becker 解隨后被不斷改進,從最初的只適用于輸運系數(shù)為常數(shù)的氣體,到適用于硬球和Maxwell 分子模型氣體[16]、軟球分子模型氣體[17]和van der Waals 氣體[18]等.盡管適用范圍不斷擴大,這些解析解對Prandtl 數(shù)的要求卻從未放寬.另一個被廣泛使用的解析解是Mott-Smith[19]提出的著名的雙峰分布函數(shù)解,然而它只是一定假設(shè)下的近似解,且從原理上看不能很好地描述弱激波.實際上,即使對強激波,Mott-Smith 解也不被實驗和直接模擬Monte Carlo(DSMC) 計算的結(jié)果支持[10].總而言之,目前的解析方法仍存在較多局限性,只能解決一些特殊條件下的簡單問題,而大部分激波結(jié)構(gòu)問題較為復(fù)雜,需借助數(shù)值方法求解.

隨著計算科學(xué)的發(fā)展,數(shù)值技術(shù)在流體力學(xué)中發(fā)揮著越來越重要的作用,激波結(jié)構(gòu)問題也受益于此.按所采用的數(shù)學(xué)物理模型的類型來劃分,激波結(jié)構(gòu)的數(shù)值模擬方法主要有3 種:

(1) 微觀方法,即基于還原論觀點、從第一性原理出發(fā)直接模擬分子運動和碰撞的粒子方法,主要包括以分子動力學(xué)[20](MD) 為代表的確定論方法和以DSMC[21-22]為代表的概率方法;

(2) 介觀方法,求解基于速度分布函數(shù)描述的Boltzmann 方程[23-24],或其簡化方程/模型[25],包括著名的BGK 模型方程[26]和離散速度模型[27-28](DVM) 等;

(3) 宏觀方法,求解基于宏觀量描述的連續(xù)介質(zhì)流體力學(xué)方程/模型,包括經(jīng)典和改進的Navier-Stokes 方程[11,29-30]、以Burnett 方程為代表的高階流體力學(xué)方程[31-32]、Grad 矩方程及其變體[14,33]、廣義流體力學(xué)[34](generalized hydrodynamics) 方程、非線性耦合本構(gòu)關(guān)系[35](NCCR) 模型以及擴展/廣延熱力學(xué)[3](extended thermodynamics) 模型等.

相對于Boltzmann 方程,控制氣體宏觀運動的流體力學(xué)方程是一種降階的描述,自由度更少,因而更簡單.例如在Chapman-Enskog 理論[36]中,Navier-Stokes 方程和Burnett 方程分別是對Boltzmann 方程的一階和二階近似.因此,第3 類方法原則上應(yīng)歸入第2 類.但從介觀Boltzmann 方程到宏觀流體力學(xué)方程的粗化并不容易,兩者之間的關(guān)系在數(shù)學(xué)上也遠非顯而易見 (詳見Hilbert 第6 問題[37]).更重要的是,流體力學(xué)方法基于與分布函數(shù)描述迥然不同的局部統(tǒng)計平均量 (分布函數(shù)的矩) 描述,體系相對獨立而內(nèi)容又極其豐富,故單獨列為一類.

前述3 類方法中,基于流體力學(xué)方程的數(shù)值模擬最受青睞,因為它在具有良好精度的同時還有較高的計算效率,其宏觀描述也為人們樂于接受,因此相關(guān)算法得到高度發(fā)展.具體到激波結(jié)構(gòu)的仿真計算中,有限差分法常被使用,這類通用的計算流體力學(xué)技術(shù)十分成熟,此處不再贅述.

由于以下3 個原因,絕大多數(shù)激波結(jié)構(gòu)問題可歸結(jié)為一維定常正激波的結(jié)構(gòu)問題.

(1) 激波的典型厚度非常小,流體穿過激波的時間極短,一般僅需數(shù)次分子碰撞即可完成從波前到波后狀態(tài)的過渡,弛豫速度遠快于其他宏觀非平衡流動現(xiàn)象.因此,對于大部分激波結(jié)構(gòu)問題,運動隨時間的發(fā)展不重要,非定常效應(yīng)可以忽略.

(2) 平面激波流動可以被視為法向的正激波流動疊加切向的平行流動.因此,二維或三維平面激波的結(jié)構(gòu)問題可被簡化為法向上的一維正激波結(jié)構(gòu)問題.

(3) 曲面激波波面的最小曲率半徑一般遠大于激波厚度,所以曲率對激波內(nèi)部結(jié)構(gòu)的影響通?？梢院雎?曲面激波的局部可被近似為平面激波,進而分解為法向上的正激波和切向的平行流.

毫無疑問,一維定常正激波是檢驗氣體動力學(xué)模型 (尤其是稀薄氣體動力學(xué)模型) 的理想算例: 極小的厚度意味著Knudsen 數(shù)較大,間斷分子效應(yīng)顯著;強耗散意味著遠離熱力學(xué)平衡,非平衡效應(yīng)顯著;一維定常的特性和簡單的邊界條件還降低了問題的復(fù)雜性.定常正激波的結(jié)構(gòu)問題,其本質(zhì)是求解一組流體力學(xué)常微分方程.對這類問題,直接數(shù)值積分比有限差分法更簡單高效.然而對激波結(jié)構(gòu)應(yīng)用空間推進的數(shù)值積分仍存在一些問題,例如沿流向推進計算的發(fā)散問題,本文的目的是解決這些問題,并建立在流體力學(xué)框架內(nèi)數(shù)值求解激波結(jié)構(gòu)的一般程序.

由于關(guān)注的問題具有普遍性,本文將只討論最基礎(chǔ)的情況,即單組分簡單氣體的情況,多組分氣體和等離子體的情況暫不考慮,以避免不必要的細節(jié)把問題復(fù)雜化.控制方程將采用簡單的Navier-Stokes 方程,這不會導(dǎo)致本文的討論失去一般性,因為宏觀連續(xù)的流體系統(tǒng)都具有Navier-Stokes 形式的守恒方程,區(qū)別僅在于用怎樣的本構(gòu)方程封閉它們,而這不影響本文的主要結(jié)論.

1 控制方程與邊界條件

一維氣體流動的質(zhì)量、動量和能量守恒方程分別可寫為

其中,t,x,ρ,V,p,u,ht,τ 和q分別為時間、流向坐標(biāo)、密度、速度、壓強、比內(nèi)能、比總焓、黏性應(yīng)力以及熱流密度.需要說明的是,為保持與氣體運動論領(lǐng)域的習(xí)慣一致,本文約定當(dāng) τ 表現(xiàn)為壓應(yīng)力時其符號為正,這與流體力學(xué)中法向應(yīng)力的定義相反.式(1)～式(3) 直接來源于三大守恒定律,也可從Boltzmann 方程的矩方程——Maxwell 輸運方程導(dǎo)出,未對 τ 和q作任何假設(shè),因此被所有連續(xù)介質(zhì)氣體動力學(xué)模型的一維形式所遵守.

若采用Navier-Stokes 和Fourier 線性本構(gòu)關(guān)系,有

式中,T是溫度,μ 和 κ 分別為動力黏度和導(dǎo)熱系數(shù).將式 (4) 代入式(1)～式(3) 即得到一維可壓Navier-Stokes 方程組.

考慮如圖1 所示的坐標(biāo)系固定在波面上的駐定正激波情況,由于流動是定常的,時間偏導(dǎo)項為0,所有參數(shù)僅與位置x有關(guān),式(1)～式(3) 可簡化為常微分方程,將所得常微分方程對x積分,得到一維定常正激波內(nèi)部流動的控制方程

式中,Jm,JP和JE是積分常數(shù),不隨x改變.實際上它們分別是質(zhì)量、動量和能量的通量,并且因為在波前波后無窮遠處 τ 和q趨于0,根據(jù)式(5)～式(7)顯然有

下標(biāo)1 和2 分別表示上游和下游無窮遠處的狀態(tài).

進一步引入理想氣體狀態(tài)方程

其中R為氣體常數(shù).式(4)～式(7) 與式(11) 組成一個封閉的微分方程組,從中消去 τ,q,p以及 ρ 可得

該系統(tǒng)滿足以下漸近邊界條件

其中V2和T2可根據(jù)來流參數(shù)V1,T1以及正激波前后參數(shù)的關(guān)系求出.

2 打靶法

式 (12)～式(13) 構(gòu)成一個封閉的一階常微分方程組,沿x對其積分可獲得V和T的沿程分布.然而邊界條件 (14) 給在無窮遠處,無法在數(shù)值積分中直接使用,人們很自然地想到將無限域截斷成有限的計算域,再用打靶法[6,30,38]將邊值問題轉(zhuǎn)化為初值問題求解,其具體步驟如下.

(1) 選一略小于V1的速度值,將其對應(yīng)的位置作為積分的起點,坐標(biāo)指定為x=0,該速度值記為V(0).

(2) 選一略大于T1的溫度值作為T(0) 的猜測值.

(3) 以V(0) 和T(0) 為初值,對式(12)～式(13) 構(gòu)成的系統(tǒng)進行數(shù)值積分,積分沿著流動的方向進行,向下游推進足夠遠后截斷.

(4) 若第3 步積分出的V和T的曲線像圖1 給出的分布曲線那樣逐漸趨平,則結(jié)束計算并輸出結(jié)果;若積分過程中發(fā)散或出現(xiàn)參數(shù)明顯遠離波后值的趨勢,則中斷計算、修正T(0) 的值并回到第3 步.

第1 步的目的是對計算域的上游邊界進行人工截斷;第3 步可采用各種數(shù)值積分算法,包括但不限于Runge-Kutta 法;第4 步的目的是對與V(0) 相對應(yīng)的初值T(0) 進行迭代修正.不同文獻使用的算法在細節(jié)上可能有所差異,例如文獻 [6,30,38] 直接求解的變量是p和 ρ,積分的起始點則選在p=(p1+p2)/2處,并從該點分別向上下游推進,這些文獻的方法與本節(jié)的方法沒有本質(zhì)區(qū)別.

2.1 打靶法的失效

用前述的打靶法對單原子氣體中Ma1=2 的正激波結(jié)構(gòu)進行了求解,數(shù)值積分采用自適應(yīng)步長的4 級Runge-Kutta 法,計算結(jié)果如圖2 所示.結(jié)果表明,如果計算程序輸入的參數(shù)組合經(jīng)過精細調(diào)整,打靶法能給出看似合理的解,然而這依賴于對T(0) 值和下游截斷位置極其謹(jǐn)慎的選擇和微調(diào),T(0) 稍有改變或計算域繼續(xù)向下游延伸就可能導(dǎo)致發(fā)散.具體到圖2 的算例,在大約x/λ1=12 (λ 表示分子平均自由程) 的位置之前,隨著計算的推進,各參數(shù)快速趨近于波后值,并在波后值附近保持了一段距離的“穩(wěn)定”,此時參數(shù)變化的趨勢以及參數(shù)分布曲線的形態(tài)是正常的,在該位置終止計算可以得到看似收斂的解.但如果繼續(xù)向下游積分,原本將要趨平的參數(shù)會快速偏離波后值,出現(xiàn)發(fā)散.

圖2 打靶法計算的單原子氣體激波結(jié)構(gòu) (Ma1=2)Fig.2 Profile of a Ma1=2 shock in a monoatomic gas calculated with the shooting method

合格的算法不應(yīng)過于依賴人工介入,然而在激波結(jié)構(gòu)問題中,打靶法很難像求解一般的兩點邊值問題那樣使用二分法、割線法之類的迭代方法自動修正初值,因為第4 步極其依賴人工判斷.此外,自動修正T(0) 需要第3 步中有一個明確的截斷標(biāo)準(zhǔn),但實際上計算域的長度很難確定,因為在計算中發(fā)現(xiàn),無論T(0) 如何接近其準(zhǔn)確值,只要向下游推進足夠遠,發(fā)散總是不可避免;過短的計算域無法保證覆蓋整個激波結(jié)構(gòu),過長的計算域則導(dǎo)致初值必須非常準(zhǔn)確才能保證收斂,因此對下游邊界的截斷和T(0)的迭代都依賴手動操作.高度的敏感性暗示著不穩(wěn)定性的存在,需要在深入分析系統(tǒng)的動力學(xué)性質(zhì)的基礎(chǔ)上提出改進算法.

2.2 積分方向的重要性

激波內(nèi)任意位置氣體的狀態(tài)可以表示為相空間中的一個點,那么激波壓縮過程可以表示為相空間中一條光滑軌線.分別記波前和波后狀態(tài)點為1 點和2 點,求解激波結(jié)構(gòu)實際就是要找到一條連接1,2 兩點的積分曲線,將V-T相平面中這樣一條曲線叫做激波過程線 (即圖3 中的曲線S).已知上下游無窮遠處速度和溫度梯度趨于0,即1,2 兩點使式 (12)～式(13) 的右側(cè)均為0,因此它們都是奇點,激波過程線的斜率 dT/dV在這兩點處是 0/0 型未定式.

圖3 V-T 平面中的相軌跡Fig.3 Phase portrait in the V-T plane

約定x增加的方向/流動方向/壓縮進行的方向為激波過程的正向.為了更直觀地了解系統(tǒng)的動力學(xué)特性,根據(jù)式 (12)～式(13) 所確定的方向場/斜率場,畫出V-T平面中的相軌跡圖,如圖3 所示,其中箭頭表示正向.相軌跡即圖中帶箭頭的實線,是隨著x的增加、相點的運動軌跡,每一條相軌跡都代表系統(tǒng)一個可能的解;相軌跡上每一點處的斜率與斜率場給定的一致,即相點的運動方向由方向場決定.相軌跡圖非常有用,它為微分方程建立了一種類似于流線圖的圖案,揭示了系統(tǒng)的定性性質(zhì),通過它可以對解的行為形成圖形化的認(rèn)識.

從圖3 可以看出,1 點是不穩(wěn)定結(jié)點,2 點是鞍點,激波過程線S是從1 點發(fā)出并最終抵達2 點的異宿軌.在相平面中創(chuàng)建曲線M和N分別代表μ=0和 κ=0 情況下的激波過程線.在V-T相平面中,曲線M的方程可通過令式 (12) 右側(cè)中括號內(nèi)的項等于0 獲得;同理曲線N的方程也可通過令式 (13) 右側(cè)中括號內(nèi)的項等于0 獲得.由各自的曲線方程易知,M是拋物線;如果氣體是定比熱的,N也是拋物線.將M和N圍成的區(qū)域命名為L.文獻 [39] 證明了 ρ-p圖中的激波過程線位于曲線M和N之間;類似地,很容易證明V-T圖中的激波過程線S也落在M和N之間,即區(qū)域L中.

如圖3 所示,從區(qū)域L內(nèi)任意一點出發(fā)進行正向積分,除非該點嚴(yán)格位于所要求解的目標(biāo)軌線S上,否則相點會先靠近2 點,但隨后又迅速遠離,始終無法抵達2 點,這也與圖2 算例發(fā)散的具體表現(xiàn)(流動參數(shù)先趨近波后值,隨后又迅速偏離) 吻合.實際計算中,即使初值點剛好落在目標(biāo)軌線上,不可避免的舍入誤差也會導(dǎo)致積分曲線偏離目標(biāo)軌線,并被導(dǎo)出L區(qū)域.波后點是鞍點,這是所有正向求解隨著數(shù)值積分向x=+∞ 推進最終都會發(fā)散的根本原因.有理由認(rèn)為,所有包含正向推進的激波結(jié)構(gòu)計算方法所獲得的“收斂”結(jié)果,都是在發(fā)散之前的“合適”位置強行終止計算所產(chǎn)生的假象,這種人為截斷既不合理,也具有很大的隨意性.

3 逆向推進計算方法

前已述及,正向推進必然導(dǎo)致數(shù)值積分曲線偏離激波過程線、相點偏離2 點.這個問題可用逆向推進解決,即從區(qū)域L內(nèi)、2 點附近選擇一初值點,從該點出發(fā)向上游積分至流動參數(shù)幾乎不再變化為止.圖3 顯示逆向積分曲線必然匯聚于1 點,所以該求解思路是可行的,但仍存在一個問題: 對初值不敏感固然是優(yōu)點,卻也導(dǎo)致逆向求解無法與打靶法搭配使用,因為打靶法靠迭代確定合適的初值,從而逼近精確解,而逆向求解會使初值誤差衰減、相點無限趨近于1 點,無法為初值的修正提供反饋信息.因此需要一種新的初值確定方法.

3.1 基于L’H?pital 法則的初值確定方法

雖然2 點的參數(shù)無法直接作為初值使用,但如果已知激波過程線的終點斜率 (dT/dV)2,就能以較高的精度在2 點附近確定一初值點 (V0,T0),其中V0的值由人為指定且略大于V2,T0可按Euler 格式給出

將式 (13) 和式 (12) 兩端分別相除,整理后得激波過程線在V-T圖中斜率的表達式

式 中JPm和JEm是常數(shù):JPm≡JP/Jm,JEm≡JE/Jm.注意理想氣體的比內(nèi)能u是且僅是溫度T的函數(shù),且du=cvdT,其中cv是定容比熱;μ 和 κ 一般也是溫度T的函數(shù).

2.2 節(jié)提到,1,2 兩點是奇點,直接將波前或波后參數(shù)代入式 (16) 會使f/g成為 0/0 型未定式,給兩端點處激波過程線斜率的確定帶來困難,這個問題可用L’H?pital 法則解決.首先,激波過程線S位于曲線M和N之間,所以其端點斜率也介于M和N的端點斜率之間;又因為曲線M和N為或者近似為拋物線,它們的端點斜率為有限值,所以S的端點斜率存在且為有限值,進而未定式f/g的極限也存在且為有限值,因此L’H?pital 法則適用.

將式 (16) 中的 ζ 乘到等號左邊,然后對等式兩端同時取極限——沿激波過程線向端點i(i=1,2)逼近,在這個過程中T可以視為V的函數(shù) (反之亦然).將f/g的分子分母同時對V求全導(dǎo),其極限的值保持不變

其中下標(biāo)V和T表示偏導(dǎo),fV=JPm-V,fT=cv,gV=2V-JPm,gT=R.值得一提的是,這里允許氣體是變比熱的,即不要求cv是常數(shù),它完全可以是T的函數(shù).

式 (20) 可整理成關(guān)于端點斜率 (dT/dV)i的一元二次方程

其具有一對相異實根

其中 γ 為比熱比,Pr為Prandtl 數(shù).不難看出這對根的符號相反.1 點是不穩(wěn)定結(jié)點,有無數(shù)對積分曲線從該點發(fā)出,其中一對在1 點處的斜率取正根;包括激波過程線在內(nèi)的其他積分曲線在1 點處相切,斜率取負根.2 點是鞍點,從該點發(fā)出和抵達該點的積分曲線各有一對,兩個根分別為這兩對積分曲線在2 點處的斜率,其中負根是所需要的激波過程線的終點斜率.

3.2 一般求解程序

綜上,激波內(nèi)部結(jié)構(gòu)逆向求解的一般程序可總結(jié)為以下4 步:

(1) 根據(jù)波前參數(shù)計算出波后參數(shù),尤其是波后馬赫數(shù)Ma2、波后速度V2以及波后溫度T2;

(2) 將波后參數(shù)代入式 (22) 中負根的表達式,獲得激波過程線的終點斜率 (dT/dV)2;

(3) 選一略大于V2的速度值V0,將其與其他所需參數(shù)代入式 (15),計算出相應(yīng)的T0;

(4) 以 (V0,T0) 為初值點,使用負步長 (即逆流向推進) 對式(12)～式(13) 構(gòu)成的系統(tǒng)進行數(shù)值積分;當(dāng)V與V1、T與T1的差值滿足精度要求時結(jié)束計算.

通過這4 步可獲得V和T的沿程分布,其他參數(shù)可基于它們進一步計算獲得而無需參與數(shù)值積分,例如將V和T的值代入式(5) 和式(11) 可求出ρ和p的值.

為了在精度和效率之間取得最佳平衡,建議在激波結(jié)構(gòu)問題中使用變步長的數(shù)值積分方法.這不僅是因為激波厚度隨激波強度的變化很大 (從無限弱激波的無窮大厚度到強激波僅數(shù)倍分子平均自由程的厚度),還因為速度梯度和溫度梯度在激波內(nèi)部不同位置差異巨大.自適應(yīng)步長可避免在參數(shù)變化劇烈的區(qū)域網(wǎng)格過疏或參數(shù)變化緩慢的地方網(wǎng)格過密.

3.3 算例驗證

為驗證逆向推進法的有效性,對單原子氣體中Ma1=1.01～100的正激波進行了求解,數(shù)值積分方法為步長自適應(yīng)的4 級Runge-Kutta 法.不失一般性地,假設(shè)氣體是定比熱的,且其輸運系數(shù)隨溫度的變化規(guī)律符合硬球分子模型的描述,即 μ 和κ 正比于T1/2.圖4 給出了代表性工況的沿程參數(shù)分布.

圖4 逆向推進法計算的單原子氣體激波結(jié)構(gòu)Fig.4 Shock structures in a monatomic gas calculated with the backward marching method

圖4 中展示的流動參數(shù)已經(jīng)過歸一化處理,x=0 的位置被統(tǒng)一規(guī)定在V=(V1+V2)/2 處,坐標(biāo)無量綱化所采用的分子平均自由程 λ 由硬球模型給定[1]

進一步地,將打靶法 (shooting method,SM) 與逆向推進法 (backward marching method,BMM) 計算的單原子氣體中Ma1=2 的正激波內(nèi)部的速度與溫度分布進行了對比,如圖5 所示.

圖5 打靶法與逆向推進法對激波結(jié)構(gòu)的計算結(jié)果的對比Fig.5 Comparison between the shock profiles calculated with the shooting method and the backward marching method

以上結(jié)果表明,逆向推進法可以正確可靠地求解激波內(nèi)部流動結(jié)構(gòu),明顯優(yōu)于打靶法.

4 結(jié)論

針對激波結(jié)構(gòu)問題,根據(jù)一維定常正激波控制方程的方向場畫出了激波內(nèi)部流動的相軌跡圖,基于其拓撲結(jié)構(gòu)分析了系統(tǒng)的動力學(xué)性質(zhì),闡述了積分方向?qū)η蠼庠搯栴}的重要性,并指出在激波結(jié)構(gòu)問題中,采用正向推進的傳統(tǒng)打靶法無法避免發(fā)散.為解決該問題,提出采用逆向推進計算方法,并給出與之配合使用的初值確定方法.在非常寬的馬赫數(shù)范圍 (1.01～100) 內(nèi)測試了逆向推進法,結(jié)果表明,該方法可以正確有效地求解激波結(jié)構(gòu)問題,且與傳統(tǒng)的打靶法相比,其具有以下優(yōu)勢:

(1) 計算無條件收斂;

(2) 對數(shù)值計算中不可避免的初值誤差、舍入誤差等不敏感,且誤差隨著計算的進行不斷降低;

(3) 無需對初值進行迭代,積分一次即獲得最終解,計算量小.

應(yīng)當(dāng)指出,黏性應(yīng)力與熱流密度的表達式雖然會影響討論的定量細節(jié),但不會使激波系統(tǒng)內(nèi)在的定性性質(zhì)發(fā)生根本變化,因此本文的求解方法及主要結(jié)論的適用范圍不限于只具有線性本構(gòu)的經(jīng)典Navier-Stokes 方程.原則上,只要本構(gòu)關(guān)系不違背熱力學(xué)第二定律,即等效黏度和等效導(dǎo)熱系數(shù)在任意位置處都非負,逆向推進法可以向各種具有Navier-Stokes 形式守恒方程 (式(1)～式(3)) 的復(fù)雜流體系統(tǒng)推廣,如改進的Navier-Stokes 方程、以Burnett 方程為代表的高階流體力學(xué)方程和Grad 矩方程等.對多原子氣體中考慮體積黏度或不透明氣體中考慮輻射傳熱的激波結(jié)構(gòu)問題,新方法亦適用.