李 芳， 劉 鑫， 尹萬旺， 張 娟， 陸林生

（江南計算技術(shù)研究所，江蘇 無錫214083）

0 引 言

本文研究一種適合于求解非定?；瘜W非平衡流動的數(shù)值計算方法。該方法可應用于超聲速非平衡流中的一些復雜物理現(xiàn)象的研究，如非定常激波／邊界層相互作用、超聲速邊界層穩(wěn)定性和轉(zhuǎn)捩和超聲速湍流流動等。目前普遍應用的非定常計算方法分為兩大 類［1］：雙時間步方法［2］（pseudo－time subiteration）和物理時間迭代方法［2］（physical time subiteration）。雙時間步方法利用到兩個時間：一個偽時間，另一個是物理時間；而物理時間迭代方法只利用一個時間，即物理時間。

目前應用較多的是雙時間步方法。雙時間步方法只需在定常程序中作少量改動即可轉(zhuǎn)化為非定常程序，因此很多軟件采用雙時間步方法將定常程序和非定常程序合并在一起。但對于一些對時間精度要求較高的復雜非定常問題，雙時間步方法暴露出兩個問題：首先偽時間步迭代通常不易收斂，因而時間精度常常低于理論值2階，有時很難反映出非定常效應，其次雙時間步方法的偽時間步迭代次數(shù)較多，一般20幾步，甚至50步以上，計算效率不高，計算大型復雜非定常問題耗時很長。龍格－庫塔（Runge－Kutta）方法是一種精度較高的物理時間迭代方法，每個物理時間迭代中，X階精度只需計算X步，因而計算效率較高。目前使用較多的是顯式 R－K方法［7－9］，顯式R－K方法時間步長小，適合求解時間分辨率很高的非定常問題，如湍流直接數(shù)值模擬。

計算非定?；瘜W非平衡流動的主要困難在于控制方程的剛性［3－4］。由于化學反應特征時間遠小于宏觀非定常流動特征時間，因而化學反應非平衡源相會帶來剛性。有時邊界層內(nèi)的細網(wǎng)格也會使粘性張力和熱流源相帶來剛性。求解剛性問題，需用隱式格式代替顯式格式。目前普遍應用的隱式方法稱為半隱式法，即將控制方程的剛性項和非剛性項分開，分別采用隱式和顯式求解的方法。

Strang［5－6］在時間分裂方法中引入半隱式法，一部分時間推進采用隱式格式求解剛性部分，另一部分時間推進顯式求解對流項，并用到上一時間步的隱式結(jié)果。這種方法消除了很多燃燒問題中遇到的隱式方法、二階時間精度和強剛性問題的耦合，缺點是時間精度只有二階。Bussing［7］等采用半隱式 Mac－Cormach法計算可壓縮化學反應流，化學反應項隱式，流動項顯式。Engquist和Sjogreen［8］采用 TVD格式和ENO格式離散對流項，采用半隱式法將源相合并到Runge－Kutta時間推進中求解了爆炸波。Zhong［9－11］以泰勒展開為基礎(chǔ)，非剛性項顯式處理，剛性項隱式處理，推導了一種滿足高精度和強穩(wěn)定性條件的半隱式Runge－Kutta方法，但這種方法系數(shù)繁多，計算過程比較復雜。

本文提出一種新的隱式Runge－Kutta方法，中間步的系數(shù)與經(jīng)典顯式Runge－Kutta方法相同，隱式體現(xiàn)在右端項的計算。這種隱式R－K的方法的時間精度和同階顯式R－K方法相當，時間步長可以取的較大，因而計算效率較高，同時計算公式簡單，實現(xiàn)方便，可以很容易地從定常程序轉(zhuǎn)化為非定常程序。

1 求解NS方程組的隱式龍格－庫塔法

計算坐標系下，三維非定常化學反應非平衡流動的無量綱NS方程為：

R－K法本是求解常微分方程的一類高精度方法，為在偏微分方程（1）中應用R－K法，需將方程（1）變?yōu)槿缦滦纬桑?/p>

