Godunov型顯式大時(shí)間步長(zhǎng)格式研究進(jìn)展

錢(qián)戰(zhàn)森

1. 中國(guó)航空工業(yè)空氣動(dòng)力研究院，沈陽(yáng) 110034 2. 高速高雷諾數(shù)氣動(dòng)力航空科技重點(diǎn)實(shí)驗(yàn)室，沈陽(yáng) 110034

自20世紀(jì)七八十年代高性能計(jì)算機(jī)出現(xiàn)以來(lái)，數(shù)值計(jì)算方法、網(wǎng)格生成技術(shù)、并行算法及圖像處理等得到了蓬勃發(fā)展，數(shù)值模擬的應(yīng)用范圍不斷拓寬，深度和復(fù)雜程度不斷加大，解決實(shí)際問(wèn)題的能力也不斷提高，使計(jì)算流體力學(xué)(CFD)逐步成為與理論分析和實(shí)驗(yàn)相并列的研究流體力學(xué)問(wèn)題的三大手段之一[1-12]。計(jì)算流體力學(xué)的應(yīng)用遍及航空航天、天氣預(yù)報(bào)、海洋與船舶、汽車(chē)和列車(chē)外型的空氣動(dòng)力學(xué)優(yōu)化、風(fēng)扇設(shè)計(jì)、降噪技術(shù)、換熱器設(shè)計(jì)等基礎(chǔ)工業(yè)的各個(gè)方面，以及生物工程、土木工程、水利工程、環(huán)境保護(hù)及風(fēng)沙治理等工程應(yīng)用領(lǐng)域[13-16]。特別是在航空航天領(lǐng)域，一方面，計(jì)算流體力學(xué)的應(yīng)用可以拓寬和深化流體力學(xué)的基礎(chǔ)研究，揭示經(jīng)典方法無(wú)法獲取的清晰圖像，例如飛機(jī)在大攻角飛行或機(jī)動(dòng)飛行的情況下，傳統(tǒng)的風(fēng)洞試驗(yàn)手段一般難以取得相關(guān)復(fù)雜現(xiàn)象的清晰信息；另一方面，相對(duì)于風(fēng)洞試驗(yàn)存在的試驗(yàn)周期長(zhǎng)、價(jià)格昂貴的缺點(diǎn)，計(jì)算流體力學(xué)快捷、低成本的優(yōu)點(diǎn)更為突顯。

然而航空航天工程問(wèn)題數(shù)值模擬的計(jì)算量是十分巨大的，特別是在大規(guī)模工程計(jì)算中。例如對(duì)復(fù)雜飛行器氣動(dòng)特性數(shù)值模擬，需要在邊界層和分離區(qū)等流動(dòng)結(jié)構(gòu)復(fù)雜的區(qū)域內(nèi)密布網(wǎng)格,而在三維計(jì)算中，網(wǎng)格尺度每增加一倍的密度，計(jì)算量便要增加約16倍。在湍流問(wèn)題的直接數(shù)值模擬中，需要真實(shí)雷諾數(shù)下的計(jì)算，其計(jì)算量將和雷諾數(shù)的9/4次方成比例。在飛行器氣動(dòng)布局優(yōu)化過(guò)程中，更是需要計(jì)算成百上千個(gè)樣本狀態(tài)，才能獲得最終的優(yōu)化氣動(dòng)布局形式。根據(jù)目前工程使用經(jīng)驗(yàn)，計(jì)算速度仍是制約當(dāng)前計(jì)算流體力學(xué)應(yīng)用的一大瓶頸，在計(jì)算機(jī)性能和資源相對(duì)有限的情況下，建立高效的數(shù)值算法以提高計(jì)算速度對(duì)開(kāi)展大規(guī)模科學(xué)與工程計(jì)算具有重大的現(xiàn)實(shí)意義。

提高計(jì)算速度的主要方法包括優(yōu)化程序設(shè)計(jì)、并行計(jì)算、區(qū)域分解、改善數(shù)值格式等。優(yōu)化程序設(shè)計(jì)主要考慮在空間循環(huán)中減少耗費(fèi)CPU時(shí)間的運(yùn)算過(guò)程，但對(duì)已開(kāi)發(fā)的格式，通過(guò)優(yōu)化程序設(shè)計(jì)所能節(jié)約的計(jì)算量是有限的。隨著計(jì)算機(jī)資源的不斷發(fā)展，并行計(jì)算已成為助力復(fù)雜飛行器數(shù)值模擬分析的重要手段，但在一定時(shí)間內(nèi)計(jì)算機(jī)資源總是有限的，模擬規(guī)模和數(shù)量往往受限。區(qū)域分解方法[17]可以根據(jù)局部流場(chǎng)的主控尺度特征對(duì)不同區(qū)域采用不同主控方程進(jìn)行模擬，也可在一定程度上明顯提高數(shù)值模擬的效率。改善數(shù)值格式，即從數(shù)值格式構(gòu)造方面提高計(jì)算速度的方法，是提高數(shù)值模擬效率最直接的途徑。格式速度每提高一倍，相當(dāng)于計(jì)算機(jī)CPU主頻提高一倍或并行計(jì)算機(jī)數(shù)目減少一半，這是相當(dāng)可觀的收益。該方法也可以和前述方法耦合使用，獲得更大的收益。因此，提高流場(chǎng)解算器數(shù)值格式的計(jì)算速度對(duì)于大規(guī)模計(jì)算具有特殊的意義。

放寬顯式計(jì)算格式的時(shí)間步長(zhǎng)限制，提高計(jì)算的CFL(Courant-Friedrichs-Lewy)數(shù)，從而提高計(jì)算效率，是一種有意義的探索，即發(fā)展所謂的顯式大時(shí)間步長(zhǎng)格式。這里CFL數(shù)的定義一般可寫(xiě)為CFL=aΔt/Δx，或Δt=CFLΔx/a，其中a為雙曲型守恒律方程的Jacobi矩陣的最大特征值(通常表征特征波的最大傳播速度)，Δt為離散的時(shí)間步長(zhǎng)，Δx為離散的空間網(wǎng)格尺度。對(duì)于顯式格式而言，一般情況下要滿(mǎn)足CFL≤1(稱(chēng)為CFL條件)，這導(dǎo)致時(shí)間步長(zhǎng)受限，從而計(jì)算效率受限。這里的顯式大時(shí)間步長(zhǎng)(Large Time Step， LTS)格式是指一種CFL數(shù)突破1限制的全離散的、單步的、顯式的計(jì)算格式，仍滿(mǎn)足CFL條件，但CFL條件需要推廣，這將在下文中詳述。

本文首先簡(jiǎn)要介紹了顯式大時(shí)間步長(zhǎng)格式的概念、分類(lèi)和優(yōu)勢(shì)，然后重點(diǎn)回顧了針對(duì)Godunov型顯式大時(shí)間步長(zhǎng)格式的發(fā)展，并給出了該類(lèi)格式今后的發(fā)展方向。本項(xiàng)工作旨在為促進(jìn)大時(shí)間步長(zhǎng)格式領(lǐng)域的交流提供參考。

1 顯式大時(shí)間步長(zhǎng)格式簡(jiǎn)介

1.1 數(shù)值格式的發(fā)展概述

1) 雙曲型守恒律的時(shí)間相關(guān)方法

在航空航天計(jì)算流體力學(xué)中廣泛采用時(shí)間相關(guān)方法[18-20]數(shù)值求解Euler或Navier-Stokes(N-S)方程。按照時(shí)間推進(jìn)方式，求解流體力學(xué)方程的時(shí)間相關(guān)格式可以歸納為顯式格式和隱式格式；按照離散方式可以分為全離散格式和半離散格式。

全離散格式指的是在離散控制方程時(shí)同時(shí)考慮時(shí)間項(xiàng)和空間項(xiàng)的導(dǎo)數(shù)近似，以取得時(shí)間和空間方向所具有的設(shè)定精度。全離散格式基本都是顯式的，一般可通過(guò)原偏微分方程以空間導(dǎo)數(shù)來(lái)表達(dá)其時(shí)間導(dǎo)數(shù)。對(duì)于線(xiàn)性常系數(shù)方程，更高階的場(chǎng)變量時(shí)間導(dǎo)數(shù)離散可容易地以其等價(jià)空間導(dǎo)數(shù)表達(dá)，但對(duì)于非線(xiàn)性方程，該變換實(shí)現(xiàn)較為復(fù)雜，甚至無(wú)法顯式表達(dá)。

半離散格式是指對(duì)控制方程進(jìn)行離散時(shí)將時(shí)間方向和空間方向離散解耦的方法。該思路可通過(guò)分別離散空間與時(shí)間導(dǎo)數(shù)較容易地導(dǎo)出設(shè)定精度的數(shù)值格式。通常可先對(duì)空間方向作離散，保持時(shí)間方向的連續(xù)性，導(dǎo)出關(guān)于時(shí)間變量的常微分方程，稱(chēng)為原控制方程的半離散方程。而導(dǎo)出的常微分方程的離散可以借助已有的標(biāo)準(zhǔn)方法給出。

顯式格式在工程實(shí)際中，特別是大規(guī)模計(jì)算中被廣泛應(yīng)用。該類(lèi)格式的優(yōu)點(diǎn)是每推進(jìn)一個(gè)時(shí)間步長(zhǎng)計(jì)算量和存儲(chǔ)量較少，程序簡(jiǎn)單，缺點(diǎn)是時(shí)間步長(zhǎng)受到穩(wěn)定性條件的限制，導(dǎo)致總體計(jì)算步數(shù)往往較多。目前已發(fā)展了大量可以求解常微分方程組初值問(wèn)題的顯式方法[21-22]，如Euler折線(xiàn)法、Adams外插法、Adams內(nèi)插法、Milne方法、Hamming算法、Adams三步四階預(yù)估-校正算法及Runge-Kutta方法等。航空航天流體力學(xué)計(jì)算中使用較多的顯式格式為Runge-Kutta方法，其時(shí)間精度可達(dá)到二或更高階，時(shí)間步長(zhǎng)可以適當(dāng)放寬。但Runge-Kutta方法屬于多步法，所有多步法均以增加計(jì)算量為代價(jià)，每多一層子循環(huán)，計(jì)算量便增加一倍，不符合高效計(jì)算的要求。

