基于BVD原理的高保真空间重构方法

肖 锋

(东京工业大学, 日本 东京 152-8550)

0 引 言

可压缩流体的运动变化规律是空气动力学的主要研究内容。伴随高马赫数流动产生的激波等不连续现象使得可压缩流动的流场结构更为复杂，并同时包含连续与间断解，给理论和实验研究带来很多本质性的困难。在过去的数十年中，由于计算机硬件的飞速发展，以及相关领域实际应用的巨大需求，计算流体力学作为最活跃的科学研究领域之一，取得了重大进展， 并已成为空气动力学非常重要的研究手段。

针对描述可压缩流体的欧拉方程，至今已提出了各种数值算法。其中，以具有严格守恒特性的Godunov格式[1]为基础开发的各类高阶守恒格式已成为求解可压缩流体运动的主流算法。基于Godunov格式的高阶守恒格式一般包括重构和演化两个重要步骤。其中，演化是根据重构得到的离散变量，利用控制方程的物理性质计算相应时间步长内的数值通量，进而预报下一时刻的物理量变化。关于欧拉方程及其黏性扩展情形下的数值通量计算，可参阅有关专著及论文[2-5]。本文集中讨论如何进行物理变量的空间重构。

作为构造高阶Godunov方法的关键，空间重构一直是关注的焦点。至今已提出了针对不同空间离散框架的各种重构方法[6]，包括有限体积或有限差分格式[7-19]、紧致格式[20-22]等基于单个单元自由度的方法， 也有间断伽列金方法[23-24]、间断伽列金/有限体积混合方法[25-27]、通量重构方法[28-30]、多矩有限体积方法[31-34]等包含多个局地自由度的高精度方法，并已被用于求解各类实际问题。这些方法通常采用高阶多项式进行重构，在用于计算解较为光滑的问题时，能够得到较高收敛率和较精确的数值结果。然而，当求解的物理问题含有激波或多相复杂介质界面等强间断时，必须使用非线性限制映射或限制器来抑制虚假的数值振荡。其中，使用WENO(weighted essentially non-oscillatory)[12-13]思想构造的高分辨格式既可有效地抑制数值振荡又能对光滑解获得较高的收敛精度，在可压缩流体的数值模拟中得到广泛应用。现有算法中设计非线性限制器的基本做法是，在解较光滑的区域尽量保持高阶多项式的性质，而在解出现间断处采用更平直或较低阶的重构函数。这类做法通常在抑制数值振荡与控制数值耗散之间，在计算鲁棒性和求解精度等方面无法兼顾，很难对光滑和间断解同时给出令人满意的计算结果。在实际应用中主要表现为：(1)存在显著的数值耗散，会抹平接触间断和较小尺度的流场结构；(2)多数非线性权函数无法使重构函数在光滑情形下完全恢复到原来的常系数高阶多项式插值函数，且依赖于一些人为确定的参数，给实际应用带来不便。作为相关领域长期以来一个未能很好解决的问题，如何构造能同时精确捕捉光滑和间断解的高保真空间重构方法，一直是可压缩流体数值方法研究的热点。

为了更好地模拟包含间断解的可压缩复杂流动，我们近年提出了以减小空间重构在网格单元边界上的变差为原则设计数值格式的思想，即，BVD(Boundary Variation Diminishing)原理。根据此原理，我们发展了一类可用于求解单相和多相可压缩流体的高保真数值算法[35-46]。这些算法能有效克服其他现有高分辨格式存在的问题，可以同时精确捕捉光滑和间断解。各种标准算例结果表明，BVD格式在模拟各种具有强间断的可压缩流体运动方面，与其他的现有格式相比具有明显的优势。 利用BVD原理进行空间重构可望成为构建可压缩流体新型高保真数值方法的有效途径。

本文对BVD格式研发的已有工作做一个简要综述，内容包括BVD的基本思想介绍、几个具有实用意义BVD算法、有代表性的标准算例结果，以及简短总结和展望。

1 BVD的基本思想

1.1 Godunov方法概要

我们考虑双曲标量守恒律:

通过以上步骤(I)和(II)获得边界数值通量后，一般可以通过线方法，采用常微分方程的数值解法对半离散方程(2)进行时间积分,得到下一时刻的数值解。

至今的高分辨数值方法研究，基本集中于开发具备更好数值性能的空间重构(I)和黎曼算子(II)。

1.2 BVD方法的基本思想

