在包含运动边界的非定常流动模拟中，流场物理空间随时间变化，边界运动导致网格也随时间变化，典型应用包括多体分离、机翼抖振颤振、舱门启闭等^[1]。CFD计算一般采用运动嵌套网格、笛卡尔自适应网格和变形动网格等技术处理网格变化。在采用上述技术处理网格变化时，均涉及流场信息的插值问题^[1-2]，即利用旧网格上的流动信息插值获得重构区域流动参数。此外，在爆轰流场模拟中，爆轰计算网格与冲击波传播计算网格差两个量级^[3-4]，为提高计算效率，往往需要将物理量由爆轰计算时使用的密网格插值传递到冲击波传播模拟中的疏网格。

插值过程一般应具备守恒性和单调性等基本特性。对爆轰计算等守恒性要求高的问题，一般的插值不能保证计算区域内物理量的守恒，会降低计算精度，且无法捕捉激波的正确位置^[5]，需要发展守恒插值方法。守恒插值方法主要可分为两类，基于面通量(SFB/DC，simplified face-based donor-cell)的方法和基于单元相交(CIB/DC, cell-intersection-based donor-cell)的方法^[6]。

SFB/DC方法^[6]效率较高，但要求新、旧网格具有相同的拓扑结构，不适用于网格重构等拓扑发生变化的情况。CIB/DC方法思想简单直接，对计算网格的类型、结构不进行任何假设，适用范围很广，并且具有保号性，插值过程中不会产生新的极值^[6]。但该方法推广至三维情况时，由于三维网格单元重叠区域切割计算非常复杂，给这种方法的应用带来了很大困难。为了实现严格意义下的三维CIB/DC算法，P.E.Farrell等^[7]进行了一系列卓有成效的工作，在2011年进一步发展为“局部超网格”(Local Supermesh)方法，成功实现了非结构网格下的三维CIB/DC算法^[8]。

本文中在P.E.Farrell等和S.Menon等^[9]工作的基础上，构造了应用于二维/三维混合网格的CIB/DC算法，实现了基于精确网格切割的积分守恒插值。验证结果表明，该方法能够保证插值过程中计算域内物理量的严格守恒，具有比常规二阶插值更高的精度。

1.   守恒插值基本思想

守恒插值要求在插值过程中保证守恒量的积分值保持不变，假设计算域为Ω，网格T^a是计算域的一种划分，则在网格T^a中，Ω被划分为N_a个网格单元T^a=(T₁^a, T₁^a, …, T_Na^a)，并满足如下条件：

在网格T_i^a中，物理量ϕ的分布函数为：

式中：ϕ_i^a(x)可以指定多种形式，通常为多项式分布。

网格T^b是计算域Ω的另一种形式的划分，划分形式与要求同T^a，则同样有T^b上的物理量分布ϕ^b(x)。如果对于Ω的任意子区域D，满足如下条件：

则称ϕ^b(x)是ϕ^a(x)的守恒插值结果。

对于一般的情形，要使式(1)满足，必须有T^b≡T^a，ϕ^b(x)≡ϕ^a(x)，但这就失去了在两套网格间进行插值的意义。为此，对D进行限制，不能是Ω的任意子区域，而必须是网格T^b中的某个网格：

上式就是一般意义上的守恒插值关系式。

对式(2)，按照基于单元相交的思想，对D进行分解：

式中：V_ij=T_i^b∩T_j^a。将上式代入式(2)，则有：

式(3)就是基于单元相交守恒插值方法(CIB/DC)的基本依据。从式(2)到式(3)未引入任何近似。因此，在给定分布ϕ(x)的情况下，CIB/DC方法是严格守恒的。

从式(3)还可以看出，如果已知网格T^a上的物理量分布ϕ^a(x)，如何计算网格单元间的相交区域V_ij便成为CIB/DC方法的关键。

2.   “超网格”基本思想

图 1为文献[8]给出的超网格示意图，超网格是由T^a和T^b所有结点以及它们边界交点构成的网格，具有如下特性：

图  1  超网格示意图

Figure  1.  Schematic of the supermesh

下载: 全尺寸图片 幻灯片

(1) 每一个超网格不再和T^a或T^b中任何网格单元相交，即它总可以在T^a中找到一个网格包含或等于它，也总可以在T^b中找到一个网格包含或等于它；

