一种移动粒子半隐式法的对称边界条件的构建方法

摘要

本发明公开了一种移动粒子半隐式法的对称边界条件的构建方法，针对移动粒子半隐式法基于高斯分布的加权平均思想求取物理量的计算特点，本发明利用了“镜像法”对靠近对称边界的粒子进行物理量的数值补偿，以保证对称边界附近的粒子在动力学和热力学求解上的正确性；考虑到某些三维模拟对象具有一组正交的对称面，希望在计算的初始布置时只保留其正交对称面的一个象限内的流域(即1/4流域)，以达到最大程度的计算简化，该“镜像法”对有一组正交对称面的算例依然使用。本发明是对移动粒子半隐式法各类边界条件研究的补充，保证计算三维对称问题时可以简化计算域，降低多余的计算消耗，提高计算效率。

说明书

技术领域

本发明属于计算流体动力学粒子法技术领域，涉及一种移动粒子半隐式法的对称边界条件的构建方法。

背景技术

自计算流体动力学(CFD)诞生以来，以欧拉视角对控制方程进行离散的方法便成为数值研究的主要研究方法并一直发展至今。而基于欧拉框架下以网格节点为离散单元的控制体对控制方程(动量方程和热平衡方程)进行离散时，对流项的存在无疑会引入额外的计算误差，这就造成了计算结果相比理论结果会产生数值耗散或是数值震荡。此外，固定于计算域内的网格也为计算带来了较大的限制，一方面是因为网格划分的质量极大地影响着计算的收敛性，另一方面，当模拟大变形流动问题时，固定网格的对流动的适应性很难满足准确模拟流动变形的过程，尽管随后产生了动网格技术，但依然只适合于计算按照一定可预见规律运动的流体。

与欧拉法相对应的拉格朗日法则是以流体单元(微团)为研究对象。这种研究方法的优点在于避免和对流项的计算，物理量的变化只和时间有关。而在流场的离散上，拉格朗日计算流体动力学以散点代替网格节点，两者的区别在于，这些散点之间只存在稳定的间距，而不存在固定的拓扑结构，在计算过程中，这些散点的空间位置会根据控制方程的计算而不断更新。由于散点的两两间距离固定且统一(针对不可压缩流体)，而整体上又随流场的流动而运动，类似若干光滑的粒子在流动，因此对应于欧拉框架下的“网格法”，拉格朗日计算流动动力学方法亦有“粒子法”这一称呼。

基于宏观流体力学基本控制方程(N-S方程)的粒子法以光滑流体动力学 (SPH)法和移动粒子半隐式(MPS)法为主。这两种方法都在模拟大变形流动和自由表面捕捉上具有很大优势。其中，由Koshizuka和Oka于1986年提出的移动粒子半隐式法主要模拟不可压缩流体。该方法对速度等物理量进行显示求解，并依据不可压缩性建立压力泊松方程，对压力进行隐式求解。

在进行MPS法数值模拟时，每一个流体粒子的物理量是对以该粒子为圆心、影响距离为半径的圆形影响域内的其他流体粒子的物理量进行加权平均得到的。粒子与粒子之间不建立稳定的拓扑关系，而是在每一个时间步内对所有流体粒子进行遍历搜索，确定粒子间的瞬时位置关系，这也是MPS法易于模拟大变形流动的基本条件。然而，大量的遍历搜索计算十分消耗计算资源，同时，影响域内粒子的加权平均算法也使得MPS法在边界设置上具有一定难度。

优化边界条件一直以来是MPS法重点关注的问题之一。对称边界条件是数值计算常用的一种边界条件，这种计算可以简化具有对称性的流场，通过对计算区域沿对称面进行截断，可以有效降低计算消耗。但针对MPS法的对称边界条件的设置方法尚不广泛。

发明内容

移动粒子半隐式法的诸多模型都涉及粒子数密度，同时，粒子数密度也是判断粒子是否为表面粒子的准则。在压力求解时，需要依据该准则将非表面流体粒子筛选出来建立压力泊松方程。当依据对称边界条件对流域进行截断时，对称面上的粒子就会被错误地判断为表面粒子，这样既不能正确计算压力，也会因为其影响域内缺少足够的含有物理量的流体粒子作为加权计算的参考，导致这些粒子的物理量计算失准。本发明的目的在于解决现有技术中的问题，提供一种移动粒子半隐式法的对称边界条件的构建方法

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

