Godunov型顯式大時間步長格式研究進(jìn)展

錢戰(zhàn)森

1. 中國航空工業(yè)空氣動力研究院，沈陽 110034 2. 高速高雷諾數(shù)氣動力航空科技重點實驗室，沈陽 110034

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

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

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

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

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

1 顯式大時間步長格式簡介

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

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

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

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

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

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

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

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

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

Godunov于1959年提出的間斷分解法被認(rèn)為是最早的激波捕捉格式[31]，即通過Riemann問題間斷分解的計算來構(gòu)造差分格式，這類格式統(tǒng)稱作Godunov型格式，這個思路直到今天依然是CFD差分格式研究的熱點，并且一直處于不斷完善和發(fā)展之中。近30年來，Godunov型格式的研究和發(fā)展主要有兩個方向：一是對格式的高階推廣，提高格式的精度以提高對各種間斷的分辨率，代表性工作包括Boris和Book[32]、Harten等[33-34]通過對一階差分格式的數(shù)值通量進(jìn)行修正而獲得二階精度的格式，Jameson等[35]在中心差分格式(JST格式)基礎(chǔ)上采用人工黏性法獲得二階精度的格式，van Leer[36-39]及Collela和Woodward[40]通過對初始時刻的解采用分段線性或拋物線重構(gòu)函數(shù)代替原Gudonov格式中的分段常數(shù)重構(gòu)以獲得二階或三階精度的格式，Harten等[41]提出的基本無震蕩(ENO)格式及Liu[42]、Shu[43]等提出并改進(jìn)的高效WENO格式，通過采用隱式差分可以在較小的模板點上獲得更高的插值精度的緊致格式及其改進(jìn)形式[44-53]；二是對近似Riemann解的研究，用各種近似Riemann求解器來代替原始Godunov格式采用的Riemann問題精確解以減少計量，代表性工作包括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解算器?？傮w來看，對近似Riemann解的研究近年來活躍度不如高階推廣高，其中較有代表性的研究工作是針對多維問題發(fā)展的混合型近似Riemann解[76-79]。

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

1.2 顯式大時間步長格式的概念

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

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

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

1.3 大時間步長格式的分類

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

1.4 大時間步長格式的優(yōu)勢

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

1) 定常問題計算效率高

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

2) 可直接應(yīng)用于非定常問題計算

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

3) 對間斷分辨率高

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

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

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

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

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

2 Godunov型大時間步長格式的提出

2.1 Godunov型格式

考慮一維雙曲型守恒律組的初值問題：

(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)行分段積分。這種方法的實現(xiàn)較煩瑣，且時間步長受到CFL≤0.5的嚴(yán)格限制，以避免左右兩相鄰間斷波發(fā)生相互干擾。因此，在實際計算中很少采用此方法。

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

(3)

式中：

(4)

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

2.2 Godunov型大時間步長格式

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

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

(5)

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

(6)

的解可表達(dá)為

(7)

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

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

(8)

同樣，對于sa<0的情況也可作類似的處理。

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

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

(9)

(10)

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

對于守恒律系統(tǒng)，與單方程不同的是間斷分解比較復(fù)雜，每一間斷可以分解出多個波，如圖1所示?？扇约僭O(shè)不同波之間相互自由穿越，互不干擾。這個假設(shè)對非線性方程也屬于線化近似，在一定CFL數(shù)范圍內(nèi)該近似足夠準(zhǔn)確，當(dāng)CFL數(shù)非常大時，它會導(dǎo)致解的扭曲，此時可以通過加密網(wǎng)格來彌補(可以證明[80]，當(dāng)網(wǎng)格加密時，該方法對所有CFL數(shù)都是收斂的)。圖4給出了Euler方程(每個間斷包含三波束)的多波束疊加示意圖[86]，這時波束雖然變得復(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型大時間步長格式的改進(jìn)

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

3.1 時間線插值技術(shù)

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

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

(13)

式中：

(14)

變量u可以寫為如下的線性組合：

(15)

(16)

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

(17)

考慮界面上的等價Riemann問題，有

(18)

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

(19)

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

(20)

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

在CFL數(shù)超過1的條件下，xA點將落在xj單元左側(cè)或xj+1單元右側(cè)。因為界面xj-1/2或xj+3/2處可能有激波存在，特征線將是不可逆的，基于逆向追蹤法的式(20)將無法直接使用。此時，界面通量可以改寫成

(21)

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

圖5 時間線插值法示意圖[83]Fig.5 Schematic of time-line interpolation method

(22)

(23)

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

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

(24)

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

(25)

(26)

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

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

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

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

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

圖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ù)

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

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

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

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

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

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

