任玉新，王乾，潘建華，章雨思，黃乾旻

清華大學(xué) 航天航空學(xué)院，北京 100084

湍流與轉(zhuǎn)捩(尤其是高超湍流與轉(zhuǎn)捩)、湍流分離流、激波/邊界層干擾、氣動聲學(xué)、湍流燃燒等問題是飛行器設(shè)計(jì)中的關(guān)鍵科學(xué)問題。這類問題具有多尺度、非線性、非定常、存在流場間斷等特點(diǎn)，對高可信度數(shù)值模擬提出了很高的要求。以湍流模擬為例，以目前計(jì)算機(jī)的能力已可以很好地滿足RANS(Reynolds-Averaged Navier-Stokes)模擬的需求，也可以開展部分RANS/LES(Large Eddy Simulation)混合模擬的工作；但還遠(yuǎn)遠(yuǎn)不能滿足LES和DNS(Direct Numerical Simulation) 的需求。未來先進(jìn)飛行器將具有復(fù)雜構(gòu)型、高機(jī)動、先進(jìn)動力系統(tǒng)、多變飛行環(huán)境等特征，對數(shù)值模擬的要求更高，RANS模擬在很多情況下不能滿足需求，預(yù)計(jì)LES將成為研究復(fù)雜湍流的主要手段。此時(shí)，提高計(jì)算效率、降低計(jì)算成本將成為計(jì)算流體力學(xué)(Computational Fluid Dynamics，CFD)發(fā)展的驅(qū)動力之一[1]。

計(jì)算成本包括時(shí)間成本和費(fèi)用成本(二者相關(guān)，但其關(guān)系依賴于計(jì)算機(jī)架構(gòu)及模擬策略)。其中時(shí)間成本可通過發(fā)展高速超算設(shè)備控制，而費(fèi)用成本則與總計(jì)算量有關(guān)。即使未來計(jì)算機(jī)的發(fā)展使高雷諾數(shù)LES和DNS成為可能，但按現(xiàn)在超算收費(fèi)標(biāo)準(zhǔn)，其費(fèi)用仍然難以承受。發(fā)展高精度格式的學(xué)術(shù)意義在于有可能實(shí)現(xiàn)在一定計(jì)算誤差的前提下使計(jì)算量更小、計(jì)算時(shí)間更短；或在同樣的計(jì)算機(jī)機(jī)時(shí)消耗下得到的數(shù)值解精度更高。因此，發(fā)展高精度格式是提高計(jì)算效率、降低計(jì)算費(fèi)用成本的有效途徑之一。為進(jìn)一步闡明這一判斷，定義數(shù)值方法的計(jì)算效率為

(1)

式中：ε為在某個(gè)特定網(wǎng)格下數(shù)值解與真值的誤差；τ為開展相應(yīng)計(jì)算所用時(shí)間；h為相應(yīng)網(wǎng)格的特征尺度。顯然，計(jì)算效率是依賴具體問題和計(jì)算網(wǎng)格的。對于k階精度格式，當(dāng)網(wǎng)格尺度變?yōu)樵瓉淼?/2，則誤差變?yōu)樵瓉淼?/2k。同時(shí)，網(wǎng)格尺度變?yōu)樵瓉淼?/2后，三維問題網(wǎng)格數(shù)增加到原來的8倍?？紤]顯式格式時(shí)間步長也要減半，其總計(jì)算量至少增加到原來的16即24倍。對于隱式格式，網(wǎng)格增多后，其迭代收斂速度通常會變慢，計(jì)算量至少增加到原來的16倍也是比較合理的估計(jì)。因此有ε(h/2)=ε(h)/2k，τ(h/2)=24τ(h)，即η(h/2)=2k-4η(h)。

顯然，對小于四階精度的格式，其計(jì)算效率隨網(wǎng)格加密而下降。例如，二階格式當(dāng)網(wǎng)格間距變?yōu)樵瓉淼?/2時(shí)，計(jì)算效率變?yōu)樵瓉淼?/4。當(dāng)進(jìn)行RANS或無黏流計(jì)算時(shí)，對網(wǎng)格數(shù)的要求相對較低，二階格式簡單、易行，可能是較好的選擇；當(dāng)開展LES或者DNS計(jì)算時(shí)，二階格式的計(jì)算效率會顯著下降，發(fā)展高精度格式的必要性得到凸顯。

由于非結(jié)構(gòu)網(wǎng)格生成具有靈活、高效、自動化程度高的優(yōu)點(diǎn)，基于非結(jié)構(gòu)網(wǎng)格的數(shù)值方法在CFD軟件中占據(jù)主流[2]?？梢灶A(yù)見，非結(jié)構(gòu)網(wǎng)格高精度數(shù)值方法也將是新一代CFD軟件的主導(dǎo)算法[2]。因此，非結(jié)構(gòu)網(wǎng)格高精度數(shù)值方法最近得到了高度重視，已經(jīng)發(fā)展出很多有效算法。其中，高精度有限體積方法是最早得到發(fā)展的方法，為二階精度有限體積方法的直接發(fā)展。為實(shí)現(xiàn)高階精度，必須在每個(gè)控制體(單元)內(nèi)部獲得物理量的高階分布,這個(gè)分布通常用多項(xiàng)式表示。由于有限體積方法直接求解的是守恒變量在單元上的平均值,為確定多項(xiàng)式的待定系數(shù)，只依靠本單元的信息是不夠的，必須使用本單元及其足夠數(shù)量相鄰單元的平均值反推多項(xiàng)式的待定系數(shù),這個(gè)過程稱為重構(gòu)(Reconstruction)[3]。參與重構(gòu)的所有單元構(gòu)成的集合稱為重構(gòu)的模板(Stencil)。即使獲得了高階的重構(gòu)多項(xiàng)式，通常也不能直接用于計(jì)算包含激波間斷的可壓縮流動，因?yàn)闀陂g斷附近產(chǎn)生非物理的數(shù)值振蕩。為抑制在間斷附近的振蕩，必須采用合理的限制措施。限制過程是某種非線性濾波運(yùn)算，其中的非線性平均算子稱為限制器(Limiter)[4]。重構(gòu)和限制過程是高精度有限體積方法的關(guān)鍵難點(diǎn)問題。

在20世紀(jì)90年代，Barth和Frederickson[5]提出了基于最小二乘法數(shù)據(jù)重構(gòu)的k-exact高精度有限體積方法。此方法的優(yōu)點(diǎn)是可以較簡單、系統(tǒng)地確定重構(gòu)模板，每個(gè)單元只需進(jìn)行一次重構(gòu)；其缺點(diǎn)是沒有很好地解決激波捕捉的梯度限制問題，對激波的計(jì)算能力較差。通常采用的Barth[6]及Venkatakrishnan[7]限制器耗散很大，對激波和光滑流動結(jié)構(gòu)的分辨能力都較差。

WENO(Weighted Essentially Non-Oscillatory)方法[8-9]是構(gòu)造非結(jié)構(gòu)網(wǎng)格高精度有限體積格式的有效途徑。WENO方法通過幾個(gè)候選重構(gòu)多項(xiàng)式的非線性加權(quán)平均得到最終的重構(gòu)多項(xiàng)式。此方法把重構(gòu)和限制過程有機(jī)地結(jié)合起來，因此具有很好的激波捕捉能力，但其代價(jià)是需要確定多個(gè)重構(gòu)模板、進(jìn)行多次重構(gòu)，因而算法復(fù)雜、計(jì)算量很大。

在高階重構(gòu)過程中需要很大的重構(gòu)模板，是高精度有限體積方法始終無法回避的問題。造成重構(gòu)模板巨大的本質(zhì)原因是對于某一個(gè)單元的多項(xiàng)式重構(gòu)，其重構(gòu)模板中的每一個(gè)單元只有一個(gè)可以利用的自由度，即物理量的單元平均值。如在k-exact重構(gòu)中，一個(gè)模板單元只能提供一個(gè)線性方程。為使得k-exact重構(gòu)穩(wěn)定，一般模板中單元數(shù)目取值應(yīng)顯著大于待定系數(shù)數(shù)目[10]，并使用最小二乘方法求解重構(gòu)方程組。因此隨著重構(gòu)多項(xiàng)式階次的提高，待定系數(shù)數(shù)目迅速增大，重構(gòu)模板也相應(yīng)地變得巨大。重構(gòu)模板巨大會導(dǎo)致一系列影響計(jì)算效率和魯棒性的問題，如要利用較遠(yuǎn)的模板單元而使計(jì)算緩存利用率低、分區(qū)交界面上要傳遞的數(shù)據(jù)多而使并行計(jì)算效率降低、物理邊界處模板偏心化而加大重構(gòu)奇異性等。另外，對于程序設(shè)計(jì)，重構(gòu)模板巨大的一個(gè)弊端是其不確定性。在搜索模板時(shí)，常常是逐層向外擴(kuò)展，但經(jīng)常出現(xiàn)某一層模板單元太多的情況，全部采用則計(jì)算量和存儲量增加過多，只用部分單元的話又需要設(shè)計(jì)取舍的準(zhǔn)則，這給程序設(shè)計(jì)帶來了很大不便。

