计算流体力学(CFD)已成为空气动力学研究和应用的主要手段之一。影响CFD精度的两个主要因素是数值方法和计算网格。数值方法研究较多且应用经验丰富，但由于复杂工程外形的网格生成强烈依赖用户经验，计算网格已成为影响CFD精准度的关键因素^[1]。近年来，耦合流场信息、自动调整网格分布的自适应技术得到了广泛的研究，且在NASA报告“2030年CFD展望研究”中将其列为亟需发展的一种网格技术^[2]。

结构网格通常只能采用移动网格点的r型自适应技术^[3-4]，优化网格的能力有限。h型自适应技术可动态增减网格单元数量，能有效地优化网格分布，包括笛卡尔网格、四面体网格和混合网格3种类型。笛卡尔网格可基于八叉树的网格生成算法进一步开展网格自适应^[5]，但模拟附面层难度很大。四面体自适应技术发展较为完备^[6]，然而也存在附面层模拟精度不高的问题。

非结构混合网格兼顾了黏性模拟精度和易于生成的要求，广泛应用于黏性流动的模拟，但临近物面的层结构网格是自适应技术的难点。基于多面体单元类型的自适应方法^[7-9]具备同时改变物面法向和切向的网格分布能力，并且有效消除了相邻网格单元间的相互影响，鲁棒性较高。

为快速解决实际飞行器的气动问题，CFD流场求解通常都采用大规模的并行技术。网格自适应也须相应地采用匹配的并行技术，否则将成为自适应CFD模拟流程中的性能瓶颈。

CFD流场求解的并行算法较为成熟，较多研究机构^[10-11]和超算中心^[12-13]都已实现了CFD流场求解的大规模并行。流场求解并行实现通常采用区域分解策略，即在负载平衡的约束下进行网格分区。为了减少内存的开销，流场求解时每个并行进程只存储局部区域的网格，各并行进程在分区交界面附近进行流场数据的通信^[10]。

当前，网格自适应并行技术的研究主要集中于h型自适应技术需要的动态负载平衡技术，用以解决因不同进程网格单元增减数量差异大造成的后续流场计算负载不平衡的问题^[14-15]。动态负载平衡算法可大致分为扩散法^[16-19]和两步法^[20-26]两种。扩散方法采用迭代思想，每步只针对并行边界两侧的网格单元实施迁移操作，直到达到全局负载平衡的指标，一般耗时较长。两步法则是通过重分区和数据迁移先后两个过程实现。根据重分区是否并行，可分为管理进程重分区^[20-22]和并行重分区^[23-26]。前者须首先汇总所有进程的网格数据，继而基于串行算法进行分区，最后将分区后的网格分发到各自进程。该方法可达到很高的负载平衡性能，但由于涉及大量网格数据的汇总和分发，并行效率通常较低且时间开销较大。并行重分区则是所有进程都采用并行分区算法(如ParMetis程序^[27]等)进行协同分区，然后将少量的网格数据传递到目标进程，虽然负载平衡性能略差，但通常效率更高和计算时间更少。

网格自适应技术一般包括误差估计、网格单元分布优化、表面网格几何投影、空间网格协调匹配和流场插值等主要步骤^[9]。为了提高整个自适应系统的并行效率，需要针对每个步骤涉及的自适应技术设计相应的并行算法。

本文基于分布式并行构架及消息传递接口(MPI)，针对混合网格自适应系统中的关键技术开展并行算法研究，主要包括网格单元分布优化并行技术、表面网格投影并行技术和空间网格协调并行技术，以及动态负载平衡技术。首先，在网格单元分布优化方面，提出了“先唯一后同一”的两步法策略实现了并行增添网格对象，采用预留编号区间方法保证了新增网格对象编号的唯一性，并设计了网格单元和网格点的特征数据恢复了网格对象的同一性。其次，采用局部曲面拟合技术实现了曲面重构和新增物理网格点投影的完全并行，基于数据通信实现了新增虚拟网格点的投影。第三，提出了半并行算法实现了空间网格匹配过程的并行化，快速解决了网格投影造成的网格单元交错问题。第四，发展了适用于基于单元剖分自适应方法的动态负载平衡技术，采用基于根节点网格单元及子单元加权的方法实现了并行重分区，采用数据打包/解包技术提高了网格数据迁移的并行性能，继而采用圆柱激波流场自适应模拟验证了动态负载平衡技术。最后，采用三角翼网格自适应加密测试了本文建立的自适应系统的并行效率。

1 网格自适应系统的并行技术