采用隱式方法，需分別對無粘對流項、粘性項和源項進行線化處理，線化方法與定常程序類似，對流項采用Jacobian矩陣分裂法線化，粘性項采用近似特征值法線化，源項采用對角化法線化［12］，線化后方程（2）的右端項可表示為：

式中，α＝0為空間顯式格式，α＝1為空間隱式格式；βD＝0時，化學源項不做隱式線化，βD≠0時，化學源項隱式線化；ρ（）、ρ（）、ρ（）分別為矩陣、、的譜半徑

采用δQ的形式，2階精度隱式R－K法可表示為：

采用δQ的形式，3階精度隱式R－K法可表示為：

由于采用相同的線化方法，因此隱式R－K方法的右端項計算方法與定常程序相同，可以方便地從定常程序轉(zhuǎn)化為非定常程序。

2 算例計算結(jié)果與分析

2.1 隱式R－K方法正確性驗證

為驗證算法正確性，選取了六個典型定常算例：湍流平板、二維燃燒室、超聲速壓縮拐角、三維進氣道、球頭激波誘導燃燒、非對稱交叉激波。分別采用隱式LU－SGS方法［13］和2階、3階隱式 R－K方法進行數(shù)值模擬，并將計算結(jié)果與實驗結(jié)果進行比較，結(jié)果表明三種隱式方法計算結(jié)果一致，且與試驗吻合。由于篇幅有限，這里僅給出球頭激波誘導燃燒計算5000步時的溫度場，如圖1所示。從圖中可看出，計算5000步時，三種時間推進方法的溫度場均已收斂，激波形狀幾乎完全相同，與試驗值吻合；3階隱式RK法的CFL數(shù)可取到60，2階隱式R－K法和隱式LU－SGS方法的CFL數(shù)較小，可見3階隱式R－K法穩(wěn)定性較好。

2.2 隱式R－K方法在非定常問題中的效率

非定常問題極其復雜，且計算量巨大，即使采用成千上萬個CPU并行計算，模擬一個大型實際問題也需要幾天甚至十幾天時間，因此求解非定常問題，除了考慮數(shù)值格式的穩(wěn)定性和精度以外，計算效率也非常重要。

20世紀60～70年代，Lehr［14］通過試驗發(fā)現(xiàn)球頭激波誘導燃燒存在非定?，F(xiàn)象，之后Jeong［15］等針對該現(xiàn)象開展了數(shù)值模擬研究，并取得較滿意的結(jié)果。由于化學反應非?？欤⑶以诤芏痰木嚯x內(nèi)伴隨著大量能量釋放，因此對數(shù)值格式的健壯性和精度要求較高。本文通過對球頭激波點火試驗進行數(shù)值模擬來校驗隱式R－K方法模擬非定?，F(xiàn)象的精度和效率。

來流條件為馬赫數(shù)4.48，壓強42663.16Pa，溫度293K，來流質(zhì)量分數(shù)為CH2＝0.0283，CO2＝0.227，CN2＝0.74507。計算采用二維軸對稱歐拉方程，無粘通量采用S－W分裂，化學反應為9組分19方程的反應模型，使用256個CPU并行計算。

時間推進分別采用2階隱式、3階隱式R－K法和2階雙時間步方法（偽時間迭代為20步），當時間步長均為10－9s時，點火過程中駐點壓強隨時間的變化曲線如圖2所示?？梢钥闯?，三種方法的駐點壓強均在點火后出現(xiàn)了一定頻率的振蕩，表明流場中產(chǎn)生了爆轟波現(xiàn)象，與實驗結(jié)果吻合。三種方法的振蕩頻率和振幅接近，均能正確描述該非定常問題。

采用不同的時間步長，3階隱式R－K法和2階雙時間步方法得到的振蕩頻率如表1所示。從表1可看出，隨著時間步長的縮小，兩種方法得到的振蕩頻率均增加，并逐漸靠近參考值：426kHz［15］。

表1 不同時間步長得到的振蕩頻率（kHz）Table 1 Vibration frequency with different time step（unit：kHz）