(2) T^a的每一个网格总可以明确表示为一个或多个超网格的集合，T^b的每一个网格总可以明确表示为一个或多个超网格的集合。

这样一来，对于T^b任意网格D可以写出它所包含的超网格，每一个超网格就是V_ij=T_i^b∩T_j^a，可以在T^a上找到具体位置和对应的物理量，采用ϕ^b(x)≡ϕ^a(x)，按照D所包含的超网格进行求和，就得到D的物理量，从而实现了从T^a到T^b的信息传递。由于T^a和T^b在超网格的交集完全等价，使传递的物理量严格守恒，实现精确守恒插值。

3.   单元相交算法

由上节可知，获得超网格是实现精确守恒插值的关键。下面给出计算平面凸多边形和三维凸多面体交集的具体算法。

3.1   平面凸多边形相交算法

如图 2所示，将目标三角形A₁A₂A₃视为切割多边形，使用切割三角形B₁B₂B₃对目标三角形进行切割，每次切割时只保留目标多边形中与切割多边形在切割线同侧的部分。图中依次使用△B₁B₂B₃的3条边B₁B₂、B₂B₃、B₃B₁切割△A₁A₂A₃，切割后保留的多边形依次为CA₃A₂G、CDEA₂G、CDEFB₁，多边形CDEFB₁即为△A₁A₂A₃和△B₁B₂B₃的交集。

图  2  二维网格切割示意图

Figure  2.  Schematic of the two-dimensional grid cutting

下载: 全尺寸图片 幻灯片

在使用直线去切割平面凸多边形时，以图 3为例，切割直线经过点O，其法向单位矢量为n，切割直线将平面划分为两部分，定义法向n所指的半平面为正半平面，反方向的半平面为负半平面。如果切割之后目标多边形仅保留正半平面的部分，那么切割操作等效于求目标多边形和正半平面的交集。定义平面上任意一点P到切割直线的距离为：

图  3  直线切割平面凸多边形

Figure  3.  Schematic of a straight line cuttingplane convex polygon

下载: 全尺寸图片 幻灯片

式中：x_P和x_P分别为点P和点O坐标。

根据上述定义，正半平面的点到切割直线的距离均为正数，负半平面的点到切割直线的距离均为负数，切割直线上的点到切割直线的距离为零。如果凸多边形各顶点到切割直线的距离均为非负数(图 3(a))，则切割结果即为原多边形；如果凸多边形各顶点到切割直线的距离均为非正数，则切割结果为空集；除上述两种情况外，凸多边形与切割直线必然有两个交点(图 3(b))。具体如下：

(1) 设目标多边形为N，其顶点依次为N₁, N₂, …, N_N；切割直线为l，切割结果为R；

(2) 计算N的各个顶点到l的距离，记顶点N_i到l的距离为d_i；

(3) 初始化一个空的有序点集P；

(4) 从顶点N₁起，依次对多边形N的边N₁N₂, N₂N₃, …, N_NN₁进行判定。

对线段N_iN_j的判定过程为：如果d_id_j＜0，则线段与切割直线有新的交点，记交点为N_new；依次将N_i、N_new、N_j中距离为非负数的顶点加入点集P；如果d_id_j≥0，则依次将N_i、N_j中距离为非负数的顶点加入点集P；

(5) 完成所有边的判定后，如果P为空集，则R为空集；若P为非空集，则依次连接P中的点形成的多边形就是所求的切割结果R。

3.2   三维凸多面体相交算法

在三维流体计算中，常用的网格包括四面体、四棱锥、三棱柱、六面体等多种类型。其中，四棱锥、三棱柱、六面体这3种网格都具有四边形表面，网格生成过程中不能严格保证四边形表面的四个顶点位于同一平面内，会给多面体表面的定义带来歧义。

为了简化算法的复杂性，同时消除歧义，在计算两个任意凸多面体的交集时，首先将凸多面体分解为多个小四面体，如图 4所示，其中点N₆为点N₂、N₃、N₄、N₅的重心；然后通过计算这些小四面体的交集，来获得多面体的交集。计算多面体P_a和P_b的交集的具体步骤为：

图  4  三维网格分解示意图

Figure  4.  Schematic of the three dimensional mesh decomposition

