一种涉及热传导与动态黏度的真实感流体仿真方法

摘要

本发明公开了一种涉及热传导与动态黏度的真实感流体仿真方法，其步骤为：1）基于光滑粒子流体动力学（SPH）模型，对流体内部及流体与外界的热传导过程进行离散建模，并根据相变焓对相变温度的影响以模拟相变过程；2）将动态黏度的计算引入以展示流体运动过程中的细节；3）调用PCISPH算法完成剩余的流体运动模拟过程；4）利用GPU加速的算法，在通用并行计算架构CUDA上将热传导、相变、流体黏度变化等过程并行处理，实现真实感相变流体的快速仿真。应用本发明能够真实、高效地模拟不同流体与外界的热传导及黏度变化过程，增强了现有方法中的模拟细节，提升了流体仿真的真实感。

说明书

技术领域

本发明属于计算机图形学领域，具体地说是一种涉及热传导与动态黏度的真实感流体仿真方法，其部分技术包括SPH流体模拟方法，热传导与相变原理以及GPU并行加速等。

背景技术

流体运动作为自然界普遍的物理现象，其仿真技术一直以来都是计算机图形学领域一个重要的研究方向，并且在影视制作、游戏开发、虚拟现实等领域都有着广泛的应用前景。基于物理的流体模拟主要通过求解流体力学中经典的纳维-斯托克斯(N-S)方程来进行，目前较流行的两种方法分别为拉格朗日法及欧拉法。其中拉格朗日法使用粒子系统模拟流体，概念易于理解，并且在处理细节及较大形变等方面表现突出，是目前应用较为广泛的一种模拟方法。

光滑粒子流体动力学(SPH)方法作为一种拉格朗日方法常用来模拟流体的运动规律。考虑到多数流体的不可压缩性，SPH的许多改进算法中逐渐完善了这一问题，如WCSPH、PCISPH方法等。其中PCISPH方法是目前应用较多的一种方法，其利用了预测-校正技术能够取得时间步长较大且每步迭代开销较低的优点，大幅提升了模拟效率，是一种非常适合模拟不可压缩流体的方法。

然而随着对真实感的进一步需求，仅仅考察流体的速度、位移已满足不了大众的视觉体验。因而流体的温度、状态等其他物理属性也必须考虑，以展示流体运动的细节。同时，当某些流体的温度发生改变时，其自身黏度往往也随之改变，尤其在诸如火山爆发、蜂蜜溶解、巧克力熔化等场景中体现得尤为明显。所以如何真实地模拟同种流体在不同温度下的黏性差异也是影响上述场景真实感仿真的一大因素，如果能很好地改善上述流体模拟过程中的细节，则对大众的视觉体验来说将会是一个巨大的提升。

发明内容

本发明的目的在于将现实生活中存在于流体中的热传导过程及黏度变化引入到基于物理的流体仿真过程中，提出了一种涉及热传导与动态黏度的真实感流体仿真方法，该方法引入了温度对流体黏度的影响，对流体热传导过程中的细节进行完善，获得了更加逼真的视觉效果。

实现本发明目的的具体技术方案是：

一种涉及热传导与动态黏度的真实感流体仿真方法，包括以下步骤：

a)基于光滑粒子流体动力学(SPH)模型模拟流体及热传导过程，具体包括：

ⅰ)邻居查找

首先将所有流体及固体用离散粒子表示，并查找得到与每个粒子i距离在光滑核半径h内的邻居粒子集合N_i；

ⅱ)各粒子导热系数计算

利用SPH算法的插值方式，考虑距离不同的邻居粒子对某个粒子i的影响，计算不同状态下粒子i的导热系数k_i，其具体公式为：

其中Phase_i表示粒子i的当前状态，m_j、ρ_j、k_j分别表示粒子i的邻居粒子j的质量、密度、导热系数，W_ij为形如的光滑核函数，它是关于粒子i与j的位置x_i与x_j之间距离||x_i-x_j||及光滑核半径h的函数；k_fluid与k_solid分别为该种流体在液态及固态下实际的导热系数；

