CN113361217B - 一种高效的两相流无网格数值模型实施方法、装置 - Google Patents

Abstract

本发明提供了一种高效的两相流无网格数值模型实施方法、装置，属于数值模拟技术领域。所述两相流无网格数值模型实施方法的主要过程为：当需要空气粒子时，在计算域中生成空气粒子，减小时间步长；当不需要空气粒子时，将计算域内的空气粒子删除，同时增大时间步长；通过上述方法，可以将两相流SPH模型的计算效率提高8倍以上。本发明能够显著提升两相流SPH模型的计算效率，节约计算成本。

Description

一种高效的两相流无网格数值模型实施方法、装置

技术领域

本发明属于数值模拟技术领域，具体涉及一种高效的两相流无网格数值模型实施方法、装置。

背景技术

随着社会经济的发展，人类为了开发利用海洋资源，在海洋中修建了大量的海洋结构物。波浪作用下海洋结构物的水动力特性一直是工程师们所关注的重要问题。近年来，随着计算流体力学数值模拟技术的不断成熟，数值模拟逐渐成为分析波浪与海洋结构物相互作用的有效手段。

光滑粒子流体动力学方法(Smoothed Particle Hydrodynamics，SPH)是一种新兴的拉格朗日形式的无网格粒子法，其最初主要用于解决三维开放空间的天体物理学问题，现已被应用到海洋工程问题。SPH方法较之传统的基于网格的数值算法，由于其粒子性和拉格朗日性，在求解大变形问题上有巨大的优势。对于一般模拟波浪与结构物的相互作用而言，具有单相流(水)SPH模型和两相流(空气+水)SPH模型。单相流SPH模型不考虑空气的作用，其计算效率高、稳健性好，能较好的预测波浪作用下结构物的水动力特性。但是对于带有腔室的海洋结构物，例如开孔沉箱的消浪室结构、高桩码头底部的空腔结构等，波浪在与这些结构物相互作用时，需要考虑腔室内空气变化对水体运动的影响，此时需要采用两相流SPH模型。两相流SPH模型虽然能够考虑腔室内空气变化对水体运动的影响，但是由于加入空气粒子，导致计算节点显著增多，同时由于空气密度很小，导致计算时间步长极大缩短，这些因素都明显降低了两相流SPH数值模型的计算效率。

发明内容

针对现有技术中存在不足，本发明提供了一种高效的两相流无网格数值模型实施方法、装置，显著提高计算效率，节约计算成本。

本发明是通过以下技术手段实现上述技术目的的。

一种高效的两相流无网格数值模型实施方法，包括步骤：

S1，建立数值波浪水槽，在流体区域布置水粒子；

S2，当结构物空腔处于封闭状态时，封闭次数k＝k+1，执行S3；若结构物空腔不满足封闭状态，则封闭次数k＝0，执行S6；

S3，若k>1，执行S4；若k＝1，执行S5；

S4，水粒子和空气粒子均进行控制方程计算，并使用四阶龙格库塔法对控制方程计算结果进行积分，其中用于获取时间步长的库朗数CFL值取0.25；

S5，生成空气粒子，水粒子和空气粒子均进行控制方程计算，并使用四阶龙格库塔法对控制方程计算结果进行积分，其中用于获取时间步长的库朗数CFL值取0.25；

S6，仅水粒子进行控制方程计算，并使用四阶龙格库塔法对控制方程计算结果进行积分，其中用于获取时间步长的库朗数CFL值取2.2；

S7，结构物腔室满足封闭条件时，封闭次数k＝k+1，执行S4；结构物腔室不满足封闭条件时，k＝0，删除所有空气粒子，水粒子进行控制方程计算，并使用四阶龙格库塔法对控制方程计算结果进行积分，其中用于获取时间步长的库朗数CFL值调整为2.2；

S8，判断是否满足时间终止条件：t>tend，若满足，则结束计算；不满足时间终止条件时，t＝t+△t，返回S2。

