CN106650064B - 一种基于粒子模型的凝结现象仿真方法 - Google Patents

Abstract

本发明公开了一种基于粒子模型的凝结现象仿真方法，该方法包括:1)对气体在固体表面凝结的现象离散建模，考虑露点温度以及湿度扩散模型中相对湿度属性对于凝结现象细节的影响；2）提出固液耦合的边界密度校正算法，弥补了边界处光滑粒子流体动力学模型粒子缺失的问题，同时在不降低时间步长的情况下防止粒子渗透；3）使用粒子系统对空气进行离散建模,可更直观地表现仿真过程，易于调整和设计。应用本发明能够便捷地调整空气的流动、热传导属性等，解决了现有空气热交换现象仿真效果不直观和细节不足的问题，具有一定的实用价值。

Description

一种基于粒子模型的凝结现象仿真方法

技术领域

本发明属于物理仿真领域，具体地说是一种基于粒子模型的凝结现象的仿真方法，其涉及到流体动力学、GPU并行、空气扩散与热交换理论以及计算机图形学等。

背景技术

空气的热交换现象在虚拟现实、医学、制造行业等领域都处在十分重要的地位。仿真学中经典的建模方法有网格法和粒子法两种。在三维中，网格法将物体分成块状，根据块与块之间的力交互以及自身的形变物理规律来离散化，对较小形变以及细节不敏感的场景仿真有较大的优势。然而网格法难以描述局部细节，无法精确计算出用来描述细节的物理量。与之相对，粒子法是对材质进行采样得到一个点集合，直接在点集合上计算力的交互，直观而且容易理解，对于细节上的表述十分清晰。在凝结现象中，网格化的气体无法具体计算出凝结液滴位置点的坐标，只能指出大概的凝结区域，在精度上不足以描述凝结现象。

光滑粒子流体动力学(SPH)是一种流体模拟算法，基本思想是将流体场离散成粒子系统进行计算，在计算单个粒子的属性时，根据属性连续的假设，对邻居粒子的同一属性根据距离进行插值计算。插值计算中的权重函数称之为核函数，常见的核函数有Poly6函数和Spiky函数。SPH方法能够直观地把物理规律施加到粒子系统，还可以用在游戏的实时的交互软件中。

基于物理的固液耦合边界处理中经典方法有：Monaghan等人提出的虚拟力法、Libersky等人提出镜像粒子法、Schechter等人提出的Ghost粒子法等等。这些方法旨在防止不同物体粒子之间的渗透问题，生成的仿真现象不够逼真。

发明内容

本发明的目的在于克服现有空气热交换仿真不直观，不注重细节的问题，提供了一种基于粒子模型的凝结现象的仿真方法，该方法使用光滑粒子流体动力学对凝结现象进行仿真，保留了细节，构建了气体在固体表面凝结现象的离散模型。

本发明的目的是这样实现的：

1、一种基于粒子模型的凝结现象的仿真方法，其特征在于该方法包括以下步骤：

a)基于光滑粒子流体动力学(SPH)模型模拟流体的运动，具体包括：

i)区分边界粒子

对所有粒子增加额外的属性，对气体粒子赋值为0，其他粒子赋值为1，构成标量场；对每个粒子计算该标量场的梯度，梯度不为零的固体粒子归到边界粒子；

ii)各粒子的密度计算

每个粒子的密度通过周围相同类型粒子的密度插值得到；在插值计算边界上气体粒子的密度时，引入密度校正算法：在插值计算过程中，把边界粒子的密度加入到插值计算中，其中边界粒子密度根据采样的疏密性，设定成气体密度的常数倍；

iii)各粒子的加速度计算

根据粒子密度，通过运动控制Navier-Stokes方程计算粒子的加速度：

其中u为流体速度，t为时间，ρ为粒子密度，p为压强，μ为黏性系数，是一阶微分算子，是二阶微分算子，f^ext为外力；方程等号右边依次为压力项，黏性力项和外力项；