iii)各粒子温度变化率计算

首先计算流体内部的热传导，对每个粒子i，其温度随时间变化率计算公式为：

$\frac{{\mathrm{dT}}_i}{dt}=\Delta \left(k_i{T}_i\right)$

其中T_i及k_i分别为粒子i的温度及导热系数，Δ为拉普拉斯算子；同时Δ(k_iT_i)的计算可通过公式Δ(k_iT_i)＝及得到，其中表示梯度；

其次计算流体与外界的热传导，对每个粒子i，其温度随时间变化率计算公式为：

${\left(\frac{{\mathrm{dT}}_i}{dt}\right)}_{ext}=\underset{j}{\Sigma }{k}_i\frac{m_J{R_i}^2}{{\rho }_i{\rho }_j}(T_b-T_i)W_{ij}$

其中T_b为固体粒子温度，R_i可通过计算得到，ρ₀为流体的静态密度；

然后更新粒子i的温度T_i，具体公式为：

$T_i(t+\Delta t)=T_i\left(t\right)+[{\left(\frac{{\mathrm{dT}}_i}{dt}\right)}_{ext}+\frac{{\mathrm{dT}}_i}{dt}]\Delta t$

其中Δt为时间步长，T_i(t)与T_i(t+Δt)分别为粒子i在时刻t与t+Δt的温度值；

iv)相变处理

由于相变焓的存在，分别设置两个温度阈值T_melt和T_solid；当粒子温度T_i大于T_melt时，粒子转化为液体粒子，当T_i小于T_solid时，粒子转化为固体粒子；当T_i处于T_solid与T_melt之间时，则将其看作为临界状态的流体粒子；

b)动态黏度计算的引入，具体包括：

在上述过程得到每个粒子i的最新时刻温度T_i后，由于不同温度对流体黏性系数存在一定的影响，在SPH模型中引入动态黏度的计算，实现流体温度变化过程中细节的逼真模拟，其中温度T_i与流体黏性系数μ_i的变化关系具体为：

log(log(μ_i+γ))＝q-y>i)

其中γ为常数，通常取值在0.6到0.9之间，q与y则为与流体自身有关的特定参数，为定值；

c)PCISPH算法调用

利用PCISPH算法，对每个粒子i分别计算其所受到压力、黏性力以及外力而产生的加速度a_i^pre、a_i^vis及a_i^other，其计算公式分别为：

${a_i}^{pre}=-\frac{1}{m_i}\underset{j}{\Sigma }(m_j\frac{p_i+p_j}{2{\rho }_j}▿W_{ij})$

${a_i}^{vis}=\frac{1}{2m_i}\underset{j}{\Sigma }\left(m_j\right({\mu }_i+{\mu }_j\left)\frac{v_j-v_i}{{\rho }_j}{▿}^2{W}_{ij}\right)$

a_i^other＝g

其中，其中m、p、ρ、μ、v分别表示粒子的质量、压强、密度、黏性系数及速度，g为重力加速度；然后通过粒子i所受的加速度计算得到其在下一时刻的位置以模拟流体运动变化规律；

d)GPU并行加速，具体包括：

利用通用并行计算架构CUDA，将所有粒子的物理属性传输至GPU，并对粒子进行排序，加速邻居查找过程；对于每个粒子，为其开辟一个独立的线程以计算其密度、温度、导热系数、加速度，并更新流体粒子的速度与位置，从而并行地提升流体仿真的效率。

本发明的有益效果：

现有的流体模拟方法对流体温度变化之后相应的黏度变化关注度并不高，尤其是利用SPH方法的流体模拟方法中。本发明充分借鉴了现有的物理原理，引入了温度对流体黏度的影响，对流体热传导过程中的细节进行完善，获得了更加逼真的视觉效果。

