黎曼边界条件在高阶精度非结构有限体积方法中的验证与应用

2021-06-24孔令发董义道

空气动力学学报 2021年3期

孔令发，刘伟，董义道

(国防科技大学空天科学学院，长沙 410073)

0 引言

相比结构化网格和笛卡尔网格，非结构网格因其对复杂外形的适应性、网格生成的自动化程度[1-3]以及对运动边界、多体分离、流场自适应等问题的模拟能力，逐渐成为学术研究和应用的重点[4]。

基于非结构网格的数值格式通常包含有限元（Finite Element, FE）[5]、间断伽辽金（Discontinuous Galerkin, DG）[6-7]和有限体积（Finite Volume, FV）[8-13]等方法，其中，具有二阶精度的FV方法因其在鲁棒性和计算开销较小等方面的优势而被广泛应用于当前主流的CFD商业与开源软件。这也一定程度上极大促进了二阶FV格式的发展，使其趋于成熟。

在二阶精度的基础上，近年来具有高阶精度的FV离散格式同样得到了较快发展，以满足高效计算与高保真数值模拟的需求[12-15]。Ollivier-Gooch与Jalali等[16-18]系统研究了高阶非结构FV格式中的k-exact重构[19,20]、曲边界处理与通量、源项积分等一系列问题。Li[21]在高阶非结构数值格式的框架内探索了新的限制器构造策略，以及考虑黏性边界条件的约束k-exact重构。Kong等基于k-exact方法，在二阶与高阶FV离散格式的框架内提出了一种新的重构模板[22]与单元参考点选择策略[23]，有效改善了原有格式在各向异性网格上的数值表现。

相比重构与离散算法研究，在非结构FV算法，尤其是高阶精度算法框架内，针对边界条件的细致研究相对较少，然而实践表明，采用不同的边界条件设置方法可导致流场求解效果的差异[24]。如果边界条件设置恰当，求解过程的鲁棒性以及计算效果能够得到明显改善；反之，若边界条件设置不当，则有可能收敛至错误的流场，甚至导致计算发散。

对于真实流动，通常可将边界条件的施加分为强施加与弱施加两类。强施加是一种变量施加方法，即通过边界条件得到边界位置的变量；而弱施加是一种通量施加方法，即通过边界条件得到边界位置的数值通量。强施加方法一般基于特征理论，认为双曲守恒律方程代表特征波的传播，在对应的边界位置存在传入或者传出计算域的特征波[25-26]。其中，外传波的特征一般取决于内部流场；而内传波需要基于无黏关系式，在界面构造演化方程[27]。

对基于非结构网格的数值格式，当前较为常用的强施加边界条件是引入虚拟单元来设置特征边界条件[28-29]，该方法能够准确反映边界附近特征波传播情况，但此过程需要人为判断特征波传播方向以得到时变幅值，对于多维效应明显的流动，通常还需要引入横向修正[29]。此外，对边界面处特征虚拟单元的变量赋值过程引入了一阶误差，严重破坏了流场求解精度[30-31]。而本文的研究目标是探讨边界条件在高阶精度有限体积方法中的表现效果，因此不再考虑这种强施加边界条件在高阶精度求解器上的应用。

不同于强施加方法，弱施加边界条件旨在得到边界面的数值通量。对于内部单元面，可根据单元面处的变量左右值采用近似黎曼解构造数值通量，而对于边界面的通量构造同样可采用此方法。Kopriva等[32-33]指出，FUN3D的边界条件处理也采用了这种方法。但当前对亚声速入口、出口以及远场这类边界条件的主要处理方式是基于一维黎曼不变量来得到边界面处的右状态矢量，这种方法仍需采用特征关系式，并基于黎曼不变量进行繁冗推导。为简化这类弱施加边界条件的处理方法，Atkins和Casper[34]针对结构网格上的高精度格式，提出了一种基于有限波模型的精度保持边界条件施加方法来统一处理包括亚声速入口、出口以及远场等边界，该方法基于有限波模型得到的精确解，将内部流场与外部瞬态条件有效联系起来，得到了高阶收敛精度。Dong等[30-31]将这种边界条件由基于结构网格的有限差分方法推广至二阶精度非结构有限体积方法，并将其与常用的特征虚拟单元边界条件对比后发现，采用黎曼边界条件能够得到较高的计算精度与收敛效率，同时也可显著提高求解器的计算稳定性。

