CN110956005A - 一种基于粒子聚类和拉格朗日方法的高压消防水枪仿真方法 - Google Patents

Abstract

本发明涉及一种基于粒子聚类和拉格朗日方法的高压消防水枪仿真方法，以实现高效且逼真的高速流体仿真为目的，并使其能用于具体的消防仿真场景。本发明使用拉格朗日法求解N‑S方程进行流体仿真的核心思想，在SPH流体仿真基础上，围绕流体粒子的速度域聚类结果展开研究。首先利用k‑means聚类算法根据粒子的速度进行聚类，得到高速流体粒子与慢速流体粒子分组，然后利用双尺度时间步长框架分别对高速流体粒子和慢速流体粒子进行流体动力学计算，并通过设计耦合算法整合高速流体粒子与普通流体粒子，实现逼真、快速、兼顾细节的大规模流体仿真，并使其能够在消防场景中得到相应的运用。

Description

一种基于粒子聚类和拉格朗日方法的高压消防水枪仿真方法

技术领域

本发明涉及一种基于粒子聚类和拉格朗日方法的高压消防水枪仿真方法，属于基于计算机图形学的流体动画仿真领域。

背景技术

目前，流体仿真是一个跨学科、跨平台的研究方向及工程领域，但是受限于现有硬件水平的计算能力的限制和人们对物理世界的不完全了解所导致的仿真算法的局限性，对流体模型进行实时高效、逼真且稳定的仿真模拟还存在诸多困难。近年来，满足工程和科研要求的真实流体仿真备受国内外学者和科研机构的关注。基于物理的真实感流体仿真研究对于从事流体相关的科研工作者和计算机图形学的研究者来说无疑面临着巨大的挑战，也具有巨大的吸引力，其在工程实践领域也具有较大的应用价值。

基于物理的流体仿真方法主要分为三种：基于网格思想的欧拉方法和基于粒子思想的拉格朗日方法以及两类思想结合的混合方法。拉格朗日方法中的SPH方法因为其丰富的粒子细节表现能力及良好的整体场景表现能力而常被应用于流体仿真的工程实践当中。受限于当前计算能力的限制，对于流体细节的精细刻画、扩大仿真场景的规模的同时保证流体仿真的效率这三者具有不可调和的矛盾，在进行精细流体细节渲染的同时很难保证实现大规模流体仿真。

在消防演练相关的仿真场景中，流体仿真是其中的一项重要的支撑内容。消防场景下的流体仿真，与一般场景中的流体仿真既有诸多共性，也有很多的细节的区别，最大的不同体现在在消防场景中的消防设施所生成的水流一般是高压强的，高压强的流体一般在进入大气环境后伴随有较高的流速，即粒子的初始速度较高。当水流撞击在火源及地面时所产生的流体细节比较丰富，同时，高速的流体粒子也给仿真系统的稳定性和实时性要求提出了较高的要求。

为了解决上述问题，本发明提出一种基于粒子聚类和拉格朗日方法的高压消防水枪仿真方法，以实现兼顾速度、规模和细节的SPH流体仿真为目标，围绕基于物理的流体仿真方法展开研究，首先用k-means算法对流体进行聚类，然后在粒子聚类结果的基础上，对高速流体粒子进行细致的流体计算，同时减少普通流体粒子的费时计算过程，从而实现了兼顾细节的大规模交互场景的实时逼真绘制，使其能在消防仿真场景中得到应用。与现有的SPH 流体模型(如PCISPH、DFSPH)相比，本方法的方法由于在慢速流体计算过程中省略了对总体流体行为无影响的计算步骤，因此显著提高了计算效率。

发明内容

本发明解决的技术问题是：克服了现有的SPH流体仿真方法的效率问题，提供一种基于粒子聚类和拉格朗日方法的高压消防水枪仿真方法，解决了流体实时性与仿真规模及细节之间的矛盾，满足了对消防仿真过程中进行高效逼真模拟的需要。

本发明采用的技术方案为：一种基于粒子聚类和拉格朗日方法的高压消防水枪仿真方法，包括以下五个步骤：

