CN104318599A - 一种基于几何特征的高精度流体动画建模方法 - Google Patents

Abstract

本发明公开了一种基于几何特征的高精度流体动画建模方法，包括以下步骤：对流体物理模型中的初始条件和边界条件进行初始化；从流体的物理模型中获取几何特征计算所需要的中间速度场、当前时间步中的流体速度场和符号距离场函数参数；构建流体动画的几何模型，对流体动画中的细节敏感区域进行自动化定位、量化和跟踪；根据几何特征的特征函数，控制物理模型的粒子集，耦合几何模型和物理模型；提取表示流体表面和流体水花的三角形网格表面，并对其进行真实感渲染，生成一帧流体动画；本发明能有效地产生高精度的流体水面和水花等细节效果，实现高精度流体动画的自动化生成，具有精度好、效率高的特点。

Description

一种基于几何特征的高精度流体动画建模方法

技术领域

本发明涉及一种基于几何特征的高精度流体动画建模方法。

背景技术

流体动画在虚拟现实、动漫特效、互动游戏、建筑设计、产品展示、医学、军事模拟仿真等方面有着广泛的应用。在这些流体动画中，视觉逼真性是最基本的追求。然而，目前很多流体动画都是动画师在工具软件中手工制作。由于流体运动非常复杂，即使是具有丰富经验的动画师，也很难随心所欲的设计想要的动画。此外流体动画制作速度很慢、成本高，是典型的“劳动密集型”产业。因此，如何自动生成具有丰富细节的流体动画是一项很有挑战性的问题。传统的流体动画方法存在计算代价高、效率低且流体细节容易丢失的问题，难以生成视觉逼真的高精度流体动画。

发明内容

本发明为了解决上述问题，提出了一种基于几何特征的高精度流体动画建模方法，本方法该方法通过对视觉逼真的高精度流体动画自动化生成，有效提高其在计算机图形学、计算机动画、特效制作、游戏娱乐等应用领域中的处理效率，为进一步推动动漫行业从“劳动密集型”到“技术密集型”转型，拓展其应用领域和范围奠定理论基础。

为了实现上述目的，本发明采用如下技术方案：

一种基于几何特征的高精度流体动画建模方法，包括以下步骤:

(1)对流体物理模型中的初始条件和边界条件进行初始化；

(2)从流体的物理模型中获取几何特征计算所需要的中间速度场、当前时间步中的流体速度场和符号距离场函数参数；

(3)构建流体动画的几何模型，对流体动画中的细节敏感区域进行自动化定位、量化和跟踪；

(4)根据几何特征的特征函数，控制物理模型的粒子集，耦合几何模型和物理模型；

(5)提取表示流体表面和流体水花的三角形网格表面，并对其进行真实感渲染，生成一帧流体动画；

(6)查看是否完成全部流体动画的制作，如果未完成，则返回步骤(2)，继续进行制作。

所述步骤(1)中，物理模型初始化的步骤包括：

(1-1)初始化背景网格中的流体速度场u；

(1-2)初始化表示流体表面的符号距离场函数φ；

(1-3)初始化对流体表面进行检测校正的粒子集P_c和表示流体水花的粒子集P_s；

(1-4)初始化边界，即使用标记位表示流体计算域中的每个网格单元的类型；为此需要将固体障碍物按计算域网格单元尺寸离散化，离散化后被固体占据的网格单元标记位设置为1，否则设置为0。

所述步骤(2)中，从流体的物理模型中获取几何特征的具体方法如下：

(2-1)根据初始化背景网格中的流体速度场u，粒子集P_s中的速度计算背景网格中的中间速度场；

(2-2)根据中间速度场，通过求解耦合欧拉和拉格朗日两种模型的压强泊松方程，获取当前时间步中的流体速度场uⁿ⁺¹；

(2-3)根据所获得的当前时间步流体速度场，计算更新后的符号距离场函数φⁿ⁺¹。

所述步骤(3)中，建立流体动画的几何模型，对流体动画中的细节敏感区域进行自动化定位、量化和跟踪的基本步骤如下：

(3-1)构建封闭流形，由封闭流体界面或者非封闭的流体界面与计算域边界形成，封闭流形边界及其周围3层网格称为窄带区域；

(3-2)定义局部中轴(Local Medial Axis,LMA)为封闭流形中位于窄带Ω_b上的内接圆圆心的轨迹；