在此基础上，本文将这种弱施加边界条件推广至高阶精度非结构有限体积方法，以检验其在高阶精度求解器中的表现效果。此外，在高阶有限体积方法的框架内将这种边界条件推广至基于制造解方法（MMS）[35-38]的流动。这种特殊的数值流动当前被广泛用来检验算法的有效性，但通常在基于MMS流动研究某些问题时，一般未考虑对其施加边界条件。较常见的处理方式是所有边界面处的变量左值由内部流场插值得到，而右值均赋值为当前位置的精确解。Folkner和Katz[35]以及Wang等[36]在二阶精度非结构有限体积方法的框架内验证了施加边界条件的MMS流动的数值表现，经检验，对MMS流动施加边界条件后可达到二阶精度。但在高精度数值方法中，对此类流动施加边界条件，尤其是黎曼边界条件的相关研究还未开展。因此，本文计划在高阶精度非结构有限体积方法的框架内，验证对MMS流动和部分真实流动施加黎曼边界条件的合理性，并对比其与常用无反射边界条件[35]间的差异，以进一步推广此类弱施加边界条件在非结构网格数值格式中的应用。

1 非结构有限体积格式

1.1 控制方程

对无黏流动探讨弱施加黎曼边界条件，需要求解的控制方程为：

其中，u为待求解变量，Fc为无黏通量。该方程可在散度定理的基础上，通过FV格式离散为：

式中Nf与NG分别为网格面与通量点的数量， Φjk为通量点处的数值通量函数，本文采用Roe格式[39]来计算该函数。 ωk为对应通量点的权系数[23]，njSj为对应单元面的面积矢量。

对于高阶方法，一个通量点无法满足积分精度要求，因此需要多个通量点（如图1所示）。对应不同精度下的通量点的数目可参考文献[23]。在基于Roe格式构造通量点数值通量时，需要知道该点的变量左右值，并且变量左右值的获取依赖于从单元参考点到面通量点的泰勒展开，因此需重构控制体单元的变量梯度及其高阶导数。

图1 数值通量构造示意图Fig. 1 Schematic diagram of the numerical flux evaluation

1.2 k-exact重构算法

对于二阶格式，常采用最小二乘（Least Square, LSQ）法与格林高斯（Green Gauss, GG）方法来获得变量梯度。但对基于积分型控制方程的高阶FV格式而言，传统的LSQ或GG方法均无法满足高于二阶的重构精度要求，因此本文采用k-exact方法[16-18]计算变量的梯度与高阶导数。该方法的详细过程本文不再重述，最终需要求解的方程组形式如下：

同时，权系数 ωij为参考点间距离的倒数，以强调与中心控制体相邻单元对重构过程的影响。

2 弱施加边界条件

由于强施加方法多依赖于特征虚拟单元来处理相应的边界条件，该方法在对无黏边界处的虚拟单元赋值过程引入了一阶误差。因此，本文不再讨论这种边界条件强施加方式。

为了在边界面施加弱边界条件，根据近似黎曼公式，需要给定边界面处的两个状态，如图2所示，其中左状态由内部流场插值得到，右状态依赖于相应的边界条件。

图2 重构模板与边界条件弱施加示意图Fig. 2 Schematic diagram for the reconstruction stencil and weak-imposition boundary condition

在基于内部流场插值时需要注意避免引入一阶误差。本文研究具有三阶精度的有限体积方法，因此首先可通过内部流场信息得到边界面左值，即：

式中WC为边界单元参考点处的变量值，b为边界面通量点。因此，通过方程（5）可得到边界面高斯积分点处的左状态矢量。边界面右值需要根据相应的边界条件获得。下面将针对亚声速入口、出口以及远场等不同边界，给出其边界条件的具体施加方式。

2.1 无黏壁面边界

对于无黏壁面，不同弱施加方法的处理过程基本一致。在施加边界条件时，需要保证速度的左右状态关于壁面边界对称。为满足这个条件，右状态矢量可定义为：

其中逆变速度由左值计算得到：

压力和温度满足零梯度条件，即：

2.2 超声速入口与出口边界

由于在超声速条件下特征波传播方向唯一，因此超声速入口、出口的处理相对简单。其中超声速入口右值直接取来流值，即：

超声速出口右值与左值相等：

2.3 亚声速入口与出口边界