从近似黎曼算子的基本形式(3)可以看出，数值通量一般由中心格式(右端第一项)和数值耗散/数值黏性(右端第二项)两部分组成。如果在构造数值格式时，能够减小数值耗散项，将有效降低格式的数值耗散，从而改善计算结果。

由数值耗散项的表达式不难看出，减小数值耗散可从两方面入手，即:

(B) 减小空间重构得到的边界变差值

各类黎曼算子的研究开发在一定意义上可以看作是为实现目标(A)所做的努力。围绕在保证计算稳定的条件下，如何控制减小数值耗散系数，提出了各类近似黎曼算子。这方面的研究可参阅相关专著[2]。

为了在空间重构过程中有效地减小数值耗散，我们提出以下设计空间重构格式的一般性原理。

如果在开发空间重构算法时能有意识地根据BVD原理，尽可能地减小边界变差值，则可以有效抑制格式的数值耗散，改进数值计算结果。我们注意到，文献[47]通过边界变差最小化的变分约束关系构造了紧致有限体积方法。

2 BVD混合格式的构建

基于BVD原理，我们提出了构建高保真数值格式的基本框架，并发展了一类基于BVD原理的混合格式。BVD格式主要包括BVD许容重构函数集和BVD选择算法两个部分。

2.1 BVD许容重构函数集

流体力学问题解的空间分布极为复杂，既包括各种尺度的涡和波动等相对光滑的流场结构，也包括激波、接触间断、不同流体间的物质界面等不连续结构。很难用一类插值函数准确描述各种不同的流场结构。在以往的研究中，各种重构函数和重构方法被相继开发出来，并用于针对不同应用问题的高分辨格式。我们把一些有代表性的重构函数分为以下三类BVD许容函数。

2.1.1 第一类：常系数多项式

这类重构函数用于逼近光滑解。包括各阶泰勒多项式、紧致格式[48-49]以及各种频散/耗散优化格式[50-54]等。作为例子，我们具体给出针对网格单元Ωi的四阶多项式:

其中系数ck(k=0,1,…,4)通过以下约束条件求出:

(5)

为书写简便，我们在此采用等网格距，Δxi=Δx。由此可以直接得到:

2.1.2 第二类：非线性系数多项式

在高分辨激波捕捉格式的发展过程中，为避免高阶多项式在间断解附近出现的非物理振荡，提出了许多非线性限制投影算法。例如，采用TVD限制器的2阶MUSCL格式[7]，可将网格单元边界左右两侧的变量值写为：

如果采用WENO类算法，则可将插值函数写成非线性加权系数多项式的形式。以下给出5阶WENO格式计算网格单元边界左右两侧变量值的算式。

其中,

关于非线性权重系数ωj(j=0,1,2)的计算方法，至今已提出了各种方案[12-19,24]。

2.1.3 第三类:非多项式类插值函数

采用某些具有单调特性的非多项式类插值函数能够避免空间重构中的数值振荡。例如，采用双曲函数、有理函数或对数函数[55-57]可以有效地抑制数值振荡。文献[58-61]使用双曲正切函数拟合跃阶型间断，并发展了一类可用于捕捉移动界面的代数VOF(Volume of Fluid)方法[62]，称为THINC(Tangent of Hyperbola for INterface Capturing)方法。和其他Sigmoid函数一样，THINC重构函数具有严格的单调性，并能很好地模拟激波、接触间断以及多相流中的物质界面等具有跃阶分布的变量。

由此，可以将通过THINC重构得到的网格单元左右两端的值写为:

这里，

B=eθβ(2C-1)

在实际应用中，可以通过调节陡度参数β的值，使THINC重构函数适用于各类数值解分布。数值频散/耗散分析表明，如果β取值于1和1.3之间，则THINC重构可以表示一类MUSCL格式[37]。取较大β值，则可以很好地拟合跃阶间断分布，并有效地抑制数值扩散。

对于多维情形，可将THINC重构函数表示为:

其中，Pi(x,y,z)+di=0为界面方程。di通过与式(9)类似的多维体积分约束关系计算。Pi(x,y,z)可采用高阶多项式，以更准确地描述界面的几何形状。关于多维高阶THINC重构可参阅相关文献[60-61]。

从以上三类许容函数中选取适当的候补函数构成BVD许容函数集，记为

ξ=1,2,…,Ξ

这里Ξ表示BVD许容函数集里所包含的候补函数的个数，通常Ξ≥2。