為克服高階有限體積格式的這些缺點(diǎn)，最近一些新型的、適用于非結(jié)構(gòu)網(wǎng)格的高精度算法得到了高度重視。它們和有限體積方法的最著區(qū)別是在確定近似解的分片光滑多項(xiàng)式過程中，主要不是從外部(相鄰單元)獲取信息，而是試圖從單元內(nèi)部得到更多的額外信息，從而保證格式模板的緊湊。為闡述方便，統(tǒng)稱這些方法為非結(jié)構(gòu)網(wǎng)格上的亞網(wǎng)格自由度高階格式。其代表性的工作有間斷迦遼金有限元方法(Discontinuous Galerkin method，DG)[11-13]、譜體積方法(Spectral Volume method，SV)[14]、譜差分方法(Spectral Difference method，SD)[15]。最近，Huynh[16]和王志堅(jiān)[17]等提出了FR(Flux Reconstruction)/CPR(Correction Procedure via Reconstruction)方法，可以把SD、DG等方法歸納到一個(gè)統(tǒng)一框架中。此外，亞網(wǎng)格自由度高階算法和有限體積(Finite Volume，F(xiàn)V)方法重構(gòu)過程的結(jié)合也得到了重視，發(fā)展了PNPM方法[18-19]、rDG(reconstructed DG)方法[20-21]、DG/FV混合算法[22]、PNPM-CPR 方法[23]、multi-moment方法[24]等。

亞網(wǎng)格自由度高階算法的共同優(yōu)點(diǎn)是計(jì)算模板緊湊。這個(gè)特點(diǎn)使得這些方法在大規(guī)模計(jì)算中具有內(nèi)存尋址消耗小、對計(jì)算機(jī)緩存的利用率高、并行計(jì)算的效率較高、與一些新HPC(High Performance Computing)體系結(jié)構(gòu)(如基于GPU(Graphics Processing Unit)的并行計(jì)算)的相容性較好等優(yōu)點(diǎn)，顯著優(yōu)于傳統(tǒng)高階有限體積格式。此外，這些方法對光滑流場的計(jì)算精度和分辨率很高。然而，這些方法還存在一些尚未完全解決的問題，主要有計(jì)算包含激波間斷的流場時(shí)缺乏高效、保持精度的限制器或非線性人工黏性方案；計(jì)算含強(qiáng)激波、強(qiáng)膨脹波和幾何突變的流場時(shí)魯棒性較差；隱式格式規(guī)模很大，需要盡量精確地計(jì)算雅可比矩陣并加以存儲，對存儲量有很高要求，目前迫切需要發(fā)展高效、低存儲的隱式時(shí)間積分技術(shù)；高階格式需要基于高階網(wǎng)格，目前的網(wǎng)格生成軟件還不能完全滿足需求。由于這些原因，基于這些方法的商業(yè)化CFD軟件發(fā)展還相當(dāng)滯后。

相反，有限體積方法發(fā)展相對成熟，在保精度限制器、高效隱式積分等方面已有較好的解決方案，對計(jì)算量和存儲量的需求也相對較低。如果能夠解決重構(gòu)模板巨大這一瓶頸問題，應(yīng)該具有很好的發(fā)展前景。發(fā)展基于緊致模板的重構(gòu)方法是近期工作的出發(fā)點(diǎn)。

所謂緊致模板，即只包含當(dāng)前單元及其有公共面(三維)或公共邊(二維)的單元的模板。通過仔細(xì)研究有限體積方法的重構(gòu)過程發(fā)現(xiàn)，由于在重構(gòu)模板中的每一個(gè)單元上只能提供一個(gè)平均值信息，要獲得高階重構(gòu)多項(xiàng)式，必然要使用很大的模板或者提供額外的信息。為解決這一難題，提出操作緊致性這一重要概念。所謂操作緊致性，即不要求實(shí)際上采用的模板是緊致的，只要求在數(shù)值求解過程中所有的數(shù)據(jù)結(jié)構(gòu)和計(jì)算操作都只用到緊致模板中的信息。操作緊致性的概念為實(shí)現(xiàn)基于緊致模板的重構(gòu)提供了廣闊的空間?；谶@一概念，提出了幾種緊致重構(gòu)方法。這些方法可以分為兩類，即隱式重構(gòu)和顯式重構(gòu)。

在隱式重構(gòu)方面提出了緊致最小二乘重構(gòu)[25-26]和變分重構(gòu)[27]方法。目前廣泛使用的高精度重構(gòu)方法對于每一個(gè)模板單元都只利用了一個(gè)信息，即物理量的單元平均值。這樣一來要確定單元上的重構(gòu)多項(xiàng)式系數(shù)，就需要個(gè)數(shù)不少于重構(gòu)多項(xiàng)式系數(shù)的模板單元，這就是現(xiàn)有重構(gòu)方案模板巨大的根本原因。為發(fā)展緊致高精度重構(gòu)，必須突破每個(gè)模板單元只貢獻(xiàn)一個(gè)信息的思路，采取“挖掘更多信息，使用更小模板”的策略。為在只包括面相鄰單元的模板上實(shí)施高階重構(gòu)，不僅要利用面相鄰單元的物理量平均值，還要利用面相鄰單元物理量的導(dǎo)數(shù)。由于面相鄰單元物理量的導(dǎo)數(shù)也是待定的，所以和結(jié)構(gòu)網(wǎng)格上的緊致差分方法類似，這類重構(gòu)過程是隱式的。隱式重構(gòu)面臨的主要困難是需要求解大型稀疏矩陣對應(yīng)的線性方程組，因而計(jì)算量很大。通過兩個(gè)手段可以顯著提高計(jì)算效率：其一是通過迭代法求解重構(gòu)過程；其二是與隱式時(shí)間推進(jìn)相結(jié)合，將重構(gòu)的迭代和隱式格式迭代求解的過程同步進(jìn)行。這些方法可使基于隱式重構(gòu)的有限體積格式計(jì)算效率與隱式k-exact格式相當(dāng)。

在顯式重構(gòu)方面提出了多步重構(gòu)方法[28]。多步重構(gòu)通過不斷擴(kuò)展原始重構(gòu)關(guān)系，每一步一般可提高一階精度，從而逐步達(dá)到高階精度。這本質(zhì)上是一個(gè)逐步擴(kuò)大模板的過程，但在操作層面只需處理緊致模板的信息。顯式重構(gòu)的優(yōu)點(diǎn)是與顯式和隱式時(shí)間推進(jìn)方法相結(jié)合都可以達(dá)到較高計(jì)算效率；而隱式重構(gòu)只有采用隱式時(shí)間積分時(shí)才可能有較高計(jì)算效率。因此，顯式重構(gòu)比隱式重構(gòu)在時(shí)間推進(jìn)方面有更大的靈活性。

除重構(gòu)過程外，高精度有限體積方法還必須解決邊界條件、激波捕捉、時(shí)間推進(jìn)、高階網(wǎng)格計(jì)算等一系列問題。本文對此也將做簡要介紹。

本文首先概要介紹有限體積方法的基本理論，重點(diǎn)對常見的k-exact重構(gòu)進(jìn)行介紹，分析大模板問題的來源；然后詳細(xì)介紹3種緊致重構(gòu)方法，即緊致最小二乘重構(gòu)、變分重構(gòu)和多步重構(gòu)；接著進(jìn)一步介紹緊致重構(gòu)有限體積方法的實(shí)施、效率提升、邊界處理、激波捕捉、高階網(wǎng)格、時(shí)間推進(jìn)、方法性能等問題，并簡述近期的發(fā)展；最后給出結(jié)論。

1 有限體積方法基本原理

原則上，有限體積方法可求解雙曲型、拋物型、橢圓型及其他偏微分方程對應(yīng)的積分型方程。其中最為重要的應(yīng)用是求解流體力學(xué)中的守恒律方程。常見的Euler方程和Navier-Stokes方程都可以寫成守恒律。以二維三角形網(wǎng)格上求解守恒律的有限體積格式為例做進(jìn)一步介紹。一般的二維守恒型方程組可以寫為

(2)

式中：U為守恒變量，對于二維Euler方程和Navier-Stokes方程而言，為有4個(gè)分量的列向量；F和G為通量函數(shù)；t、x和y分別為時(shí)間變量及空間的兩個(gè)笛卡爾坐標(biāo)。一般而言，一階守恒律中F和G為U的函數(shù)(如Euler方程)；二階守恒律中F和G為U及其梯度的函數(shù)(如Navier-Stokes方程)。在圖1所示控制體Ωi上對式(2)積分，并利用Gauss公式得到

(3)

圖1 有限體積方法的控制體Fig.1 Control volume of finite volume method

外法向量；S為控制體界面面積。

式(3)即為有限體積方法求解的積分型方程。事實(shí)上，式(3)不一定依賴于式(2)，常?？梢詮幕臼睾愣芍苯訉?dǎo)出。

直接對式(3)進(jìn)行空間離散的有限體積方法稱為半離散型有限體積方法。另一種方案即所謂全離散有限體積方法，其基本方程可以通過對式(3) 在一個(gè)時(shí)間步內(nèi)積分得到，即

(4)

式中：n為數(shù)值計(jì)算的時(shí)間步；Δt為時(shí)間步長；tn為n時(shí)間步對應(yīng)的時(shí)刻。

無論是全離散還是半離散有限體積方法，其因變量均為物理量在控制體上的平均值。而式(3)、式(4)的空間項(xiàng)涉及物理量在界面上的分布及其演化。因此，為用數(shù)值方法求解式(3)、式(4)，必須引入各種近似。這些近似構(gòu)成了有限體積方法的基本要素。為更好地解釋這些要素，將有限體積方法抽象地寫為如式(5)和式(6)的形式：

全離散有限體積格式為

(5)

半離散有限體積格式為

(6)

式中：MS為控制體邊數(shù)，對于圖1所示三角形控制體，MS=3；RΩ、EΔt(或?yàn)镋0)分別為重構(gòu)和演化積分算子；Im為沿控制體第m個(gè)邊的數(shù)值通量積分算子；Si為重構(gòu)模板，其具體含義將在1.1節(jié)中介紹。Im的構(gòu)造比較直接，其意義為

(7)

1.1 重 構(gòu)

