基于链接方向人工压缩的快速湍流数值模拟方法和装置

摘要

本申请涉及一种基于链接方向人工压缩的快速湍流数值模拟方法、装置、计算机设备和存储介质。所述方法包括：通过合成涡算法得到生成的涡的尺寸信息、数量信息和位置信息，根据尺寸信息、数量信息和位置信息得到管道输入端的波动速度，根据波动速度，通过合成湍流输入生成算法得到形成湍流所需要的外力项，将外力项加入到链接方向人工压缩方程的分布函数中，得到在管道流中加入湍流所需要增加的外力项函数，设置管道流分布函数的边界条件，在GPU中对分布函数进行迭代计算，得到网格模型中流体的密度和速度。采用本方法可以大大减少湍流数值模拟计算对全局内存的占用，提高湍流模拟在GPU上运行的效率。

说明书

技术领域

本申请涉及计算流体力学技术领域，特别是涉及一种基于链接方向人工压缩的快速湍流数值模拟方法、装置、计算机设备和存储介质。

背景技术

计算流体力学已经广泛运用在不论是工业界还是学术界的各个领域中。从航空航天领域到船舶制造业，从医学血流模拟到电影制造业，甚至芯片制造工艺等等方面，计算流体力学都发挥着至关重要的作用。计算流体力学不仅可以大幅降低物理实验的花销成本，也可以给实验的方向提供有力的证据，大大提高了实验的效率。而湍流作为流体最为重要的一种流体形式，在计算流体力学的绝大部分领域是不可被避免的。

近几年来，由于格子玻尔兹曼方法具有非常良好的局部性和天然的并行性，受到大家的青睐。而原本作为图像处理的协处理器GPU，由于其具有数量庞大的并行处理单元，已经越来越多地被运用到大规模科学与工程计算中。大量工作开始将格子玻尔兹曼方法在GPU上进行加速，并得到了非常好的效果，大大加快了模拟的速度。也有很多工作采用了多GPU的方式来增加提高模拟规模。

然而，随着求解问题的规模逐渐增大，格子玻尔兹曼方法也出现了瓶颈。格子玻尔兹曼方法一般分为两部分：碰撞步和迁移步。碰撞步在本地发生，而迁移步需要周围各自的密度分布函数。因此格子玻尔兹曼方法在每一个迭代步中每一个格子都需要保存密度分布函数，这将会占据非常大的内存。

湍流模拟占据着计算流体力学领域的重要一部分，现有的湍流模拟方法随着计算流体力学模拟规模的逐渐增大，存在计算速度慢、占用内存大的问题。

发明内容

基于此，有必要针对上述技术问题，提供一种能够提高湍流模拟计算速度、减少内存占用的基于链接方向人工压缩的快速湍流数值模拟方法、装置、计算机设备和存储介质。

一种基于链接方向人工压缩的快速湍流数值模拟方法，所述方法包括：

构建流体模拟的网格模型，获取初始化输入信息，根据初始化输入信息中包含的网格模型的坐标信息，得到链接方向人工压缩方程的管道流分布函数；

通过合成涡算法得到生成的涡的尺寸信息、数量信息和位置信息，根据所述尺寸信息、所述数量信息和所述位置信息得到管道输入端的波动速度；

根据所述波动速度，通过合成湍流输入生成算法得到形成湍流所需要的外力项；

通过将所述外力项加入到所述链接方向人工压缩方程的管道流分布函数中，得到在管道流中加入湍流所需要增加的外力项函数；

设置所述管道流分布函数的边界条件，根据所述边界条件和所述外力项函数，在GPU中对分布函数进行迭代计算，得到所述网格模型中流体的密度和速度。

在其中一个实施例中，还包括：根据所述波动速度，通过合成湍流输入生成算法得到形成湍流所需要的外力项为：

其中，

在其中一个实施例中，还包括：所述网格模型为D3Q19模型；通过将所述外力项加入到链接方向人工压缩方程的管道流分布函数中，得到在管道流中加入湍流所需要增加的外力项函数为：

其中，

在其中一个实施例中，还包括：设置所述管道流分布函数的输入端边界条件、输出端边界条件和墙边界条件，根据所述边界条件和所述外力项函数，在GPU中对分布函数进行迭代计算，得到所述网格模型中流体的密度和速度。