进一步的技术方案，所述控制方程为：

其中：ρ、u、p、g分别为粒子的密度、速度、压强和重力加速度，▽表示导数，W表示光滑核函数，r_i表示粒子i的位置，V表示体积，c₀表示声速，ρ₀表示水的初始密度，m为质量，j表示粒子j，F(p_i，p_j)是压力梯度项▽p的离散形式。

进一步的技术方案，所述采用四阶龙格库塔对控制方程计算结果进行积分计算，时间步长Δt按下式计算：

其中：h为光滑长度。

更进一步的技术方案，将流体分为三部分：Ω₁表示所有的空气粒子，Ω₂表示与空气粒子相互作用的水粒子，Ω₃表示与空气粒子不相互作用的水粒子，当粒子i和粒子j相互作用时，压力梯度项具体为：

(1)若粒子i或粒子j属于Ω₁，或者，粒子i和粒子j均属于Ω₂时：

(2)若粒子i属于Ω₂且粒子j属于Ω₃时：

(3)若i粒子属于Ω₃且j粒子属于Ω₂时：

(4)若i粒子和j粒子均属于Ω₃时：

其中，p_b为背景压强。

进一步的技术方案，所述生成空气粒子的具体过程为：

在空腔内按照初始粒子分布间距Δx布满空气粒子，定义空气粒子的C值为1，计算每个空气粒子的λ值：

此时j粒子仅为水粒子，随后删除λ_i≥1.5的空气粒子；

计算保留下来的每个空气粒子的CT值：

此时j粒子为水粒子和空气粒子，随后删除CT<0的空气粒子；

对最后保留下来的空气粒子的位置进行如下修正：

其中：β为常数，粒子i仅为空气粒子，粒子j为空气粒子和水粒子；

上述空气粒子修正过程使用四阶龙格库塔法进行数值积分，时间步长采用

当修正过程中的空气粒子的最大速度小于0.01时，修正过程完成。

一种高效的两相流无网格数值模型实施装置，其特征在于，包括：

模型建立模块，用于建立数值波浪水槽；

压强梯度离散模块，用于水粒子或空气粒子进行控制方程的计算；

积分模块，用于使用四阶龙格库塔法对控制方程计算结果进行积分；

空气粒子生成模块，用于生成空气粒子。

一种电子设备，包括存储器和处理器；

所述存储器用于存储计算机程序；

所述处理器用于执行所述计算机程序并在执行所述计算机程序时实现上述两相流无网格数值模型实施方法。

一种存储介质，所述存储介质存储有计算机程序，所述计算机程序被处理器执行时使所述处理器执行上述两相流无网格数值模型实施方法。

本发明的有益效果为：

本发明通过耦合单相流SPH模型与两相流SPH模型，提出一种高效的两相流无网格数值模型实施方法，即：当需要空气粒子时，在计算域中生成空气粒子，减小时间步长；当不需要空气粒子时，将计算域内的空气粒子删除，同时增大时间步长；通过这种方式，可以将两相流SPH模型的计算效率提高8倍以上。本发明的实施方法不仅适用于SPH方法，也适用于其他无网格方法，例如：移动粒子半隐式法(Moving Particle Semi-Implicit，MPS)等。

附图说明

图1为本发明所述两相流无网格数值模型实施方法流程图；

图2为本发明局部粒子背景压强计算结果图；

图3为本发明全体粒子背景压强计算结果图；

图4为本发明无背景压强计算结果图；

图5为本发明插值点转化为空气粒子图；

图6为本发明删除λ_i≥1.5的空气粒子后的粒子分布图；

图7为本发明删除CT<0的空气粒子后的粒子分布图；

图8为本发明空气粒子位置修正后的粒子分布图。

具体实施方式

下面结合附图以及具体实施例对本发明作进一步的说明，但本发明的保护范围并不限于此。