2階隱式、3階隱式R－K法和2階雙時間步方法均為空間導數(shù)離散的隱式格式，時間步長可以取得相當。一個物理時間步長內(nèi)，2階隱式R－K法計算2次，3階隱式R－K法計算3次，雙時間步方法計算20（偽時間迭代數(shù)）次，因此取相同的物理步長和總計算步數(shù)，雙時間步方法的總計算時間分別為2階、3階R－K方法的10倍和6.67倍，由此可見，隱式R－K方法在計算效率方面有絕對優(yōu)勢。

3 結(jié) 論

本文以經(jīng)典的顯式Runge－Kutta方法為基礎(chǔ)，推導了一種求解NS方程組的隱式R－K方法。該方法通過右端項隱式處理，解決了化學反應非平衡流動帶來的剛性問題，同時計算公式簡單，實現(xiàn)方便，可以很容易地從定常程序轉(zhuǎn)化為非定常程序。最后通過定常算例和非定常算例對隱式R－K方法進行了驗算，分析表明該方法計算結(jié)果正確，精度較高；能夠正確求解復雜非定?；瘜W反應非平衡流動，且算法健壯性好、精度高、效率高。

致謝：衷心感謝中國空氣動力研究與發(fā)展中心楊順華博士、趙慧勇博士和王西耀博士提供的測試算例，以及對計算結(jié)果的分析與處理。

［1］RUMSEY C L，SANETRIK M D，BIEDRON R T，et al.Efficiency and accuracy of time－accurate turbulent Navier－Stokes compuations［R］.AIAA paper 95－1835，1995.

［2］PULLIAM T H，Time accuracy and the use of implicit methods［R］.AIAA paper 93－3360，1993.

［3］ZHONG X.Direct numerical simulation of hypersonic boundary layer transition over blunt leading edges，Part I：a new numerical method and validation［R］.AIAA Paper 97－0755，1997.

［4］ZHONG X.Direct numerical simulation of hypersonic boundary layer transition over blunt leading edges，Part II：receptivity to sound［R］.AIAA Paper 97－0756，1997.

［5］STRANG G.On the construction and comparison of difference schemes［J］.SIAMJournalofNumerical Analysis，1968，5：506－517.

［6］LEVEQUE R J，YEE H C.A study of numerical methods for hyperbolic conservation laws with stiff source terms［J］.JournalofComputationalPhysics，1990，68：187－210.

［7］BUSSING T R A，MURMAN E M.A finite volume method for the calculation of compressible chemically reacting flows［R］.AIAA Paper 85－0296，1985.

［8］ENQUIST B，SJOGREEN B.Robust difference approximations of stiff inviscid detonation waves［R］.Report 91－03，CAM，Department of Mathematics，University of California，1991.

［9］ZHONG X.Additive semi－implicit Runge－Kutta methods for computing high－speed nonequilibrium reactive flows［J］.JournalofComputationalPhysics，1996，128：19－31.

［10］JACK JAI－ICK YOH.ZHONG X.Low－storage semiimplicit Runge－Kutta methods for reactive flow computations［R］.AIAA Paper 98－0130，1998.

［11］ZHONG X.New high－order semi－implicit Runge－Kutta Schemes for computing transient nonequilibrium hypersonic flow［R］.AIAA Paper 95－2007，1995.

［12］趙惠勇.超燃沖壓整體發(fā)動機并行數(shù)值研究［D］.四川綿陽：中國空氣動力研究與發(fā)展中心，2005.

［13］JAMESON A，YOON S.LU implicit scheme with multiple grid for the Euler equations［R］.AIAA paper 86－0105，1986.

［14］LEHR H F.Experimenta on shock－induced combustion［J］.AstronauticaActa，1972，17（4－5）：195－203.

［15］CHOI J Y，IN－SEUCK JEUNG，YOUNGBIN YOON.Computational fluid dynamics algorithms for unsteady shock－induced combustion，Part I：validation ［R］.AIAA paper 98－3217.