CN101484922A - 使用导数质点来模拟流体详细运动的方法 - Google Patents

Abstract

使用导数质点来模拟流体详细运动的方法。本发明提供了一种方法，该方法提供了一种使用质点和导数信息来显著降低速度的非物理耗散的新流体模拟技术。在求解传统纳维尔－斯托克斯方程的过程中，所述方法用质点模拟代替平流部分。当在基于网格和基于质点的模拟器之间切换时，必须转换诸如水平集和速度的物理量。开发了一种利用存储在质点中以及网格点中的导数信息的新型耗散抑制转换程序。通过多次实验，提出的技术可以再现诸如液滴/泡、薄水膜以及漩涡等高雷诺数流体的详细的运动。平流中增加的精度还能够用来在更大尺度的流体模拟中产生更好的结果，所述平流中增加的精度形成提出的技术的基础。

Description

使用导数质点来模拟流体详细运动的方法

技术领域

本发明涉及一种使用导数质点来模拟流体的详细的运动的方法。

背景技术

1 介绍

当水猛烈地与固体、空气、或水本身相互作用时，其呈现各种结构，包括液滴/泡、薄水膜、以及漩涡，如图7中所示。当流体经历以高雷诺数为特征的运动时能够出现这样的特征，其中，雷诺数表示惯性与流体的粘度之比的相对大小。高速移动的水是典型的高雷诺数流体。本发明探索这样的现象的基于物理的模拟。

假设纳维尔-斯托克斯方程正确地模拟流体的物理运动，看似真实的模拟应能够再现水的高雷诺数行为。虽然本工作主要讨论水，但也可应用于任何流体。然而，迄今为止，尚未令人满意地再现这些现象。本发明认识到导致这种失败的因素，并提出一种允许高雷诺数液体运动的真实的模拟的方法。

看似不真实的高雷诺数流体的模拟与数值耗散有关。具体地说，纳维尔-斯托克斯方程的离散化模拟不可避免地需要估计非网格点处的物理量。在大多数方法中，通过网格点处的物理量的值的插值来计算这样的值。然而，由此近似引入的误差起到类似于非物理粘度或数值耗散的作用。能够通过使用更细的网格来降低这种耗散；然而，增加网格分辨率会将计算时间和存储器需求增加至不切实际的水平。在过去的几年中，已提出多种优秀的技术来解决这个问题。将质点引入欧拉方案可以帮助捕捉流体的惯性运动，并增加模拟的有效分辨率。Enright等人[2002b]介绍了质点水平集方法，其通过在界面周围散布质点来增加表面追踪的准确性。Losasso等人[2004]提出了一种基于八叉树的多水平流体解算器，其允许在诸如水表面的更感兴趣的区域中的更精细分辨率的模拟。

不幸的是，以上技术不能产生具有充分细节和真实性的高雷诺数液体。被模拟的流体表现出比在真实物理中更大的粘度；流体常常看起来厚/粘，而且通常在复杂流动中不表现出小尺度特征的运动。申请人发现以上耗散抑制技术产生这种人为瑕疵的原因是它们不能有效地抑制速度耗散，即使它们显著地降低了质量耗散。

在本发明中，申请人介绍了一种称为导数质点(derivative particle)的新概念，并开发了一种基于该概念的流体模拟器。申请人利用用于平流部分的模拟的拉尔朗日方案的无耗散性质；对于需要模拟细节的流体区域，申请人使用质点来求解平流步骤。基于网格与基于质点的模拟器之间的切换可能引入数值耗散。本发明开发了一种允许详细流体运动的再现的新的转换程序。此工作的一个创新方面在于，除了存储物理量(速度和水平集值)之外，所述导数质点还存储那些量的空间导数，这使得能够实现非网格/非质点位置处的物理量的更精确评估。本工作的质点的使用与质点水平集方法[Enright等人2002b]的不同之处在于导数质点携带流体速度以及水平集值。

实验显示提出的模拟器准确地追踪界面。更有趣的是，提出的方法被证实显著降低非物理阻尼，允许在真实的高雷诺数流体中发生的小尺度特征的宏观运动的再现。

本说明书的其余部分如下组织：第2部分回顾先前的工作；第3部分给出所述模拟器的概述；第4部分介绍基于八叉树的约束插值剖面(CIP)解算器；第5部分介绍导数质点模型；第6部分报告我们的实验结果；第7部分总结本说明书。

2 先前的工作