SPH方法主要通过光滑核函数建立了粒子与粒子之间的联系，并引进积分表达式。可见光滑函数对SPH方法的实现起到了至关重要的作用。本发明采用的光滑核函数为高斯光滑核函数：

其中：W表示光滑核函数；两个粒子间的距离s＝||r-r'||，r表示粒子i的位置，r′表示粒子j的位置；γ为截断半径系数，通常γ＝3；h为光滑长度，对于高斯光滑核函数

Δx为粒子的初始间距。

高斯光滑核函数一个重要的性质就是它的导数可以直接计算：

利用纳维尔-斯托克斯方程(N-S方程)描述流体运动：

其中：ρ、u、p、g分别为粒子的密度、速度、压强和重力加速度，▽·表示散度，▽表示导数。

SPH方法包括单相流SPH模型和两相流SPH模型，下面就这两种模型分别进行介绍。

单相流SPH模型

使用SPH方法对式(3)进行离散，同时使用状态方程将式(3)闭合，最终得到的控制方程为：

其中：V表示体积；c₀表示声速；ρ₀表示水的初始密度，本实施例取1000kg/m³；m为质量。

采用四阶龙格库塔对式(4)进行积分计算，时间步长Δt按下式计算：

其中：库朗数CFL＝2.2，c＝c₀。

两相流SPH模型

两相流SPH模型与单相流SPH模型在控制方程上的唯一区别在于：在两相流SPH模型中，为了避免空气中极易产生的拉伸不稳定性，引入了背景压强。两相流SPH模型的控制方程为：

其中：p_b为背景压强，一般取值为1000-2000；ρ_0i表示粒子i的初始密度，若粒子为水粒子，则取值1000kg/m³，若粒子为气体粒子，则取值1.29kg/m³。

采用四阶龙格库塔法对式(6)进行积分计算，时间步长Δt按下式计算：

其中：CFL＝0.25。

由上述可知：(1)两相流SPH模型的时间步长几乎为单相流SPH模型时间步长的1/10；(2)两相流SPH模型由于引入了空气粒子，增加了计算节点。以上两点使得两相流SPH模型的计算效率远低于单相流SPH模型。为了克服这一问题，本发明将单相流SPH模型与两相流SPH模型进行耦合，提出了一种高效的两相流无网格数值模型实施方法，如图1所示，具体包括如下步骤：

步骤(1)，建立数值波浪水槽，在流体区域布置水粒子，水粒子初始间距为Δx，定义水粒子和固壁边界粒子的C值为-1。

步骤(2)，根据结构物的结构特征，判断结构物空腔是否封闭，当结构物空腔处于封闭状态时，封闭次数k＝k+1，执行步骤(3)；结构物空腔不满足封闭状态时，封闭次数k＝0，执行步骤(6)。

步骤(3)，判断k是否大于1，若k>1，执行步骤(4)；若k＝1，执行步骤(5)。

步骤(4)，水粒子和空气粒子均进行控制方程(11)计算，再使用四阶龙格库塔法对控制方程计算结果进行积分，其中用于获取时间步长的库朗数CFL值取0.25。

步骤(5)，生成空气粒子，水粒子和空气粒子均进行控制方程(11)计算，再使用四阶龙格库塔法对控制方程计算结果进行积分，其中用于获取时间步长的库朗数CFL值取0.25。

步骤(6)，仅水粒子进行控制方程(11)计算，再使用四阶龙格库塔法对控制方程计算结果进行积分，其中用于获取时间步长的库朗数CFL值取2.2。

步骤(7)，再次判断结构物空腔是否封闭，结构物腔室满足封闭条件时，封闭次数k＝k+1，执行步骤(4)；结构物腔室不满足封闭条件时，k＝0，此时删除所有空气粒子，对水粒子进行控制方程(11)计算，再使用四阶龙格库塔法对控制方程计算结果进行积分，其中用于获取时间步长的库朗数CFL值取2.2。

步骤(8)，判断是否满足时间终止条件t>tend，若满足，则结束计算；不满足时间终止条件时，t＝t+△t，返回步骤(2)。

