摘 要：文章中提出了降阶的Spectral Volume 方法 （SV方法） 的概念，对方法的公式进行推导，并成功将其应用与一维线性传递方程。对降阶的SV方法的稳定性进行了Fourier分析，并与传统的SV方法进行了对比。实验证明，降阶的SV方法在初值连续情况下，在线性传递方程上可以达到所期望的数值精度，并且在数值色散和数值耗散特性上较传统SV方法有显著提高。
关键词：SV方法 一维守恒律 降阶
中途分类号：R445文献标识码：A文章编号：1674-098X（2014）10（b）-0211-04
1 研究背景
计算流体力学（CFD）从创立至今得到了飞速发展，在工业界起着不可替代的作用。然而目前工业界所用到的CFD计算通常只能达到二阶精度，远远达不到人们的需要，并且对网格以及内存和CPU的要求超出了当今计算机所能承受的范围。另外，运算区域通常包含复杂单光滑的流体结构，比如涡旋与剪切流层的相互作用区，以及同时包含激波与其他复杂流动的区域等。因此，高精度高分辨率格式的研究一直是CFD领域中非常活跃的一个课题，发展高精度高分辨率的计算格式成为CFD工作者的一大研究方向。为采用高精度格式不但可以降低对网格规模的苛刻要求，而且能够正确分辨其中复杂的流动现象。对网格点数目的要求低，具有减少计算机内存需求和CPU时间消耗的优势。其中比较著名的是Harten提出的高阶精度的ENO（essentially non-oscillatory） 格式，随后Liu，Osher和Chen发展了原有的ENO格式，提出了有显著改进的WENO（weightedessentially non-oscillatory）格式，并获得了迅速的发展，并得到了一系列的研究。ENO/WENO格式善于捕捉物理间断，其应用迅速发展到多相流领域。
新型格式发展的目标之一是发展可以应用于非结构网格的可靠的高阶格式。理论和实践都表明，大多数高精度、高分辨率格式的性能都严重地依赖于网格的光滑性，而工业界所需要模拟的物体通常具有复杂的几何外形（图1），CFD应用最多的航空领域所需要的几何外形更是相当复杂，而这给结构网格的生成带来很大困难，一种直接有效的解决办法是采用非结构网格。然而由于稳定性等因素的影响，对于三维有限体积法，非结构算法模板所需单元往往较多，而过多的模板单元会带来边界处理、内存占用、编程复杂等多方面的困难。因此实际应用中大都只使用二阶精度的非结构算法，而这远远满足不了一些工业领域对计算精度的要求。可以应用与非结构网格的高精度、高分辨率、低耗散的实用数值格式的发展成为当今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.