网格自适应技术通常包括网格单元误差估计、自适应网格优化和网格流场插值3个部分。图 1给出了结合网格自适应的CFD流场并行模拟的执行流程，并行算法涉及各模块的数据交换。其中，网格误差估计和流场插值都是基于当前网格和进程间并行对应关系开展几何量和流场数据的代数运算，其并行实现可直接使用流场求解的并行算法。自适应网格优化可细分为网格单元分布优化、表面网格几何投影和空间网格协调匹配3个部分，是自适应系统并行算法设计的核心。

在基于网格分区的CFD并行实现中，虚拟网格技术常用来存储并行流场数据，以实现物理边界条件和并行边界条件的统一处理^[10]。本文虚拟网格单元的模板扩展为与网格点相连的网格单元，如图 2所示，实线围成的网格单元为真实物理单元，虚线相关的网格单元为虚拟单元。

通常情况下，并行算法与自适应系统采用的具体方法有关。本文针对非结构混合网格，且主要基于文献[9]中h型混合网格自适应实现方法，分别设计了自适应网格优化3项技术的并行算法。另外，针对网格加密或稀疏造成的负载不平衡，本文在两步法的基础上设计了适合多级网格自适应的动态负载平衡算法。

1.1 网格单元分布优化并行技术

采用基于单元剖分的h型自适应方法实现网格单元分布优化时，各进程并行地剖分网格单元会同时新增网格单元和网格点。为了准确建立进程间的并行连接关系，新增的网格对象(包括网格单元和网格点)需满足如下两条并行相容性准则：

1) 唯一性准则：真实物理网格对象的编号和并行进程编号在所有进程中具有唯一性。

2) 同一性准则：虚拟网格对象的编号和并行进程编号与存在于其他进程上的真实网格对象的编号和并行进程编号相同。

并行进程编号的并行相容性实现相对简单，网格单元和网格点分别按如下方法设置：通过剖分得到的新增子网格单元的进程编号与其父级网格单元相等；新增网格点的进程编号等于与该网格点相连的所有网格单元的进程编号的最小值。

本文虚拟网格单元也一并开展剖分或粗化操作，以减少子网格信息的并行通信。由于新增网格对象的同一性准则只能通过并行通信实现，同时考虑唯一性和同一性准则会由于并行通信将并行操作退化为串行。因此，本文采用两步法策略先后实现网格对象编号的并行相容性，即在新增网格对象时先保证唯一性准则，完成后再统一恢复新增网格对象的同一性。

唯一性准则要求为新增的网格对象分配全局唯一的编号。串行环境下，可以采用计数器为每个新增网格对象依次分配唯一编号。然而，并行环境下计数器分配编号将会出现冲突。本文采用预留编号区间方法实现并行编号。主要的思想是每个并行进程计数器的取值限定在互不重叠的编号区间内，只在某个并行进程的计数器取值超出了可用区间时才需要通过并行通信确定新的可用区间。并行进程p的可用编号区间[I_min, I_max]可确定为

$ \left\{ {\begin{array}{*{20}{l}} {{I_{\min }} = {N_{\max }} + (p - 1) \times {N_{\rm{T}}} + 1}\\ {{I_{\max }} = {I_{\min }} + {N_{\rm{T}}}} \end{array}} \right. $	（1）

式中：N_max为当前所有进程的计数器的最大值；n_p为并行进程总数，p=1, 2, …, n_p；N_T为计数器区间跨度，本文取512。新增网格点也可按类似方法并行地进行编号分配。

各并行进程采用预留编号区间方法独立编号新增网格对象后，每个网格单元或网格点的编号都满足全并行域内的唯一性，如图 3所示，包括虚拟网格单元(图中用字母c和d表示)和虚拟网格点(图中用字母p表示)。

本文同一性的恢复采用基于网格对象特征数据的并行通信来实现，包括“询问”和“返回”两个主要步骤。网格特征数据是区分不同网格对象的唯一性标识，在“询问”步骤通过特征数据在真实物理网格对象所在的进程查询对象编号，然后在“返回”步骤设置对应的对象编号。

对于新增网格单元的特征数据D_e，取为父级单元编号I_p和子单元序号S_c组成的数对：

$ {D_{\rm{e}}} \buildrel \Delta \over = {\mathop{\rm pair}\nolimits} \left\langle {{I_{\rm{p}}}, {S_{\rm{c}}}} \right\rangle $	（2）

对于新增网格点的特征数据D_p，取如下四元素：

$ {D_{\rm{p}}} \buildrel \Delta \over = \left\{ {{C_1}, \quad {C_2}, \quad {C_3}, \quad {C_4}} \right\} $	（3）

新增网格点位于原网格单元的相对位置包括体心、面心和网格边3种情况，分别设置为

1) 体心：C₁为原网格单元编号，C₂=C₃=C₄=Big，Big为程序语言允许的最大正整数，下同。

2) 面心：若为三角形面，C₁、C₂和C₃分别取为该三角形3个顶点的编号，C₄=Big；若为四边形面，C₁、C₂、C₃和C₄分别取为该四边形4个顶点的编号。