步骤(1)、应用粒子聚类算法，根据粒子的速度域数值对粒子进行聚类，在粒子聚类算法收敛时构造两个聚类的分组结果，并统计两个聚类的分组结果中的速度分布，求得阈值，根据阈值将流体粒子分为高速流体粒子和普通流体粒子两个集合，大于或等于速度阈值的为高速流体粒子，低于阈值的作为普通流体粒子；

步骤(2)、结合聚类算法所得到的分类结果，根据Courant–Friedrichs–Lewy条件计算得到两种时间步长：主步骤时间步长和子步骤时间步长，每个主步骤包含若干个子步骤，普通流体粒子由于运动速度缓慢，粒子属性变化量小，每个主步骤时间步长更新一次；高速粒子运动剧烈，为保证流体细节不丢失，每个子步骤时间步长更新一次，即假设每个主步骤中包含n个子步骤，则高速粒子没更新n次时，普通流体粒子更新1次；

步骤(3)、通过SPH(Smoothed Particle Hydrodynamics)方法对粒子系统进行流体动力学计算：首先导入场景的静态边界模型，建立流体仿真所需的场景边界，使用粒子生成导入的边界模型的的表面结构，并使用流体的交互方式使固体对流体产生力的交互作用；然后确定初始粒子域的位置以及后续粒子的生成位置，求解Navier-Stocks(N-S)方程，通过离散的加权求和方式解算出粒子的各物理量；再以步骤(2)中的时间步长为积分量，对高速流体粒子和普通流体粒子分别进行更新，计算得到所有粒子在下一时刻的位置，更新粒子的速度和位置信息，生成新的粒子位置信息及结构，并更新相关的物理量；

步骤(4)、把高速流体粒子和普通流体粒子进行集成，利用耦合算法实现双时间步长模型整合，计算并更新流体粒子的物理量，物理量包括密度、速度和位置；根据耦合算法实现流体仿真过程中高速流体粒子和普通流体粒子在同一仿真场景中的正确耦合；

步骤(5)、为实现真实感仿真效果，将仿真场景中的固体物体使用片面结构进行渲染呈现，而渲染的流体粒子细节直接用OpenGL绘制成粒子球在场景中进行表示，或使用Marching Cube算法生成片面加构，加以渲染呈现。

步骤(1)中，所述的粒子聚类算法采用k-means聚类算法，根据粒子速度值将粒子聚类成两个分组，该过程使用流体的中间速度值作为聚类算法的输入内容，并设定初始的平均值使其更快地实现算法的收敛，根据k-means聚类算法的原理，收敛的判断条件表示为：

上式中，S_i表示所得到的第i个聚类的粒子集，v_j为第j个粒子x_j的速度值，

为第i 个粒子集内粒子的速度平均值。

步骤(2)中，所述主步骤时间步长和子步骤时间步长的计算，对于子步骤的时间步长大小的推导，使用所有粒子的最大速度值代入CFL收敛公式进行求解，而低速流体粒子组的最大速度则被用于代入到CFL公式中求解出主步骤的时间步长大小；所有流体粒子的最大速度用v_max，all表示，则子步骤时间步长Δt_mini表示为：

低速流体集合中的速度最大值为v_max，low，主步骤时间步长Δt_major大小表示为：

上述两式中，参数a∈[0，0.5]，d为粒子半径，同时设定阈值N_tres，若主步骤时间步长大小和子步骤时间步长大小的比例低于这一阈值，一个主步骤将会只包含一个标准的SPH 模拟步骤；另外，为防止这一比例过大而导致模拟效果失真，设置一个最大的子步骤数量的上限N_max，一个主步骤内部所含有的子步骤的数量关系如下：

为保证算法稳定，设定N_tres的值不超过1.5。

本发明的原理在于：

(1)通过k-means对流体进行聚类分析，得到高速流体粒子组和普通流体粒子组，从而为后续的步骤中有针对性地对分属不同粒子组的粒子采用不同的仿真重力、结合双尺度时间步长框架提供了数据支持。

(2)为了准确地渲染流体仿真的细节表现，使用SPH方法对流体粒子进行流体动力学计算，包括在重力、空气、固体边界共同作用下的外力以及流体粒子内部相互作用下的内力的计算，通过进行隐式的循环迭代的压强项求解计算，将粒子运动计算过程中所产生的密度误差减小到阈值满园内，从而保证了流体粒子在密度上保持了不可压缩的特性，达到了仿真效果的稳定性和真实性。