(3-3)根据Step2中获得的符号距离场函数φⁿ⁺¹，使用散度值判断法计算LMA；

(3-4)根据局部中轴，定义定位量化函数(Location and Quantization Function,LF)，对于流形边界上的网格点，其函数值为给定阈值λ与该点到局部中轴距离值的较小值；窄带区域上其他网格点的定位量化函数值为该点到流形边界的最短距离加上流形边界对应点的定位量化函数值；对于计算域上其他网格点的定位量化函数值将其设置为一个大常数，定位量化函数能够实现对大曲率、细长和边界细节敏感区域的光滑定位和量化；

(3-5)根据定位量化函数的定义，通过使用KD检索的方法，计算流体表面上网格点的定位量化函数值；

(3-6)根据定位量化函数的定义，通过使用快速行进算法对流形边界点的定位量化函数值外推，计算窄带域上其他网格点的定位量化函数值；

(3-7)根据定位量化函数的定义，设定背景网格中窄带域以外的定位量化函数值为大常数，指定为N_max×Δh，其中N_max为计算域在x、y、z三个轴向上网格个数的最大值，Δh为网格单元的尺寸。

(3-8)在每个时间步中累积计算流体表面的变形程度，如果变形程度小于设定阈值，则需要对定位量化函数进行对流更新；否则需要对定位量化函数进行重新初始化。

所述步骤(4)中，其具体方法为：

(4-1)将当前时间步n的定位量化函数值LFⁿ映射为粒子的分布密度pⁿ(x)；

(4-2)根据pⁿ(x)计算需要增/删的粒子数

(4-3)在网格单元x中增加或删除个粒子，得到更新后的粒子集和

(4-4)更新中粒子的速度和位置，然后利用它们检测和校正的误差，得符号距离场函数φⁿ⁺¹；

(4-5)更新中粒子的速度和位置，得下一个时间步n+1的水花粒子集

所述步骤(5)中，具体方法包括：

(5-1)使用移动立方体算法，从表示流体表面的符号距离场函数φⁿ⁺¹中提取零等值面，获得流体表面的三角形网格表面；

(5-2)使用构建粒子局部水平集算法，从粒子集中获取表示水花粒子的符号距离场函数然后再使用移动立方体算法，获得表示水花的零等值面。

本发明的有益效果为：

(1)采用基于几何特征的方法实现对流体动画的高精度建模，通过获取几何参数、建立几何模型、几何模型与物理模型耦合等过程，充分挖掘并使用流体细节之间的相关性以及流体表面不同区域的差异性等几何特征，实现高精度流体动画的自动化生成；

(2)利用定位量化函数能够有效的反映水面运动湍急、水花飞溅等细节产生的位置及细节的多少；

(3)该方法具有精度高、高效性的特点。

附图说明

图1为本发明的流程示意图。

具体实施方式：

下面结合附图与实施例对本发明作进一步说明。

如图1所示，基于几何特征的高精度流体动画建模方法主要包括以下几个过程：

过程1：初始化；

第一步：初始化液体表面。令Ω表示计算域，计算域网格记为G。令符号距离场函数φ(x,t)表示液体表面，其函数值为t时刻网格中心x到液体表面的符号距离。显然，在液体表面上，满足φ(x,t)＝0，记为液体内部满足φ(x,t)＜0，记为Ω^-；在液体外部满足φ(x,t)＞0，记为Ω⁺。

第二步：初始化流体速度场u。指定计算域网格G上的三个三维数值分别代表速度场u的x、y、z分量，在液体外部Ω⁺区域上，三个分量的值均为0；而在液体表面及内部区域上根据具体算例指定为非零值。

第三步：初始化检测校正粒子集P_c。在流体表面内部及外部三层网格上初始化检测校正粒子集P_c。对于P_c中的每个粒子：其位置是在内外三层网格内随机生成，其半径为0.1-0.5个网格单元的大小，其速度是根据粒子所在位置从初始化的速度场u中插值得到。

第四步：初始化表示水花的粒子集P_s为空。

第五步：初始化边界条件。将流体运动所在的环境模型利用体素化方法转化为方程计算时可识别的边界。它包含两个部分：固体与液体边界、空气和液体边界。本文使用flag标记位标识。flag＝0表示被固体占据的单元网格；flag＝1表示液体单元；flag＝2表示空气单元。