在外力项中加入大气压带来的液滴阻力；计算液滴覆盖在固体表面的面积，配合外部大气压强和固体的摩擦系数，计算出由于液滴内外气压差而产生的液滴阻力：

f^friction＝μ_f×(p·s)

其中f^friction为液滴阻力，μ_f为摩擦系数，p为压强，s为液滴填充的面积；

由此得到的控制方程，用来计算粒子的加速度：

b)基于SPH模型模拟凝结现象，具体包括：

i)各粒子的湿度计算

在SPH模型上实现菲克定律，通过计算扩散系数和粒子之间的相对湿度差得到每一帧湿度的变化值：

其中RH代表相对湿度，t代表时间，c_diff代表扩散系数，N_j是粒子j的邻居粒子集合，m是质量，是二阶微分算子，r是粒子的位置，W是光滑核函数即插值权重函数，h是光滑核函数W影响区域的半径；

把计算得到的湿度变化值加上该粒子湿度，得到该粒子处新的湿度属性；

ii)各粒子的温度计算

粒子的温度变化量取决于邻居粒子的温度差：

其中T代表温度，t代表时间，c_d代表热传导系数，N_j是粒子j的邻居粒子集合，m是质量，是二阶微分算子，r是粒子的位置，W是光滑核函数即插值权重函数，h是光滑核函数W影响区域的半径；当相邻粒子材质不同时，为了保持边界处属性的连续性，热传导系数定为两种材质热传导系数的算术平均值；

把计算得到的温度变化值加上该粒子温度，得到该粒子处新的温度属性；

iii)各粒子的凝结判断

使用August-Roche-Magnus方程对Clausius-Clapeyron方程近似：

式中T代表温度，T_dp为露点温度，RH代表相对湿度，常量系数：b＝17.67，c＝243.5℃；计算出对于每个气体粒子的露点温度；当气体粒子温度低于该粒子处的露点温度时，判断该气体粒子发生凝结现象，完成粒子模型的凝结现象的仿真。

本发明在完成一帧凝结现象的仿真后，判断气体粒子距离固体-气体交界面的距离，如果该距离远大于固体边界粒子的核函数半径h，那么这种粒子只进行速度和位移的计算；当粒子远离固体边界粒子的时间超过时间阈值，删除该粒子。这一筛选可为下一帧运算加速。

本发明的有益效果：

本发明使用光滑粒子流体动力学对凝结现象进行仿真，保留了细节，构建了气体在固体表面凝结现象的离散模型。在进行凝结模拟时还考虑露点温度以及湿度扩散模型中相对湿度属性对于凝结现象的影响，通过湿度来控制空气粒子凝结出的液滴尺寸，同时可以有效避免局部粒子反复凝结造成的非真实性和粒子数量过多的问题。对于渗透问题，本发明提出气-固耦合的边界密度校正算法，弥补了边界处SPH参与计算粒子缺失的问题，在不降低时间步长的情况下防止粒子渗透。

总之，应用本发明可以精细、快速地对气-固耦合场景进行仿真。在模型的直观性、计算的高效性和仿真结果的真实感方面，都有显著提高。

附图说明

图1为气体-固体交互区间示意图；

图2为三种边界处理方法比较的示意图；

图3为本发明液滴凝结细节的效果图；

图4为本发明真实凝结场景仿真效果图。

具体实施方式

下面结合附图及实施例对本发明详细描述。

本发明的具体内容：

1)构建了气体在固体表面凝结现象的离散模型：考虑露点温度以及湿度扩散模型中相对湿度属性对于凝结现象的影响。

2)提出气-固耦合的边界密度校正算法：弥补了边界处SPH参与计算粒子缺失的问题，在不降低时间步长的情况下防止粒子渗透。

3)引入空气粒子建模：对气-固交互界面的仿真更加精细，同时易于调整和设计。

SPH的基本框架：