有限體積方法能夠得到的是物理量在各個(gè)控制體上的平均值，而有限體積方法的核心要素——數(shù)值通量的計(jì)算涉及物理量在界面上的分布及其演化。為計(jì)算數(shù)值通量，需要利用物理量在各個(gè)控制體上的平均值構(gòu)造出物理量的分布。這一過程稱為重構(gòu)，式(5)、式(6)中的RΩ稱為重構(gòu)算子。重構(gòu)函數(shù)是分別定義在各個(gè)控制體上的分片光滑函數(shù)，一般為多項(xiàng)式。

在有限體積方法的實(shí)施中，對守恒變量U的分量u逐一進(jìn)行重構(gòu)。控制體Ωi上，守恒變量的重構(gòu)函數(shù)記為Ui(x,t)，對于二維情形，x=xi+yj。只考慮Ui(x,t)的一個(gè)分量ui(x,t)的重構(gòu)問題。由于整個(gè)重構(gòu)過程是在同一個(gè)時(shí)刻進(jìn)行的，因此省略時(shí)間變量t，將重構(gòu)多項(xiàng)式記為ui(x)。一般而言，基于k次多項(xiàng)式重構(gòu)的有限體積方法，其空間精度為k+1階。其中零階重構(gòu)最為簡單，假定每個(gè)控制體內(nèi)ui(x)為常數(shù)(零次多項(xiàng)式)，即為單元平均值。Barth和Frederickson[5]提出的k-exact重構(gòu)方法可以進(jìn)行任意次數(shù)的多項(xiàng)式重構(gòu)，從而可以構(gòu)造任意階空間精度的有限體積方法；他們還提出了對重構(gòu)過程的3個(gè)基本要求：

2)k-exact性質(zhì)，即k次多項(xiàng)式重構(gòu)可以準(zhǔn)確復(fù)現(xiàn)次數(shù)不高于k的多項(xiàng)式精確解。

3) 緊湊性。有限體積方法只能提供每個(gè)控制體物理量的平均值，而當(dāng)k>0時(shí)，每個(gè)單元重構(gòu)多項(xiàng)式的待定系數(shù)個(gè)數(shù)大于1。因此，在Ωi單元進(jìn)行重構(gòu)時(shí)，必須利用其相鄰的多個(gè)控制體上平均值的信息。參與Ωi單元重構(gòu)的單元集合稱為Ωi單元的重構(gòu)模板。所謂緊湊性，即要求重構(gòu)模板中的單元應(yīng)優(yōu)先選取鄰近Ωi的控制體。

k-exact重構(gòu)的主要困難在于隨重構(gòu)多項(xiàng)式次數(shù)的提高，重構(gòu)模板中所需控制體的個(gè)數(shù)顯著增加，這就是所謂大模板問題。大模板問題增加了編程復(fù)雜度，降低了并行計(jì)算的效率，且難以以高精度實(shí)施邊界條件。

為更清楚地說明大模板問題的來源，簡要回顧k-exact重構(gòu)。設(shè)當(dāng)前時(shí)刻計(jì)算域上所有單元的平均值已知，重構(gòu)算法需要給出u在計(jì)算域上的一個(gè)分片連續(xù)的k次多項(xiàng)式分布，即在任意一個(gè)單元Ωi上，用k次多項(xiàng)式式(8)逼近u在單元i中的分布：

(8)

(9)

式中：p和q均為基函數(shù)的次數(shù)；dΩ為體積微元；xi和yi為控制體的形心坐標(biāo)；Δxi和Δyi為用來對基函數(shù)進(jìn)行無量綱化的特征尺度，可以有多種選擇，這里采用Luo等[29]的方案：

(10)

式中：xmax、ymax和xmin、ymin分別單元i的x、y坐標(biāo)最大值和最小值。基函數(shù)的無量綱化能夠克服重構(gòu)矩陣條件數(shù)隨網(wǎng)格加密而增大的問題[9]。用k=2的重構(gòu)多項(xiàng)式舉例說明指標(biāo)l與p、q之間的關(guān)系：

(11)

(12)

給定階次k時(shí)，重構(gòu)多項(xiàng)式可用式(8)～式(10) 表示。重構(gòu)基于事先選定的重構(gòu)模板，例圖2是在二維三角形單元上進(jìn)行k=3的三次多項(xiàng)式重構(gòu)時(shí)，Ωi單元的一種可能的重構(gòu)模板，記為Si。一般來講，模板中單元的個(gè)數(shù)ni應(yīng)滿足ni≥Nc(k)。

圖2 二維三角形單元的三次k-exact重構(gòu)模板示意圖Fig.2 Schematic diagram of third order k-exact reconstruction on triangular cells in 2D

k-exact重構(gòu)的基本原理是利用k-exact性質(zhì)導(dǎo)出重構(gòu)關(guān)系。對于充分光滑的函數(shù)u(x)(精確解)，在Ωi的不小于重構(gòu)模板Si的鄰域內(nèi)，總是可以找到一個(gè)k次多項(xiàng)式p(x)，使其逼近精確解u(x) 到k+1階精度，即

p(x)=u(x)+O(hk+1)

(13)

式中：O(hk+1)為無窮小量。

(14)

與式(14)相類比，k-exact重構(gòu)的重構(gòu)關(guān)系可以取

(15)

將式(8)代入式(15)，得

(16)

(17)

在實(shí)際計(jì)算中，由于重構(gòu)多項(xiàng)式采用式(9)的零均值基，當(dāng)前單元的重構(gòu)關(guān)系即守恒條件式(12) 是自動滿足的，因此不能再作為一個(gè)重構(gòu)關(guān)系。所以重構(gòu)關(guān)系式(16)可以進(jìn)一步明確為

(18)

根據(jù)對k-exact重構(gòu)的介紹和分析，可以進(jìn)一步總結(jié)如下：

1) 重構(gòu)操作是針對充分光滑函數(shù)進(jìn)行的，這是式(13)成立的前提，在這一前提下才能估計(jì)重構(gòu)精度。所述緊致重構(gòu)也是針對充分光滑函數(shù)的。由于可壓縮流動可能存在流場間斷，在有限體積方法中采用除零階重構(gòu)以外的較高精度重構(gòu)時(shí)，數(shù)值解會在間斷附近產(chǎn)生非物理振蕩。由于間斷總是局限在某些線或面上，流場大部分區(qū)域還是可以采用對光滑函數(shù)的重構(gòu)算法。在間斷附近，需要引入某種限制器對重構(gòu)多項(xiàng)式進(jìn)行局部修正。Van Leer的MUSCL(Monotonic Upwind centered Scheme for Conservation Laws)格式[4]是較早期的限制器構(gòu)造方案。限制器的構(gòu)造還可以基于TVD(Total Variation Diminishing)條件[30]、極大值原理(如Barth限制器)[6]或者TVB(Total Variation Bounded)條件(如著名的ENO(Essentially Non-Oscillatory)和WENO格式[31]等。

2) 在求解方程組式(18)時(shí)，要求線性方程組系數(shù)矩陣非奇異，或在最小二乘意義下非奇異；而且矩陣條件數(shù)對重構(gòu)精度有直接影響。為滿足這些要求，在非結(jié)構(gòu)網(wǎng)格上一般采用最小二乘法求解式(18)。經(jīng)驗(yàn)表明，比較合適的模板中單元數(shù)量為待定系數(shù)數(shù)目的1.5～2.0倍。這就是k-exact重構(gòu)方法大模板問題的由來。按照這個(gè)要求，模板中單元個(gè)數(shù)隨重構(gòu)多項(xiàng)式次數(shù)提高迅速增長。對于三維的三次(k=3)多項(xiàng)式重構(gòu)，有19個(gè)待定系數(shù)，需要約30個(gè)模板單元，而四次(k=4)多項(xiàng)式重構(gòu)則需要約50個(gè)模板單元。模板巨大會造成一系列的問題:模板單元和中心單元的編號差距大，造成CPU緩存利用率低；并行計(jì)算時(shí)虛擬網(wǎng)格量大，進(jìn)程之間需要傳遞的數(shù)據(jù)量大，通信時(shí)間長，降低并行效率；在邊界處，由于需要較多的鄰單元，只能向計(jì)算域內(nèi)部擴(kuò)展模板，造成模板的偏心化，加大重構(gòu)的奇異性。因此，模板巨大會降低格式的效率和魯棒性，已成為高精度有限體積方法發(fā)展的瓶頸問題。在1.2節(jié)將會給出解決這一問題的幾種有效方案。

1.2 演 化

雙曲型守恒律(如Euler方程)的數(shù)值通量或Navier-Stokes方程的無黏部分?jǐn)?shù)值通量計(jì)算中，演化過程的研究是計(jì)算流體力學(xué)發(fā)展過程中非常重要的研究內(nèi)容。著名的Godunov格式[32]是早期的代表性工作。在采用零階重構(gòu)的情況下，Godunov格式的演化過程通過求解黎曼問題實(shí)現(xiàn)。黎曼問題為雙曲型守恒律在分段常數(shù)初始條件下的初值問題，對于一維Euler方程存在解析解。由于精確求解黎曼問題計(jì)算量較大，還發(fā)展了多種黎曼問題的近似解法，代表性的有Roe[33]、Osher[34]和Toro[35]等的近似黎曼求解器,這些黎曼求解器可以通過維數(shù)分裂技術(shù)推廣到多維情況。

