崔鹏程，唐静，李彬，马明生，邓有奇

中国空气动力研究与发展中心 计算空气动力研究所， 绵阳 621000

在多体分离、舱门开启、返回舱降落与螺旋桨滑流等非定常流动模拟中，常采用重叠网格方法，将复杂的流动区域划分为多个简单的子区域，各子区域内的网格独立生成，数值计算在各个子网格上独立进行，各子网格之间存在重叠关系，流场信息通过插值在重叠区域进行匹配和耦合[1]。重叠网格方法降低了网格生成难度，在多体相对运动数值计算方面具有很多优势，但各子网格间的数据插值很难保证守恒性[2]。

重叠网格间高精度插值是保证流场计算正确的重要基础，目前主要有两种插值方法：守恒插值和非守恒插值。非守恒插值方法实现简单，处理无间断流动时插值精度较高，对此国内外很多学者开展了大量的研究工作。田书玲等[3-4]对物理空间的曲面单元进行了标准化处理，实现了基于非结构动态网格的Lagrange插值方法；Marstin和Mcconnaughey[5]利用有限元思想构造空间标准单元的线性插值函数，研究了双线性、三线性重叠网格插值方法；康忠良等[6]基于线性重构思想，发展了一种适用于任意单元的重叠网格插值策略，黄宇等[7]按单元尺度区分和扩充贡献单元模板，改善了网格尺寸相差较大时的插值精度；周乃春等[8]采用逆向距离权方法对重叠网格边界的物理量进行插值。

但非守恒插值方法对网格的相对尺寸依赖性较强，在网格尺寸相差较大的区域插值精度较低，且不能实现物理量在插值区域的守恒性，在非定常数值计算中可能会引起非物理扰动，此外，非守恒方法即使计算网格质量很好，也会使速度的计算值在重叠网格边界之间产生错位[1,9]。

守恒插值方法精度较高，能实现物理量在插值区域的守恒性，对网格相对尺寸依赖性较小，但过程处理复杂，实现相对困难。常见的守恒插值方法有嵌入网格方法[10]、分段交界面方法[11]、ENO(Essentially Non Oscillatory)守恒重映方法[12]和总体守恒隐式插值方法[13]等。其中嵌入网格方法和分段交界面方法都是在几何上将不重合边界转化为重合边界，在三维情况下较难实现，ENO守恒重映方法需分解几何单元，构造几何单元之间的面相交问题比较复杂，总体守恒隐式插值方法虽然计算精度较高，但需要求解附加方程组，给并行计算增添了难度。

Farrell和Maddison[14]发展了一种新的插值方法，提出“超网格”思想，构造了二维网格的积分守恒插值，2011年文献[15]进一步发展了三维多面体的切割算法。徐春光等[16]发展了二维/三维混合网格下的单元相交算法，实现了二阶精度的守恒插值方法。

本文在Farrell[14]、Menon[15]和徐春光[16]等的基础上，基于“超网格”思想与单元相交算法，发展了一种适用于非定常数值计算的隐式并行混合重叠网格守恒插值方法，在每一个非定常时间步内构造局部超网格，利用局部超网格实现背景网格向目标网格的守恒插值。经过数值试验验证，对于二阶分布的流场变量，该方法在重叠网格区域可以保证严格守恒插值，加快残差收敛速度，提高流场变量的插值精度。

1 数值方法

1.1 控制方程

雷诺平均Navier-Stokes方程守恒形式的积分方程为

(1)

式中:W=[ρρuρvρwρE]T为守恒流场变量,ρ为密度,u、v、w为速度分量，E为能量；V为体积分;S为面积分；Ω为重叠区域;Hc为对流通量；Hv为黏性通量；n为流动区域边界上指向外侧的单位法向量。

