一种无网格数值模型的模拟方法

摘要

本发明公开了一种无网格数值模型的模拟方法、计算机设备和计算机可读存储介质。本发明所提供的无网格数值模型的模拟方法构建了改进黎曼解无网格算法的浮体、流体相互作用的数值模型，相较于相关技术中的无网格方法，充分考虑了水体发生翻卷、破碎、飞溅等强非线性因素对流场的影响，改善了应力不稳定性、粒子压力场震荡、自由面表达精度不高的问题，提高了流场压力计算精度，在造波问题中达到了较高的精度，强冲击性的固‑液相互作用时的有效性及精度，使得本发明发能适用于模拟非线性波浪与浮体相互作用的场景，具有更加广泛的适用环境，从而在船舶与海洋工程、近海工程等领域均能应用，具有较高的使用价值。

说明书

技术领域

本发明涉及船舶力学技术领域，具体涉及一种无网格数值模型的模拟方法、一种计算机设备和一种计算机可读存储介质。

背景技术

流体与浮体的相互作用是一个非常复杂的船舶水动力学问题，长期以来，船舶与海洋结构物在流体中的载荷与响应一直是舰船总体工程领域所关注的核心研究内容，经过各国科技工作者半个多世纪的不断探索，流体与结构物相互作用的数值模拟和模型试验研究都有了长足的发展，从二维到三维，从频域分析到时域模拟，从线性假设到非线性理论，结构物所处的波浪环境也由最初的规则线性小振幅波逐渐发展到不规则波浪。现如今，伴随着世界各国以探索海洋为标志的“蓝色革命”的兴起，民用海洋工程及军用作战舰艇的作业环境均开始向深海恶劣海域扩张，强非线性水体与结构物相互作用机理成为海洋武器装备所必须考量的设计要素。关于该类问题的研究目前尚处于起步阶段，还未有学术界公认的行之有效的数值方法，各类方法仍处于发展阶段，距实际工程应用仍有一定距离。

目前在计算流体力学领域占据主导地位的方法主要有势流理论和粘流理论。势流理论中，通常假定浮体所处的海洋环境是均匀、不可压缩、无粘和无旋的理想流场，因而流场的速度势满足拉普拉斯方程，并满足流场自由表面和底部条件、浮体物面条件及初始条件与辐射条件，求解浮体水动力的困难在于如何满足非线性定解条件，为此国内外诸多专家开展了大量研究，建立了二阶或更高阶的数值模型，对非线性波浪与结构物相互作用所产生的绕射和辐射问题进行了研究；粘流理论中，流固耦合运动则直接应用Navier-Stokes方程和特定的边界条件进行问题求解，通过直接求解N-S方程模拟非线性波浪问题，不仅可以充分考查由于粘性作用而导致的各种复杂流动特征，使数值模型更逼近真实流场，而且可以处理船舶与海洋结构物在波浪激励下的大幅运动，其难点在于求解过程中需要不断修正或迭代压力场，需要高精度的离散格式，在有网格模型中还要处理复杂的网格运动及变形，计算难度较大，对自由面的表达精度不高。综上，采用传统的数值方法在处理强非线性流固耦合问题时所存在的问题可总结为以下几个方面：

(1)不能充分考虑水体发生翻卷、破碎、飞溅等强非线性因素对流场的影响；

(2)浮体大幅运动使流固交互界面的相互作用描述成为难点；

(3)自由液面的大变形使水体与浮式结构物相互作用的时域求解存在困难。

发明内容

为改善上述技术问题至少之一，本发明第一方面的技术方案提供了一种无网格数值模型的模拟方法。

本发明第一方面技术方案所提供的无网格数值模型的模拟方法，用于对非线性波浪与浮体相互作用进行模拟，包括如下步骤：

步骤S1，初始化流场，对每个流体粒子赋予各物理量初值，物理量初值包括位置、密度、速度、压力；

步骤S2，在流场内布置固壁边界粒子；

步骤S3，进行NNPS粒子搜索，采用高阶MUSCL格式对粒子接触间断处的黎曼问题初值进行重构(ρ，u，w)