一种移动粒子半隐式法的对称边界条件的构建方法，包括以下步骤：

步骤1，统计计算对象的对称面影响域内的流体粒子数，为存储临时粒子的数组分配空间；

步骤2，创造临时镜像粒子，临时镜像粒子的空间坐标与其对应的流体粒子根据对称面对称，物理量相同；

步骤3，通过求解宏观流体力学基本控制方程中的粘性项和体积力项，进行速度和坐标临时量的显示计算；

步骤4，进行表面张力的计算并进行镜像粒子位置的更新；

步骤5，进行压力泊松方程的建立和求解；

步骤6，给镜像粒子按照镜像关系赋压力值；

步骤7，通过求解压力梯度对速度临时量进行修正；

步骤8，释放存储临时粒子的数组空间。

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

所述步骤1中，对称面的影响域范围指的是对称面向流体一侧4.1l

所述步骤2中，镜像粒子所具有的标量物理量与对应的流体粒子相等；矢量物理量在与对称面平行的面上的分量对应相等，在对称面法向上的分量互为相反数。

所述步骤3中，粘性项中的拉普拉斯项只计算流体粒子，但每一个流体粒子进行遍历搜索时，需要包括镜像粒子；在计算完成后，需要对镜像粒子的位置和速度进行更新。

所述步骤4中，表面张力计算分为流体粒子间的受力计算和流体粒子与界面之间的受力计算，流体粒子与界面之间的受力计算体现的是液体在固体表面的浸润性，计算模型为基于杨氏方程的界面浸润模型；按照式(1)计算流体粒子间的受力，按照式(2)计算流体粒子与固体粒子间的受力，最后按照式(3)将两个力进行加和：

所述步骤6中，镜像粒子的压力按照镜像关系从与之对应的流体粒子获得。

所述步骤7中，在速度修正后对镜像粒子的空间和速度按照和流体粒子的镜像关系进行调整。

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

本发明提出一种适用于三维含表面张力的MPS法的对称边界条件——镜像粒子对称补偿法，该方法依据对称条件对接近称面的流体粒子进行物理量的补偿，而不是被当作表面粒子处理，保证即使截断流场，计算也能够正确进行。

附图说明

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

图1为本发明算法的计算流程；

图2为依据镜像粒子对靠近对称边界的流体粒子进行数值补偿的示意图；

图3为1/2方形流域和完整方形流域的核函数对比；

图4展示了图2方形流域中y＝0和y＝1两条数轴上的粒子数密度值，已进行对比验证；

图5为四分之一液滴、二分之一液滴和完整液滴润湿固壁表面的过程；

图6为四分之一液滴、二分之一液滴和完整液滴浸润过程中液滴的高度和初始高度之比，横轴表示时间；

图7为液滴接触角模拟，图中理论接触角为120°，测量X-Z平面和Y-Z平面两个角度取平均值，表征最终的模拟结果；

图8为三种算例的计算CPU时间对比，横轴表示时间，以展示减小粒子数对计算消耗的减少。

具体实施方式

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

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

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

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

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

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

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

参见图1，本发明移动粒子半隐式法的对称边界条件的构建方法，包括以下步骤：

步骤1，统计计算对象的对称面影响域内的流体粒子数，为存储临时粒子的数组分配空间。

对称面的影响域范围指的是对称面向流体一侧4.1l

步骤2，创造临时镜像粒子，临时镜像粒子的空间坐标与其对应的流体粒子根据对称面对称，物理量相同。

因为布置的临时镜像粒子和与之对应的流体粒子在空间位置上是根据对称面是对称的，这些镜像粒子所具有的标量物理量与对应的流体粒子相等。矢量物理量在与对称面平行的面上的分量同样对应相等，而在对称面法向上的分量互为相反数。

步骤3，进行速度和坐标临时量的显示求解。

镜像粒子不参与作为体积力的重力项计算(若有其他体积力项，同样不需要考虑镜像粒子)，粘性项中的拉普拉斯项同样只计算流体粒子，但每一个流体粒子进行遍历搜索时，需要包括镜像粒子。在该步计算完成后，需要对镜像粒子的位置和速度进行更新。

步骤4，进行表面张力的计算并进行镜像粒子位置的更新。