本文使用的流场解算器为自主研发的MFlow软件[17]，对流通量和黏性通量的具体内容可以参考文献[18]，控制方程的离散采用格心格式的有限体积法，时间项采用LU-SGS(Lower-Upper Symmetric Gauss-Seidel)迭代[19]，湍流模型为一方程Spalart-Allmaras(S-A)模型[20]。黏性项采用中心格式，对流项采用Roe格式[21]，非定常计算过程采用多重网格加速收敛[22]。

1.2 积分型守恒插值条件

本文以两套重叠子网格为例，简要说明重叠网格守恒插值方法。在守恒插值过程中要求守恒量的积分值保持恒定。

重叠子网格为ΓA、ΓB，重叠区域为Ω，设ΓA为背景网格，ΓB为插值网格，qA(x)为守恒变量在ΓA的分布函数，qB(x)为守恒变量在ΓB的分布函数，ΠAB为ΓA插值到ΓB的插值算子，如果满足方程：

(2)

则称qB(x)是qA(x)的守恒插值，其中qB(x)=qA(x)ΠAB。

2 超网格构造方法

2.1 超网格定义

“超网格”概念最早由Farrell和Maddison提出[14]。对于任意两套网格ΓA、ΓB，重叠区域为Ω，NA、NB分别为其网格节点，EA、EB分别为网格边，K∈Γ为网格单元。由{ΓA,ΓB}组成的超网格ΓC满足以下条件：

1)NC⊇NA∪NB。

2)V(KC∩K)∈{0,V(KC)}∀KC∈ΓC,

K∈Γ,Γ∈{ΓA,ΓB}。

即在ΓA、ΓB重叠区域内，超网格的网格节点必须包含ΓA、ΓB所有的节点，且超网格ΓC的每一个网格单元KC与ΓA、ΓB网格单元的相交体体积为零或为KC的体积，超网格的每个网格单元总是包含于ΓA或ΓB的网格单元。

2.2 局部超网格构造

按照2.1节中超网格定义，任意两套网格ΓA、ΓB的所有网格节点及边可组成它们的超网格ΓC[15]。图1为文献[15]给出的超网格示意图，图中a～e为网格节点。将目标网格ΓB的每一个网格单元与背景网格ΓA的网格单元求交，求交得到的网格单元可构成ΓA、ΓB的超网格，母网格ΓA、ΓB的每一个网格单元总是可以拆分为一个或多个超网格单元的集合。

文献[15]给出了全局超网格的构造方法，即背景网格ΓA与目标网格ΓB完全重合，所得超网格ΓC也与ΓA、ΓB完全重合。在重叠网格中各子网格区域不可能完全重合，且在非定常数值模拟中各子网格位置会相对变化。本文提出一种适用于重叠网格非定常数值模拟的局部超网格构造方法，在每一个非定常时间步内构造局部超网格，利用局部超网格实现背景网格向目标网格的守恒插值。局部超网格构造方法如下：

1) 在重叠网格各子网格之间进行挖洞操作，标识出各子网内不需要计算的网格、计算网格和插值网格等。

2) 用ADT(Alternating Digital Tree)搜索方法求出子网格每一个插值单元的所有贡献单元[23]。普通的插值方法中一个插值单元只需用到一个贡献单元信息，但在基于超网格的守恒插值方法中，一个插值单元需用到多个贡献单元的信息。

构造局部超网格需搜索出所有与插值单元相交的贡献单元，图2给出了搜索贡献单元的示意图，首先用ADT搜索方法快速定位插值单元的体心在背景网格的位置，体心所在的单元为1个贡献单元，如图2(a)所示；然后以此贡献单元为基础，检索所有与此单元共点的背景网格，搜索出与插值网格相交的网格单元，这些单元也是贡献单元；重复检索与贡献单元共点的网格单元，当所有与贡献单元共点的网格都不与插值单元相交时，终止搜索，如图2(b)所示。

3) 求出重叠网格子网格的插值单元与其每一个贡献单元的相交多面体，相交多面体可构成局部超网格。