隱式格式在模型方程的線(xiàn)性穩(wěn)定性分析中一般是呈無(wú)條件穩(wěn)定的，在數(shù)值求解實(shí)際問(wèn)題時(shí)也可取較大的時(shí)間步長(zhǎng)，目前也得到了較為廣泛的應(yīng)用。航空航天流體力學(xué)計(jì)算中具代表性的隱式方法有Beam-Warming隱式格式[23]、MacCormack隱式方法[24]、近似因式分解方法[25]、改進(jìn)的對(duì)角化近似因式分解方法[26]、LU-SGS方法[27]以及GMRES方法[28-29]等。這些隱式格式在定常流場(chǎng)計(jì)算中得到了廣泛應(yīng)用。然而，一般提高計(jì)算效率的近似處理均會(huì)導(dǎo)致迭代過(guò)程中時(shí)間精度的降階,而低階精度特別是一階精度的時(shí)間格式，不適用于非定常流的模擬。為了克服這一難題，20世紀(jì)90年代發(fā)展了模擬非定常問(wèn)題的雙時(shí)間步迭代方法[30]。但是目前工程實(shí)踐經(jīng)驗(yàn)表明，雙時(shí)間步迭代方法由于其內(nèi)迭代常常無(wú)法在短時(shí)間內(nèi)收斂或者甚至不收斂，需要通過(guò)設(shè)定內(nèi)迭代的次數(shù)來(lái)降低計(jì)算量，該措施必然引入附加誤差，甚至可能導(dǎo)致非定常流動(dòng)的模擬失真。

2) 激波捕捉格式的發(fā)展

航空航天領(lǐng)域廣泛關(guān)注可壓縮流動(dòng)的求解問(wèn)題，可壓縮流動(dòng)中的激波模擬給數(shù)值格式帶來(lái)了強(qiáng)非線(xiàn)性的重要特征，常用的激波模擬方法包括激波裝配法和激波捕捉法。相對(duì)于激波裝配法，激波捕捉法因其操作簡(jiǎn)單、實(shí)用而成為計(jì)算流體力學(xué)廣泛采用的方法。

Godunov于1959年提出的間斷分解法被認(rèn)為是最早的激波捕捉格式[31]，即通過(guò)Riemann問(wèn)題間斷分解的計(jì)算來(lái)構(gòu)造差分格式，這類(lèi)格式統(tǒng)稱(chēng)作Godunov型格式，這個(gè)思路直到今天依然是CFD差分格式研究的熱點(diǎn)，并且一直處于不斷完善和發(fā)展之中。近30年來(lái)，Godunov型格式的研究和發(fā)展主要有兩個(gè)方向：一是對(duì)格式的高階推廣，提高格式的精度以提高對(duì)各種間斷的分辨率，代表性工作包括Boris和Book[32]、Harten等[33-34]通過(guò)對(duì)一階差分格式的數(shù)值通量進(jìn)行修正而獲得二階精度的格式，Jameson等[35]在中心差分格式(JST格式)基礎(chǔ)上采用人工黏性法獲得二階精度的格式，van Leer[36-39]及Collela和Woodward[40]通過(guò)對(duì)初始時(shí)刻的解采用分段線(xiàn)性或拋物線(xiàn)重構(gòu)函數(shù)代替原Gudonov格式中的分段常數(shù)重構(gòu)以獲得二階或三階精度的格式，Harten等[41]提出的基本無(wú)震蕩(ENO)格式及Liu[42]、Shu[43]等提出并改進(jìn)的高效WENO格式，通過(guò)采用隱式差分可以在較小的模板點(diǎn)上獲得更高的插值精度的緊致格式及其改進(jìn)形式[44-53]；二是對(duì)近似Riemann解的研究，用各種近似Riemann求解器來(lái)代替原始Godunov格式采用的Riemann問(wèn)題精確解以減少計(jì)量，代表性工作包括Roe[54]、Osher[55-56]、Harten[57]、Toro[58]等分別發(fā)展出了各種通量差分分裂形式的近似Riemann解算器；Steger[59]、van Leer[60]、Liou[61-65]、Jameson[66]、張涵信[67-68]等分別發(fā)展出了各種矢通量分裂形式的近似Riemann解算器；Lax[69]、MacCormack[70]、Jameson[35]、Tadmor[71-72]、Levy[73-74]等發(fā)展的各種以Lax-Friedrichs格式[75]為基礎(chǔ)的中心型Riemann解算器。總體來(lái)看，對(duì)近似Riemann解的研究近年來(lái)活躍度不如高階推廣高，其中較有代表性的研究工作是針對(duì)多維問(wèn)題發(fā)展的混合型近似Riemann解[76-79]。

當(dāng)前廣泛采用的各種CFD求解方法，都可以看作是Godunov格式的高階推廣或間斷分解方法的改進(jìn)，即都屬于Godunov型格式，因而都需要聯(lián)合采用高階精度的插值方法和類(lèi)似Godunov格式的精確或近似的Riemann問(wèn)題解的數(shù)值通量，兩者共同組成了可壓縮流動(dòng)的激波捕捉格式。

1.2 顯式大時(shí)間步長(zhǎng)格式的概念

航空航天領(lǐng)域空氣動(dòng)力學(xué)數(shù)值模擬的主要控制方程為N-S方程，一個(gè)完整的N-S方程可以通過(guò)算子分裂[18]分解成無(wú)黏部分、有黏部分以及源項(xiàng)部分。其中無(wú)黏部分即Euler方程，屬雙曲型；有黏部分屬拋物型；源項(xiàng)部分可以是線(xiàn)性或非線(xiàn)性代數(shù)函數(shù)、常微分或偏微分形式。拋物型方程由于不存在占主導(dǎo)地位的信息傳播方向，最佳離散格式為中心差分格式。常微分方程的離散可以借助已有的線(xiàn)性多步法等標(biāo)準(zhǔn)方法處理。

目前CFD領(lǐng)域有關(guān)差分格式的研究主要集中在Euler方程相關(guān)項(xiàng)，其為典型的雙曲型守恒律方程，顯式格式和隱式格式均有所發(fā)展，增大數(shù)值格式的時(shí)間步長(zhǎng)是提高格式計(jì)算速度最有效的途徑之一。但是對(duì)于顯式格式而言，一般情況下時(shí)間步長(zhǎng)要受到穩(wěn)定性要求(CFL條件)的限制。隱式格式因其理論上沒(méi)有時(shí)間步長(zhǎng)限制，在實(shí)際模擬中也可以取較大的時(shí)間步長(zhǎng)，在工程計(jì)算領(lǐng)域得到了更廣泛的應(yīng)用。然而，Euler方程由于其雙曲型性質(zhì)，即解的依賴(lài)域是有限的，理論上不宜采用其數(shù)值解的依賴(lài)域，往往是全場(chǎng)相關(guān)的隱式格式求解。

發(fā)展顯式大時(shí)間步長(zhǎng)格式，突破顯式格式的時(shí)間步長(zhǎng)限制是另一種提高計(jì)算效率的途徑。如前文所述，這里的顯式大時(shí)間步長(zhǎng)格式其CFL數(shù)可突破1的限制，但仍保持單步、顯式、全離散的特征，并且滿(mǎn)足穩(wěn)定性和相容性等基本特性。本文重點(diǎn)關(guān)注該類(lèi)格式的研究進(jìn)展。

1.3 大時(shí)間步長(zhǎng)格式的分類(lèi)

根據(jù)目前就雙曲型守恒律方程的求解格式開(kāi)展的研究工作來(lái)看，顯式大時(shí)間步長(zhǎng)格式按照構(gòu)造方法可以分為兩類(lèi)：第1類(lèi)為幾何型構(gòu)造，第2類(lèi)為代數(shù)型構(gòu)造。第1類(lèi)主要基于Godunov型格式構(gòu)造所得，其中代表性的為L(zhǎng)eVeque提出的行波法[80-82]，針對(duì)該格式的不足，Guinot[83]提出了時(shí)間線(xiàn)插值法，本文作者[84-85]提出了膨脹波多波近似法和近距雙波干擾法，Dong和Liu[86]提出了膨脹波單元分解法，分別對(duì)其進(jìn)行了改進(jìn)，增強(qiáng)了格式的分辨率和魯棒性；第2類(lèi)主要采用代數(shù)運(yùn)算來(lái)實(shí)現(xiàn)，其中具有代表性的為Harten[87]提出的TVD大時(shí)間步長(zhǎng)格式和董海濤、李椿萱[88]提出的熵條件大時(shí)間步長(zhǎng)格式。Harten通過(guò)算子合并將原TVD格式構(gòu)造為單步大時(shí)間格式，本文作者[89]對(duì)其作了改進(jìn)并將其推廣至多維流動(dòng)的求解。董海濤和李椿萱[88]利用Hamilton-Jacobi方程精確解構(gòu)造了雙曲型守恒律方程的大時(shí)間步長(zhǎng)格式。第2類(lèi)方法與第1類(lèi)方法得到的格式性質(zhì)差異較大，本文主要關(guān)注第1類(lèi)大時(shí)間步長(zhǎng)格式，即Godunov型大時(shí)間步長(zhǎng)格式。

1.4 大時(shí)間步長(zhǎng)格式的優(yōu)勢(shì)

與常規(guī)顯式格式和隱式格式相比，Godunov型顯式大時(shí)間步長(zhǎng)格式具有以下幾個(gè)方面的突出優(yōu)點(diǎn)：

1) 定常問(wèn)題計(jì)算效率高

時(shí)間相關(guān)法是求解雙曲型守恒律方程的典型方法，但時(shí)間推進(jìn)的步長(zhǎng)直接決定著流場(chǎng)信息的推進(jìn)速度。常規(guī)顯式格式的CFL數(shù)不能超過(guò)1，大時(shí)間步長(zhǎng)格式則可以突破這一限制，實(shí)際計(jì)算結(jié)果[82-84]表明，后者計(jì)算效率也確實(shí)得到了充分的提高。隱式格式的CFL數(shù)雖然理論上不受限，但在實(shí)際復(fù)雜問(wèn)題應(yīng)用中往往不能取得過(guò)大。

2) 可直接應(yīng)用于非定常問(wèn)題計(jì)算

顯式大時(shí)間步長(zhǎng)格式可直接應(yīng)用于非定常問(wèn)題計(jì)算，且因其屬于全離散格式，具有時(shí)空一致的代數(shù)精度。與常規(guī)顯式格式和以Runge-Kutta方法為代表的多步法相比計(jì)算效率明顯提高；與隱式格式相比，其不需要雙時(shí)間步迭代，具有精度更高、魯棒性更強(qiáng)的優(yōu)勢(shì)。

3) 對(duì)間斷分辨率高

多個(gè)研究結(jié)果[82-84]表明，在同等代數(shù)精度下，隨著CFL數(shù)的增大，Godunov型大時(shí)間步長(zhǎng)格式分辨率逐漸提高。特別是對(duì)于激波等強(qiáng)非線(xiàn)性間斷和接觸間斷等線(xiàn)性或擬線(xiàn)性間斷均具有很高的分辨率，適合于空氣動(dòng)力學(xué)的激波干擾、混合層、渦結(jié)構(gòu)、聲波等復(fù)雜流動(dòng)現(xiàn)象的求解。