對于全離散的有限體積格式，由于要計(jì)算數(shù)值通量在一個(gè)時(shí)間步內(nèi)的演化，其演化過程依賴于重構(gòu)運(yùn)算。例如，Godunov格式由于其中的黎曼求解器只能處理分段常數(shù)的初值，即只能作用于零階重構(gòu)，因而格式只有一階精度。如果要達(dá)到二階精度，則需求解分段線性函數(shù)初值的黎曼問題，稱為廣義黎曼問題。廣義黎曼問題比常規(guī)黎曼問題更為復(fù)雜，一般沒有解析解。Matania和Joseph[36]研究了廣義黎曼問題的求解，使得全離散有限體積格式達(dá)到了二階精度。近年來，Toro和Titarev[37]研究了對應(yīng)更高階重構(gòu)的線性化廣義黎曼問題或?qū)?shù)黎曼問題，發(fā)展了更高階精度的全離散有限體積格式——ADER格式。徐昆[38]提出的氣體動理學(xué)格式采用全離散有限體積方法；其演化過程基于簡化的玻爾茲曼方程，且可以和任意重構(gòu)算法相結(jié)合。李杰權(quán)和杜知方[39]發(fā)展的二步四階格式也可以用于全離散有限體積方法。

半離散有限體積格式由于只需計(jì)算數(shù)值通量的無窮小時(shí)間演化，其演化算子和重構(gòu)算子是解耦的，因此Godunov[32]、Roe[33]的黎曼求解器均可以和任意重構(gòu)算法結(jié)合，構(gòu)造出任意空間精度的有限體積格式，而無需求解廣義黎曼問題。Van Leer[4]首先認(rèn)識到半離散有限體積格式的這一特性，從而把Godunov格式推廣到高階精度。

可以證明，基于黎曼求解器的Godunov類格式是一類迎風(fēng)型格式。對于半離散的有限體積格式，除Godunov類格式外還發(fā)展了其他多種迎風(fēng)型演化算法，如矢通量分裂(Flux Vector Splitting，F(xiàn)VS)格式[40]、AUSM(Advection Upstream Splitting Method)格式[41]等。此外，中心型演化算法也得到了發(fā)展。最簡單的中心型演化算法是直接對數(shù)值求積點(diǎn)兩側(cè)物理量取算術(shù)平均，但這種方法一般不穩(wěn)定；Jameson[42]通過添加混合二階和四階人工黏性項(xiàng)得到了穩(wěn)定的有限體積格式；Kurganov和Tadmor[43]通過在控制體界面附近構(gòu)造輔助控制體、在輔助控制體上求解原方程的方法構(gòu)造演化算子，也得到了穩(wěn)定的有限體積格式。

1.3 投影及時(shí)間推進(jìn)

通過重構(gòu)及演化過程得到了控制體界面處逐點(diǎn)的數(shù)值通量。顯然，數(shù)值通量的計(jì)算基于重構(gòu)函數(shù)空間，而有限體積方法的因變量是物理量的單元平均值。從重構(gòu)函數(shù)空間轉(zhuǎn)換到單元平均值的過程稱為投影，其操作為如式(7)所示的通量積分過程。對于全離散有限體積方法，通過投影過程可直接得到下一時(shí)間步物理量的單元平均值。而對于半離散有限體積方法，投影過程的結(jié)果是關(guān)于單元平均值時(shí)間導(dǎo)數(shù)的常微分方程組。為計(jì)算下一時(shí)間步物理量的單元平均值，還需確定具體的時(shí)間推進(jìn)格式。一般來說，所有求解常微分方程組的數(shù)值方法都可作為時(shí)間推進(jìn)格式，如各種Runge-Kutta方法、多步方法及隱式方法等。應(yīng)該指出，有限體積方法的因變量是物理量的單元平均值，但其最終數(shù)值解應(yīng)視為經(jīng)重構(gòu)過程得到的物理量在控制體內(nèi)部的分布函數(shù)；這一點(diǎn)對于正確理解高精度有限體積方法尤為重要。

2 緊致重構(gòu)

以半離散有限體積方法為例，介紹最近發(fā)展的非結(jié)構(gòu)網(wǎng)格高精度有限體積方法的緊致重構(gòu)算法。重構(gòu)問題的基本解法在第1節(jié)已做介紹。為在緊致模板(即圖1所示模板)上實(shí)現(xiàn)任意高階重構(gòu)，必須有效利用重構(gòu)過程的信息。對于高階重構(gòu)，重構(gòu)多項(xiàng)式中待定系數(shù)的數(shù)量很大，但已知信息僅為單元均值，從而造成了所謂大模板問題。操作緊致性是為解決這個(gè)問題提出的新思路。從數(shù)值計(jì)算角度看，操作緊致性對于數(shù)據(jù)結(jié)構(gòu)設(shè)定、程序設(shè)計(jì)、并行計(jì)算而言是和實(shí)際模板是緊致的情況完全一樣的。事實(shí)上，絕對緊致的格式是不存在的。例如DG等方法即使在一個(gè)時(shí)間步內(nèi)只用到緊致模板的信息，但還會隱含地用到更早的時(shí)間步內(nèi)、更大范圍內(nèi)網(wǎng)格上的信息。緊致重構(gòu)方法在一個(gè)時(shí)間步內(nèi)就會隱含地用到緊致模板以外的信息，但所有操作都是局部的，只用到緊致模板內(nèi)的信息。操作緊致性這一概念創(chuàng)新極大擴(kuò)展了構(gòu)造重構(gòu)算法的多樣性，使發(fā)展高階緊致重構(gòu)成為可能。同時(shí)，也使高精度緊致有限體積方法可以采用二階格式的數(shù)據(jù)結(jié)構(gòu)、計(jì)算操作和存儲策略，實(shí)現(xiàn)高精度計(jì)算，從而有效解決了高階有限體積方法模板過大這一瓶頸問題。下文將逐一介紹提出的緊致重構(gòu)方法，簡單起見，只討論一維或二維緊致重構(gòu)。重構(gòu)均基于緊致模板，即圖1中的單元Ωi及其面相鄰單元Ωj1、Ωj2、Ωj3。當(dāng)Ωj單元與Ωi相鄰時(shí)，則必為Ωj1、Ωj2、Ωj3之一。

2.1 緊致最小二乘重構(gòu)

在k-exact重構(gòu)中采用的重構(gòu)關(guān)系如式(18)所示，即

(19)

此時(shí)，每個(gè)模板中的單元只能提供一個(gè)關(guān)系，從而造成了所謂大模板問題。注意到由重構(gòu)的守恒性，式(19)也可以寫為

?Ωj∈Si,j≠i

(20)

即要求ui(x)在Ωj上的平均值與uj(x)在Ωj上的平均值相同。按照“挖掘更多信息，使用更小模板”的思路，對式(20)進(jìn)行擴(kuò)展，即不僅要求重構(gòu)多項(xiàng)式本身滿足式(20)，同時(shí)要求其各階導(dǎo)數(shù)也滿足類似關(guān)系，稱為導(dǎo)數(shù)守恒關(guān)系，從而有[25-26]

?Ωj∈Si,j≠i, 0≤m+n≤M≤k

(21)

式中：M為導(dǎo)數(shù)守恒關(guān)系中最高階導(dǎo)數(shù)的階數(shù)。取不同的m和n時(shí)，式(21)各方程的量綱不同，為此對這些方程進(jìn)行無量綱化，引入權(quán)函數(shù)wi,m+n=(whi)m+n，其中w為自由參數(shù)，用于控制各個(gè)階次的導(dǎo)數(shù)的守恒方程的權(quán)重；hi為單元特征長度尺度，一般取為控制體外接圓半徑。合理地調(diào)節(jié)w可以改善格式精度和計(jì)算效率，具體可見文獻(xiàn)[44]。考慮權(quán)函數(shù)后的導(dǎo)數(shù)守恒關(guān)系為

(22)

導(dǎo)數(shù)守恒關(guān)系中導(dǎo)數(shù)階數(shù)為0時(shí)，式(21)退化為式(20)。因此，導(dǎo)數(shù)守恒關(guān)系包含了比常規(guī)k-exact重構(gòu)更多的條件，從而可以縮小重構(gòu)模板。

將重構(gòu)多項(xiàng)式式(8)代入式(22)，得到線性方程組：

(23)

(24)

式中：矩陣和向量的定義為

(25)

當(dāng)將j取遍重構(gòu)模板中的面相鄰單元{j1,j2,j3}并對其對應(yīng)的線性方程組式(24)進(jìn)行組裝，得到單元i的重構(gòu)線性方程組：

(26)

式中：

(27)

在二維情況下，一般取導(dǎo)數(shù)守恒方程的最高階次M=k-1。這樣一來，重構(gòu)線性方程組式(26) 的方程個(gè)數(shù)為3[Nc(k-1)+1]，未知量個(gè)數(shù)為Nc(k)，可以推導(dǎo)出3[Nc(k-1)+1]-Nc(k)=k2-k+1>0，即方程個(gè)數(shù)總是大于未知數(shù)的個(gè)數(shù)，二者的數(shù)量對比見表1。所以重構(gòu)線性方程組式(26)是超定的。

表1 2～4階緊致重構(gòu)線性方程組的方程和未知量個(gè)數(shù)

用最小二乘法求解超定線性方程組式(26)，因此這種重構(gòu)方法稱為緊致最小二乘(Compact Least Squares，CLS)重構(gòu)方法。在一維情況下，全場單元的法方程組成一個(gè)塊三對角方程組，可以用直接法快速求解；但在多維情況下，非結(jié)構(gòu)網(wǎng)格單元編號是不連續(xù)的，中心單元和面相鄰單元的編號可能相差很大，因此全場單元的法方程組成一個(gè)大型稀疏線性方程組。全場的法方程聯(lián)立組成一個(gè)大型稀疏線性方程組，所以每個(gè)單元實(shí)際依賴的模板都是全場單元的集合。如用直接法求解，那么求解過程無論如何都不能做到緊致。因此，感興趣的是能夠保持重構(gòu)方程組求解過程緊致性的迭代方法。如第2節(jié)第1段所述，緊致是指求解過程的操作緊致性，即在求解過程的任何一個(gè)步驟只需使用面相鄰模板單元的信息。程序?qū)嵤r(shí)只需定義基于緊致模板的數(shù)據(jù)結(jié)構(gòu)，操作上緊致這一性質(zhì)已足以克服重構(gòu)模板巨大造成的緩存命中率低及并行計(jì)算時(shí)進(jìn)程間的數(shù)據(jù)交換量大等問題。