(3)为了对流体进行逼真的细节刻画，本发明提出了一种基于粒子聚类和拉格朗日方法的流体仿真技术。该方法在SPH流体仿真的基础上，使用双尺度时间步长框架，保证了流体仿真细节同时减少了部分费时流体计算，从而达到了良好的流体真实视觉表现和高效的流体仿真过程。

本发明与现有技术相比的优点在于：

(1)本发明提出的基于粒子速度域的聚类过程，易于在工程应用中实现并高效运行，同时为后续的流体仿真步骤提供了重要的预处理数据。

(2)高精度的拉格朗日方法，本发明利用SPH对流体进行流体动力学计算，能够呈现丰富的细节表现，特别是在高速流体仿真中能够呈现出足够的流体细节。

(3)结合基于粒子速度域的粒子聚类算法以及配套的双尺度时间步长框架，既发挥了 SPH方法逼真的渲染效果的优势，又发挥了对不同速度大小的粒子进行不同处理所带来的效率优势。

附图说明

图1为本发明的一种基于粒子聚类和拉格朗日方法的高压消防水枪仿真方法流程图；

图2为本发明中基于双尺度时间步长框架的SPH改进算法示意图；

图3为本发明中高速流体冲击普通流体效果渲染图，其中流体由不同速度集合的粒子球渲染；

图4本发明中高速流体冲击普通流体效果渲染图，其中流体由网格模型进行渲染。

具体实施方式

图1给出了基于粒子聚类和拉格朗日方法的高压消防水枪仿真方法的总体处理流程，下面结合其他附图及具体实施方式进一步说明本发明。

本发明提供一种基于粒子聚类和拉格朗日方法的高压消防水枪仿真方法，步骤如下：

步骤(1)、基于粒子速度域的流体粒子聚类：应用k-means聚类算法，根据粒子的速度域数值对粒子进行聚类，这是流体仿真的预处理部分。

步骤(2)、流体自由表面追踪：应用PLSM方法构建流体自由表面。输出步骤(1)中得到的流体网格，用PLSM方法标记流体运动界面在LBM网格内部的位置，追踪流体界面，以此作为流体的近似表面。

步骤(3)、流体细节仿真：通过SPH方法来对流体细节进行物理仿真，模拟流体运动过程中产生的诸如浪花、水珠等细节。

步骤(4)、集成仿真场景、进行算法耦合；把步骤(1)、(2)、(3)中粒子聚类过程、双尺度时间步长结构和SPH仿真的流体集成到同一个场景中，利用耦合算法进行流体动力学计算过程，控制流体仿真过程中高速流体粒子和普通流体粒子的全局耦合计算时机。

具体实施为SPH方法和k-means聚类方法的结合，主要步骤介绍如下：

1、基于粒子速度域的粒子聚类方法

本发明采用k-means方法对消防仿真场景中的流体粒子进行聚类计算，k-means是一个常用于数据分析的算法，其算法目的是将一组样本根据其的大小或位置，聚类到若干个互不相交的分组当中，并使得每个分组中的样本值离该分组中的平均值相比其它分组更小。由于 k-means算法具有易于实现，且能够将样本值快速地划分为指定数量的分组，因此在各领域被广泛应用。算法的收敛条件是达到组内平方和的最小值，将聚类算法与流体仿真相结合，首先需要提供进行聚类的数据，具体的数据为粒子在外力作用下所得到的中间速度值，其形式如下所示：

其中，v_t为上一时间节点上的粒子速度值，v_t+Δt是当前时刻的粒子速度值，f为粒子在当前状态下所受的外力大小，m为粒子的质量，Δt为仿真过程中的时间步长大小。

由于k-means算法本身具有一定的缺陷，在仿真的过程中，其聚类的结果与之前一段时间内所得到的结果缺乏时间和空间上的联系，从而会导致处于中间速度值的粒子会在两个聚类中间反复摇摆，从而产生了效果上的抖动。为了防止这一抖动所造成的仿真效率上稍许降低，本发明额外引入了滑动窗口算法对所得到的聚结果进行平滑处理，从而增强聚类结果在时空上的相互联系。