本发明在对热传导的离散建模中考虑了相变焓对相变过程的影响，并且对不同状态的流体粒子采用各自的导热系数，较现有方法更为严谨，尊于事实，保证了模拟过程的真实感。同时，借助GPU加速技术本发明提高了计算效率，在实验过程中较之CPU方法取得了最高接近16倍的加速比。

总之，应用本发明可以快速地模拟流体的热传导及黏度变化过程，在细节、真实感体验及时间效率上都有所提升。

附图说明

图1为本发明方法模拟牛奶状液体内部热传导示意图；

图2为本发明方法模拟蜡烛熔化及凝固过程的示意图；

图3为本发明方法模拟蜡烛熔化及凝固过程中各粒子数随时间的变化关系示意图；

图4为传统恒定黏度方法与本发明方法模拟巧克力加热熔化过程的对比图；

图5为本发明方法模拟不同温度火山熔岩从山坡上下滑的对比图。

具体实施方式

本发明包括以下步骤：

1)基于光滑粒子流体动力学(SPH)模型模拟流体及热传导过程：

由于相变在一个较小的温度范围内进行，设置相变过程中的临界状态粒子，从而设置流体粒子在液态、临界态及固态之间的转换条件。为不同状态的流体粒子设置其导热系数，对同种流体内部及流体与外界接触面之间的热传导过程进行离散建模。

2)动态黏度计算的引入：

由于不同温度对流体黏性系数存在影响，在SPH模型中引入动态黏度的计算，实现流体温度变化过程中细节的逼真模拟。

3)PCISPH算法调用：

利用PCISPH算法保证流体的不可压缩性，并以此为框架进行流体模拟过程。

4)GPU并行加速，具体包括：

利用通用并行计算架构CUDA，发挥GPU加速的优势，将流体模拟过程中的热传导、相变、流体黏度变化等并行处理，从而提升整体运行效率。

SPH基本框架：

SPH方法将流体看成是由许多粒子构成的系统，每个粒子i有各自的物理属性，包括质量m_i、密度ρ_i、体积V_i、压强p_i、速度v_i、位置x_i等。该方法利用邻居粒子属性插值的方法计算粒子i的各物理属性，具体形式如下：

$A_i=\underset{j}{\Sigma }\frac{m_j}{{\rho }_j}{A}_j{W}_{ij}---\left(1\right)$

其中A_i表示粒子i的某一物理属性值，j表示粒子i支持域范围内的粒子。W_ij为形如的光滑核函数，它是关于粒子i与粒子j之间距离||x_i-x_j||及光滑核半径h的函数。同样地，及也可用类似的方法插值得到。

随着时间t的推移，流体粒子相应的物理属性随之改变，根据公式(2)更新流体粒子的位置：

$\frac{{\mathrm{dx}}_i}{dt}=v_i---\left(2\right)$

其中粒子i的速度v_i体现了流体运动的规律，其变化率即由拉格朗日形式的N-S方程决定，具体公式为：

$\frac{{\mathrm{dv}}_i}{dt}=-\frac{1}{{\rho }_i}▿p_i+\mu {▿}^2{v}_i+\frac{F_i^{other}}{m_i}---\left(3\right)$

其中μ表示流体的整体黏度系数，表示粒子i受到的外力。公式(3)右端三项分别代表粒子i所受到压力、黏性力以及外力而产生的加速度，即：

${a_i}^{pre}=-\frac{1}{m_i}\underset{j}{\Sigma }(m_j\frac{p_i+p_j}{2{\rho }_j}▿W_{ij})---\left(4\right)$

${a_i}^{vis}=\frac{\mu }{m_i}\underset{j}{\Sigma }\left(m_j\frac{v_j-v_i}{{\rho }_j}{▿}^2{W}_{ij}\right)---\left(5\right)$

a_i^other＝g(6)

其中g为重力加速度。