4) 以局部超网格为媒介，在每一步内迭代计算中实现贡献单元向插值单元的守恒插值。当前时间步计算结束时，释放超网格。下一个时间步开始计算时，根据当前子网格的相对位置重新建立超网格。

2.3 网格求交算法

由2.2节可知，插值单元与贡献单元求交是获取超网格的关键，不同的网格单元有不同的网格求交算法。对凸多边形或多面体而言，二维和三维网格求交方法已经较为成熟[24]，本文采用快捷的Sutherland-Hodgman 逐边切割方法[16,25]。

Sutherland-Hodgman方法的基本思想为用切割多边形的各边界线逐次对目标多边形进行切割，从而得到一个新的多边形。图3为Suther-land-Hodgman方法示意图，以两个二维多边形求交为例，三角形A1A2A3与B1B2B3求交。把A1A2A3当做目标三角形，B1B2B3当做切割三角形。用B1B2B3的各条边依次对A1A2A3进行切割，得到相交多边形B1CDEF。

三维情况下Sutherland-Hodgman方法的思想为用切割多面体的各个面逐次对目标多面体进行切割，从而形成一个新的多面体。在三维数值模拟中，常见的网格单元包括四面体、三棱柱、金字塔和六面体等单元。其中三棱柱、金字塔和六面体单元具有四边形表面。在用常见的Pointwise、Icem等商业软件生成网格的过程中，为了提高物面网格投影的精度和空间网格的光滑性，会对网格进行扭转或弯曲变形，从而导致四边形的4个顶点不能严格地分布在同一平面内，在多面体切割的过程中会产生歧义[26]。为了消除歧义，在多面体求交时，首先将三棱柱、金字塔和六面体单元分解为多个小四面体单元。通过小四面体相交体的组合来实现多面体求交。图4所示为四面体网格求交算法，四面体P1与P2求交，将P2当做目标四面体，P1当做切割四面体，用P1的各个面对P2进行切割，保留切割面内侧的交点与交线，从而获得P1、P2的相交体P3。相比平面多边形而言，四面体求交得到的多面体构造比较复杂，图4中四面体P1、P2的相交体P3是具有10个顶点的7面体。

图5所示为用网格求交方法构造三维超网格的方法，具体步骤如下：

1) 在获取重叠网格局部超网格时，先用ADT搜索方法求出子网格每一个插值单元的所有贡献单元。在三维情况下，一个插值单元可能有多个对应的贡献单元，图5(a)为一个插值单元与其相交的背景单元，这个插值单元有2个贡献单元。

(3)

(4)

则插值单元PC与贡献单元PD的相交体PI可表示为小四面体交集的并集：

(5)

其体积VI与体心XI分别为

(6)

3) 插值单元与贡献单元的相交体可组成局部超网格，图5(b)为一个插值单元与其超网格，每一个超网格单元都由插值单元与其贡献单元相交获得。

3 基于超网格的守恒插值方法

3.1 积分转换

本文使用的流场计算器MFlow软件控制方程的离散采用格心格式的有限体积法，流场数据储存在网格体心处，网格K的体积为VK，其体心XK可以表示为

(7)

与式(6)类似，对于任意的多面体网格，其体积和体心都可以分解为四面体来计算，即

VK=∑Vi

(8)

式中：Vi和Xi分别为四面体的体积与体心。

一般流体计算中认为流场变量呈二阶分布，将流场变量Taylor展开：

(9)

可以证明，流场变量在K单元的平均值就等于体心值，流场变量在K单元的积分值等于体心值与体积的乘积,即

qK

(10)

3.2 守恒插值方法

以超网格为媒介，可以实现重叠网格子网格之间的守恒插值，具体方法如下。

式(2)是重叠网格守恒插值条件，将其在背景网格ΓA、插值网格ΓB上离散化，即

(11)

根据超网格的特点，插值网格ΓA的每一个网格单元KB总是可以拆分为一个或多个超网格单元的集合，设ΓC为其超网格，网格单元为KC，则式(11)可写为[15]