表面张力计算分为流体粒子间的受力计算和流体粒子与界面之间的受力计算，流体粒子与界面之间的受力计算体现的是液体在固体表面的浸润性，计算模型为基于杨氏方程的界面浸润模型，因此在进行表面张力的镜像数值补偿计算时，因注意将流体粒子间的计算和与固体间的计算分开进行。下式(1)为流体粒子间的受力计算，式(2)为流体粒子与固体粒子间的受力计算，最后将两个力进行加和(如式(3))，这符合。

步骤5，进行压力泊松方程的建立和求解。

镜像粒子不参与压力泊松方程的构造，但参与其他流体粒子的粒子数密度计算。

步骤6，给镜像粒子按照镜像关系赋压力值。

镜像粒子的压力依然按照镜像关系从与之对应的流体粒子获得，这是为了在下一步根据压力梯度修正时对称边界范围内的流体粒子可以得到正确的压力修正结果。

步骤7，通过求解压力梯度对速度临时量进行修正。

在速度修正后同样需要对镜像粒子的空间和速度按照和流体粒子的镜像关系进行调整，以便于之后其他物理量的求解(比如温度的求解)。

步骤8，释放存储临时粒子的数组空间。

释放存储临时粒子的数组空间是为了保证在下一时间步依照更新后的流场获得全新的镜像粒子的总数。

本发明的原理如下：

本发明在计算具有对称特性的MPS三维算例时，可以按照对称面对流场进行截断，只布置对称面一面的流场，之后实际计算时也只计算该部分流体，起到节省计算资源、减小计算消耗的目的，具体方案如下：

确定计算对象是否具有流动对称性，判断算例是否适用对称边界条件。

当可用对称边界条件时，仅布置对称面一侧的流场为计算初场。

在显示部分的计算步骤中，为对称面影响域内的流体粒子根据“镜像法”进行物理量的数值补偿。

在构建压力泊松方程进行隐式求解时，对系数矩阵对角元的数值根据“镜像法”进行数值补偿。

通过压力梯度进行粒子速度和坐标的修正时，对对称面影响域内的流体粒子根据“镜像法”进行压力补偿，计算正确的压力梯度。

本发明适用于二维移动粒子半隐式法和三维移动粒子半隐式法。对于一切具有对称特点的算例均可使用，且当算例具有一组正交的对称面时，该方法依然可以使用。镜像粒子在每一时间步计算之前都会根据当前的流体粒子信息重新排布，每一个镜像粒子都有一个与之根据对称面对称的流体粒子，物理量也对称。以动态镜像粒子补充的方法保证靠近对称边界的粒子的物理量基于影响域内其他粒子加权平均的计算方法可以正确地进行，补偿计算依据下以下公式：

实施例：

如图1所示，本发明对称边界条件的算法，包括以下步骤：

初始布置需要按照对称边界条件建立，通常将对称边界条件设为x＝0、y＝0或 z＝0的平面，将被被截断为1/2或者1/4的流体域布置在坐标轴的正方向，而之后的镜像粒子将会被布置在对称轴的负方向，镜像粒子与对称面影响范围内的流体粒子是一一对应的。

为初步验证这种计算的正确性，计算如图2和图3所示的二分之一方形流域和完整方形流域临时粒子数密度的对比。临时粒子数密度n*是MPS法中判定表面粒子的关键依据，也是参与压力泊松方程求解、保证流体不可压缩的基本条件。粒子数密度n*计算公式如下：

图2以云图的形式展示了完整方形流域的右半部分和1/2方形流域的粒子数密度n*的计算结果比较。1/2方形流域在对称补偿模型的修正下得到的粒子数密度与完整方形流域相同。图3重点考察了两种流域中y＝0、y＝1两个轴上各个粒子的粒子数密度，通过线图定量地验证了对补偿模型的正确性。

初始条件的设置无需针对对称边界条件做特殊处理。

统计对称面影响域内流体粒子数的个数Nm和与流体粒子的应对关系。其中，对称面影响域是指对称面向流体所在的方向4.1l

按照所有粒子(包括本身参与物理量求解的流体粒子和第一层壁面粒子、幽灵粒子、镜像粒子)重新为存储计算粒子的数组

进入到MPS法对NS方程的求解中，首先进行的是速度临时量的求解，这里涉及到了体积力、粘性项和表面张力的计算。体力直接施加在具体的流体粒子上，但粘性项涉及速度扩散项的计算，需要考虑粒子与粒子之间的空间关系对物理量造成的影响。

