CN113761812B - 基于无网格法的单层壁面复杂几何流域的求解方法及系统 - Google Patents

Abstract

本发明公开了一种基于无网格法的单层壁面复杂几何流域的求解方法及系统，通过将流体粒子进行划分为近壁面流体粒子与内部流体粒子，结合光滑通用壁面模型与无表面粒子判定技术，采用虚拟粒子自动补偿近壁面流体粒子数密度，接触模型提供壁面作用力，利用二阶龙格库塔法解决压缩不稳定性，判断运行时刻是否超出了终止时刻，如果超出，输出瞬时粒子位置、压力和速度参数，并绘制近壁面流体粒子的流线、计算图以及压力、速度分布图。本发明通过采用无网格法，流体粒子依据粘性流体的运动方程进行流动，并与周围的流体和固体粒子建立力学关系，用于模拟计算区域存在大变形的流动问题，对于现代复杂工程问题数值仿真技术的提升有重大意义。

Description

基于无网格法的单层壁面复杂几何流域的求解方法及系统

技术领域

本发明属于仿真技术领域，涉及一种基于无网格法的单层壁面复杂几何流域的求解方法及系统。

背景技术

随着电子计算机的迅速发展，数值模拟技术在工程领域的应用越来越广泛。目前主流仿真技术如有限容积法等，需要通过网格剖分的方式计算流场，对于流通域形状剧烈变化的问题，常采用网格调整和重构技术来处理，但过程复杂且常面临网格质量下降和收敛性变差的风险。无网格粒子法如移动粒子半隐式法基于流体质点概念，在拉格朗日框架下采用没有固定拓扑关系的离散点代替网格及节点，避免了所有基于网格的限制，适用于计算任意非定常大变形流动。

然而，算法的精度、稳定性以及边界条件的适用性问题一直是无网格粒子法应用于复杂工程问题内流计算的难点；原始算法由于表面粒子误判问题容易导致较大误差，使得计算发散；传统的壁面边界条件难以用于存在复杂曲面的模型，因为大曲率表面附近的粒子数密度误差较大，而且大量参与压力计算的壁面粒子会明显降低计算效率；目前已有的改进表面粒子判定法有数量判断法、光照法、几何判定法，但以上方法计算复杂，对于存在复杂流动的三维问题，仍较难完全消除表面粒子误判。国际上主流的粒子法壁面边界可分为两类，一是基于壁面权函数的多边形壁面模型，另一类基于边界力，例如统一半解析壁面(USAW)模型；在多边形壁面模型中，网格的质量直接影响模型的精度，对于复杂的几何体，粒子数密度补偿仍然是一个问题；基于边界力的壁面模型受力的数学处理较复杂，目前大多应用于二维问题。

发明内容

本发明的目的在于解决现有技术中的问题，提供一种基于无网格法的单层壁面复杂几何流域的求解方法及系统，通过将流体粒子进行划分为近壁面流体粒子与内部流体粒子，结合光滑通用壁面模型与无表面粒子判定技术，采用虚拟粒子自动补偿近壁面流体粒子数密度，壁面粒子仅有单层且不参与压力计算，壁面作用力由接触模型来提供，压缩不稳定性利用二阶龙格库塔法解决，判断运行时刻是否超出了终止时刻，如果超出，输出瞬时粒子位置、压力和速度参数，并绘制近壁面流体粒子的流线、计算图以及压力、速度分布图。

为达到上述目的，本发明采用以下技术方案予以实现：

基于无网格法的单层壁面复杂几何流域的求解方法，包括：

S1：对参数和变量进行初始化设定，配置初始环境；

S2：以流体粒子到最近的壁面粒子的法向量投影距离为依据，区分内部流体粒子和近壁面流体粒子；

S3：对流体粒子进行拉普拉斯算子计算，并对近壁面流体粒子进行补偿；

S4：判断当前时间步是否为虚拟时间步，若是，执行S5；若否，则初始化虚拟时间步内的变量，并进入虚拟时间步，重复S2、S3和S4；

S5：将虚拟时间步的压力梯度值与当前时间步的压力梯度值进行计算以修正速度参数，以修正后的速度参数更新位移参数；

S6：判断当前时刻是否超出终止时刻，若是，执行S7；若否，输出瞬时流体粒子位置、压力和速度参数，重复S2、S3、S4、S5和S6；