过程2：根据当前流体速度场u和粒子集P_s获取计算域网格中的中间速度场u^*。

第一步：令当前时间步为n，计算域网格G中的流体速度场为uⁿ，SPH粒子集P_s中的粒子速度为为了计算动态流体速度场uⁿ⁺¹，需要求解流体运动的NS方程

$\frac{1}{\mathrm{\ρ}}\frac{\mathrm{D\ρ}}{\mathrm{Dt}}+\▿\·u=0---\left(1\right)$

$\frac{\∂u}{\∂t}=-u\·(\▿u)-\frac{1}{\mathrm{\ρ}}\▿p+f---\left(2\right)$

其中，ρ为密度，u为流体速度,p为压强，f为外力，▽·为散度，▽为梯度。

第二步：对于区域上的每个网格单元，其中间速度场u^*可通过执行对流和添加外力两个步骤得到中间速度场u^*。其中，对流步可以使用半拉格朗日法求解，即添加外力步可表示为u^*＝u(x)+Δtf(x,t)，Δt为时间步长。

第三步：对于Ω⁺区域需计算P_s粒子集所在的每个网格单元中的u^*。为此，首先添加外力即然后再将其投影到网格G中，即其中，ω_p(x)为粒子对网格的贡献大小

${\mathrm{\ω}}_p\left(x\right)=\left\{\begin{array}{cc}c(1-{\left|\right|x-x_p\left|\right|}^2/r_p^2)& \mathrm{if}{\left|\right|x-x_p\left|\right|}^2\≤r_p^2\\ 0& \mathrm{else}\end{array}\right.$

这里，c为控制常数,x_p为粒子位置，r_p为粒子半径。

过程3：根据中间速度场，通过求解耦合欧拉和拉格朗日两种模型的压强泊松方程，获取计算域中的流体速度场uⁿ⁺¹。

第一步：变换方程(1)(2)，得到压强泊松方程。根据过程2方法获得u^*后，方程(2)右侧只剩下-1/ρ▽p，利用有限差分展开，得

uⁿ⁺¹＝u^*-▽pΔt/ρ   (3)

对方程(3)两侧同时取散度，可得压强泊松方程

▽·(▽pΔt/ρ)＝▽·u^*-▽·uⁿ⁺¹   (4)

第二步：计算公式(4)的右侧数值。在公式(4)中，▽·uⁿ⁺¹在计算域的不同位置取值不同。其中，区域属于液体内部，具有严格的不可压缩性，故▽·uⁿ⁺¹＝0。而在Ω⁺区域，液态粒子与空气混合，具有一定可压缩性，▽·uⁿ⁺¹≠0，其值可以根据网格单元中P_s粒子个数与目标粒子数估计。根据质量守恒方程(1)，ρ可以表示为ρ＝N_pm_p/V_p，其中，Ν_p为网格单元内的P_s粒子数目，m_p为粒子的质量，V_p为粒子所占的体积。故方程(1)可变换为

$\▿\·u=-\frac{1}{N_p}\frac{{\mathrm{DN}}_p}{\mathrm{Dt}}---\left(5\right)$

为了计算稳定，本文使用时间步上的细分积分形式来计算。其中Δτ为细分后的时间步长，为n时刻每个网格单元的SPH粒子个数，为目标粒子的个数。

第三步：求解压强p。根据第一步和第二步，方程(4)的右侧为已知，则方程(4)转化为关于压强p的线性方程组。使用预条件共轭梯度法高效求解该线性方程组即可获得压强解p。

第四步：获取n+1时刻的计算域G中的流体速度场uⁿ⁺¹。将p代入公式(3)，可得流体区域中新的速度场。最后用速度外推方法将流体区域中的速度外推到空气网格中获取计算域G中的速度场uⁿ⁺¹。

过程4：根据所获得计算域G中的速度场uⁿ⁺¹，计算更新后的符号距离场函数φⁿ⁺¹。

第一步：令当前时间步n对应的符号距离场函数为φⁿ，利用半拉格朗日法求解LS方程

$\frac{\∂\mathrm{\φ}}{\∂t}+u\·\▿\mathrm{\φ}=0---\left(6\right)$

得未校正的符号距离场函数这里，u＝uⁿ⁺¹为更新后的流体速度场。

第二步：使用检测校正粒子集P_c中的粒子对进行误差检测和校正，得符号距离场函数φⁿ⁺¹。具体的，令s_p为粒子符号，为n+1时间步的粒子位置，如果存在

$s_p{\mathrm{\&phi;}}_m^{n+1}\left(x_c^{n+1}\right)<0$

则说明产生粒子逃逸，即检测到流体表面计算存在误差。此时，需要利用逃逸粒子定义局部LS函数对流体表面进行校正。校正时，需要首先计算逃逸正负粒子集E⁺和E^-网格点的局部LS的最大值和最小值然后在逃逸粒子集所在网格点x中将两个距离场进行融合。如果|φ⁺(x)|≤|φ^-(x)|，则φ_c(x)＝φ⁺(x)，否则φ_c(x)＝φ^-(x)。将所得φ_c(x)值替换x位置处的距离场函数值便可以得到误差校正后的流体表面φⁿ⁺¹。

过程5：定义和计算局部中轴LMA。

第一步：构建封闭流形，由封闭流体界面或者非封闭的流体界面与计算域边界形成。封闭流形边界及其周围3层网格称为窄带区域。窄带定义为其中，为液体表面外侧三层网格区域，为液体表面内侧三层网格区域，表示液体表面。

第二步：定义局部中轴LMA为上述封闭流形位于窄带Ω_b上的内接圆圆心的轨迹。LMA为一个点集。

第三步：初始化LMA点集为空。

第三步：估计LMA集中的点p^M。为此，计算符号距离场函数φ的梯度的散度值。表示为

$D=\&dtri;\&CenterDot;(\&dtri;\mathrm{\&phi;})=\frac{{\&PartialD;}^2\mathrm{\&phi;}}{{\&PartialD;x}^2}+\frac{{\&PartialD;}^2\mathrm{\&phi;}}{{\&PartialD;y}^2}+\frac{{\&PartialD;}^2\mathrm{\&phi;}}{{\&PartialD;z}^2}=\frac{1}{\mathrm{\&Delta;}{h}^2}({\mathrm{\&phi;}}_{i-1,j,k}+{\mathrm{\&phi;}}_{i,j-1,k}+{\mathrm{\&phi;}}_{i,j,k-1}-6{\mathrm{\&phi;}}_{i,j,k}+{\mathrm{\&phi;}}_{i+1,j,k}+{\mathrm{\&phi;}}_{i,j+1,k}+{\mathrm{\&phi;}}_{i,j,k+1})$

其中，Δh为空间步长，一般定义为空间网格的大小。如果即网格单元的散度值大于一定阈值，说明该单元内含LMA点。然后用该网格中心点(网格点)x近似LMA上的点p^M，并入到已有的LMA点集中，即LMA＝LMA∪x。

过程6：定义定位量化函数(Location and Quantization Function,LF)。

第一步：定义界面上的点x到LMA的距离。表示为：

$L(x,\mathrm{LMA})=\min (\underset{p^M\&Element;\mathrm{LMA}}{\inf }\left(\right||x-p^M|\left|\right),\mathrm{\&lambda;})---\left(7\right)$

其中，inf为下界值，λ＝10Δh为阈值。

第二步：定义LF函数。令x为计算域上的任意一点，y为界面上的一点，则函数LF定义为：

$\mathrm{LF}\left(x\right)=\left\{\begin{array}{cc}L(x,\mathrm{LMA})& x\&Element;\&PartialD;\mathrm{\&Omega;}\\ k\left|\right|x-y\left|\right|+\mathrm{LF}\left(y\right)& x\&Element;{\mathrm{\&Omega;}}_b-\&PartialD;\mathrm{\&Omega;}\\ \&infin;& x\&Element;\mathrm{\&Omega;}-{\mathrm{\&Omega;}}_b\end{array}\right.---\left(8\right)$

其中，∞表示大常数，k为控制参数，控制界面内LF的值随距离变化的剧烈程度。y为界面上距离x的最近点，L(x,LMA)表示界面上的点x到LMA的距离。

过程7：计算表面上网格点的LF函数值。

第一步：构建KD树存储LMA点集。LMA的KD树构建使用LMA点集中具有最大跨度值的维度作为切分维度、用该维度上的均值作为切分值。

第二步：使用KD检索的方法，计算流体表面上网格点的LF函数值。从KD树中可以快速从LMA点集中获取距离x最近的点p^M。最后，根据公式(7)计算表面上网格点x的LF值。

过程8：计算区域上的LF函数值。

第一步：初始化网格点集合：建立三个集合Far，Close和Accepted；将所在的网格点置入Accepted中，并将这些网格点的位于上的邻居网格点置入Close中；将上的剩余的网格点置入Far中。

第二步：初始化LF值：Accepted中所有网格点的LF值为算法的输入，即上的LF值；Close集合和Far集合中所有网格点的LF值初始化为一个大常数，即LF＝∞。

第三步：建立以LF值为指标的最小堆，用于保存Close集合。

第四步：更新Close集合中所有网格点x的LF值：令Y＝{y}为x在Accepted集合中的邻居网格点集合，则x的LF值更新为

LF(x)＝min(LF(x),k||x-y||+LF(y))

第五步：根据重新计算的LF值更新保存Close集的最小堆。

第六步：选择最小堆顶元素Trial(即Close集合中具有最小LF值的网格点)，将其置入Accepted集合；将Trail中的所有位于Far集合中的邻居点转移到Close集合中。

第七步：转到第四步，迭代运行直到最小堆为空，转第八步。

第八步：上述Accepted集合中所有上的网格点的LF值即为所求。

过程9：计算区域上的LF函数值。

第一步：初始化网格点集合：建立三个集合Far，Close和Accepted；将所在的网格点置入Accepted中，并将这些网格点的位于上的邻居网格点置入Close中；将上的剩余的网格点置入Far中。

第二步：初始化LF值：Accepted中所有网格点的LF值为算法的输入，即上的LF值；Close集合和Far集合中所有网格点的LF值初始化为一个大常数，即LF＝∞。

第三步：建立以LF值为指标的最小堆，用于保存Close集合。

第四步：更新Close集合中所有网格点x的LF值：令Y＝{y}为x在Accepted集合中的邻居网格点集合，则x的LF值更新为

LF(x)＝min(LF(x),k||x-y||+LF(y))

第五步：根据重新计算的LF值更新保存Close集的最小堆。

第六步：选择最小堆顶元素Trial(即Close集合中具有最小LF值的网格点)，将其置入Accepted集合；将Trail中的所有位于Far集合中的邻居点转移到Close集合中。

第七步：转到第四步，迭代运行直到最小堆为空，转第八步。

第八步：上述Accepted集合中所有上的网格点的LF值即为所求。

过程10：计算窄带域以外的LF函数值。根据公式(8)，将其设置为一个大常数。

过程11：计算流体表面的变形程度值S。

第一步：定义表面变形函数S。流体表面变形函数S定义为一段时间流体表面所包围的流体体积的变化率

$S(t_i,t_{i-1})=\underset{t=t_{i-1},\mathrm{\&Delta;t}}{\overset{t_i}{\mathrm{\&Sigma;}}}\left(\right|A(t+\mathrm{\&Delta;t})-A\left(t\right)\left|\right)/A\left(t_{i-1}\right)---\left(9\right)$

其中，t_i-1和t_i分别上次重采样时刻和当前时刻，Δt为时间步长，A(t)为t时刻被流体界面包围的流体体积，表示为

$A\left(t\right)=\underset{{\mathrm{\&Omega;}}_b}{\&Integral;}\mathrm{\&delta;}\left({\mathrm{\&phi;}}^t\right)|\&dtri;{\mathrm{\&phi;}}^t|\mathrm{dx}---\left(10\right)$

其中，δ函数为Smeared-out Heaviside函数的导数

$\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)=\left\{\begin{array}{cc}0& \mathrm{\&phi;}<-\mathrm{\&epsiv;}\\ \frac{1}{2\mathrm{\&epsiv;}}+\frac{1}{2\mathrm{\&epsiv;}}\cos \left(\frac{\mathrm{\&pi;\&phi;}}{\mathrm{\&epsiv;}}\right)& -\mathrm{\&epsiv;}\&le;\mathrm{\&phi;}\&le;\mathrm{\&epsiv;}\\ 0& \mathrm{\&epsiv;}<\mathrm{\&phi;}\end{array}\right..$

这里，ε＝1.5Δh，Δh为空间步长。

第二步：计算的S值。在计算中，公式(10)可转化为

$A\left(t\right)=\underset{(i,j,k)\&Element;{\mathrm{\&Omega;}}_b}{\mathrm{\&Sigma;}}\frac{1}{\mathrm{\&Delta;h}}\mathrm{\&delta;}\left(\mathrm{\&phi;}(i,j,k,t)\right)\sqrt{\mathrm{\&Phi;}{\left(I\right)}^2+\mathrm{\&Phi;}{\left(J\right)}^2+\mathrm{\&Phi;}{\left(K\right)}^2}$

Φ(I)＝φ(i+1,j,k,t)-φ(i,j,k,t)，

Φ(J)＝φ(i,j+1,k,t)-φ(i,j,k,t)，

Φ(K)＝φ(i,j,k+1,t)-φ(i,j,k,t)，

其中(i,j,k)为窄带Ω_b上的网格点。这样，根据过程4符号距离场函数φⁿ⁺¹获取的参数，带入公式(9)(10)即可计算出表面变形程度S。

过程12：更新LF函数值。

第一步：确定当前时间步的性质。如果S的值小于阈值θ，说明表面变形不大，无需对几何特征进行重新计算，该时间步属于对流步；否则该时间步属于重新初始化步。

第二步：如果是对流步，则只需像流体表面一样，对其进行对流更新即可。此时LF的更新表示为

$\frac{\&PartialD;F}{\&PartialD;t}=-u\&CenterDot;\&dtri;F---\left(11\right)$

其中，F＝LF(x,t)，u＝uⁿ⁺¹为当前时刻已获取的网格G中的流体速度场。这样根据当前时刻的LF值LFⁿ和更新后的G中速度场uⁿ⁺¹，利用半拉格朗日法求解公式(11)，即可获得更新后的LF值LFⁿ⁺¹。

第三步：如果是重新初始化步，则需要按照过程(7)-过程(10)的方法重新计算LF函数值。

过程13：根据从几何模型计算中获得的特征函数，增加或删除粒子集中的粒子。

第一步：将当前时间步n的LF函数值LFⁿ映射为粒子的分布密度pⁿ(x)。由于LF在大曲率、细长区域以及边界区域的函数值较小，而这些区域所需的粒子反而较多，因此，先对LF进行变换

$p^n\left(x\right)=1.0-\frac{{\mathrm{LF}}^n\left(x\right)-\min \left({\mathrm{LF}}^n\right)}{\max \left({\mathrm{LF}}^n\right)-\min \left({\mathrm{LF}}^n\right)}.$

这样，LF值较小的网格(即大曲率、细长区域和界面边界等)pⁿ(x)值较大，从而满足在这些区域分布更多粒子的需求。然后将归一化，获取粒子的分布密度

$p^n\left(x\right)=\frac{p^n\left(x\right)}{\underset{x}{\mathrm{\&Sigma;}}{p}^n\left(x\right)}$

第二步：根据pⁿ(x)计算需要增/删的粒子数为

$N_{\mathrm{del}}^n\left(x\right)=\mathrm{\&kappa;}(\underset{x}{\mathrm{\&Sigma;}}{N}_p^n\left(x\right)\&times;p^n\left(x\right)-N_p^n\left(x\right))$

其中，为当前网格单元x中的粒子数，为计算域中粒子总个数，κ为控制参数，控制粒子个数变化的剧烈程度。

第三步：在网格单元x中增加或删除个粒子，得到更新后的粒子集和在每个网格单元x中，如果则在网格内随机删除个粒子。否则，随机生成个粒子，粒子位置在网格单元x内部随机产生，速度初始化为零。

第四步：确定当前时间步n下的新粒子的分类并分别置入和中置入和中以获得和对于新粒子p，如果即该粒子位于区域，则属于其符号为负。否则，则该粒子位于区域。这里，为过程(4)中获取的未校正的符号距离场函数。如果粒子位于界面的ω层(ω＝3)网格内，即则该粒子属于其符号为正。如果则粒子属于

过程14：更新中粒子的速度和位置，对实施误差检测校正获取符号距离场函数φⁿ⁺¹。

第一步：更新中粒子的速度。对中的粒子，其速度可以利用三次线性插值的方法从G中速度场uⁿ⁺¹获取。即其中为粒子的位置。

第二步：更新中粒子的位置。以作为输入，通过求解方程即可获得更新后的粒子位置如果中的负粒子逃逸到液体表面三层网格之外的空气网格中，则将其从中移除，添加到中。

第三步：获取n+1时间步中新的粒子集上述中所有粒子速度和位置更新后，即可得新的粒子集

第四步：利用新的粒子集检测和校正获得由于为过程(4)中获取的未校正的符号距离场函数。由于本发明中粒子集的分布是根据几何特征计算所得的，因此它对中的误差可以进行更精确的检测和更准确的校正，提高流体表面的精度。具体的，令s_p为粒子符号，如果存在

$s_p{\mathrm{\&phi;}}_m^{n+1}\left(x_c^{n+1}\right)<0;$

则说明粒子逃逸，即检测到流体表面计算存在误差。此时，需要利用逃逸粒子定义局部LS函数对流体表面进行校正。校正时，需要首先计算逃逸正负粒子集E⁺和E^-网格点的局部LS的最大值和最小值然后在逃逸粒子集所在网格点x中将两个距离场进行融合。如果|φ⁺(x)|≤|φ^-(x)|，则φ_c(x)＝φ⁺(x)，否则φ_c(x)＝φ^-(x)。将所得φ_c(x)值替换x位置处的距离场函数值便可以得到误差校正后的流体表面φⁿ⁺¹。

过程15：更新中粒子的速度和位置，得下一个时间步n+1的水花粒子集

第一步：更新中粒子的速度由于uⁿ⁺¹的计算是将粒子速度投影到网格中并求解网格和粒子耦合的压强方程获取的。因此，只需将uⁿ⁺¹反投影到粒子上即可获得新的粒子速度。具体的，首先把中间速度场u^*到uⁿ⁺¹的增量反投影到粒子上，即Δu_p＝Σ(uⁿ⁺¹-u^*)ω_p(x)，则更新后的粒子速度为参数s表示粒子速度受压强的影响程度。

第二步：更新中粒子位置根据求解方程可得更新后的粒子位置如果粒子运动到流体表面以内，即则从中移除该粒子。

第三步：获取n+1时间步中新的粒子集上述中粒子速度和位置更新后，即得新的粒子集

过程16：从上述计算结果φⁿ⁺¹和中提取表示流体表面和流体水花的三角形网格表面，并对其进行真实感渲染后即可生成流体动画。

第一步：使用移动立方体算法，从表示流体表面的符号距离场函数φⁿ⁺¹中提取零等值面，获得流体表面的三角形网格表面。

第二步：使用构建粒子局部水平集算法，从粒子集中获取表示水花粒子的符号距离场函数定义为：

${\mathrm{\&phi;}}_p^{n+1}\left(x\right)=|x-\stackrel{\&OverBar;}{x}|-\stackrel{\&OverBar;}{r}$

其中，为粒子周围邻居粒子的平均位置，记为为粒子周围邻居粒子的平均半径，记为这里i为当前粒子的邻居粒子标号。w_i为权重值。

$w_i=\frac{k\left(\right|x-x_i|/R)}{\underset{i}{\mathrm{\&Sigma;}}k\left(\right|x-x_j|/R)}$

这里，k为核函数定义为k(d)＝max(0,(1-d²)³)。

获取水花粒子的符号距离场函数之后，再使用移动立方体算法，获得零等值面，该等值面为水花的三角形网格表面。

第三步：使用全局光照算法，对获取的流体表面和流体水花的三角形网格表面进行真实感渲染，即可获得流体动画。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。

Claims (7)

1.一种基于几何特征的高精度流体动画建模方法，其特征是：包括以下步骤:

(1)对流体物理模型中的初始条件和边界条件进行初始化；

(2)从流体的物理模型中获取几何特征计算所需要的中间速度场、当前时间步中的流体速度场和符号距离场函数参数；

(3)构建流体动画的几何模型，对流体动画中的细节敏感区域进行自动化定位、量化和跟踪；

(4)根据几何特征的特征函数，控制物理模型的粒子集，耦合几何模型和物理模型；

(5)提取表示流体表面和流体水花的三角形网格表面，并对其进行真实感渲染，生成一帧流体动画；

(6)查看是否完成全部流体动画的制作，如果未完成，则返回步骤(2)，继续进行制作。

2.如权利要求1所述的一种基于几何特征的高精度流体动画建模方法，其特征是：所述步骤(1)中，物理模型初始化的步骤包括：

(1-1)初始化背景网格中的流体速度场u；

(1-2)初始化表示流体表面的符号距离场函数φ；

(1-3)初始化对流体表面进行检测校正的粒子集P_c和表示流体水花的粒子集P_s；

(1-4)初始化边界，即使用标记位表示流体计算域中的每个网格单元的类型；为此需要将固体障碍物按计算域网格单元尺寸离散化，离散化后被固体占据的网格单元标记位设置为1，否则设置为0。

3.如权利要求1所述的一种基于几何特征的高精度流体动画建模方法，其特征是：所述步骤(2)中，从流体的物理模型中获取几何特征的具体方法如下：

(2-1)根据初始化背景网格中的流体速度场u，粒子集P_s中的速度计算背景网格中的中间速度场；

(2-2)根据中间速度场，通过求解耦合欧拉和拉格朗日两种模型的压强泊松方程，获取当前时间步中的流体速度场uⁿ⁺¹；

(2-3)根据所获得的当前时间步流体速度场，计算更新后的符号距离场函数φⁿ⁺¹。

4.如权利要求1所述的一种基于几何特征的高精度流体动画建模方法，其特征是：所述步骤(3)中，建立流体动画的几何模型，对流体动画中的细节敏感区域进行自动化定位、量化和跟踪的基本步骤如下：

(3-1)构建封闭流形，由封闭流体界面或者非封闭的流体界面与计算域边界形成，封闭流形边界及其周围3层网格称为窄带区域；

(3-2)定义局部中轴为封闭流形中位于窄带Ω_b上的内接圆圆心的轨迹；

(3-3)根据Step2中获得的符号距离场函数φⁿ⁺¹，使用散度值判断法计算LMA；

(3-4)根据局部中轴，定义定位量化函数；

(3-5)根据定位量化函数的定义，通过使用KD检索的方法，计算流体表面上网格点的定位量化函数值；

(3-6)根据定位量化函数的定义，通过使用快速行进算法对流形边界点的定位量化函数值外推，计算窄带域上其他网格点的定位量化函数值；

(3-7)根据定位量化函数的定义，设定背景网格中窄带域以外的定位量化函数值为大常数，指定为N_max×Δh，其中N_max为计算域在x、y、z三个轴向上网格个数的最大值，Δh为网格单元的尺寸；

(3-8)在每个时间步中累积计算流体表面的变形程度，如果变形程度小于设定阈值，则需要对定位量化函数进行对流更新；否则需要对定位量化函数进行重新初始化。

5.如权利要求4所述的一种基于几何特征的高精度流体动画建模方法，其特征是：所述步骤(3-4)中，具体方法为：根据局部中轴，定义定位量化函数，对于流形边界上的网格点，其函数值为给定阈值λ与该点到局部中轴距离值的较小值；窄带区域上其他网格点的定位量化函数值为该点到流形边界的最短距离加上流形边界对应点的定位量化函数值；对于计算域上其他网格点的定位量化函数值将其设置为一个大常数，定位量化函数能够实现对大曲率、细长和边界细节敏感区域的光滑定位和量化。

6.如权利要求1所述的一种基于几何特征的高精度流体动画建模方法，其特征是：所述步骤(4)中，其具体方法为：

(4-1)将当前时间步n的定位量化函数值LFⁿ映射为粒子的分布密度pⁿ(x)；

(4-2)根据pⁿ(x)计算需要增/删的粒子数

(4-3)在网格单元x中增加或删除个粒子，得到更新后的粒子集和

(4-4)更新中粒子的速度和位置，然后利用它们检测和校正的误差，得符号距离场函数φⁿ⁺¹；

(4-5)更新中粒子的速度和位置，得下一个时间步n+1的水花粒子集

7.如权利要求1所述的一种基于几何特征的高精度流体动画建模方法，其特征是：所述步骤(5)中，具体方法包括：

(5-1)使用移动立方体算法，从表示流体表面的符号距离场函数φⁿ⁺¹中提取零等值面，获得流体表面的三角形网格表面；

(5-2)使用构建粒子局部水平集算法，从粒子集中获取表示水花粒子的符号距离场函数然后再使用移动立方体算法，获得表示水花的零等值面。