步骤S4，采用时间步进方法，解出守恒型变量(ω，ωρ，ωρv)，相除后得到第 N+1步粒子的速度、密度、体积，将密度代入状态方程得到压力，结合XSPH法修正速度，更新粒子位置，若达到输出条件即输出数据；

步骤S5，之后将时间步向后推进，并返回步骤S3，直至完成设定时间内全部时间步的计算。

在上述技术方案的基础上，进一步地，步骤S2具体为：设置固壁边界粒子模型，包括：

步骤S21，设置公式，用于计算边界粒子j对水粒子i施加的排斥力权值A：

其中v

步骤S22，获取水粒子i到边界粒子j的切向距离x及法向距离y，当q＜1 且x＜Δp时，计算边界粒子j对水粒子i施加排斥力f

f

其中，q＝y/2Δp，Δp为水粒子i与边界粒子j间的初始间距。

在上述技术方案的基础上，进一步地，步骤S3包括：搜索中心粒子影响范围内的相邻粒子，例如选取a、b作为相互作用的两个水粒子(简称两粒子)；将两粒子的接触间断物理量采用黎曼求解器近似，并引入高阶MUSCL格式对HLLC黎曼求解器左右两端的输入值进行重构，得到具有高阶精度的光滑解输入，代入控制方程组得到黎曼解输出值；根据输出值(ρ，u，w，p)

在上述技术方案的基础上，进一步地，定义物理量Φ＝[ρ，u，w]，定义粒子a 与粒子b接触断面处的黎曼求解器输入初值：

其中，

其中，β为限制因子；r

将HLLC黎曼求解器的输入通量定义为Q＝[ρ，ρu，ρw]

根据控制方程组中所定义的迁移速度v

u

其中，Q

引入粒子接触算法，以确定速度投影关系式：

a、b为相互作用的两粒子，x-z为原坐标系，X-T为重构坐标系，u

其中，θ为a、b两粒子连线与x轴的夹角；

得到迁移速度v

u

左右行波及星区波速表达式如下：

其中，γ＝7，p

采用SPH方法对欧拉运动微分方程进行离散，得到基于黎曼解的守恒型控制方程组为：

其中，

优选地，限制因子β＝1.3。

在上述技术方案的基础上，进一步地，步骤S4包括：

采用Symplectic时间步进法进行计算，定义粒子a的物理量为 F

其中，n为时间步，ω为粒子体积，ρ为密度，v为粒子速度，x为粒子位置，Δt为时间步，Δt满足CFL条件，即Δt≤0.4h/c，c为声速，h为光滑长度。

在上述任一技术方案的基础上，进一步地，在步骤S5之后还包括：步骤S6，对单个粒子i的位置求散度，表达式为：

在上述任一技术方案的基础上，进一步地，在步骤S5之后还包括：步骤S7，计算边界的压力值：

综上，本发明所提供的无网格数值模型的模拟方法，构建了改进黎曼解无网格算法的浮体、流体相互作用的数值模型，通过该数值模型，提高了流场压力计算精度，在造波问题中达到了较高的精度，提高了在模拟具有强冲击性的固-液相互作用问题时的有效性及精度。同时，提出了一种适用于改进黎曼解无网格算法的“数值波高仪”后处理方法，用于流体自由表面粒子捕捉及浪高计算；提出了固壁法向排斥力结合一阶核函数修正的固壁处理算法，用于模拟固壁及壁面压力估算，相较于相关技术中的多种无网格方法，充分考虑了水体发生翻卷、破碎、飞溅等强非线性因素对流场的影响，改善了应力不稳定性、粒子压力场震荡、自由面表达精度不高的问题，使得本发明发能适用于模拟非线性波浪与浮体相互作用的场景，具有更加广泛的适用环境，从而在船舶与海洋工程、近海工程等领域均能应用，具有较高的使用价值。

本发明第二方面的技术方案提供了一种计算机设备，包括存储器、处理器和存数在存储器上并可在处理器上运行的计算机程序，处理器执行计算机程序时，能够实现本发明第一方面技术方案中任一项的无网格数值模型的模拟方法的步骤。因此，本发明第二方面的技术方案提供的计算机设备具有上述本发明第一方面任一项技术方案中任一项的全部有益效果，在此不再赘述。