具体的滑动窗口算法是：以仿真过程中的一段给定长度的时间为单位，统计这一段时间内各个主步骤帧内的粒子聚类结果及其比例，判断其在一段时间内的趋势。如果趋势持续地朝向其中一侧变化，或者在一定时间内，上下浮动的比例不超过给定的阈值，则不对粒子聚类的结果进行整流优化；如果在一定的时间内，粒子聚类结果的比例反复上下波动比例高于给定的阈值，则根据其波动的方向，将波动的比例强制缩小到阈值的范围之内，使其波动的幅度变小，从而减少仿真场景中受到波动的粒子数。

2、基于SPH方法的流体仿真方法

基于SPH的流体仿真方法主要包含以下四个子过程：邻域搜索过程、外力项计算过程、压强项计算过程、速度及位置更新过程。

因为SPH是一个以一定范围内的粒子相互作用产生力的关系为基础的流体仿真算法，一个粒子周围的粒子情况会在计算中多次被使用，因此邻域搜索过程是为了寻找并记录仿真场景中的流体粒子以供后续使用，以减少后续的计算工作。邻域搜索的范围是有一定限制的，通常设定为粒子半径的若干倍。在算法的实现过程，需要将粒子根据其位置划分到其所对应的网格中，而网格的半径则设定为SPH方法中内核函数的支持半径的大小。在搜索的过程中，只针对粒子所在的网格本身及与直接相邻的周边26个网格进行搜索即可，从而减少了搜索过程中所需要的遍历的粒子数量，加快了整个搜索过程所需要的时间长度。

外力项的计算过程包含流体的粘性项、空气的作用项等计算。对于消防场景中的高压水枪所对应的流体类型，属于弱粘性的流体类型，因此在仿真的时候直接使用XSPH中的粘性计算方法对流体之间所表现的粘性作用；而对于空气对于水流的作用，在高压消防水枪的场景中，空气对于水流的作用较弱，因此使用了人工的空气拉力作用来计算空气对水流所产生的作用力，以减少在仿真计算过程中的计算过程，使仿真的过程速度更快，从而达到实时性所需的要求。

压强项计算过程是流体仿真过程中最为费时的一个项目，即使是高压水枪，水流内部存在高压，且在从水管中输出时水的流速较大，但流体内部依然保持着固有的密度不变，流体在仿真的过程中依然需要满足不可压缩性质的要求，但进行解算的时候会带来一些困难，因为速度的增大减少了仿真过程中所使用的时间步长的大小，从而增加了仿真相同时间迭代收敛所需要的循环次数，进而出于满足实时性要求的需求，限制了仿真时所能计算的流体粒子数量的最大值。为了加快流体动力学的计算过程，在这里应用了目前国内外所研究和应用的 SPH方法中最先进且效率最高的DFSPH作为压强项求解的基本方法。DFSPH分为两个部分对压强项进行计算，即粒子速度的无旋度条件的满足和流体密度固定不变条件的满足，从而加快了流体不可压缩性质的满足过程。

经过外力的计算和压强项的计算之后，可以综合外力和内力的矢量，结合粒子本身的质量大小，计算得到作用于粒子上的加速度。同时结合时间步长大小，以及粒子在上一时间步骤结束时的速度值，可以得到当前步骤结束时的速度值，同时将速度值所产生的位移更新到粒子的位置变化上。

3、双尺度时间步长框架和SPH方法的耦合

在得到聚类结果之后，需要针对所属不同分组的流体粒子使用不同的仿真和调度策略，本发明采用了一种双尺度的时间步长框架，使属于高速流体粒子组的流体粒子得到充足的流体动力学计算机会，从而使消防场景中的高压水枪喷射出的流体得到足够的计算机会，能够产生足够真实的流体细节，使得场景渲染的效果更为真实；而属于普通流体粒子组的流体则省略了一部分费时的计算过程，如压强项计算的过程，使得整体的流体仿真效率得到提升。同时为了保证流体整体行为的正确性，会在一定仿真步骤进行之后，对普通流体的时间步长大小进行弥补，从而弥补普通流体粒子在此之前所缺失的流体动力学计算，两者的耦合计算使得场景中的粒子成为一个紧密耦合的整体。这两者的耦合过程是通过子步骤和主步骤来进行的。