与壁面、超声速入口和出口的处理方式不同，在处理亚声速入口、出口以及远场时，本文采用黎曼边界条件。该方法通过求解对应的黎曼问题，在边界面位置可得到物理上相容的左右状态，因此避免了采用特征关系式与黎曼不变量对这类边界条件的判断与推导[35]。

如图3所示，有限波模型假设内部流场与外部瞬态状态通过两道由接触间断分隔的膨胀波进行过渡，图中标有u的虚线代表接触间断，另外两条红色直线对应膨胀波，膨胀波之间的区域将其称为星区。

现构造如下黎曼问题，无穷远处的自由来流参数作为右状态矢量，记为WR，而左状态矢量WL由内部流场插值得到。通过精确求解此黎曼问题，可以获得星区的左右状态矢量，分别记为与。星区的压力可通过如下关系式获得：

其中，c为局部声速，γ为比热比。逆变速度的定义为：un=u·n，其中u为笛卡尔坐标系下的速度矢量，n为边界面处的单位法向矢量。在得到星区压力的基础上，可采用等熵关系式得到星区的密度：

以及星区的逆变速度：

在此基础上，边界右状态可通过星区状态矢量即给定。

图3 有限波模型示意图Fig. 3 Schematic diagram of the finite wave model

相比常用的无反射边界条件[35]，采用黎曼边界条件处理亚声速入口、出口与远场边界时，均可通过精确求解对应的黎曼问题来给定边界右值，无需再基于法向速度与特征关系来判断不同边界，并单独处理。因此，这种边界条件施加方式可有效简化对亚声速入口、出口以及远场边界的处理过程。

3 MMS流动中的精度验证

当前，MMS方法广泛应用于对算法计算精度的检验，但在基于该方法验证精度时，一般未对其施加任何物理边界条件，即边界面的右状态均以参考点处的精确解赋值。本节将对MMS流动采用黎曼方法处理亚声速入口、出口边界，以验证对MMS流动施加黎曼边界条件的合理性。同时为了验证工作的完整性，还测试了对MMS流动施加无黏壁面以及超声速入口、出口这类具有真实物理意义边界条件后的数值表现，以期在高阶精度非结构有限体积框架内，为MMS流动提供一种便捷合理的边界处理方式。

测试网格分别采用了六套由疏到密的规则三角形网格与扰动三角形网格。网格示意图与网格点分布情况见图4与表1。

图4 规则与扰动三角形网格示意图Fig. 4 Schematics of regular and randomly perturbed grids

表1 背景四边形网格在x与y方向的分布量Table 1 Distribution of background quadrilateral grid in x- and y-directions

当来流速度为亚声速时，计算域的左右边界分别设置为亚声速入口与出口；对于超声速流动，左右边界分别对应超声速入口和出口。此外，设置壁面边界时，考虑到特征波的传播方向，统一将计算域上边界设置为壁面，示意图见图5。

图5 边界条件设置示意图Fig. 5 Schematic diagram of the corresponding boundary conditions

3.1 无边界条件

在研究不同边界条件的影响之前，为了验证求解器的求解精度，首先测试了不添加边界条件的MMS流动计算效果。在相应的边界高斯积分点处，左值均由内部流场插值得到，而右值由MMS制造解给定。其形式如下：

流动控制方程如式（1），因此可首先根据控制方程推导出源项s的统一表达式为：

基于该解析解，得到的流场如图6所示。

计算收敛后，统计了每套网格误差的L1、L2以及L∞范数来充分测试高阶求解器的计算精度。图7中Slope3表示三阶精度标准线。从图7中可以看出，不添加边界条件时，在疏网格上精度一开始并未达到三阶，但随着网格加密，计算结果逐渐趋于设计精度。因此，在不施加边界条件的情况下，首先验证了求解器具有三阶精度。接下来对计算域的不同位置单独施加黎曼边界条件，来验证在高阶方法中对此类流动施加边界条件的合理性。

图6 当前制造解流场Fig. 6 Flow fields of current manufactured solutions

图7 无边界条件的误差曲线Fig. 7 Error curves without any boundary condition

3.2 亚声速入口、出口边界

对于真实流动的亚声速入口与出口，可通过精确求解黎曼问题来获得对应的边界条件。在求解黎曼问题时，WR通常用无穷远处的来流参数赋值。而MMS流动与真实流动不同，一般不包含无穷远处的参数信息。因此尝试在MMS流动中，WR统一用边界处的制造解提供，并在得到右状态的基础上，通过求解黎曼问题来得到星区状态矢量最终由该矢量来获得亚声速入口、出口处的右状态参数。