4) 大規(guī)模并行計(jì)算可擴(kuò)展性好

大時(shí)間步長(zhǎng)格式因是顯式格式，對(duì)計(jì)算機(jī)內(nèi)存要求較低。同時(shí)因其信息依賴(lài)域有限，故而在分區(qū)并行中具有較大優(yōu)勢(shì)，尤其適合于當(dāng)前CFD領(lǐng)域普遍采用的幾何分區(qū)并行計(jì)算，一方面不同幾何區(qū)域界面之間的數(shù)據(jù)通訊量較小，另一方面收斂效率基本不受幾何分區(qū)的影響。故大時(shí)間步長(zhǎng)格式有著非常優(yōu)越的大規(guī)模并行計(jì)算可擴(kuò)展性。與之相比較，隱式格式則在大規(guī)模并行計(jì)算方面的擴(kuò)展性稍欠。由于其計(jì)算往往是全場(chǎng)相關(guān)的，如果要嚴(yán)格按照原始方式推進(jìn)，則需要在不同幾何分區(qū)之間進(jìn)行大量的數(shù)據(jù)傳遞，這將占用龐大的時(shí)間；若簡(jiǎn)單地在每一個(gè)分區(qū)內(nèi)單獨(dú)使用諸如LU-SGS等形式的隱式推進(jìn)格式，則將導(dǎo)致算法的魯棒性和收斂效率大大降低。

5) 可與現(xiàn)有加速收斂的算法良好兼容

長(zhǎng)期以來(lái)已發(fā)展了不少行之有效的基于時(shí)間相關(guān)法求解定常問(wèn)題的快速收斂算法，如當(dāng)?shù)貢r(shí)間步長(zhǎng)[90]、殘值光順[91-93]、焓阻尼技術(shù)[94]和多重網(wǎng)格技術(shù)[30，95]等。大時(shí)間步長(zhǎng)格式因其顯式特點(diǎn)，可與這些加速收斂的算法很好地兼容，大大改善對(duì)于定常問(wèn)題的求解效率，這對(duì)于大時(shí)間步長(zhǎng)格式的工程應(yīng)用推廣具有十分重要的意義。

2 Godunov型大時(shí)間步長(zhǎng)格式的提出

2.1 Godunov型格式

考慮一維雙曲型守恒律組的初值問(wèn)題：

(1)

圖1 局部Godunov平均示意圖Fig.1 Schematic of local Godunov averaging

(2)

由于函數(shù)u(x,tn+1)在區(qū)間[xj-1/2,xj]和[xj,xj+1/2]上由不同的間斷分解產(chǎn)生，需要進(jìn)行分段積分。這種方法的實(shí)現(xiàn)較煩瑣，且時(shí)間步長(zhǎng)受到CFL≤0.5的嚴(yán)格限制，以避免左右兩相鄰間斷波發(fā)生相互干擾。因此，在實(shí)際計(jì)算中很少采用此方法。

2) 第2種方法 采用守恒型格式：

(3)

式中：

(4)

顯然，以上兩種方法的CFL數(shù)都要受到限制，即不能超過(guò)1，從而給時(shí)間推進(jìn)步長(zhǎng)帶來(lái)了嚴(yán)格的限制。

2.2 Godunov型大時(shí)間步長(zhǎng)格式

沿著原始Godunov型格式的思想，LeVeque最早提出了一種用于計(jì)算雙曲型守恒律的大時(shí)間步長(zhǎng)格式[80-82]。該格式的CFL數(shù)可突破1的限制，其核心思想是：沿用分段常數(shù)重構(gòu)函數(shù)來(lái)代替網(wǎng)格點(diǎn)上的離散值，形成一系列的Riemann問(wèn)題。在對(duì)該系列的Riemann問(wèn)題作間斷分解后，將得到一系列的波(包括激波、膨脹波、接觸間斷等)。于是可對(duì)其進(jìn)行時(shí)間積分求出數(shù)值解。當(dāng)CFL≤0.5時(shí)，波系間顯然不會(huì)發(fā)生相互干擾，可直接采用原始Godunov型格式的第1種方法或第2種方法求解。當(dāng)CFL≤1.0時(shí)，采用第2種平均方法也不致出現(xiàn)實(shí)質(zhì)性困難。為了使CFL數(shù)進(jìn)一步提高，LeVeque對(duì)波系之間的相互干擾作如下的線(xiàn)性假設(shè)：假設(shè)兩束波相交時(shí)互不影響，即兩束波的傳播速度和強(qiáng)度均保持不變，也無(wú)新的波產(chǎn)生，只是簡(jiǎn)單的相互穿越。

由于大時(shí)間步長(zhǎng)格式允許CFL>1，因此需面對(duì)間斷分解后的強(qiáng)波與強(qiáng)波之間的干擾問(wèn)題：強(qiáng)波相交問(wèn)題。對(duì)于線(xiàn)性方程組，異族強(qiáng)波相互穿越并且穿越后波速與波強(qiáng)不變(直線(xiàn)穿越)，同族強(qiáng)波則變?yōu)橄嗷テ叫械闹本€(xiàn)，故LeVeque的上述假設(shè)在常系數(shù)線(xiàn)性雙曲型守恒律是精確成立的；對(duì)于非線(xiàn)性方程，異族強(qiáng)波相互穿越并且穿越后波速與波強(qiáng)會(huì)有改變，同族強(qiáng)波一般不易相交，如相交則合并成一個(gè)強(qiáng)波。可見(jiàn)對(duì)非線(xiàn)性方程，LeVeque假設(shè)會(huì)帶來(lái)一定誤差，但是實(shí)際計(jì)算結(jié)果表明，其仍可獲得較好的計(jì)算結(jié)果。

(5)

則式(5)的系列分片常數(shù)Riemann問(wèn)題

(6)

的解可表達(dá)為

(7)

首先考慮單束波的簡(jiǎn)單問(wèn)題，即由標(biāo)量雙曲型守恒律主控的問(wèn)題，如圖2所示[84]。設(shè)在點(diǎn)(xj+1/2,tn)處有Δu=uj-uj+1的間斷，sa為該間斷的傳播速度。當(dāng)sa>0時(shí)，對(duì)于x∈[xj+1/2,xj+1/2+saΔt]，un(x,tn+1)的值相對(duì)于un(x,tn)存在Δu的變化，即un+1(x,tn+1)=un(x,tn)+Δu。若該束波的傳播已越過(guò)[xj+1/2,xj+3/2]的整個(gè)區(qū)間，則對(duì)區(qū)間[xj+1/2,xj+3/2]內(nèi)的每一點(diǎn)x，un(x,tn+1)在時(shí)間段[tn,tn+1]均增加了Δu，于是整個(gè)區(qū)間的平均值也增加了Δu。若在時(shí)間段[tn,tn+1]內(nèi)波束未穿過(guò)整個(gè)區(qū)間，僅到達(dá)區(qū)間[xj+1/2,xj+3/2]內(nèi)某一點(diǎn)ε(xj+1/2<ε

圖2 單波束在網(wǎng)格上的傳播[84]Fig.2 A single wave propagation on grid[84]

(8)

同樣，對(duì)于sa<0的情況也可作類(lèi)似的處理。

采用顯式格式數(shù)值求解該問(wèn)題，由穩(wěn)定性條件(CFL條件)知微分方程數(shù)值解的依賴(lài)域必須包含精確解的依賴(lài)域，因此CFL數(shù)取K時(shí)，至少需2K+1個(gè)模板單元，此時(shí)推進(jìn)時(shí)間步長(zhǎng)可以達(dá)到Δt=KΔx/sa，即可能有2K個(gè)波束(左右各K個(gè))對(duì)目標(biāo)單元有貢獻(xiàn)，圖3僅給出了特征波速為正的波束示意圖[86]。根據(jù)疊波法思路，CFL=K的大時(shí)間步長(zhǎng)格式可寫(xiě)成

圖3 多波束在網(wǎng)格上的傳播[86]Fig.3 Multi-wave propagation on grid[86]

(9)

(10)

LeVeque[80，82]的數(shù)值實(shí)驗(yàn)表明，以上格式可以突破CFL≤1的限制，且在理論上CFL數(shù)可以取無(wú)限大。但實(shí)際上，由于線(xiàn)化假設(shè)等原因，CFL數(shù)只能取到有限值。

對(duì)于守恒律系統(tǒng)，與單方程不同的是間斷分解比較復(fù)雜，每一間斷可以分解出多個(gè)波，如圖1所示。可仍假設(shè)不同波之間相互自由穿越，互不干擾。這個(gè)假設(shè)對(duì)非線(xiàn)性方程也屬于線(xiàn)化近似，在一定CFL數(shù)范圍內(nèi)該近似足夠準(zhǔn)確，當(dāng)CFL數(shù)非常大時(shí)，它會(huì)導(dǎo)致解的扭曲，此時(shí)可以通過(guò)加密網(wǎng)格來(lái)彌補(bǔ)(可以證明[80]，當(dāng)網(wǎng)格加密時(shí)，該方法對(duì)所有CFL數(shù)都是收斂的)。圖4給出了Euler方程(每個(gè)間斷包含三波束)的多波束疊加示意圖[86]，這時(shí)波束雖然變得復(fù)雜，但在標(biāo)量情形下發(fā)展的掃描方法仍可直接推廣應(yīng)用。

圖4 Euler方程多波束在網(wǎng)格上的傳播[86]Fig.4 Multi-wave propagation on grid for Euler equations[86]

在方程組的情形下，式(10)可以推廣為

(11)

(12)

3 Godunov型大時(shí)間步長(zhǎng)格式的改進(jìn)

LeVeque[80-81]提出的大時(shí)間步長(zhǎng)格式主要針對(duì)標(biāo)量雙曲型守恒律方程，自其提出該思路以來(lái)，研究者一直試圖將其拓展到復(fù)雜問(wèn)題的求解中，主要的應(yīng)用包括氣體動(dòng)力學(xué)Euler方程和水動(dòng)力學(xué)淺水波方程兩個(gè)方面。主要的改進(jìn)也涉及兩個(gè)方面，一是對(duì)非線(xiàn)性方程的改進(jìn)處理方法，二是從標(biāo)量方程到方程組的推廣。但目前仍存在一些問(wèn)題有待解決，包括強(qiáng)非線(xiàn)性、膨脹波、同族波系干擾等帶來(lái)的問(wèn)題。

3.1 時(shí)間線(xiàn)插值技術(shù)

時(shí)間線(xiàn)(Time-line)插值技術(shù)是Guinot[83]提出的一種改進(jìn)方法，旨在減小數(shù)值通量計(jì)算的模板點(diǎn)數(shù)，從而進(jìn)一步降低計(jì)算量，其主要思路如下。

因u(xj+1/2,t)為處于xj+1/2界面上的Riemann問(wèn)題的解，則界面數(shù)值通量也可表達(dá)為

(13)

式中：

(14)

變量u可以寫(xiě)為如下的線(xiàn)性組合：

(15)

(16)

按照特征值的正負(fù)，可將式(16)改寫(xiě)為

(17)

考慮界面上的等價(jià)Riemann問(wèn)題，有

(18)

在特征變量空間內(nèi)，按照特征線(xiàn)法，有

(19)

考慮λ(m)為正值的情形，在CFL數(shù)不超過(guò)1的條件下，根據(jù)特征線(xiàn)法可將時(shí)間積分轉(zhuǎn)換為空間積分：

(20)

式中：xA點(diǎn)是特征線(xiàn)的起始點(diǎn)，在CFL≤1條件約束下，一般落于xj或xj+1單元。

在CFL數(shù)超過(guò)1的條件下，xA點(diǎn)將落在xj單元左側(cè)或xj+1單元右側(cè)。因?yàn)榻缑鎥j-1/2或xj+3/2處可能有激波存在，特征線(xiàn)將是不可逆的，基于逆向追蹤法的式(20)將無(wú)法直接使用。此時(shí)，界面通量可以改寫(xiě)成