S7：处理瞬时参数，得到流线、计算图以及压力、速度分布图。

本发明的进一步改进在于：

对参数和变量进行初始化设定，配置初始环境具体为：参数包括流体物性、接触模型系数和重力加速度，变量包括速度、压力以及时间；

配置初始环境包括特殊壁面粒子的标定和流体粒子的初始化数据；其中，初始数据包括初始流体粒子坐标、壁面粒子坐标以及壁面粒子法向量；特殊壁面粒子是指在尖角和拐点处不存在法向量的壁面粒子。

以流体粒子到最近壁面粒子的法向量投影距离为依据，区分内部流体粒子和近壁面流体粒子，具体为：

以流体粒子到最近壁面粒子的法向量投影距离为判断依据，投影距离由公式(1)计算得到，判据如公式(2)所示；

其中，|r_iw|为投影距离；|r_ij|为流体粒子i与最近壁面粒子j距离；|n_j|为最近壁面粒子法向量；l₀为流体粒子直径；S为特殊壁面粒子集合；N₁、N₂和N₃分别为近壁面、次近壁面和内部流体粒子集合。

对近壁面流体粒子进行补偿具体为：利用壁面函数进行补偿，粘性补偿数值由公式(3)计算：

其中，为粘性项补偿值；d、λ、n₀、l₀分别为维度、粘性修正系数、粒子数密度常数和粒子直径；a为粘性补偿修正系数。

补偿近壁面流体粒子还包括：通过显式计算得到近壁面流体粒子的速度；采用无表面粒子判定技术计算近壁面流体粒子的密度；根据接触模型计算壁面对近壁面流体粒子的作用力；利用无表面粒子判定技术计算压力泊松方程系数；基于不完全Cholesky分解的共轭梯度法求解压力泊松方程并求解压力梯度，在压力泊松方程中引入速度散度，提高近壁面粒子的承压性；并通过近壁面压力补偿得到近壁面流体粒子的压力梯度。

采用无表面粒子判定技术计算近壁面流体粒子数密度的计算式如公式(4)所示：

其中，<n">_i为补全后粒子数密度；〈n>_i为原始方法计算的粒子数密度；〈n'>_i为考虑虚拟粒子后的粒子数密度；Nr为发生挤压的粒子数个；无表面粒子判定技术通过在粒子数密度缺失的自由表面或近壁面流体粒子附近产生虚拟粒子，从而补足其粒子数密度；

所述接触模型由弹簧和阻尼器组成，在壁面法向和切向两个方位各有一个弹簧和阻尼器；当近壁面流体粒子即将穿透壁面时，壁面通过接触模型向该类流体粒子施加斥力，防止其穿透；壁面对近壁面流体粒子的作用力的计算式如公式(5)所示：

其中，f_iw为壁面作用力；上标n和t表示法向与切向；k、ζ和μ为弹簧、阻尼和摩擦系数；

所述利用无表面粒子判定技术计算压力泊松方程系数的计算式如公式(6)所示：

其中，M为流体粒子个数；为显式计算后的速度；d为维度；λ为修正系数；n⁰为粒子数密度常数；n"为补全后粒子数密度；γ为伪可压系数；

所述近壁面压力补偿值由公式(7)得到：

其中，为近壁面流体粒子的压力梯度补偿值；m为流体粒子的质量；dr为近壁面流体粒子到平衡面的距离。

S5中修正速度具体为：基于二阶龙格库塔的修正速度计算公式(8)所示：

其中，为t时间步的压力梯度；/>为虚拟时间步的压力梯度，u^t为t时间步的修正速度，u^t+1为t+1时间步的修正速度。

进入虚拟时间步为当前时间步加上虚拟时间步长；虚拟时间步为当前时间步的下一个时间步；虚拟时间步长为当前时间步进入虚拟时间步的推进步长。

处理瞬时参数，得到流线、计算图以及压力、速度分布图，具体为：压力和速度分布图为将压力、速度数据导入TECPLOT软件生成，将速度矢量数据导入PARAVIEW软件即生成流线图；将多个时刻某一流体微团的位置坐标导入TECPLOT软件即生成该微团的流动迹线图。

基于无网格法的单层壁面复杂几何流域的求解系统，包括：