在其中一个实施例中，所述输入端边界条件包括输入端速度和输入端密度；所述输出端边界条件包括输出端速度和输出端密度，通过所述初始化输入信息中包含的平均速度和所述波动速度得到所述输入端速度；

通过二阶精度的外插值方法得到所述输入端密度；

设置所述输出端密度为1；

通过二阶精度的外插值方法得到所述输出端速度；

通过格子玻尔兹曼方法中的反弹格式墙边界条件得到所述墙边界条件；

根据所述边界条件和所述外力项函数，在GPU中对分布函数进行迭代计算，得到所述网格模型中流体的密度和速度。

在其中一个实施例中，在设置所述管道流分布函数的边界条件，根据所述边界条件和所述外力项函数，在GPU中对分布函数进行迭代计算，得到所述网格模型中流体的密度和速度之前，还包括：

在所述GPU的内存中，分配两份全局内存，所述全局内存用于保存计算得到的所述密度和所述速度；

在所述GPU的内存中，分配一份共享内存，所述共享内存被分为多个8*8*8的区块，所述共享内存用于保存所述密度和所述速度的备份。

一种基于链接方向人工压缩的快速湍流数值模拟装置，所述装置包括：

管道流分布函数获取模块，用于构建流体模拟的网格模型，获取初始化输入信息，根据初始化输入信息中包含的网格模型的坐标信息，得到链接方向人工压缩方程的管道流分布函数；

波动速度获取模块，用于通过合成涡算法得到生成的涡的尺寸信息、数量信息和位置信息，根据所述尺寸信息、所述数量信息和所述位置信息得到管道输入端的波动速度；

外力项获取模块，用于根据所述波动速度，通过合成湍流输入生成算法得到形成湍流所需要的外力项；

分布函数获取模块，用于通过将所述外力项加入到链接方向人工压缩方程的管道流分布函数中，得到在管道流中加入湍流所需要增加的外力项函数；

密度和速度计算模块，用于设置所述管道流分布函数的边界条件，根据所述边界条件和所述外力项函数，在GPU中对分布函数进行迭代计算，得到所述网格模型中流体的密度和速度。

一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现以下步骤：

构建流体模拟的网格模型，获取初始化输入信息，根据初始化输入信息中包含的网格模型的坐标信息，得到链接方向人工压缩方程的管道流分布函数；

通过合成涡算法得到生成的涡的尺寸信息、数量信息和位置信息，根据所述尺寸信息、所述数量信息和所述位置信息得到管道输入端的波动速度；

根据所述波动速度，通过合成湍流输入生成算法得到形成湍流所需要的外力项；

通过将所述外力项加入到所述链接方向人工压缩方程的管道流分布函数中，得到在管道流中加入湍流所需要增加的外力项函数；

设置所述管道流分布函数的边界条件，根据所述边界条件和所述外力项函数，在GPU中对分布函数进行迭代计算，得到所述网格模型中流体的密度和速度。

一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现以下步骤：

构建流体模拟的网格模型，获取初始化输入信息，根据初始化输入信息中包含的网格模型的坐标信息，得到链接方向人工压缩方程的管道流分布函数；

通过合成涡算法得到生成的涡的尺寸信息、数量信息和位置信息，根据所述尺寸信息、所述数量信息和所述位置信息得到管道输入端的波动速度；

根据所述波动速度，通过合成湍流输入生成算法得到形成湍流所需要的外力项；

通过将所述外力项加入到所述链接方向人工压缩方程的管道流分布函数中，得到在管道流中加入湍流所需要增加的外力项函数；

设置所述管道流分布函数的边界条件，根据所述边界条件和所述外力项函数，在GPU中对分布函数进行迭代计算，得到所述网格模型中流体的密度和速度。

上述基于链接方向人工压缩的快速湍流数值模拟方法、装置、计算机设备和存储介质，通过合成涡算法得到生成的涡的尺寸信息、数量信息和位置信息，根据尺寸信息、数量信息和位置信息得到管道输入端的波动速度，根据所述波动速度，通过合成湍流输入生成算法得到形成湍流所需要的外力项，将所述外力项加入到链接方向人工压缩方程的分布函数中，得到在管道流中加入湍流所需要增加的外力项函数，设置管道流分布函数的边界条件，在GPU中对分布函数进行迭代计算，得到所述网格模型中流体的密度和速度。本发明的方法大大减少了湍流数值模拟计算对全局内存的占用，提高了湍流模拟在GPU上运行的效率。