为了实现单相流SPH模型与两相流SPH模型进行耦合，本发明主要通过以下两项技术：(1)局部背景压强的施加，(2)空气粒子的生成。

(1)局部背景压强的施加

在上述单相流SPH模型和两相流SPH模型中，压力梯度项▽p_i的离散形式F(p_i，p_j)可以写为：

对式(8)变形得到：

上式右端第一项为压力梯度的精确粒子插值，以下主要分析公式(9)右端的第二项。

单相流SPH模型中，当粒子位于流体内部时，即光滑核函数不被截断时，公式(9)右端第二项理论上为0，但是，由于粒子分布可能没有达到理想均匀程度，所以该项实际并不为0，此时若p_i>0，该项会使粒子向周围粒子缺失的地方移动，即会使粒子分布更加均匀；如果p_i<0时，会使粒子分布更加不均匀。所幸，对于大部分以重力为主要作用的波浪问题，p_i>0。当流体粒子位于自由液面区域时，即便光滑核函数被截断，由于此时p_i较小，公式(9)中的第二项也近似为0，因此该项只会造成很小的自由液面漂移，对计算结果影响不大。

两相流SPH模型中，由于引入了背景压强，则公式(8)变为：

p_b的一般取值为1000-2000，由于两相流SPH模型中没有自由液面，所有粒子的光滑核函数均不会被截断，因此背景压强的引入能够使得粒子分布均匀，避免拉伸不稳定性。

在单相流SPH与两相流SPH的耦合模型中，存在自由液面，即光滑核函数被截断，如果所有流体粒子均引入背景压强，会造成强烈的自由液面漂移，严重影响计算。针对这一问题，给出一套局部背景压强引入方法：

其控制方程为：

将流体分为三部分，Ω₁表示所有的空气粒子，Ω₂表示与空气粒子相互作用的水粒子，Ω₃表示与空气粒子不相互作用的水粒子。当粒子i和粒子j相互作用时：

若粒子i或粒子j属于Ω₁，或者，粒子i和粒子j均属于Ω₂时：

若粒子i属于Ω₂且粒子j属于Ω₃时：

若i粒子属于Ω₃且j粒子属于Ω₂时：

若i粒子和j粒子均属于Ω₃时：

图2、图3和图4分别给出了局部粒子背景压强的计算结果、全体粒子背景压强的计算结果和无背景压强的计算结果。可以看出：局部粒子背景压强既能达到光滑的压强场，也能保证空气粒子的均匀分布。

(2)空气粒子的生成

在SPH方法中，粒子分布的均匀程度对数值计算结果的精确度至关重要，因此提出一套生成均匀分布的空气粒子的流程显得十分重要。

首先在结构物的空腔内按照初始粒子分布间距Δx布置一系列插值点，当需要生成空气粒子时，这些插值点转化为空气粒子，如图5。可以看出此时新生成的粒子与其他粒子相互渗透，因此下一步将相互渗透的粒子删除。

由于自由液面处流体粒子的位置散度小于1.5。基于此原理，计算每个空气粒子i在水中插值的位置散度λ的值，删除λ_i≥1.5的空气粒子，参见图6。

此时j粒子仅为水(fluid)粒子。

删除上述部分粒子后，在水体自由液面处，存在着一系列与水粒子十分接近的空气粒子，为了删除这部分粒子，定义空气粒子的C值为1，水粒子和数值波浪水槽固壁边界粒子的C值为-1。计算每个空气粒子i的CT值：

此时j粒子为水(fluid)粒子和空气(gas)粒子，且i与j粒子不为同一粒子。随后删除CT<0的空气粒子，参见图7。

可以看到，最后保留下来的空气粒子与水粒子间存在空隙，因此需要对这部分空气粒子的位置进行修正，修正公式如下：