设定模块，所述设定模块用于对参数和变量进行初始化设定，配置初始环境；

划分模块，所述划分模块用于以流体粒子到最近的壁面粒子的法向量投影距离为依据，区分内部流体粒子和近壁面流体粒子；

补偿模块，所述补偿模块用于对流体粒子进行拉普拉斯算子计算，并对近壁面流体粒子进行补偿；

第一判断模块，所述第一判断模块用于判断当前时间步是否为虚拟时间步，若否，则初始化虚拟时间步的变量，并进入虚拟时间步；

修正模块，所述修正模块用于将虚拟时间步的压力梯度值与当前时间步长的压力梯度值进行计算以修正速度参数，以修正后的速度参数更新位移参数；

第二判断模块，所述第二判断模块用于判断当前时刻是否超出终止时刻，若否，输出瞬时流体粒子位置、压力和速度参数；

图像绘制模块，所述图像绘制模块用于处理瞬时参数，得到流线、计算图以及压力、速度分布图。

与现有技术相比，本发明具有以下有益效果：

本发明提出了一种基于无网格法的单层壁面复杂几何流域的求解方法及系统，与传统基于网格剖分的方法不同，本发明通过采用无网格粒子法，计算节点不再是有固定编号和固定拓扑结构的网格节点，而是没有相互固定关联的流体粒子，流体粒子依据粘性流体的运动方程流动到任何地方，再与当地周围的流体和固体粒子建立力学关系，特别适用于模拟计算区域存在大变形的流动问题，如扭曲破碎、融合等。因此，无网格法在处理包含动静相干现象的复杂问题时，流域是整体离散和整体求解，没有动静交界面的任何假设处理和信息交换；无需分块处理流体域，对于流体机械内部形状复杂的流通域，均可利用粒子的自适应性填充布置，根据近壁面流体粒子的参数输出流体微团运动轨迹，根据输出的运动迹线研究流动流体粒子的运动机理。

附图说明

为了更清楚的说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1为本发明的基于无网格法的单层壁面复杂几何流域的求解方法流程图；

图2为初始的壁面粒子数据图；

图3为本发明计算出的复杂静水蘑菇算例的压力分布图；

其中，a为本发明计算出的复杂静水蘑菇算例的压力分布三维图；

b为本发明计算出的复杂静水蘑菇算例的压力分布前视图；

图4为复杂静水蘑菇算例底部压力数值随时间变化图；

图5为复杂静水蘑菇算例水平面高度随时间变化图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。

因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。

在本发明实施例的描述中，需要说明的是，若出现术语“上”、“下”、“水平”、“内”等指示的方位或位置关系为基于附图所示的方位或位置关系，或者是该发明产品使用时惯常摆放的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”等仅用于区分描述，而不能理解为指示或暗示相对重要性。

此外，若出现术语“水平”，并不表示要求部件绝对水平，而是可以稍微倾斜。如“水平”仅仅是指其方向相对“竖直”而言更加水平，并不是表示该结构一定要完全水平，而是可以稍微倾斜。

在本发明实施例的描述中，还需要说明的是，除非另有明确的规定和限定，若出现术语“设置”、“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。

下面结合附图对本发明做进一步详细描述：

参见图1，图1公布了一种基于无网格法的单层壁面复杂几何流域的求解方法，包括：

步骤1：对参数和变量进行初始化设定，配置初始环境。

对参数和变量进行初始化设定，配置初始环境具体为：参数包括流体物性、接触模型系数和重力加速度，变量包括速度、压力以及时间，将速度、压力和时间变量置0；配置初始环境包括特殊壁面粒子的标定和流体粒子的初始化数据；其中，初始数据包括初始流体粒子坐标、壁面粒子坐标以及壁面粒子法向量；特殊壁面粒子是指在尖角和拐点处不存在法向量的壁面粒子。接触模型系数包括弹簧刚度、阻尼系数和摩擦系数。

将初始流体粒子坐标、高精度壁面数据(图2)、壁面法向量数据和特殊壁面粒子数据导入程序用于计算。

步骤2：以流体粒子到最近的壁面粒子的法向量投影距离为依据，区分内部流体粒子和近壁面流体粒子；对流体粒子进行拉普拉斯算子计算，并对近壁面流体粒子进行补偿。