下载: 全尺寸图片 幻灯片

(1) 将P_a分解为一组四面体{T_aⁱ}，将P_b分解为一组四面体{T_b^j}，并且满足：

(2) 计算{T_aⁱ}和{T_bⁱ}中任意两个四面体的交集：

并记录交集T_c^ij的体积和重心分别为V_c^ij，x_c^ij。

(3) 多面体P_a和P_b交集可表示为：

交集的体积和重心分别是：

该方法的核心在于如何计算两个四面体的交集T_c^ij。相对于平面多边形的情形，即使四面体是最简单的多面体，其交集的计算也是十分复杂的：两个四面体相交，最多可能产生一个12顶点的8面体。本文中使用切割法计算两个四面体的交集：将一个四面体看做目标四面体，另一个四面体为切割四面体，依次使用切割四面体的4个表面去切割目标四面体，最终留下的多面体就是这两个四面体的交集。

4.   相交单元的搜索算法

相对于通常的流场插值方法，守恒插值的计算效率较低，主要有两个原因：(1)计算两个单元交集的算法较为复杂；(2)对同一个新网格单元，普通插值方法只需用到一个旧网格单元(称为贡献单元)的信息，而守恒插值需要用到多个贡献单元的信息。前文已对单元相交算法进行了详细的描述，本节将介绍守恒插值中多个贡献单元的搜索算法。

在贡献单元的搜索算法中，首先需要快速定位出新网格单元在旧网格中的位置。对此，可采用四叉树^[11]或八叉树^[12](三维)方法加速搜寻，如图 5(a)所示，完全填充网格为旧网格，其上的非填充网格为新网格中的某个单元，通过四叉树方法，可快速定位出新网格单元的格心所在的旧网格单元(图中填充单元)，这是新网格单元的第1个贡献单元；然后，逐个搜索各个贡献单元的直接相邻单元，如果相邻单元与新网格单元相交，则将该单元加入贡献单元的列表，如图 5(b)中浅色填充单元所示；当所有贡献单元的相邻单元(图 5(c)与(b)中不同的填充单元)都不与新网格单元相交时，搜索结束。

图  5  相交单元的搜索方法

Figure  5.  Search method of intersection elements

下载: 全尺寸图片 幻灯片

5.   守恒插值验证算例

设计了3个算例来验证本文发展的守恒插值方法。第1个算例主要考核相同拓扑结构下网格点位置变化的影响；第2个算例主要考核网格形状的影响；第3个算例主要考核三维情况下的应用情况。

算例 1

考虑如下的测试函数：

计算区域为[0, 1]×[0, 1]，计算网格为一组拓扑结构相同的非结构三角形网格，分别表示为M₀、M₁、…、M_N，其中M₀如图 6所示。网格M_N中的第(i, j)个网格点的坐标由下式给出：

图  6  数值计算建模图

Figure  6.  Geometry of the specimenused for simulation

下载: 全尺寸图片 幻灯片

式中：α(t)=0.5sin(4πt), ξ=(i-1)/(i_max-1), η=(j-1)/(j_max-1)t=n/T, 其中i=1, 2, …, i_max; j=1, 2, …, j_max; n=0, 1, …, N。

测试过程中，在网格M₀上给定测试函数的准确分布，然后依次插值到网格M₁、M₂、…、M_N上，共进行N次插值，最后比较网格M_N上的函数分布与准确值之间的误差。在本文计算中，取i_max=j_max=65，N=200，T=20。

图 7给出了计算误差随插值次数的变化情况，其中纵坐标代表计算区域内各单元函数分布插值后计算值与理论值之差的1范数。两种插值方法中，插值误差都随插值次数的增加而增大，但守恒型插值方法的误差始终较小，不到常规二阶插值误差的一半。

图  7  计算结果的误差比较

Figure  7.  Different errors of computational results

下载: 全尺寸图片 幻灯片

图 8为两种方法的物理量积分值比较，其中纵坐标代表每个单元计算值总和与理论值的比，从图中明显可以看出，守恒插值保证了计算域内物理量的守恒。

图  8  计算结果的积分比较

Figure  8.  Different integrals of computational results

下载: 全尺寸图片 幻灯片

算例2

考虑如下的测试函数：

