本发明涉及计算机图形学领域，具体涉及一种基于SPH方法近似求解的粘性流体建模方法。
背景技术：
粘性流体现象广泛存在于自然界、日常生活中、计算机游戏动画中和工业生产中。例如，蜂蜜、岩浆、泥土、油漆、蜡等。在计算机流体模拟的众多方法中，基于物理的流体模拟使用流体力学方程来描述水体内部速度、压力、密度等物理量随时间的变化，从而能够真实呈现粘性流体的运动。
计算机模拟粘性流体的最普遍的方法是基于著名的纳维-斯托克方程(参考文件1：M. Griebel,T.Dornsheifer,T.Neunhoeffer.Numerical Simulation in Fluid Dynamics:A Practical Introduction[M].SIAM:Society for Industrial and Applied Mathematics,1997.)，简称N-S方程。但是N-S方程是一个非线性偏微分方程，在求解过程中需要大量地计算节点，导致计算复杂度非常高。同时，为了更逼真地模拟粘性流体，还需要在N-S方程中加入弹性应力项或弹簧模型，从而导致计算复杂度更高。因此采用传统的N-S方程来对粘性流体建模，从计算机应用的角度出发是不可行的。但是可以通过对N-S方程简化近似求解，减少计算的时间复杂度。光滑粒子动态学方法(Smooth Particle Hydrodynamics)，简称SPH方法，是目前对N-S方程近似简化比较成熟和流行的方案。SPH是一种用于求解连续介质动力学方程的无网格粒子方法，其通过引入空间场函数和核函数的概念，将流体控制方程离散化。流体的问题域被离散化成有限的粒子，每个粒子带有自身材料性质和其他物理属性，如：速度、压力、密度等。每个粒子的每一时刻物理属性可以由上一时刻的邻近粒子加权求得。通过粒子物理属性的动态更新，就可以模拟出复杂的粘性流体的运动现象，如：冰块熔化、岩浆流动等。
Stora等人(参考文件2：Stora D,Agliati P O,Cani M P.Animating lava flows[C].Proceedings of the Graphics Interface,1999:203-210.)采用SPH方法模拟火山熔岩粘性运动的过程。其在状态方程添加温度项并计算热传递量来控制流体的粘度，用数字高程模型模拟地形，减少了计算粒子投影到水平面的高度计算。由于作者只考虑了流体内部的作用力以及流体与不同地形之间的作用力，所以并没有模拟出与不规则岩石碰撞时流体折叠翻卷等细节表现。Steele等人(参考文件3：Steele K,Cline D,Egbert P K,Dinerstein J.Modeling and rendering viscous liquids[J].Computer Animation and Virtual Worlds,2004,15(3-4):183-192.)基于SPH方法，利用粘附矩阵来表示流体与物体间的作用力，同时用粒子表示物体表面，较快捷地描绘出粘性流体与固体交互。由于该方法在每个时间步长内都要进行重复的自适应粒子位置更新来维持体积守恒，从而导致计算量较大，耗损时间严重。常元章等人(参考文件4：Chang Yuan-zhang, Bao Kai,Liu You-quan,et al.A Unified Particle-Based Method for Newtonian and Viscoelastic Fluids Animation[J].Chinese Journal of Computers,2010,33(7):1286-1294.)通过基于SPH方法在传统的N-S方程中引入弹性应力项，统一模拟了牛顿流体和粘弹性流体。此方法有效且可控性强，通过调节参数可以实现不同粘弹性的流体。但是因为求解粘弹性方程过程复杂，导致模拟速度不够快，实时性不佳。
上述基于传统的SPH方法，都是在SPH框架下进行数值插值来计算各种力和其他物理属性，此种计算模式需要严格的时间步长和大量的计算时间。其次，由于方程的稳定性可能会导致粒子发生聚簇现象，导致无法真实地模拟粘性流体的细节。
技术实现要素：
本发明的目的是为了克服上述现有基于SPH方法无法逼真模拟粘性流体，具有严格的时间步长约束，以及计算中高数值耗散带来的计算速度慢等缺点。本发明提出一种基于SPH方法近似求解的粘性流体建模方法，旨在提高粘性流体的建模过程计算效率和粘性流体模拟效果的逼真度。
本发明提出的一种基于SPH方法近似求解的粘性流体建模方法，具体包括如下步骤：
步骤一：根据SPH方法的特点，简化N-S方程构建模拟粘性流体基本框架。
SPH方法是粒子系统的一种插值方法，场量被离散化定义到各个处于空间具体位置的粒子上，并可以在空间的任何位置被计算。在SPH方法中，每个粒子的物理属性是通过邻近粒子以对称光滑核函数为权的物理属性累加求和得到。根据SPH方法，一个标量A插值计算的表达式如式(1)所示。
其中j表示粒子编号，mj表示粒子的质量，r表示粒子的位置矢量，W(r,h)表示支持域为h的光滑核函数。W(r,h)是偶函数，满足正则化条件等性质。
而A的梯度和拉普拉斯算子可以表示为式(2)和式(3)。
N-S方程是描述粘性不可压缩流体质量守恒和动量守恒的运动方程。其分别由式(4)和式 (5)表示：
其中v是速度矢量，ρ是密度，t是时间，p是压力，a是外力的密度场矢量，即ρa表示外力矢量，μ是粘度参数。
由于在基于SPH方法的粒子系统中，粒子的总数不变且每个粒子的质量相同，故式(4) 所描述的质量守恒是满足的，即式(4)可以被忽略。由于流体的移动是由粒子随时间的变化来表示，所以对于粒子系统而言，速度场的导数实质上可以简单表示为Dv/Dt，至于对流项粒子系统是可以省略的。故式(5)可以简化为式(6)。
力场的总和f可由式(7)表示。
由式(6)和式(7)可得每个粒子i的加速度gi为：
其中vi表示粒子i的速度，fi表示粒子i的力，ρi表示粒子i的密度。
由式(6)可知，外力ρa、模型压力和粘度决定粒子动量ρ(Dv/Dt)的变化，故此三个物理量是粘性流体模拟框架的主要求解对象。传统SPH方法构建粘性流体，根据式(8) 通过描述力密度场，构建统一的基于SPH方法求解物理模型的框架。即外力fe，模型压力fp和粘力fv都统一由SPH方法求解，且根据牛顿第二定律描述粒子的运动状态。这种简单显式对粒子各种力的数值求解，需要严格的时间步长约束和大量的计算时间。fp对粒子运动状态的影响，通过压力-密度松弛方法计算得到；而fv对粒子运动的影响，即粒子运动表现出粘弹性，通过在粒子间添加、修改和删除弹簧质子体现。
上述粒子系统的模拟粘性流体的框架中，根据粒子系统的特点，满足质量守恒和省略对流项的计算，极大简化了N-S方程，将有效提高粘性流体模拟的计算效率。同时，采用压力- 密度松弛方法和弹簧质子可以避免时间步长的约束和粒子聚簇现象，并能更加真实模拟粘性流体。
步骤二：基于步骤一构建的粘性流体模拟框架求解粒子物理属性模型。
步骤2.1，计算密度与压力，求解压力-密度松弛。由步骤一中式(1)可得，粒子的密度计算公式如式(9)所示。
由于考虑到所有粒子有相同的质量，故省略质量属性，则得式(10)。
其中Wq(r-rj,h)＝(1-rij/h)2为二次光滑核函数。
本发明定义伪压力Pi(由于省略质量属性)，其与当前粒子密度ρi和静止密度ρ0的差成正比，比例系数为k，其中k为刚度参数的常数值，Pi如式(11)所示。
Pi＝k(ρi-ρ0) (11)
则粒子之间的压力-密度松弛值Dij可由式(12)计算得到。
其中为粒子i到粒子j的单位法向量，Wl(r-rj,h)＝(1-rij/h)为线性光滑核函数。
Dij用于调整粒子j在时间步长△t内到达的预测位置。同时，根据牛顿第三定律，粒子i同样需要一个数值相等、方向相反的Dij来调整预测位置。
步骤2.2，计算邻近密度与邻近压力，改进压力-密度松弛。由式(11)可知，当粒子的密度ρi达到静止密度ρ0时，则粒子将聚集少数粒子形成聚簇现象，导致无法真实地模拟流体。为了解决聚簇现象，本发明定义依赖于不同的光滑核函数计算得到的新密度，即邻近密度如式(13)所示。
其中Wc(r-rj,h)＝(1-rij/h)3为三次光滑核函数。
与邻近密度对应的邻近压力则可以将式(12)求解的Dij改进为，如式(14) 所示。
步骤2.3，构建弹簧模型，求解弹簧松弛。流体的粘弹性主要由弹性、塑性和粘性三部分构成。其中粘性通过基于传统的SPH方法框架计算得到。弹性与塑性则通过在邻近粒子间增加、修改或删除弹簧模拟。塑性表示材料在外力作用下，材料发生永久变形而不被破坏其完整性的能力。在本发明提出的弹簧模型中可定义塑性为粒子间弹簧长度的变化并修正粒子间的弹簧长度。可由式(15)求得。
其中△t为时间步长，α为塑性系数，r为粒子间的距离，L为静止弹簧长度。
然后，定义材料容忍变化长度d＝γLij，其中γ为容忍变化率(一般在0-0.2范围取值)，Lij为粒子间的弹簧长度。当粒子间的距离rij大于L+d时(拉伸状态)，则如式(16)所示修改Lij； rij小于L-d时(压缩状态)，则如式(17)所示修改Lij。
Lij＝Lij+△tα(rij-L-d) (14)
Lij＝Lij-△tα(L-d-rij) (15)
其中弹性通过在邻近粒子对之间添加一个弹簧，当邻近粒子对中粒子发生位移时，粒子位置通过弹簧松弛量修正，以表现弹性。其中由式(18)计算求得。
其中ks是弹簧模型的一个常数系数，Lij是粒子对ij之间的弹簧长度，rij表示粒子对之间的距离，表示粒子i到粒子j的单位向量。则对每对邻近粒子对的粒子位置进行弹簧松弛修正，如式(19)所示。
其中xi和xj分别表示粒子i和j的当前位置。
步骤2.4，求解热方程，计算温度变化率dT/dt。本发明为了模拟材料熔化过程，丰富粘性流体运动场景，引入温度变化对流体粘性的影响。若不模拟材料熔化场景时，不需进行该步骤。温度随时间的变化遵守热方程式(其中λ为热扩散系数)规律。由步骤一中式(3)可计算温度拉普拉斯算子如式(20)所示。
步骤2.5，对粒子速度进行修正，描述流体粘性。流体粘性具有平滑速度场的效果。故可在邻近粒子对之间，引入速度松弛量I，修正粒子的速度，描述流体粘性。粒子之间的粘性与粒子对之间的相对速度相关，故定义粒子间相对径向速度速度松弛量I可由式(21)求得。
其中△t为时间步长，β为与温度相关的参数，由式(23)求得。u为相对径向速度，r为粒子间距离，h为支持域，表示粒子i到粒子j的单位向量。
则根据式(22)对粒子i和j的速度进行修正。
其中βmin和βmax分别是β取值的最小值和最大值，其取值与特定的液体材料属性相关； Tmin和Tmax分别是对应温度T的最小值与最大值。由式(22)可知当T越大时，β越小，速度松弛量I越小，描述的粘性就越小。
步骤三：根据利用步骤二所求解的物理模型，构建粘性流体运动场景，输出结果并显示。步骤二中求解的粘性流体模型，仅仅是基于数学物理求解得出的结果，从视觉上描述，是由一系列粒子构成，并不能逼真地呈现现实中的粘性流体场景。因此还需要对上述求解模型，进行相应的渲染操作，才能真实呈现粘性流体模型。本发明中采用OpenGL提供的库函数进行编程，对粘性流体进行渲染，从而展示逼真的流体效果。
本发明方法的优点和积极效果在于：本发明利用粒子系统的特点有效地简化N-S方程，构建基于SPH方法框架求解粘性流体物理模型参数，提高计算效率；然后，用压力-密度松弛方法对粒子位置修正，以代替传统SPH方法显式求解压力项对粒子运动的影响，减少了计算量，解决了粒子聚簇现象；同时，通过弹簧模型描述邻近粒子对之间的弹性和塑性，引入粘性相关速度松弛量矫正粒子速度，描述流体粘性，从而多维度地描述了流体的粘弹性，更加逼真地模拟粘性流体运动场景；最后，引入热力方程，描述了温度对材料的影响，实现了材料熔化的场景，丰富了粘性流体的运动场景。
附图说明
图1是SPH光滑核函数粒子近似求解原理图；
图2是本发明粒子空间位移变化原理简图；
图3是本发明一种基于SPH方法近似求解的粘性流体建模方法的整体步骤流程图；
图4是本发明粘性材料受巨大吸引力粒子运动状态模拟序列图；
图5是本发明冰块熔化过程中粘性流体运动状态模拟序列图；
图6是本发明大型粘性液滴各种粒子运动状态和以及表面渲染后状态图。
具体实施方式
下面将结合附图和实施例对本发明的技术方案作进一步的详细说明。
本发明提出一种基于SPH方法近似求解的粘性流体建模方法，利用粒子系统流体建模的特点极大简化N-S方程，构建粘性流体模拟的基本框架；并通过压力-密度松弛模型，修正粒子的运动位置，解决了采用传统SPH方法产生的粒子聚簇现象的问题；然后，引入弹簧模型描述流体弹性与塑性，对粒子对速度进行修正，并引入温度因素调节修正值大小，真实地描述流体粘性，从而多维度地描述了流体的粘弹性；此外，引入热力方程，描述材料熔化场景，丰富粘性流体运动场景。
本发明提出一种基于SPH方法近似求解的粘性流体建模方法，旨在提高传统SPH方法模拟粘性流体建模的性能，下面具体对实现步骤进行说明。
步骤一：根据SPH方法的特点，简化N-S方程构建模拟粘性流体基本框架。
步骤1.1，根据SPH方法的特点，简化N-S方程。N-S方程可由式(4)和式(5)表示，由于在基于SPH方法的粒子系统中，粒子的总数不变且每个粒子的质量相同，故式(4)所描述的质量守恒是满足的，即式(4)可忽略不计。由于流体的移动是由粒子位置随时间的变化来表示，所以对于粒子系统而言，速度场的导数实质上可以简单表示为Dv/Dt，至于对流项粒子系统是可以省略的。故根据N-S方程描述粘性流体的运动状态，最终可简化近似，根据式 (6)描述粘性流体的运动状态。
步骤1.2，基于SPH方法，构建基本求解框架。SPH方法是粒子系统的一种插值方法，物理量被离散化定义到各个处于空间具体位置的粒子上，并可以在空间的任何位置被计算。在SPH方法中，每个粒子的物理属性是通过邻近粒子以对称光滑核函数为权的物理属性累加求和近似得到。其中光滑核函数W(r-rj,h)是偶函数，满足正则化条件等性质。且邻近粒子是指处在光滑核函数域内的粒子。图1“SPH光滑核函数粒子近似求解原理图”展示了粒子i的物理属性是由支持域h内的粒子j的物理属性以光滑核函数为权，累加求和所得。
故由式(1)、式(2)或式(3)可以计算各种类型的物理量。由于假定粒子系统中粒子的质量是相同的，故本发明在计算过程中省略粒子的质量。且不同的物理量，需要根据不同的光滑核函数计算。
步骤1.3，改进传统SPH方法描述粘性流体运动的方式，确定粘性流体模拟框架。由式 (6)可知外力ρa、模型压力和粘度决定粒子动量ρ(Dv/Dt)的变化，故此三个物理模型是粘性流体模拟框架的主要求解对象。传统SPH方法构建粘性流体，根据式(8)通过描述力密度场，构建统一的基于SPH方法求解物理模型的框架。即外力fe、模型压力fp和粘力 fv都统一由SPH方法求解，且根据牛顿第二定律表现粒子的运动状态。这种简单显式对粒子各种力的数值插值求解，需要严格的时间步长约束和大量的计算时间。本发明提出的粒子系统中，外力fe简单包括重力和其他可以显式地通过对应的加速度值给出的外力(如吸引力或拉力等)；然后模型压力fp对粒子运动的影响，通过压力-密度松弛调节粒子位移来描述；而fv对粒子运动的影响，即粒子运动表现出粘弹性，通过弹簧模型来描述。图2“粒子空间位移变化原理简图”展示了本发明提出的粘性流体模拟方法中粒子运动的原理。首先，粒子i 在外力作用下运动到位置其次，经过压力-密度松弛Dij调节粒子的位置；然后，经过弹簧模型计算弹簧松弛量调节粒子的位置。故本发明提出模拟粘性流体的方法步骤，如图3 “一种基于SPH方法近似求解的粘性流体建模方法的整体步骤流程图”所示。
步骤二：基于步骤一构建的粘性流体模拟框架求解粒子物理属性模型。
步骤2.1，计算密度与压力，求解压力-密度松弛。粒子的密度ρ可以根据SPH方法求得，故由式(10)可得式(24)，来实际计算粒子i密度ρi。
根据式(11)计算出粒子的压力Pi，根据式(12)计算出压力-密度松弛值Dij。
步骤2.2，计算邻近密度与邻近压力，改进压力-密度松弛。由式(11)可知，当粒子的密度ρi达到静止密度ρ0时，则粒子将聚集少数粒子形成聚簇现象，导致无法真实地模拟流体。为了解决聚簇现象，本发明定义依赖于不同的光滑核函数计算得到的新密度，即邻近密度通过邻近密度计算对应的邻近压力最后利用邻近压力修正压力-密度松弛值Dij，如式(14)所示。Dij用于调整粒子j在时间步长△t内到达的预测位置。同时，根据牛顿第三定律，粒子i同样需要一个数值相等、方向相反的Dij去调整预测位置。如式(25)所示。
步骤2.3，构建弹簧模型，求解弹簧松弛。本发明为了实现粘性流体的粘弹性效果，引入弹簧模型。通过在邻近粒子对之间添加、修改和删除弹簧Lij构建弹簧模型，并依据粒子间的弹簧长度计算出弹簧松弛量由去修正粒子的位置，实现粘性效果，如式(19)所示。
由于互为邻近粒子的粒子对之间只需经过本步骤一次计算，故要求计算过程中粒子编号 i<j。如果邻近粒子对之间，没有弹簧，则增加一个初始值为Lij＝h(h为粒子系统中光滑函数的支持域)的弹簧；根据特定材料设定其对应容忍变化的长度d＝γLij，其中γ为容忍变化率(一般在0-0.2范围取值)。当粒子间的距离rij大于L+d时(拉伸状态)，则如式(16)所示修改Lij；当rij小于L-d时(压缩状态)，则如式(17)所示修改Lij。
步骤2.4，求解热方程，计算温度变化率dT/dt。本发明为了模拟材料熔化过程，丰富粘性流体运动场景，引入温度的变化对流体粘性的影响。但由于本步骤计算量较大，若不模拟熔化场景时，可将本步骤省略，直接给出每个粒子恒定的温度。温度随时间的变化遵守热方程式(其中λ为热扩散系数)规律。根据SPH方法，如式(20)所示，求解则根据热方程可求解粒子i下一时刻的温度Ti，步骤2.5中将引入与温度相关的参数，调节粒子之间的速度，从而描述温度对流体粘性的影响。
步骤2.5，对粒子速度进行修正，描述流体粘性。首先，经过上述步骤对粒子的位置进行调节后，根据计算粒子此时的速度。粒子之间的粘性与粒子对之间的相对速度有关，故引入相对径向速度并由式(21)求出速度松弛量I。其中参数β由式(23)求出，式(23)中的温度T可以由粒子i和j的温度平均值求出。然后根据式(22)对粒子i和 j的速度进行修正。由于粒子i和j互为邻近粒子，为避免重复对速度修正计算，而导致抵消速度修正最终结果，故要求在计算过程中粒子编号i<j。
步骤三：根据利用步骤二所求解的物理模型，构建粘性流体运动场景，输出结果并显示。步骤二中求解的粘性流体模型，仅仅是基于数学物理求解得出的结果，从视觉上描述，是由一系列粒子构成，并不能逼真地呈现现实中的粘性流体场景。因此还需要对上述求解模型，进行相应的渲染操作，才能真实呈现粘性流体模型。本发明中采用OpenGL提供的库函数进行编程，对粘性流体进行渲染，从而展示逼真的流体效果。
本发明所有测试都是在戴尔台式机上进行的，CPU型号是英特尔Core i7-4790@3.6GHz，显卡型号AMD RadeonTM R5 240。
图4展示了粘性流体材料在巨大吸引力的场景下的运动状态的实验效果，特别突出地展示了粘性流体运动时的粘弹性效果。其中，图序列从上到下展示了粘性流体材料在右端上方受到巨大外部吸引力情况下粒子的运动状态。
图5展示了冰块熔化的模拟过程，特别突出地展示了材料熔化时，粘性流体的运动状态场景。其中，(a)是冰块的初始状态，(b)、(c)和(d)是冰块随时间变化而发生熔化过程中的序列图。
图6中上面一行的各子图展示了大型粘性液滴在不受外力(最左边的子图)和受各种外力条件下的不同粒子运动状态，下面一行的各子图展示了对上面的子图进行表面渲染后的相应运动状态。