(21)

式中：τ表示從界面xj-1/2發(fā)出的特征線(xiàn)與界面xj+1/2的交點(diǎn)，如圖5所示。仍考慮λ(m)為正值的情形，進(jìn)一步可有

圖5 時(shí)間線(xiàn)插值法示意圖[83]Fig.5 Schematic of time-line interpolation method

(22)

(23)

式中：τ′表示從界面xj+1/2逆向發(fā)出的特征線(xiàn)與界面xj-1/2的交點(diǎn)。式(23)的表達(dá)相當(dāng)于把特征線(xiàn)不可跨域的激波區(qū)采用相鄰的界面值來(lái)代替。

對(duì)于界面xj-1/2的變量采用線(xiàn)性重構(gòu)，則有

(24)

式中：中點(diǎn)值和斜率的表達(dá)式分別為

(25)

(26)

其中：s為線(xiàn)性重構(gòu)的斜率，可以對(duì)s引入限制器以防數(shù)值解振蕩，具體可參見(jiàn)文獻(xiàn)[83]。

該格式需要從左向右依次推進(jìn)，以得到所有波速為正的界面通量貢獻(xiàn)，從右向左依次推進(jìn)，以得到所有波速為負(fù)的界面通量貢獻(xiàn)。所以，數(shù)值通量表現(xiàn)出很強(qiáng)的邊界條件依賴(lài)性，需要妥善處理邊界條件，詳情可見(jiàn)文獻(xiàn)[83]。本質(zhì)上來(lái)講，因?yàn)闊o(wú)法實(shí)現(xiàn)真正意義上的局部推進(jìn)，該格式僅是一種半顯式格式，可以看作是對(duì)隱式格式的一種簡(jiǎn)化。

3.2 膨脹波的處理技術(shù)

2.2節(jié)提到的疊波法只能直接用于波系均為間斷(激波或接觸間斷)的情形，在非線(xiàn)性條件下會(huì)出現(xiàn)膨脹波，而膨脹波在tn+1時(shí)刻有可能跨越若干個(gè)區(qū)間。LeVeque[82]曾提出采用一束非物理膨脹間斷來(lái)代替膨脹波，并采用此法求解了一維等熵Euler方程。該方法計(jì)算效率高，但存在違背熵條件的非物理解的可能性。在疊波法中，如果對(duì)Riemann問(wèn)題分解后出現(xiàn)的膨脹波處理不當(dāng)，則在遇到較強(qiáng)膨脹波時(shí)會(huì)導(dǎo)致數(shù)值解惡化，甚至出現(xiàn)包含膨脹間斷的非物理解。

針對(duì)Euler 方程的特點(diǎn)，本文作者與Lee[84]提出采用若干個(gè)膨脹間斷波來(lái)近似間斷分解中出現(xiàn)的膨脹波，可稱(chēng)之為L(zhǎng)TS-Godunov-MW格式。如圖6所示，各個(gè)膨脹間斷波之間的狀態(tài)由膨脹波波頭和波尾的狀態(tài)通過(guò)線(xiàn)性插值得到。這里的膨脹間斷波嚴(yán)格意義上也是違反熵條件的，但是因?yàn)榕c壓力相比，熵變隨著近似波數(shù)目的增多衰減要快得多，所以使得最終數(shù)值解出現(xiàn)違反熵條件的概率急劇下降。嘗試了分別采取一束膨脹間斷波(其本質(zhì)就是一束膨脹激波，簡(jiǎn)稱(chēng)單波近似)、兩束膨脹間斷波(簡(jiǎn)稱(chēng)雙波近似)和三束膨脹間斷波(簡(jiǎn)稱(chēng)三波近似)的處理。理論上可以采用更多波束(多于三波)來(lái)代替膨脹波，且采用越多波越逼近于原膨脹波，但也相應(yīng)加大了計(jì)算量。數(shù)值實(shí)驗(yàn)表明，在絕大數(shù)情況下，兩波近似已可提供足夠的精度，既能避免非物理解的出現(xiàn)，又不致過(guò)大地增加計(jì)算量。

圖6 膨脹波的多波束近似示意圖[84]Fig.6 Schematic of multi-wave approximation of expansion wave[84]

圖7 膨脹波網(wǎng)格單元分解法示意圖[86]Fig.7 Schematic of grid cell decomposition method for expansion wave[86]

3.3 波系干擾的處理技術(shù)

對(duì)于非線(xiàn)性雙曲型守恒律組，如一維Euler方程，在CFL>0.5時(shí)，如圖8所示[84]，兩個(gè)相鄰間斷發(fā)出的波即可能發(fā)生干擾。對(duì)于諸如Euler方程的非線(xiàn)性方程，異族強(qiáng)波相互穿越并且穿越后波速與波強(qiáng)均會(huì)有改變，同族強(qiáng)波不易相交，如相交則合并成為一個(gè)強(qiáng)波。在實(shí)際問(wèn)題中，非線(xiàn)性波系的干擾非常復(fù)雜。相鄰的兩個(gè)間斷分解后，各發(fā)出3個(gè)波——左傳波(可以是左向激波或左向膨脹波)、接觸間斷、右傳波(可以是右向激波或右向膨脹波)，其中激波和膨脹波為非線(xiàn)性波，接觸間斷為線(xiàn)性波。

圖8 波系干擾示意圖[84]Fig.8 Schematic of interaction of wave systems[84]

從圖8可以看到，在給定的時(shí)間步Δt內(nèi)僅兩個(gè)相鄰間斷就可能出現(xiàn)多個(gè)波對(duì)干擾問(wèn)題，對(duì)于N個(gè)間斷的情形，將出現(xiàn)更多個(gè)波對(duì)的相互干擾。文獻(xiàn)[96]對(duì)兩族波之間的相互作用，包括激波、膨脹波和接觸間斷的兩兩相互作用進(jìn)行了詳細(xì)論述，給出了解析解。但是，在實(shí)際應(yīng)用中，考慮完整的波系干擾(包括不相鄰的兩個(gè)或多個(gè)間斷)，在編程和計(jì)算等實(shí)際操作上是難以實(shí)現(xiàn)的。為簡(jiǎn)化計(jì)算，文獻(xiàn)[84]提出僅考慮相鄰的兩個(gè)間斷發(fā)出的相鄰波——左間斷發(fā)出的右傳波(圖8中波a)和右間斷發(fā)出的左傳波(圖8中波b)——之間的碰撞問(wèn)題。經(jīng)以上簡(jiǎn)化后，只需處理兩個(gè)異族波的對(duì)撞問(wèn)題，如圖9所示。其對(duì)撞存在4種可能性：兩個(gè)激波的對(duì)撞、左膨脹波和右激波的對(duì)撞、左激波和右膨脹波的對(duì)撞、兩個(gè)膨脹波的對(duì)撞，可采用文獻(xiàn)[96]給出的解析解來(lái)實(shí)現(xiàn)。波a和波b碰撞后必然產(chǎn)生兩個(gè)新的非線(xiàn)性波和一個(gè)接觸間斷，即兩波之間的狀態(tài)由單一狀態(tài)M變?yōu)閮蓚€(gè)由接觸間斷所間隔的狀態(tài)M*(若波a和波b均為膨脹波時(shí)則不存在接觸間斷)。數(shù)值計(jì)算結(jié)果表明，通過(guò)采用相鄰間斷干擾精確處理的方法，使得格式的魯棒性大大提高，穩(wěn)定計(jì)算的CFL數(shù)范圍得到有效擴(kuò)展。

圖9 兩波碰撞示意圖[84]Fig.9 Schematic of interaction of two discontinuities[84]

3.4 近似Riemann求解器的應(yīng)用

在大時(shí)間步長(zhǎng)格式的研究中，文獻(xiàn)[80, 84, 86]的工作均采用了Riemann問(wèn)題的精確解來(lái)進(jìn)行間斷分解。近期Lindqvist[97]、Prebeg[98]等嘗試了采用近似Riemann求解器來(lái)進(jìn)行間斷分解。

1) LTS-Roe格式

Lindqvist等[97]針對(duì)標(biāo)量雙曲型方程，采用單參數(shù)(One Parameter)格式，描述了多種大時(shí)間步長(zhǎng)格式，其討論了數(shù)值黏性和通量差分分裂兩種表達(dá)形式的大時(shí)間步長(zhǎng)格式，并探討了其TVD性質(zhì)。基于數(shù)值黏性形式，標(biāo)量雙曲型方程的大時(shí)間步長(zhǎng)格式的數(shù)值通量可寫(xiě)為

(27)

基于通量差分分裂形式，標(biāo)量雙曲型方程的大時(shí)間步長(zhǎng)格式可寫(xiě)為

(28)

Lindqvist[97]等指出數(shù)值黏性和通量差分分裂兩種形式的大時(shí)間步長(zhǎng)格式是一一對(duì)應(yīng)的，具有如下關(guān)系：

(29a)

(29b)

式中:i∈{1,2,…,K-1}。

(30)

2K+1點(diǎn)的LTS-Lax-Friedrichs格式的數(shù)值黏性系數(shù)為

(31a)

(31b)