计算区域是边长为1的正方形，L为正方形的对角线长度，r是与正方形中心的距离。测试中使用两种不同的网格，网格1是均匀分布的结构网格，网格2是非结构网格，并与网格1具有相同的边界网格点分布，如图 9所示。测试中，首先在网格1上给出测试函数的准确分布，然后将其插值到网格2上，再插值回网格1，如此反复，进行200次插值后，将物理量分布与初始的准确分布进行比较，统计插值过程中的误差。分别采用常规二阶插值和守恒型二阶插值进行测试，并使用不同疏密的网格以分析插值方法的精度。测试中使用的4组网格的网格量如表 1所示。

图  9  数值计算建模图

Figure  9.  Geometry of the specimen used for simulation

下载: 全尺寸图片 幻灯片

表  1  计算网格参数

Table  1.  Parameters of the computational grid

序号	结构网格	非结构网格
		节点数	单元数
1	33×33	695	1 260
2	65×65	2 593	4 928
3	129×129	10 044	19 574
4	257×257	39 687	78 348

下载: 导出CSV 

| 显示表格

图 10是计算误差与网格步长之间的关系，其中横坐标是表 1中4个不同结构网格计算单元的边长，纵坐标代表计算区域内各单元函数分布插值后计算值与理论值之差的1范数。可以看出在相同的网格步长下，守恒型插值的误差更小，而且误差与网格步长之间的变化斜率更大，说明其具有更高阶的精度。计算结果表明，常规二阶插值精度约为1.4阶，而守恒型二阶插值精度约为1.87阶。表 2是不同算例的守恒误差比较，守恒插值方法的守恒误差接近机器零，明显优于普通二阶插值。

图  10  插值方法的精度比较

Figure  10.  Comparison of precision with different interpolation methods

下载: 全尺寸图片 幻灯片

表  2  不同算例的守恒误差

Table  2.  Conservation errors indifferent numerical cases

序号	二阶插值	守恒型二阶插值
1	3.724×10-3	0
2	-1.447×10-3	0
3	-3.985×10-4	1×10-14
4	5.636×10-5	-2×10-14

下载: 导出CSV 

| 显示表格

算例3

考虑如下的测试函数：

计算区域是边长为1的正方体。每次测试使用两套网格，首先在网格1上给出测试函数的准确分布，然后将其插值到网格2上，再插值回网格1，如此反复，进行10次插值后，将物理量分布与初始的准确分布进行比较，统计插值过程中的误差。分别采用普通二阶插值和守恒型二阶插值进行测试，并使用不同疏密的网格以分析插值方法的精度。测试中使用的4组网格的网格量如表 3所示，其中网格步长由下式得到：

表  3  计算网格参数

Table  3.  Parameters of the computational grid

序号	网格1			网格2		dx
	节点数	单元数		节点数	单元数
1	914	4 033		676	2 711	0.067
2	2 222	10 632		1 424	6 075	0.050
3	6 527	32 913		3 733	16 681	0.035
4	26 099	138 542		12 600	59 347	0.022

下载: 导出CSV 

| 显示表格

式中：N_E1和N_E2分别为两套网格的单元个数，N_D=3为网格的维数。

图 11是测试函数在三维空间中的分布云图，图 12是计算误差与网格步长之间的关系，其中横坐标是表 3中dx的值，纵坐标代表计算区域内各单元函数分布插值后计算值与理论值之差的1范数。在三维条件下，守恒型插值依然具有更小的计算误差和更高的计算精度。

图  11  测试函数分布云图

Figure  11.  Nephogram of different test functions

下载: 全尺寸图片 幻灯片

图  12  插值方法的精度比较

Figure  12.  Comparison of precision with different interpolation methods

下载: 全尺寸图片 幻灯片

6.   结论

采用在新旧网格相交单元剖分构建出来的超级网格作为网格之间信息传递的基础，采用四叉树(二维)或八叉树(三维)方法搜寻算法提高计算效率，基于以上方法提出了网格之间流场物理量的守恒插值算法。该算法理论上可以实现网格之间物理量的精确守恒，二维和三维算例也表明产生的误差接近机器零，明显优于常用的二阶插值。本文中提出的方法可以解决爆轰流场模拟中网格尺度和物理量相差较大情况下的插值难题。