采用奇異值分解(Singular Value Decomposition，SVD)方法計(jì)算式(26)的最小二乘解，并采用Block Gauss-Seidel方法逐個(gè)單元進(jìn)行迭代求解以保持求解過程的操作緊致性。式(26)的Block Gauss-Seidel迭代公式為

(28)

2.2 變分重構(gòu)

變分重構(gòu)[26]是一類模板緊致的高精度重構(gòu)方法的總稱。這類重構(gòu)方法的原理是定義一個(gè)泛函，需要確定的重構(gòu)關(guān)系恰好使泛函取極小值。要求泛函重構(gòu)基于緊致模板，且具有k-exact性質(zhì)。為滿足這個(gè)要求，可以在任意兩個(gè)控制體的交界面上定義界面跳躍積分(Interfacial Jump Integration，IJI)。即對于一個(gè)單元交界面f，令其左、右兩側(cè)的單元編號分別是L、R，那么界面f上的IJI定義是

(29)

(30)

式中：Nf為計(jì)算域上的單元交界面的總數(shù)。事實(shí)上，當(dāng)所有控制體的平均值均由同一個(gè)次數(shù)不高于k的多項(xiàng)式計(jì)算，且每個(gè)單元的重構(gòu)多項(xiàng)式均取為這個(gè)多項(xiàng)式時(shí)，I取其極小值0。因此，這個(gè)泛函具有k-exact性質(zhì)，可以保證重構(gòu)的精度。

一般情況下，變分重構(gòu)用來求解重構(gòu)多項(xiàng)式待定系數(shù)的基本關(guān)系式是令泛函I取極小值導(dǎo)出的：

(31)

式中:N為控制體總數(shù)。

l=1,2,…,Nc(k);i=1,2,…,N

(32)

式中：dij為界面兩側(cè)單元形心之間的距離。

式(32)的矩陣形式為

(33)

式中：

(34)

可用精度足夠高的高斯數(shù)值積分精確地計(jì)算重構(gòu)矩陣表達(dá)式中的面積分。

組裝所有單元的重構(gòu)線性方程組式(33)，得到以全場單元的重構(gòu)多項(xiàng)式待定系數(shù)為未知數(shù)的大型稀疏線性方程組：

Au=b

(35)

式中：

(36)

可以證明矩陣A具有如下性質(zhì)：①A是對稱的；②A是正定的；③ 2D-A是對稱正定的。根據(jù)性質(zhì) ① 和性質(zhì) ②，矩陣A是對稱正定的，那么矩陣A一定是可逆的，變分重構(gòu)的線性方程組式(35) 有唯一解。重構(gòu)矩陣非奇異是變分重構(gòu)相對于CLS和k-exact重構(gòu)的一個(gè)巨大的優(yōu)勢，因?yàn)楹髢烧叨紵o法保證重構(gòu)矩陣的非奇異性。

線性方程組式(35)是全場單元的重構(gòu)線性方程組聯(lián)立組成的一個(gè)大型稀疏線性方程組，每個(gè)單元的實(shí)際依賴模板都是全場單元的集合。如果用直接法求解，那么無論如何都不能做到緊致。事實(shí)上，引入式(35)的唯一目的是證明重構(gòu)矩陣的正定性。在實(shí)際求解過程中，無需組裝式(35)所示的大型稀疏矩陣線性方程組。根據(jù)操作緊致性的要求，仍然采用迭代方法逐個(gè)單元求解式(33)。由于矩陣A和2D-A均是對稱正定的，因此根據(jù)數(shù)值分析理論[45]，變分重構(gòu)線性方程組式(33)的Block Jacobi、Block Gauss-Seidel和Block SOR(Successive Over-Relaxation)迭代都是收斂的。變分重構(gòu)矩陣的優(yōu)良特性不僅保證了重構(gòu)矩陣的非奇異性，還保證了迭代方法的收斂性。在實(shí)際計(jì)算中，一般采用Block SOR方法求解，其迭代公式為

i=1,2,…,N

(37)

式中：ω為松弛因子。進(jìn)一步研究發(fā)現(xiàn)，松弛因子ω取值在1.3左右的Block SOR方法譜半徑最小，收斂速度最快。因此求解變分重構(gòu)的迭代方法是Block SOR方法，松弛因子ω=1.3。

變分重構(gòu)是一類模板緊致高精度重構(gòu)方法的總稱，式(29)只給出了一個(gè)具體的例子。一般來講，構(gòu)造變分重構(gòu)的泛函要遵循以下原則：

1) 泛函對應(yīng)的變分重構(gòu)必須具有k-exact特性，這樣便能夠保證變分重構(gòu)的精度。

2) 泛函對應(yīng)的變分重構(gòu)必須是緊致的，即中心單元在泛函中只與面相鄰單元相關(guān)。

3) 泛函對應(yīng)的變分重構(gòu)矩陣必須是對稱正定的。

滿足條件1)～3)的重構(gòu)泛函可以有無窮多種。給出一種更為通用的泛函形式，即在式(30)中，If的定義為

(38)

式中：ωv、ωs、ωv,p和ωs,p為可人為選定的權(quán)函數(shù)；ΩL和ΩR分別為左、右兩側(cè)控制體。

式(38)不僅包含IJI，還包括界面兩側(cè)控制體上的重構(gòu)多項(xiàng)式及其各階導(dǎo)數(shù)相互之間的差在控制體上的積分。通過合理確定權(quán)函數(shù)可以進(jìn)一步提高變分重構(gòu)的精度和效率，這是目前正在開展的研究。

2.3 多步重構(gòu)

CLS重構(gòu)和變分重構(gòu)都是隱式重構(gòu)，即使采用迭代方法求解計(jì)算量仍然很大。在第3節(jié)將介紹與隱式時(shí)間積分結(jié)合時(shí)，進(jìn)一步提高隱式重構(gòu)計(jì)算效率的措施。當(dāng)與顯式時(shí)間積分結(jié)合時(shí)，隱式重構(gòu)不適于計(jì)算非定常流動。為改善緊致重構(gòu)在顯式時(shí)間積分情況下的計(jì)算效率，提出了一種多步重構(gòu)方法。這種方法大量使用矩陣操作，不容易直觀理解。簡單起見只介紹一維重構(gòu)，高維重構(gòu)算法詳見文獻(xiàn)[28]。

以三次多項(xiàng)式重構(gòu)為例說明多步重構(gòu)的思路和實(shí)施過程。一維情況下的緊致模板為

(39)

式中：φl,i為零均值基，φ1,i=(x-xi)/h，φ2,i=(x-xi)2/h2-1/12，φ3,i=(x-xi)3/h3，其中h=Δx,Δx為網(wǎng)格間距。

均勻網(wǎng)格上單元Ωi的三次重構(gòu)多項(xiàng)式表達(dá)式為

(40)

顯然，重構(gòu)多項(xiàng)式中有3個(gè)待定參數(shù)，而把式(18)的重構(gòu)關(guān)系應(yīng)用于一維情形，只能得到兩個(gè)獨(dú)立的重構(gòu)關(guān)系：

(41)

因此，重構(gòu)問題是欠定的。在緊致模板下實(shí)現(xiàn)高階顯式重構(gòu)可通過多步方法實(shí)現(xiàn)。一般三次多項(xiàng)式重構(gòu)需要分3步實(shí)現(xiàn)，逐一展開介紹。

1) 第1步重構(gòu)

(42)

(43)

(44)

(45)

2) 第2步重構(gòu)

經(jīng)過第1步重構(gòu)得到了每個(gè)單元對應(yīng)的規(guī)范化重構(gòu)關(guān)系，即式(46)。由于重構(gòu)模板的緊致性要求，在第2步重構(gòu)中只列出單元Ωi的緊致模板上的規(guī)范化重構(gòu)關(guān)系，即

(46)

(47)

進(jìn)一步，假定各單元的重構(gòu)函數(shù)之間是光滑過渡的，因此根據(jù)式(47)有

(48)

以及

(49)

通過式(48)和式(49)將緊致模板中Ωi-1、Ωi+1單元上重構(gòu)多項(xiàng)式的待定系數(shù)表示為Ωi上重構(gòu)多項(xiàng)式待定系數(shù)的線性組合。這個(gè)過程稱為延拓算法。利用這個(gè)算法，式(46)可以改寫為

(50)

經(jīng)過延拓運(yùn)算得到的式(50)稱為第2步重構(gòu)的原始重構(gòu)關(guān)系，式(50)的3個(gè)方程均與Ωi上的重構(gòu)多項(xiàng)式ui(x)相關(guān)。對于一維問題，式(50)中有3個(gè)方程，原則上可以確定三次重構(gòu)多項(xiàng)式的3個(gè)待定系數(shù)。但是對于二維和三維非結(jié)構(gòu)網(wǎng)格而言，第2步重構(gòu)的原始重構(gòu)關(guān)系可能不足以確定三次重構(gòu)多項(xiàng)式的待定系數(shù)。一般而言，無論是幾維問題，多步重構(gòu)的每一步比上一步提高一階精度是可行的。在第1步重構(gòu)中介紹過，第1步的重構(gòu)關(guān)系足以進(jìn)行線性重構(gòu)，則第2步總是可以進(jìn)行二次多項(xiàng)式重構(gòu)的。為此，將式(50)改寫為

(51)

(52)

(53)