2K+1點(diǎn)的LTS-Roe格式的數(shù)值黏性系數(shù)為

(32a)

(32b)

Lindqvist[97]等指出LTS-Roe格式和LTS-Lax-Friedrichs格式分別為耗散最小和最大的具有TVD性質(zhì)的大時(shí)間步長(zhǎng)格式。LTS-Lax-Friedrichs格式數(shù)值耗散過(guò)大，在工程計(jì)算中不實(shí)用；LTS-Roe格式數(shù)值耗散在某些情形下不足，導(dǎo)致數(shù)值解在間斷波附近出現(xiàn)非物理振蕩，并有可能出現(xiàn)違背熵條件的可能，得到非物理的解。Lindqvist等結(jié)合了兩種格式的優(yōu)劣，提出了混合型的LTS-RoeLxF(β)格式：

(33)

該格式屬于LTS-Roe格式和LTS-Lax-Friedrichs格式的線(xiàn)性組合，通過(guò)LTS-Lax-Friedrichs格式引入了更多的數(shù)值黏性，消弱了數(shù)值解在激波附近的非物理振蕩，并使得數(shù)值解朝熵解靠攏。其在數(shù)值算例中探索了參數(shù)β對(duì)計(jì)算結(jié)果的影響。

Lindqvist等[97]采用局部線(xiàn)性近似方法，將基于標(biāo)量雙曲型方程建立的大時(shí)間步長(zhǎng)格式推廣至方程組情形，但是由于方程組的復(fù)雜性，在標(biāo)量情形所討論的有關(guān)格式的各種性質(zhì)將不再全部有效。盡管如此，Lindqvist等仍采用一維Euler方程激波管問(wèn)題作為數(shù)值算例，初步驗(yàn)證了所建立的方法。

2) LTS-HLL格式

Prebeg等[98]在Lindqvist等[97]研究工作的基礎(chǔ)上提出了基于HLL[57]格式和HLLC格式[58]的大時(shí)間步長(zhǎng)格式(LTS-HLL/HLLC)。其采用更廣義的雙參數(shù)(Two Parameters)描述形式，表達(dá)了包含Lindqvist等[97]提出的多種格式在內(nèi)的大時(shí)間步長(zhǎng)格式，重點(diǎn)探討了基于HLL和HLLC格式的大時(shí)間步長(zhǎng)格式LTS-HLL和LTS-HLLC格式。

HLL格式是Harten等[57]提出的一個(gè)新穎的方法，其通過(guò)近似求解Riemann問(wèn)題來(lái)給出數(shù)值通量，該格式假設(shè)Riemann問(wèn)題的解由兩道波組成，并且通過(guò)估算給出兩道波各自的傳播速度，可快速得到數(shù)值通量，不像Roe格式和Osher格式那樣詳細(xì)給出Riemann問(wèn)題解的波系結(jié)構(gòu)。該格式比Roe格式計(jì)算效率更高，但是由于對(duì)Riemann問(wèn)題的間斷分解后的物理解結(jié)構(gòu)采用了兩波近似，使得格式對(duì)接觸間斷的分辨率下降。針對(duì)這一問(wèn)題，Toro等[58]通過(guò)添加接觸間斷的影響，擴(kuò)展兩波近似為三波形式，給出了改近的HLL格式，稱(chēng)為HLLC格式，提高了對(duì)接觸間斷的分辨率。

Prebeg等[98]通過(guò)對(duì)標(biāo)量雙曲型守恒律引入兩參數(shù)形式的數(shù)值解，得到了更廣義的方法。LTS-HLL格式寫(xiě)成數(shù)值黏性形式為

(34a)

(34b)

式中：SL和SR分別為HLL格式中近似Riemann問(wèn)題的解中兩道波的波速。

LTS-HLL格式也可以改寫(xiě)成疊波形式

(35)

式中：w為波強(qiáng)，即分解后所得的波兩側(cè)狀態(tài)的差量。其中S1=SL和S2=SR，并且有

(36)

LTS-HLLC格式直接針對(duì)方程組構(gòu)造，以一維Euler方程為例，這里波的個(gè)數(shù)變?yōu)?個(gè)，寫(xiě)成疊波形式為

(37)

式中：

S1=SL,S2=SC,S3=SR

(38)

其中：SC為HLLC格式中接觸間斷的傳播速度，進(jìn)一步有

(39)

容易理解，LTS-HLL格式耗散較大，LTS-HLLC格式耗散較小，分辨率有明顯提高。但是從數(shù)值實(shí)驗(yàn)也發(fā)現(xiàn)LTS-HLLC格式在間斷附近容易滋生非物理的數(shù)值振蕩。為改善該類(lèi)格式的數(shù)值特性，參考HLLE格式[99]的波速選取，Prebeg建立了LTS-HLLE格式，為進(jìn)一步改進(jìn)格式，Prebeg等[98]提出了LTS-HLLE格式與LTS-Lax-Friedrichs格式的混合格式LTS-HLLELxF(β)格式：

(40)

4 Godunov型大時(shí)間步長(zhǎng)格式的高階精度推廣方法

自LeVeque[80-81]提出大時(shí)間步長(zhǎng)格式以來(lái)，所發(fā)展的絕大多數(shù)格式均只具有一階代數(shù)精度。雖然在同等代數(shù)精度下，隨著CFL數(shù)的增大，Godunov型大時(shí)間步長(zhǎng)格式分辨率逐漸提高，甚至超過(guò)了常規(guī)二價(jià)代數(shù)精度的格式，但是構(gòu)造更高階的格式對(duì)于提高格式分辨率還是大有裨益的。本文作者[85]首次嘗試了構(gòu)造二階精度大時(shí)間步長(zhǎng)的工作，其綜合采用了Billett和Toro[100]提出的加權(quán)平均狀態(tài)(Weighted Average State，WAS)方法和LeVeque的疊波法[80]，建立了具有二階精度的Godunov型大時(shí)間步長(zhǎng)格式，但目前CFL數(shù)僅能達(dá)到2.0。

4.1 加權(quán)平均狀態(tài)方法

WAS方法是Toro等[100]提出的用于求解雙曲型守恒律的一種二階精度的數(shù)值格式。其仍然按照數(shù)據(jù)重構(gòu)、間斷分解和通量計(jì)算等3步求解數(shù)值通量，在數(shù)據(jù)重構(gòu)步仍然采用分段常數(shù)重構(gòu)，間斷分解仍采用精確或近似Riemann解進(jìn)行，Toro方法的核心思想在通量計(jì)算步，其將原始的Godunov格式的通量計(jì)算方法進(jìn)行了推廣，即

(41)

(42)

式中：u*(x,t)為tn時(shí)刻間斷分解后的數(shù)值解。式(42)是一種廣義的數(shù)值通量計(jì)算方法，本節(jié)的討論僅局限于

(43)

的局域。與原始Godunov型格式僅有時(shí)間方向的積分不同，改進(jìn)的數(shù)值通量的計(jì)算引入了空間方向的積分。為獲得二階精度格式，這里對(duì)時(shí)間方向采用中點(diǎn)積分方式，則有

(44)

如圖10所示[11]，對(duì)區(qū)域[-Δx/2,Δx/2]進(jìn)行積分，可得

圖10 WAS方法通量計(jì)算示意圖[11]Fig.10 Schematic of flux calculation of WAS method[11]

(45)

(46)

式中：Ck為波速為Sk的第k個(gè)波對(duì)應(yīng)的CFL數(shù)。式(45)也可以改寫(xiě)為

(47)

Toro推導(dǎo)出了具有TVD性質(zhì)的WAS格式來(lái)保證數(shù)值計(jì)算的穩(wěn)定性。引入TVD約束條件后，WAS方法的數(shù)值通量式(47)可改寫(xiě)為

(48)

式中：Φ為限制器函數(shù)，具體形式見(jiàn)文獻(xiàn)[100]。

4.2 WAS型二階精度大時(shí)間步長(zhǎng)Godunov格式

采用WAS方法的起點(diǎn)是引入狀態(tài)變量的加權(quán)平均式(42)，然后采用式(41)給出數(shù)值通量。對(duì)于Euler方程，式(42)中的變量狀態(tài)可以選擇守恒變量、原始變量或其他狀態(tài)參數(shù)，均不影響格式的守恒性。為了同原有的一階Godunov型大時(shí)間步長(zhǎng)格式相容，文獻(xiàn)[85]仍采用守恒變量，具體構(gòu)造思路如下。

圖11 WAS方法中單束波的傳播示意圖[85]Fig.11 Illustration of single wave propagation for WAS method[85]

步驟1設(shè)置

U(j,1)=U(j,2)=Uj

(49)

步驟3由間斷xj+1/2處兩側(cè)的兩個(gè)半?yún)^(qū)間值(U(j,1)和U(j,2))的平均值給出計(jì)算數(shù)值通量Fj+1/2的Uj+1/2，即

(50)

步驟4采用守恒格式更新變量值，即

(51)

(52)

(53)

Toro[11]指出，所選擇的變量q只需滿(mǎn)足穿過(guò)每一個(gè)波的跳躍不為0即可。通常可選擇密度ρ或內(nèi)能E。本文作者[85]選擇了密度ρ。J的取值與迎風(fēng)方向有關(guān)，在此取

(54)

由于在時(shí)間方向采用了中點(diǎn)積分，故而該格式的CFL數(shù)理論上僅能達(dá)到2.0，可以考慮構(gòu)造使用更多點(diǎn)的積分進(jìn)一步擴(kuò)展CFL數(shù)，但算法的復(fù)雜度將大幅增加。

5 多維問(wèn)題推廣

目前提出的大時(shí)間步長(zhǎng)格式都是針對(duì)一維問(wèn)題的，向多維情形的推廣主要采用維數(shù)分裂法[101]。維數(shù)分裂法雖然會(huì)引入一定誤差，但是在一定精度范圍內(nèi)還是可以滿(mǎn)足要求。

5.1 維數(shù)分裂

考慮多維Euler方程[85]，其離散形式可寫(xiě)成

(55)

式中：λ=Δt/Δξ=Δt/Δη=Δt/Δζ，這里ξ、η、ζ為曲線(xiàn)貼體坐標(biāo)系下的廣義坐標(biāo)。記算子δi(·)=δi+1/2(·)-δi-1/2(·)，則式(55)可以寫(xiě)成算子形式

(56)

若給式(56)右端添加一個(gè)二階小量:

(57)

則可得到

(58)

以上形式的全離散型顯式差分形式可以在i、j、k3個(gè)方向上分別進(jìn)行準(zhǔn)一維的計(jì)算，程序設(shè)計(jì)也大大簡(jiǎn)化。

