0 引言

随着计算机科学技术的不断发展，计算流体力学(Computational Fluid Dynamics，CFD)已经成为与风洞试验、飞行试验互为补充的三大空气动力学研究手段之一，被广泛应用于航空航天飞行器设计过程之中。

长期以来，作为经典物理中的世纪性难题，湍流被认为是一个重大的基础科学问题，其流动机理至今仍是一个未解之谜^[1]。利用CFD模拟层流流动已经较为成熟，但是面对实际流动中的湍流问题，由于其物理本身的高度复杂性，只能做简化处理。当前，湍流模拟主要有三类方法：雷诺平均N-S方程(RANS)、大涡模拟(LES)、直接数值模拟(DNS)。DNS能直接进行流动的数值模拟，但是由于需要极其庞大的计算资源，在现在和可预见的未来都难以用于模拟真实外形；LES可精确模拟流动分离，但是其计算资源仍然耗费巨大。RANS是目前唯一能够直接用于实际工程计算的方法，且对常规流动问题的模拟能取得令人满意的结果，是当前CFD工程应用中大量采用的方法。

目前的工程应用中大部分RANS数值模拟都是采用SA和k-ω SST两类湍流模型，这两类模型方程中均涉及壁面距离，即空间点到物体表面的最近距离。在计算前的预处理过程中，需要首先计算壁面距离。壁面距离计算方法主要有两类：一类是通过几何方法直接计算，另一类是通过求解偏微分方程的形式求解^[2-4]。对于求解偏微分方程的方法，在简单几何外形上速度较快，但是在真实复杂外形上难以应用，主要原因有：(1) CFD代码开发者除了要开发流动方程的求解器外，还要开发一个额外的求解器，带来额外的工作负担；(2)只要是偏微分方程求解类问题，都存在稳定性问题，而在复杂外形上常常难以收敛，容易算发散；(3)偏微分方程求解得到的是近似解，而壁面距离的精确性直接影响到CFD求解器的稳定性、收敛性和计算结果的精度(下文将详细阐述)。因此，目前实际应用中基本上都是用几何方法直接搜索计算。

直接搜索法，即对于空间中的任意一个点，搜索壁面上的离散单元(如三角形、四边形)离其最近的距离。搜索过程中，需要一一比较与壁面上的每个单元的距离以找到最近点。对于N个空间点、M个壁面面元，搜索比较次数为N×M次，对于二维问题来说，这种搜索方式尚且能忍受，但是对于三维真实航空航天飞行器外形，N值一般为千万、上亿量级，M值一般为百万、千万量级，在复杂应用中甚至到达数十亿，直接搜索将耗费大量时间，是整个CFD模拟所无法承受的。针对上述效率问题，研究者提出了几种高效的计算方法。

NSMB^[5]、Spalart^[6]等介绍了一种沿网格线找到离其最近的壁面点的方法，但是这种方法只适应于结构网格计算，无法应用到非结构网格中，算法不具备通用性。王刚等^[7]提出了一种逐层推进的方法，基本原理是从壁面开始，向流场逐层推进搜索最近距离，该方法适用于结构网格和非结构网格，但是实现过程比较繁琐，编程困难，而且难以并行化。

赵慧勇等^[8]介绍了一种基于二分法将壁面循环划分为若干个盒子的方法，即循环盒子法，首先计算离空间点最近的若干个盒子，再在这些盒子中搜索最近的壁面点。循环盒子法的本质是将整个壁面用八叉树划分，每个方盒子相当于一个八叉树节点，从而实现了串行直接搜索方法的大幅加速。这种方法适用于任意网格，但是这种方法仍存在两方面的不足：(1)算法难以并行化; (2)当网格量大幅增加后，效率将大幅下降。这是因为，如果盒子数过少，则每个盒子中的壁面面元数将大幅增加，每个盒子中的搜索相当于直接搜索，导致效率大幅下降；相反，如果盒子数过多，在计算离空间点最近的盒子过程中需要首先进行盒子距离排序，这一排序过程将耗费大量时间，甚至与盒子内单元直接搜索时间相当。上述两个问题导致循环盒子法的扩展性受限。杨永国等^[9]、李广宁等^[4]对这一方法中的盒子划分方法、搜索方法进行了改进提高，但算法上并无本质区别。