一般根据流场结构的主要特征以及希望改善的数值解性质选取候补函数。例如，针对以间断解为主的应用问题，可以选取非线性限制多项式(第二类)和THINC函数(第三类)构成BVD许容函数集[35-36,43]。或选取采用不同β值的THINC函数构成BVD许容函数集[37]。针对以声场或旋涡为主的应用问题，为了充分利用高阶常系数多项式的良好性质，可以避开非线性限制多项式，只选取高阶常系数多项式(第一类)和THINC函数(第三类)构成BVD许容函数集[38-42,44-46]。研究表明，通过构造合适的BVD算法，即使不采用非线性限制多项式，也能有效地抑制数值振荡，如后叙的PnTm-BVD格式。

2.2 BVD选择算法

我们对各类BVD许容函数用于连续和间断分布重构时产生的BV或TBV值进行了分析及数值实验，得到以下主要结论：

(1) 对于光滑分布，采用第一类许容重构函数(常系数多项式)能够得到更小的BV或TBV值。且多项式阶数越高，边界变差值越小。与Weierstrass多项式逼近定理一致。

(2) 对于跃阶型间断解，在第二类许容重构函数(带限制器的非线性系数多项式)中，高阶格式的TBV值小于低阶格式的TBV值。数值耗散较大的格式产生较大的TBV值。与基于非线性系数多项式重构相比，采用较大β值的THINC重构可以获得更小的TBV值。

(3) 对间断进行插值时，THINC函数产生的BV只出现在空间较窄区域，而基于多项式的重构函数则会导致较宽区域内出现明显的BV。

由此，我们提出了各种选择算法。下面列举几个有较强实用意义的BVD算法。

2.2.1 BVD算法Ⅰ

BVD许容函数集包括两个候补函数，即

2.2.2 BVD算法Ⅱ

采用包含两个候补函数的BVD许容函数集，通过以下步骤确定重构函数。

2.2.3 BVD算法Ⅲ

针对包含两个候补函数的BVD许容函数集，通过以下步骤确定重构函数。

(2) 选取使得边界总变差最小的函数作为最终的重构函数，

针对多个候补函数构成的BVD许容函数集(Ξ≥3)，算法III可写成更一般的形式:

文献[43]采用三个候补函数建立了非结构网格的BVD格式，能很好地捕捉间断和涡结构。

2.2.4 BVD算法Ⅳ

当BVD许容函数集包括3个以上候补函数时，可通过多步BVD算法确定重构函数。其一般步骤如下:

for (ξ=2;ξ≤Ξ;ξ++)

{

}

其中，BVDξ(·,·)表示基于BVD原理，针对两个候补函数的选择算法，如上述的算法Ⅰ至算法Ⅲ。

(2) 对ξ=1，…,m-1循环计算

和

j=i-1,i,i+1。

文献[38-40]分别给出了P4T2-BVD、P6T3-BVD、P8T3-BVD、P10T3-BVD、P12T3-BVD和P14T3-BVD等高阶BVD格式。这些格式对光滑解可以得到对应的常系数多项式的最高收敛阶数，亦即，PnTm-BVD格式具有n+1阶。同时，能有效避免数值振荡。与其他现有的高阶方法不同，这类BVD格式完全不需要使用传统的非线性限制映射来消除数值振荡。

使用ADR(Approximate Dispersion Relation)方法[63]可以分析数值格式的频散和耗散性质。图1显示各阶PnTm-BVD格式的数值频散和耗散特征。为了比较，这里还显示了WENO-JS[13]、WENO-M[14]、WENO-Z[15]以及基于各阶格式对应的常系数多项式的线性格式的频散/耗散关系。对于所有波数，PnTm-BVD可以回复其对应高阶线性格式的频散/耗散关系，而采用非线性限制多项式的高阶格式则在高波数区明显偏离线性格式。

(a) 频散特征

(b) 耗散特征

此外，文献[44]还提出了以4次常系数多项式和采用可变β值的THINC作为候补函数的BVD格式。其中，网格单元Ωi的THINC函数的β值通过以下算式给出:

j=i-1,i,i+1

该格式称为P4Tβv-BVD格式。

2.3 多维BVD格式的构建

上述构建一维BVD格式的思路原则上都可以推广到多维情形。对于结构网格可以采用方向分离方法，通过一维格式进行空间重构。对于非结构网格，文献[43]用重构函数沿网格单元各边界的积分值计算BV。针对网格单元Ωi和与其相邻的单元Ωij，沿单元边界Γij的BV值由下式计算，

需要指出，BVD算法在同一网格单元上要对所有候补函数进行重构，并储存相应的BV或TBV值，会增加计算量和存储量。然而，由于各网格单元间的计算过程相互独立，BVD算法不会给并行处理带来不利影响。