(12)

这样就把ΓA到ΓB的守恒插值转换到ΓA到ΓC的守恒插值。根据超网格的特点，背景网格ΓA的每一个网格单元KA总是可以拆分为一个或多个超网格单元的集合。超网格单元KC的体积为VKC，则其体心XKC可表示为

(13)

设KC为贡献单元KA所包含的一个超网格单元，可以利用二阶Taylor重构[6,27-28]将任意变量插值从贡献单元体心XKA插值到超网格体心XKC处，即

(14)

设贡献单元KA包含m个超网格，可以证明通过式(9)的插值方式，变量q在KA上的积分可写为在KA所包含的m个超网格上的积分之和，即

(15)

综合式(11)、式(12) 、式(15)，背景网格ΓA的流场变量可通过超网格ΓC插值到插值网格ΓB，即

(16)

式(15)整个插值过程对二阶分布的流场变量严格守恒，具体实现过程总结如下：

1) 将流场变量从贡献单元KA插值到它的超网格单元KC，即

(17)

2) 将超网格ΓC的流场变量采用体积权插值到插值单元KB，设KB包含n个超网格单元KC，则其流场变量q(XKB)可写为

(18)

这样以超网格为媒介，将流场变量从贡献单元KA插值到插值单元KB。如果流场变量q(x)二阶分布，则整个插值过程严格守恒，这种插值方法不会改变原先流场的二阶分布，实现了超网格单元和插值单元流场变量的精确守恒插值。

3.3 隐式并行插值方法

在重叠网格插值过程中，根据贡献单元是否同时也是插值单元，插值方法可分为隐式插值和显式插值。显式插值方法需区分背景网格与插值网格，插值速度较慢，不适合复杂构型的大规模重叠网格。本文发展了一种隐式并行插值方法，不严格区分背景网格与插值网格，如图6所示，每一套网格既向其他网格发送数据，同时又从其他网格接收数据。

这种隐式并行的重叠网格插值方法插值速度较快，适用于多个相对运动的子网格间的流场数据插值，能处理复杂构型的大规模重叠网格。

4 算例分析

本文首先采用不同类型的流场方程对所建立的插值方法的守恒性和插值精度进行验证。然后将该方法用于CFD重叠网格标模的数值模拟，通过跟试验数据对比，考核该方法的鲁棒性和工程实用性。

4.1 守恒性验证

4.1.1 线性方程

假设物理量φ在网格区域线性分布，满足方程：

φ(x,y,z)=2x+3y+z

(19)

如图7所示，ΓA、ΓB分别为结构与非结构四面体网格，网格区域一致。在网格上设定物理量φ的分布满足式(19)，然后ΓA、ΓB间进行100次相互插值，统计φ在ΓB网格区域内的积分值与每个单元的平均插值误差绝对值。

图8对比了本文插值方法与三线性插值方法、逆向距离权插值方法在守恒性和插值精度上的异同。图8(a)为φ在ΓB网格区域内的积分值，可以看出用本文守恒插值方法计算得出的积分值在插值过程中保持恒定，按照式(11)守恒插值的定义，该方法在计算过程中可保持物理量积分的严格守恒。而三线性插值方法和逆向距离权插值方法不具备守恒性，且随着插值次数的增加，物理量积分值与理论值的误差越来越大。

图8(b)为φ在ΓB网格区域内每个单元的平均插值误差绝对值，可以看出与三线性插值方法和逆向距离权插值方法相比，本文所建立的插值方法插值精度更高。

4.1.2 三角函数

假设物理量φ在网格区域呈三角函数分布，满足方程：

φ(x,y,z)=1+sin(2πx)sin(2πy)sin(2πz)

(20)

插值网格与计算条件同4.1.1小节。图9对比了三角函数的插值过程中，本文插值方法与三线性插值方法、逆向距离权插值方法在守恒性和插值精度的异同。可以看出，本文发展的守恒插值方法对三角函数的插值计算具有较强的守恒性(见图9(a))，且插值精度明显高于三线性插值方法与逆向距离权插值方法(见图9(b))。