Boger^[10]采用了一种高效的ADT数据结构，利用三角不等式，设计了一种壁面距离的快速计算方法，可以有效利用空间中点的坐标信息，实现对直接搜索法的大幅加速。但是正如李广宁^[4]所述，这种方法存在两个缺点，一是搜索半径难以确定，二是不平衡的ADT树导致外流计算问题中壁面距离计算效率大幅下降。

上述方法在几何直接搜索法上进行了不同的改进，相对于原始方法而言，计算效率有了大幅的提升，对于一般工程应用中的百万、千万量级计算网格是适用的。但是，随着设计人员对模拟分辨率要求的不断提高、以及如DES等CFD计算方法不断进入工程应用，目前工程应用中的计算网格量不断增加，航空航天飞行器设计型号任务的网格量接近上亿量级，甚至更多。另外，以“神威·太湖之光”、“天河二号”等为代表的高性能计算机硬件不断发展，目前CFD软件也在不断进行硬件适配，以适应于超大规模并行计算，计算网格量可达十亿、百亿量级。使用大规模网格带来的问题是，上述壁面距离计算方法已完全无法满足大规模并行计算的预处理需求。

针对上述问题，本文提出了一种基于ADT数据结构的湍流模拟壁面距离并行计算方法，适用于十亿、百亿量级网格的大规模数值模拟问题。

1 壁面距离精确计算的重要性 1.1 湍流模型方程基础

CFD工程应用中常用的湍流模型方程有一方程SA模型和两方程k-ω SST模型。这里以SA模型为例，说明湍流壁面距离对计算结果的影响。SA湍流模型方程为^[11-12]：

$ \begin{aligned} \frac{\partial \widetilde{\nu}}{\partial t}=&-u_{j} \frac{\partial \widetilde{\nu}}{\partial x_{j}}+\frac{1}{\sigma}\left[\nabla \cdot(\nu+\tilde{\nu}) \nabla \tilde{\nu}+c_{b 2}(\nabla \tilde{\nu})^{2}\right]+\\ & C_{b 1} \widetilde{S}_{\nu}-C_{w 1} f_{w}\left(\frac{\tilde{\nu}}{d}\right)^{2} \end{aligned} $

(1)

上式写成非守恒物质导数形式为：

$ \begin{aligned} \frac{\mathrm{D} \tilde{\nu}}{\mathrm{D} t}=& C_{b 1} \widetilde{S} \tilde{\nu}+\frac{1}{\sigma}\left[\nabla \cdot(\nu+\tilde{\nu}) \nabla \tilde{\nu}+c_{b 2}(\nabla \tilde{\nu})^{2}\right]-\\ & C_{w 1} f_{w}\left(\frac{\tilde{\nu}}{d}\right)^{2} \end{aligned} $

(2)

这里，$\widetilde{v}$为SA模型方程变量，d为空间点到壁面的最小距离，也就是本文所要计算的对象。在式(2)右端，第一项为湍流生成项，第二项为扩散项，第三项为壁面破坏项。其中，扩散项又由非守恒扩散项和守恒扩散项组成。该模型的壁面边界条件为：

$ \begin{aligned} &\widetilde{\nu}_{\text {wall }}=0\\ &\widetilde{\nu}_{\text {farfield }}=3 \nu_{\infty} \sim 5 \nu_{\infty} \end{aligned} $

(3)

这里，下标wall、farfield分别指代固体壁面(即物面)和远场。式(3)求在壁面上的湍流变量为0。SA方程是一个精巧的数学模型，通过式(2)中右端三项相互调节，使整个流场内的湍流黏性达到平衡状态。尤其是，为了与式(3)“壁面上湍流变量为0”这一边界条件相匹配，该方程通过右端三项驱动，使得空间湍流变量沿壁面法向，从空间到壁面逐渐衰减至0，在紧邻壁面的一层网格的湍流变量值是一个接近0的小值。此外，定常状态计算收敛时，要求湍流变量不再随时间变化，即要求式(1)右端几项之和为0。在未收敛前，右端几项之和往往不等于0，将其定义为“残差”。

1.2 壁面距离精确性的影响