光滑粒子流体力学(smoothed-particle hydrodynamics，简称SPH)将流体场离散成粒子系统进行计算；计算单个粒子的属性时，根据属性连续的假设，使用邻居粒子的同一属性根据距离进行插值计算：

其中A_j是粒子j的属性值，N_j是粒子j的邻居粒子集合，m是质量，r_j是粒子j的位置，W是光滑核函数即插值权重函数，h是光滑核函数W影响区域的半径。

粒子的运动控制方程是经典的动量守恒方程Navier-Stokes方程：

其中u为流体速度，t为时间，ρ为粒子密度，p为压强，μ为黏性系数，是一阶微分算子，是二阶微分算子，f^ext为外力；方程等号右边依次为压力项，黏性力项和外力项；方程等号右边依次为压力项，黏性力项和外力项。在每一帧对每个点计算以上三部分力，汇总得到该点的加速度，更新点的位置。

本发明所述的考虑热传导的SPH流体模型，其具体为：

(1)粒子的温度变化量取决于邻居粒子的温度差，而改变速率则由热传导系数决定。热传导过程可通过SPH离散计算：

其中T代表温度，c_d代表热传导系数。当相邻粒子材质不同时，为了保持边界处属性的连续性，热传导系数定为两种材质热传导系数的算术平均值。

湿度扩散模型根据菲克定律，其中扩散系数定义为：沿扩散方向，在单位时间每单位浓度梯度的条件下，垂直通过单位面积所扩散某物质的质量或摩尔数。扩散系数的大小主要取决于扩散物质和扩散介质的种类及其温度和压力。空气相对湿度的扩散模型：

其中RH代表相对湿度，c_diff代表扩散系数。

当空气发生凝结现象，凝结所得到的水滴体积取决于空气的相关湿度。在完成凝结后，空气的相对湿度降低，通过湿度扩散模型平衡周围空气的湿度。

(2)发生热交换的区域主要是在具有较大温度差的不同材质的交界处，当空气粒子到固体边界的距离远大于边界固体粒子的核函数半径，那么该空气粒子对于交界面处固体粒子属性的贡献很小，所以迭代计算该粒子的属性价值很低，在本发明中称为低效用粒子。

在交界面处的粒子，离开交界面作用域后，根据速度和受力判断出其持续远离的趋势，那么这种粒子也判定为低效用粒子。对于低效用的粒子只考虑速度、位移，不考虑邻居粒子之间的热交换，力交互以及相变。当粒子跑出场景，或者变成低效用粒子的时间超过时间阈值后，进行删除操作。

本发明所述的凝结现象的仿真，其具体为：

(1)露点温度指空气在水汽含量和气压都不改变的条件下，冷却到饱和时的温度。在离散系统中相变过程的压强和温度的关系可以通过Clausius-Clapeyron方程来表示：

其中L是单位质量的物质所含有的潜伏热能，v是单位质量的物质所占有的体积。

August-Roche-Magnus对公式(5)给出了一个近似估计：

式中有几个常量系数，根据用途取不同的值，本发明中常量系数取：b＝17.67，c＝243.5℃出自David Bolton的论文。根据公式(6)可以得到在环境温度T时，相对湿度为RH的空气粒子对应的露点温度T_dp。

(2)在气体-固体的交界面，由于SPH计算气体粒子密度时，核半径内的固体粒子对密度计算没有贡献，同时占据气体粒子的空间，导致气体-固体交界面附近的气体粒子密度小于空气粒子静止密度，且由于密度差生成指向固体内部的压力，会导致渗透现象。另外较大的时间步长也会导致气体粒子渗入固体内部。一旦气体粒子进入固体内部，最直接的影响就是气体-固体的热交换发生在固体内部，而不是固体表面，出现诸如内部先凝结出液滴这种不符合物理规律的相变情况。