4.2 定常数值计算

4.2.1 ONERA-M6机翼

将重叠网格数值计算结果与单套整体网格数值计算结果对比是一种常用的检验重叠网格技术优劣的方法[7,23]。为了考察本文重叠网格插值方法的插值精度，对ONERA-M6机翼分别采用非结构重叠网格与单套整体网格进行数值模拟，计算来流条件为:马赫数Ma=0.84, 迎角α=3.06°, 来流温度T∞=255.56 K, 雷诺数Re=1.172×107。

重叠网格由两部分组成：背景网格和机翼网格，机翼表面采用各向异性网格对机翼前后缘进行加密，机翼近壁面采用三棱柱网格以较好地模拟黏性，远壁面采用四面体、金字塔和六面体的混合网格，机翼网格量约为60.3万。背景网格全部采用四面体网格单元，网格量约为6.1万。两套子网格的重叠效果如图10所示，总网格量为66.4万。

单套全机网格的近壁面网格分布与重叠网格完全一致，空间网格尺度和网格数量与两套重叠网格相当，总网格量为54.8万。

如图11(a)和图11(b)所示，本文将不同重叠网格插值方法的压力系数Cp计算结果与相同网格尺度单套整体网格的计算结果进行对比，可以看出相比三线性插值方法与逆向距离权插值方法，本文所提出的插值方法计算得出的压力系数和单套整体网格的计算结果更为符合,其中L为机翼展长。由于本文建立的插值方法提高了重叠网格交界面的插值精度，使得激波附近的压力系数计算更为精确。

从图11(a)和图11(b)可以发现CFD计算结果与风洞试验结果符合良好，前缘吸力峰和激波位置与风洞试验数据较为一致。但仔细对比发现，整体网格与守恒插值方法计算得出的Cp与试验值略有偏差。ONERA-M6机翼是一个跨声速机翼，如图11(c)和图11(d)所示，大量的数值计算[29-33]表明，对于ONERA-M6机翼跨声速状态，全局加密网格、S-A模型、DES和DDES方法的数值计算结果都无法与风洞试验完全符合，这是数值计算和风洞试验的固有偏差[29-30]。而三线性与距离权插值方法在激波后部与风洞试验完全符合，这是完全不合理的，也恰好说明这两种重叠网格插值方法数值误差较大。本文建立的守恒插值方法计算结果与整体网格、文献[29-33]一致，说明该方法插值精度较高，产生的数值误差较小。

图12给出了用不同重叠网格插值方法计算的残差收敛曲线，可以看出由于本文守恒插值方法在重叠网格区域插值精度较高，引入误差较小，使得收敛速度较快。本文插值方法的收敛速度分别比逆向距离权方法和三线性方法快20.0%和10.8%。

4.2.2 NACA0012翼型激波干扰

取2个并列错位的NACA0012翼型进行超声速扰流计算，来流条件为:Ma=2,α=0°,Re=1.812×107。

在此计算条件下，由于翼型并列错位，2个翼型之间会产生强烈的激波干扰，流场结构比较复杂，本文采用重叠网格对此算例进行数值模拟，如图13所示,重叠网格由两部分组成：上翼型网格和下翼型网格。

为了考核插值方法在重叠网格边界的插值精度，人为地将计算网格分布得较为稀疏，在重叠网格边界,网格的尺度相差较大。