根據(jù)式(52)，式(53)可以進(jìn)一步記為

(54)

式(54)是第2步重構(gòu)得到的規(guī)范化重構(gòu)關(guān)系，有兩個(gè)獨(dú)立方程。

3) 第3步重構(gòu)

首先列出Ωi的緊致模板上第2步重構(gòu)得到的規(guī)范化重構(gòu)關(guān)系，共有6個(gè)獨(dú)立方程。然后采用延拓算法將緊致模板中Ωi-1、Ωi+1單元上重構(gòu)多項(xiàng)式的待定系數(shù)表示為Ωi上重構(gòu)多項(xiàng)式待定系數(shù)的線性組合。這樣就得到了計(jì)算三次重構(gòu)多項(xiàng)式ui(x)的6個(gè)新的原始重構(gòu)關(guān)系。這6個(gè)關(guān)系足以用最小二乘法確定ui(x)中的3個(gè)待定系數(shù)，從而可以得到三次重構(gòu)多項(xiàng)式。其具體操作與第2步類似，不再贅述。

綜上所述，多步重構(gòu)主要算法為部分求逆算法和延拓算法。部分求逆算法在保留高階待定系數(shù)的情況下對低階系數(shù)做最小二乘，可以得到新的重構(gòu)關(guān)系即規(guī)范化重構(gòu)關(guān)系。一般第1步重構(gòu)對重構(gòu)多項(xiàng)式一次項(xiàng)做部分求逆，第2步對一次和二次項(xiàng)做部分求逆，并以此類推到更多重構(gòu)步。延拓算法延拓緊致模板上的規(guī)范化重構(gòu)關(guān)系，把緊致模板中Ωi-1、Ωi+1單元上重構(gòu)多項(xiàng)式的待定系數(shù)表示為Ωi上重構(gòu)多項(xiàng)式待定系數(shù)的線性組合，從而將緊致模板上的規(guī)范化重構(gòu)關(guān)系轉(zhuǎn)化為可以計(jì)算ui(x)待定系數(shù)的新原始重構(gòu)關(guān)系。多步重構(gòu)實(shí)際上是在每步隱含地?cái)U(kuò)充模板的過程。例如，在第2步重構(gòu)使用的規(guī)范化重構(gòu)關(guān)系(即式(45))涉及Si-1、Si和Si+13個(gè)緊致模板，包含的單元為{Ωi-2,Ωi-1,Ωi,Ωi+1,Ωi+2}。然而，在每一重構(gòu)步的操作遍歷所有單元時(shí)，僅對每個(gè)單元緊致模板的信息進(jìn)行操作，從而保證了算法的操作緊致性。

顯然，3個(gè)重構(gòu)步中部分求逆和延拓的操作可以在任意網(wǎng)格上對任意維數(shù)的問題進(jìn)行。因此，多步重構(gòu)的思路可以直接用于多維任意非結(jié)構(gòu)網(wǎng)格。對于多維問題，其計(jì)算量比k-exact重構(gòu)略大。如多維非結(jié)構(gòu)網(wǎng)格的三次多項(xiàng)式重構(gòu)需要分3步進(jìn)行，其中第3步的計(jì)算量和k-exact重構(gòu)類似。但是由于多步重構(gòu)前兩步的計(jì)算量較小，所以總的計(jì)算量比k-exact重構(gòu)增加并不明顯。與k-exact重構(gòu)相比，多步重構(gòu)的優(yōu)點(diǎn)是基于緊致模板；與CLS重構(gòu)和變分重構(gòu)相比，多步重構(gòu)的優(yōu)勢是可以靈活地采用顯式和隱式時(shí)間積分。

3 緊致高精度有限體積方法其他問題的討論

3.1 重構(gòu)的實(shí)施

第2節(jié)所述隱式和顯式重構(gòu)過程的具體實(shí)施在于確定每個(gè)單元上的重構(gòu)關(guān)系。如在開展CLS重構(gòu)或變分重構(gòu)且用迭代法求解時(shí)，需計(jì)算式(28)或式(37)中涉及的所有矩陣；采用顯式多步重構(gòu)時(shí)，也要進(jìn)行類似的矩陣運(yùn)算，但每步可以直接計(jì)算，無需迭代求解。無論采用哪種緊致重構(gòu)，都無需組裝全場的大矩陣，只需針對每個(gè)單元的重構(gòu)關(guān)系進(jìn)行操作。每個(gè)單元的重構(gòu)關(guān)系涉及矩陣規(guī)模都比較小，對守恒變量的每個(gè)分量都是相同的，因此總體計(jì)算量不大。此外，由式(25)和式(34)可以看出，涉及的矩陣均只與網(wǎng)格的幾何參數(shù)相關(guān)，在不進(jìn)行網(wǎng)格自適應(yīng)時(shí)，可以只在計(jì)算啟動后計(jì)算一次并存儲下來，從而顯著提高計(jì)算效率。在計(jì)算矩陣中的元素時(shí)，一般需要進(jìn)行數(shù)值積分運(yùn)算，如式(25)中的

(55)

或式(34)中的

(56)

在程序設(shè)計(jì)方面，必須記錄非結(jié)構(gòu)網(wǎng)格的基本數(shù)據(jù)結(jié)構(gòu)，如控制體頂點(diǎn)編號及坐標(biāo)、每個(gè)控制體的編號及由哪些頂點(diǎn)組成等。此外，由于緊致重構(gòu)的操作緊致性，在程序設(shè)計(jì)中只需記錄和當(dāng)前控制體有公共面的控制體；為方便通量計(jì)算，還需記錄所有界面編號及其兩側(cè)的控制體編號。可見，在數(shù)據(jù)結(jié)構(gòu)方面，緊致有限體積方法和常規(guī)的二階精度有限體積方法完全相同。這就使可在現(xiàn)有二階格心型有限體積程序基礎(chǔ)上，于不改變數(shù)據(jù)結(jié)構(gòu)的情況下引入高階緊致重構(gòu)。這給程序的發(fā)展帶來了很大方便。

3.2 隱式重構(gòu)提高計(jì)算效率的措施

即使采用迭代法求解，如果每次迭代到收斂，CLS重構(gòu)和變分重構(gòu)的計(jì)算量仍然是非常大的。因此，如何進(jìn)一步提高計(jì)算效率是基于隱式重構(gòu)的有限體積方法能夠成功應(yīng)用的關(guān)鍵。對于定常問題的計(jì)算，由于不關(guān)心時(shí)間推進(jìn)方法收斂過程中的時(shí)間精度，可在每個(gè)時(shí)間步只執(zhí)行迭代過程式(28)或式(37)一次，并讓迭代和時(shí)間推進(jìn)同時(shí)收斂到定常解。由于執(zhí)行迭代過程式(28)或式(37) 一次的計(jì)算量和進(jìn)行一次k-exact重構(gòu)基本相同，因此在定常計(jì)算中，隱式重構(gòu)的緊致有限體積方法計(jì)算量和基于傳統(tǒng)k-exact重構(gòu)的有限體積方法基本相同。圖3為定常問題基于CLS重構(gòu)和k-exact重構(gòu)的有限體積方法收斂速度比較，CLSFV代表基于CLS重構(gòu)的FV方法，可見二者基本相同[44]。

圖3 定常問題基于CLS重構(gòu)和k-exact重構(gòu)的FV方法收斂速度比較[44]Fig.3 Comparison of convergence speeds of CLS reconstruction and k-exact reconstruction FV schemes in steady flow[44]

在求解非定常問題時(shí)，為提高計(jì)算效率，隱式重構(gòu)一般和雙時(shí)間步(Dual Time-stepping)隱式格式聯(lián)合使用。由于流體力學(xué)基本方程常常是高度非線性的，在非定常時(shí)間推進(jìn)過程中必須引入迭代求解過程，即使使用常規(guī)重構(gòu)方法時(shí)也是如此。雙時(shí)間步隱式方法就是常見的一種迭代過程,此方法把一個(gè)物理時(shí)間步的推進(jìn)過程分為虛擬時(shí)間步的多次內(nèi)迭代，使迭代解逐步趨近非定常真解。利用這個(gè)性質(zhì)提出了一種重構(gòu)和雙時(shí)間步耦合迭代解法[25]，在每次內(nèi)迭代過程中，迭代過程式(28) 或式(37)也只執(zhí)行一次。隨著內(nèi)迭代的進(jìn)行，重構(gòu)迭代和時(shí)間推進(jìn)迭代同步收斂。顯然，這種方法的計(jì)算量也和基于傳統(tǒng)k-exact重構(gòu)的雙時(shí)間隱式有限體積方法相當(dāng)。通過上述分析可以確信，隱式重構(gòu)有限體積方法在一定條件下仍然可以達(dá)到很高的計(jì)算效率。事實(shí)上，緊致重構(gòu)求解的是一個(gè)線性方程組，而虛擬時(shí)間步內(nèi)迭代求解的是一個(gè)非線性方程組。數(shù)值實(shí)驗(yàn)表明，緊致重構(gòu)的收斂速度會比虛擬時(shí)間步內(nèi)迭代的收斂速度快得多。圖4為靜態(tài)CLS重構(gòu)收斂速度和時(shí)間推進(jìn)殘差收斂速度的比較，可以清楚地看到這一現(xiàn)象[44]。

圖4 靜態(tài)CLS重構(gòu)和雙時(shí)間步迭代殘差收斂速度比較[44]Fig.4 Comparison of residual convergence speeds of static CLS reconstruction and dual-time stepping iteration[44]

3.3 邊界處理

邊界處理是影響數(shù)值計(jì)算精度的重要因素，主要討論重構(gòu)過程中的邊界處理。物理邊界條件的實(shí)施精度依賴于邊界附近重構(gòu)的精度，一般直接作用于數(shù)值通量，和常規(guī)有限體積方法相同。3種重構(gòu)方法的邊界處理有一些不同的特點(diǎn)。在邊界附近的典型緊致重構(gòu)模板如圖5所示。