对于防止粒子渗透的方法，基于物理的经典方法有虚拟力法以及Ghost粒子法。虚拟力法是通过对粒子小半径范围内的其他粒子，施加与相对距离相关较大的力，迫使其离开受力范围。这种处理方法有明显的人为痕迹，以气-固为例，气体和固体之间会形成明显的真空层，大大降低了热传导的能力。Ghost粒子法的主要思想是在固体表面和内部随机采样更密集的粒子并通过松弛算法使粒子分布均匀化，采样粒子用来补足交界处SPH计算的有效粒子数，从而消去由于交界处密度差而产生的压力。但是Ghost粒子法算法十分耗时，需要额外的采样和松弛步骤，有很高的复杂度。

在时间复杂度和基于物理的要求下，本发明提出了一种基于SPH边界密度的校正算法。以气体-固体耦合为例：构造标量场，对于固体粒子赋值为1，其余赋值为0，类似2D下的高度场。对每个固体粒子计算标量值c_s：

在这之后对c_s标量场求梯度，那么对于固体内部，c_s梯度都是零，而对于边界粒子则梯度非零。同时梯度是一个从固体内部指向外部，沿表面法向的矢量。设定该矢量模的阈值c_grad，如果某粒子处c_s梯度值大于c_grad那么判定这个固体粒子是表面粒子。在使用SPH计算空气密度的时候，范围内如果存在这样的边界固体粒子，在计算该空气粒子密度时，把边界固体粒子的密度设为空气密度的常数倍，累加到密度核函数。由于固体粒子原本的分布要比Ghost粒子稀疏，为了能达到相同的贡献，需要增加一个常数项作为倍数。而当固体粒子局部数量足够多时，只需要把固体粒子密度设为距离最近的空气粒子密度即可。这一算法完全内嵌在现有的SPH里，十分易于实现和移植，而且时间复杂度等于SPH算法的复杂度，在GPU并行的情况下几乎没有额外的时间复杂度。

(3)在凝结过程中，空气粒子通过与较冷固体表面进行热交换，凝结成液态水，并依附在固体表面。附着的液滴很好地覆盖了固体表面，造成局部真空，产生内外气压差。而当液滴体积较小时，液滴受到由大气压作为正压力产生的摩擦力，平衡了重力，使液滴静止在固体表面。然而现有的SPH中，粘滞力无法提供这个摩擦力，粘滞力公式为：

力的大小取决于离散系统相邻粒子的速度差，当液滴静止时力为零，显然无法达到液滴静止在固体表面的现象。

本发明基于上述原理引入液滴阻力：

f^friction＝μ_f×(p·s) (9)

其中f^friction为液滴阻力，μ_f为摩擦系数，p为压强，s为液滴填充的面积。

实施例

图2展示了三种不同的边界处理方法，依次为虚拟力方法，Ghost粒子法和本发明方法。图中中间圆形粒子代表边界处理方法所涉及到的固体粒子，包围在外部的其余粒子为空气粒子。空气粒子的运动趋势是从左向右。

左边图片表示的虚拟力方法对于近距离的粒子施加较大的力，在图片左半部可以看到空气和固体之间有一条真空带，而在图片的右半部分可以看到空气粒子有明显的远离固体粒子的趋势。中间的Ghost粒子法涉及了大量的粒子，能很好处理边界粒子，然而需要频繁对固体的表面和体积进行采样和松弛操作，时间复杂度高且并行难度大。右边图片展示的本发明方法，仅需要修正空气粒子密度计算时固体表面粒子的贡献，就能够达到防渗透同时不推开空气的效果，由于其容易嵌入现有SPH框架，可以并行加速。实验中的Ghost粒子最大尝试数为30。由于算法中每采样一个新点都需要比较邻居区域中其他采样点的距离，并行的难度较大，选择在CPU上运行。需要注意的是Ghost粒子方法需要对物体进行体采样，所以物体的体积大小对于时间复杂度影响较大。