本发明第三方面的技术方案提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被执行时能够实现本发明第一方面技术方案中任一项的无网格数值模型的模拟方法的步骤。因此，本发明第三方面的技术方案提供的计算机可读存储介质具有上述本发明第一方面任一项技术方案的全部有益效果，在此不再赘述。

附图说明

图1为本发明一个实施例中的无网格数值模型的模拟方法的流程框图。

图2为本发明一个实施例中的HLLC黎曼解结构图。

图3为本发明一个实施例中的粒子接触算法示意图。

图4为本发明一个实施例中的捕捉到的表面粒子与原图的自由面的对比示意图。

图5为本发明一个实施例中的法向排斥力模型。

图6为本发明一个实施例中的浮体受力示意图。

图7为本发明一个实施例中的无网格数值模型与传统SPH在同一时刻流场压力云图的对比示意图。

图8为本发明一个实施例中的弱非线性规则波流场压力云图。

图9为本发明一个实施例中的弱非线性规则波浪高对比图。

图10为本发明一个实施例中的楔形体入水压力分布示意图。

图11为本发明一个实施例中的楔形体入水与试验值、理论值的对比图。

图12为本发明一个实施例中的楔形体完全入水对比图。

图13、图14、图15为本发明一些实施例中的浮体与非线性波浪相互作用自由液面变化对比图。

图16为本发明一个实施例中的相邻单元格搜索方法示意图。

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明的保护范围。

目前，现有的一些相关技术提供了一种称为无网格的数值方法，与传统基于网格的CFD计算方法相比，无网格法对于整个计算域上的各个节点，不需要采用预定义的网格信息对其离散，而是通过一系列任意分布的粒子或节点来求解具有不同类型边界条件的偏微分方程组(PDEs)，从而得到精确稳定的数值解，所有粒子或节点不需要网格对其进行连接，各个点在控制方程的作用下可以自由运动，因此特别适用于处理具有极端状态的高压、高速、高能的强冲击、大变形问题。截至目前为止，已相继提出了多种无网格方法，如移动粒子半隐式法(MPS)、无网格局部Petrov-Galerkin法(MLPG)、光滑粒子流体动力学方法(SPH)、重构核粒子法(RKPM)、弥散元法(DEM)等。由于无网格法起步较晚，均存在一定的不足和局限性，需进一步优化和完善。

但是，随着SPH方法在流体力学领域应用的不断深入，其内在缺陷也逐渐暴露出来，如应力不稳定性、粒子压力场震荡、自由面表达精度不高等。为了拓展SPH方法在工业领域的应用，提高算法的计算精度及效率，本发明的一些实施例提供了一种无网格数值模型的模拟方法。

实施例1

如图1所示，本发明的一个实施例提供了一种无网格数值模型的模拟方法，包括：

Step1：初始化流场，对每个流体粒子赋予各物理量初值，如位置、密度、速度、压力。

Step2：布置固壁边界粒子。

Step3：进行NNPS粒子搜索，采用高阶MUSCL格式对粒子接触间断处黎曼问题初值进行重构(ρ，u，w)

Step4：采用时间步进方法，解出上述守恒型变量(ω，ωρ，ωρv)，相除后得到 N+1步粒子的速度、密度、体积，将密度代入状态方程得到压力，结合XSPH法修正速度，更新粒子位置，若达到输出条件即输出数据，之后将时间步向后推进，再次循环Step3。

Step5：若计算未结束，继续循环Step3、Step4。

实施例2

实施例1中的无网格数值模型的模拟方法具体提供了一种改进黎曼解无网格算法的流体控制方程。

与传统SPH方法不同，改进黎曼解无网格数值模型不再将相互作用的两粒子的物理量(密度、压力、速度等)直接代入控制方程进行求解，而是认为粒子间的相互作用存在于接触断面处，在求解过程中引入粒子接触算法的思想，将两粒子的接触间断物理量采用黎曼求解器近似，同时引入具有二阶精度的 TVD-MUSCL格式对HLLC黎曼求解器左右两端的输入值进行重构，得到具有高阶精度的光滑解输入，同时采用一种不降低光滑解区域精度的通量限制器，用以抑制间断附近的数值波动。