其中：常数β＝0.002，粒子i仅为空气粒子，粒子j为空气粒子和水粒子，i与j不为同一粒子。上述空气粒子修正过程使用四阶龙格库塔法进行数值积分，时间步长采用

当修正过程中的空气粒子的最大速度小于0.01时，修正过程完成，参见图8。

空气粒子的压强和速度通过下式从水粒子中插值得到，为使用所有空气粒子获得压强和速度值，插值过程进行20次循环：

当计算中需要空气粒子时，在需要区域生成空气粒子，此时CFL＝0.25，当不需要空气粒子时，删除空气粒子，此时CFL＝2.2，通过这种方式，可极大提高两相流SPH模型的计算效率，根据发明人的数值算例分析，在相同计算条件下，采用本发明方法的估计计算效率为传统两相流SPH模型计算效率的8倍。

本发明高效的两相流无网格数值模型实施装置，包括模型建立模块、压强梯度离散模块、积分模块和空气粒子生成模块；

模型建立模块，用于建立数值波浪水槽；压强梯度离散模块，用于水粒子或空气粒子进行控制方程的计算；积分模块，用于使用四阶龙格库塔法对控制方程计算结果进行积分；空气粒子生成模块，用于生成空气粒子。

基于与高效的两相流无网格数值模型实施方法相同的发明构思，本申请还提供了一种电子设备，该电子设备包括一个或多个处理器和一个或多个存储器，存储器中存储了计算机可读代码，其中，计算机可读代码当由一个或多个处理器执行时，进行高效的两相流无网格数值模型实施。其中，存储器可以包括非易失性存储介质和内存储器；非易失性存储介质可存储操作系统和计算机可读代码。该计算机可读代码包括程序指令，该程序指令被执行时，可使得处理器执行任意一种高效的两相流无网格数值模型实施方法。处理器用于提供计算和控制能力，支撑整个电子设备的运行。存储器为非易失性存储介质中的计算机可读代码的运行提供环境，该计算机可读代码被处理器执行时，可使得处理器执行任意一种高效的两相流无网格数值模型实施方法。