三维情况下，壁面被离散为三角形或四边形网格集，空间点的精确壁面距离是点到三角形面片或四边形面片集的最近距离。一般情况下，为了编程方便、提高计算速度，会采用近似方法计算壁面距离，例如上文提到的求解偏微分方程的方法，或者几何上的近似计算，即将点到三角形/四边形面心的最近距离来代替壁面距离，从而带来计算误差。

壁面距离计算的精确性对湍流模拟的健壮性、计算结果精度有重要的影响。这里以二维湍流平板算例为例说明。图 1(a)是算例边界条件，下侧x>0为壁面、x < 0为对称面，左右侧为压力入/出口，上侧为远场。图 1(b)是计算网格，流向和法向网格维度为35×25，为了说明壁面距离精确性，将紧邻壁面的第一层网格划分为各向异性三角形。要注意的是，为了观察方便，图 1(b)中的x和y方向采用不同的长宽比可视化，真实的第一层三角形网格的长宽比为1×10⁷量级的高度拉伸的各向异性三角形。

对于图 1(b)中的每一个三角形或四边形计算网格，精确壁面距离是其y坐标。这里，首先采用精确的壁面距离将该算例计算到收敛状态，使得湍流方程(1)的右端项即残差接近于0，然后在第一层的三角形网格中人为引入壁面距离误差，观察对湍流方程的影响。图 2是给定第一层三角形网格不同的误差后，对右端各项以及残差的影响，可以看到，当误差增加后，与壁面距离相关的破坏项迅速变化，残差从0迅速增加到1×10³量级。从局部放大图看，即使误差有小幅增加，也将使得残差被快速放大。

SA方程中通过壁面距离来判断空间网格离壁面的相对远近，第一层网格的壁面距离本应该相对小，残差应该是接近于0的小值，其壁面距离值被人为增大后，使得其本来具备的“近壁面”属性丧失而被“误认为”是“远壁面”单元，从而使得湍流模型方程中的壁面破坏效应迅速丧失，将严重影响CFD计算的收敛性，甚至导致计算发散。为了解决这一问题，在一些CFD代码中，通过强行将壁面上的湍流方程通量修正为0的方式来提高计算鲁棒性，然而，这样可能会带来另一个严重的问题：计算无法完全收敛。图 3中显示了人为修正壁面通量对平板湍流模拟的影响，图 3(a)是x=1.05处边界层内沿y向的湍流方程右端项分布。如果从整体看分布图没有任何异常，但是，当仅观察壁面附近的几层网格(图 3b)时，可以看到，即使对已经收敛后的流场进行残差修正，也会使得壁面附近的1~2层网格上的湍流项失真，湍流场在迭代计算过程中不断呈现“收敛-修正-收敛-修正”的循环，使得湍流方程求解无法完全收敛。而当去掉残差修正后(图 3c)，湍流方程收敛时其右端三项之和为0，符合SA模型方程的控制机理。

总的来说，壁面距离的精确与否，将直接影响计算的稳定性和计算精度。对于不精确的壁面距离，如果在湍流方程求解过程中采用类似残差修正等手段，虽然能保证计算顺利进行，但是会影响计算收敛性，从而影响计算结果的准确性。

2 计算几何基础

一般情况下，真实的壁面外形在非间断处是光滑连续的，在CAE或者计算机辅助几何设计里，一般采用参数曲线、曲面描述。在CFD计算中生成计算网格时，会将这些连续的参数化曲面离散为三角形/四边形面片，对于壁面距离计算问题，将转化为计算空间点到三角形/四边形面片最近距离的求解问题。一些CFD代码为了降低实现难度，提高计算效率，往往会采用近似的几何方法来计算壁面距离，最常用的是直接用空间单元中心到壁面面元面心的距离来求得最近壁面距离。本文尝试采用较为精确的计算方法。几何方法计算壁面距离的基础是一些基本的几何代数运算，这里先简要介绍本文所述方法中采用的计算几何基础^[13]。

壁面距离计算的通用情况是计算空间点到凸N边形的距离，为了提高计算精度，将凸N边形的每条边和多边形中心连接，形成N个三角形，最终转换为计算点到N个三角形的距离。分两步计算点到三角形的距离：首先将空间点投影到三角形所在平面，然后判断投影点是否在三角形内部。