首先需要计算出双尺度时间步长框架中的主步骤和子步骤的时间步长大小，根据CFL 条件的约束，

其中参数a∈[0，0.5]，d为粒子半径。将上式分别代入高速流体粒子组的最大速度值以及普通流体粒子组的最大速度值，得到主步骤的时间步长大小和子步骤的时间步长大小Δt_mini以及Δt_major，并根据其求出比值向下取整，得到主步骤中的子步骤数量。同时，当比值较小时，即仿真场景中粒子的速度较为接近时，将此数量设置为1，而当比值较大时，则设置为一个固定的最大值，以防导致流体仿真所可能导致的失真现象，具体表达形式如下所示：

每个主步骤中含有N_step数量的子步骤，当N_step＝1时，则所对应的主步骤中的子步骤与常规的SPH步骤一致。如果含有多个子步骤，则除了第一个子步骤中含有粒子聚类的预处理步骤以及最后一个子步骤中含有对全体流体粒子进行流体动力学计算的过程之外，其余的子步骤中均只对场景中的高速流体粒子进行流体动力学计算。而在最后一个子步骤中，为了强化流体的细节效果，针对普通流体粒子在之前的子步骤静止等待高速流体粒子运动的情况，对全局统一的时间步长的情况进行了相应的改变，具体表现为增加了普通流体粒子在最后一个子步骤中的时间步长大小：

其中μ为由人工所指定的参数，一般选取范围不应当使普通流体粒子违反CFL条件的约束。具体主步骤与子步骤之间的关系如图2所示。

算法表示的步骤为：

1)计算得到全局粒子的中间速度值，并进行暂存；

2)通过k-means聚类算法对全体流体粒子进行聚类分析，结果保存为两个聚类，其中一个聚类中的粒子速度值比另外一个聚类中的粒子速度值大；

3)根据所得到的聚类结果，分别得到两个聚类的最大速度值，根据最大速度值以及CFL条件推导到对应的两个时间步长大小，其中大的时间步长大小设定为主步骤时间步长大小，较小的时间步长大小设定为子步骤时间步长大小；

4)根据两个时间步长大小计算其比值，得到主步骤所包含的子步骤的数量，以及各子步骤的结构成分。当比值较小时，则一个主步骤只包含一个子步骤，当比值较大时，则设定一个最大值，一个主步骤最多包含该最大值数量的子步骤；

5)根据主步骤中所包含的子步骤的数量，使用SPH方法计算仿真场景中粒子的流体动力学算式；

6)在除了主步骤的最后一个子步骤之外的子步骤中，只针对高速流体进行流体动力学计算，更新其位置以及速度值，普通流体粒子静止等待高速流体粒子的运动过程，直接结束；

7)在主步骤的最后一个子步骤中，针对仿真场景中的所有流体粒子进行仿真，但为了使仿真所得到的结果更为精确，在这里对普通的流体粒子进行时间步长大小方面的补偿，使得普通流体粒子能够更好地与高速流体粒子进行耦合和碰撞互动，使之成为一个有机的整体；

8)当主步骤完成最后一个子步骤后，一个流体仿真的阶段完成，将进入下一个仿真步骤的模拟计算过程。

4、渲染绘制

最后，对于消防仿真场景中的固定物体，直接根据其面片结构进行相应的渲染及呈现，而场景中的流体粒子则用OpenGL绘制成粒子球，或者使用Marching Cube算法将场景中的粒子球进行采样，生成流体表面的三角面片结构，从而使仿真效果更具真实感。如图3中所示，图中龙模型为固定物体，则根据龙模型的三角面片结构进行渲染，冲击固定物体的流体粒子则绘制为粒子球，其中深色粒子球表示普通流体粒子，浅色粒子球表示高速流体粒子。如图4中所示，图中龙模型为固定物体，则根据龙模型的三角面片结构进行渲染，冲击固定物体的流体粒子使用Marching Cube算法进行采样，生成流体表面的三角面片结构，渲染为连续的流体表面模型。