根据公式(2)～(6)即可求得流体粒子在每一时刻的速度v_i(t)与位置x_i(t)。

传统的SPH方法通常利用理想气体状态方程(7)计算粒子i所受压强。其中ρ₀代表流体的静态密度，K为和流体相关的常数，一般只与温度有关。

p_i＝K(ρ_i-ρ₀)(7)

然而这种计算压强的方法无法保证流体的不可压缩性，因而在流体模拟过程中常会出现一些失真的现象。目前已有几种高效的不可压缩方法可以很好地改善这个问题。

PCISPH作为一种不可压缩的流体模拟方法，主要运用了预测-校正的思想方法。对于粒子i，该方法首先计算在除压力以外其他所有力作用下产生的加速度以预测粒子的位置与速度，并利用公式(1)计算粒子i的预测密度值从而得到其与静态密度之间的误差并以此更新压强其中δ的计算方式如下：

$\delta =\frac{-1}{\beta (-{\Sigma }_j▿W_{ij}·{\Sigma }_j▿W_{ij}-{\Sigma }_j(▿W_{ij}·▿W_{ij}))}---\left(8\right)$

$\beta ={\mathrm{\Delta t}}^2{m}^2\frac{2}{{{\rho }_0}^2}---\left(9\right)$

其中Δt表示时间步长，m为流体粒子质量。

然后利用更新后的压强重新计算压力，并作用在粒子i上产生位移。以上预测-校正的步骤会循环进行直至密度误差小于事先给定的某个阈值，即认为流体在模拟过程中保证了其不可压缩的性质。除此以外，PCISPH方法可以使用比传统SPH更大的时间步长，因而能够大大提高模拟的速率。因此，本发明选择PCISPH作为基础框架。

本发明的相变模拟过程具体为：

通常情况下相变被认为发生在单一温度下，即理想相变，此时相变焓为常数。然而由于相变焓的变化，导致相变其实是在一个小的温度范围内进行。由于这一情况的存在，本发明在模拟相变过程中根据温度变化设置了临界状态粒子，并且采用如下方法来处理粒子状态间的转化。

首先分别设置两个温度阈值T_melt和T_solid。当粒子温度T_i大于T_melt时，粒子转化为液体粒子，当T_i小于T_solid时，粒子转化为固体粒子。当T_i处于T_solid与T_melt之间时，则将其看作为临界状态的流体粒子。具体过程如下算法1所示，其中Phase_i表示粒子i的状态。

算法1.相变处理算法

若T_i＞T_melt

则Phase_i＝液态

否则若T_i＜T_solid

则Phase_i＝固态

否则Phase_i＝临界态

本发明对热传导过程的离散建模具体为：

对于发生在流体内部的热传导，一般方法中计算单个粒子温度变化率的方式为：

$\frac{{\mathrm{dT}}_i}{dt}={\mathrm{k\Delta T}}_i---\left(10\right)$

其中k为该物质整体的导热系数。

然而，由于导热系数与物质的状态有关，本发明针对每个流体粒子增加了独立的导热系数k_i。由于处于相变临界状态的流体导热系数很难通过具体实验测量得到，在此利用SPH插值计算的思想近似处于临界状态的粒子i的导热系数为：

$k_i=\underset{j}{\Sigma }\frac{m_j}{{\rho }_j}{k}_j{W}_{ij}---\left(11\right)$

对于临界状态的粒子，考虑光滑核半径内其它粒子对它的影响。从而对于不同状态的粒子，其导热系数可由如下公式得到：

其中k_fluid与k_solid分别为液态与固态下的导热系数。于是将公式(10)转化为：

$\frac{{\mathrm{dT}}_i}{dt}=\Delta \left(k_i{T}_i\right)---\left(13\right)$

其次将公式(13)表示成离散形式，得到Δ(k_iT_i)的计算公式为：

