多面体网格原位分布式并行边界层物理量获取方法和装置与流程

1.本发明涉及计算流体力学等多物理场数值模拟的工程应用的技术领域，尤其是涉及一种多面体网格原位分布式并行边界层物理量获取方法和装置。


背景技术：

2.流体力学中环绕物体的流动区域可分为两个区域进行处理：一个是靠近物面薄的粘性层，称为边界层，边界层内粘性作用非常重要，往往具有强烈的剪切现象；二是边界层以外的区域，一般可忽略粘性并按照无粘流动处理。prandtl基于上述假设提出的边界层理论已被称为现代流体力学的基石之一。边界层转捩表征边界层从稳定的层流状态发展为紊乱的湍流状态的连续物理过程，是现代流体力学研究中的前沿问题。
3.当前计算流体力学中的边界层转捩预测方法如转捩准则方法、转捩模型方法等，需要获取壁面附近的边界层内物理量的值。而主流方法的计算效率和计算消耗是不可接受的，只能采用间隔较多迭代步数下进行采样计算边界层物理量降低计算消耗，而不是实时更新边界层物理量，这极大限制了转捩模式、转捩准则方法在大规模网格下边界层转捩预测中的应用，导致边界层转捩预测精度降低，计算成本变大，边界层转捩相关工程设计周期变长。


技术实现要素：