采用SPH方法对欧拉运动微分方程进行离散，得到基于黎曼解的守恒型控制方程为：

其中：

带有张量积的项可展开如下：

实施例3

实施例1中的无网格数值模型的模拟方法具体提供了一种HLLC黎曼求解器。

间断处的黎曼精确解求解非常困难，一般采用近似求解。HLLC黎曼求解器为HLL求解器的改进，是一种近似求解黎曼解的方法，将波解由双波结构推演为三波结构，间断处的波解结构如图2所示。

将HLLC黎曼求解器的输入通量定义为Q＝[ρ，ρu，ρw]

u

其中：Q

其中：θ为两粒子连线与x轴的夹角。迁移速度v

u

左右行波及星区波速表达式如下：

其中：γ＝7，p

以上方程中的变量，均为高阶MUSCL格式重构后的量。

实施例4

实施例1中的无网格数值模型的模拟方法具体提供了一种高阶MUSCL格式。

MUSCL为高阶精度的Godunov格式，也被称为限制器插值方法。本实施例将MUSCL格式引入到改进黎曼解无网格数值模型算法中，粒子的物理量不再是恒定的常数分布，而被认为是线性分布，通过线性插值计算，从而得到两粒子接触断面处的值，作为黎曼问题的初始间断值，也可以认为是对黎曼求解器左右两边的输入值进行重构，以满足更高精度的计算要求。

定义物理量Φ＝[ρ，u，w]，粒子a与b接触断面处的黎曼求解器输入初值由下式定义：

其中：

其中：β为限制因子，该值对会对格式所带来的数值粘性产生影响，建议取值为1.3。

实施例5

实施例1中的无网格数值模型的模拟方法具体提供了一种时间步进步骤。

采用Symplectic时间步进法进行计算，定义粒子a的物理量为F

其中：n为时间步，ω为粒子体积，ρ为密度，v为粒子速度，x为粒子位置，Δt为时间步。Δt应满足CFL条件，即Δt≤0.4h/c，c为声速，h为光滑长度。

实施例6

实施例1中的无网格数值模型的模拟方法具体提供了一种自由液面捕捉及“数值波高仪”。

在推板造波问题中，需要在固定测点处设置虚拟的“数值波高仪”，从而得到固定测点处的波高时历变化曲线，本实施例所提出的“数值波高仪”后处理计算方法主要分两步进行：首先对位于自由表面的粒子进行捕捉；之后对得到的自由面粒子进行插值计算，得到与固定测点处最为接近的自由面粒子，其在z方向的坐标值即为该位置的波面高度。由此可知，如何能够精确的捕捉自由面粒子成为 SPH数值造波中的关键技术。

本实施例提出一种高精度的自由面捕捉算法，该算法的核心思想是对单个粒子的位置求散度，表达式如下：

其中：r为粒子位置，

本实施例建议当

实施例7

实施例1中的无网格数值模型的模拟方法具体提供了一种固壁处理步骤。

固壁边界的施加方法一直以来是SPH方法中的重点和难点，选取方法的好坏对计算精度有直接的影响，截至目前为止还未有一种公认合理的方法出现，这成为阻碍SPH方法发展的一个重要原因。目前常用的方法主要有排斥力法、镜像粒子法、虚粒子法，这三种方法各有优缺点，其中排斥力法的边界粒子设置较为简单，可以处理具有复杂形状的问题，而且能够有效的防止粒子穿越边界，其缺点在于临近边界域粒子存在积分域缺失，守恒性较差，而且边界力公式由人工给定，受经验因素影响较大，在计算初期会对流体域产生扰动；镜像粒子法与虚粒子法能提供边界域外的积分粒子，守恒性好但是均难以处理复杂边界问题，同时在计算时容易出现粒子穿透边界的现象。