5.2 對(duì)稱(chēng)維數(shù)分裂

式(58)可以改寫(xiě)成

(59)

式中：T為推進(jìn)算子。將T分解為3個(gè)方向的子推進(jìn)算子，則有

(60)

常用的方法是如式(60)所示的簡(jiǎn)單的維數(shù)分裂，但對(duì)于某些情形該分裂方法可能會(huì)導(dǎo)致一定的數(shù)值振蕩，尤其在大CFL數(shù)時(shí)問(wèn)題可能更突出。Dong和Liu[86]發(fā)現(xiàn)這是由于非對(duì)稱(chēng)的維數(shù)分裂造成的，他們建議使用所謂的對(duì)稱(chēng)算子分裂。以二維情形為例：

(61)

圖12[86]給出了算子分裂的示意，使用3個(gè)間斷的解為例來(lái)分析，3個(gè)間斷形成一個(gè)三波點(diǎn)，可用于考察簡(jiǎn)單分裂和對(duì)稱(chēng)分裂對(duì)于斜間斷數(shù)值解的影響。從圖中可以看出，對(duì)于斜間斷解，簡(jiǎn)單分裂與分裂的順序有關(guān)，一個(gè)在左上，一個(gè)在右下，兩者間存在一個(gè)如下的非對(duì)稱(chēng)小量差異：

圖12 簡(jiǎn)單分裂和對(duì)稱(chēng)分裂的對(duì)比[86]Fig.12 Comparison of simple split and symmetric split[86]

(62)

式中：Tξ=I-τTξ；Tη=I-τTη。即便是這個(gè)小量的存在，也會(huì)改變斜間斷的形狀，使之發(fā)生畸變。在多維情形，這個(gè)畸變可能會(huì)導(dǎo)致間斷附近出現(xiàn)非物理的振蕩。使用對(duì)稱(chēng)分裂方法可以去除該非對(duì)稱(chēng)差異，可以抑制斜間斷的變形及其附近的數(shù)值解振蕩。

5.3 分片Riemann問(wèn)題的引入

一維情形下，Euler方程含有3個(gè)變量，但對(duì)于三維情形，每一方向均需要處理5個(gè)變量的問(wèn)題。文獻(xiàn)[85]參考了文獻(xiàn)[102]的方法，采用當(dāng)?shù)匦D(zhuǎn)矩陣 將其變換至局部笛卡爾坐標(biāo)系。該變換矩陣的形式為

(63)

(64)

式中：

(65)

(66)

經(jīng)過(guò)坐標(biāo)變換，物理域中位于(x,y,z)點(diǎn)處ξ、η或ζ方向的間斷分解問(wèn)題可轉(zhuǎn)換為局部笛卡爾坐標(biāo)系下的一維Riemann問(wèn)題。該問(wèn)題和原始笛卡爾坐標(biāo)系下的一維問(wèn)題具有相同的形式，但多出了兩個(gè)速度分量，Toro[11]稱(chēng)之為分片Riemann問(wèn)題。以沿ξ方向的三維分片Riemann問(wèn)題為例，其變換后的形式可寫(xiě)為

(67)

其初值記為

(68)

η和ζ方向的分析類(lèi)同。完成坐標(biāo)變換后即可直接進(jìn)行計(jì)算，給出每一步的數(shù)值解后再采用當(dāng)?shù)匦D(zhuǎn)矩陣T的逆矩陣T-1:

(69)

(70)

5.4 邊界條件處理

在絕大多數(shù)算法研究的論文中僅考慮了一維激波管算例，在該算例中，邊界條件均采取一階外推給定，處理較為容易。在數(shù)值模擬多維問(wèn)題時(shí)，邊界條件的設(shè)置則相對(duì)復(fù)雜。本文作者[84-85]通過(guò)在計(jì)算域邊界外添加虛網(wǎng)格的方式簡(jiǎn)化邊界條件的設(shè)置，時(shí)間步長(zhǎng)越大則所需的虛網(wǎng)格數(shù)量越多，如圖13所示。多維計(jì)算中，入流遠(yuǎn)場(chǎng)邊界虛網(wǎng)格上的值采用入流無(wú)反射條件給定；出流邊界采用出流無(wú)反射外推邊界條件；壁面邊界虛網(wǎng)格上的值由鏡像反射條件[103]給定。對(duì)于Euler方程，其壁面切向流動(dòng)不受約束，可使用鏡像反射法給出邊界外的虛流場(chǎng)。對(duì)于N-S方程，壁面滿(mǎn)足無(wú)滑移條件，虛網(wǎng)格上速度可由反對(duì)稱(chēng)鏡像法給定。

圖13 邊界附近虛網(wǎng)格示意圖[84]Fig.13 Schematic of ghost cells near boundary[84]

6 Godunov型大時(shí)間步長(zhǎng)格式的性能分析

自LeVeque提出大時(shí)間步長(zhǎng)格式以來(lái)，有多個(gè)工作研究了該類(lèi)格式的性能，但絕大多數(shù)工作集中在針對(duì)標(biāo)量雙曲型守恒律的穩(wěn)定性條件、TVD特性、熵穩(wěn)定性、分辨率等性能的探討，對(duì)于方程組的情形，目前開(kāi)展的理論分析仍很少見(jiàn)。當(dāng)然，計(jì)算效率作為大時(shí)間步長(zhǎng)格式的主要研究出發(fā)點(diǎn)，也得到了較為廣泛的關(guān)注。本節(jié)簡(jiǎn)要介紹有關(guān)的代表性工作。

6.1 收斂特性

1) 穩(wěn)定性條件

直觀上看，該類(lèi)格式似乎違背了CFL條件，無(wú)法實(shí)現(xiàn)穩(wěn)定計(jì)算，其實(shí)即使在CFL>1時(shí)，其依然滿(mǎn)足穩(wěn)定性條件。CFL條件指出數(shù)值解的依賴(lài)域不得小于物理?xiàng)l件下的依賴(lài)域，其實(shí)該條件并未限定具體的CFL數(shù)。根據(jù)CFL條件，當(dāng)格式的模版點(diǎn)為2K+1時(shí)，其CFL數(shù)的上限可達(dá)到K。文獻(xiàn)[85]指出Godunov型大時(shí)間步長(zhǎng)格式正是通過(guò)增加CFL數(shù)來(lái)自動(dòng)增加模版點(diǎn)的數(shù)目以保證滿(mǎn)足CFL條件。該格式采用波的傳播速度來(lái)完成模版點(diǎn)選擇的自適應(yīng)性，即CFL越大，每一時(shí)間步長(zhǎng)中涉及的模版點(diǎn)數(shù)就越多。

為何普通多點(diǎn)格式如二階TVD格式、MUSCL格式、WENO格式等在時(shí)間方向顯式離散時(shí)，其CFL數(shù)仍不能超過(guò)1。這是因?yàn)樵诙ATVD、MUSCL、WENO等格式中，除了3 個(gè)主要模版點(diǎn)之外，其他的模版點(diǎn)僅存在于系數(shù)之中(主要體現(xiàn)在限制器中)，處于次要地位，因而不能獲得更大的CFL數(shù)。但該處理可帶來(lái)光滑區(qū)格式精度的提高。

2) TVD性質(zhì)

TVD是標(biāo)量雙曲型守恒律方程的一個(gè)重要特性，根據(jù)Lax-Wendroff定理，只有具有TVD性質(zhì)的守恒型相容格式才能收斂于弱解，并且可以保證數(shù)值解的無(wú)非物理振蕩。對(duì)于3點(diǎn)格式的TVD特性，Harten等開(kāi)展了充分的研究，取得了很大的成功。但是對(duì)于多點(diǎn)格式的TVD特性，長(zhǎng)期以來(lái)開(kāi)展的研究并不多，根據(jù)第5節(jié)的討論，大時(shí)間步長(zhǎng)格式就是典型的多點(diǎn)格式。

LeVeque[104]證明了對(duì)于標(biāo)量雙曲型守恒律方程，LTS-Godunov格式是具有TVD性質(zhì)的，因而可以收斂到弱解。Jameson和Lax[105-106]最早研究了廣義2K+1點(diǎn)格式的TVD特性，證明了該類(lèi)多點(diǎn)格式具有TVD格式的充分條件。Lindqvist等[97]擴(kuò)展了Jameson和Lax的結(jié)論，給出了標(biāo)量雙曲型方程的2K+1點(diǎn)格式具有TVD性質(zhì)的充分必要條件，給出了2K+1點(diǎn)TVD格式的數(shù)值黏性系數(shù)需滿(mǎn)足的上下邊界條件，指出LTS-Roe格式和LTS-Lax-Friedrichs格式分別為耗散最小和最大的具有TVD性質(zhì)的大時(shí)間步長(zhǎng)格式，據(jù)此發(fā)展了3.4節(jié)的LTS-RoeLxF(β)大時(shí)間步長(zhǎng)格式。

應(yīng)該說(shuō)明的是，上述工作都是針對(duì)標(biāo)量方程開(kāi)展的研究，在推廣至方程組的過(guò)程中采用了特征場(chǎng)分解的方法，對(duì)于分解后的單個(gè)特征場(chǎng)則分別采用標(biāo)量情形發(fā)展的計(jì)算格式。

3) 熵穩(wěn)定性

在標(biāo)量雙曲型守恒律方程情形下，具有TVD性質(zhì)的守恒型格式收斂于弱解，但是弱解不一定是熵解，只有具有熵穩(wěn)定性的格式才能收斂至熵解。對(duì)于常規(guī)格式的熵穩(wěn)定性，目前已有較多的結(jié)論，但是對(duì)于大時(shí)間步長(zhǎng)格式的熵穩(wěn)定性的研究尚很不充分。LeVeque[104]曾推測(cè)，只要間斷分解過(guò)程中使用的Riemann求解器可提供熵解，則采用疊波法的LTS-Godunov格式就可以收斂到熵解，但時(shí)至今日，該推測(cè)依然沒(méi)有得到嚴(yán)格證明。Wang和Warnecke[107-108]曾證明：在CFL≤2的條件下，如果通量函數(shù)是單調(diào)的，則LTS-Godunov格式是熵穩(wěn)定的；如果初始值是單調(diào)的，則LTS-Godunov格式對(duì)于任意CFL數(shù)都是熵穩(wěn)定的。Wang等[109]證明了在某些特定的初值條件下，LTS-Godunov格式對(duì)于任意CFL數(shù)都是熵穩(wěn)定的。Tang和Warnecke[110]研究了LTS-Lax-Friedrichs格式，證明了因其為單調(diào)格式，故而是熵穩(wěn)定的。