以流体粒子到最近壁面粒子的法向量投影距离为判断依据，投影距离由公式(1)计算得到，判据如公式(2)所示；

其中，|r_iw|为投影距离；|r_ij|为流体粒子i与最近壁面粒子j距离；|n_j|为最近壁面粒子法向量；l₀为流体粒子直径；S为特殊壁面粒子集合；N₁、N₂和N₃分别为近壁面、次近壁面和内部流体粒子集合。

计算流体粒子粘性项算子并补偿近壁面粘性项；其中，对近壁面流体粒子进行补偿具体为：利用壁面函数进行补偿，粘性补偿数值由公式(3)计算：

其中，为粘性项补偿值；d、λ、n₀、l₀分别为维度、粘性修正系数、粒子数密度常数和粒子直径；a为粘性补偿修正系数。

补偿近壁面流体粒子还包括：通过显式计算得到近壁面流体粒子的速度；采用无表面粒子判定技术计算近壁面流体粒子数密度；根据接触模型计算壁面对近壁面流体粒子的作用力。

采用无表面粒子判定技术计算近壁面流体粒子数密度的计算式如公式(4)所示：

其中，〈n">_i为补全后粒子数密度；<n>_i为原始方法计算的粒子数密度；<n'>_i为考虑虚拟粒子后的粒子数密度；Nr为发生挤压的粒子数个；无表面粒子判定技术通过在粒子数密度缺失的自由表面或近壁面流体粒子附近产生虚拟粒子，从而补足其粒子数密度；

接触模型由弹簧和阻尼器组成，在壁面法向和切向两个方位各有一个弹簧和阻尼器；当近壁面流体粒子即将穿透壁面时，壁面通过接触模型向该类流体粒子施加斥力，防止其穿透；弹簧系统用于提供排斥力，阻尼系统持续衰减被弹开流体粒子的速度，防止被弹开的粒子产生过大的速度而对压力场施加非物理振荡；壁面对近壁面流体粒子的作用力的计算式如公式(5)所示：

其中，f_iw为壁面作用力；上标n和t表示法向与切向；k、ζ和μ为弹簧、阻尼和摩擦系数；

步骤3，利用基于不完全Cholesky分解的共轭梯度法求解压力泊松方程；计算压力梯度；补偿近壁面压力梯度数值；

利用无表面粒子判定技术计算压力泊松方程系数；基于不完全Cholesky分解的共轭梯度法求解压力泊松方程并求解压力梯度，在压力泊松方程中引入速度散度，提高近壁面粒子的承压性；并通过近壁面压力补偿得到近壁面流体粒子的压力梯度。

利用无表面粒子判定技术计算压力泊松方程系数的计算式如公式(6)所示：

其中，M为流体粒子个数；为显式计算后的速度；d为维度；λ为修正系数；n⁰为粒子数密度常数；n"为补全后粒子数密度；γ为伪可压系数；

近壁面压力补偿值由公式(7)得到：

其中，为近壁面流体粒子的压力梯度补偿值；m为流体粒子的质量；dr为近壁面流体粒子到平衡面的距离，平衡面与壁面距离为0.5l₀。

步骤4：判断当前时刻是否虚拟时间步，若否，则初始化虚拟时间步的变量，并进入虚拟时间步，重复步骤2、步骤3和步骤4；若是，将计算出的虚拟时间步的压力梯度值与当前时间步的压力梯度值按公式(8)计算修正速度，以此速度更新位移。

进入虚拟时间步为当前时间步加上虚拟时间步长Δt；虚拟时间步为当前时间步的下一个时间步；虚拟时间步长为当前时间步为了进入虚拟时间步的推进步长。

修正速度具体为：基于二阶龙格库塔的修正速度计算公式(8)所示：

其中，为t时间步长的压力梯度；/>为虚拟时间步的压力梯度，u^t为t时间步长的修正速度，u^t+1为t+1时间步长的修正速度。

步骤5：判断当前时刻是否超出终止时刻，若否，输出瞬时粒子位置、压力和速度等参数，以Δt'为时间步长，推进时间层，重复步骤2、步骤3、步骤4和步骤5，若是，计算出的压力场如图3所示，图3为本发明计算出的复杂静水蘑菇算例的压力分布图；底部压力值与理论解对比结果如图4所示，图4为复杂静水蘑菇算例底部压力数值随时间变化图；自由面位置与理论解的对比如图5，图5为复杂静水蘑菇算例水平面高度随时间变化图。