由于Foster和Metaxas[1996；1997]首先介绍了基于完全3D纳维尔-斯托克斯模拟的流体动画技术，所以该方法已在计算机图形学界盛行。Jos Stam[1999]介绍了一种称为稳定流体(Stable Fluids)的稳定流体解算器。该解算器的平流步骤使用半拉格朗日方法[Staniforth和C^ot`e 1991]来实施，其即使在使用大时步时也保持稳定。从那时起，已有计算机图形学中基于半拉格朗日方法的快速流体动画技术的积极开发。Rasmussen等人[2003]提出了一种使用2D半拉格朗日解算器来产生气体的3D大尺度动画的技术，且Feldman等人[2003]提出了一种将基于质点的燃烧模型并入稳定流体解算器的爆炸模型。Treuille等人[2003；2004]介绍了一种生成满足指定的关键帧约束的流体流动的优化技术。另外，Feldman等人提出了结合交错网格和非结构网格[Feldman等人2005a]的混合网格的使用，而且还提出了可变形网格的使用[Feldman等人2005b]，对于该方法他们必须修改半拉格朗日方法。

为了模拟流体，除稳定的纳维尔-斯托克斯解算器之外，还需要用于追踪液体表面的模型。为此，Foster和Fedkiw[2001]提出了一种将无质量质点引入到水平集场的混合表面模型。这种模型推动了质点水平集方法[Enright等人2002b]的开发，所述质点水平集方法能够以显著的准确率捕捉流体表面的动态运动。Enright等人[2005]证明该质点水平集方法允许半拉格朗日方案中的大时步的使用。该方法已被用于建模流体与刚性体之间的相互作用[Carlson等人2004]以及流体与可变性薄壳物体之间的相互作用[Guendelman等人2005]，而且被用于模拟粘弹性流体[Goktekin等人2004]、沙[Zhu和Bridson 2005]、多相流体[Hong和Kim 2005]、以及水滴[Wang等人2005]。

作为基于网格的方法的替代，还已研究了纯粹基于质点的方法。Stam和Fiume[1995]使用光滑质点流体动力学(SPH)来建模气态现象。在SPH中，通过颗粒的集合来表现流体，且模拟器通过计算纳维尔-斯托克斯方程的每一项来计算它们的运动。M¨uller等人[2003]使用SPH模型来模拟流体，并且在[M¨uller等人2005]中，他们使用该技术来模拟多相流体相互作用。Premoˇze等人[2003]将移动质点半隐式(MPS)法引入到图形学界，这是一种在模拟流体的不可压缩运动的过程中显示出比SPH更好的性能的技术。然而，纯粹基于质点的方法的基本缺陷是它们缺少表面追踪能力；如果使用不适当的质点数目，则它们可能产生粒状或过度光滑表面。

削弱模拟结果的视觉真实性(以及物理准确性)的主要因素是速度场的数值耗散。为了降低模拟气态现象时的耗散，Fedkiw等人[2001]使用三次插值。它们还包括称为涡度约束(vorticity confinement)的附加步骤，其将速度场的漩度(curl)放大，产生烟气运动中的真实的涡漩形部分。通过采用涡流(vortex)质点法来执行涡度耗散的更基于物理的预防[Selle等人2005；Park和Kim 2005]。这种方法用漩度型纳维尔-斯托克斯方程来工作并利用质点来求解(solve)平流项，这导致涡度的有效保持。然而，在液体的情况下，粘度表现出更突出的作用；特别是，涡漩形运动存在时间短而且较少观察到。因此，建模高雷诺数液体的行为需要在更一般的情况下工作的速度耗散抑制技术。为了降低液体模拟中的耗散，Song等人[2005]提出了一种基于CIP平流方法的技术。他们的技术用三阶精度来求解速度平流。然而，通过这种方法，不能避免来自基于网格的平流的数值粘度。Zhu和Bridson[2005]提出了一种基于流体隐式质点(FLuid-Implicit-Particle)(FLIP)方法的流体模拟技术。该FLIP方法[Brackbill和Ruppel 1986]用质点来求解平流步骤，而不包括网格的所有其它步骤。虽然这种方法中的质点平流没有数值耗散，但当在网格与质点之间来回传递速度值时，引入了插值误差。

因此，长期以来存在对使用导数质点来模拟流体的详细运动的方法的需要。本发明意在解决这些问题并满足期盼已久的需要。

发明内容

本发明设法解决现有技术的缺点。

本发明的目的是提供一种使用导数质点来模拟流体的详细运动的方法。

本发明的另一目的是提供一种使用导数质点来模拟流体的详细运动的方法，其中，用拉格朗日方案来实施平流部分。

本发明的又一目的是提供一种使用导数质点来模拟流体的详细运动的方法，其提供在平流与非平流部分的切换处转换物理量的程序，抑制不必要的数值耗散的介入。

3 概述

在多相流体解算器的开发中，申请人假设空气和水两者都是不可压缩的。不可压缩流体的纳维尔-斯托克斯方程可以写成

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

$\▿\·u=0---\left(2\right)$

其中u是速度，p是压力，f是外力，μ是粘度系数，ρ是流体密度。为了准确地建模跨越介质之间的界面的密度和粘度的不连续性，申请人采用具有可变密度的虚拟流体方法[Liu等人2000；Kang等人2000；Hong和Kim 2005]。在我们的解算器中，表面追踪是基于水平集方法。根据下式而更新水平集场φ

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

符号距离条件

$|\▿\mathrm{\φ}|=1---\left(4\right)$

对于纳维尔-斯托克斯方程的时间积分，申请人采用分步法[Stam1999]，其递增地计算(accounts for)等式(1)中的项的效果和等式(2)中的质量守恒条件。申请人开发的技术在这里关注于平流部分的增加。因此，如图8所示，申请人将分步(fractional steps)分组成两部分，非平流部分和平流部分。等式(3)的时间积分仅涉及平流部分。

平流部分由网格平流部分和质点平流部分组成。所述网格平流部分平流输送根据当前速度场从非平流部分计算的速度和水平集场。我们的方法与纯粹欧拉模拟的不同之处在于我们的模拟器包括质点平流部分，其用于模拟需要产生细节的区域。为此，申请人在空气-水界面中引入了大量质点。

在网格平流部分中，平流输送速度和平流输送水平集是欧拉平流步骤。这些步骤的精确处理形成高雷诺数行为的成功模拟的基础。对于欧拉平流，申请人开发了一种基于八叉树的CIP解算器(第4部分)，其将速度和质量耗散降低到显著的水平。

当流体区域的平流需要使用质点平流部分来模拟时，当前在网格上限定的物理量(例如水平集和速度)被输送到质点。在该输送之后，质点能够被直接地平流输送。质点平流的结果然后被传回到网格。在程序中的这一点上，不管模拟是经由网格平流部分还是质点平流部分进行的，最终速度和水平集值都存储在网格点上。申请人在第5部分中介绍了网格至质点和质点至网格速度/水平集转换程序。基于质点的平流的使用和转换程序的开发是再现流体运动的细节所必不可少的。

一种使用导数质点来模拟流体的详细运动的方法包括以下步骤：a)基于八叉树数据结构用自适应网格来建模流体；b)对于网格点处的流体速度求解不可压缩流体的纳维尔-斯托克斯方程；c)在纳维尔-斯托克斯方程的时间积分中使用分步法；d)将分步分组成非平流部分和平流部分；e)根据模拟的详细程度来选择网格平流部分或质点平流部分；f)对于网格平流部分使用基于八叉树的CIP方法；g)对于质点平流部分使用导数质点模型。

当模拟的详细程度高时，使用对于质点平流部分使用导数质点模型的步骤。

对于高度感兴趣的区域的细节自适应网格增加了网格的数目。

八叉树数据结构包括树形数据结构，每个内部结点具有多达八个孩子。

不可压缩流体的纳维尔-斯托克斯方程可以写成

$\frac{\∂u}{\∂t}=-(u\·\▿)u-\frac{\▿p}{\mathrm{\ρ}}+\▿\·(\mathrm{\μ}\▿u)+\frac{f}{\mathrm{\ρ}},$

$\▿\·u=0,$

其中u是速度，p是压力，f是外力，μ是粘度系数，ρ是流体密度。

该方法还可以包括使用水平集方法来追踪流体的表面的步骤。所述水平集方法包括水平集场φ。根据下式来更新所述水平集场φ：

$\frac{\∂\mathrm{\φ}}{\∂t}+u\·\▿\mathrm{\φ}=0$

符号距离 $|\▿\mathrm{\φ}|=1.$

$\frac{\∂\mathrm{\φ}}{\∂t}+u\·\▿\mathrm{\φ}=0$ 的时间积分仅涉及平流部分。网格平流部分平流输送根据当前速度场从非平流部分计算的速度和水平集场。

网格平流部分使用欧拉平流步骤来平流输送速度和水平集。通过基于八叉树的CIP(约束插值剖面)法来执行欧拉平流。

对于网格平流部分使用基于八叉树的CIP方法的步骤包括用基于八叉树的CIP方法来平流输送速度和水平集及其导数的步骤。直接从方程 $\frac{\∂\mathrm{\φ}}{\∂t}+u\·\▿\mathrm{\φ}=0$ 获得平流输送导数的方程，其中，φ是被平流输送的物理量，且通过对x微分该方程来获得空间变量x的导数的平流输送方程，即：

$\frac{\∂{\mathrm{\φ}}_x}{\∂t}+u\·\▿{\mathrm{\φ}}_x=-u_x\·\▿\mathrm{\φ}.$

求解x的导数的平流输送方程的步骤包括采用使用分步法的步骤，其使用分步法的步骤包括以下步骤：a)使用有限差分法来求解非平流部分 $\∂{\mathrm{\φ}}_x/\∂t=-u_x\·\▿\mathrm{\φ};$ b)平流输送根据 $\frac{\∂{\mathrm{\φ}}_x}{\∂t}+u\·\▿{\mathrm{\φ}}_x=0$ 的结果；以及c)在由网格限定的单元(cell)的中心处进行压力值取样并在结点处进行包括速度、水平集及其导数的其它值的取样。

该方法可以进一步包括以下步骤：a)将基于八叉树的CIP方法与自适应网格相结合；以及b)通过用基于八叉树的CIP方法模拟平流来使由八叉树数据结构的网格分辨率中的区域性变化引起的八叉树人为瑕疵最小化。

导数质点模型使用基于质点的拉格朗日方案。用 $\mathrm{\φ}\left(x\right)={\mathrm{\φ}}_p+\frac{\▿{\mathrm{\φ}}_p}{|\▿{\mathrm{\φ}}_p|}\·(x-x_p)$ 和 $\▿\mathrm{\φ}\left(x\right)=\▿{\mathrm{\φ}}_p$ 来计算网格点处的水平集及其导数，其中p是由导数质点携带的梯度。

该方法可以进一步包括以下步骤：a)在对于质点平流部分使用导数质点模型的步骤之前将网格上限定的速度和水平集转换成质点上限定的速度和水平集(网格至质点)；以及b)在返回到对于网格平流部分使用基于八叉树的CIP模型的步骤之前将质点上限定的速度和水平集转换成网格上限定的速度和水平集(质点至网格)。

转换速度和水平集网格至质点的步骤包括下述步骤，即对于速度u^P的转换使用具有用单调CIP插值来代替线性插值的FLIP方法的步骤，且 $u^P=u_{-}^P+{\mathrm{\Σ}}_{i=1}^4{w}_i\mathrm{\Δ}{u}_i^G,$ 其中

是就在非平流部分开始之前的质点速度，且

是由于非平流部分而引起的G_i处的速度变化。

对于速度u^P、网格点G、选自G的每个象限的最近质点P₁、P₂、P₃和P₄、P_i的速度以及

和

的x和y分量的导数，转换速度和水平集质点至网格的步骤包括以下步骤：a)将

的导数坐标变换成

从而表示沿着方向

的分量且表示垂直于

的分量，其中，类似地获得

的坐标变换的导数

b)沿着方向对步骤(a)中获得的结果进行一维单调CIP插值以便计算位置A处的u_A及其

方向导数M；c)通过应用上述坐标变换的反转(inverse)来从

获得(u_Ax，u_Ay)；d)对质点P₃和P₄执行步骤(a)～(c)以计算B的u_B和(u_Bx，u_By)；以及e)对步骤(c)和(d)中获得的结果执行y方向单调CIP插值，其给出G处的速度和导数。

对于质点平流部分使用导数质点模型的步骤包括以下步骤：a)调整质点的水平集值；以及b)执行质点重新播种。

本发明的优点是(1)所述使用导数质点来模拟流体的详细运动的方法降低质量和速度两者的耗散，导致动态流体的真实再现；(2)所述使用导数质点来模拟流体的详细运动的方法用拉格朗日方案来实施平流部分；以及(3)所述使用导数质点来模拟流体的详细运动的方法提供质点至网格和网格至质点程序。

虽然简要地概述了本发明，但通过以下附图、详细说明及权利要求书，能够获得对本发明的更全面的理解。

附图说明

参照附图将更透彻地理解本发明的这些及其它特征、方面和优点。

图1是示出了根据本发明的使用导数质点来模拟流体的详细运动的方法的流程图；

图2是示出了根据本发明的方法的另一流程图；

图3是示出了根据本发明的方法的又一流程图；

图4是示出了网格至质点和质点至网格转换的一部分流程图；

图5是示出了网格至质点和质点至网格转换的另一部分流程图；

图6是示出了使用导数质点模型的步骤的细节的一部分流程图；

图7示出了由球的冲击而产生的水的模拟结果；

图8示出了根据本发明的方法的体系结构；

图9示出了从质点速度到网格速度的计算；

图10示出了从2D溃坝(breaking-dam)模拟截取的截图：上部和下部序列分别示出了用质点水平集方法的传统线性模型和导数质点模型产生的结果；

图11示出了从水滴模拟截取的截图：左图像和右图像分别示出了用质点水平集方法的传统线性模型和导数质点模型产生的结果；

图12示出了从旋转的水箱的模拟截取的截图；以及

图13示出了从洪水模拟截取的截图：底端行的三个图像示出了洪水的推进(顶视图)；顶部图像示出了侧视图。

具体实施方式

4 开发基于八叉树CIP解算器

本部分描述通过将CIP方法与八叉树数据结构相结合而进行的模拟器的网格平流部分的开发。

4.1 CIP方法介绍

在用分步法来求解纳维尔-斯托克斯方程(第3部分)的过程中，平流项

和

由于其双曲线性质而需要特别注意。半拉格朗日方法提供了一种稳定方案，在该方案内处理以上双曲线方程[Stam 1999]。然而，不幸的是，这种方法遭遇严重的数值耗散，该数值耗散源于用来根据相邻网格点处的物理量而确定回推(backtracked)的非网格点处的物理量的线性插值。

同时CIP方法是最初由Yabe和Aoki[1991a；1991b]提出并随后由Yabe等人[2001]改进的三阶技术。这种方法的关键思想是不仅平流输送其物理量，而且还平流输送其导数。这里，问题将是如何平流输送导数。Yabe和Aoki[1991a]的感兴趣的观察结果是能够直接从初始双曲线方程获得用于平流输送导数的方程

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

其中，φ是正在被平流输送的物理量。对于空间变量x微分方程(5)，申请人得到

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

其可以用来预测x随时间的演变。在求解方程(6)的过程中，申请人再次应用分步法：首先，申请人使用有限差分来求解非平流部分 $\∂{\mathrm{\φ}}_x/\∂t=-u_x\·\▿\mathrm{\φ};$ 然后，申请人平流输送根据下式的结果

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

对CIP平流的详细介绍可以在[Song等人2005]中找到。以上方法的二维和三维情况的扩展存在于[Yabe和Aoki 1991b]。Song等人[2005]发现[1991b]中给出的推广可能引起不稳定，并提出解决这个问题的修改。在本工作中，为了求解平流部分，申请人采用带有二阶Runge-Kutta时间积分的CIP方法。

4.2 将CIP与自适应八叉树网格相结合

由Losasso等人[2004]介绍的基于八叉树数据结构的自适应网格的使用允许在全部流体中以非齐次精度来执行模拟。从实践观点看，这种方法是有用的，因为其允许通过引入少量的附加计算和存储器来以更高的精度模拟诸如表面的更感兴趣的区域。这里，申请人采用这种八叉树数据结构，并提出这种自适应方法的实际值能够通过将其与CIP方法相结合而得到进一步改善。

Losasso等人[2004]对于半拉格朗日平流步骤使用线性插值。然而，如果申请人用CIP插值来代替线性插值，则平流将具有三阶精度。如在Guendelman等人[2005]中一样，申请人在单元中心处对压力值进行取样，但是在结点处对所有其它量—速度、水平集及其导数—进行取样。当单元被细化时，通过参照所有导数值来对新网格点处的值进行CIP插值。

我们注意到CIP方法非常好地与八叉树数据结构相配合，因为虽然其具有三阶空间精度，但其被限定在单一网格单元模板(stencil)而不是多个模板上。由于这种紧凑性，申请人能够在没有任何主要修改的情况下将CIP方法用于模拟自适应网格。相反，对于被限定在多个模板上的高阶方案，难以将该方法从一阶扩展到三阶：由于模拟自适应地细化网格，所以将产生具有不同尺寸的单元，这使得多模板高阶方案的开发成为难题。

我们还将注意到一种申请人称为八叉树人为瑕疵的特征。八叉树数据结构的网格分辨率的区域性变化产生耗散的量的变化。质量耗散的区域性变化不显著。然而，对于速度耗散，区域性差异可能显著，尤其是对于快速流体运动。例如，在溃坝的模拟中，水经历快速的横向移动。接近分辨率高的表面的流体经历少量的数值扩散并因此而进行快速移动。相反，底部的流体由于低的网格分辨率而以更粘流体的特征的方式移动。当这两类运动一起出现时，水的上部好像是在下部上爬行。包括CIP在内的高阶方案避免不了这种人为瑕疵。然后，当使用本工作中提出的基于八叉树的CIP解算器来模拟平流时，所述人为瑕疵较不显著；申请人将这种改善归因于数值耗散的总体降低，尤其是在低分辨率区域中。

5 导数质点模型

在动量/质量守恒方面，基于质点的拉格朗日方案比基于网格的欧拉方案更精确。因此，申请人将拉格朗日方案应用于需要更详细地模拟的区域。然而，这种方法在表面追踪和压力/粘度计算方面具有其本身的局限性。因此，申请人仅将该方法应用于平流部分的模拟；模拟的所有其它部分均基于欧拉方案。拉格朗日与欧拉方案之间的切换要求物理量(速度和水平集值)的网格至质点和质点至网格转换。应小心地开发用于执行这些转换的程序以便它们不引入不必要的数值扩散。在本部分中，申请人介绍了不削弱拉格朗日平流的需要的特性的新型转换程序，其是使得高雷诺数流体中的小尺度特征能够再现的必要的组件。为简单起见，对2D进行本部分的说明。

5.1 网格至质点速度转换

在图8中所示的平流部分的末端，(1)需要被平流输送的速度存储在网格点处，而且(2)大量质点散布在网格上。质点平流部分必须找出质点的速度从而其速度和水平集值产生拉格朗日运动。假设质点P处于由四个网格点G₁、G₂、G₃与G₄限定的单元中。对于i∈{1，2，3，4}，使u^G _i＝(u^G _i，v^G _i)为G_i处的网格速度。申请人正在寻找能够用于P的速度u^P＝(u^P，v^P)的公式。

在单元中质点方法[Harlow 1963]中使用的一种可能方法是使用双线性插值 $u^P={\mathrm{\Σ}}_{i=1}^4{w}_i{u}_i^G,$ 其中，w_i是基于质点位置P而确定的双线性加权。不幸的是，该转换引入严重的数值扩散。在FLIP方法[Brackbill和Ruppel 1986]中，仅网格速度的递增部分被传递质点速度。即， $u^P={u_{-}^P+\mathrm{\Σ}}_{i=1}^4{w}_i{\mathrm{\Δu}}_i^G,$ 其中，是就在非平流部分开始之前的质点速度，且

是由于非平流部分而引起的G_i处的速度变化。这种方法已被示出显著降低数值扩散[Zhu和Bridson 2005]。因此，在网格至质点速度转换程序的开发中，申请人使用FLIP方法，但是用单调CIP插值[Song等人2005]来代替线性插值。这不仅给出u^P，而且给出

$(\frac{\∂}{\∂x}{u}^P,\frac{\∂}{\∂y}{u}^P)\mathrm{and}(\frac{\∂}{\∂x}{v}^P,\frac{\∂}{\∂y}{v}^P)$

5.2 质点至网格速度转换

假设每个质点根据速度u^P移动，携带速度的空间导数以及速度本身。申请人现在必须开发将质点平流的结果传送到网格的程序。质点至网格速度转换比网格至质点转换更复杂。

参照图9，使G为网格点。申请人从G的每个象限选择最近的质点，即P₁、P₂、P₃和P₄。使为P_i的速度，并使

和

分别为

的x和y分量的导数。申请人必须找到G处的速度u^G＝(u^G，v^G)及其空间导数的公式。由于所述四个质点没有成矩形地放置，所以申请人不能使用传统的基于网格的CIP插值。

我们的质点至网格速度转换由以下步骤组成。为简单起见，申请人仅示出x分量的计算。

(1)申请人将

的导数

坐标变换成，从而，

表示沿方向的分量，且u_1⊥|表示垂直于

的分量。同样地，申请人得到

的坐标变换的导数

(2)申请人沿着方向

对步骤(1)中得到的结果执行一维单调CIP插值以便计算位置A处的u_A及其

方向的导数

(3)申请人通过以上坐标变换的反转来从

获得(u_Ax，u_Ay)。

(4)同样地，申请人对质点P₃和P₄执行步骤(1)～(3)以计算B的u_B和(u_Bx，u_By)。

(5)申请人对步骤(3)和(4)中得到的结果执行y方向单调CIP插值，这给出G处的速度和导数。

这种单调插值法能够直接扩展到3D情况。申请人注意到本工作未采用一般径向基插值方案，因为(1)他们不利用导数信息，而且(2)已证明难以修改方案以便它们体现导数但保持单调性。

5.3 水平集转换

应不同于速度转换来执行水平集转换；水平集值表示到表面的最短距离，所以，该换算将不基于插值法。这里，申请人使用广泛使用的水平集转换程序[Enright等人2002a；Enright等人2002b]的修改的且更准确的形式。

在初始的质点水平集方法中，使用球隐函数φ(x)＝s_p(|φ_p|-|x-x_p|)来计算网格点x处的水平集值，其中，p是质点的水平集，且对于正(负)质点，s_p＝+1(-1)。在本工作中，利用存储在导数质点中的导数信息，申请人用下式来计算该网格点处的水平集及其导数

$\mathrm{\φ}\left(x\right)={\mathrm{\φ}}_p+\frac{\▿{\mathrm{\φ}}_p}{|\▿{\mathrm{\φ}}_p|}\·(x-x_p)---\left(8\right)$

$\▿\mathrm{\φ}\left(x\right)=\▿{\mathrm{\φ}}_p,---\left(9\right)$

其中p是由导数质点携带的梯度。由于方程(8)是基于梯度方向的距离而不是欧氏距离，所以其产生比球隐函数更准确的结果。除以上修改之外，该转换程序与Enright等人[2002a]中提出的相同。

6 试验结果

本发明中提出的技术在具有双G5 2.5GHz处理器和5.5GB的存储器的Power Mac上实施。申请人使用程序来模拟真实世界中产生高雷诺数流体行为的多种实验情况。在实验中，申请人用以下常数来求解由水和空气组成的两相流体：g＝9.8m/sec²、ρ_water＝1000kg/m³、μ_water＝1.137×10³kg/ms、ρ_air＝1.226kg/m³、以及μ_air＝1.78×10⁵kg/ms，其中g是重力。使用移动立方体算法来进行水表面的提取，而且通过智能射

来进行渲染。

溃坝：为了比较导数质点模型与用质点水平集方法的线性模型之间的数值扩散，申请人用有效分辨率128²来执行2D溃坝测试，其中，在重力场释放0.2×0.4m的水柱。图10示出了结果，申请人能够看到导数质点产生扩散较少的结果：波的断裂更剧烈，而且很好地保存涡度。

落到浅水上的大块水：图11示出了从下述模拟截取的截图：即其中大块水落到0.05m深的蓄水池上。申请人对于蓄水池的底表面使用无滑动边界条件，这引起水在顶部快速地移动而在底部慢速地移动，导致皇冠状飞溅。这项实验用256³的有效网格分辨率来执行。还用使用质点水平集方法的传统线性模型来模拟相同的情景。比较证明所提出的技术产生更真实的结果。导数质点模型平均每时间步长用60秒，而线性模型平均每时间步长用30秒，均包括文件输出时间。

固体球的冲击：在这项实验中，具有半径0.15m的固体球以速度(5.0，-3.0，0.0)m/sec冲击到水中。该冲击产生图7中所示的复杂结构。这项实验表明所提出的技术能够产生猛烈的固体-水-空气相互作用中发生的详细的流体运动和表面特征。这项实验的有效分辨率是192×128×128。

旋转水箱：图12示出了具有半满的水的0.9×0.9×0.3m旋转箱。在实验中，离心力产生将水推向水箱侧的效果。这种模拟以96×96×32的有效网格分辨率来执行。这项实验证明即使在粗分辨率的网格中，导数质点模型也能够模拟流体的复杂运动。

洪水：在这项实验中，申请人模拟在6m×2m×4m域上发生的中等尺度的洪水。图13示出了从这项实验截取的截图。当移除堤坝时，水开始滑落台阶。由于沿路的阻碍，水进行产生复杂几何形状的猛烈运动。这项实验的有限分辨率是384×128×256。

图1～6示出了本发明的流程图。

如图1～6中所示，一种用于使用导数质点来模拟流体的详细运动的方法包括以下步骤：a)用基于八叉树数据结构的自适应网格建模流体(S100)；b)对于网格点处的流体速度求解不可压缩流体的纳维尔-斯托克斯方程(S200)；c)在纳维尔-斯托克斯方程的时间积分中使用分步法(S300)；d)将分步分组成非平流部分和平流部分(S400)；e)根据模拟的详细程度来选择网格平流部分或质点平流部分(S500)；f)对于网格平流部分使用基于八叉树的CIP方法(S600)；以及g)对于质点平流部分使用导数质点模型(S700)。

当模拟的详细程度高时使用对于质点平流部分使用导数质点模型的步骤(S700)。

对于高度感兴趣的区域的细节自适应网格增加网格的数目。

八叉树数据结构包括树形数据结构，每个内部结点具有多达八个孩子。

不可压缩流体的纳维尔-斯托克斯方程可以写成

$\frac{\∂u}{\∂t}=-(u\·\▿)u-\frac{\▿p}{\mathrm{\ρ}}+\▿\·(\mathrm{\μ}\▿u)+\frac{f}{\mathrm{\ρ}},$

$\▿\·u=0,$

其中u是速度，p是压力，f是外力，μ是粘度系数，ρ是流体密度。

该方法可以进一步包括使用水平集方法来追踪流体的表面的步骤(S800)。

如图1至3中所示，所述水平集方法包括水平集场φ。根据下式来更新所述水平集场φ：

$\frac{\∂\mathrm{\φ}}{\∂t}+u\·\▿\mathrm{\φ}=0$

符号距离 $|\▿\mathrm{\φ}|=1.$

$\frac{\∂\mathrm{\φ}}{\∂t}+u\·\▿\mathrm{\φ}=0$ 的时间积分仅涉及平流部分。网格平流部分平流输送根据当前速度场从非平流部分计算的速度和水平集场。

网格平流部分使用欧拉平流步骤来平流输送速度和水平集。通过基于八叉树的CIP(约束插值剖面)方法来执行所述欧拉平流。

对于网格平流部分使用基于八叉树的CIP方法的步骤(S600)包括用基于八叉树CIP方法来平流输送速度和水平集及其导数的步骤(S610)。直接从方程 $\frac{\∂\mathrm{\φ}}{\∂t}+u\·\▿\mathrm{\φ}=0$ 获得用于平流输送导数的方程，其中，φ是被平流输送的物理量，且通过对于x求方程的微分来获得空间变量x的导数的平流输送方程，即：

$\frac{\∂{\mathrm{\φ}}_x}{\∂t}+u\·\▿{\mathrm{\φ}}_x=-u_x\·\▿\mathrm{\φ}.$

如图2和图3中所示，求解x的导数的平流输送方程的步骤(S610)包括采用使用分步法的步骤，且所述使用分步法的步骤包括以下步骤：a)使用有限差分法来求解非平流部分 $\∂{\mathrm{\φ}}_x/\∂t=-u_x\·\▿\mathrm{\φ}$ (S620)；b)平流输送根据 $\frac{\∂{\mathrm{\φ}}_x}{\∂t}+u\·\▿{\mathrm{\φ}}_x=0$ 的结果(S630)；以及c)在由网格限定的单元的中心处进行压力值取样并在结点处进行包括速度、水平集及其导数的其它值进行取样(S640)。

如图3中所示，该方法可以进一步包括以下步骤：a)将基于八叉树的CIP法与自适应网格相结合(S650)；以及b)通过用基于八叉树的CIP法模拟平流来使由八叉树数据结构的网格分辨率中的区域性变化引起的八叉树人为瑕疵最小化(S660)。

导数质点模型使用基于质点的拉格朗日方案。用 $\mathrm{\φ}\left(x\right)={\mathrm{\φ}}_p+\frac{\▿{\mathrm{\φ}}_p}{|\▿{\mathrm{\φ}}_p|}\·(x-x_p)$ 和 $\▿\mathrm{\φ}\left(x\right)=\▿{\mathrm{\φ}}_p$ 来计算网格点处的水平集及其导数，其中p是由导数质点携带的梯度。

该方法可以进一步包括以下步骤：a)在对于质点平流部分使用导数质点模型的步骤(S700)之前将网格上限定的速度和水平集转换成质点上限定的速度和水平集(S690)(网格至质点)；以及b)在返回到对于网格平流部分使用基于八叉树的CIP法的步骤(S600)之前将质点上限定的速度和水平集转换成网格上限定的速度和水平集(S710)(质点至网格)。

转换速度和水平集的步骤(S690)网格至质点包括将具有用单调CIP插值代替线性插值的FLIP方法用于速度 $u^P=u_{-}^P+{\mathrm{\Σ}}_{i=1}^4{w}_i\mathrm{\Δ}{u}_i^G$ 和uP的转换的步骤(S692)，其中，

是就在非平流部分开始之前的质点速度，且

是由于非平流部分而引起的Gi处的速度变化。

如图5中所示，对于速度u^P、网格点G、选自G的每个象限的最近质点P₁、P₂、P₃和P₄、P_i的速度、以及

和

的x和y分量的导数，转换速度和水平集的步骤(S710)质点至网格包括以下步骤：a)将的导数

坐标变换成

因此

表示沿着方向

的分量且u_1⊥表示垂直于

的分量，其中，类似地获得

的坐标变换的导数

(S712)；b)沿着方向

对步骤(a)中获得的结果进行一维单调CIP插值以便计算位置A处的u_A及其

方向的导数

(S714)；c)通过应用以上坐标变换的反转来从

获得(u_Ax，u_Ay)(S716)；d)对质点P₃和P₄执行步骤(a)～(c)以计算B的u_B和(u_Bx，u_By)(S718)；以及e)对步骤(c)和(d)中获得的结果执行y方向单调CIP插值，其给出G处的速度和导数(S720)。

如图6中所示，对于质点平流部分使用导数质点模型的步骤包括以下步骤：a)调整质点的水平集值(S702)；以及b)执行质点重新播种(S704)。

7 结论

在本发明中，申请人已介绍了一种称为导数质点的新概念，并提出了增加平流部分的精度的流体模拟技术。使用欧拉方案来实施模拟器的非平流部分，而使用拉格朗日方案来实施平流部分。在开发以这种方式结合的模拟器时的重要问题是如何转换平流与非平流部分的切换处的物理量。申请人通过开发有效地抑制不必要的数值耗散的转换程序而成功地克服了这个问题。在不需要执行质点平流的流体区域中，使用基于八叉树的CIP解算器来模拟平流，这是申请人通过将CIP插值与八叉树数据结构相结合而在本工作中开发的。这种方法还用来降低数值耗散。多次实验的结果证明提出的技术显著降低了速度耗散，导致动态流体的真实再现。

没有先前的研究人员的开拓性工作，本工作将是不可能的。完全使用质点来求解平流部分的思想来自PIC[Harlow 1963]和FLIP[Brackbill和Ruppel 1986]方法，但是对于质点至网格和网格至质点速度换算，申请人使用基于导数的三次插值来代替线性插值。另外，让质点携带水平集值的想法来自质点水平集方法[Enright等人2002b]。这种模型在本工作中得到扩展，从而，质点还携带速度信息。此外，利用导数信息来自CIP方法[Yabe等人2001]。这里，申请人将这种方法的应用扩展至质点，引起质点至网格速度转换程序的开发，该转换程序包含CIP方法本身的扩展：非矩形构造中的质点的CIP插值。

在实验中，强调再现细节的能力。但是，提出的技术还能够应用于大尺度的流体运动以产生准确的结果。该技术包括增加的用于存储导数信息的存储器的量。然而，申请人证明能够在当前的PC上模拟2563网格，这已经是用于描绘感兴趣的流体情景的可用的分辨率。

虽然已参照不同的实施例示出并描述了本发明，但本领域的技术人员将认识到在不脱离权利要求书所限定的本发明的精神和范围的情况下可以进行形式、细节、组成以及操作上修改。

参考文献

BRACKBILL，J.U.，AND RUPPEL，H.M.1986.Flip：A method foradaptively zoned，particle-in-cell calculations of fluid flows in twodimensions.J.Comp.Phys.65，314-343.

CARLSON，M.，MUCHA，R.J.，AND TURK，G.2004.Rigid fluid：Animating the interplay between rigid bodies and fluid.ACM Transactionson Graphics 23，3，377-384.

ENRIGHT，D.，FEDKIW，R.，FERZIGER，J.，AND MITCHELL，I.2002.A hybrid particle level set method for improved interface capturing.J.Comp.Phys.183，83-116.

ENRIGHT，D.，MARSCHNER，S.，AND FEDKIW，R.2002.Animation and rendering of complex water surfaces.ACM Transactionson Graphics 21，3，736-744.

ENRIGHT，D.，LOSASSO，F.，AND FEDKIW，R.2005.A fast andaccurate semi-lagrangian particle level′set method.Computers andStructures 83，479-490.

FEDKIW，R.，STAM，J.，AND JENSEN，H.W.2001.Visualsimulation of smoke.Computer Graphics(Proc.ACM SIGGRAPH 2001)35，15-22.

FELDMAN，B.E.，O′BRIEN，J.F.，AND ARIKAN，0.2003.Animating suspended particle explosions.ACM Transactions on Graphics22，3，708-715.

FELDMAN，B.E.，O′BRIEN，J.F.，AND KLINGNER，B.M.2005.Animating gases with hybrid meshes.ACM Transactions on Graphics 24，3，904-909.

FELDMAN，B.E.，O′BRIEN，J.F.，KLINGNER，B.M.，ANDGOKTEKIN，T.G.2005.Fluids in deforming meshes.In Proceedings ofEurographics/ACM SIGGRAPH Symposium on Computer Animation2005，255-259.

FOSTER，N.，AND FEDKIW，R.2001.Practical animation of liquids.Computer Graphics(Proc.ACM SIGGRAPH 2001)35，23-30.

FOSTER，N.，AND METAXAS，D.1996.Realistic animation ofliquids.Graphical models and image processing：GMIP 58，5，471-483.

FOSTER，N.，AND METAXAS，D.1997.Controlling fluidanimation.In Computer Graphics International 97，178-188.

GOKTEKIN，T.G.，BARGTEIL，A.W.，AND O′BRIEN，J.F.2004.A method for animating viscoelastic fluids.ACMT ransactions onGraphics 23，3，463-468.

GUENDELMAN，E.，SELLE，A.，LOSASSO，F.，AND FEDKIW，R.2005.Coupling water and smoke to thin deformable and rigid shells.ACMTransactions on Graphics 24，3，973-981.

HARLOW，F.H.1963.The particle-in-cell method for numericalsolution of problems in fluid dynamics.In Experimental arithmetic，high-speed computations and mathematics.

HONG，J.-M.，AND KIM，C-H.2005.Discontinuous fluids.ACMTransactions on Graphics 24，3，915-920.

KANG，M.，FEDKIW，R.，AND LIU，X.-D.2000.A boundary-condition capturing method for multiphase incompressible flow.J.Sci.Comput.15，323-360.

LIU，X.-D.，FEDKIW，R.，AND KANG，M.2000.A boundarycondition capturing method for poisson′s equation on irregular domains.J.Comp.Phys.160，151-178.

LOSASSO，F.，GIBOU，F.，AND FEDKIW，R.2004.Simulatingwater and smoke with an octree data structure.ACM Transactions onGraphics 23，3，457-462.

MCNAMARA，A.，TREUILLE，A.，POPOVI′C，Z.，AND STAM，J.2004.Fluid control using the adjoint method.ACM Transactions onGraphics 23，3，449-456.

M＂ULLER，M.，CHARYPAR，D.，AND GROSS，M.2003.Particlebased fluid simulation for interactive applications.In Proceedingsof ACM SIGGRAPH/Eurographics Symposium on Computer Animation2003，154-159.

M＂ULLER，M.，SOLENTHALER，B.，KEISER，R.，AND GROSS，M.2005.Particle-based fluid-fluid interaction.In Proceedings ofEurographics/ACM SIGGRAPH Symposium on Computer Animation2005，237-244.

PARK，S.I.，AND KIM，M.J.2005.Vortex fluid for gaseousphenomena.In Proceedings of Eurographics/ACM SIGGRAPHSymposium on Computer Animation 2005，261-270.

PREMCTZ E，S.，TASDIZEN，T.，BIGLER，J.，LEFOHN，A.，ANDWHITAKER，R.T.2003.Particle-based simulation of fluids.InEurographics 2003 proceedings，Blackwell Publishers，401-410.

RASMUSSEN，N.，NGUYEN，D.Q.，GEIGER，W.，AND FEDKIW，R.2003.Smoke simulation for large scale phenomena.ACM Transactionson Graphics 22，3，703-707.

SELLE，A.，RASMUSSEN，N.，AND FEDKIW，R.2005.A vortexparticle method for smoke，water and explosions.ACM Transactions onGraphics 24，3，910-914.

SONG，O.-Y.，SHIN，H.，AND KO，H.-S.2005.Stable butnondissipative water.ACM Transactions on Graphics 24，1，81-97.

STAM，J.，AND FIUME，E.1995.Depicting fire and other gaseousphenomena using diffusion processes.Computer Graphics(Proc.ACMSIGGRAPH′95)29，Annual Conference Series，129-136.

STAM，J.1999.Stable fluids.Computer Graphics(Proc.ACMSIGGRAPH′99)33，Annual Conference Series，121-128.

STANIFORTH，A.，AND CλO TTE，J.1991.Semi-lagrangianintegration scheme for atmospheric model-a review.Mon.Weather Rev.119，12，2206-2223.

TREUILLE，A.，MCNAMARA，A.，POPOVI′C，Z.，AND STAM，J.2003.Keyframe control of smoke simulations.ACM Transactions onGraphics 22，3，716-723.

WANG，H.，MUCHA，P.J.，AND TURK，G.2005.Water drops onsurfaces.ACM Transactions on Graphics 24，3，921-929.

YABE，T.，AND AOKI，T.1991.A universal solver for hyperbolicequations by cubic-polynomial interpolation i.one-dimensional solver.Comp.Phys.Comm.66，219-232.

YABE，T.，AND AOKI，T.1991.A universal solver for hyperbolicequations by cubic-polynomial interpolation ii.two-and three dimensionalsolvers.Comp.Phys.Comm.66，233-242.

YABE，T.，XIAO，F.，AND UTSUMI，T.2001.The constrainedinterpolation profile method for multiphase analysis.J.Comp.Phys.169，556-593.

ZHU，Y.，AND BRIDSON，R.2005.Animating sand as a fluid.ACMTransactions on Graphics 24，3，965-972.

Claims (19)

1.一种使用导数质点来模拟流体的详细运动的方法，包括以下步骤：

a)用基于八叉树数据结构的自适应网格来建模流体；

b)对于网格点处的流体速度求解用于不可压缩流体的纳维尔-斯托克斯方程；

c)在纳维尔-斯托克斯方程的时间积分中使用分步法；

d)将各分步分组成非平流部分和平流部分；

e)根据模拟的详细程度来选择网格平流部分或质点平流部分；

f)对网格平流部分使用基于八叉树的CIP方法；以及

g)对质点平流部分使用导数质点模型，

其中，当模拟的详细程度高时，使用所述的对质点平流部分使用导数质点模型的步骤。

2.根据权利要求1所述的方法，其中，对于在细节上高度关注的区域，所述自适应网格增加网格的数目。

3.根据权利要求1所述的方法，其中，所述八叉树数据结构包括树形数据结构，其中，每个内部结点具有多达八个的孩子。

4.根据权利要求1所述的方法，其中，所述用于不可压缩流体的纳维尔-斯托克斯方程能够写成

$\frac{\∂u}{\∂t}=-(u\·\▿)u-\frac{\▿p}{\mathrm{\ρ}}+\▿\·(\mathrm{\μ}\▿u)+\frac{f}{\mathrm{\ρ}},$

$\▿\·u=0$

其中，u是速度，p是压力，f是外力，μ是粘度系数，ρ是流体密度。

5.根据权利要求1所述的方法，还包括使用水平集方法来追踪流体的表面的步骤，其中，所述水平集方法包括水平集场φ。

6.根据权利要求5所述的方法，其中，根据下式来更新所述水平集场φ|

$\frac{\∂\mathrm{\φ}}{\∂t}+u\·\▿\mathrm{\φ}=0$

带有符号距离 $|\▿\mathrm{\φ}|=1.$

7.根据权利要求6所述的方法，其中，所述 $\frac{\∂\mathrm{\φ}}{\∂t}+u\·\▿\mathrm{\φ}=0$ 的时间积分仅涉及平流部分。

8.根据权利要求6所述的方法，其中，所述网格平流部分根据当前的速度场来平流输送根据非平流部分计算的速度和水平集场。

9.根据权利要求6所述的方法，其中，所述网格平流部分使用欧拉平流步骤来平流输送速度和水平集。

10.根据权利要求9所述的方法，其中，通过基于八叉树的CIP(约束插值剖面)方法来执行所述欧拉平流。

11.根据权利要求1所述的方法，其中，所述的对网格平流部分使用基于八叉树的CIP方法的步骤包括用基于八叉树的CIP方法来平流输送速度和水平集和它们的导数的步骤。

12.根据权利要求11所述的方法，其中，直接从方程 $\frac{\∂\mathrm{\φ}}{\∂t}+u\·\▿\mathrm{\φ}=0$ 获得用于平流输送导数的方程，其中，φ是被平流输送的物理量，且通过关于x求方程的微分来获得空间变量x的导数的平流输送方程，即：

$\frac{{\∂\mathrm{\φ}}_x}{\∂t}+u\·{\▿\mathrm{\φ}}_x={-u}_x\·\▿\mathrm{\φ}.$

13.根据权利要求12所述的方法，其中，求解x的导数的平流输送方程的步骤包括使用分步法的步骤，其中，所述使用分步法的步骤包括以下步骤：

a)使用有限差分法来求解非平流部分 ${\∂\mathrm{\φ}}_x/\∂t=-u_x\·\▿\mathrm{\φ};$

b)平流输送根据 $\frac{{\∂\mathrm{\φ}}_x}{\∂t}+u\·\▿{\mathrm{\φ}}_x=0$ 的结果；以及

c)在由网格所限定单元的中心处对压力值取样，并在各结点处对包括速度、水平集和它们的导数的其它值取样。

14.根据权利要求1所述的方法，进一步包括步骤：

a)将基于八叉树的CIP方法与自适应网格相结合；以及

b)通过用基于八叉树的CIP方法模拟平流来最小化八叉树数据结构的网格分辨率中的区域性变化所引起的八叉树人为瑕疵。

15.根据权利要求1所述的方法，其中，所述导数质点模型使用基于质点的拉格朗日方案，其中，用下式来计算网格点处的水平集及其导数

$\mathrm{\φ}\left(x\right)={\mathrm{\φ}}_p+\frac{{\▿\mathrm{\φ}}_p}{\left({\▿\mathrm{\φ}}_p\right)}\·(x-x_p)$

以及

$\▿\mathrm{\φ}\left(x\right)={\▿\mathrm{\φ}}_p,$

其中，p是导数质点携带的梯度。

16.根据权利要求1所述的方法，进一步包括以下步骤：

a)在所述的对质点平流部分使用导数质点模型的步骤之前将各网格上定义的速度和水平集转换成各质点上定义的速度和水平集(网格到质点)；以及

b)在返回到所述的对网格平流部分使用基于八叉树的CIP方法的步骤之前将各质点上定义的速度和水平集转换成各网格上定义的速度和水平集(质点到网格)。

17.根据权利要求16所述的方法，其中，网格到质点的转换速度和水平集的步骤包括对于速度u^P的转换使用具有用单调CIP插值代替线性插值的FLIP方法的步骤，其中 $u^P=u_{-}^P+{\mathrm{\Σ}}_{i-1}^4{w}_i{\mathrm{\Δu}}_i^G,$ 其中是恰好在非平流部分开端之前的质点速度，且

是由于非平流部分引起的G_i处的速度变化。

18.根据权利要求16所述的方法，其中，对于速度u^P、网格点G、选自G的每个象限的最近质点P₁、P₂、P₃和P₄、Pi的速度、以及

和

的x分量和y分量的导数，

质点到网格的转换速度和水平集的步骤包括以下步骤：

a)将

的导数

坐标变换成(u_1‖u_1⊥)，使得u_1‖表示沿着方向

的分量且u_1⊥表示垂直于

的分量，其中，类似地获得

的坐标变换的导数(u_2‖，u_2⊥)；

b)沿着方向

对步骤(a)中获得的结果进行一维单调CIP插值以计算位置A处的u_A及其

方向的导数(u_A‖，u_A⊥)；

c)通过应用上述坐标变换的逆变换，从(u_A‖，u_A⊥)获得(u_Ax，u_Ay)；

d)对质点P₃和P₄执行步骤(a)～(c)以计算B的u_B和(u_Bx，u_By)，B为网格线和连接P₃和P₄的线段之间的交点；以及

e)对步骤(c)和(d)中获得的结果执行y方向单调CIP插值，给出G处的速度和导数。

19.根据权利要求15所述的方法，其中，所述的对质点平流部分使用导数质点模型的步骤包括以下步骤：

a)调整质点的水平集值；以及

b)执行质点重新播种。

Images