本实施例的数值模拟系统采用法向排斥力固壁处理方法，该模型认为固壁边界粒子对水粒子只提供法向力作用，而没有切向力，示意图如图5所示。表达式如下：

f

其中：n

实施例8

实施例1中的无网格数值模型的模拟方法具体提供了一种浮体运动控制方程。

本实施例采用法向排斥力模型作为浮体的移动边界条件，当水粒子位置处于边界粒子的影响域范围内时，边界粒子对其沿固壁法向方向施加排斥力，从而防止粒子穿越边界，浮体的受力由边界粒子作用于水粒子的排斥力间接得到，在动量守恒的考虑下，边界粒子给予水粒子的排斥力的同时，水粒子也对边界粒子产生大小相同、方向相反的反作用力，将浮体上各边界粒子的受力进行求和从而得到浮体所受的总力，以此计算浮体的运动，作用方式如图6所示，表达式如下：

其中：j为浮体边界固壁粒子，i为j影响域内的水粒子，N为影响域内的水粒子数量，m为粒子质量，

实施例9

实施例1中的无网格数值模型的模拟方法具体提供了一种通过一阶核函数修正边界压力估算的步骤。

采用法向排斥力模型作为浮体动边界条件时，不能直接得到边界粒子的压力值，因此需对边界粒子影响域内的水粒子压力进行积分插值运算，间接估算边界处压力，但是固壁边界处不可避免的会出现核函数截断，不能满足归一性，故需要进行一阶核函数修正，边界压力计算采用下式：

实施例10

实施例1中的无网格数值模型的模拟方法具体提供了一种相邻单元格搜索方法。

光滑核函数具有一个紧支域，在每个中心计算粒子以κh支持域半径内包含有限个的粒子数目，而所包含的这些粒子将会在粒子近似计算过程中被使用，由于SPH为无网格拉格朗日算法，粒子位置即时变化，不同时刻下粒子支持域内所包含的粒子都有所不同，故在进行积分插值前，首要任务便是确定中心粒子影响范围内的其它粒子，这些粒子通常被称为相关粒子的最近相邻粒子(Nearest Neighboring Particle，简称NNP)，搜索相邻粒子的过程被称为最新相邻粒子搜索 (Nearest Neighboring Particle Search，简称NNPS)。在计算过程中，NNPS的工作量往往要占到整个SPH计算所耗时的60％～80％，因此找到一个高效的搜索方法十分必要，本实施例构建的数值模型采用相邻单元格搜索法(neighboringcell search)。

在光滑长度恒定的情况下，采用相邻单元格搜索法十分有效。核心思想是将整个计算域划分为固定尺寸的网格单元，具体的尺寸由所采用核函数的影响半径κh给定，首先对每个网格单元进行坐标式定义，如图16所示，其中数字1、2、 3.....表示位置坐标定义，计算过程中位置坐标始终保持不变，此时域内的粒子分布在各个单元内，在每一时间步计算初期，通过对粒子位置的计算，将每一粒子标定进对应坐标位置的网格单元，又由于核函数的影响半径与网格尺寸相等，故该粒子影响域内的其它粒子必定只存在于标定网格周边的网格单元以及该粒子本身被标定的网格单元内，在一维、二维、三维空间的搜索范围分别是3、9、 27个单元，在NNPS计算时只需对相邻坐标的单元格内所标定的粒子进行遍历搜索，其本质与全配对搜索类似，但极大的缩小了目标搜索范围，因此极大的缩短了计算时间，相邻单元格搜索方法的复杂度阶数为O(N log N)。

在对计算域内的所有粒子进行单元格位置标定后，其后的搜索步骤为：首先由计算域底部最左边的单元格开始，沿X轴向右依次扫描，之后沿Z轴方向向上推进一行，再向右依次扫描，如此重复，如图13所示，按位置标号依次为： (1、2...6)--(7、8...12)...--(19、20...24)。在此过程中仅对每个单元格本身，以及与其相邻的东(E)、北(N)、东北(NE)、西北(NW)四个方向的单元格内的粒子进行计算，用以避免粒子相互作用的重复性计算，如图16所示，以图中23号单元格为例，仅需对24(E)、31(NW)、32(N)、33(NE)上述四个方向的单元格及其自身内部的粒子进行关联计算，因为与其相邻的13、14、15、22号单元格，在中心单元格为其自身时，已经与23号单元格关联完毕。