图 4是点投影示意图，将空间点P投影到△V₀V₁V₂所在的平面$\mathscr{P}$得到投影点P_I。图中，投影距离d = r · n，投影点P_I=P-d· n，即将空间点沿逆法向方向投影到平面$\mathscr{P}$。

图 5中，采用参数化坐标方法判断投影点P_I是否在三角形内^[14]。参数化的平面$\mathscr{P}$可以表示为式(4)。

$ \begin{aligned} V(s, t) &=V_{0}+s\left(V_{1}-V_{0}\right)+t\left(V_{2}-V_{0}\right) \\ &=V_{0}+s \mathit{\boldsymbol{u}}+t \mathit{\boldsymbol{v}} \end{aligned} $

(4)

上式中，u =V₁-V₀和v =V₂-V₀分别是三角形的两条边向量。s和t是参数化坐标，对于任意点P_I=V(s, t)，若位于三角形中，必满足下式：

$ s \geqslant 0, \quad t \geqslant 0, s+t \leqslant 1.0 $

(5)

因此，如要判断投影点P_I是否在三角形内，只需要找到其参数化坐标s和t即可。设w=V₁-V₀，P_I的参数化坐标s和t可通过下式计算：

$ \begin{aligned} &s_{I}=\frac{(\boldsymbol{u} \cdot \boldsymbol{v})(\boldsymbol{w} \cdot \boldsymbol{v})-(\boldsymbol{v} \cdot \boldsymbol{v})(\boldsymbol{w} \cdot \boldsymbol{u})}{(\boldsymbol{u} \cdot \boldsymbol{v})^{2}-(\boldsymbol{u} \cdot \boldsymbol{u})(\boldsymbol{v} \cdot \boldsymbol{v})}\\ &t_{I}=\frac{(\boldsymbol{u} \cdot \boldsymbol{v})(\boldsymbol{w} \cdot \boldsymbol{u})-(\boldsymbol{u} \cdot \boldsymbol{u})(\boldsymbol{w} \cdot \boldsymbol{v})}{(\boldsymbol{u} \cdot \boldsymbol{v})^{2}-(\boldsymbol{u} \cdot \boldsymbol{u})(\boldsymbol{v} \cdot \boldsymbol{v})} \end{aligned} $

(6)

通过式(6)计算得到参数化坐标，然后通过式(5)判断投影点P_I是否在三角形内，若在三角形内，则投影距离d是空间点P投影到△V₀V₁V₂的距离。特别地，若投影点不在所有物面三角形内部，则计算空间单元中心到物面点/三角形中心的最近距离来作为壁面距离，这种情况极少。

3 大规模壁面距离快速计算

本文的大规模壁面距离计算方法几个基本要素包括并行计算方法、基于ADT数据结构的搜索、MPI/OpenMP混合并行方法。

3.1 并行计算壁面距离

现有的壁面距离计算方法基本上都是串行的，一些学者提出的并行计算方法只是在某个步骤、分任务式的并行计算，其本质还是串行方式，可扩展性不强。随着计算网格规模的增加，当达到上亿甚至十亿、百亿量级后，现有串行计算方法都会失效，必须采用并行计算。

CFD中的并行计算基本上都是基于区域分解的，即将离散后整个块计算区域分解为不同的空间子分区，这些子分区再被分配给不同的进程。非结构网格上的一个分区是指一个网格块，而结构网格的分区可能包括多个网格块。本文的分区包括两个概念——“空间网格分区”和“壁面网格分区”。空间网格分区，是采用开源分区软件METIS将原始的整个计算域划分后的子分区，分区内的单元是四面体/六面体/金字塔等三维拓扑单元。壁面网格分区，是指空间网格分区中包含的壁面面元构成的子分区，分区内的单元是三角形/四边形等二维拓扑单元。图 6是球体外形计算网格的空间网格分区和壁面网格分区。空间网格分区用于CFD模拟过程中分配给每个进程，而壁面网格分区天然地成为并行壁面计算时的分区。

