一种基于SPH方法近似求解的粘性流体建模方法技术

本发明专利技术是一种基于SPH方法近似求解的粘性流体建模方法。本方法首先根据粒子系统方法的特点对N‑S方程进行简化；其次根据简化后的方程，构建模拟粘性流体的框架；然后根据模拟框架基于SPH方法对各项物理量进行近似求解。通过压力‑密度法近似求解N‑S方程中对应的压力项，解决了SPH方法模拟粘性流体受严格的时间步长约束的问题；同时，解决了流体运动过程中粒子聚簇的现象。通过构建弹簧模型，即在邻近粒子对之间添加、修改和删除弹簧，实现流体运动的粘弹性；并通过引入温度对粒子相对速度的影响，从而描述温度对流体粘性的影响。本发明专利技术与传统方法相比，计算简单，基本上能够实现实时模拟，且能够更加逼真地模拟粘性流体运动场景。

【技术实现步骤摘要】
一种基于SPH方法近似求解的粘性流体建模方法
本专利技术涉及计算机图形学领域，具体涉及一种基于SPH方法近似求解的粘性流体建模方法。
技术介绍
粘性流体现象广泛存在于自然界、日常生活中、计算机游戏动画中和工业生产中。例如，蜂蜜、岩浆、泥土、油漆、蜡等。在计算机流体模拟的众多方法中，基于物理的流体模拟使用流体力学方程来描述水体内部速度、压力、密度等物理量随时间的变化，从而能够真实呈现粘性流体的运动。计算机模拟粘性流体的最普遍的方法是基于著名的纳维-斯托克方程(参考文件1：M.Griebel,T.Dornsheifer,T.Neunhoeffer.NumericalSimulationinFluidDynamics:APracticalIntroduction[M].SIAM:SocietyforIndustrialandAppliedMathematics,1997.)，简称N-S方程。但是N-S方程是一个非线性偏微分方程，在求解过程中需要大量地计算节点，导致计算复杂度非常高。同时，为了更逼真地模拟粘性流体，还需要在N-S方程中加入弹性应力项或弹簧模型，从而导致计算复杂度更高。因此采用传统的N-S方程来对粘性流体建模，从计算机应用的角度出发是不可行的。但是可以通过对N-S方程简化近似求解，减少计算的时间复杂度。光滑粒子动态学方法(SmoothParticleHydrodynamics)，简称SPH方法，是目前对N-S方程近似简化比较成熟和流行的方案。SPH是一种用于求解连续介质动力学方程的无网格粒子方法，其通过引入空间场函数和核函数的概念，将流体控制方程离散化。流体的问题域被离散化成有限的粒子，每个粒子带有自身材料性质和其他物理属性，如：速度、压力、密度等。每个粒子的每一时刻物理属性可以由上一时刻的邻近粒子加权求得。通过粒子物理属性的动态更新，就可以模拟出复杂的粘性流体的运动现象，如：冰块熔化、岩浆流动等。Stora等人(参考文件2：StoraD,AgliatiPO,CaniMP.Animatinglavaflows[C].ProceedingsoftheGraphicsInterface,1999:203-210.)采用SPH方法模拟火山熔岩粘性运动的过程。其在状态方程添加温度项并计算热传递量来控制流体的粘度，用数字高程模型模拟地形，减少了计算粒子投影到水平面的高度计算。由于作者只考虑了流体内部的作用力以及流体与不同地形之间的作用力，所以并没有模拟出与不规则岩石碰撞时流体折叠翻卷等细节表现。Steele等人(参考文件3：SteeleK,ClineD,EgbertPK,DinersteinJ.Modelingandrenderingviscousliquids[J].ComputerAnimationandVirtualWorlds,2004,15(3-4):183-192.)基于SPH方法，利用粘附矩阵来表示流体与物体间的作用力，同时用粒子表示物体表面，较快捷地描绘出粘性流体与固体交互。由于该方法在每个时间步长内都要进行重复的自适应粒子位置更新来维持体积守恒，从而导致计算量较大，耗损时间严重。常元章等人(参考文件4：ChangYuan-zhang,BaoKai,LiuYou-quan,etal.AUnifiedParticle-BasedMethodforNewtonianandViscoelasticFluidsAnimation[J].ChineseJournalofComputers,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为本文档来自技高网...

【技术保护点】
1.一种基于SPH方法近似求解的粘性流体建模方法，其特征在于，实现步骤如下：步骤一：根据SPH方法的特点，简化N‑S方程构建模拟粘性流体基本框架；首先，根据SPH方法粒子系统的特点，简化N‑S方程：

【技术特征摘要】
1.一种基于SPH方法近似求解的粘性流体建模方法，其特征在于，实现步骤如下：步骤一：根据SPH方法的特点，简化N-S方程构建模拟粘性流体基本框架；首先，根据SPH方法粒子系统的特点，简化N-S方程：其中ρ是密度，t是时间，p是压力，a是外力密度场矢量，即ρa是外力矢量，μ是粘度参数。其次，基于SPH方法，构建物理量基本求解框架。然后，根据简化N-S方程，对方程中对应项基于SPH方法近似求解，构建模拟粘性流体的框架。即，外力项ρa一般情况只包括重力，直接显式给出；模型压力项根据压力-密度松弛法对粒子位置的影响近似求解得到；粘性项通过构建弹簧模型以及对粒子速度根据温度进行矫正来近似求解得到。步骤二：基于步骤一构建的粘性流体模拟框架求解粒子物理属性模型。步骤2.1，计算密度与压力，求解压力-密度松弛。基于SPH计算框架，求解粒子i的密度；根据理想气体状态方程改进求得粒子i的压强Pij，进而在粒子i到粒子j位置矢量方向上求取压力-密度松弛值Dij。该步骤整体求解过程作如下形式计算：其中rij为粒子i和j之间的距离，h为光滑核函数的支持域，ρ0为静止密度常数，k为一常数值，为粒子i到j方向上的单位向量。步骤2.2，计算邻近密度与邻近压力，改进压力-密度松弛。当粒子i密度ρi达到静止密度ρ0时，则粒子i将聚集少数粒子形成聚簇现象，导致无法真实的模拟流体。故本发明引入邻近密度和邻近压力Pin，作如下形式计算：其中kn为一常数值。根据求解的邻近密度和邻近压力值，对步骤2.1求解的压力-密度值作如下形式矫正：对于粒子i的位置和其邻近粒子j作如下松弛计算：xj←xj+Dij/2xi←xi-Dij/2步骤2.3，构建弹簧模型，求解弹簧松弛。通过构建弹簧模型，实现粘性流体的粘弹性效果。在邻近粒子对之间添加、修改和删除弹簧Lij，构建弹簧模型，并依据粒子间的弹簧长度计算出弹簧松弛量以修正粒子i和j位置。由于互为邻近粒子的粒子对之间只需经过本步骤一次计算，故要求计算过程中粒子编号i＜j。如果邻近粒子对之间，没有弹簧，则增加一个初始值为Lij＝h(h为粒子系统中光滑核函数的支持域常数)的弹簧；根据特定材料设定其对应容忍变化的长度d＝γLij，其中γ为容忍变化率(一般在0-0.2范围取值)。当粒子间的距离rij大于L+d(拉伸状态)，则Lij作如下修改：Lij＝Lij+Δtα(rij-L-d)当粒子间的距离rij小于L-d时(压缩状态)，则Lij作如下修改：Lij＝Lij-Δtα(L-d-rij)则粒子对之间弹簧松弛量可以作如下计算：其中α为塑性系数，ks是弹簧模型的一个常数系数。然后对粒子i和j的位置作如下松弛计算：步骤2.4，求解热方程，计算温度变化率dT/dt。本...

【专利技术属性】
技术研发人员：艾明晶，杜安定，陈保贺，
申请(专利权)人：北京航空航天大学，
类型：发明
国别省市：北京,11