1) LTS-Roe格式

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

(27)

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

(28)

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

(29a)

(29b)

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

(30)

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

(31a)

(31b)

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

(32a)

(32b)

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

(33)

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

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

2) LTS-HLL格式

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

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

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

(34a)

(34b)

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

LTS-HLL格式也可以改寫成疊波形式

(35)

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

(36)

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

(37)

式中：

S1=SL,S2=SC,S3=SR

(38)

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

(39)

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

(40)

4 Godunov型大時間步長格式的高階精度推廣方法

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

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

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

(41)

(42)

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

(43)

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

(44)

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

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

(45)

(46)

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

(47)

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

(48)

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

4.2 WAS型二階精度大時間步長Godunov格式

采用WAS方法的起點是引入狀態(tài)變量的加權(quán)平均式(42)，然后采用式(41)給出數(shù)值通量。對于Euler方程，式(42)中的變量狀態(tài)可以選擇守恒變量、原始變量或其他狀態(tài)參數(shù)，均不影響格式的守恒性。為了同原有的一階Godunov型大時間步長格式相容，文獻(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è)的兩個半?yún)^(qū)間值(U(j,1)和U(j,2))的平均值給出計算數(shù)值通量Fj+1/2的Uj+1/2，即

(50)

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

(51)

(52)

(53)

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

(54)

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

5 多維問題推廣

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

5.1 維數(shù)分裂

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

(55)

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

(56)

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

(57)

則可得到

(58)

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

5.2 對稱維數(shù)分裂

式(58)可以改寫成

(59)

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

(60)

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

(61)

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

圖12 簡單分裂和對稱分裂的對比[86]Fig.12 Comparison of simple split and symmetric split[86]

(62)

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

5.3 分片Riemann問題的引入

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

(63)

(64)

式中：

(65)

(66)

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

(67)

其初值記為

(68)

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

(69)

(70)

5.4 邊界條件處理

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

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

6 Godunov型大時間步長格式的性能分析

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

6.1 收斂特性

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

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

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

2) TVD性質(zhì)

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

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

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

3) 熵穩(wěn)定性

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

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

6.2 分辨率

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

1) 誤差分析

文獻(xiàn)[86]給出標(biāo)量雙曲型守恒律方程的大時間步長格式的截斷誤差系數(shù)為

(71)

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

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

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

3) 數(shù)值實驗驗證

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

6.3 計算效率

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

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

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

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

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

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

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

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

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

7.1 一維Sod激波管問題

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

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

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

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

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

7.2 二維翼型繞流

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

圖22 NACA0012翼型繞流計算結(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 三維機翼繞流

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

圖23 文獻(xiàn)[84]M6機翼繞流計算結(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機翼繞流計算結(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格式對高超聲速雙橢球繞流的計算結(jié)果，來流馬赫數(shù)Ma為4.0、攻角α為0°，計算結(jié)果如圖25所示。其指出大時間步長Godunov格式在不同CFL數(shù)下的計算結(jié)果和Haten-Yee TVD格式[118]所得到的結(jié)果相符，都能準(zhǔn)確地捕捉到橢球頭部的脫體弓形激波和兩橢球鑲嵌部位附近的二次激波，圖中給出了CFL數(shù)為1.0，3.0和5.0時雙橢球物面及流場對稱面上的等壓線分布和對稱面上的上下壁面壓力分布。

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

8 大時間步長格式研究的發(fā)展方向

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

1) 高階精度Godunov型大時間步長格式

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

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

目前基于疊波法所發(fā)展的大時間步長格式具有無自由參數(shù)的優(yōu)點，但是從數(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(β)格式。其目的都是通過LTS-LxF引入更多的耗散，但是其引入的混合參數(shù)β為固定值，對于復(fù)雜流動計算難以得到表現(xiàn)優(yōu)異的全場統(tǒng)一的混合參數(shù)。故而，發(fā)展自適應(yīng)形式的混合參數(shù)，使其根據(jù)流場結(jié)構(gòu)的特點自適應(yīng)地調(diào)節(jié)，是一種優(yōu)化的策略。對于疊波法而言，該思路也是可以嘗試的，在疊波過程引入自適應(yīng)調(diào)節(jié)參數(shù)，也可方便調(diào)節(jié)格式的耗散特性。

3) 方程源項的大時間步長處理

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

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

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

致 謝

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