图14(a)～图14(c)给出了不同插值方法计算得到的马赫数等值线分布，白色、黑色分别为上、下翼型网格的马赫数等值线，图14(d)为细网格计算结果，3种插值方法都正确模拟了流场结构。可以看出，在上下翼型重叠网格交界面附近，用三线性和逆向距离权插值方法计算得到的马赫数等值线有明显的错位，而守恒插值方法得到的马赫数等值线较为光滑。在下面翼型的尾流区域，人为设置了两套子网格的尺度，使得在重叠网格边界两套网格的尺度相差较大，在这种情况下，三线性和逆向距离权插值方法计算得到的等值线出现了明显的间断，而守恒插值方法得到的马赫数等值线仍然较为连续。此外，2个NACA0012翼型并列错位，产生了复杂的激波干扰流场，相比三线性和逆向距离权方法，用守恒插值方法计算得到的激波分布更为精细、光滑。

4.3 三维机翼分离投放模型

选择机翼/挂架/带舵外挂物(Wing/Pylon/Finned-Store，WPFS)模型，对本文所建立的重叠网格守恒插值方法进行验证。该模型是美国进行分离投放CFD数值验证时使用的一个机翼和导弹的简化模型，有详实可靠的风洞试验数据，因此常被当做验证重叠网格分离投放计算的标准算例。外挂物共有4片尾翼，呈X形分布。在该试验中，还使用了弹射力，使外挂物在分离初期保持抬头趋势，避免其尾部和挂架碰撞，具体数据详见文献[34]。

计算状态取和文献[34]一致：来流马赫数Ma=0.95，初始迎角α=0°，初始侧滑角为0°，本文采用非定常雷诺平均Navier-Stokes方程、S-A湍流方程、6自由度(DOF)方程来模拟计算该模型的外挂物分离历程。采用6套重叠子网格进行计算，分别为：用于模拟计算的包裹机翼和挂架的背景网格、包裹带挂物的子网格以及包裹每一个舵的子网格，总网格量为553万。WPFS模型子网格的重叠效果如图15(a)所示，表面网格如图15(b)所示。

图16(a)是初始状态在执行“挖洞”操作后得到的重叠网格计算单元的剖面图，物理量通过相互重叠的网格部分进行插值传递。图16(b)是时间t=0.3 s时各个子网格相互重叠示意图，可以看到在计算过程中，子网格跟随物体位置的变化而变化。图17给出了t=0 s时刻与t=0.3 s时刻的实时压力分布云图，计算得出的流场结构清晰合理。

将计算得到的外挂物分离动态特性同试验数据进行比较，图18(a)和图18(b)分别为外挂物的质心位移随时间的变化量以及对应时刻的姿态角。可以看出，计算得到的轨迹同试验数据吻合较好，证明本文建立的重叠网格守恒插值方法准确有效，鲁棒性强，工程实用。

5 结 论

1) 针对二阶线性分布的流场变量，采用局部超网格方法在重叠网格区域可实现精确守恒插值，对非线性分布的流场变量，本文插值方法守恒性更高。

2) 结合局部超网格与ADT搜索方法，重叠网格插值单元能够快速搜索到多个贡献单元，合理地扩大了贡献单元模板，提高了插值精度和流场光滑性，加快了残差收敛速度。

3) 本文重叠网格守恒插值方法适用于非结构网格与结构网格的各类型网格单元，工程实用、鲁棒性强。插值时采用体积加权方法，当子网格尺度相差较大时，仍然能够保持较好的流场光滑性与连续性。

本文发展了一种重叠网格守恒插值新方法，下一步将开展重叠网格洞面优化技术与嵌入网格方法，进一步提高重叠子网格尺度相差较大时流场的光滑性。

参 考 文 献

[1] 李鹏, 高振勋, 蒋崇文. 重叠网格方法的研究进展[J]. 力学与实践, 2014, 36(5):551-565.

LI P, GAO Z X, JIANG C W. The progress of the overlapping grid techniques[J]. Mechanics in Engineering, 2014, 36(5):551-565(in Chinese).

[2] SHIH T I P. Overset grids: Fundamentals and practical issues: AIAA-2002-3259[R]. Reston, VA: AIAA, 2002.

[3] 田书玲, 伍贻兆, 夏健. 用动态非结构重叠网格法模拟三维多体相对运动绕流[J]. 航空学报, 2007, 28(1):46-51.