對于CLS重構(gòu)而言，物理邊界(非周期邊界)單元要進(jìn)行特殊處理。這是因?yàn)槲锢磉吔鐔卧拿嫦噜弳卧獋€(gè)數(shù)至少比內(nèi)點(diǎn)單元少一個(gè)，造成重構(gòu)線性方程個(gè)數(shù)的減少。如對于圖5所示的邊界單元，有2個(gè)面相鄰單元，如果仿照內(nèi)點(diǎn)單元取M=k-1，對于CLS重構(gòu)，式(27)中與Ωj3有關(guān)的項(xiàng)不再存在。因此，二維三次多項(xiàng)式重構(gòu)時(shí)，有9個(gè) 待定系數(shù)，重構(gòu)關(guān)系中有12個(gè)方程。雖然方程個(gè)數(shù)仍然大于未知量個(gè)數(shù)，但數(shù)值測試表明這種邊界重構(gòu)格式不穩(wěn)定。即使將守恒導(dǎo)數(shù)的階次調(diào)到最大，即M=k，獲得了20個(gè)方程，也很難獲得穩(wěn)定的邊界重構(gòu)格式。但數(shù)值測試結(jié)果表明，將邊界單元的重構(gòu)格式精度降低一階便可以得到穩(wěn)定的邊界重構(gòu)格式。根據(jù)Gustafsson[46]的分析，邊界單元重構(gòu)格式精度降低一階不會造成整體格式精度階數(shù)的降低。測試結(jié)果表明，邊界單元重構(gòu)精度降低一階時(shí)，用L1誤差范數(shù)衡量，整體格式勉強(qiáng)達(dá)到理論精度階數(shù)；用L∞誤差范數(shù)衡量，整體格式精度會降低一階。

圖5 邊界單元的緊致重構(gòu)模板Fig.5 Compact reconstruction stencil of boundary cells

多步重構(gòu)和CLS重構(gòu)面臨的問題基本相同。特別是當(dāng)邊界單元只有一個(gè)面相鄰內(nèi)部單元的情況，有可能出現(xiàn)多步重構(gòu)的第1或第2步原始重構(gòu)關(guān)系數(shù)量不夠的情況。為解決這個(gè)問題，提出一種邊界單元的異步重構(gòu)方案，即第m步重構(gòu)先對內(nèi)部單元進(jìn)行，在得到內(nèi)部單元的規(guī)范化重構(gòu)關(guān)系后，再把其規(guī)范重構(gòu)關(guān)系延拓到面相鄰邊界單元；最后利用這些規(guī)范化重構(gòu)關(guān)系進(jìn)行邊界單元的第m步重構(gòu)。對于多步重構(gòu)，邊界單元重構(gòu)階數(shù)也要比內(nèi)點(diǎn)低一階。

在處理黏性流動的無滑移邊界條件時(shí)，CLS重構(gòu)和多步重構(gòu)可把無滑移條件作為一個(gè)重構(gòu)關(guān)系，與其他重構(gòu)關(guān)系同時(shí)進(jìn)行最小二乘計(jì)算；也可把無滑移條件當(dāng)作約束條件，采用約束最小二乘方法求解。

在邊界處理方面，變分重構(gòu)有明顯優(yōu)勢。其最突出的優(yōu)勢是可以達(dá)到與內(nèi)點(diǎn)一致的精度，無需降階。另一個(gè)優(yōu)勢是可以通過引入邊界泛函考慮物理邊界條件的影響。以式(30)的泛函為例，當(dāng)考慮邊界條件時(shí)，式(30)改寫為

(57)

式中：Ibf為表征邊界條件的IJI；Nbf為物理邊界面的總數(shù)。以無滑移等溫固壁邊界條件和對稱邊界條件為例介紹邊界面的IJI構(gòu)造方法。對于無滑移等溫固壁，守恒變量U的邊界IJI為

(58)

式中：dbf為邊界面bf所屬的單元中心到邊界面bf的中心的距離；下標(biāo)L和bf分別代表邊界界面內(nèi)側(cè)和外側(cè)的物理量，外側(cè)的物理量由邊界條件確定。對于壁面溫度為Twall的靜止固壁，守恒變量在邊界面上的值計(jì)算公式為

(59)

式中：ρ、E、cp、u、γ和v分別為密度、總能、定壓比熱、比熱比和兩個(gè)方向上的速度分量。

3.4 激波捕捉

良好的激波捕捉性能是對可壓縮流動求解器的基本要求，對于高階格式一般需要引入限制器或非線性人工黏性抑制激波附近的振蕩。對有限體積和DG等方法的限制器開展了比較系統(tǒng)的研究[47-50]。對于緊致高精度有限體積方法，推薦的限制器是WBAP(Weighted Biased Averaging Procedure)限制器[48]。WBAP限制器具有如下7個(gè)優(yōu)點(diǎn)：① 光滑區(qū)保精度；② 基于緊致模板；③ 直 接作用于重構(gòu)多項(xiàng)式系數(shù)，無需數(shù)值積分，計(jì)算效率較高；④ 通過采用二次重構(gòu)技術(shù)，無需在多個(gè)模板上開展重構(gòu)，計(jì)算效率顯著優(yōu)于WENO格式；⑤ 適用于結(jié)構(gòu)和非結(jié)構(gòu)網(wǎng)格；⑥ 通過采用逐次限制的方法可對任意階精度重構(gòu)進(jìn)行限制；⑦ 在一定條件下可退化為TVD限制器或滿足極大值原理，具有良好的抑制振蕩的能力。其具體原理和實(shí)施過程見文獻(xiàn)[48]，不再贅述。WBAP限制器和這3種緊致重構(gòu)算法都可以較好地結(jié)合，得到的有限體積格式具有很好的激波捕捉能力。圖6為三次CLS重構(gòu)有限體積格式(四階精度)計(jì)算復(fù)雜激波繞射的結(jié)果。

圖6 正激波繞過復(fù)雜幾何形狀障礙物(t=1.5時(shí)刻障礙物附近的密度分布)Fig.6 Normal shock wave flows past obstacle with complex geometry (density distribution near obstacle at t=1.5)

3.5 高階網(wǎng)格

復(fù)雜幾何外形(如飛機(jī)機(jī)身)通常包含曲面邊界，如用直邊網(wǎng)格表示曲面邊界，會過濾掉高階曲面信息，從而降低高階數(shù)值方法的精度。對于DG等方法，可能計(jì)算出非物理解。有限體積方法雖然一般不會出現(xiàn)非物理解，但整體精度會下降。處理曲面邊界的實(shí)用方法是采用高階曲面或曲線網(wǎng)格，目前已有一些高階網(wǎng)格的生成方法[51]。當(dāng)采用高階網(wǎng)格時(shí)，必須采用參數(shù)變換的方法計(jì)算式(55)、式(56)的數(shù)值積分。發(fā)展了直角參數(shù)坐標(biāo)系下建立插值基函數(shù)的一般方法，并給出高階三角形、四邊形、四面體、三棱柱和六面體單元的參數(shù)變換公式。借助于參數(shù)變換公式可求出數(shù)值積分點(diǎn)的物理坐標(biāo)和Jacobi值進(jìn)行數(shù)值積分。計(jì)算金字塔單元的體積分時(shí)，最佳方案是將金字塔單元剖成兩個(gè)四面體，金字塔單元的體積分就等于這兩個(gè)四面體單元的體積分之和。表2給出了參數(shù)變換對應(yīng)的網(wǎng)格類型及總控制點(diǎn)數(shù)。圖7為二維三角形二次網(wǎng)格及其6個(gè)控制點(diǎn)(黑圓點(diǎn))，圖8為三維四面體二次網(wǎng)格及其10個(gè)控制點(diǎn)。對于坐標(biāo)x，參數(shù)變換公式可以一般地寫為

圖7 二次三角形單元TRI_6Fig.7 Quadratic triangular cell TRI_6

圖8 二次四面體單元TETRA_10Fig.8 Quadratic tetrahedron cell TETRA_10

表2 網(wǎng)格單元類型及其控制點(diǎn)數(shù)Table 2 Types of grid units and number of reference points

(60)

式中：ξ、η和ζ為參數(shù)坐標(biāo)，取值范圍均為[0,1]；xi為控制點(diǎn)笛卡爾坐標(biāo)；vi(ξ,η,ζ)為第i個(gè)控制點(diǎn)的Lagrange插值基函數(shù)；nc為控制點(diǎn)的總數(shù)。

舉例說明基函數(shù)vi(ξ,η,ζ)的表達(dá)式。對于二維三角形二次網(wǎng)格，有v1=(ξ+η-1)(2ξ+2η-1)，v2=ξ(2ξ-1)，v3=η(2η-1)，v4=-4ξ(ξ+η-1)，v5=4ξη，v6=-4η(ξ+η-1)。

對于三維四面體二次網(wǎng)格，有v1=(-1+ζ+η+ξ)(-1+2ζ+2η+2ξ)，v2=ξ(-1+2ξ)，v3=η(-1+2η)，v4=ζ(-1+2ζ)，v5=-4ξ(-1+ζ+η+ξ)，v6=4ηξ，v7=-4η(-1+ζ+η+ξ)，v8=-4ζ(-1+ζ+η+ξ)，v9=4ζξ，v10=4ζη。

其他情形見文獻(xiàn)[44]。

3.6 時(shí)間推進(jìn)