Lindqvist等[97]針對(duì)標(biāo)量雙曲型方程指出，LTS-Roe格式的數(shù)值耗散小于LTS-Godunov格式的耗散，因而可能出現(xiàn)不滿(mǎn)足熵條件的情形，并且指出Riemann解中膨脹波的處理是違背熵條件的主要因素。因?yàn)樵诏B波法中，對(duì)膨脹波作了多波束近似，故必然存在出現(xiàn)膨脹間斷的風(fēng)險(xiǎn)。針對(duì)膨脹波處理，Lindqvist等提出了兩種方法：對(duì)于超聲速膨脹波，其建議采用隨機(jī)選取法來(lái)改變時(shí)間步長(zhǎng)，進(jìn)而竭力避免膨脹間斷的出現(xiàn)；對(duì)于跨聲速膨脹波，其沿用了Harten[34]提出的熵修正函數(shù)。采用這兩種策略改進(jìn)后的格式稱(chēng)為L(zhǎng)TS-Roe*，數(shù)值實(shí)驗(yàn)表明其解的熵穩(wěn)定性有了明顯改善。

6.2 分辨率

多個(gè)研究均表明Godunov型大時(shí)間步長(zhǎng)格式的分辨率與時(shí)間步長(zhǎng)相關(guān)，且時(shí)間步長(zhǎng)越大，則格式的分辨率越高，即引入的耗散越小。本節(jié)從誤差分析、數(shù)值耗散機(jī)制分析和數(shù)值實(shí)驗(yàn)驗(yàn)證3個(gè) 角度簡(jiǎn)述對(duì)該特性的研究。

1) 誤差分析