$▿\left(k_i{T}_i\right)=\underset{j}{\Sigma }{m}_j\frac{k_j{T}_j-k_i{T}_i}{{\rho }_j}▿W_{ij}---\left(14\right)$

$\Delta \left(k_i{T}_i\right)=\underset{j}{\Sigma }{m}_j\frac{▿\left(k_j{T}_j\right)}{{\rho }_j}▿W_{ij}---\left(15\right)$

至此即完成了流体内部的热传导的离散建模。

对于流体与外界的热传导，一般主要发生在固液交互过程中，因而本发明忽略了空气对流体的影响。对于固液交界面处的热传导，一般方法将固体作为一个整体，通过如下两个公式进行计算：

${\left(\frac{{\mathrm{dT}}_i}{dt}\right)}_{ext}=k_i(T_b-T_i)\frac{{R_i}^2}{{\rho }_i}---\left(16\right)$

$\frac{4}{3}{{\mathrm{\pi R}}_i}^3=\frac{m_i}{{\rho }_0}---\left(17\right)$

其中表示流体粒子i与外界热传导而造成的温度变化率，T_b表示外界固体粒子的温度。由此可知流体粒子与外界热传导而产生的温度变化率与两种物质的温度差及接触面积成正比，与自身密度成反比。

在本发明中，首先将固体离散成粒子的形式，考虑局部固体粒子b对流体粒子i温度的单独影响以及每个粒子各自不同的导热系数k_i，并将公式(16)离散为：

${\left(\frac{{\mathrm{dT}}_i}{dt}\right)}_{ext}=\underset{j}{\Sigma }{k}_i\frac{m_j{R_i}^2}{{\rho }_i{\rho }_j}(T_b-T_i)W_{ij}---\left(18\right)$

类似地，固体粒子温度T_b的变化率也可通过以上公式变换求得。因而至此完成了流体与外界的热传导离散建模过程。

本发明的动态黏度计算具体为：

对于同种流体，其黏度往往随着温度的升高而降低，当温度较低时往往表现得较为黏稠，这一现象在一些温度变化范围较大的流体中更为明显，例如火山熔岩。针对这一现象，在SPH方法中引入了动态黏度的概念，利用温度与流体黏性力之间的物理关系，将其应用到单个粒子上得到：

log(log(μ_i+γ))＝q-y>i)(19)

其中μ_i表示粒子i的黏度系数，γ为常数，通常取值在0.6到0.9之间，q与y则为与流体自身有关的特定参数，为定值，一般通过具体实验测得，本发明中分别取值为2.5左右与1.0左右。

由于粒子的黏度系数发生变化，相应的由黏度力产生的加速度公式(5)也可修改为：

${a_i}^{vis}=\frac{1}{2m_i}\underset{j}{\Sigma }\left(m_j\right({\mu }_i+{\mu }_j\left)\frac{v_j-v_i}{{\rho }_j}{▿}^2{W}_{ij}\right)---\left(20\right)$

本发明的具体实现过程为：

将所有的流体与外界固体边界都用粒子离散表示，并为每个粒子附上初始位置、速度、温度等属性。本发明使用蛙跳法更新流体粒子速度与位置，对于时间步长则使用CFL条件确定。

利用通用并行计算架构CUDA，将所有粒子的物理属性传输至GPU，并对粒子进行排序，加速邻居查找过程；对于每个粒子，为其开辟一个独立的线程以计算其密度、温度、导热系数、加速度，并更新流体粒子的速度与位置。

同时，由于上述物理属性的计算过程与SPH方法中其他问题如不可压缩性、边界问题等的处理过程相互独立，因而可以很方便地移植到其他SPH改进算法中。基于PCISPH算法的具体流程如下算法2所示。

算法2.涉及热传导与动态黏度的流体仿真算法

流体模拟过程中：

对所有粒子i

查找邻居集合N_i

对所有粒子i

利用公式(12)计算其导热系数k_i

