一种涉及热传导与动态黏度的真实感流体仿真方法
【专利摘要】本发明公开了一种涉及热传导与动态黏度的真实感流体仿真方法，其步骤为：1）基于光滑粒子流体动力学（SPH）模型，对流体内部及流体与外界的热传导过程进行离散建模，并根据相变焓对相变温度的影响以模拟相变过程；2）将动态黏度的计算引入以展示流体运动过程中的细节；3）调用PCISPH算法完成剩余的流体运动模拟过程；4）利用GPU加速的算法，在通用并行计算架构CUDA上将热传导、相变、流体黏度变化等过程并行处理，实现真实感相变流体的快速仿真。应用本发明能够真实、高效地模拟不同流体与外界的热传导及黏度变化过程，增强了现有方法中的模拟细节，提升了流体仿真的真实感。
【专利说明】
一种涉及热传导与动态黏度的真实感流体仿真方法
技术领域
[0001] 本发明属于计算机图形学领域，具体地说是一种涉及热传导与动态黏度的真实感 流体仿真方法，其部分技术包括SPH流体模拟方法，热传导与相变原理以及GPU并行加速等。
【背景技术】
[0002] 流体运动作为自然界普遍的物理现象，其仿真技术一直以来都是计算机图形学领 域一个重要的研究方向，并且在影视制作、游戏开发、虚拟现实等领域都有着广泛的应用前 景。基于物理的流体模拟主要通过求解流体力学中经典的纳维-斯托克斯(N-S)方程来进 行，目前较流行的两种方法分别为拉格朗日法及欧拉法。其中拉格朗日法使用粒子系统模 拟流体，概念易于理解，并且在处理细节及较大形变等方面表现突出，是目前应用较为广泛 的一种模拟方法。
[0003] 光滑粒子流体动力学(SPH)方法作为一种拉格朗日方法常用来模拟流体的运动规 律。考虑到多数流体的不可压缩性，SPH的许多改进算法中逐渐完善了这一问题，如WCSPH、 PCISPH方法等。其中PCISPH方法是目前应用较多的一种方法，其利用了预测-校正技术能够 取得时间步长较大且每步迭代开销较低的优点，大幅提升了模拟效率，是一种非常适合模 拟不可压缩流体的方法。
[0004] 然而随着对真实感的进一步需求，仅仅考察流体的速度、位移已满足不了大众的 视觉体验。因而流体的温度、状态等其他物理属性也必须考虑，以展示流体运动的细节。同 时，当某些流体的温度发生改变时，其自身黏度往往也随之改变，尤其在诸如火山爆发、蜂 蜜溶解、巧克力熔化等场景中体现得尤为明显。所以如何真实地模拟同种流体在不同温度 下的黏性差异也是影响上述场景真实感仿真的一大因素，如果能很好地改善上述流体模拟 过程中的细节，则对大众的视觉体验来说将会是一个巨大的提升。
【发明内容】
[0005] 本发明的目的在于将现实生活中存在于流体中的热传导过程及黏度变化引入到 基于物理的流体仿真过程中，提出了一种涉及热传导与动态黏度的真实感流体仿真方法， 该方法引入了温度对流体黏度的影响，对流体热传导过程中的细节进行完善，获得了更加 逼真的视觉效果。
[0006] 实现本发明目的的具体技术方案是：
[0007] -种涉及热传导与动态黏度的真实感流体仿真方法，包括以下步骤：
[0008] a)基于光滑粒子流体动力学(SPH)模型模拟流体及热传导过程，具体包括：
[0009] i)邻居查找
[0010] 首先将所有流体及固体用离散粒子表示，并查找得到与每个粒子i距离在光滑核 半径h内的邻居粒子集合Ni;
[0011] ii)各粒子导热系数计算
[0012]利用SPH算法的插值方式，考虑距离不同的邻居粒子对某个粒子i的影响，计算不 同状态下粒子i的导热系数lu，其具体公式为：
[0014]其中Phasei表示粒子i的当前状态，mj、Pj、kj分别表示粒子i的邻居粒子j的质量、 密度、导热系数：
i的光滑核函数，它是关于粒子i与j的位置^与&之 间距离I |Xl_Xj| I及光滑核半径h的函数;kfluid与kscllld分别为该种流体在液态及固态下实际 的导热系数；
[0015] iii)各粒子温度变化率计算
[0016] 首先计算流体内部的热传导，对每个粒子i，其温度随时间变化率计算公式为：
[0018] 其中^及匕分别为粒子i的温度及导热系数，Δ为拉普拉斯算子；同时Δ (luTO的 计算可通过公式
得到，其中I7表示梯 度；
[0019] 其次计算流体与外界的热传导，对每个粒子i，其温度随时间变化率计算公式为：
[0021] 其中Tb为固体粒子温度，心可通过
计算得到，Po为流体的静态密度；
[0022] 然后更新粒子i的温度，具体公式为：
[0024] 其中Δ t为时间步长，Ti⑴与Ti(t+ Δ t)分别为粒子i在时亥[|t与t+ Δ t的温度值；
[0025] iv)相变处理
[0026] 由于相变焓的存在，分别设置两个温度阈值Tmelt和Tscllid;当粒子温度Ti大于T melt 时，粒子转化为液体粒子，当Ti小于Tsolid时，粒子转化为固体粒子；当Ti处于Tsolid与Tmelt之 间时，则将其看作为临界状态的流体粒子；
[0027] b)动态黏度计算的引入，具体包括：
[0028] 在上述过程得到每个粒子i的最新时刻温度1\后，由于不同温度对流体黏性系数 存在一定的影响，在SPH模型中引入动态黏度的计算，实现流体温度变化过程中细节的逼真 模拟，其中温度h与流体黏性系数&的变化关系具体为：
[0029] log(log(yi+y )) = q-y log(Ti)
[0030] 其中γ为常数，通常取值在0.6到0.9之间，q与y则为与流体自身有关的特定参数， 为定值；
[0031] c)PCISPH 算法调用
[0032]利用PCISPH算法，对每个粒子i分别计算其所受到压力、黏性力以及外力而产生的 加速度31@、&1'^及&14'其计算公式分别为 ：
[0035] aiother = g
[0036] 其中，其中!11、？、0^^分别表示粒子的质量、压强、密度、黏性系数及速度4为重力 加速度;然后通过粒子i所受的加速度计算得到其在下一时刻的位置以模拟流体运动变化 规律；
[0037] d)GPU并行加速，具体包括：
[0038]利用通用并行计算架构⑶DA，将所有粒子的物理属性传输至GPU，并对粒子进行排 序，加速邻居查找过程;对于每个粒子，为其开辟一个独立的线程以计算其密度、温度、导热 系数、加速度，并更新流体粒子的速度与位置，从而并行地提升流体仿真的效率。
[0039]本发明的有益效果：
[0040] 现有的流体模拟方法对流体温度变化之后相应的黏度变化关注度并不高，尤其是 利用SPH方法的流体模拟方法中。本发明充分借鉴了现有的物理原理，引入了温度对流体黏 度的影响，对流体热传导过程中的细节进行完善，获得了更加逼真的视觉效果。
[0041] 本发明在对热传导的离散建模中考虑了相变焓对相变过程的影响，并且对不同状 态的流体粒子采用各自的导热系数，较现有方法更为严谨，尊于事实，保证了模拟过程的真 实感。同时，借助GPU加速技术本发明提高了计算效率，在实验过程中较之CPU方法取得了最 高接近16倍的加速比。
[0042]总之，应用本发明可以快速地模拟流体的热传导及黏度变化过程，在细节、真实感 体验及时间效率上都有所提升。
【附图说明】
[0043] 图1为本发明方法模拟牛奶状液体内部热传导示意图；
[0044] 图2为本发明方法模拟蜡烛熔化及凝固过程的示意图；
[0045] 图3为本发明方法模拟蜡烛熔化及凝固过程中各粒子数随时间的变化关系示意 图；
[0046] 图4为传统恒定黏度方法与本发明方法模拟巧克力加热熔化过程的对比图；
[0047] 图5为本发明方法模拟不同温度火山熔岩从山坡上下滑的对比图。
【具体实施方式】
[0048]本发明包括以下步骤：
[0049] 1)基于光滑粒子流体动力学(SPH)模型模拟流体及热传导过程：
[0050] 由于相变在一个较小的温度范围内进行，设置相变过程中的临界状态粒子，从而 设置流体粒子在液态、临界态及固态之间的转换条件。为不同状态的流体粒子设置其导热 系数，对同种流体内部及流体与外界接触面之间的热传导过程进行离散建模。
[0051] 2)动态黏度计算的引入：
[0052]由于不同温度对流体黏性系数存在影响，在SPH模型中引入动态黏度的计算，实现 流体温度变化过程中细节的逼真模拟。
[0053] 3)PCISPH 算法调用：
[0054]利用PCISPH算法保证流体的不可压缩性，并以此为框架进行流体模拟过程。
[0055] 4)GPU并行加速，具体包括：
[0056]利用通用并行计算架构CUDA，发挥GPU加速的优势，将流体模拟过程中的热传导、 相变、流体黏度变化等并行处理，从而提升整体运行效率。
[0057] SPH基本框架：
[0058] SPH方法将流体看成是由许多粒子构成的系统，每个粒子i有各自的物理属性，包 括质量nu、密度Pl、体积I、压强Pl、速度Vl、位置^等。该方法利用邻居粒子属性插值的方法 计算粒子i的各物理属性，具体形式如下：
[0060]其中Ai表示粒子i的某一物理属性值，j表示粒子i支持域范围内的粒子。Wij为形如 i的光滑核函数，它是关于粒子i与粒子j之间距离I |Xl-x」I及光滑核半径h的 函数。同样地，t7也及户4也可用类似的方法插值得到。
[0061 ]随着时间t的推移，流体粒子相应的物理属性随之改变，根据公式(2)更新流体粒 子的位置：
[0063]其中粒子i的速度Vi体现了流体运动的规律，其变化率即由拉格朗日形式的N-S方 程决定，具体公式为：
[0065]其中μ表示流体的整体黏度系数，?表示粒子i受到的外力。公式（3)右端三项 分别代表粒子i所受到压力、黏性力以及外力而产生的加速度，即：
[0068] aiother = g (6)
[0069] 其中g为重力加速度。
[0070]根据公式(2)~(6)即可求得流体粒子在每一时刻的速度^(〇与位置Xl(t)。
[0071] 传统的SPH方法通常利用理想气体状态方程(7)计算粒子i所受压强。其中Po代表 流体的静态密度，K为和流体相关的常数，一般只与温度有关。
[0072] pi = K(Pi-p〇) (7)
[0073] 然而这种计算压强的方法无法保证流体的不可压缩性，因而在流体模拟过程中常 会出现一些失真的现象。目前已有几种高效的不可压缩方法可以很好地改善这个问题。 [0074] PCISPH作为一种不可压缩的流体模拟方法，主要运用了预测-校正的思想方法。对 于粒子i，该方法首先计算在除压力以外其他所有力作用下产生的加速度以预测粒子的位 置与速度，并利用公式（1)计算粒子i的预测密度值P〖，从而得到其与静态密度之间的误差 一 flo*并以此更新压强Pf+= 5p;rori，其中δ的计算方式如下：
[0077] 其中△ t表示时间步长，m为流体粒子质量。
[0078] 然后利用更新后的压强重新计算压力，并作用在粒子i上产生位移。以上预测-校 正的步骤会循环进行直至密度误差小于事先给定的某个阈值，即认为流体在模拟过程中保 证了其不可压缩的性质。除此以外，PCISPH方法可以使用比传统SPH更大的时间步长，因而 能够大大提高模拟的速率。因此，本发明选择PCISPH作为基础框架。
[0079]本发明的相变模拟过程具体为：
[0080]通常情况下相变被认为发生在单一温度下，即理想相变，此时相变焓为常数。然而 由于相变焓的变化，导致相变其实是在一个小的温度范围内进行。由于这一情况的存在，本 发明在模拟相变过程中根据温度变化设置了临界状态粒子，并且采用如下方法来处理粒子 状态间的转化。
[0081 ] 首先分别设置两个温度阈值Tmelt和Tscllid。当粒子温度Ti大于T melt时，粒子转化为液 体粒子，当Ti小于Tsolid时，粒子转化为固体粒子。当Ti处于Tsolid与Imelt之间时，则将其看作 为临界状态的流体粒子。具体过程如下算法1所示，其中示粒子i的状态。
[0082]算法1.相变处理算法
[0083] 若 Ti>Tmeit
[0084] 则 Phasei =液态
[0085] 否则若Ti<Ts〇iid
[0086] 则 Phasei=固态
[0087] 否则phasei =临界态
[0088] 本发明对热传导过程的离散建模具体为：
[0089] 对于发生在流体内部的热传导，一般方法中计算单个粒子温度变化率的方式为：
[0091] 其中k为该物质整体的导热系数。
[0092] 然而，由于导热系数与物质的状态有关，本发明针对每个流体粒子增加了独立的 导热系数1^。由于处于相变临界状态的流体导热系数很难通过具体实验测量得到，在此利 用SPH插值计算的思想近似处于临界状态的粒子i的导热系数为：
[0094]对于临界状态的粒子，考虑光滑核半径内其它粒子对它的影响。从而对于不同状 态的粒子，其导热系数可由如下公式得到：
[0096]其中kfluid与kscllld分别为液态与固态下的导热系数。于是将公式(10)转化为：
[0098]其次将公式(13)表示成离散形式，得到△(kdO的计算公式为：
[0101]至此即完成了流体内部的热传导的离散建模。
[0102]对于流体与外界的热传导，一般主要发生在固液交互过程中，因而本发明忽略了 空气对流体的影响。对于固液交界面处的热传导，一般方法将固体作为一个整体，通过如下 两个公式进行计算：
[0105]
表示流体粒子i与外界热传导而造成的温度变化率，Tb表示外界固体粒 子的温度。由此可知流体粒子与外界热传导而产生的温度变化率与两种物质的温度差及接 触面积成正比，与自身密度成反比。
[0106] 在本发明中，首先将固体离散成粒子的形式，考虑局部固体粒子b对流体粒子i温 度的单独影响以及每个粒子各自不同的导热系数lu，并将公式(16)离散为：
[0108] 类似地，固体粒子温度Tb的变化率也可通过以上公式变换求得。因而至此完成了 流体与外界的热传导离散建模过程。
[0109] 本发明的动态黏度计算具体为：
[0110] 对于同种流体，其黏度往往随着温度的升高而降低，当温度较低时往往表现得较 为黏稠，这一现象在一些温度变化范围较大的流体中更为明显，例如火山熔岩。针对这一现 象，在SPH方法中引入了动态黏度的概念，利用温度与流体黏性力之间的物理关系，将其应 用到单个粒子上得到：
[0111] log(log(yi+y )) = q-y log(Ti) (19)
[0112] 其中示粒子i的黏度系数，γ为常数，通常取值在0.6到0.9之间，q与y则为与 流体自身有关的特定参数，为定值，一般通过具体实验测得，本发明中分别取值为2.5左右 与1.0左右。
[0113]由于粒子的黏度系数发生变化，相应的由黏度力产生的加速度公式(5)也可修改 为：
[0115] 本发明的具体实现过程为：
[0116] 将所有的流体与外界固体边界都用粒子离散表示，并为每个粒子附上初始位置、 速度、温度等属性。本发明使用蛙跳法更新流体粒子速度与位置，对于时间步长则使用CFL 条件确定。
[0117]利用通用并行计算架构⑶DA，将所有粒子的物理属性传输至GPU，并对粒子进行排 序，加速邻居查找过程;对于每个粒子，为其开辟一个独立的线程以计算其密度、温度、导热 系数、加速度，并更新流体粒子的速度与位置。
[0118] 同时，由于上述物理属性的计算过程与SPH方法中其他问题如不可压缩性、边界问 题等的处理过程相互独立，因而可以很方便地移植到其他SPH改进算法中。基于PCISPH算法 的具体流程如下算法2所示。
[0119] 算法2.涉及热传导与动态黏度的流体仿真算法
[0120] 流体模拟过程中：
[0121] 对所有粒子i
[0122] 查找邻居集合Ni
[0123] 对所有粒子i
[0124] 利用公式(12)计算其导热系数匕
[0125] 对所有粒子i
[0126] 利用公式（16)计算对&拟
[0127] 对所有粒子i
[0128] 利用公式
[0129] 利用公式（18)计算(dVdtUt
[0131] 调用算法1相变处理
[0132] 若Phasei辛固态
[0133] 则利用公式(19)计算黏度系数yi
[0134] 执行PCISPH算法，其中黏力加速度计算采用公式(20)
[0135] 在每一个时间步长的计算中，热传导模型与动态黏度计算模块在查找完邻居粒子 之后进行，继而执行PCISPH算法，值得注意的是其中由黏度力产生的加速度需要利用公式 (20)计算得到。最后，根据模拟得到的粒子位置信息采用Marching Cubes算法计算流体的 表面，然后进行后期渲染工作。
[0136] 实施例
[0137] 本发明涉及热传导与动态黏度的真实感流体仿真方法，效果展示说明如下：
[0138] 图1展示了温度较高的兔子形状的牛奶落入温度较低的牛奶中的效果。该场景共 使用了87164个流体粒子，从左至右每列分别为模拟的第0、80、160帧，其中上排为温度场 图，下排为效果图。从图中可以直观地看到热传导的整个过程，温度较高的牛奶在模拟之初 全部集中于容器中央，而后随时间的推移热量逐渐向外围传递，整体颜色趋于均匀。此场景 主要呈现了同种流体内部的热传导过程。
[0139] 图2展示了蜡烛熔化及凝固的过程，从左至右分别为模拟的第200、800、1400、 2000、2600、3200帧。从刚点燃蜡烛起，其上部逐渐熔化形成液态的蜡油。在蜡油累积到一定 程度之后，开始沿着蜡烛表面向下流动，最终遇到底部温度较低的大理石桌面进行热传导 导致蜡油温度降低从而发生凝固。
[0140] 为了更具体地反映该场景中存在的熔化再凝固过程，在此记录了场景中各状态粒 子数随模拟时间的变化情况，如图3所示。从图中可以看到在模拟之初，流体粒子及临界状 态粒子数由于热传导的发生而显著增加，而固态粒子数明显减少，即主要为蜡烛的熔化过 程。从第3000帧左右开始，流体粒子及临界状态粒子数增长趋于平缓，最后呈下降趋势，固 态粒子数逐渐增加，其主要原因即为蜡油与温度较低的固体边界接触而发生了凝固现象。 结合效果图与具体统计数据，可以清晰地展现出蜡烛熔化及凝固的整个过程。
[0141] 图4对比了使用恒定黏度系数与本发明涉及动态黏度两者模拟结果之间的差异， 从左至右分别为模拟的第20、70、120、190帧。在持续加热的情况下，上排使用恒定的黏度系 数，在不同温度下巧克力始终处于黏度较大的状态，表现得较为黏稠，不易流动。而下排使 用了本文动态黏度的计算方法，巧克力的黏度在熔化之初与上排黏度相仿，在达到较高温 度后表现得较稀，表面积增大，最终几乎占满了整个锅底。这与真实情况下巧克力加热熔化 过程中由黏稠变为稀薄的细节更为接近。从该场景中可以看到本发明的动态黏度计算可以 真实地反映流体运动过程中的细节，达到更加逼真的效果。
[0142] 为了体现不同温度流体的黏度差异，图5分别比较了初始温度为700°C (上排）及 1200°C (下排）的火山熔岩从同一山坡上下滑的过程。根据动态黏度计算方法，由于初始温 度的不同，两种情况下火山熔岩的平均黏度会有显著差异。
[0143] 因此，从图5中可以观察到:具有较高温度的熔岩下滑速度较快，距离较远，而温度 较低的熔岩则由于其黏度较大下滑明显较慢。此实验反映了实际情况下流体黏度与温度之 间的关系，验证了流体动态黏度计算的合理性。
[0144] 最后，对各场景分别利用GPU及CPU模拟，在下表中给出了各自每帧所需的平均时 间，其中最后一列为本发明在GPU上获得的加速效果。可以看到通过并行化处理，本发明获 得了最高接近16倍的加速比，提升了模拟的效率。
[0145]
[0146] 以上列举的仅是本发明的具体实施例。显然，本发明不限于以上实施例，还可以有 许多变形。本领域的普通技术人员能从本发明公开的内容直接导出或联想到的所有变形， 均应认为是本发明的保护范围。
【主权项】
1. 一种设及热传导与动态黏度的真实感流体仿真方法，其特征在于该方法包括W下步 骤： a) 基于光滑粒子流体动力学(SPH)模型模拟流体及热传导过程，具体包括： 0邻居查找 首先将所有流体及固体用离散粒子表示，并查找得到与每个粒子i距离在光滑核半径h 内的邻居粒子集合Ni; ii )各粒子导热系数计算 利用SPH算法的插值方式，考虑距离不同的邻居粒子对某个粒子i的影响，计算不同状 态下粒子i的导热系数ki，其具体公式为：其中化asei表示粒子i的当前状态，mj、Pj、kj分别表示粒子i的邻居粒子j的质量、密度、 导热系数，Wij为形如％=:W(fc^)的光滑核函数，它是关于粒子i与j的位置XI与xj之间距 离II义广刮II及光滑核半径h的函数;kfluid与ksolid分别为该种流体在液态及固态下实际的导 热系数； 化）各粒子溫度变化率计算 首先计算流体内部的热传导，对每个粒子i，其溫度随时间变化率计算公式为：其中T汲ki分别为粒子i的溫度及导热系数，A为拉普拉斯算子;同时Δ AiTi)的计算可 通过公?得到，其中17表示梯度； 其次计算流体与外界的热传导，对每个粒子i，其溫度随时间变化率计算公式为：其中Tb为固体粒子溫度，Ri可通过2计算得到，P0为流体的静态密度； 然后更新粒子i的溫度Τι，具体公式为：其中A t为时间步长，Ti(t)与Ti(t+Δ t)分别为粒子i在时亥Ijt与t+A t的溫度值； iv)相变处理 由于相变洽的存在，分别设置两个溫度阔值Tmelt和TsDlid ;当粒子溫度Ti大于Tmelt时，粒 子转化为液体粒子，当Ti小于TsDlid时，粒子转化为固体粒子；当Ti处于TsDlid与Tmelt之间时， 则将其看作为临界状态的流体粒子； b) 动态黏度计算的引入，具体包括： 在上述过程得到每个粒子i的最新时刻溫度Τι后，由于不同溫度对流体黏性系数存在一 定的影响，在SKI模型中引入动态黏度的计算，实现流体溫度变化过程中细节的逼真模拟， 其中溫度Τι与流体黏性系数μι的变化关系具体为： log(log(yi+丫））=q-y log(Ti) 其中丫为常数，通常取值在0.6到0.9之间，q与y则为与流体自身有关的特定参数，为定 值； c) PCISPH算法调用 利用PCISPH算法，对每个粒子i分别计算其所受到压力、黏性力W及外力而产生的加速 度曰严6、日严3及ai°theT，其计算公式分别为：其中，其中111、口、口、4、乂分别表示粒子的质量、压强、密度、黏性系数及速度，旨为重力加速 度;然后通过粒子i所受的加速度计算得到其在下一时刻的位置W模拟流体运动变化规律； d) GPU并行加速，具体包括： 利用通用并行计算架构CUDA，将所有粒子的物理属性传输至GPU，并对粒子进行排序， 加速邻居查找过程;对于每个粒子，为其开辟一个独立的线程W计算其密度、溫度、导热系 数、加速度，并更新流体粒子的速度与位置，从而并行地提升流体仿真的效率。
【文档编号】G06F17/50GK106096215SQ201610602958
【公开日】2016年11月9日
【申请日】2016年7月28日 公开号201610602958.7, CN 106096215 A, CN 106096215A, CN 201610602958, CN-A-106096215, CN106096215 A, CN106096215A, CN201610602958, CN201610602958.7
【发明人】王长波, 张申帆, 孔凡龙, 李晨
【申请人】华东师范大学