处理瞬时参数，得到流线、计算图以及压力、速度分布图，具体为：压力和速度分布图为将压力、速度数据导入TECPLOT软件生成，将速度矢量数据导入PARAVIEW软件即生成流线图；将多个时刻某一流体微团的位置坐标导入TECPLOT软件即生成该微团的流动迹线图。

基于无网格法的单层壁面复杂几何流域的求解系统，包括：

设定模块，所述设定模块用于对参数和变量进行初始化设定，配置初始环境；

划分模块，所述划分模块用于以流体粒子到最近的壁面粒子的法向量投影距离为依据，区分内部流体粒子和近壁面流体粒子；

补偿模块，所述补偿模块用于对流体粒子进行拉普拉斯算子计算，并对近壁面流体粒子进行补偿；

第一判断模块，所述第一判断模块用于判断当前时间步是否为虚拟时间步，若否，则初始化虚拟时间步的变量，并进入虚拟时间步；

修正模块，所述修正模块用于将虚拟时间步的压力梯度值与当前时间步长的压力梯度值进行计算以修正速度参数，以修正后的速度参数更新位移参数；

第二判断模块，所述第二判断模块用于判断当前时刻是否超出终止时刻，若否，输出瞬时流体粒子位置、压力和速度参数；

图像绘制模块，所述图像绘制模块用于处理瞬时参数，得到流线、计算图以及压力、速度分布图。

以上仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.基于无网格法的单层壁面复杂几何流域的求解方法，其特征在于，包括：

S1：对参数和变量进行初始化设定，配置初始环境；

S2：以流体粒子到最近的壁面粒子的法向量投影距离为依据，区分内部流体粒子和近壁面流体粒子；

S3：对流体粒子进行拉普拉斯算子计算，并对近壁面流体粒子进行补偿；

S4：判断当前时间步是否为虚拟时间步，若是，执行S5；若否，则初始化虚拟时间步内的变量，并进入虚拟时间步，重复S2、S3和S4；

S5：将虚拟时间步的压力梯度值与当前时间步的压力梯度值进行计算以修正速度参数，以修正后的速度参数更新位移参数；

S6：判断当前时刻是否超出终止时刻，若是，执行S7；若否，输出瞬时流体粒子位置、压力和速度参数，重复S2、S3、S4、S5和S6；

S7：处理瞬时参数，得到流线、计算图以及压力、速度分布图。

2.根据权利要求1所述的基于无网格法的单层壁面复杂几何流域的求解方法，其特征在于，所述对参数和变量进行初始化设定，配置初始环境具体为：参数包括流体物性、接触模型系数和重力加速度，变量包括速度、压力以及时间；配置初始环境包括特殊壁面粒子的标定和流体粒子的初始化数据；其中，初始数据包括初始流体粒子坐标、壁面粒子坐标以及壁面粒子法向量；特殊壁面粒子是指在尖角和拐点处不存在法向量的壁面粒子。

3.根据权利要求1所述的基于无网格法的单层壁面复杂几何流域的求解方法，其特征在于，所述以流体粒子到最近壁面粒子的法向量投影距离为依据，区分内部流体粒子和近壁面流体粒子，具体为：

以流体粒子到最近壁面粒子的法向量投影距离为判断依据，投影距离由公式(1)计算得到，判据如公式(2)所示；

其中，|r_iw|为投影距离；|r_ij|为流体粒子i与最近壁面粒子j距离；|n_j|为最近壁面粒子法向量；l₀为流体粒子直径；S为特殊壁面粒子集合；N₁、N₂和N₃分别为近壁面、次近壁面和内部流体粒子集合。

4.根据权利要求1所述的基于无网格法的单层壁面复杂几何流域的求解方法，其特征在于，所述对近壁面流体粒子进行补偿具体为：利用壁面函数进行补偿，粘性补偿数值由公式(3)计算：

其中，为粘性项补偿值；d、λ、n₀、l₀分别为维度、粘性修正系数、粒子数密度常数和粒子直径；a为粘性补偿修正系数。