验证这两种边界条件时，需要保证在对应的入口与出口位置速度为亚声速，因此给定解的形式如下：

对应该制造解，得到的流场如图8所示。

图8 当前制造解流场Fig. 8 Flow fields of current manufactured solutions

对入口与出口条件单独验证，在验证其中一个边界条件时，其余边界右值均使用MMS精确解赋值，这样有利于研究不同边界对流动求解精度的直接影响。统计误差如图9与图10所示。

从图9与图10反映的结果可以看出，对MMS流动单独施加亚声速入口、出口边界条件时，不论是在规则网格还是扰动网格上，均达到了三阶精度。因此，对MMS流动的亚声速入口与出口施加黎曼边界条件不会影响流动的求解精度。

图9 施加亚声速入口边界条件的误差曲线Fig. 9 Error curves with the subsonic inlet boundary condition

图10 施加亚声速出口边界条件的误差曲线Fig. 10 Error curves with the subsonic outlet boundary condition

3.3 无黏壁面边界

不同的弱施加边界条件对无黏壁面的处理方式基本一致。壁面法向速度为0，压力和温度的右值与左值相等。为保证速度的法向分量为0，验证壁面时给定解的形式为：

Folkner和Katz等[35]在二阶精度有限体积方法中验证了对MMS流动施加此类含有真实物理意义边界条件的可行性。本节将其推广至高阶精度非结构有限体积方法，以检验在高阶方法中对MMS流动施加此类边界的合理性。基于方程（27）得到的流场如图11所示。

图11 MMS流场信息Fig. 11 Flow fields of the MMS flow

在验证壁面边界时，除壁面外其余边界右值均由MMS精确解赋值。在规则与扰动网格下的统计误差如图12所示。由图12可以看出，对MMS流动施加相应的壁面边界条件后，不论在规则网格还是扰动网格上，计算结果仍可达到三阶精度。由此验证了对MMS流动施加壁面边界条件的可行性，同时也将文献[35]中提出的壁面边界验证方式有效推广至高阶精度非结构有限体积算法。

图12 施加无黏壁面边界条件的误差曲线Fig. 12 Error curves with the inviscid wall boundary condition

3.4 超声速入口、出口边界

相比亚声速入口、出口的边界条件施加方式，超声速入口与出口的边界条件施加相对简便—入口采用MMS精确解赋值，出口的右值等于左值。为保证入口、出口处为超声速，给定解的形式为：

图13为该制造解对应的流场云图。与3.2和3.3节的验证方式相同，对入口与出口条件单独验证。图14与15分别给出了施加超声速入口与出口条件得到的误差。从图14与图15反映的结果来看，单独对MMS流动施加超声速入口或出口条件，并未破坏计算精度，L1、L2与L∞误差均达到了三阶精度。我们将三类边界条件在规则密网格上得到的计算精度汇总，结果如表2所示。

从表2可以看出，对MMS流动的不同边界单独施加相应的边界条件后，在密网格上L1与L2误差均达到了三阶精度，L∞误差略低于但接近三阶精度。此外，从误差统计曲线可以明显看出，对MMS流动施加合理的边界条件后，随着网格加密，求解精度的变化较小，即使在稀疏网格上结果同样接近三阶精度，并趋于一致收敛。因此，这也从一个层面论证了在高阶方法中对MMS流动施加具有真实物理意义边界条件的合理性。

综上，黎曼边界条件在MMS流动中具有较好的数值表现，尤其在处理亚声速入口、出口时，为MMS流动的高精度数值模拟提供了一种简单有效的边界条件施加方式。

图13 MMS流场信息Fig. 13 Flow fields of the MMS flow

图14 施加超声速入口边界条件的误差曲线Fig. 14 Error curves with the supersonic inlet boundary condition

图15 施加超声速出口边界条件的误差曲线Fig. 15 Error curves with the supersonic outlet boundary condition

表2 不同弱施加边界条件在密网格之间的计算精度Table 2 Accuracy of different weak boundary conditions between vfin and vvfin grids

4 亚声速无黏圆柱绕流验证