并行计算壁面距离的流程如图 7。对于N个进程并行的情况，进程i的分区如果含有壁面面元，则将这些壁面面元压缩到数据容器后发送到服务器进程(一般为0号进程)。服务器进程将所有进程的壁面面元收集后，将全局壁面面元发送给其他的所有进程，从而使得每个进程中都有全局壁面面元信息。图 8是宽体客机的壁面信息图。该飞行器的计算网格被分为2048个子区，其中的若干个分区含有壁面面元，将这些壁面面元由服务器进程收集后发送给每个进程，使得每个进程都含有图中的壁面面元信息的备份。最后，每个进程中的计算网格根据其上的全局壁面面元信息，计算每个空间网格单元中心到全局壁面的最近距离。

这种并行计算流程的优点是，每个进程在进行壁面距离计算时，由于其中含有全局的壁面信息，因此可以采用几乎现有的所有方法计算壁面距离，如直接搜索法、循环盒子法等，改造简单，适应性和扩展性好。

3.2 基于ADT数据结构的搜索算法

图 7的并行壁面距离计算流程中，每个进程需要计算其中的单元中心到壁面面元的最近距离，如果采用直接搜索法，计算时间将取决于网格单元数最多的进程，大规模计算的每个分区单元数在数十万量级，而壁面面元数可能达到千万甚至亿量级，直接搜索法的遍历次数是二者的乘积，计算时间可能长达数天，效率存在瓶颈。本文采用了一种基于ADT(Alternating Digital Tree)数据结构的搜索算法。

ADT数据结构实际上是一种特殊的二叉树^[15]，其基本搜索原理是：对于每一个空间多边形，将其几何坐标的最小盒子，依次按照笛卡尔坐标三个方向插入到空间中。搜索时，在空间点的几何坐标盒子中查找含有的空间多边形，这种搜索方式是二分法在多维情况下的扩展。该方法由于效率极高而在计算几何方法中得到广泛的应用。例如，在网格生成、几何求交等算法中，能实现给定空间范围中元素的快速查找。

如果直接利用ADT数据结构搜索，会出现类似文献[4]中提到的问题——当搜索对象离搜索目标距离较远时(如外流问题的远场区域)，采用ADT搜索将面临极为尴尬的问题，即在给定的搜索半径下要么搜索不到目标，要么搜索到的目标过多而退化为直接搜索法。针对这一问题，提出一种迭代搜索方法。首先，建立全局壁面面元的ADT树，建立过程是将壁面面元的坐标方盒逐个插入到ADT树的过程，详细可参考文献[15]。然后，遍历当前进程所含分区中的所有单元，基于ADT查找计算最近的壁面距离。图 9是ADT搜索过程，其中的N_max为人为给定的最大ADT节点数(其中存储的数据是壁面面元)，其目的是限制搜索到的节点数，使其不至于过多而影响效率。

表 1是ADT迭代搜索算法描述。ADT查找的核心是确定搜索半径，对于复杂外形，很难通过几何外形直接计算得到搜索半径。本文采用的这种迭代搜索方法，只需要在最开始给定一个初始搜索半径，迭代过程中根据搜索到的ADT节点数量不断地二分调整搜索半径，直到搜索到理想的节点数为止。以图 8的搜索为例，第1次搜索时由于搜索半径过小，没有搜索到壁面面元；将其搜索半径扩大2倍后进行第2次搜索，但导致搜索的壁面面元数超过了节点限制数N_max；第3次搜索时，搜索半径设定为前两次搜索半径之二分之一后，搜索到的壁面面元数满足需求，然后在这些搜索到的节点内计算最近的壁面距离。

3.3 MPI/OpenMP混合并行计算

上面介绍了基于区域分解的并行壁面距离计算方法。一般来说，每个进程上分配一个网格分区，在这种条件下，要求每个进程中都必须包含全局壁面的数据备份。这里的“全局壁面数据”包含两类数据，一是所有壁面上点的坐标、单元与点的拓扑关系，二是壁面面元构建的ADT树数据。

当进行大规模网格计算时，壁面面元数可能达到上千万，每个进程中都要存储全局壁面数据将会导致计算节点的内存不够。可以采用的解决方法是，在一个节点内只存储一个全局壁面信息，节点中的网格分区逐个、串行地利用全局壁面信息计算壁面距离，但如此一来将导致一个节点内只有一个进程在进行壁面距离计算，其他进程处于空闲状态。