综上，本发明所提供的无网格数值模型的模拟方法至少具有以下有益效果：

(1)提高流场压力计算精度

以具有强冲击、大变形的溃坝算例为例，分别给出传统SPH及改进黎曼解无网格数值模型在同一时刻流场压力云图对比，如图7所示，可以看出，改进黎曼解无网格数值模型相较传统SPH算法，流场压力计算精度有质的提升。

(2)提高数值波高仪算法有效性

为了考察改进黎曼解无网格数值模型“数值波高仪”的计算有效性，模拟生成非线性规则波，与二阶斯托克斯波(Stokes)理论值进行对比，图8为流场压力云图，图9为x＝8m处捕捉的浪高值对比图，可以看出，本发明提出的“数值波高仪”能够以较高的精度对自由表面粒子进行捕获并计算出该位置的浪高变化过程，同时，改进黎曼解无网格数值模型算法模型在造波问题中达到了较高的精度，而传统SPH面临波浪衰减的内在问题。

(3)提高流固耦合计算精度

本发明以两个经典流固耦合算例进行算法精度的验证，分别为物体入水及浮体与波浪的相互作用。

(3.1)能够用于描述楔形体入水问题

入水问题异常复杂，涉及到流、固、气三相耦合等强非线性问题，同时砰击的过程中会产生剧烈的自由面变形，这些都导致传统的网格方法难以处理该类问题。本发明构建改进黎曼解无网格数值模型对二维楔形体入水进行模拟，首先针对入水初期的砰击过程进行计算，研究了楔形体在此过程中的受力情况及物面压力分布，并与相关的理论值及试验值进行对比，以期验证耦合边界处理方法在模拟具有强冲击性的固-液相互作用问题时的有效性及精度。

图10给出楔形体入水流场压力云图，且图11给出了与理论解及试验值的对比，改进黎曼解无网格数值模型与理论解吻合较好，与试验值的误差主要来自于未考虑空气的影响。

图12给出楔形体完全入水的自由液面形态(右图)，并与试验高频照片(左图)进行了对比，吻合度较高。

(3.2)能够用于描述浮体与波浪的相互作用

流体和浮体的相互作用是一个复杂的问题，在船舶与海洋工程、近海工程等领域都是非常重要的研究方向。本发明构建改进黎曼解无网格数值模型对非线性波浪与浮体相互作用进行模拟，并与试验值进行对比，验证算法精度及有效性。

图13、图14、图15给出了各个时刻的改进黎曼解无网格数值模型模拟结果，并与试验高频照片进行了对比，可以看出，各个时刻的自由液面吻合情况良好，算法精度较高。

本发明的一个实施例提供了一种计算机设备，包括存储器、处理器和存数在存储器上并可在处理器上运行的计算机程序，处理器执行计算机程序时，能够实现上述实施例中任一项的无网格数值模型的模拟方法的步骤。

本发明的一个实施例提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被执行时能够实现上述实施例中任一项的无网格数值模型的模拟方法的步骤。

需要说明的是，本领域技术人员可以理解上述实施例的各种方法中的全部或部分步骤是可以通过程序来指令相关的硬件完成的，改程序可存储于计算机的可读存储介质中，存储介质包括只读存储器(ROM，Read-Only Memory)、随机存储器(RAM，Random AccessMemory)、可编程只读存储器(PROM，Programmable Read-Only Memory)、只读光盘(CD-ROM，Compact Disc Read-Only Memory) 或其他光盘存储器、磁带存储器，或者能够用于携带或存储数据的计算机刻度的任何其他介质。

在本说明书的描述中，具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

当然，本发明还可以有其他多种实施例，在不背离本发明精神及其实质的情况下，熟悉本领域的技术人员可以根据本发明作出各种相应的改变和变形，但这些改变和变形都应属于本发明的权利要求的保护范围。