附图说明

图1为一个实施例中基于链接方向人工压缩的快速湍流数值模拟方法的流程示意图；

图2为一个实施例中合成湍流输入生成算法的作用区域示意图；

图3为一个实施例中输入端密度的边界条件示意图；

图4为一个实施例中输出端速度的边界条件示意图；

图5为一个实施例中墙边界条件示意图；

图6为一个实施例中基于链接方向人工压缩的快速湍流数值模拟装置的结构框图；

图7为一个实施例中计算机设备的内部结构图。

具体实施方式

为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。

本申请提供的基于链接方向人工压缩的快速湍流数值模拟方法，可以应用于如下应用环境中。以CPU和GPU作为湍流数值模拟硬件环境，构建湍流模拟的网络模型，将通过合成涡算法和合成湍流输入生成算法得到的外力项加入链接方向人工压缩方程的分布函数中，再根据预先设置的边界条件，在GPU中对分布函数进行迭代计算，得到网格模型中流体的密度和速度。

在一个实施例中，如图1所示，提供了一种基于链接方向人工压缩的快速湍流数值模拟方法，包括以下步骤：

步骤102，构建流体模拟的网格模型，获取初始化输入信息，根据初始化输入信息中包含的网格模型的坐标信息，得到链接方向人工压缩方程的管道流分布函数。

链接方向人工压缩方程（Link-wise artificial compressibility method）是在“Link-wise artificial compressibility method, Journal of ComputationalPhysics 231(2012) 5109 – 5143”中提出的。

本实施例中使用的网格模型是D3Q19模型，初始化输入信息中包含D3Q19模型的坐标信息，根据初始化输入信息中包含的网格模型的坐标信息，得到链接方向人工压缩方程的管道流分布函数为：

其中，

步骤104，通过合成涡算法得到生成的涡的尺寸信息、数量信息和位置信息，根据尺寸信息、数量信息和位置信息得到管道输入端的波动速度。

合成涡算法（Synthetic eddy method）可看做是在管道流输入端形成涡的输入端边界条件。合成涡算法最初由Jarrin在“A synthetic-eddy-method for generatinginflow conditions for large-eddy simulations,International Journal of Heatand Fluid Flow 27 (2006) p.585–593”中提出，本实施例所采用的是Skillen在“Accuracy and efficiency improvements in synthetic eddy methods,International Journal of Heat and FluidFlow 62 (2016) 386 – 394”所提出的合成涡算法的改进版本，在合成涡算法中，输入端的速度分成平均速度和波动速度两个部分，平均速度是由初始化输入信息预先设置的，波动速度是根据生成的涡的尺寸信息、数量信息和位置信息计算得到的。

步骤106，根据波动速度，通过合成湍流输入生成算法得到形成湍流所需要的外力项。

合成湍流输入生成算法（Synthetic turbulence inflow generator -STIG）源于Klein在“A digital filter based generation of inflow data for spatiallydeveloping direct numerical or large eddysimulations, Journal ofComputational Physics 186 (2003) 652 – 665”中提出的数字滤波器概念。由于管道流速度可以分解成平均速度和波动速度，虽然可以将合成涡算法直接作用于管道流的输入端，但实验表明，将合成涡算法直接作用于输入端的链接方向人工压缩方程的管道流在进行数值模拟时并不稳定，程序易崩溃，因此，在本实施例中，将合成涡算法得到的波动速度作为源项加入到合成湍流输入生成算法，将合成湍流输入生成算法作为外力项作用于管道流的输入端来形成湍流。

步骤108，通过将外力项加入到链接方向人工压缩方程的管道流分布函数中，得到在管道流中加入湍流所需要增加的外力项函数。

将合成湍流输入生成算法作为外力项增加到链接方向人工压缩方程的分布函数中，如图2所示，在管道流靠近输入端的位置，形成了“合成湍流输入生成方法”的作用区域，作用区域长度为

其中，

步骤110，设置管道流分布函数的边界条件，根据边界条件和外力项函数，在GPU中对分布函数进行迭代计算，得到网格模型中流体的密度和速度。

