摘 要：文章中提出了降階的Spectral Volume 方法 （SV方法） 的概念，對方法的公式進(jìn)行推導，并成功將其應用與一維線(xiàn)性傳遞方程。對降階的SV方法的穩定性進(jìn)行了Fourier分析，并與傳統的SV方法進(jìn)行了對比。實(shí)驗證明，降階的SV方法在初值連續情況下，在線(xiàn)性傳遞方程上可以達到所期望的數值精度，并且在數值色散和數值耗散特性上較傳統SV方法有顯著(zhù)提高。
關(guān)鍵詞：SV方法 一維守恒律 降階
中途分類(lèi)號：R445文獻標識碼：A文章編號：1674-098X（2014）10（b）-0211-04
1 研究背景
計算流體力學(xué)（CFD）從創(chuàng )立至今得到了飛速發(fā)展，在工業(yè)界起著(zhù)不可替代的作用。然而目前工業(yè)界所用到的CFD計算通常只能達到二階精度，遠遠達不到人們的需要，并且對網(wǎng)格以及內存和CPU的要求超出了當今計算機所能承受的范圍。另外，運算區域通常包含復雜單光滑的流體結構，比如渦旋與剪切流層的相互作用區，以及同時(shí)包含激波與其他復雜流動(dòng)的區域等。因此，高精度高分辨率格式的研究一直是CFD領(lǐng)域中非?；钴S的一個(gè)課題，發(fā)展高精度高分辨率的計算格式成為CFD工作者的一大研究方向。為采用高精度格式不但可以降低對網(wǎng)格規模的苛刻要求，而且能夠正確分辨其中復雜的流動(dòng)現象。對網(wǎng)格點(diǎn)數目的要求低，具有減少計算機內存需求和CPU時(shí)間消耗的優(yōu)勢。其中比較著(zhù)名的是Harten提出的高階精度的ENO（essentially non-oscillatory） 格式，隨后Liu，Osher和Chen發(fā)展了原有的ENO格式，提出了有顯著(zhù)改進(jìn)的WENO（weightedessentially non-oscillatory）格式，并獲得了迅速的發(fā)展，并得到了一系列的研究。ENO/WENO格式善于捕捉物理間斷，其應用迅速發(fā)展到多相流領(lǐng)域。
新型格式發(fā)展的目標之一是發(fā)展可以應用于非結構網(wǎng)格的可靠的高階格式。理論和實(shí)踐都表明，大多數高精度、高分辨率格式的性能都嚴重地依賴(lài)于網(wǎng)格的光滑性，而工業(yè)界所需要模擬的物體通常具有復雜的幾何外形（圖1），CFD應用最多的航空領(lǐng)域所需要的幾何外形更是相當復雜，而這給結構網(wǎng)格的生成帶來(lái)很大困難，一種直接有效的解決辦法是采用非結構網(wǎng)格。然而由于穩定性等因素的影響，對于三維有限體積法，非結構算法模板所需單元往往較多，而過(guò)多的模板單元會(huì )帶來(lái)邊界處理、內存占用、編程復雜等多方面的困難。因此實(shí)際應用中大都只使用二階精度的非結構算法，而這遠遠滿(mǎn)足不了一些工業(yè)領(lǐng)域對計算精度的要求?？梢詰门c非結構網(wǎng)格的高精度、高分辨率、低耗散的實(shí)用數值格式的發(fā)展成為當今CFD的重要任務(wù)之一。以此為目標的一些高階格式已經(jīng)發(fā)明出來(lái)，其中比較著(zhù)名的是Cockburn發(fā)展的間斷Galerkin有限元方法（DG），間斷 Galerkin有限元方法的插值模板小，各單元之間僅通過(guò)界面通量計算相聯(lián)系，因而比較適合在復雜外形。而另一種成功的格式是Wang在自2002年發(fā)展的應用于非結構網(wǎng)格的算法SV方法。SV方法是一種高精度高分辨率低耗散的數值格式。同間斷Galerkin有限元方法一樣，SV方法可以達到相同的精度。然而相比之下，SV方法比間斷Galerkin有限元方法有著(zhù)更好的穩定性（可以允許更大的CFL數）。與傳統的有限體積法相比，對計算機內存和CPU的要求更低。所以SV方法有著(zhù)非常大的能在將來(lái)被應用于工業(yè)界的潛力。同時(shí)，SV方法的穩定性和可行性還需要進(jìn)一步驗證和改進(jìn)。
2 SV方法介紹
SV方法，是比較新提出的另一種可以應用于完全非結構網(wǎng)格的高精度、高分辨率、低耗散的守恒型格式。Wang于2002年提出一維SV方法的概念并隨后將其發(fā)展到二維和三維。SV方法是經(jīng)典有限體積法的延伸。它提出了控制體積（Control Volume）的概念，從而避免的復雜模板（stencil）的使用，從而使格式可以輕松應用與非結構網(wǎng)格上。經(jīng)典的有限體積法中，變量在每個(gè)網(wǎng)格中的值被用來(lái)進(jìn)行多項式插值。為了達到相應的精度，由相鄰多個(gè)單元格組成的模板被用來(lái)進(jìn)行插值。而非結構網(wǎng)格使得模板中相鄰網(wǎng)格的選取成為一大難題。不同于經(jīng)典的有限體積法，SV方法把網(wǎng)格中的每一個(gè)單元格進(jìn)一步分成若干控制體積（CV），對每個(gè)控制體積中變量的平均值加以利用來(lái)實(shí)現多項式插值從而達到所要求的精度。SV方法的優(yōu)勢在于其在保證相同精度的同時(shí)免除了對周?chē)W(wǎng)格的依賴(lài)，從而降低了對網(wǎng)格的要求。此外，由于不同單元格之間的信息交流少，SV方法非常適宜并行運算的實(shí)現。
有限體積法（包括ENO和WENO方法）利用相鄰多個(gè)自由度（模板）內狀態(tài)變量的值進(jìn)行多項式構建。而在SV方法中，我們把這個(gè)單個(gè)的單元格叫做Spectral Volume（SV）。我們只利用此SV中的信息進(jìn)行高階多項式的構建，為了得到所需要的自由度，每個(gè)SV被劃分為更小的體積單元 Control Volume（CV）。而方法的精度取決于每個(gè)SV劃分為CV的數量。
對于二階的SV方法，流量的計算是線(xiàn)性的，所以并不復雜。對于高階的情況 （k>2），重構過(guò)程如下。在下一個(gè)時(shí)間步上，每個(gè)CV中狀態(tài)變量的平均值各自更新，在SV內部，狀態(tài)變量在床CV邊界時(shí)時(shí)連續的，所以我們并不需要Riemann求解器，流量的值我們可以通過(guò)插值多項式直接到的解析值。當兩個(gè)相鄰CV的邊界也是 SV邊界時(shí)，我們需要通過(guò)插值得到邊界上Gauss求積點(diǎn)處狀態(tài)變量的值，然后通過(guò) Riemann求解器來(lái)計算流量的值。這樣方法保證了計算的簡(jiǎn)潔以及高階特性。
總的來(lái)說(shuō)，假設我們重構出至多k-1 階多項式，k階SV方法計算的一般過(guò)程是通過(guò)插值計算每個(gè)Gauss積分點(diǎn)處狀態(tài)量的值；利用k階高斯求積形式和Riemann 求解器計算單元邊界上面通量的積分。在單元格內部CV邊界上，由于我們假定狀態(tài)變量在單元格內是連續的，通量的值我們直接應用重構多項式的解析值；利用TVD的Runge-Kutta格式進(jìn)行時(shí)間推進(jìn)。
如同有限體積方法一樣，在計算過(guò)程中SV方法利用狀態(tài)變量的平均值。然而SV方法并不像有限體積法那樣利用在整個(gè)單元格內平均，它利用的是狀態(tài)變量在每個(gè)控制體積內的平均值，通過(guò)CV的劃分來(lái)增加自由度，從而達到所需要的高階。關(guān)于SV方法的具體操作以及一維和多維的各種算例可以在文章中找到。
之前提到的間斷Galerkin有限元方法和SV方法在特性上有許多相似之處，同樣是緊致類(lèi)方法，從而適于并行運算，能達到高階精度，由于在單元格邊Riemann求解器的應用，它們是完全守恒型的數值計算格式，并且伴隨著(zhù)TVD或者TVB限制器的應用，兩種方法都非常適于應用于非結構網(wǎng)格和復雜的幾何外形。自然而然，兩種方法的研究比較出現在很多文章中。由于考慮局部的細致信息，SV方法被認為相對應間斷Galerkin有限元方法可以更加細致的捕捉間斷。此外間斷Galerkin有限元方法會(huì )隨著(zhù)計算階數的升高而導致相應的CFL數減小，穩定性會(huì )隨之下。而Zhang和Shu通過(guò)理論以及算例指出兩種方法可以達到相同的精度，SV方法可以允許更大的CFL數，然而在相同的網(wǎng)格上SV方法的誤差大于間斷Galerkin有限元方法。
3 降階SV方法
SV方法的穩定性受很多因素影響，我們提出的一種對其可能的改進(jìn)是降階的SV 方法。由于SV方法是緊致型格式，所有信息的處理幾乎都在單元格自身內部進(jìn)行，這使得其格式的改造變得相對容易。一種可能的構想是減小SV方法基底的維度。以基本的一維情況為例，理論上k階的SV方法需要將SV劃分成k個(gè)CV并以k-1階多項式進(jìn)行近似
（1）
在傳統的SV方法中，基底的階數與 CV的個(gè)數是相等的，從而保證方法的精度。在這里我們提出，可以嘗試以降低多項式基底的維度來(lái)改善SV方法在穩定性方面的表現。
假如SV被分為k個(gè)CV并以k-1-n階多項式（n≤k-1）近似，其理論精度值為 k-n。
（2）
如圖2給出了一維情況下以二階多項式描繪三個(gè)CV中狀態(tài)變量值的示意圖。它的實(shí)際精度以及在具體算例中的表現需要進(jìn)行驗證。
另外可以期待的一點(diǎn)是此方法對格式穩定性的影響，可以通過(guò)Fourier分析法研究格式的數值色散和數值耗散特性，k個(gè)CV并以k-1-n階多項式近似的降階的SV方法被期待比傳統的k-n階SV方法有更好的穩定性，即在數值色散和數值耗散方面有更好的表現。
SV方法插值邊界處Gauss積分點(diǎn)上狀態(tài)變量的值的關(guān)系式是：
（3）
其中向量指的是SV中每個(gè)CV上狀態(tài)變量的平均值，其維度是此單元格內CV的數量，向量q表示邊界上每個(gè)的Gauss求積點(diǎn)處狀態(tài)變量的值，其維度是此單元格內 Gauss 求積點(diǎn)的數量。
相應的構建方法是：
（4）
其中為的維度，為SV內部近似狀態(tài)變量的重構多項式，它是基底多項式的線(xiàn)性組合
（5）
其中在某一固定的時(shí)間步上，為常數，代表基底線(xiàn)性組合的系數。是多項式基底的維度。而q的每個(gè)元素則就是重構多項式在CV邊界處對應Gauss積分點(diǎn)上的取值。
（6）
指的是Gauss積分點(diǎn)的數目，其中指第k個(gè)Gauss積分點(diǎn)的坐標。從而矩陣運算可以從中提取出。我們得到多項式取值的矩陣P和基底在每個(gè)CV上的平均值矩陣L。
（7）
（8）
分析一下維度我們得到
（9）
對于不同的劃分方式，我們只需要根據基底和CV的幾何分布得到矩陣P和L便可以進(jìn)行運算。但這里我們假設，L是可逆方陣。如果我們縮減基底的維數，L便不再可逆，q的表達式會(huì )有相應變化，這也將是新的自適應階數SV方法的創(chuàng )新之處。
（10）
此時(shí)我們在計算公式中以L(fǎng)的假逆矩陣代替其逆矩陣。而我們定義其假逆矩陣：
（11）
這里M是一個(gè)正定矩陣，它的階數等于每個(gè)SV中CV的數量。在這種定義之下，我們重新定義了一種基于矩陣M的數量積： 對于任意向量a和b，我們有
（12）
與此同時(shí)，利用幾何中投影的定義，我們也定義了一種基于矩陣 M 的投影：
（13）
在具體操作中我們發(fā)現，M矩陣的選取要滿(mǎn)足許多要求。當我們利用降階的SV方法時(shí)其中非常重要的一點(diǎn)是如下守恒定律必須要滿(mǎn)足：
（14）
其中是保存狀態(tài)變量平均值的向量，其中的每個(gè)元素保存了同一個(gè)SV中指定CV上狀態(tài)變量的平均值，是所研究 SV中CV的數量。這個(gè)條件的物理意義是，對一個(gè)SV中每個(gè)CV上的狀態(tài)變量平均值投影到矩陣L的像空間，投影前后保證整個(gè)SV中狀態(tài)變量的積分總量不變。事實(shí)上，如果這一關(guān)系沒(méi)有滿(mǎn)足，則原方程的守恒形式將被打破。之后的數值測試結果表明，如果以上提到的守恒定律沒(méi)有被滿(mǎn)足，計算的結果是完全錯誤的。滿(mǎn)足這一條件的矩陣M也不是唯一的，利用簡(jiǎn)單的矩陣運算，我們得到一個(gè)符合守恒定律條件的如下矩陣M：
（15）
數值測試
在測試的初期，我們通過(guò)編寫(xiě)如下基本的一維線(xiàn)性輸運方程 Matlab 程序驗證 SV 方法精度：
（16）
方程的輸運速度設為a=1。
為了方便驗證且不失一般性，我們提出T=60s以保證。
我們首先驗證初始條件連續
的情況下SV方法的精度。我們采用初始條件：
（17）
測試結果顯示在表1中。其中自由度指的是區間[0，3]被劃分成的CV的數量。在命名方法上，我們采用：SV+每個(gè)SV中CV數+模擬多項式的階數。如SV43表示每個(gè)SV劃分為四個(gè)CV，但降一階，采用三階精度的多項式近似。若每個(gè)SV中CV數=模擬多項式的階數，則格式為傳統的SV方法。在測試中我們采用三階Runge-Kutta時(shí)間積分方法。
我們可以從表1中看出，理論上達到二階的數值格式，即SV22，SV32及SV42方法，在低自由度的網(wǎng)格上由于誤差太大，格式很難達到所期望的精度，而在過(guò)高自由度的網(wǎng)格上，由于數值計算誤差，所得到的精度可能與理論值有所偏差。但總體來(lái)說(shuō)，降階的SV方法在允許的誤差范圍內可以達到理論預期值。
4 格式穩定性研究
我們通過(guò)Fourier分析法對格式的穩定性進(jìn)行研究自適應SV方法較之于傳統SV方法的優(yōu)勢。這里我們之研究格式空間離散帶來(lái)的誤差而不考慮時(shí)間推進(jìn)格式。所以這里我們研究計算中的半離散的數值格式：
（18）
其中即我們所研究的狀態(tài)變量。通過(guò)對方程系統的線(xiàn)性化我們得到矩陣A。這里我們研究最基本的一屆線(xiàn)性傳遞方程，從而矩陣A為常數矩陣。我們提取矩陣A的特征值λ，其特征值組成互為共軛的數對，如圖3所示。λ的實(shí)部和虛部分別對應格式的數值色散和數值耗散特性。如Van den Abeele在[5]中指出，我們找到不同波數下對應的λ值，從而得出半離散系統的數值耗散律和數值色散率。而我們知道，連續系統的理想數值耗散律恒為零，而數值色散律為波數的線(xiàn)性函數。通過(guò)與連續系統的理想值比較我們便可得出格式在數值色散律和數值耗散率對應于不同波數的誤差，圖4和圖5給出了誤差計算方法的示意。圖中橫坐標K為波數，而格式波數的范圍為[0，+H]，其中為每個(gè)SV中CV的數量，而H為負值，表明降階SV方法所降低的階數，而ε表示數值格式與理想情況相比較的誤差。結果比較在圖6-9中給出。
我們可以看出對于大波數情況下，降階SV方法相較傳統SV方法在數值耗散和數值色散特性上均有優(yōu)勢。對于二階情況SV42與SV32相較傳統的SV22方法，在小波數時(shí)誤差甚至小兩個(gè)數量級以上。而對于三階情況SV43也比SV33方法在數值耗散方面對于大波數成分有較大優(yōu)勢。
5 結語(yǔ)
我們成功發(fā)展了降階的SV格式，推導了降階SV格式在流量計算中為得到Gauss積分點(diǎn)處狀態(tài)變量的插值公式。隨后我們在文章中驗證了降階SV格式與傳統SV格式的精度，證實(shí)了其可以達到所期望的計算精度。隨后我們對其穩定性做了分析。研究證明，新格式在數值色散和數值耗散方面較傳統SV格式有明顯優(yōu)勢。我們仍然會(huì )在今后對新格式特性的進(jìn)一步更深入的研究。
參考文獻
[1]閻超，于劍，徐晶磊，等.CFD模擬方法的發(fā)展成就與展望[J].力學(xué)進(jìn)展，2011（5）：562-589.
[2]Z.J.Wang.Spectral（Finite）Volume Method for Conservation Laws on Unstructured Grids： Basic Formulation.J.Comput[J].Phys，2001，178：210-251.
[3]Z.J.Wang，Y.Liu.Spectral（Finite）Volume Method for Conservation Laws on Unstructured Grids II： Extension to Two-Dimensional Scalar Equation.J.Comput[J]. Phys，2002，179：665-697.
[4]Y.Liu，M.Vinokur，Z.J.Wang.Spectral（Finite）Volume Method for Conservation Laws on Unstructured Grids V：Extension to three-dimensional systems. J.Comput[J].Phys，2006，212：454-472.
[5]K.Van den Abeele.Development of high-order accurate schemes for unstructured grids[J].Phd thesis in Vrije Universiteit Brussel，2009（3）.
[6]K.Van den Abeele，C.Lacor.An accuracy and stability study of the 2D spectral volume method.J.Comput[J].Phys.，2007，226（1）：1007-1026.
[7]M.Zhang，C.-W.Shu.An analysis of and a comparison between the discontinuous Galerkin and the spectral finite volume methods.Comput[J].Fluids，2005（34）：581-592.
[8]Z.J.Wang，Y.Liu：Extension of the spectral volume method to high-order boundary representation.J.Comput[J].Phys，2006（211）：154-178.
[9]H.T.Huynh：A flux reconstruction approach to high-order schemes including discontinuous Galerkin methods[J].AIAA Paper，2007，10（105）：219.
[10]閆超.計算流體力學(xué)方法及應用[M].北京：北京航空航天大學(xué)出版社，2006.
[11]A.Harten.High resolution schemes for hyperbolic conservation laws，J.Comput[J]. Phys，1983，49（3）：339-357.
[12]A.Harten，B.Engquist，S.Osher， et al Uniformly high order essentially non-oscillatory schemes III，J. Comput[J].Phys，1987，71（2）：231-303.
[13]M.Zhang，C.-W.Shu：An analysis of and a comparison between the discontinuous Galerkin and the spectral finite volume methods.Comput[J].Fluids，2005，34（4-5）：581-592.
[14]J.A.Ekaterinaris.High-order accurate，low numerical diffusion methods for aerodynamics， Progress in Aerospace Science，2005，41（3）：192-300.