4.本发明的目的在于提供一种多面体网格原位分布式并行边界层物理量获取方法和装置，解决基于多面体网格的计算流体力学方法中，难以高效率、低内存地获取边界层参数，导致所构造的转捩模式和转捩准则方法应用受限、计算成本较高等技术问题。
5.第一方面，本发明提供了一种多面体网格原位分布式并行边界层物理量获取方法，包括：将每个处理器的壁面网格单元广播到其余处理器中，以使所述每个处理器基于直线距离建立壁面网格单元的k-d树数据结构，并进行邻近点搜索；基于所述邻近点确定出的网格体单元的最近距离壁面网格单元，将每个处理器中的每个网格单元对应的网格体单元编号和重叠区网格体单元，按照法向距离插入法向网格体单元链表中；每个所述处理器根据物理场解算器依次求解无粘欧拉(euler)方程与粘性纳维斯托克斯（navier-stokes）方程，获取流场物理解并基于所述流场物理解更新区域分解重叠区域的物理场信息；各个所述处理器根据所述区域分解重叠区域的物理场信息计算当地处理器内的边界层物理量特征量，通过并行计算函数全局规约计算全局的壁面网格单元边界层特征物理量。
6.结合第一方面，本发明提供了第一方面的第一种可能的实施方式，其中，在将每个处理器的壁面网格单元广播到其余处理器中，以使所述每个处理器基于直线距离建立壁面
网格单元的k-d树数据结构，并进行邻近点搜索的步骤之前，所述方法还包括：初始化并行计算环境，对输入的网格进行分布式并行区域分解，创建多层重叠的并行计算通信数据结构。
7.结合第一方面，本发明提供了第一方面的第二种可能的实施方式，其中，将每个处理器的壁面网格单元广播到其余处理器中，以使所述每个处理器基于直线距离建立壁面网格单元的k-d树数据结构，并进行邻近点搜索的步骤，包括：每个处理器将收集到壁面网格单元信息广播到其他处理器；其中，所述壁面网格单元信息包括：壁面网格单元的法向量、壁面网格面心坐标、组成壁面网格的网格点坐标和全局网格面编号；根据壁面在区域分解后所处的处理器位置，对全局壁面网格单元赋予局部编号；基于直线距离，以所述全局壁面网格单元的面心坐标为数据，建立用于表征壁面网格单元的k-d树数据结构；基于所述k-d树数据结构，输入点的最邻近搜索，确定所述壁面网格单元中与所述输入点最接近的点。
8.结合第一方面，本发明提供了第一方面的第三种可能的实施方式，其中，基于所述邻近点确定出的网格体单元的最近距离壁面网格单元，将每个处理器中的每个网格单元对应的网格体单元编号和重叠区网格体单元，按照法向距离插入法向网格体单元链表中的步骤，包括：对每个处理器中的每个网格单元重复执行以下步骤，直至每个所述处理器中的每个所述网格单元均被遍历：查找每个所述网格体单元对应的最近距离壁面网格单元，记录最近距离壁面网格单元的网格面编号，以及所述网格体单元和所述最近距离壁面网格单元之间的直线距离；根据所述直线距离，求解所述网格体单元到所述壁面网格单元的法向距离，并将该网格体单元编号按法向距离大小排列插入到所述壁面网格单元的法向网格体单元链表中；启动通信接口，交换区域分解重叠区网格体单元的距离信息，并将所述重叠区网格体单元按法向距离大小排列插入到与所述重叠区网格体单元对应的目标壁面网格单元的法向网格体单元链表中。
9.结合第一方面，本发明提供了第一方面的第四种可能的实施方式，其中，每个所述处理器根据物理场解算器依次求解无粘欧拉方程与粘性纳维斯托克斯方程，获取流场物理解并基于所述流场物理解更新区域分解重叠区域的物理场信息的步骤，包括：每个处理器计算当地网格单元的无粘欧拉方程物理场；每个处理器计算当地网格单元的粘性纳维斯托克斯方程物理场；根据预先建立多层重叠的并行计算通信数据结构，确定接口交换区域分解的重叠区网格体单元的无粘以及粘性物理场变量；基于所述接口交换区域分解的重叠区网格体单元的无粘以及粘性物理场变量，更新区域分解重叠区域的物理场信息。
10.结合第一方面，本发明提供了第一方面的第五种可能的实施方式，其中，各个所述处理器根据所述区域分解重叠区域的物理场信息计算当地处理器内的边界层物理量特征
量，通过并行计算函数全局规约计算全局的壁面网格单元边界层特征物理量的步骤，包括：对各个所述处理器中每个全局壁面网格单元重复执行以下步骤，直至每个所述处理器中每个全局壁面网格单元均被遍历：根据所述壁面网格单元的网格体单元链表，收集流动物理量中的无粘及粘性物理量；根据壁面网格单元的法向网格体单元的无粘和粘性物理量，计算当地边界层特征量；针对所述当地边界层特征量，确定每个所述壁面网格单元的特征点；利用并行计算函数全局规约每个所述壁面网格单元的第一特征点，以及每个所述壁面网格单元的第一特征量与法向距离，并规约每个所述壁面网格单元的第二特征点以及每个所述壁面网格单元的第二特征量与法向距离。
11.结合第一方面，本发明提供了第一方面的第六种可能的实施方式，其中，所述方法还包括：根据各个所述处理器计算出的全局的壁面网格单元边界层特征物理量以及当地处理器对应的物理场信息，计算出所述壁面网格单元的其它边界层物理量；根据所述其它边界层物理量计算目标物理场求解值。
12.第二方面，本发明还提供一种多面体网格原位分布式并行边界层物理量获取装置，包括：搜索模块，将每个处理器的壁面网格单元广播到其余处理器中，以使所述每个处理器基于直线距离建立壁面网格单元的k-d树数据结构，并进行邻近点搜索；插表模块，基于所述邻近点确定出的网格体单元的最近距离壁面网格单元，将每个处理器中的每个网格单元对应的网格体单元编号和重叠区网格体单元，按照法向距离插入法向网格体单元链表中；获取模块，每个所述处理器根据物理场解算器依次求解无粘欧拉方程与粘性纳维斯托克斯方程，获取流场物理解并基于所述流场物理解更新区域分解重叠区域的物理场信息；计算模块，各个所述处理器根据所述区域分解重叠区域的物理场信息计算当地处理器内的边界层物理量特征量，通过并行计算函数全局规约计算全局的壁面网格单元边界层特征物理量。
13.第三方面，本发明提供一种电子设备，包括存储器、处理器，所述存储器中存储有可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现前述实施方式任一项所述的方法的步骤。
14.第四方面，本发明提供一种机器可读存储介质，所述机器可读存储介质存储有机器可执行指令，所述机器可执行指令在被处理器调用和执行时，机器可执行指令促使处理器实现前述实施方式任一项所述的方法的步骤。
15.本发明实施例带来了以下有益效果：1 .通过本方法可原位分布式并行计算壁面网格单元对应的边界层物理量，无需移动复制全场数据，计算效率高，通信次数及通信量少，内存占用少。
16.2 .本方法根据边界层理论，使用无粘流场和粘性流场计算出判定边界层的特征
量1（第一特征量），该特征物理量计算简单可靠，可有效降低了复杂外形下的特征物理量判定难度。
17.3 .本方法使用调优的k-d树数据结构处理最近点搜索问题，所建立的k-d树效率高，作为独立模块嵌入到现有程序中。
18.本发明的其他特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点在说明书以及附图中所特别指出的结构来实现和获得。
19.为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。
附图说明
20.为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
21.图1为是本发明所述的多面体网格原位分布式并行边界层物理量获取方法的步骤示意图；图2为本发明的输入接口对应的二维翼型多面体网格计算区域示意图；图3为本发明所对应的图剖分无重叠区域分解及人工边界示意图；图4为本发明所形成的三层重叠区域分解并行计算区域示意图；图5为本发明各处理器获取所有壁面网格单元示意图；图6为本发明构造k-d树插入根节点后示意图；图7为本发明所求解的壁面法向距离示意图；图8为本发明所构造的当地壁面法向单元列表示意图；图9为本发明计算特征量所使用的局部分解坐标系示意图；图10为本发明计算特征量所使用的当地-全局归并示意图；图11为本发明计算所得的边界层厚度的示意图；图12为本发明提供的一种多面体网格原位分布式并行边界层物理量获取装置的功能模块示意图；图13为本发明提供的电子设备的硬件架构示意图。
具体实施方式
22.为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
23.边界层转捩是一种十分常见的流动现象，对航空飞行器、涡轮机械及航天飞行器等不同领域中的各类流动有重要影响。对于航空飞行器流动，飞行器湍流流动下的摩擦阻力一般为层流状态下的数倍，所以期望通过层流化提高航行经济性，但在某些情况下又希
望边界层尽快转捩，进而抑制边界层流动分离。对于航空涡轮机械流动，边界层转捩直接影响总压恢复系数、推进效率等发动机性能参数。对于航天飞行器，高速飞行的气动加热问题十分突出重要，边界层转捩过程伴随气动加热的剧烈变化，预测边界层转捩有益于航天飞行器的热防护设计工作。
24.计算流体力学中的边界层转捩预测方法如转捩准则方法、转捩模型方法等，需要获取壁面附近的边界层内物理量的值。例如在可压缩流动的转捩预测中，边界层的速度剖面具有重要作用，边界层内出现的垂直于边界层外无粘流动的速度分量，称之为横流。横流由相对于边界层外的无粘流速度偏转角和流动速度综合决定，具体表现为：在边界层外，由于偏转角为零，横流速度为零；在壁面处由于流动速度为零，横流速度也为零，边界层内的横流速度存在拐点。由于速度拐点的存在，横流流动具有不稳定性，在一定外界扰动的条件下，容易诱发边界层转捩。由此可见，边界层内物理量对边界层转捩预测意义极大。目前的计算流体力学方法中，大多基于某种网格划分进行分布式计算，壁面附近的边界层网格单元可能并不在当地数据结构中，在原位数据结构上求解边界层十分困难，因此边界层内物理量采用的主要获取方法是重构一套边界层法向网格数据结构（或称为法向积分路径），对整个计算区域的物理量进行复制并分发到新重构的数据结构中，然后在新建立的数据结构中进行边界层物理量计算。这个过程不是原位并行的——即没有使用物理场求解器的原始分布式并行数据结构，而是重构了新的专门用于边界层物理量计算的全局数据结构。由于不是原位并行的，整个过程需要对整体流场进行收集复制和分发计算，计算后需要再次集中归并到原始数据结构中。
25.如上所述，这个计算过程具有两次必须的全局串行操作，由于目前计算流体力学的典型计算网格在几千万到几亿的量级，该操作的通信量、内存占用均极大。实际计算中，往往需要频繁获取边界层物理量，目前主流方法的计算效率和计算消耗是不可接受的，只能采用间隔较多迭代步数下进行采样计算边界层物理量降低计算消耗，而不是实时更新边界层物理量，这极大限制了转捩模式、转捩准则方法在大规模网格下边界层转捩预测中的应用，导致边界层转捩预测精度降低，计算成本变大，边界层转捩相关工程设计周期变长。
26.基于此，本发明实施例提供的一种多面体网格原位分布式并行边界层物理量获取方法，可以高效率、低内存地获取边界层参数，使所构造的转捩模式和转捩准则方法能够广泛应用，降低计算成本。
27.为便于对本实施例进行理解，首先对本发明实施例所公开的一种多面体网格原位分布式并行边界层物理量获取方法进行详细介绍，该方法可应用于服务器、上位机等智能控制设备中。
28.图1为本发明实施例提供的一种多面体网格原位分布式并行边界层物理量获取方法流程图。
29.参照图1，该方法主要包括以下步骤：步骤s102，将每个处理器的壁面网格单元广播到其余处理器中，以使每个处理器基于直线距离建立壁面网格单元的k-d树数据结构，并进行邻近点搜索。
30.步骤s104，基于邻近点确定出的网格体单元的最近距离壁面网格单元，将每个处理器中的每个网格单元对应的网格体单元编号和重叠区网格体单元，按照法向距离插入法向网格体单元链表中。
31.步骤s106，每个处理器根据物理场解算器依次求解无粘欧拉方程与粘性纳维斯托克斯方程，获取流场物理解并基于流场物理解更新区域分解重叠区域的物理场信息。
32.步骤s108，各个处理器根据区域分解重叠区域的物理场信息计算当地处理器内的边界层物理量特征量，通过并行计算函数全局规约计算全局的壁面网格单元边界层特征物理量。
33.在步骤s102之前，该方法还包括步骤s101：初始化并行计算环境，对输入的网格进行分布式并行区域分解，创建多层重叠的并行计算通信数据结构。
34.其中，初始化后消息传递接口（messagepassinginterface，mpi）环境，输入计算流体力学的计算区域网格，使用一种针对多面体网格的多层区域分解并行计算方法及系统对网格进行预处理，对输入的网格进行分布式并行区域分解，创建多层重叠的并行计算通信数据结构，考虑到边界层物理量计算的需要，区域分解重叠的层次为典型的三层；如图2所示，为二维翼型的多面体网格计算区域示意图，其中粗线为翼型壁面；图3为初次区域分解后的无重叠计算区域示意图，计算区域被分解为两部分，分配到两个处理器中，其壁面部分也被分解成两部分；图4为创建多层重叠的区域分级并行计算通信数据结构后的各处理器网格示意图，图中的阴影部分为从其他处理器中移动的重叠区域网格单元，各处理器通过通用接口交换阴影部分的物理场数据，如图4中箭头所示。
35.在一些实施例中，各处理器可广播本处理器的壁面网格单元到所有处理器，每个处理器均获取了计算区域的全部壁面网格单元列表，对全部壁面网格单元列表内所有的网格面单元建立基于直线距离的k-d树（k-dimensional树）数据结构，用于单元的最近点搜索，作为一种示例，该步骤s102，可包括：步骤1.1），每个处理器将收集到壁面网格单元信息广播到其他处理器；其中，壁面网格单元信息包括：壁面网格单元的法向量、壁面网格面心坐标、组成壁面网格的网格点坐标和全局网格面编号；其中，各处理器收集当地的壁面网格单元的法向量、壁面网格面心坐标、组成壁面网格的网格点坐标和全局网格面编号，各处理器将收集的壁面网格单元信息广播到所有处理器的全局壁面网格单元数据存储向量中；步骤1.2），根据壁面在区域分解后所处的处理器位置，对全局壁面网格单元赋予局部编号；具体地，各个处理器单元对全局壁面网格单元赋予局部编号，局部编号顺序按照各处理器的壁面网格顺序排列依次递增，先给处理器0编号，然后在处理器0基础上给处理器1编号；由于局部编号按处理器排列，可根据局部编号快速确定该壁面在区域分解后实际的处理器位置；图5为各处理器收集壁面信息并编号后的示意图，从图5中可以看出，每个处理器中的同一壁面具有相同的编号，如图5左右两侧处理器的相同位置壁面编号相同，处理器0和处理器1的最左端壁面均编号f10，最右端壁面均编号f34；同一处理器中的壁面网格单元具有连续编号并随着处理器的编号递增，可以理解的是，壁面网格单元从f0到f49进行连续递增编号，图5中仅示出若干位置所对应的壁面网格单元编号；步骤1.3），基于直线距离，以全局壁面网格单元的面心坐标为数据，建立用于表征壁面网格单元的k-d树数据结构；基于k-d树数据结构，输入点的最邻近搜索，确定壁面网格单元中与输入点最接近的点。
36.其中，以全局壁面网格单元的面心坐标为数据建立k-d树，k-d树是一种每个树节点都为k维点的二叉树，在本发明的方法中，k等于3。建立k-d树的过程为逐点插入编号为wi的壁面网格单元的面心坐标，使用壁面网格单元wi的面心坐标x轴作为分割超平面，把空间分割成两个半空间。在分割超平面左边的点位于节点左边的子树，在分割超平面右边的点位于节点右边的子树。图6为插入头节点a之后的空间分割示意图，整个空间被a点的x轴分割为两个区域。
37.随后依次用下一个插入壁面网格单元的面心坐标y轴，面心坐标z轴作为分割超平面建立子树。该k-d树建立完成后将被用于输入点的最邻近搜索，用来找出在k-d树中（即壁面网格单元中）与输入点最接近的点。k-d树建立过程中的关键参数为最大叶子节点数量maximum_leaf_number，距离函数容差dist_eps，网格单元包围盒尺寸bounding_box。最大叶子节点数量maximum_leaf_number决定了k-d树的高度和宽度，由于搜索算法需要遍历树的高度，然后回溯，因此该参数极大影响搜索算法的效率，在本发明的方法中，其最优值为30≤maximum_leaf_number≤50。距离函数容差dist_eps影响最近点查找的精度，计算方法如下：根据壁面网格单元的构成网格点坐标计算所有壁面网格单元的边长，所有壁面网格单元的最小边长即为resolution，利用dist_eps=resolution*beta得到距离函数容差dist_eps，其中1e-8≤beta≤1e-6，beta是一个可调参数，表示对长度缩放一定的量级。网格单元包围盒尺寸bounding_box用于快速定位可能最近点，是搜索算法高效率的关键，包括coord_min[wi][j]和coord_max[wi][j]，其中，wi为壁面网格单元编号，j为坐标轴序号；coord_min为壁面网格单元的构成网格点坐标分量下确界，coord_max为壁面网格单元的构成网格点坐标分量上确界。图6中绘制了壁面网格单元f10的网格单元包围盒尺寸bounding_box。
[0038]
在实际应用中，步骤s104可通过对当地处理器中的所有当地网格体单元进行循环，在k-d树中使用最近查找算法查找该编号为ci的网格体单元对应的编号为wj的最近距离壁面网格单元，记录该网格体单元的最近距离壁面网格单元的编号wj及其直线距离变量distance，根据直线距离求解网格体单元到壁面网格单元的法向距离distance_normal，同时将该网格体单元编号按法向距离大小排列插入到编号为wj的最近距离壁面网格单元的法向网格体单元链表wall_norm_cell_list[wj]中；启动通信接口，交换区域分解重叠区网格体单元的距离信息，并将重叠区网格体单元按法向距离大小排列插入到对应编号为wj的最近距离壁面网格单元（目标壁面网格单元）的法向网格体单元链表wall_norm_cell_list[wj]中；作为一种可选的实施例，可通过以下步骤实现，包括：对每个处理器中的每个网格单元重复执行以下步骤，直至每个处理器中的每个网格单元均被遍历：步骤2.1），查找每个网格体单元对应的最近距离壁面网格单元，记录最近距离壁面网格单元的网格面编号和网格体单元和最近距离壁面网格单元之间的直线距离；其中，对当地处理器中的所有当地网格体单元进行循环，在k-d树中使用最近查找算法查找该编号为ci的网格体单元对应的最近距离壁面网格单元，搜索结果返回的是最近距离壁面网格单元的编号wj及所计算的直线距离变量distance，将直线距离变量distance存储于同样用于表征直线距离的数组distance_cell[ci]，编号wj存储于最近网格壁面nearst_wall[ci]；其中，distance_cell[ci]表示第ci个网格体单元cell的距离是直线距
离变量distance存储的值，每个网格体单元cell都有一个值；第ci个网格体单元cell的最近距离壁面网格单元的编号是wj，每个网格体单元cell都有一个最近距离的壁面网格单元，全部最近距离的壁面网格单元所构成的最近网格壁面nearst_wall[ci]。
[0039]
步骤2.2），根据直线距离，求解网格体单元到壁面网格单元的法向距离，并将该网格体单元编号按法向距离大小排列插入到壁面网格单元的法向网格体单元链表中；需要说明的是，由于最近点算法计算的距离是编号为ci的网格体单元格心到编号为wj的最近距离壁面网格单元面心的直线距离，而不是法向距离，在近壁区网格倾斜扭曲时两者相差较大，但在远壁区取直线距离结果更光滑，因此采用下列算法进行法向距离计算：如果所计算的直线距离distance_cell[ci]大于特征距离l，法向距离distance_normal_cell[ci]直接取直线距离distance_cell[ci]的值；如果所计算的直线距离distance_cell[ci]小于等于特征距离l，法向距离distance_normal_cell[ci]等于直线距离distance_cell[ci]与编号为wj的最近距离壁面网格单元的面法向wall_norm[wj]的内积，即；图7为对网格体c的网格体格心所计算的直线距离与计算法向距离的示意图，图7展示了头部区域网格体直接采用直线距离作为法向距离的示意图；其中，图7中m代表直线尽头的圆点，即该直线所指向的网格的格心，同理c点该直线所指向网格的格心。
[0040]
步骤2.3），启动通信接口，交换区域分解重叠区网格体单元的距离信息，并将重叠区网格体单元按法向距离大小排列插入到与重叠区网格体单元对应的目标壁面网格单元的法向网格体单元链表中。
[0041]
示例性地，将网格体单元编号ci插入到编号为wj的最近距离壁面网格单元的法向网格体单元链表wall_norm_cell_list[wj]中，使用的算法为插入排序算法，按照法向距离distance_normal_cell[ci]为插入排序算法的关键字，找到第一个大于该值的链表节点，将该单元编号插入到该链表节点之前；图10展示了壁面f22所对应的法向网格体单元wall_norm_cell[f22]，处理器0包括当地的网格体单元a，b，c，处理器1包括当地的网格体单元d，e；各处理器启动前述步骤s101所形成的通信接口，交换区域分解重叠区体网格的法向距离distance_normal_cell[ci]与最近网格壁面nearst_wall[ci]值，将并将重叠区网格体单元按法向距离大小排列插入到对应编号为wj的最近距离壁面网格单元wj构成的最近网格壁面nearst_wall[ci]的法向网格体单元链表wall_norm_cell_list[wj]中；图8箭头指示器展示了壁面f22的区域分解重叠区法向网格体单元列表，处理器0交换获得了区域分解重叠区的网格体单元a，b，c；处理器1交换获得了区域分解重叠区的网格体单元d，e。
[0042]
作为一种可选的实施例，各处理器启动物理场解算器，依次求解无粘欧拉方程，与粘性纳维斯托克斯方程，获取流场物理解，并更新区域分解重叠区域的物理场信息；该步骤s106可包括：步骤3.1)，每个处理器计算当地网格单元的无粘欧拉方程物理场；步骤3.2)，每个处理器计算当地网格单元的粘性纳维斯托克斯方程物理场；步骤3.3)，根据预先建立多层重叠的并行计算通信数据结构，确定接口交换区域分解的重叠区网格体单元的无粘以及粘性物理场变量；
这里，使用步骤s101建立的接口交换区域分解的重叠区网格体单元的无粘以及粘性物理场变量；步骤3.4)，基于接口交换区域分解的重叠区网格体单元的无粘以及粘性物理场变量，更新区域分解重叠区域的物理场信息。
[0043]
在一些实施例中，各处理器根据无粘物理场与粘性物理场计算当地处理器内的边界层物理量特征量，并将当地边界层特征量数据记录到壁面网格单元的边界层特征物理量向量容器中，通过用于并行计算的mpi函数全局规约计算全局的壁面网格单元边界层特征物理量；该步骤s108主要包括：对各个处理器中每个全局壁面网格单元重复执行以下步骤，直至每个处理器中每个全局壁面网格单元均被遍历：步骤4.1），根据壁面网格单元的网格体单元链表，收集流动物理量中的无粘及粘性物理量；这里，对全局壁面网格单元wi进行循环处理，对该壁面网格单元的wall_norm_cell[wi]链表中壁面网格单元wi的法向网格体单元cj的流动物理量进行处理，收集无粘及粘性物理量信息；其中，法向网格体单元cj是法向网格体单元wall_norm_cell[wi]里的其中一个数据，法向网格体单元wall_norm_cell[wi]中包括一串网格体单元。
[0044]
步骤4.2），根据壁面网格单元的法向网格体单元的无粘和粘性物理量，计算当地边界层特征量；具体地，根据壁面网格单元wi的法向网格体单元cj的无粘和粘性物理量，计算当地边界层特征量；首先获取当地速度分解坐标系对当地网格的粘性速度分量进行分解，分解坐标系为利用壁面网格单元wi的法向网格体单元cj速度单位矢量a1_unit，壁面单位法向量b1_unit计算得到横向单位向量，利用壁面单位法向量b1_unit、横向单位向量c1_unit计算得到流向单位向量，壁面单位法向量b1_unit、横向单位向量c1_unit、流向单位向量d1_unit组成了当地的分解坐标系，如图9所示；特征量1（第一特征量）g的计算方法如下，分解出粘性速度矢量u的粘性流向速度分量u_stream=u
•
d1_unit=ux*unitx+uy*unity+uz*unitz，计算出粘性流向速度分量u_stream与无粘速度u_inviscid的比值为g =u_stream/u_inviscid为特征量1；特征量1无需边界层内其他网格的信息，具有完全当地性的特点，计算时无需与其他网格体单元发生联系；特征量2（第二特征量）h基于粘性物理场的相关物理量，具体为，其中， norm表示法向方向。d（u_stream）/dnorm表示流向速度在法向上的偏微分；ρ为粘性物理场密度值，特征值2需使用相邻的法向网格体单元，但由于本发明构造了三层重叠的并行区域分解，因此当地处理器计算该特征值时不会出现相邻网格体单元缺失的情况，当地计算无任何困难；如图6所示，计算特征量2时，使用了单元a，单元b，单元c以及壁面等一系列法向单
元的值。
[0045]
步骤4.3），针对当地边界层特征量，确定每个壁面网格单元的特征点；具体地，针对壁面网格单元wi的法向网格体单元cj的当地边界层特征量，第一特征量g=u_stream/u_inviscid，当g》= g’， g’为给定值且法向距离distance_normal_cell[cj]最小时的cj为第一特征点wall[wi]-cj_character_local_1；当第二特征量，为给定值且法向距离distance_normal_cell[cj]最小时为第二特征点wall[wi]-cj_character_local_2；步骤4.4），利用并行计算函数全局规约每个壁面网格单元的第一特征点，以及每个壁面网格单元的第一特征量（特征量1）与第一法向距离，并规约每个壁面网格单元的第二特征点以及每个壁面网格单元的第二特征量（特征量2）与第二法向距离。
[0046]
利用mpi函数全局规约所有壁面网格单元的第一特征点wall[wi]-cj_character_local_1，及其第一特征量g与第一法向距离distance_normal_cell[cj_character_local_1]，规约所有壁面网格单元的第二特征点wall[wi]-cj_character_local_2及其第二特征量h与第二法向距离distance_normal_cell[cj_character_local_2]；对于第一特征量，其全局的壁面网格单元特征点为所有的g》=g’，g’为给定值且第一法向距离distance_normal_cell[cj_character_local_1]最小时的网格体单元，记为wall[wi]-cj_character_1；对于第二特征量，其全局的壁面网格单元特征点为所有的，为给定值且第二法向距离distance_normal_cell[cj_character_local_2]最小时的网格体单元，记为wall[wi]-cj_character_2；所述步骤中的mpi函数全局规约操作主要函数为全收集函数mpi_allgather，可参见图10。
[0047]
在一些实施例中，前述实施例的方法还包括：步骤s109，根据各个处理器计算出的全局的壁面网格单元边界层特征物理量以及当地处理器对应的物理场信息，计算出壁面网格单元的其它边界层物理量；根据其它边界层物理量计算目标物理场求解值。
[0048]
可以理解的是，根据各处理器计算出的全局特征量信息以及当地物理场，计算出其它的边界层物理量信息，根据新计算的壁面网格单元边界层物理量信息计算下一步物理场求解值。
[0049]
本实例描述了应用本发明方法在针对翼身组合体实际模型时的原位分布式并行具体实施方案。翼身组合体网格大约为651万多面体网格单元，使用本发明的边界层特征量计算方法进行单核，双核，8核的分布式并行计算，其计算时间如下表1所示：表1本实施例的原位分布式计算耗时
测试外形单核边界层特征量计算时间2核边界层特征量计算时间8核边界层特征量计算时间翼身组合体(651万网格)7.10874.20987s1.86731
从表1中可以看出，本实施例的边界层特征量计算方法具有高度的并行性，计算效率高。
[0050]
本实例描述了应用本发明方法在针对平板模型时的具体实施方案。本发明方法应用于层流平板边界层流动，来流马赫数为0.3，攻角和侧滑角均为0
°
，壁面为绝热温度条件，
这个模型在工况下有动量积分理论关系式，其中，为雷诺数，该雷诺数用于标识流动状态，如图11所示的是应用本发明方法所计算的边界层厚度信息，其中横坐标x表示平板到前缘头部的距离，纵坐标θ表示边界层的厚度，本发明计算所得规律与理论上的规律一致。
[0051]
如图12所示，本发明实施例提供一种多面体网格原位分布式并行边界层物理量获取装置，包括：搜索模块，将每个处理器的壁面网格单元广播到其余处理器中，以使所述每个处理器基于直线距离建立壁面网格单元的k-d树数据结构，并进行邻近点搜索；插表模块，基于所述邻近点确定出的网格体单元的最近距离壁面网格单元，将每个处理器中的每个网格单元对应的网格体单元编号和重叠区网格体单元，按照法向距离插入法向网格体单元链表中；获取模块，每个所述处理器根据物理场解算器依次求解无粘欧拉方程与粘性纳维斯托克斯方程，获取流场物理解并基于所述流场物理解更新区域分解重叠区域的物理场信息；计算模块，各个所述处理器根据所述区域分解重叠区域的物理场信息计算当地处理器内的边界层物理量特征量，通过并行计算函数全局规约计算全局的壁面网格单元边界层特征物理量。
[0052]
在一些实施例中，所述装置还包括，创建模块，用于初始化并行计算环境，对输入的网格进行分布式并行区域分解，创建多层重叠的并行计算通信数据结构。
[0053]
在一些实施例中，搜索模块，还具体用于，每个处理器将收集到壁面网格单元信息广播到其他处理器；其中，所述壁面网格单元信息包括：壁面网格单元的法向量、壁面网格面心坐标、组成壁面网格的网格点坐标和全局网格面编号；根据壁面在区域分解后所处的处理器位置，对全局壁面网格单元赋予局部编号；基于直线距离，以所述全局壁面网格单元的面心坐标为数据，建立用于表征壁面网格单元的k-d树数据结构；基于所述k-d树数据结构，输入点的最邻近搜索，确定所述壁面网格单元中与所述输入点最接近的点。
[0054]
在一些实施例中，插表模块，还具体用于，对每个处理器中的每个网格单元重复执行以下步骤，直至每个所述处理器中的每个所述网格单元均被遍历：查找每个所述网格体单元对应的最近距离壁面网格单元，记录最近距离壁面网格单元的网格面编号，以及所述网格体单元和所述最近距离壁面网格单元之间的直线距离；根据所述直线距离，求解所述网格体单元到所述壁面网格单元的法向距离，并将该网格体单元编号按法向距离大小排列插入到所述壁面网格单元的法向网格体单元链表中；启动通信接口，交换区域分解重叠区网格体单元的距离信息，并将所述重叠区网格体单元按法向距离大小排列插入到与所述重叠区网格体单元对应的目标壁面网格单元的法向网格体单元链表中。
[0055]
在一些实施例中，获取模块，还具体用于，每个处理器计算当地网格单元的无粘欧拉方程物理场；每个处理器计算当地网格单元的粘性纳维斯托克斯方程物理场；根据预先建立多层重叠的并行计算通信数据结构，确定接口交换区域分解的重叠区网格体单元的无粘以及粘性物理场变量；基于所述接口交换区域分解的重叠区网格体单元的无粘以及粘性物理场变量，更新区域分解重叠区域的物理场信息。
[0056]
在一些实施例中，计算模块，还具体用于，对各个所述处理器中每个全局壁面网格
单元重复执行以下步骤，直至每个所述处理器中每个全局壁面网格单元均被遍历：根据所述壁面网格单元的网格体单元链表，收集流动物理量中的无粘及粘性物理量；根据壁面网格单元的法向网格体单元的无粘和粘性物理量，计算当地边界层特征量；针对所述当地边界层特征量，确定每个所述壁面网格单元的特征点；利用并行计算函数全局规约每个所述壁面网格单元的第一特征点，以及每个所述壁面网格单元的第一特征量与法向距离，并规约每个所述壁面网格单元的第二特征点以及每个所述壁面网格单元的第二特征量与法向距离。
[0057]
在一些实施例中，计算模块，还具体用于，根据各个所述处理器计算出的全局的壁面网格单元边界层特征物理量以及当地处理器对应的物理场信息，计算出所述壁面网格单元的其它边界层物理量；根据所述其它边界层物理量计算目标物理场求解值。
[0058]
本发明实施例提供的用于实现一种电子设备，本实施例中，所述电子设备可以是，但不限于，个人电脑（personal computer，pc）、笔记本电脑、监控设备、服务器等具备分析及处理能力的计算机设备。
[0059]
作为一种示范性实施例，可参见图13，电子设备110，包括通信接口111、处理器112、存储器113以及总线114，处理器112、通信接口111和存储器113通过总线114连接；上述存储器113用于存储支持处理器112执行上述方法的计算机程序，上述处理器112被配置为用于执行该存储器113中存储的程序。
[0060]
本文中提到的机器可读存储介质可以是任何电子、磁性、光学或其它物理存储装置，可以包含或存储信息，如可执行指令、数据，等等。例如，机器可读存储介质可以是：ram（radom access memory，随机存取存储器）、易失存储器、非易失性存储器、闪存、存储驱动器（如硬盘驱动器）、任何类型的存储盘（如光盘、dvd等），或者类似的存储介质，或者它们的组合。
[0061]
非易失性介质可以是非易失性存储器、闪存、存储驱动器（如硬盘驱动器）、任何类型的存储盘（如光盘、dvd等），或者类似的非易失性存储介质，或者它们的组合。
[0062]
可以理解的是，本实施例中的各功能模块的具体操作方法可参照上述方法实施例中相应步骤的详细描述，在此不再重复赘述。
[0063]
本发明实施例所提供计算机可读存储介质，所述可读存储介质中存储有计算机程序，所述计算机程序代码被执行时可实现上述任一实施例所述的方法，具体实现可参见方法实施例，在此不再赘述。
[0064]
所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的系统和装置的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。
[0065]
另外，在本发明实施例的描述中，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。
[0066]
在本发明的描述中，需要说明的是，术语“中心”、“上”、“下”、“左”、“右”、“竖直”、“水平”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、
以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性。
[0067]
最后应说明的是：以上所述实施例，仅为本发明的具体实施方式，用以说明本发明的技术方案，而非对其限制，本发明的保护范围并不局限于此，尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化，或者对其中部分技术特征进行等同替换；而这些修改、变化或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围，都应涵盖在本发明的保护范围之内。