为此，设计了一种MPI/OpenMP混合并行的壁面距离计算模式。CFD软件通常运行在分布式集群上，采用MPI进行进程间的通信。对于多核处理器来说，节点内部的内存是共享的，如果在节点内部也采用MPI通信，可能会由于带宽的限制带来通信总量的增加，采用单一的MPI通信往往无法得到最好的效果。其中一种解决方案是采用MPI/OpenMP混合并行，即在节点内部用OpenMP并行，而在节点之间用MPI并行，实现节点间和节点内部的两级并行，充分利用消息传递和共享内存两种编程模型的优点。每个进程在进行壁面距离计算时，各进程无需通信，因此可以在节点内部对分区的循环遍历实施OpenMP并行。当节点内的各分区都计算完壁面距离后，在节点间采用MPI进行数据同步，从而实现MPI/OpenMP混合并行。表 2是混合并行计算算法。

4 计算方法测试

从精度和速度两方面，分别在国产定制集群和“天河二号”集群上，对上述的方法进行了测试。测试平台是国产通用CFD软件平台PHengLEI^[16-17](前身为HyperFLOW)，该软件是作者所在团队研发的具有完全自主知识产权的通用CFD软件平台，主要面向航空航天飞行器设计应用部门。软件采用C++语言编程、面向对象的软件设计理念。设计的体系结构和数据结构具有良好的可扩展性和通用性，以适应不同类型网格(结构网格、非结构网格、混合网格、重叠网格)的计算。在该软件平台上同时开发了结构、非结构求解器，二者可独立运行，也可同步耦合计算，即在流场中同时含有结构、非结构网格的情况下，在结构网格上调用结构求解器，在非结构网格上调用非结构求解器。

4.1 计算精度测试

为了测试不同计算方法的精度，采用了结构/非结构两套高度拉伸的各向异性网格。(1)结构网格。将图 1中的二维四边形网格直接在展向拉伸，得到三维六面体结构网格。(2)非结构网格。将图 1中的二维四边形网格三角化，然后展向拉伸，形成图 10中高度拉伸、高度扭曲的各向异性三棱柱网格。这两套网格都有很大的拉伸率，在壁面附近，长宽比达到2×10⁷。二者的不同点是网格正交性，结构网格在法向具有很好的正交性，全场网格单元的夹角都几乎为90°；而三棱柱化后的各向异性非结构网格具有很大的扭曲度，在壁面附近最小夹角仅为0.27°。采用一套正交性很好、一套正交性很差的网格，其目的是测试精确方法和近似方法的壁面距离计算精度。

在两套网格中，y=0平面即为壁面，各棱柱单元中心的y向坐标即为其精确壁面距离，可以用于测试壁面距离计算的误差。

图 11是不同方法计算的壁面距离比较，其中，“Y coordinate”是y方向的坐标也是壁面距离精确值，“Approximation”是直接用单元中心到壁面面元面心的距离来求壁面距离的近似方法，这种方法在很多代码中都用到。从图 11(a)的结构网格计算结果看，采用直接搜索方法和ADT搜索方法计算得到的壁面距离和精确值几乎完全重合，“Approximation”近似法虽然在远离壁面处有一些误差，但是在壁面附近(y接近0)和精确值符合很好。从图 11(b)的非结构网格计算结果看，同样地，直接搜索方法和ADT搜索方法计算得到的壁面距离和精确值几乎完全重合，但是“Approximation”方法计算的壁面距离在壁面附近误差很大。

正如第1.1节所述，对于RANS方程湍流模拟而言，其鲁棒性和收敛性主要受壁面附近壁面距离精确性的影响。从精度测试结果看，传统方法中常用的“Approximation”近似计算方法，对于具有良好正交性的网格，计算的壁面距离在壁面附近能满足湍流模拟的需要，而在网格正交性不好的情况下，计算精度差，将严重影响CFD模拟结果。因此，近似法受网格质量的影响很大，严重依赖于网格生成人员的经验。而对于直接搜索法和ADT搜索法，计算的壁面距离精度都是满足需要的，受网格质量的影响较小。

4.2 计算效率测试

测试算例是第三届美国AIAA高升力CFD预测活动中提供的JSM(JAXA Standard Model)标模外形。该外形是客机、运输机等大型飞行器的简化研究外形，主要用于评估CFD技术对高升力外形的预测能力，为工程应用提供经验指导。为了说明方法的可扩展性，分别采用了粗/密两套不同规模的计算网格在不同的运行环境下进行了测试。图 12是该外形的几何外形及计算网格。