对所有粒子i

利用公式(16)计算

对所有粒子i

利用公式(13)(15)计算dT_i/dt

利用公式(18)计算(dT_i/dt)_ext

更新温度

调用算法1相变处理

若Phase_i≠固态

则利用公式(19)计算黏度系数μ_i

执行PCISPH算法，其中黏力加速度计算采用公式(20)

在每一个时间步长的计算中，热传导模型与动态黏度计算模块在查找完邻居粒子之后进行，继而执行PCISPH算法，值得注意的是其中由黏度力产生的加速度需要利用公式(20)计算得到。最后，根据模拟得到的粒子位置信息采用Marching Cubes算法计算流体的表面，然后进行后期渲染工作。

实施例

本发明涉及热传导与动态黏度的真实感流体仿真方法，效果展示说明如下：

图1展示了温度较高的兔子形状的牛奶落入温度较低的牛奶中的效果。该场景共使用了87164个流体粒子，从左至右每列分别为模拟的第0、80、160帧，其中上排为温度场图，下排为效果图。从图中可以直观地看到热传导的整个过程，温度较高的牛奶在模拟之初全部集中于容器中央，而后随时间的推移热量逐渐向外围传递，整体颜色趋于均匀。此场景主要呈现了同种流体内部的热传导过程。

图2展示了蜡烛熔化及凝固的过程，从左至右分别为模拟的第200、800、1400、2000、2600、3200帧。从刚点燃蜡烛起，其上部逐渐熔化形成液态的蜡油。在蜡油累积到一定程度之后，开始沿着蜡烛表面向下流动，最终遇到底部温度较低的大理石桌面进行热传导导致蜡油温度降低从而发生凝固。

为了更具体地反映该场景中存在的熔化再凝固过程，在此记录了场景中各状态粒子数随模拟时间的变化情况，如图3所示。从图中可以看到在模拟之初，流体粒子及临界状态粒子数由于热传导的发生而显著增加，而固态粒子数明显减少，即主要为蜡烛的熔化过程。从第3000帧左右开始，流体粒子及临界状态粒子数增长趋于平缓，最后呈下降趋势，固态粒子数逐渐增加，其主要原因即为蜡油与温度较低的固体边界接触而发生了凝固现象。结合效果图与具体统计数据，可以清晰地展现出蜡烛熔化及凝固的整个过程。

图4对比了使用恒定黏度系数与本发明涉及动态黏度两者模拟结果之间的差异，从左至右分别为模拟的第20、70、120、190帧。在持续加热的情况下，上排使用恒定的黏度系数，在不同温度下巧克力始终处于黏度较大的状态，表现得较为黏稠，不易流动。而下排使用了本文动态黏度的计算方法，巧克力的黏度在熔化之初与上排黏度相仿，在达到较高温度后表现得较稀，表面积增大，最终几乎占满了整个锅底。这与真实情况下巧克力加热熔化过程中由黏稠变为稀薄的细节更为接近。从该场景中可以看到本发明的动态黏度计算可以真实地反映流体运动过程中的细节，达到更加逼真的效果。

为了体现不同温度流体的黏度差异，图5分别比较了初始温度为700℃(上排)及1200℃(下排)的火山熔岩从同一山坡上下滑的过程。根据动态黏度计算方法，由于初始温度的不同，两种情况下火山熔岩的平均黏度会有显著差异。

因此，从图5中可以观察到：具有较高温度的熔岩下滑速度较快，距离较远，而温度较低的熔岩则由于其黏度较大下滑明显较慢。此实验反映了实际情况下流体黏度与温度之间的关系，验证了流体动态黏度计算的合理性。

最后，对各场景分别利用GPU及CPU模拟，在下表中给出了各自每帧所需的平均时间，其中最后一列为本发明在GPU上获得的加速效果。可以看到通过并行化处理，本发明获得了最高接近16倍的加速比，提升了模拟的效率。

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