TIAN S L, WU Y Z, XIA J. Simulation of flows past multi-body in relative motion with dynamic unstructured overset grid method[J]. Acta Aeronautica et Astronautica Sinica, 2007, 28(1): 46-51 (in Chinese).

[4] 伍贻兆, 田书玲, 夏健. 基于非结构动网格的非定常流数值模拟方法[J]. 航空学报, 2011, 32(1):15-26.

WU Y Z, TIAN S L, XIA J. Unstructured grid methods for unsteady flow simulation[J]. Acta Aeronautica et Astronautica Sinica, 2011, 32(1): 15-26 (in Chinese).

[5] MARSTIN C, MCCONNAUGHEY H. Computational problems on composite grids: AIAA-1984-1611[R]. Reston, VA: AIAA, 1984.

[6] KANG Z L, YAN C, YU J, et al. A fast and reliable overset unstructured grids approach[J]. Acta Mechanica Sinica, 2013, 29(2): 149-157.

[7] 黄宇, 阎超, 王文, 等. 混合重叠网格插值方法的改进及应用[J]. 北京航空航天大学学报, 2017, 43(2):285-292.

HUANG Y, YAN C, WANG W, et al. An improved interpolation method for hybrid overset grid and its application[J]. Journal of Beijing University of Aeronautics and Astronautics, 2017, 43(2): 285-292 (in Chinese).

[8] 周乃春, 李彬, 郑鸣, 等. 带控制率导弹投放数值模拟[J]. 空气动力学学报, 2013, 31(3): 107-115.

ZHOU N C, LI B, ZHENG M, et al. Missile separation simulation with control laws[J]. Acta Aerodynamica Sinica, 2013, 31(3): 107-115 (in Chinese).

[9] FARRELL P E, PIGGOTT M D, PAIN C C, et al. Conservative interpolation between unstructured meshes via supermesh construction[J]. Computer Methods in Applied Mechanics and Engineering, 2009, 198(33): 2632-2642.

[10] ZHENG Y, LIOU M S. Progress in the three-dimensional DRAGON grid scheme: AIAA-2001-2540[R]. Reston, VA: AIAA, 2001.

[11] XU K, SUN G, CAI J Y. On interface conservative attributions for computations of the complex flows of high-lift system based on chimera tech-nique[C]∥Computer and Automation Engineering (ICCAE), 2010 The 2nd International Conference on IEEE. Piscataway, NJ: IEEE Press, 2010: 279-283.

[12] 张宇飞, 陈海昕, 符松. 基于高阶守恒重映对窗口嵌入技术的改进[J]. 计算物理, 2011, 28(2): 167-173.

ZHANG Y F, CHEN H X, FU S. Improvement on window embedment technology with high order conservative remapping[J]. Chinese Journal of Computational Physics, 2011, 28(2): 167-173(in Chinese).

[13] ZHAO X, GUAN H W, YANG Z, et al. An implicit and globally conservative unstructured chimera grid method: AIAA-2011-0777[R]. Reston, VA: AIAA, 2011.

[14] FARRELL P E, MADDISON J R. Conservative interpolation between volume meshes by local Galerkin projection[J]. Computer Methods in Applied Mechanics and Engineering, 2011, 200(1): 89-100.

[15] MENON S, SCHMIDT D P. Conservative interpolation on unstructured polyhedral meshes: An extension of the supermesh approach to cell-centered finite-volume variables[J]. Computer Methods in Applied Mechanics and Engineering, 2011, 200(41): 2797-2804.

[16] 徐春光, 董海波, 刘君. 基于单元相交的混合网格精确守恒插值方法[J]. 爆炸与冲击, 2016, 36(3): 305-312.

XU C G, DONG H B, LIU J. An accurate conservative interpolation method for mixed grid based on the intersection of grid cells[J]. Explosion and Shock Waves, 2016, 36(3): 305-312 (in Chinese).