图3展示的是本发明空气在玻璃表面的凝结现象。左图中空气粒子在固体表面液化凝结，受到表面张力和大气压的影响，停留或者沿着固体表面滑落。另外液滴生成有局部密集的现象，这是因为湿度扩散模型则起到湿度扩散传播的效果，使得凝结过液滴的空气粒子可以受到周围空气粒子的湿度传播，在适当的条件下再次凝结。中间图片显示了当相邻水滴靠近时，两个水滴会进行合并成为一个更大的水滴。右图中固体表面生成了大小不一的液滴粒子，这是因为空气粒子凝结时相对湿度的差异性会直接影响凝结水滴的体积。

图4展示的是本发明对桌面上的玻璃制品施加湿度较大的空气，对真实场景的凝结现象进行仿真。

以上列举的仅是本发明的具体实施例。显然，本发明不限于以上实施例，还可以有许多变形。本领域的普通技术人员能从本发明公开的内容直接导出或联想到的所有变形，均应认为是本发明的保护范围。

Claims (1)

1.一种基于粒子模型的凝结现象的仿真方法，其特征在于，该方法包括以下步骤：

a)基于光滑粒子流体动力学SPH模型模拟流体的运动，具体包括：

ⅰ)区分边界粒子

对所有粒子增加额外的属性，对气体粒子赋值为0，其他粒子赋值为1，构成标量场；对每个粒子计算该标量场的梯度，梯度不为零的固体粒子归到边界粒子；

ⅱ)各粒子的密度计算

每个粒子的密度通过周围相同类型粒子的密度插值得到；在插值计算边界上气体粒子的密度时，引入密度校正算法：在插值计算过程中，把边界粒子的密度加入到插值计算中，其中边界粒子密度根据采样的疏密性，设定成气体密度的常数倍；

ⅲ)各粒子的加速度计算

根据粒子密度，通过运动控制Navier-Stokes方程计算粒子的加速度：

其中u为流体速度，t为时间，ρ为粒子密度，p为压强，μ为黏性系数，是一阶微分算子,是二阶微分算子,f^ext为外力；方程等号右边依次为压力项，黏性力项和外力项；

在外力项中加入大气压带来的液滴阻力；计算液滴覆盖在固体表面的面积，配合外部大气压强和固体的摩擦系数，计算出由于液滴内外气压差而产生的液滴阻力：

f^friction＝μ_f×(p·s)

其中f^friction为液滴阻力，μ_f为摩擦系数，p为压强，s为液滴填充的面积；

由此得到的控制方程，用来计算粒子的加速度：

b)基于SPH模型模拟凝结现象，具体包括：

ⅰ)各粒子的湿度计算

在SPH模型上实现菲克定律，通过计算扩散系数和粒子之间的相对湿度差得到每一帧湿度的变化值：

其中RH代表相对湿度，t代表时间，c_diff代表扩散系数，N_i是粒子i的邻居粒子集合，m是质量，是二阶微分算子，r是粒子的位置，W是光滑核函数即插值权重函数，h是光滑核函数W影响区域的半径；

把计算得到的湿度变化值加上该粒子湿度，得到该粒子处新的湿度属性；

ⅱ)各粒子的温度计算

粒子的温度变化量取决于邻居粒子的温度差：

其中T代表温度，t代表时间，c_d代表热传导系数,N_i是粒子i的邻居粒子集合，m是质量，是二阶微分算子，r是粒子的位置，W是光滑核函数即插值权重函数，h是光滑核函数W影响区域的半径；当相邻粒子材质不同时，为了保持边界处属性的连续性，热传导系数定为两种材质热传导系数的算术平均值；

把计算得到的温度变化值加上该粒子温度，得到该粒子处新的温度属性；

ⅲ)各粒子的凝结判断

使用August-Roche-Magnus方程对Clausius-Clapeyron方程近似：

式中T代表温度，T_dp为露点温度，RH代表相对湿度，常量系数：b＝17.67，c＝243.5℃；计算出对于每个气体粒子的露点温度；当气体粒子温度低于该粒子处的露点温度时，判断该气体粒子发生凝结现象，完成粒子模型的凝结现象的仿真。