表 3是计算网格、测试环境和计算效率对比。计算网格为棱柱/四面体/金字塔组成的混合网格，其中，粗网格单元为0.52亿，共2048个空间网格分区，在国产定制集群上测试；密网格单元为33.2亿，共12288个空间网格分区，在“天河二号”上测试。对于两套网格，分别采用直接搜索方法和ADT搜索方法计算了壁面距离，二者都采用并行方式计算。可以看到，在两套网格、两种测试环境下，并行ADT搜索方法都有大幅的效率提升。尤其是，直接搜索法时间将随网格量的增加而指数增长，当网格量达到上亿量级后，壁面距离计算很可能耗费与CFD模拟过程差不多的时间，这将是不可接受的。

正如上文所言，由于ADT搜索过程中每个进程都需要有全局壁面信息的备份，可能带来内存开销的大幅增加，上文的MPI/OpenMP混合并行的目的正是为了解决内存开销问题。表 4中给出了粗网格在不同线程/进程数情况下的内存耗费和测试时间，可以看到，即使对于粗网格而言，如果用纯MPI并行计算(2048×1)，每个节点内的16个CPU核都拥有一个全局壁面信息备份，单节点内存耗费将达到23.2%；而在MPI/OpenMP混合并行情况下，如果1个节点内只用2个MPI进程(256×8)，单节点内存耗费降到7.8%。对于密网格，如果采用纯MPI并行计算，其内存耗费约为120 GB，将远远超过节点内存，必须采用混合并行计算。

对于直接搜索法和ADT搜索法，虽然都可以采用3.3节的MPI/OpenMP混合并行计算方案，但是二者在混合并行时有很大的差异。为了说明这一点，采用了ONERA-M6机翼外形(图 13)，在国产定制集群上进行了测试。计算网格为棱柱/四面体/金字塔组成的混合网格，共46万单元，分为了64个空间网格分区，包括43个壁面网格分区，壁面三角形总数为1.89万。之所以采用这个相对较小的算例来比较，是由于直接搜索法计算效率低，如果用大规模算例测试时间太长。表 5是不同的MPI进程数、OpenMP线程数情况下两种方法的测试结果。对于直接搜索法，随着MPI进程数减少、OpenMP线程增加，计算时间大幅增加。而对于ADT搜索法，计算时间除了4进程×16线程外，其他情况下计算时间没有明显区别。

出现以上结果的原因是：对于直接搜索法，每个空间点都要遍历所有的壁面面元，因此，如果多个线程访问同一个内存中的壁面信息备份，必然会出现数据竞争；对于ADT搜索法，由于是根据其空间坐标关系遍历其坐标周围一定搜索半径的区域，因此搜索时只和该区域中的其他对象产生竞争，而根据ADT数据结构，即使在同一个区域中的空间对象，也是依次在3个笛卡尔坐标方向轮换存储于内存中的，又进一步减少了数据竞争的机会。对于4进程×16线程情况，意味着一个节点内只有一个壁面信息备份，该壁面信息被全部16个线程访问，会增加数据竞争的概率，即使是ADT搜索方法，其计算效率也会受到影响。

5 结论

针对超大规模网格下CFD湍流模拟中壁面距离计算效率低的问题，设计了一种基于ADT搜索的并行计算方法，采用了MPI/OpenMP混合并行方法解决内存开销过大的问题，实现了复杂外形大规模网格壁面距离高精度、高效并行计算。通过33.2亿网格大规模并行测试，结果表明，即使与并行化的直接搜索方法相比，新方法在效率上也有量级的提升，计算时间完全能够满足大规模CFD湍流模拟的需求，解决了面向下一代E级CFD计算的前处理问题。

该方法还存在待改进之处。ADT搜索实际上是一种三维的二分法，其搜索“方盒”的各个面天然地与Oxyz笛卡尔坐标面平行，若几何外形的某处局部也与笛卡尔坐标面平行，将导致ADT搜索效率下降，即要么搜索不到任何点, 要么搜索到该平面中所有的点。解决方法是，将几何坐标旋转一个角度(小于90°)，从而使得壁面面元不再与笛卡尔坐标面平行。下一步工作将继续开展深化研究。