1、文章编号：0258-1825(2023)07-0093-09笛卡尔网格下精确高效的壁面距离计算方法孟爽1，2，周丹1，李雪亮2，毕林2，3，*(1.中南大学轨道交通安全教育部重点实验室，长沙410075；2.空气动力学国家重点实验室，绵阳621000；3.中国空气动力研究与发展中心计算空气动力研究所，绵阳621000)摘要：对于笛卡尔网格方法，壁面距离是采用虚拟单元法精确处理物面边界的重要参数，同时也是网格自适应后制约流动计算效率的关键因素之一。针对现有壁面距离计算方法结果不精确、效率不高的问题，引入三角形参数化方法，将空间点到三角形物面离散网格的最小距离问题转换为约束条件下一维极值问题，仅需
2、通过符号判断和少量加减乘运算，即可确定最小距离，计算精度和效率大幅提高；发展嵌套包围盒概念的 KDT（K-dimensionaltree）物面网格数据存储结构，优化 KDT 最近邻搜索算法中距物面较远数据点回溯过程，实现了最小距离对应的三角形的快速定位。运用球、导弹、DPW6 等三维几何构型对上述方法考核验证结果表明，计算得到的壁面距离与解析值的误差在百万分之一以内，十亿量级网格规模下的单核计算效率接近已有文献中的并行计算效率。关键词：笛卡尔网格；壁面距离；计算效率；KDT；回溯方法中图分类号：V211.3文献标识码：Adoi：10.7638/kqdlxxb-2021.0359Accurate
3、 and efficient wall distance calculation method for Cartesian gridsMENGShuang1，2，ZHOUDan1，LIXueliang2，BILin2，3，*(1.Key Laboratory of Traffic Safety on Track(Central South University),Ministry of Education,Changsha410075,China；2.State Key Laboratory of Aerodynamics,Mianyang621000,China；3.Computationa
4、l Aerodynamics Institute of China Aerodynamics Research and Development Center,Mianyang621000,China)Abstract:ThewalldistanceofCartesiangridsisanessentialparameterfortheproperwalltreatmentusingghostcellsandisalsooneofthecriticalfactorsgoverningtheefficiencyoftheflowfieldsimulationaftermeshadaptation.
5、Thispaperproposesatriangularparameterizationmethodthatconvertstheproblemofcomputingtheminimumdistancebetweenspatialpointsanddiscretizedtriangularmeshesonthesurfaceintoaconstrainedone-dimensionalextremumproblem.Thissimplificationonlyrequiressymbolicjudgmentsandasmallnumberofaddition,subtraction and m
6、ultiplication operations to obtain the minimal distance,yielding significantimprovementsintheaccuracyandefficiencycomparedtoexistingmethods.Meanwhile,aKDT(K-dimensionaltree)datastructurebasedonthenestedenclosingboxconceptisdevelopedtooptimizethebacktrackingofdatapointsfarfromthesurfaceintheKDTneares
7、tneighborsearchalgorithm.Theapplicationofthismethodtothree-dimensional geometries such as spheres,missiles,and DPW6 demonstrates that the error between thecomputedwalldistanceandtheresolvedoneiswithinonemillionth.Moreover,thecomputationalcostsofusingasinglecoreforbillion-scalegridsarecomparabletotho
8、seofparallelcomputationusingexistingmethods.Keywords:Cartesiangrid；walldistance；computationalefficiency；K-dimensionaltree；backtracking收稿日期：2021-11-08；修订日期：2022-03-14；录用日期：2022-04-12；网络出版时间：2023-05-11基金项目：国家数值风洞工程（NNW2018-ZT1A02）；中南大学研究生自主探索创新项目（206021722）作者简介：孟爽（1994），男，河南许昌人，博士研究生，研究方向：网格生成技术.E-mai
9、l：通信作者：毕林*，研究员，研究方向：转捩、湍流.E-mail：引用格式:孟爽,周丹,李雪亮,等.笛卡尔网格下精确高效的壁面距离计算方法J.空气动力学学报,2023,41(7):93101.MENGS,ZHOUD,LIXL,etal.AccurateandefficientwalldistancecalculationmethodforCartesiangridsJ.ActaAerodynamicaSinica,2023,41(7):93101(inChinese).doi:10.7638/kqdlxxb-2021.0359第41卷第7期空气动力学学报Vol.41，No.72023年7月AC
10、TA AERODYNAMICA SINICAJul.，2023符号说明符号说明G物面三角形单元集合NG中的元素个数i维度（本文指x、y、z三个维度）GiminGimax(,)几何的最小包围盒dKDT的深度u当前划分维（x，y，或z）Tcenter三角形中心KiminKimax(,)节点K对应的包围盒顶点KL节点K的左子节点KR节点K的右子节点KLiminKLimax(,)节点K的左子节点对应的包围盒KRiminKRimax(,)节点K的右子节点对应的包围盒Kis节点K在第i维度上剖分面对应的值lmin最小距离Knearest最小距离对应的节点m剖分轴 0 引言笛卡尔网格可以不依赖物面网格直接生
11、成，因具有自动化程度高、复杂外形适应性好、非定常/多尺度流动结构捕捉能力强等优势而广受关注1。壁面距离的精确高效计算对实现笛卡尔网格方法具有重要意义：一方面，由于笛卡尔网格的非贴体特性，通常采用虚拟单元法处理物面边界2-6，需要根据物面确定虚拟单元的壁面距离，寻找参考点并插值获得其物理量，壁面距离精确计算对虚拟单元物面处理精准度有较大影响；另一方面，对于非定常流动或运动物体，笛卡尔网格会根据流场的变化或物体的运动进行自适应7-8，自适应后的网格壁面距离需重新计算，壁面距离计算效率成为制约流动计算效率的关键因素之一。目前，壁面距离的计算方式主要分为求解偏微分方程9-12和采用几何方法直接计算13
12、-14两类。第一类方法，将最小距离转换为关于波传播问题的偏微分方程数值求解15，需要额外的计算花费，壁面距离计算精度受数值离散精度的影响，且对于复杂外形适应性较差。第二类方法是根据空间点与物面离散网格的几何关系计算壁面距离。通常为简单起见，在物面网格足够密的情况下，求解空间网格点的壁面距离时，一般采用空间点到物面网格中心的距离代替最小距离。赵钟等16发现，这种近似距离会引起计算误差，不仅影响计算结果的精度，而且对计算稳定性造成影响。为提高计算精度，有学者发展了根据空间点关于物面网格所在平面投影点与物面网格的几何关系计算壁面距离的 3D 投影法16、通过坐标变换转化为二维问题的 2D 法17等，
13、但这些方法涉及大量向量计算和关系判断，计算量较大。在几何方法求解壁面距离时，由于物面网格数目较大，计算效率的核心问题是如何快速定位到最短距离相关的物面网格，物面网格数据结构是关键。最常见的方法是将物面网格顺序存入数组，采用遍历法搜索物面点得到壁面距离。网格规模较小时，计算量尚在承受范围之内。但是对于三维复杂几何外形的流场计算，物面网格数量可达到 O(105)O(106)量级，空间网格点可达到 O(108)O(109)量级，采用遍历法计算量达到 O(1013)O(1015)量级。此种规模的计算量对整个计算周期的影响是巨大的。Boger18提出采用 ADT(alternatingdigitaltr
14、ee)二叉树数据结构存储物面网格，计算效率大幅提高，成为目前最常用的物面网格数据存储方式。但对于复杂外形，ADT 平衡性下降，计算效率仍然差强人意。郭中州等14采用 KDT(K-dimensionaltree)存储物面网格，解决了 ADT 平衡性问题，计算效率进一步提升，但对于离物面较远的空间点，在运用 KDT 最近邻搜索算法时，超球面与分割面位置关系判断失准，需要反复回溯，计算量大。本文针对上述问题，引入三角形参数化方法，将空间点到三角形物面离散网格的最小距离问题转换为约束条件下一维极值问题；同时发展嵌套包围盒概念的 KDT 物面网格数据存储结构，优化 KDT 最近邻搜索算法中距物面较远数据
15、点回溯过程，实现最小距离对应的三角形的快速定位。1 点到三角形最小距离精确算法T(s,t)=B+sE0+tE1(s,t)D=(s,t):s 0,t 0,s+t 1T(s,t)点到空间三角形最小距离的问题19可描述为点P 和三角形之间的最小距离，其中，B 为三角形的一个顶点，E0和 E1分别为此顶点对应的三角形两条边。点 P 和三角形内任一点距离的平方为椭圆函数：Q(s,t)=|T(s,t)P|2=as2+2bst+ct2+2ds+2et+f(1)(s,t)Da=E0E0b=E0E1c=E1E1d=E0(BP)e=E1(BP)f=(BP)(BP)其 中，。点 P 和三角形之间的关系如图 1 所示
16、。94空气动力学学报第41卷Q(s,t)最小距离问题转换为在 D 内求连续可微函数的极值问题。令：Q=2(as+bt+d,bs+ct+e)(2)s=(becd)/(acb2)t=(bdae)/(acb2)Q=(0,0)(s,t)DT(s,t)Q=V0s+t=1当且时，有。如果，最小距离为点 P 到的距离；否则，最小距离出现在 D 的边界上。为了找到正确的边界，采用三角形区域 D 将 s-t 平面划分为七个区域，如图 2 所示。其中区域对应三角形区域 D，表示 Q 对应的椭圆与相切。tst0s0Q=VV0Q=V0Q=011(s,t)图 2 用三角域D划分s-t平面以及不同等高线Fig.2 Par
17、tition of a s-t plane by triangle domain Dand distance contours(s,t)T(s,t)如果在区域，则与点 P 最近的点在三角形上。(s,t)s 0,1Q(s,t)s+t=1s0 0,1Q(s0,1s0)如果在区域，则与点 P 最近的点计算问题转换为在内等值线与直线的接触问题（可能是相切，也可能是相交于两个端点），即寻找，使得达到极小值，这是个一维极值问题，仅需通过符号判断和少量加减乘运算即(s,t)可确定。在区域或区域时，处理方式与区域类似。(s,t)s 0,1Q(s,t)s+t=1s 0,1Q(s,t)s=0如果在区域，则与点 P
18、 最近的点计算问题存在两种可能：1）在内等值线与直线的接触问题；2）在内等值线与直线的接触问题。如图 3 所示。ts101 Q Q图 3 以区域为中心的等高线与三角形接触的两种情况Fig.3 Two contour levels with centers in region touching the triangle在等值线上的任何一点，Q 的方向均朝向椭圆内部。为简单起见，我们根据Q 沿顶点(0,1)的两条边的方向导数的符号进行判断：s+t=1(1,1)/2)Q(0,1)s+t=1s 0,11）对于边，顶点(0,1)处的方向导数为。如果值为负，最小值出现在边上。最小值的情形对应图 3 红色曲
19、线，处理方式和在区域类似。s=0(0,1)Q(0,1)s=0(s,t)2）对 于 边，顶 点(0,1)处 的 方 向 导 数 为。如 果 值 为 负，最 小 值 出 现 在 边上。最小值的情形对应图 3 绿色曲线，处理方式和在区域类似。(s,t)(s,t)在区域或区域时，处理方式与在区域类似。stT(s,t)=B+sE0+tE1 st求得 和 之后，将其带入即可得到点到三角形最小距离的交点，从而构建虚拟单元方法中的“镜像点”。此种方法将最小距离求解问题转换为一维极值问题，在计算距离时只在求 或时进行一次浮点除法，其他均为简单的符号判断或加减乘运算，相比于 3D 投影法16或 2D 法17，在保
20、证计算精度的同时计算量相对较小。2 物面单元的存储 2.1 嵌套包围盒概念的 KDT 参数定义与文献 14 中只存储物面网格中心点不同，本文BT(s,t)PE0E1图 1 点P和三角形之间的关系Fig.1 The relationship between P and the triangle第7期孟爽等：笛卡尔网格下精确高效的壁面距离计算方法95发展了基于嵌套包围盒概念的 KDT，不仅按逆时针顺序存储物面网格顶点信息（方便笛卡尔网格物面处理），还同时包含 KDT 节点对应所有物面网格的最小包围盒，KDT 节点从上到下形成嵌套包围盒形式。嵌套包围盒的引入保证 KDT 每个节点包含当前节点对应三角
21、形的所有信息，且在最小距离查询时，仅需要通过与剖分面（剖分面的大小由包围盒确定）坐标对比，就可以快速排除与目标点距离大于当前最小距离的三角形，效率大幅提升。2.2 嵌套包围盒概念的 KDT 建立步骤步骤 1：在正式创建 KDT 之前，首先需要确定目标数据集 G 中三角形元素的数量 N。然后，需要根据这 N 个元素确定 G 的最小包围盒，包围盒的两个顶点分为(Gimin,Gimax)，其中i=1，2，3，其定义如图4 所示。zxyGimaxGimin图 4 几何外形包围盒Fig.4 Geometry bounding box步骤 2：创建一个空的 KDT 备用。步骤 3：若 N 为 1，则将此元
22、素插入到当前 KDT的节点中去，同时，将此节点的左右子树均置空。步骤 4：若 N 为 2，则将第一个元素插入到当前KDT 的节点中去，同时规定第一维为划分维。此时，只有第二个元素尚未被插入到 KDT 中。将第二个元素中心点与第一个元素的中心点在第一维的值进行比较。若第二个元素大于第一个元素，则对当前节点的右子树进行步骤 3 的操作，同时当前节点的左子树置空；否则，对当前节点的左子树进行步骤 3 的操作，同时当前节点的右子树置空。此步骤中 d 加 1。步骤 5：若 N 大于 2，则首先确定划分维。为了保证 KDT 的平衡，划分维的确定是计算 N 个三角形的中心点 Tcenter在各个维度的方差，
23、然后找出方差最大的维度即为划分维。然后将此 N 个元素中心点在方差最大的维度上按从小到大的顺序进行排序。将位于中间的元素插入到当前 KDT 的节点中去。小于此元素的将全部位于当前节点的左子树，大于此元素的将全部位于当前节点的右子树。对当前节点左子树和右子树的元素递归执行步骤 3 或步骤 4 或步骤 5，直至所有元素全部插入到 KDT 中。值得注意的是，在向 KDT 中插入节点的同时，节点对应的包围盒也在进行同步划分。例如 KDT 根节点 对 应 的 包 围 盒 即 为 数 据 集 G 的 包 围 盒(Gimin,Gimax)。假设节点 K 对应的包围盒为(Kimin,Kimax)，则K 的左子
24、节点 KL 的包围盒(KLimin,KLimax)为：KLimin=Kimin,KLimax=Kimax,i,uKiu,i=u(3)K 的右子节点 KR 的包围盒(KRimin,KRimax)为：KRimin=Kimin,i,uKiu,i=u,KRimax=Kimax(4)由于节点包围盒表示此空间区域内所有三角形元素集合的最小包围盒，而在进行空间区域剖分时的依据是三角形中心点确定的平面，所以节点包围盒与包围盒之间，剖分面与剖分面之间会存在相交的部分。这种相交只会出现在三维及以上的情况。2.3 创建 KDT 时间复杂度分析由于 KDT 是一个二叉树，所以对于给定的 N 个元素建成一个深度为 d
25、的最优 KDT 有：N=2d+2d1+20=(22d1)(5)d=log2(N+1)/2)Nlog2N则，这里的“”表示不小于其内部数值的最小整数。KDT 建立过程中需不断进行求中值和排序操作，且当树深为 d 时，不需对仅剩的一个元素进行求中值和排序等操作。对 N 个元素创建 KDT，时间复杂度为，这其中包括了求中值、排序以及构造 N 个节点的操作。若对相同的 N 个元素创建一个 ADT20，只需在三个维度交替平分搜索空间，每次划分后的两个子空间插入当前 ADT 的两个子节点，递归循环直至所有物面点插入 ADT 即可，不需要求中值和排序操作，仅需进行构造 N 个节点的操作。因此对 N 个元素创
26、建 KDT 的时间复杂度要大于 ADT。但物面网格一般是不变的，一次计算只需建立一次KDT 或者ADT，因此相对于ADT 来说，KDT牺牲少许建树时间获取的计算时间收益是可观的。3 最小距离查询计算空间某点到物面的最小距离需要两个步骤。第一，得到此点到物面最小距离的点所在的物面三角形；第二，采用点到三角形最小距离精确算法得到此点到物面的最小距离。本节基于 KDT 最近邻算法对第一步进行设计。96空气动力学学报第41卷 3.1 算法设计1）首先从 KDT 的根节点开始查询，根据目标点与根节点剖分面的位置关系，确定与目标点 T 潜在的最近点是在根节点的左子树还是右子树。2）根据 1）确定的子树，对
27、此子树的节点进行步骤 1）的操作。3）当递归访问至叶子节点时，记录 T 与当前节点之间的距离作为当前最近距离 lmin，当前节点作为与目标点 T 最近的节点记为 Knearest。4）进行回溯操作。回溯操作从当前最近的点Knearest开 始 向 上 进 一 步 查 找。构 造 以 T 为 球 心，lmin为半径的球。若球面与 Knearest父节点的剖分面相交，则查找 Knearest的兄弟节点即 Knearest父节点的另外子树的节点，同时计算 T 与 Knearest父节点之间的距离，若小于当前最小距离，则 Knearest和 lmin均进行更新；若球面与剖分面不相交，则说明 Knear
28、est的兄弟节点中不存在比 Knearest与 T 更近的节点。5）执行步骤 4）直至根节点。得到的 Knearest即为与 T 最近的点，lmin为 Knearest与 T 的最小距离。至此最小距离查询模块结束。3.2 算法优化abab在实际流场计算中，存在空间网格点与物面距离较远的情况，此时采用上述最小距离查询方法将面临回溯操作效率低下的问题。为了直观说明问题，本小节以二维数据点为例进行说明，给定随机分布的17 个点，对这些点创建 KDT，如图 5 所示。给定目标单元 T，正常查找得到的最邻近单元是 G。在回溯操作中，以 T 为圆心，T 和 G 之间的距离 lTG为半径的圆与根节点 J 的
29、剖分轴 m 相交。根据上述回溯步骤，此时需在 J 的左子树进行查询。然而，实际情况是圆与J 的左子树并不相交。因此，后续的回溯操作均为无效回溯，这将大大浪费计算资源。为此，本文在回溯操作时，首先判断以 T 为圆心，T 和 G 之间的距离lTG为半径的圆是否与线段相交（由剖分轴和当前包围盒确定）。若相交，则进行正常的回溯操作；若不相交，则直接排除对应节点另外子树存在最小距离的可能性。采用改进前后的回溯法进行的回溯操作如图 5（bc）中绿色箭头所示。从图中可以看出，改进后的回溯方法可极大减少不必要的回溯过程，进而提高最小距离查询效率。图 6 为球、导弹和 DPW6 三种几何外形及空间笛卡尔网格。为
30、了定量对比改进前后回溯算法的效率，本文对这三种不同几何外形进行了壁面距离查询测试。对几何外形进行壁面距离查找时，存储所有空间网格点的壁面距离，同时记录得到壁面距离所需的查找次数，然后将查找次数的总和除以空间网格点数得到平均查找次数，结果如图 7 所示。从图中可以看出，针对不同外形，本文改进的回溯方法均可有效减小查找次数。针对同一几何外形，物面三角形数量越多，查找次数减少的越明显。4 算例测试本小节对本文发展的壁面距离算法的准确性和效率进行了测试。测试的几何外形仍为球、导弹和DPW6。为了保证测试结果的合理性，本文所有工况ALOQP NMKGTmlTGJHFIBCDEabJONPQIKMAHGC
31、ELFDBJONPQIKMAHGCELFDB(a)目标点与 KDT 节点之间位置关系(b)算法改进前最小距离查询路径(c)算法改进后最小距离查询路径图 5 圆与剖分轴相交且回溯无效的情况Fig.5 Sketch of an optimized backtracking process when thecircle intersects the split axis第7期孟爽等：笛卡尔网格下精确高效的壁面距离计算方法97均在 AMDEPYC7702 的处理器上测试，测试平台为Window10 系统及 VisualStudio2012 版本。每个工况测试 10 次，统计平均获得程序运行的时间。首先
32、对壁面距离结果的准确性进行对比测试。对壁面距离的两种计算方式进行对比，第一种是近似壁面距离，即空间点到物面三角形中心的最小距离，第二种是精确壁面距离，即采用本文发展的算法计算得到的最小距离。计算结果如图 8 所示，图 8（a）为近似壁面距离计算的结果，图 8（b）为精确壁面距离计算所得结果。从图 8（a）中可以明显看出，球的近似壁面距离云图及等值线呈现波浪状，与实际物理特征明显不符。采用本文方法计算得出的精确壁面距离等值线光顺平滑，符合物理特征。对于几何外形较为复杂的 DPW6，在远离壁面区域，两种计算方法得出的壁面距离结果接近；而在近壁区域，采用近似壁面距离计算方法得到的壁面距离为 250
33、的等值线与实际几何特征误差较大，而精确计算结果的近壁面等值线可较好反映近壁物面特征。(a)球(b)导弹(c)DPW6图 6 三种几何外形Fig.6 Three geometries02.04.06.08.020304050607080查找次数三角形单元数量/105球-常规回溯法球-改进回溯法导弹-常规回溯法导弹-改进回溯法DPW6-常规回溯法DPW6-改进回溯法图 7 求最小距离所需查找次数Fig.7 Number of queries for finding the minimum distance2000180016001400120010008006004002000zxyzxyzxyW
34、all distanceWall distanceWall distance111098765432101000750500250700006000050000400003000020000100000(a)近似壁面距离结果98空气动力学学报第41卷此外，以球为研究对象，我们将本文方法计算得到的壁面距离与解析解进行对比。结果表明本文方法的计算结果与解析解差异在百万分之一以内。因此，本文方法得到的壁面距离精确性较高。在计算效率测试方面，对遍历法+精确距离、ADT+精确距离、KDT+精确距离，KDT+近似距离 4 种方法的计算效率进行对比。物面三角形数量、空间网格数量以及 CPU 运行时间如表 1
35、 所示。从表中可以看出，采用 KDT 方法的计算效率相比遍历法和 ADT 方法均有量级的提升。尤其对于复杂的 DPW6 来说，采用遍历法的时间为 3616.87s，KDT 只需 2.92s，效率提升超过 1000 倍。而对于简单的球来说，效率提升仅不到 300 倍。因此，壁面距离查询效率的提升程度与几何外形密切相关。精确距离的计算要比近似距离的计算多调用一个点到三角形距离的函数，此函数在回溯过程中频繁调用会增加计算耗时。从定量结果来看，精确距离计算耗时为近似距离计算耗时的 2 倍左右，但这对整个计算周期的耗时影响并不明显。因此，通过增加少量计算时间来获取计算精度收益是值得的。为了进一步测试本文
36、发展的方法在大规模网格计算时的效率，我们选用 DPW6 外形进行物面网格离散并生成空间笛卡尔网格。物面三角形网格数量跨越万、十万以及百万 3 个量级；空间笛卡尔网格数量跨越千万、亿以及十亿 3 个量级。效率测试结果如表 2 所示，结果表明，对于同一种物面离散方式，CPU运行时间与空间网格数量大致呈正比；对于同一数量空间网格（由于物面网格数量不同，加密相同层级生成的空间笛卡尔网格数量有略微差异），物面网格数量从 30021 增加到 1417181 时，CPU 运行时间仅增加为原来的 12 倍，这再次说明了本文基于 KDT 的回溯方法的高效性。为了与文献 16 中计算效率进行对比，本文选用 JSM
37、（JAXAstandardmodel）标模进行壁面距离查询测试，网格规模为 33 亿。本文方法（无并行）进行壁面距离查询所需时间大致为 0.51h，文献 16 中采用zxyzxyzxy2000180016001400120010008006004002000Wall distanceWall distance111098765432101000750500250Wall distance700006000050000400003000020000100000(b)精确壁面距离结果图 8 不同几何外形近似壁面距离与精确壁面计算结果云图及等值线Fig.8 Contours of approxima
38、te and accurate calculation resultsfor wall distances of different geometric shapes表 1 壁面距离查询所需时间Table 1 Wall-distance querying time几何三角形数量空间网格数量CPU运行时间/s遍历法+精确距离ADT+精确距离KDT+精确距离KDT+近似距离球4926194175722778.72130.309.865.08导弹3352258658335840.5994.837.672.66DPW63002142054623616.87140.492.921.89第7期孟爽等：笛卡
39、尔网格下精确高效的壁面距离计算方法99并行方法对进行壁面距离查询所需时间为 0.42h。因此，除去网格类型、测试环境带来的误差，总体来说，本文方法在大规模网格下仍表现出较高效率。表 2 DPW6 壁面距离计算效率Table 2 Wall-distance calculation time for DPW6三角形数量空间网格数量CPU运行时间/s300216855076437.64274828969147.411100690083588.371764936847963746.78274551797172.751099579946692.5514171816844956580.8327439357
40、6229.161098946203794.95 5 结论在本文工作中，我们基于笛卡尔网格发展了一种壁面距离计算方法，具备精确性、高效性以及通用性。基于本文内容，主要有以下 3 个方面结论：1）在提高壁面距离计算的准确性方面，我们将物面三角形进行参数化（s 和 t）表示，根据 s 和 t 的取值范围将 s-t 平面分为 7 个区域，对不同区域进行分情况讨论，从而将求最小距离问题巧妙地转化为求约束条件下一维极值问题，计算量大幅减少。2）在提高计算效率方面，采用平衡的 KDT 存储物面三角形所有顶点信息并建立相应的嵌套包围盒，使得在最小距离查询时尽可能快速排除在目标范围之外的单元；同时优化了当前超球
41、面与剖分平面的位置判断法则，解决了在距离物面较远处因无效回溯次数过多导致查询效率下降的问题。3）此外，本文的壁面距离查询模块是一个单独的模块，在计算壁面距离时，仅需知道空间网格点的位置信息以及物面三角形信息。因此，本文发展的壁面距离计算方法不仅对笛卡尔网格，对结构网格、非结构网格以及重叠网格等也可较好兼容，可方便程序的移植。本文采用 3 个三维几何外形对本文方法进行测试，结果表明计算得到的壁面距离与解析值的误差在百万分之一以内，有较高的精确性；十亿量级网格规模下的单核计算效率与已有文献 16 的并行计算效率相当。虽然在大规模网格测试中表现出较高的计算效率，但仍有提升空间。下一步我们计划添加Op
42、enMP/MPI 并行方法进一步提升计算效率，以满足下一代超大规模 CFD 湍流模拟的需求。参 考 文 献：CAPIZZANO F.Automatic generation of locally refined Cartesianmeshes:datamanagementandalgorithmsJ.InternationalJournalforNumericalMethodsinEngineering,2018,113(5):789813.doi:10.1002/nme.56361XINJJ,CHENZL,SHIF,etal.Aradialbasisfunction-basedghostce
43、llmethod for complex rigid or flexible moving boundary flowsJ.InternationalJournalofComputationalMethods,2021,18(1):2050025.doi:10.1142/s02198762205002552HONG S,YOON D,HA S,et al.A ghost-cell immersed boundarymethodforunifiedsimulationsofflowoverfinite-andzero-thicknessmovingbodiesatlargeCFLnumbersJ
44、.EngineeringApplicationsofComputationalFluidMechanics,2021,15(1):437461.doi:10.1080/19942060.2021.18809713PETON N,LARDJANE N.An immersed boundary method forgeometrical shock dynamicsJ.Journal of Computational Physics,2020,417:109573.doi:10.1016/j.jcp.2020.1095734AL-MAROUF M,SAMTANEY R.A versatile em
45、bedded boundaryadaptivemeshmethodforcompressibleflowincomplexgeometryJ.JournalofComputationalPhysics,2017,337:339378.doi:10.1016/j.jcp.2017.02.0445SESHADRIPK,DEA.Anovelsharpinterfaceimmersedboundaryframework for viscous flow simulations at arbitrary Mach numberinvolving complex and moving boundaries
46、J.Computers&Fluids,2020,206:104579.doi:10.1016/pfluid.2020.1045796PRONS,BENOITC.Automaticoff-bodyoversetadaptiveCartesianmeshmethodbasedonanoctreeapproachJ.JournalofComputationalPhysics,2013,232(1):153173.doi:10.1016/j.jcp.2012.07.0297GRIMBERGS,FARHATC.FastcomputationofthewalldistanceinunsteadyEuler
47、ianfluid-structurecomputationsJ.InternationalJournalforNumericalMethodsinFluids,2019,89(4-5):143161.doi:10.1002/fld.46868TUCKERPG.Differentialequation-basedwalldistancecomputationforDESandRANSJ.JournalofComputationalPhysics,2003,190(1):229248.doi:10.1016/S0021-9991(03)00272-99XUJL,YANC,FANJJ.Computa
48、tionsofwalldistancesbysolvingatransport equationJ.Applied Mathematics and Mechanics,2011,32(2):141150.doi:10.1007/s10483-011-1401-810TUCKERPG,RUMSEYCL,SPALARTPR,etal.ComputationsofwalldistancesbasedondifferentialequationsJ.AIAAJournal,2005,43(3):539549.doi:10.2514/1.862611TUCKER P G.Assessment of ge
49、ometric multilevel convergencerobustnessandawalldistancemethodforflowswithmultipleinternalboundariesJ.Applied Mathematical Modelling,1998,22(4-5):293311.doi:10.1016/S0307-904X(98)10007-012王刚,曾铮,叶正寅.混合非结构网格下壁面最短距离的快速计算方法J.西北工业大学学报,2014,32(4):511516.WANG G,ZENG Z,YE Z Y.An efficient search algorithm f
50、or13100空气动力学学报第41卷calculatingminimumwalldistanceofunstructuredmeshJ.JournalofNorthwestern Polytechnical University,2014,32(4):511516(inChinese).doi:10.3969/j.issn.1000-2758.2014.04.007郭中州,何志强,夏陈超,等.高效计算网格壁面距离的KD树方法J.国防科技大学学报,2017,39(4):2125.GUOZZ,HEZQ,XIACC,etal.KDtreemethodforefficientwalldistance