在求解定常問題時(shí)，緊致重構(gòu)可以和任意時(shí)間推進(jìn)格式相結(jié)合。當(dāng)求解非定常問題時(shí)，顯式多步重構(gòu)可以和任意顯式和隱式時(shí)間推進(jìn)格式相結(jié)合。隱式重構(gòu)即CLS重構(gòu)和變分重構(gòu)，只有和隱式雙時(shí)間步格式相結(jié)合使用才能實(shí)現(xiàn)高求解效率。此外，在求解非定常問題時(shí)，高精度有限體積方法在時(shí)間方向也要達(dá)到高階精度。例如，當(dāng)采用三次多項(xiàng)式緊致重構(gòu)時(shí)，空間為四階精度，則時(shí)間離散格式也應(yīng)達(dá)到四階精度。三次CLS和變分重構(gòu)下，采用的時(shí)間離散格式為多步隱式Runge-Kutta方法，具體為三步四階精度的SDIRK4方法[52]。其具體實(shí)施步驟如下。

(61)

三步四階精度SDIRK4方法的離散公式為

(62)

式中：

(63)

其中：bα為向量；aαβ為矩陣。矩陣aαβ和向量bα各元素的值見表3，ζ=0.128 886 400 515。

表3 三步四階SDIRK4方法的Butcher表Table 3 Butcher table of 3-stage p4 SDIRK4 method

求解式(63)時(shí)，在其左端加一個(gè)虛擬時(shí)間項(xiàng)，采用雙時(shí)間步方法迭代求解。迭代方法為LU-SGS(Lower-Upper Symmetric Gauss-Seidel)或GMRES(Generalized Minimum RESidual)+LU-SGS方法。

3.7 計(jì)算精度和效率

按照式(1)的計(jì)算效率定義，比較合理的計(jì)算精度和效率的衡量標(biāo)準(zhǔn)是在數(shù)值誤差一定時(shí)，計(jì)算時(shí)間越短越好；或計(jì)算時(shí)間一定時(shí)，數(shù)值誤差越小越好。按照這一標(biāo)準(zhǔn)，無論哪一種重構(gòu)方法，高階格式的性能都顯著優(yōu)于二階格式。對于同樣精度格式之間的比較，采用二維等熵渦問題進(jìn)行分析。等熵渦問題的計(jì)算網(wǎng)格見圖9。將網(wǎng)格逐次加密進(jìn)行求解。圖10(a)[25]表明，同階的CLS有限體積格式明顯優(yōu)于k-exact格式；圖10(b)[26]表明，在計(jì)算網(wǎng)格較密時(shí)，變分有限體積(Variational Finite Volume，VFV)格式計(jì)算效率優(yōu)于CLS有限體積格式，網(wǎng)格較少時(shí)則相反；圖10(c)[27]表明，多步重構(gòu)(MSR)有限體積格式和k-exact格式基本相當(dāng)，其主要優(yōu)勢在于模板的緊致性。

圖9 二維等熵渦問題的計(jì)算網(wǎng)格Fig.9 Computing grid of isentropic vortex problem in 2D

圖10 等熵渦問題計(jì)算效率比較Fig.10 Comparison of computational efficiency of isentropic vortex problem

對于有激波的三維問題計(jì)算，高階緊致有限體積也取得了很好的結(jié)果。圖11[44]為四階變分有限體積格式計(jì)算出ONERA-M6跨聲速繞流問題的馬赫數(shù)等值線，計(jì)算網(wǎng)格由約100萬四面體網(wǎng)格組成。將變分有限體積方法應(yīng)用于求解RANS方程，圖12[43]為四階變分有限體積格式采用SA模式計(jì)算平板湍流邊界層的速度剖面，與其進(jìn)行比較的是網(wǎng)格數(shù)為高階格式網(wǎng)格16倍時(shí)用CFL3D軟件得到的結(jié)果,可見二者精度基本相當(dāng)。

圖11 四階VFV格式計(jì)算出的ONERA-M6跨聲速繞流問題馬赫數(shù)等值線[44]Fig.11 Contours of Mach number in transonic flow past ONERA-M6 airfoil computed by fourth order VFV scheme[44]

圖12 四階VFV格式得到的湍流平板邊界層問題速度剖面[43]Fig.12 Velocity profile in turbulence boundary layer problem computed by fourth order VFV scheme[43]

3.8 發(fā)展及推廣

由于變分重構(gòu)非奇異的優(yōu)良特性，目前主要發(fā)展的是基于變分重構(gòu)的非結(jié)構(gòu)網(wǎng)格高精度緊致有限體積方法。

潘建華等[53]將自適應(yīng)網(wǎng)格技術(shù)引入變分重構(gòu)有限體積方法，顯著提高了計(jì)算效率。最近，進(jìn)一步研究了非定常流動網(wǎng)格自適應(yīng)準(zhǔn)則問題，提出了基于亞網(wǎng)格能量耗散的自適應(yīng)算法。圖13(a)和圖13(b)是采用這種準(zhǔn)則計(jì)算三維圓柱湍流的網(wǎng)格及渦量等值線(只顯示了一個(gè)橫截面)，初步結(jié)果表明，這種方法比基于渦量的網(wǎng)格自適應(yīng)準(zhǔn)則效果更好。

圖13 四階VFV格式得到的圓柱繞流問題自適應(yīng)網(wǎng)格分布和渦量等值線Fig.13 Distribution of adaptive grids and contours of vorticity in flow past cylinder computed by fourth order VFV scheme

黃乾旻等[54]將變分重構(gòu)有限體積方法推廣到求解非守恒方程，通過引入對流重構(gòu)方法克服了常規(guī)方法的數(shù)值穩(wěn)定性問題；他還研究了大長寬比網(wǎng)格上變分重構(gòu)的構(gòu)造以及RANS模式的保正問題，顯著提高了高精度有限體積方法RANS模擬的精度和魯棒性。圖14為采用SA湍流模式計(jì)算30P30N高升力構(gòu)型的湍流流動的計(jì)算結(jié)果。結(jié)果表明，通過改進(jìn)變分重構(gòu)泛函和湍流模型的保正處理，變分重構(gòu)高精度有限體積法可實(shí)現(xiàn)RANS方程魯棒穩(wěn)定的高精度求解，不會出現(xiàn)負(fù)的湍流黏性系數(shù)，且壁面摩阻系數(shù)比二階格式更為準(zhǔn)確。

圖14 四階VFV格式得到的16°攻角下30P30N高升力構(gòu)型參數(shù)Fig.14 Parameter of 30P30N high lift configuration at angle of attack of 16° computed by fourth order VFV scheme

變分重構(gòu)算法也促進(jìn)了其他非結(jié)構(gòu)網(wǎng)格高精度格式的發(fā)展。如把變分重構(gòu)作為重構(gòu)DG方法的重構(gòu)算法取得了很好的結(jié)果，有效解決了多維重構(gòu)的穩(wěn)定性問題，代表性工作見文獻(xiàn)[55-56]。又如Tamaki和Imamura[57]在CLS重構(gòu)的啟發(fā)下，提出了迭代最小二乘重構(gòu)方法。Nishikawa[58]在變分重構(gòu)的啟發(fā)下提出了隱式Green-Gauss梯度重構(gòu)方法。變分重構(gòu)方法還將作為非結(jié)構(gòu)網(wǎng)格高精度有限體積方法的方案之一，在國家數(shù)值風(fēng)洞工程中繼續(xù)得到發(fā)展及應(yīng)用。

4 結(jié) 論

緊致重構(gòu)技術(shù)是為解決非結(jié)構(gòu)網(wǎng)格高精度有限體積方法重構(gòu)模板巨大這一瓶頸問題而提出的。在概念層面，提出了操作緊致性這一重要概念，即不要求實(shí)際上采用的模板是緊致的，只要求在數(shù)值求解過程中所有的數(shù)據(jù)結(jié)構(gòu)和計(jì)算操作都只用到緊致模板中的信息。這一概念創(chuàng)新極大擴(kuò)展了構(gòu)造重構(gòu)算法的多樣性，使發(fā)展高階緊致重構(gòu)成為可能。

在方法層面，提出了CLS、變分、多步等緊致重構(gòu)方法。從實(shí)際性能角度看，變分重構(gòu)非奇異、計(jì)算效率高、邊界不降階，總體上是3種方法中最優(yōu)的。CLS重構(gòu)可以看作k-exact重構(gòu)的推廣，比較簡單，性能上不如變分重構(gòu)，但優(yōu)于多步重構(gòu)。變分和CLS重構(gòu)是隱式重構(gòu)，在計(jì)算非定常流動時(shí)必須和隱式雙時(shí)間步方法結(jié)合才能達(dá)到較高計(jì)算效率；而多步重構(gòu)無須迭代求解，可以和任意顯式和隱式時(shí)間推進(jìn)格式結(jié)合，有效計(jì)算非定常流動。因此，幾種緊致重構(gòu)方法各有特點(diǎn)，共同構(gòu)成了非結(jié)構(gòu)網(wǎng)格高精度有限體積方法的基礎(chǔ)算法。

在實(shí)施層面，綜述了在具體實(shí)施方案、邊界處理、激波捕捉、高階網(wǎng)格、時(shí)間推進(jìn)等方面一些需要注意的問題。

在應(yīng)用層面，所發(fā)展的方法已經(jīng)成功地應(yīng)用于一系列二維、三維，有黏、無黏，層流和湍流流動的模擬，在計(jì)算精度和計(jì)算效率等方面相對于二階精度有限體積方法顯示出了明顯的優(yōu)勢。

緊致重構(gòu)有限體積方法徹底解決了常規(guī)有限體積方法的大模板問題，為非結(jié)構(gòu)網(wǎng)格高精度有限體積方法的發(fā)展提供了廣闊空間。近期的一些新進(jìn)展表明，緊致高精度有限體積方法是一種有發(fā)展前途的新方法。