<ol id="ebnk9"></ol> 
摘 要：文章中提出了降階的Spectral Volume 方法 （SV方法） 的概念，對方法的公式進行推導，并成功將其應用與一維線性傳遞方程。對降階的SV方法的穩定性進行了Fourier分析，并與傳統的SV方法進行了對比。實驗證明，降階的SV方法在初值連續情況下，在線性傳遞方程上可以達到所期望的數值精度，并且在數值色散和數值耗散特性上較傳統SV方法有顯著提高。
關鍵詞：SV方法 一維守恒律 降階
中途分類號：R445文獻標識碼：A文章編號：1674-098X（2014）10（b）-0211-04
1 研究背景
計算流體力學（CFD）從創立至今得到了飛速發展，在工業界起著不可替代的作用。然而目前工業界所用到的CFD計算通常只能達到二階精度，遠遠達不到人們的需要，并且對網格以及內存和CPU的要求超出了當今計算機所能承受的范圍。另外，運算區域通常包含復雜單光滑的流體結構，比如渦旋與剪切流層的相互作用區，以及同時包含激波與其他復雜流動的區域等。因此，高精度高分辨率格式的研究一直是CFD領域中非?；钴S的一個課題，發展高精度高分辨率的計算格式成為CFD工作者的一大研究方向。為采用高精度格式不但可以降低對網格規模的苛刻要求，而且能夠正確分辨其中復雜的流動現象。對網格點數目的要求低，具有減少計算機內存需求和CPU時間消耗的優勢。其中比較著名的是Harten提出的高階精度的ENO（essentially non-oscillatory） 格式，隨后Liu，Osher和Chen發展了原有的ENO格式，提出了有顯著改進的WENO（weightedessentially non-oscillatory）格式，并獲得了迅速的發展，并得到了一系列的研究。ENO/WENO格式善于捕捉物理間斷，其應用迅速發展到多相流領域。
新型格式發展的目標之一是發展可以應用于非結構網格的可靠的高階格式。理論和實踐都表明，大多數高精度、高分辨率格式的性能都嚴重地依賴于網格的光滑性，而工業界所需要模擬的物體通常具有復雜的幾何外形（圖1），CFD應用最多的航空領域所需要的幾何外形更是相當復雜，而這給結構網格的生成帶來很大困難，一種直接有效的解決辦法是采用非結構網格。然而由于穩定性等因素的影響，對于三維有限體積法，非結構算法模板所需單元往往較多，而過多的模板單元會帶來邊界處理、內存占用、編程復雜等多方面的困難。因此實際應用中大都只使用二階精度的非結構算法，而這遠遠滿足不了一些工業領域對計算精度的要求?？梢詰门c非結構網格的高精度、高分辨率、低耗散的實用數值格式的發展成為當今CFD的重要任務之一。以此為目標的一些高階格式已經發明出來，其中比較著名的是Cockburn發展的間斷Galerkin有限元方法（DG），間斷 Galerkin有限元方法的插值模板小，各單元之間僅通過界面通量計算相聯系，因而比較適合在復雜外形。而另一種成功的格式是Wang在自2002年發展的應用于非結構網格的算法SV方法。SV方法是一種高精度高分辨率低耗散的數值格式。同間斷Galerkin有限元方法一樣，SV方法可以達到相同的精度。然而相比之下，SV方法比間斷Galerkin有限元方法有著更好的穩定性（可以允許更大的CFL數）。與傳統的有限體積法相比，對計算機內存和CPU的要求更低。所以SV方法有著非常大的能在將來被應用于工業界的潛力。同時，SV方法的穩定性和可行性還需要進一步驗證和改進。
2 SV方法介紹
SV方法，是比較新提出的另一種可以應用于完全非結構網格的高精度、高分辨率、低耗散的守恒型格式。Wang于2002年提出一維SV方法的概念并隨后將其發展到二維和三維。SV方法是經典有限體積法的延伸。它提出了控制體積（Control Volume）的概念，從而避免的復雜模板（stencil）的使用，從而使格式可以輕松應用與非結構網格上。經典的有限體積法中，變量在每個網格中的值被用來進行多項式插值。為了達到相應的精度，由相鄰多個單元格組成的模板被用來進行插值。而非結構網格使得模板中相鄰網格的選取成為一大難題。不同于經典的有限體積法，SV方法把網格中的每一個單元格進一步分成若干控制體積（CV），對每個控制體積中變量的平均值加以利用來實現多項式插值從而達到所要求的精度。SV方法的優勢在于其在保證相同精度的同時免除了對周圍網格的依賴，從而降低了對網格的要求。此外，由于不同單元格之間的信息交流少，SV方法非常適宜并行運算的實現。
有限體積法（包括ENO和WENO方法）利用相鄰多個自由度（模板）內狀態變量的值進行多項式構建。而在SV方法中，我們把這個單個的單元格叫做Spectral Volume（SV）。我們只利用此SV中的信息進行高階多項式的構建，為了得到所需要的自由度，每個SV被劃分為更小的體積單元 Control Volume（CV）。而方法的精度取決于每個SV劃分為CV的數量。
對于二階的SV方法，流量的計算是線性的，所以并不復雜。對于高階的情況 （k>2），重構過程如下。在下一個時間步上，每個CV中狀態變量的平均值各自更新，在SV內部，狀態變量在床CV邊界時時連續的，所以我們并不需要Riemann求解器，流量的值我們可以通過插值多項式直接到的解析值。當兩個相鄰CV的邊界也是 SV邊界時，我們需要通過插值得到邊界上Gauss求積點處狀態變量的值，然后通過 Riemann求解器來計算流量的值。這樣方法保證了計算的簡潔以及高階特性。
總的來說，假設我們重構出至多k-1 階多項式，k階SV方法計算的一般過程是通過插值計算每個Gauss積分點處狀態量的值；利用k階高斯求積形式和Riemann 求解器計算單元邊界上面通量的積分。在單元格內部CV邊界上，由于我們假定狀態變量在單元格內是連續的，通量的值我們直接應用重構多項式的解析值；利用TVD的Runge-Kutta格式進行時間推進。
如同有限體積方法一樣，在計算過程中SV方法利用狀態變量的平均值。然而SV方法并不像有限體積法那樣利用在整個單元格內平均，它利用的是狀態變量在每個控制體積內的平均值，通過CV的劃分來增加自由度，從而達到所需要的高階。關于SV方法的具體操作以及一維和多維的各種算例可以在文章中找到。
之前提到的間斷Galerkin有限元方法和SV方法在特性上有許多相似之處，同樣是緊致類方法，從而適于并行運算，能達到高階精度，由于在單元格邊Riemann求解器的應用，它們是完全守恒型的數值計算格式，并且伴隨著TVD或者TVB限制器的應用，兩種方法都非常適于應用于非結構網格和復雜的幾何外形。自然而然，兩種方法的研究比較出現在很多文章中。由于考慮局部的細致信息，SV方法被認為相對應間斷Galerkin有限元方法可以更加細致的捕捉間斷。此外間斷Galerkin有限元方法會隨著計算階數的升高而導致相應的CFL數減小，穩定性會隨之下。而Zhang和Shu通過理論以及算例指出兩種方法可以達到相同的精度，SV方法可以允許更大的CFL數，然而在相同的網格上SV方法的誤差大于間斷Galerkin有限元方法。
3 降階SV方法
SV方法的穩定性受很多因素影響，我們提出的一種對其可能的改進是降階的SV 方法。由于SV方法是緊致型格式，所有信息的處理幾乎都在單元格自身內部進行，這使得其格式的改造變得相對容易。一種可能的構想是減小SV方法基底的維度。以基本的一維情況為例，理論上k階的SV方法需要將SV劃分成k個CV并以k-1階多項式進行近似
（1）
在傳統的SV方法中，基底的階數與 CV的個數是相等的，從而保證方法的精度。在這里我們提出，可以嘗試以降低多項式基底的維度來改善SV方法在穩定性方面的表現。
假如SV被分為k個CV并以k-1-n階多項式（n≤k-1）近似，其理論精度值為 k-n。
（2）
如圖2給出了一維情況下以二階多項式描繪三個CV中狀態變量值的示意圖。它的實際精度以及在具體算例中的表現需要進行驗證。
另外可以期待的一點是此方法對格式穩定性的影響，可以通過Fourier分析法研究格式的數值色散和數值耗散特性，k個CV并以k-1-n階多項式近似的降階的SV方法被期待比傳統的k-n階SV方法有更好的穩定性，即在數值色散和數值耗散方面有更好的表現。
SV方法插值邊界處Gauss積分點上狀態變量的值的關系式是：
（3）
其中向量指的是SV中每個CV上狀態變量的平均值，其維度是此單元格內CV的數量，向量q表示邊界上每個的Gauss求積點處狀態變量的值，其維度是此單元格內 Gauss 求積點的數量。
相應的構建方法是：
（4）
其中為的維度，為SV內部近似狀態變量的重構多項式，它是基底多項式的線性組合
（5）
其中在某一固定的時間步上，為常數，代表基底線性組合的系數。是多項式基底的維度。而q的每個元素則就是重構多項式在CV邊界處對應Gauss積分點上的取值。
（6）
指的是Gauss積分點的數目，其中指第k個Gauss積分點的坐標。從而矩陣運算可以從中提取出。我們得到多項式取值的矩陣P和基底在每個CV上的平均值矩陣L。
（7）
（8）
分析一下維度我們得到
（9）
對于不同的劃分方式，我們只需要根據基底和CV的幾何分布得到矩陣P和L便可以進行運算。但這里我們假設，L是可逆方陣。如果我們縮減基底的維數，L便不再可逆，q的表達式會有相應變化，這也將是新的自適應階數SV方法的創新之處。
（10）
此時我們在計算公式中以L的假逆矩陣代替其逆矩陣。而我們定義其假逆矩陣：
（11）
這里M是一個正定矩陣，它的階數等于每個SV中CV的數量。在這種定義之下，我們重新定義了一種基于矩陣M的數量積： 對于任意向量a和b，我們有
（12）
與此同時，利用幾何中投影的定義，我們也定義了一種基于矩陣 M 的投影：
（13）
在具體操作中我們發現，M矩陣的選取要滿足許多要求。當我們利用降階的SV方法時其中非常重要的一點是如下守恒定律必須要滿足：
（14）
其中是保存狀態變量平均值的向量，其中的每個元素保存了同一個SV中指定CV上狀態變量的平均值，是所研究 SV中CV的數量。這個條件的物理意義是，對一個SV中每個CV上的狀態變量平均值投影到矩陣L的像空間，投影前后保證整個SV中狀態變量的積分總量不變。事實上，如果這一關系沒有滿足，則原方程的守恒形式將被打破。之后的數值測試結果表明，如果以上提到的守恒定律沒有被滿足，計算的結果是完全錯誤的。滿足這一條件的矩陣M也不是唯一的，利用簡單的矩陣運算，我們得到一個符合守恒定律條件的如下矩陣M：
（15）
數值測試
在測試的初期，我們通過編寫如下基本的一維線性輸運方程 Matlab 程序驗證 SV 方法精度：
（16）
方程的輸運速度設為a=1。
為了方便驗證且不失一般性，我們提出T=60s以保證。
我們首先驗證初始條件連續
的情況下SV方法的精度。我們采用初始條件：
（17）
測試結果顯示在表1中。其中自由度指的是區間[0，3]被劃分成的CV的數量。在命名方法上，我們采用：SV+每個SV中CV數+模擬多項式的階數。如SV43表示每個SV劃分為四個CV，但降一階，采用三階精度的多項式近似。若每個SV中CV數=模擬多項式的階數，則格式為傳統的SV方法。在測試中我們采用三階Runge-Kutta時間積分方法。
我們可以從表1中看出，理論上達到二階的數值格式，即SV22，SV32及SV42方法，在低自由度的網格上由于誤差太大，格式很難達到所期望的精度，而在過高自由度的網格上，由于數值計算誤差，所得到的精度可能與理論值有所偏差。但總體來說，降階的SV方法在允許的誤差范圍內可以達到理論預期值。
4 格式穩定性研究
我們通過Fourier分析法對格式的穩定性進行研究自適應SV方法較之于傳統SV方法的優勢。這里我們之研究格式空間離散帶來的誤差而不考慮時間推進格式。所以這里我們研究計算中的半離散的數值格式：
（18）
其中即我們所研究的狀態變量。通過對方程系統的線性化我們得到矩陣A。這里我們研究最基本的一屆線性傳遞方程，從而矩陣A為常數矩陣。我們提取矩陣A的特征值λ，其特征值組成互為共軛的數對，如圖3所示。λ的實部和虛部分別對應格式的數值色散和數值耗散特性。如Van den Abeele在[5]中指出，我們找到不同波數下對應的λ值，從而得出半離散系統的數值耗散律和數值色散率。而我們知道，連續系統的理想數值耗散律恒為零，而數值色散律為波數的線性函數。通過與連續系統的理想值比較我們便可得出格式在數值色散律和數值耗散率對應于不同波數的誤差，圖4和圖5給出了誤差計算方法的示意。圖中橫坐標K為波數，而格式波數的范圍為[0，+H]，其中為每個SV中CV的數量，而H為負值，表明降階SV方法所降低的階數，而ε表示數值格式與理想情況相比較的誤差。結果比較在圖6-9中給出。
我們可以看出對于大波數情況下，降階SV方法相較傳統SV方法在數值耗散和數值色散特性上均有優勢。對于二階情況SV42與SV32相較傳統的SV22方法，在小波數時誤差甚至小兩個數量級以上。而對于三階情況SV43也比SV33方法在數值耗散方面對于大波數成分有較大優勢。
5 結語
我們成功發展了降階的SV格式，推導了降階SV格式在流量計算中為得到Gauss積分點處狀態變量的插值公式。隨后我們在文章中驗證了降階SV格式與傳統SV格式的精度，證實了其可以達到所期望的計算精度。隨后我們對其穩定性做了分析。研究證明，新格式在數值色散和數值耗散方面較傳統SV格式有明顯優勢。我們仍然會在今后對新格式特性的進一步更深入的研究。
參考文獻
[1]閻超，于劍，徐晶磊，等.CFD模擬方法的發展成就與展望[J].力學進展，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]閆超.計算流體力學方法及應用[M].北京：北京航空航天大學出版社，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.
<ol id="ebnk9"></ol>