在本发明之前，还没有将链接方向人工压缩方法用于管道流模拟的研究，因此，本发明也首次研究将链接方向人工压缩方法用于管道流模拟时的边界条件。由于链接方向人工压缩方法既有格子玻尔兹曼方法的特点，又有有限差分方法的特点，为了达到最好的效果，本实施例中同时结合了格子玻尔兹曼方法和有限差分方法来设置管道流的边界条件。

原本作为图像处理的协处理器GPU，由于其具有数量庞大的并行处理单元，已经越来越多地被运用到大规模科学与工程计算中。本实施例将GPU作为模拟计算的硬件环境，可以提升数值模拟的速度。

上述基于链接方向人工压缩的快速湍流数值模拟方法中，通过合成涡算法得到生成的涡的尺寸信息、数量信息和位置信息，根据尺寸信息、数量信息和位置信息得到管道输入端的波动速度，根据波动速度，通过合成湍流输入生成算法得到形成湍流所需要的外力项，将外力项加入到链接方向人工压缩方程的分布函数中，得到在管道流中加入湍流所需要增加的外力项函数，设置管道流分布函数的边界条件，在GPU中对分布函数进行迭代计算，得到网格模型中流体的密度和速度。本发明的方法大大减少了湍流数值模拟计算对全局内存的占用，提高了湍流模拟在GPU上运行的效率。

在一个实施例中，根据波动速度，通过合成湍流输入生成算法得到形成湍流所需要的外力项为：

其中，

在一个实施例中，网格模型为D3Q19模型；通过将外力项加入到链接方向人工压缩方程的管道流分布函数中，得到在管道流中加入湍流所需要增加的外力项函数为：

其中，

在一个实施例中，还包括：设置管道流分布函数的输入端边界条件、输出端边界条件和墙边界条件，根据边界条件和外力项函数，在GPU中对分布函数进行迭代计算，得到网格模型中流体的密度和速度。输入端边界条件包括输入端速度和输入端密度；输出端边界条件包括输出端速度和输出端密度；通过初始化输入信息中包含的平均速度和波动速度得到输入端速度；通过二阶精度的外插值方法得到输入端密度；设置输出端密度为1；通过二阶精度的外插值方法得到输出端速度；通过格子玻尔兹曼方法中的反弹格式墙边界条件得到墙边界条件；根据边界条件和外力项函数，在GPU中对分布函数进行迭代计算，得到网格模型中流体的密度和速度。

具体的，如图3所示，位置0在输入端边界，位置1和2是相邻两个格子，则输入端的密度可以通过下式得到：

如图4所示，位置0在输出端边界，而位置1和2是相邻两个格子，输出端的速度可以通过下式得到：

如图5所示，墙边界条件采用格子玻尔兹曼方法中常见的反弹格式墙边界条件，将墙的位置假设在距离最近格子点半个格子的位置，反弹的密度分布函数通过相反方向的密度分布函数计算得到：

其中

在一个实施例中，在设置管道流分布函数的边界条件，根据边界条件和外力项函数，在GPU中对分布函数进行迭代计算，得到网格模型中流体的密度和速度之前，还包括：通过提高管道流输入端和输出端流体的粘度参数，在管道流输入端和输出端设置海绵层。

由于链接方向人工压缩方法在输入端和输出端的模拟非常不稳定，容易中途程序崩溃，因此在输入和输出的部分采用了海绵层来增加流体粘度。除了两端流体粘度不同意外，其他部分的流体粘度始终相同。

在一个实施例中，在设置管道流分布函数的边界条件，根据边界条件和外力项函数，在GPU中对分布函数进行迭代计算，得到网格模型中流体的密度和速度之前，还包括：在GPU的内存中，分配两份全局内存，全局内存用于保存计算得到的密度和速度；在GPU的内存中，分配一份共享内存，共享内存被分为多个8*8*8的区块，共享内存用于保存密度和速度的备份。

分配两份全局内存同时读写，可以提高GPU数值模拟的计算效率。由于共享内存在硬件上距离GPU的计算单元近，分配一份共享内存来保存密度和速度的备份，可以通过空间局部性来提高对全局内存的访问效率。

应该理解的是，虽然图1的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，图1中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。

在一个实施例中，如图6所示，提供了一种基于链接方向人工压缩的快速湍流数值模拟装置，包括：管道流分布函数获取模块602、波动速度获取模块604、外力项获取模块606、分布函数获取模块608和密度和速度计算模块610，其中：

