隱式Runge－Kutta方法在非定?；瘜W(xué)非平衡流動(dòng)中的應(yīng)用

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

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

0 引 言

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

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

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

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

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

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

計(jì)算坐標(biāo)系下，三維非定?；瘜W(xué)反應(yīng)非平衡流動(dòng)的無量綱NS方程為：

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

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

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

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

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

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

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

2.1 隱式R－K方法正確性驗(yàn)證

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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

致謝：衷心感謝中國空氣動(dòng)力研究與發(fā)展中心楊順華博士、趙慧勇博士和王西耀博士提供的測試算例，以及對計(jì)算結(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ā)動(dòng)機(jī)并行數(shù)值研究［D］.四川綿陽：中國空氣動(dòng)力研究與發(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.