技术特征：
1.一种多面体网格原位分布式并行边界层物理量获取方法，其特征在于，包括：将每个处理器的壁面网格单元广播到其余处理器中，以使所述每个处理器基于直线距离建立壁面网格单元的k-d树数据结构，并进行邻近点搜索；基于所述邻近点确定出的网格体单元的最近距离壁面网格单元，将每个处理器中的每个网格单元对应的网格体单元编号和重叠区网格体单元，按照法向距离插入法向网格体单元链表中；每个所述处理器根据物理场解算器依次求解无粘欧拉方程与粘性纳维斯托克斯方程，获取流场物理解并基于所述流场物理解更新区域分解重叠区域的物理场信息；各个所述处理器根据所述区域分解重叠区域的物理场信息计算当地处理器内的边界层物理量特征量，通过并行计算函数全局规约计算全局的壁面网格单元边界层特征物理量。2.根据权利要求1所述的方法，其特征在于，在将每个处理器的壁面网格单元广播到其余处理器中，以使所述每个处理器基于直线距离建立壁面网格单元的k-d树数据结构，并进行邻近点搜索的步骤之前，所述方法还包括：初始化并行计算环境，对输入的网格进行分布式并行区域分解，创建多层重叠的并行计算通信数据结构。3.根据权利要求2所述的方法，其特征在于，将每个处理器的壁面网格单元广播到其余处理器中，以使所述每个处理器基于直线距离建立壁面网格单元的k-d树数据结构，并进行邻近点搜索的步骤，包括：每个处理器将收集到壁面网格单元信息广播到其他处理器；其中，所述壁面网格单元信息包括：壁面网格单元的法向量、壁面网格面心坐标、组成壁面网格的网格点坐标和全局网格面编号；根据壁面在区域分解后所处的处理器位置，对全局壁面网格单元赋予局部编号；基于直线距离，以所述全局壁面网格单元的面心坐标为数据，建立用于表征壁面网格单元的k-d树数据结构；基于所述k-d树数据结构，输入点的最邻近搜索，确定所述壁面网格单元中与所述输入点最接近的点。4.根据权利要求2所述的方法，其特征在于，基于所述邻近点确定出的网格体单元的最近距离壁面网格单元，将每个处理器中的每个网格单元对应的网格体单元编号和重叠区网格体单元，按照法向距离插入法向网格体单元链表中的步骤，包括：对每个处理器中的每个网格单元重复执行以下步骤，直至每个所述处理器中的每个所述网格单元均被遍历：查找每个所述网格体单元对应的最近距离壁面网格单元，记录最近距离壁面网格单元的网格面编号，以及所述网格体单元和所述最近距离壁面网格单元之间的直线距离；根据所述直线距离，求解所述网格体单元到所述壁面网格单元的法向距离，并将该网格体单元编号按法向距离大小排列插入到所述壁面网格单元的法向网格体单元链表中；启动通信接口，交换区域分解重叠区网格体单元的距离信息，并将所述重叠区网格体单元按法向距离大小排列插入到与所述重叠区网格体单元对应的目标壁面网格单元的法向网格体单元链表中。
5.根据权利要求2所述的方法，其特征在于，每个所述处理器根据物理场解算器依次求解无粘欧拉方程与粘性纳维斯托克斯方程，获取流场物理解并基于所述流场物理解更新区域分解重叠区域的物理场信息的步骤，包括：每个处理器计算当地网格单元的无粘欧拉方程物理场；每个处理器计算当地网格单元的粘性纳维斯托克斯方程物理场；根据预先建立多层重叠的并行计算通信数据结构，确定接口交换区域分解的重叠区网格体单元的无粘以及粘性物理场变量；基于所述接口交换区域分解的重叠区网格体单元的无粘以及粘性物理场变量，更新区域分解重叠区域的物理场信息。6.根据权利要求2所述的方法，其特征在于，各个所述处理器根据所述区域分解重叠区域的物理场信息计算当地处理器内的边界层物理量特征量，通过并行计算函数全局规约计算全局的壁面网格单元边界层特征物理量的步骤，包括：对各个所述处理器中每个全局壁面网格单元重复执行以下步骤，直至每个所述处理器中每个全局壁面网格单元均被遍历：根据所述壁面网格单元的网格体单元链表，收集流动物理量中的无粘及粘性物理量；根据壁面网格单元的法向网格体单元的无粘和粘性物理量，计算当地边界层特征量；针对所述当地边界层特征量，确定每个所述壁面网格单元的特征点；利用并行计算函数全局规约每个所述壁面网格单元的第一特征点，以及每个所述壁面网格单元的第一特征量与法向距离，并规约每个所述壁面网格单元的第二特征点以及每个所述壁面网格单元的第二特征量与法向距离。7.根据权利要求2所述的方法，其特征在于，所述方法还包括：根据各个所述处理器计算出的全局的壁面网格单元边界层特征物理量以及当地处理器对应的物理场信息，计算出所述壁面网格单元的其它边界层物理量；根据所述其它边界层物理量计算目标物理场求解值。8.一种多面体网格原位分布式并行边界层物理量获取装置，其特征在于，包括：搜索模块，将每个处理器的壁面网格单元广播到其余处理器中，以使所述每个处理器基于直线距离建立壁面网格单元的k-d树数据结构，并进行邻近点搜索；插表模块，基于所述邻近点确定出的网格体单元的最近距离壁面网格单元，将每个处理器中的每个网格单元对应的网格体单元编号和重叠区网格体单元，按照法向距离插入法向网格体单元链表中；获取模块，每个所述处理器根据物理场解算器依次求解无粘欧拉方程与粘性纳维斯托克斯方程，获取流场物理解并基于所述流场物理解更新区域分解重叠区域的物理场信息；计算模块，各个所述处理器根据所述区域分解重叠区域的物理场信息计算当地处理器内的边界层物理量特征量，通过并行计算函数全局规约计算全局的壁面网格单元边界层特征物理量。9.一种电子设备，其特征在于，包括存储器、处理器及存储在所述存储器上并且能够在所述处理器上运行的程序，所述处理器执行所述程序时实现如权利要求1至7中任一项所述的方法。10.一种计算机可读存储介质，其特征在于，所述可读存储介质中存储有计算机程序，
所述计算机程序被执行时实现权利要求1-7中任意一项所述的方法。

技术总结
本发明公开了一种多面体网格原位分布式并行边界层物理量获取方法和装置，涉及计算流体力学等多物理场数值模拟的工程应用的技术领域，该方法基于多层相互重叠的区域分解分布式网格，可原位分布式计算当地处理器中的边界层物理量信息而不需重构整个流动物理场，解决目前方法难以高效率、低内存地获取边界层参数，导致所构造的转捩模式和转捩准则方法应用受限、计算成本较高等技术问题。计算成本较高等技术问题。计算成本较高等技术问题。


技术研发人员：林博希 张亮 初敏
受保护的技术使用者：北京凌云智擎软件有限公司
技术研发日：2023.07.12
技术公布日：2023/8/13