3) 网格边：C₁和C₂分别取为网格边两个端点的编号，C₃=C₄=Big。

不同进程上定义同一网格面或网格边所采用的网格点的顺序可能不同，因而同一网格点的特征数据包含的4个元素的顺序也可能不同。对4个元素进行升序排序，则两个特征数据相等判断只需逐一判断每个元素。另外，虽然特征数据分为3种情况分别设置，但得到的特征数据互不相同，因而在并行算法设计时可采用同一数据结构。图 3中新增的网格单元和网格点恢复编号的同一性后，如图 4所示。

以上并行相容性算法可能会造成网格对象的编号不连续，包括同一进程内编号不连续和全并行域内某些编号未使用。因此，在实现了新增网格对象的并行相容后，可基于并行对应关系重新排序所有物理网格对象，并通过并行通信同步虚拟网格对象在重排序后的新编号。

1.2 表面网格几何投影并行技术

当表面网格投影基于CAD系统时，每个进程可以并行地调用CAD系统进行投影。当采用基于曲面拟合方法重构几何信息时，曲面拟合方法的可并行性直接影响并行算法的设计。通常非局部的曲面拟合方法的可并行性较低，并行通信量大，并行效率低。曲面拟合采用模板少且紧致的方法，如文献[9]中局部Coons曲面拟合方法，有利于降低并行算法实现难度和提高并行性能。

对于包含新增网格点的表面网格单元，局部Coons曲面拟合只需要表面网格单元的顶点和顶点处的法向信息，后者可采用点相连的网格面法向信息加权得到。因此，该拟合方法只涉及当前表面网格单元和顶点相连的临近表面网格单元。

根据本文虚拟网格单元的设置方法可知，真实的物理表面网格单元，如图 5中的网格面F，需要的所有表面网格都已存在于当前的并行进程中，因而曲面拟合不需要并行通信。同时，由于重构的曲面为参数化曲面，新增表面网格点的投影目标点也不需要并行通信。因此，物理表面网格单元上新增网格点投影可以完全独立执行。

对于虚拟表面网格上的新增网格点，曲面拟合需要其他并行进程上的表面网格信息，必须通过并行通信获取。然而，考虑到该虚拟表面网格单元同时作为物理网格单元存在于其他并行进程上，且新增网格点的目标点也已经求出，因此只需通过并行通信直接同步投影的目标点即可。

1.3 空间网格匹配并行技术

空间网格匹配是用来解决因表面网格投影造成的空间网格单元交叉等问题，通常采用动网格技术将表面网格点的位移传递到内部网格点上。空间网格匹配并行技术的性能与采用的动网格技术的并行算法直接相关。

针对文献[9]中的动网格技术，可以采用半并行算法^[28]，其主要思想是每个并行进程都收集完整的表面网格点数据，独立地处理各自的内部空间网格点。

2 流场求解方法简介

本文流场求解采用课题组自主研发的MFlow软件^[29]。该软件基于积分型式的雷诺平均Navier-Stokes方程和采用格心格式的有限体积方法。文中采用Roe格式^[30]计算对流通量和采用中心差分格式计算黏性通量，采用LU-SGS方法^[31]进行方程的迭代求解，湍流模拟采用标准SA模型^[32]。

3 动态负载平衡技术

通常情况下，增添或删除网格单元的操作主要集中在少数几个并行进程中，造成后续流场求解的负载严重不平衡，将显著降低并行效率。本文基于并行重分区-网格数据迁移的两步法思想，发展了适用于多层级网格自适应^[9]的动态负载平衡技术。

在完成增添或删除网格单元后，多层级网格单元间的连接关系存在跨层级情况，构建连接关系的算法较为复杂。在并行环境下构建单元连接关系，不同的网格单元需要不同层级的虚拟网格单元。为了避免虚拟网格单元的遴选过程带来的计算开销和降低算法复杂度，本文将所有层级的虚拟网格单元都通过并行通信进行传递，以适应所有物理网格单元的连接关系的构建。

并行重分区采用ParMetis程序库^[27]，分区所需的网格连接图的构建基于最初的根节点网格单元，根节点网格单元的权重取为该单元包含活动子网格单元(流场求解需要)的总数。该方法能有效保证后续流场求解的负载平衡性能。在完成根节点网格单元的并行进程编号设置后，所有子网格单元的并行进程编号设置为对应的根节点网格单元的进程编号值。

网格重分区完成后，不再属于当前进程的网格单元，以及临近的虚拟网格单元，都需要通过并行通信迁移到目标进程。网格单元的数据包含整数和实数等多种数据类型，且通常较为分散地存储在内存中。为了减少并行通信函数的启动频率，提高通信效率，本文采用数据打包/解包的方法，如图 6所示，将多种数据类型、分散的网格数据转为连续存储，增大并行通信的数据粒度。