文獻(xiàn)[86]給出標(biāo)量雙曲型守恒律方程的大時(shí)間步長(zhǎng)格式的截?cái)嗾`差系數(shù)為

(71)

圖14 誤差系數(shù)隨CFL數(shù)變化[86]Fig.14 Error coefficients vs CFL numbers[86]

2) 數(shù)值耗散機(jī)制分析

Toro[11]提出典型的Godunov型格式的求解過(guò)程可以分為3個(gè)步驟，即重構(gòu)、演化和平均。重構(gòu)主要實(shí)現(xiàn)網(wǎng)格節(jié)點(diǎn)或單元內(nèi)平均值分布的高階近似，故一般不會(huì)引入耗散。演化時(shí)，耗散大小與所采用的Riemann求解器相關(guān)，特別是采用精確Riemann解時(shí)，可以認(rèn)為不引入耗散。平均是為了實(shí)現(xiàn)了間斷分解后的復(fù)雜波系解在下一時(shí)刻向網(wǎng)格上的投影，該步是引入耗散的主要環(huán)節(jié)。文獻(xiàn)[85]指出，大時(shí)間步長(zhǎng)格式的分辨率隨CFL數(shù)的增大而逐漸提高的原因在于， Godunov型格式求平均步會(huì)產(chǎn)生較大的耗散，導(dǎo)致間斷面的抹平，對(duì)激波和接觸間斷分辨率不高，而大時(shí)間步長(zhǎng)格式與普通Godunov型格式相比，由于采取了大的時(shí)間步長(zhǎng)，在相同的物理計(jì)算時(shí)間內(nèi)，減少了求平均的次數(shù)，從而降低了耗散，提高了對(duì)間斷的分辨率。CFL數(shù)越大，相同時(shí)間內(nèi)求平均的次數(shù)就越少，所以格式的分辨率也就越高。

3) 數(shù)值實(shí)驗(yàn)驗(yàn)證

自大時(shí)間步長(zhǎng)格式提出以來(lái)，多個(gè)研究者均對(duì)其進(jìn)行了數(shù)值算例驗(yàn)證。文獻(xiàn)[80，82]針對(duì)標(biāo)量方程和等熵Euler方程，文獻(xiàn)[83，111-115]針對(duì)淺水波方程，文獻(xiàn)[84，86，97]針對(duì)Euler方程均進(jìn)行了大量數(shù)值實(shí)驗(yàn)，計(jì)算結(jié)果均表明Godunov型大時(shí)間步長(zhǎng)格式具有隨著CFL數(shù)增大而分辨率提高的性質(zhì)。這與上述的誤差分析和數(shù)值耗散機(jī)制分析的結(jié)論是一致的。

6.3 計(jì)算效率

在針對(duì)大時(shí)間步長(zhǎng)格式的研究中，絕大多數(shù)工作僅針對(duì)格式本身的特性，并未直接給出計(jì)算CPU時(shí)間的直接對(duì)比。文獻(xiàn)[84,86,115]等各自給出了LTS-Godunov格式的CPU時(shí)間，評(píng)估了所發(fā)展的各種大時(shí)間步長(zhǎng)格式的計(jì)算效率，給出了其在典型模型問(wèn)題求解中的加速特性。

文獻(xiàn)[84]給出了采用雙波近似的LTS-Godunov-MW格式對(duì)一維Sod激波管和二維NACA0012翼型的流場(chǎng)計(jì)算的效率對(duì)比,如圖15所示[84]。在CFL數(shù)5.0以?xún)?nèi)，對(duì)一維Sod激波管問(wèn)題格式加速比隨CFL數(shù)增加而提高，在CFL數(shù)5.0時(shí)能達(dá)到23%；對(duì)二維NACA0012翼型問(wèn)題加速特性有所下降，在CFL數(shù)5.0時(shí)加速比能達(dá)到42%左右。但該格式效率還存在優(yōu)化的空間。

圖15 文獻(xiàn)[84]計(jì)算效率隨CFL數(shù)變化Fig.15 Computational efficiency vs CFL numbers in Ref.[84]

文獻(xiàn)[86]給出了改進(jìn)后的LTS-WA格式對(duì)一維Burges方程、Sod激波管問(wèn)題和二維Ringleb流動(dòng)的流場(chǎng)計(jì)算的效率對(duì)比,如圖16所示[86]。對(duì)于一維Burges方程，LTS-WA格式在CFL數(shù)9條件下的加速比約為53%；對(duì)于一維Sod激波管問(wèn)題，LTS-WA格式在CFL數(shù)9條件下的加速比能達(dá)到21%；對(duì)于二維Ringleb流動(dòng)，LTS-WA格式在CFL數(shù)8條件下的加速比能達(dá)到25%。

圖16 文獻(xiàn)[86]計(jì)算效率隨CFL數(shù)變化Fig.16 Computational efficiency vs CFL numbers in Ref.[86]

文獻(xiàn)[115]給出了LTS-Roe格式在一維淺水波方程求解中的效率對(duì)比，并研究了膨脹波多波束近似和CFL數(shù)之間的關(guān)系，多波束處理會(huì)增加計(jì)算時(shí)間，但總體上來(lái)看增大CFL數(shù)還是可明顯提高計(jì)算效率，如圖17所示[115]。在雙波近似條件下，以原始Roe格式在CFL數(shù)0.8條件下的計(jì)算時(shí)間作為參考，LTS-Roe格式在CFL數(shù)4條件下的加速比能達(dá)到45%。

圖17 文獻(xiàn)[115]計(jì)算效率隨CFL數(shù)變化Fig.17 Computational efficiency vs CFL numbers in Ref.[115]

7 典型問(wèn)題應(yīng)用示例

Godunov型顯式大時(shí)間步長(zhǎng)格式的主要應(yīng)用對(duì)象為以雙曲型守恒律方程組為控制方程的問(wèn)題，主要包含空氣動(dòng)力學(xué)領(lǐng)域的Euler方程[82，84，86，97-98]、水力學(xué)領(lǐng)域的淺水波方程[111-115]和電磁學(xué)領(lǐng)域的Maxwell方程[116]等典型問(wèn)題。本文主要關(guān)注其在空氣動(dòng)力學(xué)領(lǐng)域的應(yīng)用，給出了其在激波管、翼型、機(jī)翼等典型一維、二維和三維流動(dòng)的應(yīng)用示例。

7.1 一維Sod激波管問(wèn)題

Sod激波管問(wèn)題[117]是CFD計(jì)算格式驗(yàn)證的經(jīng)典算例，幾乎所有的格式研究都采用其作為驗(yàn)證算例。文獻(xiàn)[84，86，97-98]均采用所發(fā)展的大時(shí)間步長(zhǎng)格式對(duì)該問(wèn)題進(jìn)行了計(jì)算。文獻(xiàn)[84]的計(jì)算結(jié)果如圖18所示，其通過(guò)采用多波束對(duì)膨脹波進(jìn)行近似處理，可以很好地避免非物理解的產(chǎn)生，并且隨著CFL數(shù)的增大，格式對(duì)間斷波的分辨率顯著提高。文獻(xiàn)[86]的計(jì)算結(jié)果如圖19所示，同樣也表明隨著CFL數(shù)的增大，格式的分辨率提高。文獻(xiàn)[97]等的計(jì)算結(jié)果如圖20所示，其通過(guò)熵修正技術(shù)可以使得LTS-Roe*格式在CFL數(shù)6.0范圍內(nèi)獲得較為滿(mǎn)意的結(jié)果。文獻(xiàn)[98]的結(jié)果如圖21所示，其發(fā)展的LTS-HLL格式基本呈現(xiàn)無(wú)數(shù)值振蕩，但耗散相對(duì)較大，LTS-HLLC格式耗散明顯減小，呈現(xiàn)略微數(shù)值振蕩。

圖18 文獻(xiàn)[84]Sod激波管問(wèn)題計(jì)算結(jié)果Fig.18 Numerical results of Sod shock tube problem in Ref.[84]

圖19 文獻(xiàn)[86]Sod激波管問(wèn)題計(jì)算結(jié)果Fig.19 Numerical results of Sod shock tube problem in Ref.[86]

圖20 文獻(xiàn)[97]Sod激波管問(wèn)題計(jì)算結(jié)果Fig.20 Numerical results of Sod shock tube problem in Ref.[97]

圖21 文獻(xiàn)[98]Sod激波管問(wèn)題計(jì)算結(jié)果Fig.21 Numerical results of Sod shock tube problem in Ref.[98]

7.2 二維翼型繞流

目前大部分的大時(shí)間步長(zhǎng)格式研究工作主要集中在一維格式的構(gòu)造上，對(duì)二維問(wèn)題進(jìn)行計(jì)算的工作不多。Dong和Liu[86]曾采用LTS-WA格式計(jì)算了二維Riemann問(wèn)題，但仍采用了等距網(wǎng)格。文獻(xiàn)[84]首次將LTS-Godunov格式用于二維NACA0012翼型繞流計(jì)算。來(lái)流條件分別為馬赫數(shù)Ma=0.8、攻角α=1.25°和馬赫數(shù)Ma=1.2、攻角α=7.0°。文中兩個(gè)繞流情形均給出了CFL數(shù)從1.0～5.0的計(jì)算結(jié)果，可以發(fā)現(xiàn)隨著CFL數(shù)的增大，LTS-Godunov格式的分辨率逐漸提高，該趨勢(shì)尤其明顯地體現(xiàn)在對(duì)下翼面強(qiáng)度較弱的激波的分辨能力上，如圖22所示。

圖22 NACA0012翼型繞流計(jì)算結(jié)果(Ma=0.8、α=1.25°)[84]Fig.22 Numerical results of flow field around NACA0012 airfoil (Ma=0.8，α=1.25°)[84]

7.3 三維機(jī)翼繞流

目前僅文獻(xiàn)[84]和Dong[86]等給出了采用LTS-Godunov格式計(jì)算M6機(jī)翼三維繞流的結(jié)果，計(jì)算的來(lái)流馬赫數(shù)Ma=0.839 5、攻角α=3.06°。此算例為超臨界狀態(tài)，流場(chǎng)結(jié)構(gòu)較復(fù)雜，上翼面共有3道激波：前一道為弱激波，激波后仍為超聲速流動(dòng)；后一道激波為強(qiáng)激波，波后為亞聲速流動(dòng)；此兩道激波在翼梢附近相交，匯合成λ型第3道最強(qiáng)激波。下翼面均為亞聲速流動(dòng)。圖23 給出了文獻(xiàn)[84]的部分計(jì)算結(jié)果，顯示了CFL數(shù)1.0～5.0時(shí)的上壁面壓力云圖。從圖中可以看出，大時(shí)間步長(zhǎng)Godunov型格式能正確捕捉到上翼面由于兩道激波相交而形成的λ激波，且和一維及二維情形相似，隨著CFL數(shù)的提高，該格式的分辨率也逐步提高，主要表現(xiàn)在λ激波的捕捉上。文獻(xiàn)[84]還給出了相應(yīng)CFL數(shù)下的壓力分布及其與風(fēng)洞試驗(yàn)結(jié)果的對(duì)比，驗(yàn)證了該類(lèi)格式的計(jì)算可靠性。文獻(xiàn)[86]給出了CFL數(shù)最大到8.0 的計(jì)算結(jié)果，得到的流場(chǎng)分布如圖24所示，結(jié)論與文獻(xiàn)[84]基本一致，且文獻(xiàn)[86]針對(duì)80% 展長(zhǎng)處的壓力分布捕捉精度不足的問(wèn)題，通過(guò)對(duì)網(wǎng)格的改進(jìn)得到了更優(yōu)化的計(jì)算結(jié)果。

圖23 文獻(xiàn)[84]M6機(jī)翼繞流計(jì)算結(jié)果(Ma=0.839 5、α=3.06°)Fig.23 Numerical results of flow field around M6 wing (Ma=0.839 5, α=3.06°) in Ref.[84]

圖24 文獻(xiàn)[86]M6機(jī)翼繞流計(jì)算結(jié)果(Ma=0.839 5、α=3.06°)Fig.24 Numerical results of flow field around M6 wing (Ma=0.839 5,α=3.06°) in Ref.[86]

7.4 高超聲速雙橢球繞流

本文作者[85]還給出了LTS-Godunov格式對(duì)高超聲速雙橢球繞流的計(jì)算結(jié)果，來(lái)流馬赫數(shù)Ma為4.0、攻角α為0°，計(jì)算結(jié)果如圖25所示。其指出大時(shí)間步長(zhǎng)Godunov格式在不同CFL數(shù)下的計(jì)算結(jié)果和Haten-Yee TVD格式[118]所得到的結(jié)果相符，都能準(zhǔn)確地捕捉到橢球頭部的脫體弓形激波和兩橢球鑲嵌部位附近的二次激波，圖中給出了CFL數(shù)為1.0，3.0和5.0時(shí)雙橢球物面及流場(chǎng)對(duì)稱(chēng)面上的等壓線(xiàn)分布和對(duì)稱(chēng)面上的上下壁面壓力分布。

圖25 雙橢球高超聲速繞流計(jì)算結(jié)果(Ma=4.0、α=0°)[85]Fig.25 Numerical results of hypersonic flow field around double ellipsoid (Ma=4.0,α=0°)[85]

8 大時(shí)間步長(zhǎng)格式研究的發(fā)展方向

自LeVeque提出大時(shí)間步長(zhǎng)格式的概念以來(lái)，經(jīng)過(guò)30多年的發(fā)展，Godunov型大時(shí)間步長(zhǎng)格式引起了多個(gè)研究者的興趣。特別是在近年來(lái)隨著超級(jí)計(jì)算機(jī)系統(tǒng)的迅速發(fā)展，該類(lèi)格式因其優(yōu)異的并行特性，得到了更廣泛的關(guān)注。經(jīng)過(guò)逐漸探索和完善，愈來(lái)愈接近推向工程應(yīng)用計(jì)算的階段，但是仍有一些問(wèn)題有待解決。根據(jù)作者觀點(diǎn)，以下的發(fā)展方向值得進(jìn)一步探索。

1) 高階精度Godunov型大時(shí)間步長(zhǎng)格式

6.2節(jié)的分析表明Godunov型大時(shí)間步長(zhǎng)格式所取的時(shí)間步長(zhǎng)越大，則格式的分辨率越高，其在大時(shí)間步長(zhǎng)下的分辨率甚至高過(guò)了常規(guī)的高階精度格式。但是發(fā)展高階精度的大時(shí)間步長(zhǎng)格式仍是必要的，特別是在提高黏性分辨率方面，有更大的潛力。但是構(gòu)造高階格式并不容易，如果采用比分段常數(shù)更高精度的高階重構(gòu)，如分段線(xiàn)性重構(gòu)，則會(huì)面臨所謂的廣義Riemann問(wèn)題的求解，目前對(duì)于廣義Riemann問(wèn)題尚未有有效的精確或近似解法，這將引起疊波法的使用困難。為了避免該難題，繼續(xù)采用分段常數(shù)重構(gòu)相對(duì)容易，本文作者[85]曾結(jié)合加權(quán)平均狀態(tài)方法和疊波法，構(gòu)造了具有二階精度的Godunov型大時(shí)間步長(zhǎng)格式，但目前CFL數(shù)僅能達(dá)到2.0。如何克服廣義Riemann問(wèn)題帶來(lái)的疊波法應(yīng)用困難，應(yīng)該是下一步工作要探索的重點(diǎn)。

2) 引入自適應(yīng)參數(shù)的大時(shí)間步長(zhǎng)格式

目前基于疊波法所發(fā)展的大時(shí)間步長(zhǎng)格式具有無(wú)自由參數(shù)的優(yōu)點(diǎn)，但是從數(shù)值算例還是可以發(fā)現(xiàn)該格式在大CFL數(shù)下，在激波波后仍有微幅數(shù)值振蕩，這可能是數(shù)值耗散不足造成的。Lindqvist等[97]結(jié)合了LTS-Roe和LTS-LxF兩種格式的優(yōu)劣，提出了混合型的LTS-RoeLxF(β)格式；Prebeg等[98]提出了LTS-HLLE格式與LTS-Lax-Friedrichs格式的混合型LTS-HLLELxF(β)格式。其目的都是通過(guò)LTS-LxF引入更多的耗散，但是其引入的混合參數(shù)β為固定值，對(duì)于復(fù)雜流動(dòng)計(jì)算難以得到表現(xiàn)優(yōu)異的全場(chǎng)統(tǒng)一的混合參數(shù)。故而，發(fā)展自適應(yīng)形式的混合參數(shù)，使其根據(jù)流場(chǎng)結(jié)構(gòu)的特點(diǎn)自適應(yīng)地調(diào)節(jié)，是一種優(yōu)化的策略。對(duì)于疊波法而言，該思路也是可以嘗試的，在疊波過(guò)程引入自適應(yīng)調(diào)節(jié)參數(shù)，也可方便調(diào)節(jié)格式的耗散特性。

3) 方程源項(xiàng)的大時(shí)間步長(zhǎng)處理

在航空飛行器繞流的數(shù)值模擬中，黏性是不可忽略的重要因素。當(dāng)前飛行器空氣動(dòng)力學(xué)研究普遍采用RANS湍流模型來(lái)模擬黏性湍流流動(dòng)，湍流模型中一般存在描述湍流生成和耗散的源項(xiàng)。在考慮化學(xué)非平衡的高超聲速流動(dòng)或燃燒效應(yīng)的流動(dòng)中也存在描述化學(xué)反應(yīng)效應(yīng)的源項(xiàng)。因此，開(kāi)展帶源項(xiàng)的雙曲型守恒律的大時(shí)間步長(zhǎng)格式研究具有重要的實(shí)際工程意義。Murillo[111]、Hernandez[113]等在淺水波方程的求解中發(fā)展了處理渠底坡度和摩擦損失等帶來(lái)的源項(xiàng)處理方法，可為空氣動(dòng)力學(xué)領(lǐng)域的控制方程求解所借鑒。但空氣動(dòng)力學(xué)領(lǐng)域控制方程的源項(xiàng)一般非線(xiàn)性更強(qiáng)，故處理難度更大，是今后需要重點(diǎn)探討的方向之一。

4) 與自適應(yīng)網(wǎng)格技術(shù)結(jié)合的真正多維大時(shí)間步長(zhǎng)格式

自適應(yīng)網(wǎng)格技術(shù)是空氣動(dòng)力學(xué)數(shù)值模擬的重要手段，可隨著流場(chǎng)結(jié)構(gòu)的變化動(dòng)態(tài)調(diào)整網(wǎng)格的分布，是高精度求解復(fù)雜流場(chǎng)問(wèn)題的高效手段。目前的大時(shí)間步長(zhǎng)格式均基于一維情形構(gòu)造，然后借助維數(shù)分裂法來(lái)實(shí)現(xiàn)多維問(wèn)題推廣。但實(shí)際上Godunov型大時(shí)間步長(zhǎng)格式的疊波法具備推廣至多維情形的潛力，并且疊波法具有先天追蹤波系信號(hào)的特點(diǎn)，在采用自適應(yīng)網(wǎng)格技術(shù)方面具有獨(dú)特優(yōu)勢(shì)。如能將疊波法和自適應(yīng)網(wǎng)格技術(shù)相結(jié)合，則可能構(gòu)造出顯式的自適應(yīng)網(wǎng)格大時(shí)間步長(zhǎng)格式，具備更高效求解空氣動(dòng)力學(xué)問(wèn)題的能力，也是值得探索的方向。

致 謝

感謝清華大學(xué)陳海昕教授，本文是在他的鼓勵(lì)下完成的。