应当理解的是，处理器可以是中央处理单元(Central Processing Unit，CPU)，该处理器还可以是其他通用处理器、数字信号处理器(Digital Signal Processor，DSP)、专用集成电路(Application Specific Integrated Circuit，ASIC)、现场可编程门阵列(Field-Programmable Gate Array，FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。其中，通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。

本申请的实施例中还提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机可读代码，所述计算机可读代码中包括程序指令，所述处理器执行所述程序指令，实现本申请的高效的两相流无网格数值模型实施方法。

其中，所述计算机可读存储介质可以是前述实施例所述电子设备的内部存储单元，例如所述计算机设备的硬盘或内存。所述计算机可读存储介质也可以是所述电子设备的外部存储设备，例如所述电子设备上配备的插接式硬盘、智能存储卡(SmartMedia Card，SMC)、安全数字(Secure Digital，SD)卡、闪存卡(Flash Card)等。

所述实施例为本发明的优选的实施方式，但本发明并不限于上述实施方式，在不背离本发明的实质内容的情况下，本领域技术人员能够做出的任何显而易见的改进、替换或变型均属于本发明的保护范围。

Claims (8)

1.一种高效的两相流无网格数值模型实施方法，其特征在于，包括步骤：

S1，建立数值波浪水槽，在流体区域布置水粒子；

S2，当结构物腔室处于封闭状态时，封闭次数k＝k+1，执行S3；若结构物室腔不满足封闭状态，则封闭次数k＝0，仅水粒子进行控制方程计算，并使用四阶龙格库塔法对控制方程计算结果进行积分，其中用于获取时间步长的库朗数CFL值取2.2；

S3，若k>1，水粒子和空气粒子均进行控制方程计算，并使用四阶龙格库塔法对控制方程计算结果进行积分，其中用于获取时间步长的库朗数CFL值取0.25；若k＝1，生成空气粒子，水粒子和空气粒子均进行控制方程计算，并使用四阶龙格库塔法对控制方程计算结果进行积分，其中用于获取时间步长的库朗数CFL值取0.25；

S4，t＝t+△t，结构物腔室满足封闭条件时，封闭次数k＝k+1，水粒子和空气粒子均进行控制方程计算，并使用四阶龙格库塔法对控制方程计算结果进行积分，其中用于获取时间步长的库朗数CFL值取0.25；结构物腔室不满足封闭条件时，k＝0，删除所有空气粒子，水粒子进行控制方程计算，并使用四阶龙格库塔法对控制方程计算结果进行积分，其中用于获取时间步长的库朗数CFL值取2.2，其中△t为时间步长；

S5，判断是否满足时间终止条件：t>tend，若满足，则结束计算；不满足时间终止条件时，返回S2。

2.根据权利要求1所述的两相流无网格数值模型实施方法，其特征在于，所述控制方程为：

其中：ρ、u、p、g分别为粒子的密度、速度、压强和重力加速度，W表示光滑核函数，r_i表示粒子i的位置，V表示体积，c₀表示声速，ρ₀表示水的初始密度，m为质量，j表示粒子j，F(p_i，p_j)是压力梯度项

的离散形式。

3.根据权利要求1所述的两相流无网格数值模型实施方法，其特征在于，采用四阶龙格库塔对控制方程计算结果进行积分计算，时间步长Δt按下式计算：

其中：h为光滑长度，c＝c₀，表示声速。

4.根据权利要求2所述的两相流无网格数值模型实施方法，其特征在于，将流体分为三部分：Ω₁表示所有的空气粒子，Ω₂表示与空气粒子相互作用的水粒子，Ω₃表示与空气粒子不相互作用的水粒子，当粒子i和粒子j相互作用时，压力梯度项具体为：

(1)若粒子i或粒子j属于Ω₁，或者，粒子i和粒子j均属于Ω₂时：

(2)若粒子i属于Ω₂且粒子j属于Ω₃时：

(3)若i粒子属于Ω₃且j粒子属于Ω₂时：

(4)若i粒子和j粒子均属于Ω₃时：

其中，p_b为背景压强。

5.根据权利要求2所述的两相流无网格数值模型实施方法，其特征在于，所述生成空气粒子的具体过程为：

在腔室内按照初始粒子分布间距Δx布满空气粒子，定义空气粒子的C值为1，计算每个空气粒子的λ值：

此时j粒子仅为水粒子，定义水粒子的C值为-1，随后删除λ_i≥1.5的空气粒子；

计算保留下来的每个空气粒子的CT值：

此时j粒子为水粒子和空气粒子，随后删除CT<0的空气粒子；

对最后保留下来的空气粒子的位置进行如下修正：

其中：β为常数，粒子i仅为空气粒子，粒子j为空气粒子和水粒子，W表示光滑核函数，r_i表示粒子i的位置，V表示体积，u为粒子的速度；

上述空气粒子修正过程使用四阶龙格库塔法进行数值积分，时间步长采用

当修正过程中的空气粒子的最大速度小于0.01时，修正过程完成。

6.一种实现权利要求1-5任一项权利要求所述的两相流无网格数值模型实施方法的装置，其特征在于，包括：

模型建立模块，用于建立数值波浪水槽；

压强梯度离散模块，用于水粒子或空气粒子进行控制方程的计算；

积分模块，用于使用四阶龙格库塔法对控制方程计算结果进行积分；

空气粒子生成模块，用于生成空气粒子。

7.一种电子设备，其特征在于，包括存储器和处理器；

所述存储器用于存储计算机程序；

所述处理器用于执行所述计算机程序并在执行所述计算机程序时实现如权利要求1-5任一项所述的两相流无网格数值模型实施方法。

8.一种存储介质，其特征在于，所述存储介质存储有计算机程序，所述计算机程序被处理器执行时使所述处理器执行如权利要求1-5任一项所述的两相流无网格数值模型实施方法。

Images