此处采用二维圆柱超声速绕流算例，开展基于梯度误差估计^[9]的网格自适应流动模拟，以初步验证本文建立的动态负载平衡技术。来流条件为：马赫数为2.0，静压为30 800 Pa，静温为230 K。采用4核并行计算，网格自适应共迭代3次。图 7对比了初始网格和自适应3次后计算的流场，自适应后激波和膨胀波都捕捉得更精细。自适应网格负载平衡的动态变化过程见图 8，图中同一颜色的网格归属相同的并行进程。可以看出，为了保持后续流场求解的负载平衡，网格并行分区的交界面随着自适应迭代动态变化。另外，改变归属进程的网格单元都集中于并行交界面附近的位置，有利于降低网格单元数据迁移的通信开销。

4 网格自适应并行技术效率测试

为了测试本文自适应系统的并行技术在大规模并行计算中的并行效率，采用尖前缘三角翼模型。同时，对自适应系统3个主要模块的并行性能也进行了统计，包括网格分布优化、物面几何投影和空间网格匹配。测试来流条件为：雷诺数为1.38×10⁷，静温为322 K，马赫数为0.4，迎角为13.3°，侧滑角为0°。由于来流对称，采用包含对称面的半模网格开展流动模拟，初始网格约2 500万，网格单元包含四面体、金字塔和三棱柱3种类型。网格自适应采用均匀加密所有单元，自适应后网格单元总量约为2.0亿。测试的并行核数分别为：16、32、64、128、256和512。由于计算机内存限制，此处以16核作为并行效率计算的基准，并行加速比S和并行效率E计算方法分别为

$ \left\{ {\begin{array}{*{20}{l}} {S = \frac{{{T_0}}}{{{T_n}}}}\\ {E = \frac{S}{{n/{n_0}}} \times 100\% } \end{array}} \right. $	（4）

式中：T为一定迭代步数计算的平均时间；n为并行核数；下标0表示相应的基准值。

在128核以内，整个网格自适应系统和3个主要模块的并行加速比差别较小，如图 9所示。更大并行规模时，网格分布优化和表面网格投影两个模块的加速比优于整个自适应系统，空间网格匹配模块加速比最低。该模块加速比性能最低的原因是本文采用半并行算法，每个并行进程都针对所有表面网格点开展插值运算，降低了该并行算法在较大并行规模时的并行性能。

图 10给出了整个网格自适应系统和3个主要模块的并行效率，整体来看，128核以内自适应系统的并行效率大于80%，表明自适应系统有较高的可扩展性。更大并行规模时，空间网格匹配的并行效率下降最快，而表面网格投影模块并行效率最高。由于网格投影模块中曲面拟合和物理网格点投影都不需要并行通信，仅需要并行同步虚拟网格点，因而该模块的并行效率很高。

不同并行核数时，网格自适应系统、3个主要模块和并行通信耗费的总时间及比例见表 1，3个主要模块中，表面网格投影模块计算时间最少，小于总时间的1%，网格分布优化模块计算时间最长。图 11给出了耗时比例随并行规模增加的变化，表面网格投影模块耗时比例略有减小，网格分布优化模块耗时比例逐渐减小，空间网格匹配模块略有增加。另外，并行通信时间所占比例随并行规模较为快速地增加，表明该网格自适应系统属于通信密集型并行系统。

图 12给出了采用256核求解流场的密度残差和升力系数的收敛曲线，迭代约11 000步后残差下降约3个量级，气动力获得收敛解。此时，流场求解共需约3.2×10⁴ s(约8.9 h)。使用256核完成网格自适应共需200 s，约为流场求解的0.6%。因此，尽管256核时网格自适应并行效率不高，但从总时间开销看是完全可接受的。

5 结论

针对非结构混合网格，设计了网格单元分布优化、表面网格几何投影和空间网格匹配相关的并行算法。同时，为后续流场模拟建立了基于并行重分区-网格数据迁移的动态负载平衡技术。最后，采用三角翼模型测试了不同并行规模下3项并行算法的并行性能。

1) 提出的“先唯一后同一”两步策略有效实现了网格单元分布优化模块中网格单元和网格点的并行增添或删除，确保了并行环境下网格的并行相容性。

2) 采用的局部曲面拟合技术仅需同步虚拟网格点的信息，物理网格点投影每个进程完全独立，该方法具备很高可并行性。

3) 建立的并行自适应系统在大规模并行时并行效率有所降低，但对于定常流动模拟，与流场求解相比，总时间消耗很少，该系统的并行性能完全可以接受。

本文提出的两步法策略和采用局部曲面拟合技术的思想，对其他的网格自适应方法也具有一定的借鉴意义。