为了测试上述边界条件在真实流动中的表现，在亚声速无黏圆柱绕流上检验其数值效果。该流动是典型的亚声速算例，边界条件包括亚声速远场和无黏壁面。由之前的分析可以知，黎曼边界可有效简化对亚声速远场的处理过程。为对比黎曼边界条件在高阶求解器中的数值表现，在处理远场边界时，还测试了常用的基于一维黎曼不变量的无反射边界条件[35]。网格采用O型拓扑，壁面附近的局部放大图如图16所示。主要来流参数Ma=0.38，远场取20倍圆柱直径。网格生成可参照文献[38]。计算在三套网格上进行，密网格径向与周向分别分布32与128个网格点。网格径向分布情况为：

其中α=1.1，圆柱半径r0=0.5。对最密网格疏化可得到剩余两套网格。

图16 圆柱绕流计算网格示意图Fig. 16 Schematic diagram of grids for the cylinder flow simulation

基于两种不同的弱施加边界条件，首先测试了由两种边界条件得到的流场压力云图，如图17所示。从流场云图可以看出，由两种边界条件得到的流场均具有较好的对称性。

图17 不同边界条件下的圆柱绕流流场Fig. 17 Flow fields of the cylinder flow with different boundary conditions

为进一步量化对比两种边界条件之间的差异，统计了沿圆柱表面的总压恢复系数与熵误差。总压恢复系数的定义如下：

根据等熵关系式可得到总压pt和静压p之间的关系，

由于该算例为无黏亚声速流动，因此理想情况下总压恢复系数应当精确等于1.0。图18给出了总压恢复系数沿壁面的分布情况。从图中可明显看出，在梳网格上，总压损失较大，而随着网格的不断加密，数值耗散逐渐降低，总压恢复系数接近理想值。

图18 圆柱绕流总压恢复系数分布曲线Fig. 18 Total pressure recovery coefficients of the cylinder flow

此外，统计了由不同边界条件得到的熵误差，结果如图19所示。从图19可以看出，采用两种远场边界条件得到的计算精度基本一致，但由黎曼边界得到的误差值略低于常用的无反射边界，更利于计算结果的准确性。综上，在亚声速圆柱绕流中，黎曼边界条件的数值表现优于无反射边界条件。

图19 熵误差统计Fig. 19 Entropy errors

5 Gaussian Pulse非定常熵波验证

本节引入添加初始高斯脉冲（Gaussian Pulse）的非定常熵波传播算例，来测试黎曼边界条件在出口处的无反射特性。初始扰动函数如下[40]：

式中xc取0.5，计算域为x×y∈[0,1]×[0,1]。并且除密度外，其他变量的初始化[31]如下：

计算过程中，分别对比了由黎曼边界条件与无反射边界条件得到的出口特性。在施加相应的边界条件时，为防止y方向边界条件对计算结果的影响，计算域的上下边界均采用外推边界。具体的边界条件施加示意图如图20。从图20可以看出，黎曼边界条件针对亚声速入口与出口可通过精确求解对应的黎曼问题来统一处理，相比无反射边界条件更加简便。

图20 两种边界条件的施加示意图Fig. 20 Schematic diagram for the imposition of two boundary conditions

为对比两种边界在出口处的无反射特性，分别设置了四套由疏到密的均匀四边形网格，分布量从60×60至160×160（图21）。在四套网格下，统计了添加初始高斯扰动的ρ沿x方向的变化情况，同时参考文献[40]，分别给出了时间t= 0.9 s、1.1 s以及1.5 s时的结果。统计结果如图22所示。表3、表4为不同时刻两种边界条件得到的L2误差的具体数值。

图21 疏网格与密网格Fig. 21 Coarse and fine grids

图22 不同网格上密度在三个时刻的统计结果Fig. 22 Density on different grids at three moments

表3 无反射边界条件在不同时刻的L2误差统计Table 3 L2 errors of the non-reflective boundary condition at different moments

表4 黎曼边界条件在不同时刻的L2误差统计Table 4 L2 errors of the Riemann boundary condition at different moments

从图22中可以看出，随着网格的不断加密，两种边界条件的密度变化曲线接近，在出口边界均具有较好的无反射特性。因此说明了对亚声速非定常扰流施加黎曼边界条件具有一定的合理性。此外，结合表3的数据可以看出，黎曼边界条件与无反射边界条件在四套网格上，不同时刻的L2误差接近，并在t=0.9 s与t=1.5 s时误差均低于无反射边界条件。

最后，相比无反射边界条件，黎曼边界条件针对亚声速入口与出口，可通过精确求解黎曼问题来统一处理，有效简化了边界条件的处理过程。