管道流分布函数获取模块602，用于构建流体模拟的网格模型，获取初始化输入信息，根据初始化输入信息中包含的网格模型的坐标信息，得到链接方向人工压缩方程的管道流分布函数。

波动速度获取模块604，用于通过合成涡算法得到生成的涡的尺寸信息、数量信息和位置信息，根据尺寸信息、数量信息和位置信息得到管道输入端的波动速度。

外力项获取模块606，用于根据波动速度，通过合成湍流输入生成算法得到形成湍流所需要的外力项。

分布函数获取模块608，用于通过将外力项加入到链接方向人工压缩方程的管道流分布函数中，得到在管道流中加入湍流所需要增加的外力项函数。

密度和速度计算模块610，用于设置管道流分布函数的边界条件，根据边界条件和外力项函数，在GPU中对分布函数进行迭代计算，得到网格模型中流体的密度和速度。

密度和速度计算模块610还用于获取管道流分布函数的输入端边界条件、输出端边界条件和墙边界条件，根据边界条件和外力项函数，在GPU中对分布函数进行迭代计算，得到网格模型中流体的密度和速度。

密度和速度计算模块610还用于通过初始化输入信息中包含的平均速度和波动速度得到输入端速度；通过二阶精度的外插值方法得到输入端密度；设置输出端密度为1；通过二阶精度的外插值方法得到输出端速度；通过格子玻尔兹曼方法中的反弹格式墙边界条件得到墙边界条件；根据边界条件和外力项函数，在GPU中对分布函数进行迭代计算，得到网格模型中流体的密度和速度。

密度和速度计算模块610还用于通过提高管道流输入端和输出端流体的粘度参数，在管道流输入端和输出端设置海绵层。

密度和速度计算模块610还用于在GPU的内存中，分配两份全局内存，全局内存用于保存计算得到的密度和速度；在GPU的内存中，分配一份共享内存，共享内存被分为多个8*8*8的区块，共享内存用于保存密度和速度的备份。

关于基于链接方向人工压缩的快速湍流数值模拟装置的具体限定可以参见上文中对于基于链接方向人工压缩的快速湍流数值模拟方法的限定，在此不再赘述。上述基于链接方向人工压缩的快速湍流数值模拟装置中的各个模块可全部或部分通过软件、硬件及其组合来实现。上述各模块可以硬件形式内嵌于或独立于计算机设备中的处理器中，也可以以软件形式存储于计算机设备中的存储器中，以便于处理器调用执行以上各个模块对应的操作。

在一个实施例中，提供了一种计算机设备，该计算机设备可以是终端，其内部结构图可以如图7所示。该计算机设备包括通过系统总线连接的处理器、存储器、网络接口、显示屏和输入装置。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现一种基于链接方向人工压缩的快速湍流数值模拟方法。该计算机设备的显示屏可以是液晶显示屏或者电子墨水显示屏，该计算机设备的输入装置可以是显示屏上覆盖的触摸层，也可以是计算机设备外壳上设置的按键、轨迹球或触控板，还可以是外接的键盘、触控板或鼠标等。

本领域技术人员可以理解，图7中示出的结构，仅仅是与本申请方案相关的部分结构的框图，并不构成对本申请方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。

在一个实施例中，提供了一种计算机设备，包括存储器和处理器，该存储器存储有计算机程序，该处理器执行计算机程序时实现上述方法实施例中的步骤。

在一个实施例中，提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述方法实施例中的步骤。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器（ROM）、可编程ROM（PROM）、电可编程ROM（EPROM）、电可擦除可编程ROM（EEPROM）或闪存。易失性存储器可包括随机存取存储器（RAM）或者外部高速缓冲存储器。作为说明而非局限，RAM以多种形式可得，诸如静态RAM（SRAM）、动态RAM（DRAM）、同步DRAM（SDRAM）、双数据率SDRAM（DDRSDRAM）、增强型SDRAM（ESDRAM）、同步链路（Synchlink） DRAM（SLDRAM）、存储器总线（Rambus）直接RAM（RDRAM）、直接存储器总线动态RAM（DRDRAM）、以及存储器总线动态RAM（RDRAM）等。

以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

以上实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求为准。