MPS方法中流体粒子物理量的计算通过影响与内其他粒子物理量的加权平均获得。实际计算时只对初始布置的流体粒子进行求解计算，在计算时，一个具体的流体粒子i的核函数影响域范围内需要包括镜像粒子和壁面粒子。这里配合图1来解释：进行对称边界影响区域时内的流体粒子i的物理量计算时，除了计算粒子i的影响域内本身存在的流体粒子，还需要考虑粒子i的影响域在超过对称边界以外的部分本应包含的粒子。式(1)表示了对称补偿模型需要计算的每一项，i’表示镜像区域与当前求解的粒子i互为空间对称的镜像粒子。i与i’的速度矢量在对称边界法向方向的分量互为相反数，切向方向的分量相等，温度与压力等标量也相等。Ω

对于表面张力的计算，需要分别考虑流体粒子与流体粒子之间的受力关系和流体粒子与固体粒子之间的浸润关系。这里针对含有表面自由能模型的MPS法对称边界条件做以说明。

表面自由能模型计算流体粒子的表面张力时，是依据势能函数计算全场流体粒子的受力。当流体粒子处在流体内部时，该粒子影响域内的其他粒子是近似对称分布的，因此其受力平衡。而表面流体粒子会受到一个垂直于表面、指向流体内部的的力，即表面张力，其公式如式(4)：

式中，f

粒子最终所受表面张力项为表面张力和界面张力之和：

粘性项涉及拉普拉斯算子计算，计算公式如下：

式中：d表示维度，λ为修正量，n

由于流体粒子的临时空间坐标和临时速度矢量都发生了变化，因此镜像粒子的空间位置和速度也需要按照对应的流体粒子进行更新，以保证接下来的压力泊松方程的构造是正确的。

流体的速度场和压力场通过压力Poisson方程式得到耦合，将其左端项以扩散模型替代，可以离散得到一组线性代数方程组，求解方程组即可得到压力值P。对动量方程中压力项的计算是模拟流体不可压缩条件的重要内容，也是整个MPS 方法求解过程中的关键步骤。

将压力Poisson方程进行离散求解，将左端项进行离散变形后得到：

将式(6)写成矩阵Ap＝b的形式，系数矩阵为：

其中，N为参与压力计算的粒子数，这里只包括流体粒子和第一层壁面粒子，不包括镜像粒子。设a

a

未知数为矩阵p中的元素各粒子压力值P

压力梯度求解之后仍然要赋予镜像粒子相应的压力值以便于下一步压力梯度的计算。

压力梯度计算需要使用MPS法梯度算子，公式如下：

式中，P

所求解的压力梯度值将用于MPS法中N-S方程求解的最后一项——速度修正，然后，该时间步内的求解全部完成。以上求解过程涉及到MPS法原始算法的部分可参考MPS法的提出者Seiichi Koshizuka的论文：S.Koshizuka,Y.Oka. Moving Particle Semi-implicit Method for Fragment of Incompressible Fluid.Nuclear ScienceEngineering.1996.

当涉及到其他物理量的计算，如MPS法研究较多的温度场求解，同样按照镜像关系对镜像粒子按照镜像关系赋予物理值，以满足其在基于核函数的加权计算中可以得到真实准确的计算结果。

本发明所提出的对称边界条件一方面是为MPS法在边界条件上的研究进行补充，另外也为节省计算量、降低多余的计算消耗提供一种操作简单、效果明显的方法。图4分别利用对称边界条件模拟了1/2液滴、1/4液滴在平板上的浸润情况，和完整液滴的模拟过程进行对比，液滴的高度变化在图5中展示，三条变化去曲线重合度较好。

图6展示了1/4液滴在理论接触角为120度时的接触角验证模拟，并分别测量了X-Z平面、Y-Z平面上的接触角，分别为123.54度和123.36度，其平均值为123.45度，与理论界的误差为2.9％，具有可接受的计算精度。

图7展示了1/4液滴、1/2液滴和完整液滴完整计算过程中的CPU时间消耗情况，横轴代表模拟对象发展的时间，纵轴为相应时间的CPU时间，能够看出粒子数的减半对计算效率的提升是成指数型增长的。

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