3 数值结果

我们利用标量守恒律和欧拉方程的各种标准算例，对BVD算法进行了系统的检验，并与其他现有方法进行了比较。

3.1 精度检验

对于光滑解，BVD格式一般都能从候补函数中选择高阶重构，从而获得高收敛率。以下列举采用不同阶数常系数多项式作为候补函数的BVD格式的收敛检验算例。文献[38-40]分别针对平流和欧拉方程进行了网格收敛率检验。对于一维平流方程，选用以下具有一定光滑度的初始条件，

对于二维欧拉方程的密度扰动传播算例,初始条件设为:

(x,y)∈[-1,1]×[-1,1]

表1给出各阶PnTm-BVD格式的误差及收敛率，均能达到相对应的常系数多项式的阶数。其误差的绝对值与采用常系数多项式作为重构函数的结果基本一致。

表1 PnTm-BVD格式的数值误差及精度检验

3.2 Jiang-Shu平流实验

作为检验数值格式是否能够同时捕捉连续和间断解的算例，文献[13]设计了一个既包括光滑分布又包括间断解的平流初始条件。

图2分别显示了经过一百万步后WENO-JS[13]，WENO-Z[15]和P4T2-BVD格式[38]的计算结果。可以看出，由于数值误差的不断积累，WENO-JS和WENO-Z格式已无法保持初始分布的特征。而P4T2-BVD格式则能很好保持初始分布，对连续和间断解都能得到高保真的计算结果。P4T2-BVD格式既能有效地避免数值振荡，又可以消除数值耗散，使解的结构在长期积分后仍能可靠再现。这是现有其他方法很难做到的。

图2 Jiang-Shu平流标准算例在t=2000(1×106步)时的计算结果. 从左至右分别为五阶WENO-JS[13]格式、WENO-Z[15]格式和P4T2-BVD格式[38]的计算结果，网格单元数为400.

3.3 一维欧拉方程的数值实验

对于一维欧拉方程，我们利用各类标准算例对BVD格式进行了检验，在抑制数值振荡的同时，对间断和光滑解均能得到高保真的计算结果。以下给出几个具体算例的结果。

3.3.1 双激波相互作用[64]

本算例通过两个激波的相互作用及边界反射产生复杂的流场结构，包括激波、接触间断和膨胀波，被广泛用于检验数值方法是否产生数值振荡，以及对各类间断和连续解的捕捉能力。图3给出了P4T2-BVD以及P4Tβv-BVD格式的计算结果。可以看到BVD格式能够精确地得到各类特征的数值解。在图4中我们进一步显示9阶和11阶BVD格式的计算结果，高阶BVD格式也能有效抑制数值振荡，保持良好的鲁棒性，并得到高保真的数值计算结果。在上述结果中，BVD格式能够有效地控制数值耗散，可使最左端的接触间断厚度保持在3至4个网格之内。这是现有基于多项式重构及非线性限制映射的高分辨格式难以实现的。

图3 Two interacting blast waves标准算例在t=0.038时刻的P4T2-BVD和P4Tβv-BVD格式的计算结果(密度), 网格单元数为400.

图4 Two interacting blast waves标准算例在t=0.038时刻的P8T3-BVD(左)和P10T3-BVD(右)格式的计算结果(密度),网格单元数为400.

3.3.2 激波与密度扰动相互作用[65]

为了进一步检验BVD格式对高频周期波的捕捉能力，我们计算了激波与密度扰动相互作用的算例。图5显示不同阶数BVD格式的计算结果。可以看到，不同阶数的BVD格式都能抑制数值振荡，并在光滑区正确地选择多项式作为重构函数。随着多项式阶数的提高，BVD能有效减小数值耗散，更精确地捕捉高频波动[40]。

(a) P4T2-BVD[38]和P6T3-BVD格式[38]的计算结果

(b) P8T3-BVD[38]和P10T3-BVD格式[38]的计算结果

3.3.3 Le Blanc 激波管问题[66]

作为一个较困难的标准算例，Le Blanc激波管问题的解包括非常强的激波和稀疏波，经常用来检验算法的鲁棒性[66]。图6给出了P4T2-BVD的计算结果(密度)。没有出现明显的数值振荡。并能很好再现稀疏波、接触间断和激波[46]。数值实验表明，随着网格分辨率的提高，右行激波能够收敛到相应精确解的位置。

图6 Le Blanc一维激波管标准算例在t=6时刻的P4T2-BVD格式的计算结果(密度), 网格单元数为800.

3.3.4 刚性爆轰波[67]