[17] 崔鹏程, 邓有奇, 唐静, 等．基于伴随方程的网格自适应及误差修正技术[J]．航空学报, 2016,37(10): 2992-3002.

CUI P C, DENG Y Q, TANG J, et al. Adjoint-based grid adaptation and error correction[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(10): 2992-3002 (in Chinese).

[18] BLAZEK J. Computational fluid dynamics: principles and applications[M]. 3rd ed. Oxford: Elsevier, 2015:75-120.

[19] KIM J S, KWON O J. Improvement on block LU-SGS scheme for unstructured mesh Navier-Stokes computations: AIAA-2002-1061[R]. Reston, VA: AIAA, 2002.

[20] SPALART S R. A one-equation turbulence model for aerodynamic flows: AIAA-1992-0439[R]. Reston, VA: AIAA, 1992.

[21] ROE P L. Approximate Riemann solvers, parameter vectors, and difference schemes[J]. Journal of Computational Physics, 1981, 43(2): 357-372.

[22] 李彬, 唐静, 邓有奇, 等. 并行的多重网格方法在离散伴随优化中的应用[J]. 航空学报, 2014, 35(8): 2091-2101.

LI B, TANG J, DENG Y Q, et al. Application of parallel multigrid algorithm to discrete adjoint optimization[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(8): 2091-2101 (in Chinese).

[23] BEATRICE R, JAY S. Robust and scalable overset grid assembly for partitioned unstructured meshes: AIAA-2013-0797[R]. Reston, VA: AIAA, 2013.

[24] THEOHARIS T, IAN P. Two parallel methods for polygon clipping[J]. Computer Graphics Forum, 2010, 8(2):107-114.

[25] MARTINEZ F, RUEDA A J, FEITO F R. A new algorithm for computing Boolean operations on polygons[J]. Computers & Geosciences, 2009, 35(6): 1177-1185.

[26] 唐静, 邓有奇, 马明生, 等. 飞翼气动优化中参数化和网格变形技术研究[J]. 航空学报, 2015, 36(5): 1480-1490.

TANG J, DENG Y Q, MA M S, et al. Parametrization and grid deformation techniques for fly-wing shape optimization[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(5): 1480-1490 (in Chinese).

[27] ALAUZET F, MEHRENBERGER M. P1-conservative solution interpolation on unstructured triangular meshes[J]. International Journal for Numerical Methods in Engineering, 2010, 84(13): 1552-1588.

[28] GARIMELLA R, KUCHARIK M, SHASHKOV M. An efficient linearity and bound preserving conservative interpolation (remapping) on polyhedral meshes[J]. Computers & Fluids, 2007, 36(2): 224-237.

[29] MAYEUR J, DUMONT A, DESTARAC D, et al. RANS simulations on TMR test cases and M6 wing with the ONERA elsA flow solver: AIAA-2015-1745[R]. Reston, VA: AIAA, 2015.

[30] DURRANI N, QIN N. Comparison of RANS, DES and DDES results for ONERA M6 wing at transonic flow speed using an in-house parallel code: AIAA-2011-0190[R]. Reston, VA: AIAA, 2011.

[31] DA SILVA R G, AZEVEDO J L F, BASSO E. Simulation of ONERA M6 wing flows for assessment of turbulence modeling capabilities: AIAA-2016-0549[R]. Reston, VA: AIAA, 2016.

[32] SAXENA S K, NAIR M T. Implementation and testing of Spalart-Allmaras model in a multi-block code: AIAA-2002-0835[R]. Reston, VA: AIAA, 2002.

[33] LUO H, BAUM J D, LÖHNER R. A Hermite WENO-based limiter for discontinuous Galerkin method on unstructured grids[J]. Journal of Computational Physics, 2007, 225(1): 686-713.

[34] HEIM E R. CFD wing/pylon/finned store mutual interference wind tunnel experiment: AEDC-TSR-91-P4[R]. New York: AEDC, 1991.