本发明未详细阐述的技术内容属于本领域技术人员的公知技术。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (3)

1.一种基于粒子聚类和拉格朗日方法的高压消防水枪仿真方法，其特征在于，包括以下五个步骤：

步骤(1)、应用粒子聚类算法，根据粒子的速度域数值对粒子进行聚类，在粒子聚类算法收敛时构造两个聚类的分组结果，并统计两个聚类的分组结果中的速度分布，求得阈值，根据阈值将流体粒子分为高速流体粒子和普通流体粒子两个集合，大于或等于速度阈值的为高速流体粒子，低于阈值的作为普通流体粒子；

步骤(2)、结合聚类算法所得到的分类结果，根据Courant–Friedrichs–Lewy条件计算得到两种时间步长：主步骤时间步长和子步骤时间步长，每个主步骤包含若干个子步骤，普通流体粒子由于运动速度缓慢，粒子属性变化量小，每个主步骤时间步长更新一次；高速粒子运动剧烈，为保证流体细节不丢失，每个子步骤时间步长更新一次，即假设每个主步骤中包含n个子步骤，则高速粒子没更新n次时，普通流体粒子更新1次；

步骤(3)、通过SPH(Smoothed Particle Hydrodynamics)方法对粒子系统进行流体动力学计算：首先导入场景的静态边界模型，建立流体仿真所需的场景边界，使用粒子生成导入的边界模型的的表面结构，并使用流体的交互方式使固体对流体产生力的交互作用；然后确定初始粒子域的位置以及后续粒子的生成位置，求解Navier-Stocks(N-S)方程，通过离散的加权求和方式解算出粒子的各物理量；再以步骤(2)中的时间步长为积分量，对高速流体粒子和普通流体粒子分别进行更新，计算得到所有粒子在下一时刻的位置，更新粒子的速度和位置信息，生成新的粒子位置信息及结构，并更新相关的物理量；

步骤(4)、把高速流体粒子和普通流体粒子进行集成，利用耦合算法实现双时间步长模型整合，计算并更新流体粒子的物理量，物理量包括密度、速度和位置；根据耦合算法实现流体仿真过程中高速流体粒子和普通流体粒子在同一仿真场景中的正确耦合；

步骤(5)、为实现真实感仿真效果，将仿真场景中的固体物体使用片面结构进行渲染呈现，而渲染的流体粒子细节直接用OpenGL绘制成粒子球在场景中进行表示，或使用Marching Cube算法生成片面加构，加以渲染呈现。

2.根据权利要求1所述的一种基于粒子聚类和拉格朗日方法的高压消防水枪仿真方法，其特征在于：步骤(1)中，所述的粒子聚类算法采用k-means聚类算法，根据粒子速度值将粒子聚类成两个分组，该过程使用流体的中间速度值作为聚类算法的输入内容，并设定初始的平均值使其更快地实现算法的收敛，根据k-means聚类算法的原理，收敛的判断条件表示为：

上式中，S_i表示所得到的第i个聚类的粒子集，v_j为第j个粒子x_j的速度值，

为第i个粒子集内粒子的速度平均值。

3.根据权利要求1所述的一种基于粒子聚类和拉格朗日方法的高压消防水枪仿真方法，其特征在于：步骤(2)中，所述主步骤时间步长和子步骤时间步长的计算，对于子步骤的时间步长大小的推导，使用所有粒子的最大速度值代入CFL收敛公式进行求解，而低速流体粒子组的最大速度则被用于代入到CFL公式中求解出主步骤的时间步长大小；所有流体粒子的最大速度用v_max,all表示，则子步骤时间步长Δt_mini表示为：

低速流体集合中的速度最大值为v_max,low，主步骤时间步长大Δt_major小表示为：

上述两式中，参数a∈[0,0.5]，d为粒子半径，同时设定阈值N_tres，若主步骤时间步长大小和子步骤时间步长大小的比例低于这一阈值，一个主步骤将会只包含一个标准的SPH模拟步骤；另外，为防止这一比例过大而导致模拟效果失真，设置一个最大的子步骤数量的上限N_max，一个主步骤内部所含有的子步骤的数量关系如下：

为保证算法稳定，设定N_tres的值不超过1.5。

Images