5.根据权利要求4所述的基于无网格法的单层壁面复杂几何流域的求解方法，其特征在于，补偿近壁面流体粒子还包括：通过显式计算得到近壁面流体粒子的速度；采用无表面粒子判定技术计算近壁面流体粒子的密度；根据接触模型计算壁面对近壁面流体粒子的作用力；利用无表面粒子判定技术计算压力泊松方程系数；基于不完全Cholesky分解的共轭梯度法求解压力泊松方程并求解压力梯度，在压力泊松方程中引入速度散度，提高近壁面粒子的承压性；并通过近壁面压力补偿得到近壁面流体粒子的压力梯度。

6.根据权利要求5所述的基于无网格法的单层壁面复杂几何流域的求解方法，其特征在于，所述采用无表面粒子判定技术计算近壁面流体粒子数密度的计算式如公式(4)所示：

其中，<n">_i为补全后粒子数密度；<n>_i为原始方法计算的粒子数密度；<n'>_i为考虑虚拟粒子后的粒子数密度；Nr为发生挤压的粒子数个；无表面粒子判定技术通过在粒子数密度缺失的自由表面或近壁面流体粒子附近产生虚拟粒子，从而补足其粒子数密度；

所述接触模型由弹簧和阻尼器组成，在壁面法向和切向两个方位各有一个弹簧和阻尼器；当近壁面流体粒子即将穿透壁面时，壁面通过接触模型向该类流体粒子施加斥力，防止其穿透；壁面对近壁面流体粒子的作用力的计算式如公式(5)所示：

其中，f_iw为壁面作用力；上标n和t表示法向与切向；k、ζ和μ为弹簧、阻尼和摩擦系数；

所述利用无表面粒子判定技术计算压力泊松方程系数的计算式如公式(6)所示：

其中，M为流体粒子个数；为显式计算后的速度；d为维度；λ为修正系数；n⁰为粒子数密度常数；n"为补全后粒子数密度；γ为伪可压系数；

所述近壁面压力补偿值由公式(7)得到：

其中，为近壁面流体粒子的压力梯度补偿值；m为流体粒子的质量；dr为近壁面流体粒子到平衡面的距离。

7.根据权利要求1所述的基于无网格法的单层壁面复杂几何流域的求解方法，其特征在于，所述S5中修正速度具体为：基于二阶龙格库塔的修正速度计算公式(8)所示：

其中，为t时间步的压力梯度；/>为虚拟时间步的压力梯度，u^t为t时间步的修正速度，u^t+1为t+1时间步的修正速度。

8.根据权利要求1所述的基于无网格法的单层壁面复杂几何流域的求解方法，其特征在于，所述进入虚拟时间步为当前时间步加上虚拟时间步长；虚拟时间步为当前时间步的下一个时间步；虚拟时间步长为当前时间步进入虚拟时间步的推进步长。

9.根据权利要求1所述的基于无网格法的单层壁面复杂几何流域的求解方法，其特征在于，所述处理瞬时参数，得到流线、计算图以及压力、速度分布图，具体为：压力和速度分布图为将压力、速度数据导入TECPLOT软件生成，将速度矢量数据导入PARAVIEW软件即生成流线图；将多个时刻某一流体微团的位置坐标导入TECPLOT软件即生成该微团的流动迹线图。

10.基于无网格法的单层壁面复杂几何流域的求解系统，其特征在于，包括：

设定模块，所述设定模块用于对参数和变量进行初始化设定，配置初始环境；

划分模块，所述划分模块用于以流体粒子到最近的壁面粒子的法向量投影距离为依据，区分内部流体粒子和近壁面流体粒子；

补偿模块，所述补偿模块用于对流体粒子进行拉普拉斯算子计算，并对近壁面流体粒子进行补偿；

第一判断模块，所述第一判断模块用于判断当前时间步是否为虚拟时间步，若否，则初始化虚拟时间步的变量，并进入虚拟时间步；

修正模块，所述修正模块用于将虚拟时间步的压力梯度值与当前时间步长的压力梯度值进行计算以修正速度参数，以修正后的速度参数更新位移参数；

第二判断模块，所述第二判断模块用于判断当前时刻是否超出终止时刻，若否，输出瞬时流体粒子位置、压力和速度参数；

图像绘制模块，所述图像绘制模块用于处理瞬时参数，得到流线、计算图以及压力、速度分布图。