与现有的数值方法相比，BVD格式能够消除数值耗散，精确地捕捉不连续界面。图7显示了包括刚性化学反应的爆轰波计算结果。P4T2-BVD格式能有效控制反应面的温度跃阶厚度，从而保证反应面始终以正确的速度传播[42,46]。

图7 刚性爆轰波标准算例在t=π/5时刻的P4T2-BVD格式的计算结果(密度),网格单元数为200.

3.4 二维欧拉方程的数值实验

3.4.1 双马赫反射[64]

作为检验二维算法的标准算例，双马赫反射问题自提出以来，在相关文献中得到广泛应用。在本算例中，沿两个反射激波间的滑移线会出现开尔文-亥姆霍兹不稳定，进而产生涡列。在数值计算中，这些涡列是否出现和发展取决于数值耗散的大小。图8分别给出了P4-T2-BVD和P14-T3-BVD格式的计算结果。由于BVD原理旨在减小格式的数值耗散，以此为基础设计的格式都能更多地保留流场的细微结构，捕捉到更小尺度的流体运动。与其他五阶格式相比，P4-T2-BVD格式的结果显示了充分发展的涡列。随着精度提高，十五阶的P14-T3-BVD格式可以捕捉到更多旋涡结构[44]。

3.4.2 爆轰波扰流问题[68]

二维爆轰波的90°凸角绕流在绕角下流区域会产生密度和压力的低值区，极易在数值结果中出现负值，导致计算失败。文献[46]在P4-T2-BVD格式中引入了MOOD(Multi dimensional Optimal Order Detection)正定修正算法[69]，从而保证了热力学变量数值解的正定性。图9显示二维爆轰波90°凸角绕流的密度和温度分布，对于反应区和各种流场结构均获得较精确的计算结果。

(a) P4T2-BVD

3.4.3 空气激波与水柱相互作用[70]

采用THINC函数作为候补重构函数的BVD格式可以有效控制数值耗散，从而避免接触间断和物质界面在计算过程中被平滑扩散的现象。文献[36]采用MUSCL-THINC-BVD格式对各类两项可压缩流动进行了模拟，取得了良好的计算结果。

图10给出了空气与水柱相互作用的计算结果。BVD格式能有效地保持水/气界面的厚度，并对各种流场结构都能得到高保真的计算结果。与文献[70]及其他现有算法相比，BVD格式能在抑制数值振荡的同时，获得间断解的高保真计算结果，显示出在模拟可压缩多相介质等具有显著间断解结构特征应用问题时的明显优势和潜力。文献[43]将基于MUSCL和THINC的BVD格式扩展到非结构网格，并计算了包含移动界面的可压缩多相流算例，获得了良好的数值结果。

(a) 密度

(b) 温度

图10 空气与水柱相互作用的两相流算例在t=2.15时刻MUSCL-THINC-BVD格式的计算结果(密度梯度), 网格单元数为2000×500.

4 结 论

本文概述了BVD方法的基本思想，设计相关格式的基本思路，以及一些具有很强实用价值的BVD格式。并通过单相和两相可压缩流动的一些典型算例验证了BVD格式的特点和优势。

BVD方法把减小数值耗散作为设计空间重构算法的指导思想，具有明确的物理含义和普遍性，是一条设计和改进数值格式的新途径。针对研究对象，选用合适的BVD许容函数作为空间重构的候补选项，并设计相应的BVD算法可以构建一类全新的数值格式。这些格式与基于多项式重构和非线性限制修正的传统方法有本质区别。数值实验结果表明，BVD格式可以得到明显优于已有算法的计算结果，能同时对各种尺度的连续和间断结构都获得高保真的数值解。

采用THINC或类似的Sigmoid函数作为候补重构函数可以明显改善基于欧拉网格的数值方法对接触间断和物质界面在计算过程中被平滑扩散的问题，对解决含有化学反应面或移动物质界面流动问题的实际应用具有重要意义。

相关研究工作仍在继续，例如，可以通过BVD判据，在光滑区引入中心格式进一步控制数值耗散[71]，以及将BVD格式扩展到间断伽列金的数值框架[72]。希望通过不断努力和持续发展，能够提出具有更好性能的BVD格式，并成为解决实际应用问题的有力工具。

致谢:邓希、谢彬、孙紫尧、程李东、金鹏、陈春刚、Shimizu、Inaba、Tann、Wakimura等东京工业大学工学院肖研究室的学生以及合作者姜振华、Loubère、Abe等对本项研究工作做出了重要贡献，特此致谢。