CN107273617B - 一种利用浅水波方程获得表面流流体运动的实时模拟方法及系统 - Google Patents
一种利用浅水波方程获得表面流流体运动的实时模拟方法及系统 Download PDFInfo
- Publication number
- CN107273617B CN107273617B CN201710468199.4A CN201710468199A CN107273617B CN 107273617 B CN107273617 B CN 107273617B CN 201710468199 A CN201710468199 A CN 201710468199A CN 107273617 B CN107273617 B CN 107273617B
- Authority
- CN
- China
- Prior art keywords
- calculating
- simulator
- model
- shallow water
- fluid
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 239000012530 fluid Substances 0.000 title claims abstract description 95
- 238000000034 method Methods 0.000 title claims abstract description 73
- 238000004088 simulation Methods 0.000 title claims abstract description 65
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 title claims abstract description 43
- 238000004364 calculation method Methods 0.000 claims abstract description 38
- 238000004090 dissolution Methods 0.000 claims abstract description 14
- 239000002245 particle Substances 0.000 claims description 33
- 239000007788 liquid Substances 0.000 claims description 30
- 239000011159 matrix material Substances 0.000 claims description 20
- 230000000694 effects Effects 0.000 claims description 19
- 239000007787 solid Substances 0.000 claims description 19
- 238000009499 grossing Methods 0.000 claims description 12
- 230000008859 change Effects 0.000 claims description 11
- 230000003993 interaction Effects 0.000 claims description 10
- 230000008569 process Effects 0.000 claims description 9
- PXFBZOLANLWPMH-UHFFFAOYSA-N 16-Epiaffinine Natural products C1C(C2=CC=CC=C2N2)=C2C(=O)CC2C(=CC)CN(C)C1C2CO PXFBZOLANLWPMH-UHFFFAOYSA-N 0.000 claims description 7
- 238000001556 precipitation Methods 0.000 claims description 6
- 230000002452 interceptive effect Effects 0.000 claims description 4
- 230000009466 transformation Effects 0.000 claims description 4
- 230000005855 radiation Effects 0.000 claims description 3
- 229940050561 matrix product Drugs 0.000 claims description 2
- 239000000203 mixture Substances 0.000 claims description 2
- 238000010586 diagram Methods 0.000 description 5
- 230000008901 benefit Effects 0.000 description 2
- 238000004422 calculation algorithm Methods 0.000 description 2
- 230000005012 migration Effects 0.000 description 2
- 238000013508 migration Methods 0.000 description 2
- 101100478173 Drosophila melanogaster spen gene Proteins 0.000 description 1
- 101100513476 Mus musculus Spen gene Proteins 0.000 description 1
- 230000001133 acceleration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000007664 blowing Methods 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 230000002706 hydrostatic effect Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 238000005461 lubrication Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000009877 rendering Methods 0.000 description 1
- 230000009897 systematic effect Effects 0.000 description 1
- 239000010409 thin film Substances 0.000 description 1
- 230000016776 visual perception Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- General Physics & Mathematics (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Computer Graphics (AREA)
- Software Systems (AREA)
- Processing Or Creating Images (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
一种表面流流体运动的实时模拟方法，通过使用浅水波方程，对其中的不同项的计算提出创新性的计算方法，获得系统性，高效率，能够处理非可忽略速度流体的表面流流体解决方案。其中，本方法首次提出了基于特征的粘滞模型与可处理表面溶解的表面多组分流体运动模型。本方法对于浅水波方程中移流、表面张力、与3D模拟器结合等问题的处理方式相比前人的方法也更加高效。本方法可以完全并行在GPU上，获得实时的模拟速度。本发明具有较强的实用性，能够为表面流流体的模拟提供完整全面的高效解决方案。
Description
技术领域
本发明涉及计算机图形学流体模拟与渲染技术领域，尤其涉及一种流体模拟方法与系统。
背景技术
流体现象对于真实世界的视觉感受而言非常重要。计算机图形学领域为了重现丰富多彩的流体现象，开发了多种多样的流体模拟器。
在之前的工作中，大多数的关注重点在于如何在较少障碍物的三维空间中模拟流体。而生活中最常见的流体，如空气，水等具有较为透明的特征，其多样化的运动现象在固体表面边界附近时才在视觉上更加明显。这即是一类表面流流体问题的研究对象：即对于任意的三维模型，研究一薄层流体在模型表面的运动。
图形学中存在一些研究流体在二维几何流型中运动的工作，如Z.Fan等人2005年发表“Adapted unstructured lbm for flow simulation on curved surfaces”,S.Auer等人2012年发表“Real-time Fluid effects on surfaces using the closest pointmethod”,分别使用LBM方法和最近点方法将二维平面的零厚度流体模拟扩展到二维流型中，获得类似三维流体模拟的投影效果。而对于在几何模型表面上存在的有限厚度的薄层流体的模拟，H.Wang等人在2007年发表”Solving general shallow wave equations onsurfaces”，O.Azencot等人在2015年发表“Functional thin films on surfaces”,分别使用浅水波方程和基于润滑理论的优化方法获得了一定的表面流效果。然而，他们的效果只能适用于以可忽略速度运动的高粘滞性流体，并且无法在任意模型表面上做到高度并行，因而速度较慢，无法实时化。
在任意几何模型表面模拟流体还涉及到较为复杂的数值计算技巧。L.Shi等人在2004年发表“Inviscid and incompressible fluid simulation on triangle meshes”，S.Jeong等人在2013年发表“Combustion waves on the point set surface”,研究几何模型表面的移流计算的处理方式。然而前者的计算较为复杂，后者只涉及标量场的计算，而表面流问题需要处理速度矢量场的移流。上述H.Wang等人的文章提出了表面流中一种表面张力的计算方式，然而其亦存在近似精度低，假设较强等问题。N.Thurey等人在2006年发表“Animation of open water phenomena with coupled shallow water and freesurface simulations”,N.Chentanez等人在2015年发表“Coupleing 3d eulerian,heightfield and particle methods for interactive simulation of large scaleliquid phenomena”,将浅水波模拟器与3D流体模拟器相结合，然而他们仅能处理向上表面的互相结合问题。近年来，B.Ren等人在2014年发表“Multiple-fluid sph simulationusing a mixture model”，研究三维空间中的可混多组分流体运动，但目前图形学的相关研究中尚不涉及表面流中的多组分流体运动计算方法。
发明内容
本发明的目的是解决在表面流流体的模拟中，如何高效地处理方程中的复杂的数值计算问题，提供一种表面流流体运动的实时模拟方法，以便高效地恢复真实的表面流运动，并与三维模拟器结合，恢复包含多组分流体运动在内的多样化的表面流流体现象，做到任意模型上的实时表面流流体模拟。
本发明的技术方案：
一种利用浅水波方程获得表面流流体运动的实时模拟方法，包括：
第1:利用一种基于特征的粘滞力模型，对几何模型进行预计算处理；
包括对几何模型进行平滑化并计算平滑化三角网格的外法向，在平滑化三角网格的格点上计算控制面积，粘滞力张量矩阵与局部切平面坐标系间仿射变换矩阵；
第2：在模拟开始时载入预计算的数据；
对于可形变模型表面上的模拟，在每个时间步开始时重新计算平滑化三角网格的外法向，仿射变换矩阵与粘滞力张量矩阵；
所述的可形变模型表面上的模拟及重新计算方法进一步包括：
在可形变模型表面上的模拟中，在平滑时调整平滑网格顶点的位置，使得原网格的顶点总在平滑网格顶点沿其法向方向的直线上；之后，通过移动平滑网格顶点的位置或编辑底面高度控制原模型的形变；
第3:使用单组分浅水波方程计算单组分的表面流流体运动，或使用多组分浅水波方程计算多组分的表面流流体运动；
第4:选择性包括与3D模拟器进行交互；
对于包含交互的情况，在第1步中的预计算时根据模型生成粒子3D模拟器中的粒子边界；同时将3D模拟器进行一个时间步的计算，并在时间步的最后，计算浅水波方程中三角网格上的模拟器以及粒子3D模拟器间的物理量交互。
其中，
第1步所述预计算进一步包括：
第1.1、将原几何模型平滑化；
第1.2、对每个平滑化网格顶点，计算邻域内的一个半平整化模型(SOM)；
所述SOM的计算方法包括：
第1.2.1、对每个所考虑的平滑模型上的格点v，将其周围邻域内的点沿该考虑格点的法向投影至切平面；
第1.2.2、v沿其法向对应一个原模型上的点vp,对于任意原模型上的点vp',同样可沿同一方向获得平滑模型上的投影对应点v’；
第1.2.3、v，v’,vp,vp'一定能够构成一个平面，与v的切平面交于一直线，vv’与该直线的夹角为θ，则将vp'以vp为中心在该平面内旋转θ角度；
第1.2.4、对所有v附近的vp'作上述三步，获得v邻域内的一个“半平整化模型”(SOM)；
第1.3、在SOM上计算每个三角形的形状运算符S(Shape Operator)；
第1.4、利用第1.3步中的计算结果获得每个平滑网格顶点切平面内的平均粗糙程度最大的方向和与其垂直的方向，以及两个方向上的粗糙系数；
第1.5、利用第1.4步结果获得粘滞力张量矩阵，粘滞力张量矩阵与速度矢量的矩阵乘积结果即为表面流流体运动中所受的模型产生的粘滞力。
第3步所述的计算单组分的表面流流体运动的方法包括：
第3.1、在每个平滑网格点上存储表面流流体的速度场与深度场；
第3.2、在每个时间步开始时采用半拉格朗日法对速度场与深度场进行移流计算；
所述的在每个时间步开始时采用半拉格朗日法对速度场与深度场进行移流计算包括：
第3.2.1、在考虑顶点i时，对于所有周围的顶点j，设其法向分别为ni、nj，将顶点j的3D坐标沿方向ni×nj旋转，使ni与旋转后的nj方向相同；
第3.2.2、此时两点i、j的切平面平行，计算两点局部切平面坐标系间的放射矩阵Mij；
第3.2.3、在半拉格朗日法进行速度插值时，即能够直接使用矩阵Mij和顶点j的速度uj插值获得顶点i的速度ui；
第3.2.4、进一步对深度场进行插值，方法是直接对原深度场进行插值；
第3.3、计算流体所受粘滞力；
第3.4、计算流体所受表面张力压强；
所述计算流体所受表面张力压强的方法是：
第3.4.1、判断网格点处于液体内部还是边缘；
第3.4.2、对于液体内部的格点，利用表面能的梯度获得表面张力，进而计算出表面张力压强；
第3.4.3、对于液体边缘的格点，利用表面能沿表面切平面方向的梯度，根据接触角数据，计算出边缘的表面张力，进而计算出表面张力压强；
第3.5、在网格上计算速度的散度与流体压强的梯度；
第3.6、根据浅水波方程获得速度场与深度场的变化率。
对于固体组分和液体组分的混合组分，在第3.5步和第3.6步之间还应当包括：
(1)：计算液体与固体间的粘滞作用力；之后：
(2)：计算固体的溶解与析出及其对流体运动造成的影响；
所述计算固体的溶解与析出及其对流体运动造成影响的方法是：
A：根据用户指定的溶解度和溶解速率计算固体当前时间步的溶解量，或者当超过溶解度时的析出量；
B：根据步骤A中的数据计算液体因此损失或获得的动量，加入到液体速度变化率的计算中。
最后根据浅水波方程获得固体与液体两组分的速度场与深度场的变化率。
第4步所述与3D模拟器进行交互的方法是：
第4.1、在模拟开始前根据模型生成3D模拟器的粒子边界；
第4.2、同时将粒子3D模拟器进行一个时间步的计算；
第4.3、根据用户给定的阈值计算浅水波模拟器应向粒子3D模拟器传送的液体质量、位置和速度，并生成3D模拟器中的粒子；
第4.4、判断3D模拟器中运动轨迹与浅水波模拟器模型表面相交的粒子，根据粒子质量计算浅水波方程中的体积源项，影响下一时间步的浅水波方程模拟器计算，同时在3D模拟器中删除该粒子。
本发明同时还提供了一种利用浅水波方程获得表面流流体运动的实时模拟系统，该系统包括：
输入/输出模块：用于模拟过程中数据的输入输出，包括加载模拟过程中所需的必要物理参数与用户指定的参数，保存流体模拟器所获得的图片与数据结果等；
预计算模块：用于预计算数据与需要对预计算数据重新计算情况下的重新计算；
载入与更新模块：用于载入预计算数据，及在每时间步开始时需要对预计算数据重新计算情况下的重新计算；
移流计算模块：用于计算3D模型表面上的单组分或多组分流体的移流；
表面张力计算模块：用于计算单组分或多组分的表面流流体的表面张力；
流体速度与深度更新模块：用于计算单组分或多组分流体的表面流流体速度与深度的变化率；
与3D模拟器交互模块：用于处理与3D粒子模拟器的交互模拟。
本发明的优点和有益效果：
采用本发明方法的主要优势在于提供了一个系统性，高效率，能够处理非可忽略速度流体的表面流流体解决方案。其中，本方法首次提出了基于特征的粘滞模型与可处理表面溶解的表面多组分流体运动模型。本方法对于浅水波方程中移流、表面张力、与3D模拟器结合等问题的处理方式相比前人的方法也更加高效。本方法可以完全并行在GPU上，获得实时的模拟速度。
附图说明
图1是本发明算法流程图；
图2是复杂表面上(t＝0.4、2.7、4.5、6.4、9.0s时)沿物体特征的表面流运动效果；
图3是可形变表面的(t＝0.7、3.5、4.4、7.9、10.0s时)表面流运动效果；
图4是由艺术家根据表面颜色梯度人为指定T获得的(t＝2.0、4.0s时)艺术性效果；
图5是在上述基础上(t＝2.0、4.0s时)多组分表面流模拟效果示意图；
图6是(t＝1.8、6.5、8.3、12.3s时)与3D模拟器结合的模拟效果；
图7是(t＝4.6、6.1s时)不同风速环境下的模拟效果；
图8是获得顶点v邻域内一个“半平整化模型”(SOM)示意图；
图9是控制面积示意图，即连接各三角形质心的阴影区域。
图10是利用浅水波方程获得表面流流体运动的实时模拟系统框图。
具体实施方式
一、如图10所示利用浅水波方程获得表面流流体运动的实时模拟系统，该系统包括：
输入/输出模块：用于模拟过程中数据的输入输出，包括加载模拟过程中所需的必要物理参数与用户指定的参数，保存流体模拟器所获得的图片与数据结果等；
预计算模块：用于预计算数据与需要对预计算数据重新计算情况下的重新计算；
移流计算模块：用于计算3D模型表面上的单组分或多组分流体的移流；
表面张力计算模块：用于计算单组分或多组分的表面流流体的表面张力；
流体速度与深度更新模块：用于计算单组分或多组分流体的表面流流体速度与深度的变化率；
与3D模拟器交互模块：用于处理与3D粒子模拟器的交互模拟。
二、利用浅水波方程获得表面流流体运动的实时模拟方法
本发明实时模拟方法的算法流程(参见图1)步骤具体如下：
A：利用一种基于特征的粘滞力模型，预计算所有可以预计算的数据并存储。
其中，粘滞力张量可以由艺术家人为指定获得想要的艺术效果。对于存在与3D模拟器交互的情况，根据模型生成3D模拟器中的粒子边界。
在本发明方法中，为了成功模拟出真实世界中流体会沿着表面沟壑等细节特征流动的现象，提出了一种基于特征的粘滞力模型如下：
首先，将原模型平滑化，本方法的后续模拟将进行在该平滑后的模型上。为了恢复平滑过程中损失的原模型细节信息，在每个平滑模型的格点上将计算一个粘滞力张量T。方法如下：
1.对每个所考虑的平滑模型上的格点v，将其周围邻域内的点沿该考虑格点的法向投影至切平面；
2.v沿其法向对应一个原模型上的点vp,对于任意原模型上的点vp',同样可沿同一方向获得平滑模型上的投影对应点v’；
3.v，v’,vp,vp'一定可以构成一个平面，与v的切平面交于一直线，vv’与该直线的夹角为θ，则将vp′以vp为中心在该平面内旋转θ角度；
4.对所有v附近的vp′作上述三步，获得v邻域内的一个“半平整化模型”(SOM)。如图8所示。
5.在该SOM上计算每个三角面片的形状运算符(Shape Operator)S；
6.对每个v，计算如下等式：
其中，投影相交集Iv由所有沿v法向方向投影后，与v在平滑表面上控制面积相交的SOM上的三角面片组成。A与S为三角面片的面积与形状运算符。x为v切平面上的方向矢量。该等式的意义在于找到v切平面上“平均粗糙程度最大”的方向与该方向上的粗糙系数τmax。控制面积如图9所示，为连接各三角形质心的阴影区域。
7.计算与垂直的方向并计算下式：
8.设其中e为局部切平面坐标系中的坐标轴矢量，M为仿射矩阵。则T＝μMdiag{τmax,τ⊥}M-1
其中μ为粘滞系数。为计算方便，在实现中用户可直接指定摩擦系数γ＝μ/ρ用于代入(2)式后的计算。
9.最后，Γb＝T·u
以上基于特征的粘滞力模型将可以获得自然的流体沿表面特征运动的模拟效果。
在上面所述的模拟方法中，对于非形变模型上的模拟，原始模型的平滑化，各类仿射矩阵计算，各顶点的控制面积与粘滞力张量T均使用一个预计算步骤计算并存储，在之后复数的模拟过程中可以直接调用，无需重复计算。
B：在模拟开始时载入上述预计算的数据，对可形变模型，在每个时间步开始时重新计算必须更新的数据。具体而言，我们在预计算的平滑化时调整平滑网格顶点的位置，使得原网格的顶点总在平滑网格顶点沿其法向方向的直线上。之后，我们通过移动平滑网格顶点的位置或编辑底面高度控制原模型的形变。这样，仅需重新计算平滑网格的法向，仿射矩阵与粘滞力张量。我们的方法采用后面这种较为高效的手段模拟可形变模型上的表面流运动。
C：采用浅水波公式计算单组分或多组分的表面流流体运动，具体项(移流，表面张力，粘滞力等)的技术实现方法具体说明如下。
在本发明中，采用浅水波方程作为基本模拟方法，其基本方程为：
在上述两式中，第(1)式为表面流流体沿表面法向方向的深度的变化方程(质量方程)，第(2)式为流体的速度变化方程(运动方程)。其中d为流体沿表面法向方向的深度，u为速度，P为压强，由静水压ρgh与表面张力压强Ps相加而得，ρ为密度，Γ＝Γw-Γb为风带来的粘滞力Γw与底面带来的粘滞力Γb之差，aext为外力等导致的加速度，Q为体积源。
为风所带来的粘滞力计算公式，γw为用户设定的系数，uw为风速。
对于(1)(2)两式的移流项，本方法中采用半拉格朗日法处理。对于矢量场的移流快速插值，本方法提出使用如下方式：
(1)在考虑顶点i时，对于所有周围的顶点j，设其法向分别为ni，nj，将顶点j的3D坐标沿方向ni×nj旋转，使ni与旋转后的nj方向相同。
(2)此时两点i，j的切平面平行，可以计算两点局部切平面坐标系间的放射矩阵Mij。
(3)在半拉格朗日法进行速度插值时，即可直接使用矩阵Mij和顶点j的速度uj插值获得顶点i的速度ui。
对于顶点v处表面张力的计算，我们基于表面能观点使用如下等式：
求和对所有包括所考虑的顶点v的流体表面的三角面片进行，式中的梯度表示仅对v的坐标进行的偏导算子。σ为表面张力系数。
对液体与固体之间的接触角，设其为θc，计算
σlg＝σ,
则
其中三个求和分别对所有包含v的液体表面，液体-固体交面，固体-气体交面进行。其中的梯度算子代表沿v的切平面的二维偏导。
最后，对液体内部的顶点，使用作为表面张力压强，对液体边界上的点，使用作为表面张力压强。其中面积A的垂直角标分别代表垂直于液体表面的投影面积与垂直与v的切平面上液体边界线法向的投影面积。该表面张力压强的正负号可根据(5)(6)两式左边的力矢量与模型外法向是否同向确定：同向为负，异向为正。
本方法还提供了表面多组分流体运动的模拟方式，即液体中溶解固体的模拟，公式如下：
角标l,s分别代表液体、固体的物理量。
其中Γl＝Γbl+Γdrag，由以下两式计算：
Γbl＝λγlT·ul-κ(λ-1)γsT
Γdrag＝κ(λ-1)Cd(ul-us)|ul-us|
其中
而
定义当前溶解量Qs＝ep(dl+ds)|ul|，则液体因固体溶解的动量损失
其中ep为用户设置的溶解速率。
在模拟中，本方法设置一个溶解率上限，在模拟时检查ds与dl的比值不超过该上限。对于ds超过的部分，将其从中减去并黏附回模型表面上使得底面高度b相应增加同样的量。
D：对于包含3D模拟器交互的情况，同时将3D模拟器进行一个时间步，在时间步的最后，处理两个模拟器间的物理量交互。
本方法中，利用体积源Q处理与3D模拟器的结合问题。当某点i的液体深度d大于用户指定的∈1时，我们将d设为用户指定的∈2，并将总计体积为(d-∈2)Ai的液体粒子加入SPH三维模拟器，其中Ai为i的控制面积。该SPH粒子的初始位置为其中vi是i的世界位置坐标，b是浅水波模型中的底面高度，n为该处的外法向。该粒子的初速为其中分子上为当前时间步计算出的液体深度变化。于此同时，我们还搜索i点直接相连的顶点，若其深度大于用户指定的(稍小的)值η1，用与上面同样的后续步骤将值设置为用户指定的η2。
反过来，当一个SPH粒子运动至模型表面上时，我们认为其被表面流体吸收。此时在i点(1)式中，Q在一个时间步内临时增加同时删除该SPH粒子。
图2至图7分别给出了采用本发明方法进行实时模拟的效果，其中各参数均被设置如下：
γ＝1.3,σ＝1.0,κ＝2.0,Cd＝10.0,ep＝0.005
1＝3.0,2＝1.6,η1＝2.4,η2＝2.1。
图2是复杂表面上(t＝0.4、2.7、4.5、6.4、9.0s时)沿物体特征的表面流运动效果；所使用的模型顶点数约44000，模拟速度为每帧7.44ms。
图3是可形变表面(t＝0.7、3.5、4.4、7.9、10.0s时)的表面流运动效果；模型定点数23200，模拟速度41.36ms/帧，其中用于形变表面信息重新计算的时间为28.92ms/帧。
图4是由艺术家根据表面颜色梯度人为指定T获得的(t＝2.0、4.0s时)艺术性效果；其中底面颜色的梯度方向被指定为最大粗糙程度的方向，相应粗糙系数数值为颜色梯度值，且该粗糙系数被指定为5倍于垂直方向的粗糙系数；摩擦系数被指定为80网格顶点数52200，速度9.94ms/帧。
图5是在上述基础上(t＝2.0、4.0s时)多组分表面流模拟效果示意图；其中，固液密度比为2：1，速度15.24ms/帧。
图6是(t＝1.8、6.5、8.3、12.3s时)与3D模拟器结合的模拟效果；模型包含23200个顶点，粒子模拟器中共有129000个粒子。模拟时间为使用浅水波的三角网格上模拟器20.52ms/帧，粒子模拟器26.73ms/帧，合计47.25ms/帧。
图7是(t＝4.6、6.1s时)不同风速环境下的模拟效果；图中由上到下为：无风；风迎面吹来；风由左手方向吹来的结果；使用的风速场由另外的现有三维网格气体模拟器事先生成。
Claims (10)
1.一种利用浅水波方程获得表面流流体运动的实时模拟方法，其特征在于该方法包括：
第1:利用一种基于特征的粘滞力模型，对几何模型进行预计算处理；
包括对几何模型进行平滑化并计算平滑化三角网格的外法向，在平滑化三角网格的格点上计算控制面积，粘滞力张量矩阵与局部切平面坐标系间仿射变换矩阵；
第2：在模拟开始时载入预计算的数据；
对于可形变模型表面上的模拟，在每个时间步开始时重新计算平滑化三角网格的外法向，仿射变换矩阵与粘滞力张量矩阵；
第3:使用单组分浅水波方程计算单组分的表面流流体运动，或使用多组分浅水波方程计算多组分的表面流流体运动；
第4:选择性包括与3D模拟器进行交互；
对于包含交互的情况，在第1步中的预计算时根据模型生成粒子3D模拟器中的粒子边界；同时将3D模拟器进行一个时间步的计算，并在时间步的最后，计算浅水波方程中三角网格上的模拟器以及粒子3D模拟器间的物理量交互。
2.根据权利要求1所述的利用浅水波方程获得表面流流体运动的实时模拟方法，其特征在于第1步所述预计算进一步包括：
第1.1、将原几何模型平滑化；
第1.2、对每个平滑化网格顶点，计算邻域内的一个半平整化模型SOM；
所述SOM的计算方法包括：
第1.2.1、对每个所考虑的平滑模型上的格点v，将其周围邻域内的点沿该考虑格点的法向投影至切平面；
第1.2.2、v沿其法向对应一个原模型上的点vp,对于任意原模型上的点vp',同样能够沿同一方向获得平滑模型上的投影对应点v’；
第1.2.3、v，v’,vp,vp'一定能够构成一个平面，与v的切平面交于一直线，vv’与该直线的夹角为θ，则将vp'以vp为中心在该平面内旋转θ角度；
第1.2.4、对所有v附近的vp'作上述三步，获得v邻域内的一个“半平整化模型”SOM；
第1.3、在SOM上计算每个三角形的形状运算符S(Shape Operator)；
第1.4、利用第1.3步中的计算结果获得每个平滑网格顶点切平面内的平均粗糙程度最大的方向和与其垂直的方向，以及两个方向上的粗糙系数；
第1.5、利用第1.4步结果获得粘滞力张量矩阵，粘滞力张量矩阵与速度矢量的矩阵乘积结果即为表面流流体运动中所受的模型产生的粘滞力。
3.根据权利要求2所述的利用浅水波方程获得表面流流体运动的实时模拟方法，其特征在于，所述的可形变模型表面上的模拟及重新计算方法进一步包括：
在可形变模型表面上的模拟中，在平滑时调整平滑网格顶点的位置，使得原网格的顶点总在平滑网格顶点沿其法向方向的直线上；之后，通过移动平滑网格顶点的位置或编辑底面高度控制原模型的形变。
4.根据权利要求1所述的利用浅水波方程获得表面流流体运动的实时模拟方法，其特征在于第3步所述的计算单组分的表面流流体运动的方法包括：
第3.1、在每个平滑网格点上存储表面流流体的速度场与深度场；
第3.2、在每个时间步开始时采用半拉格朗日法对速度场与深度场进行移流计算；
第3.3、计算流体所受粘滞力；
第3.4、计算流体所受表面张力压强；
第3.5、在网格上计算速度的散度与流体压强的梯度；
第3.6、根据浅水波方程获得速度场与深度场的变化率。
5.根据权利要求4所述的利用浅水波方程获得表面流流体运动的实时模拟方法，其特征在于，对于固体组分和液体组分的混合组分，在第3.5步和第3.6步之间还应当包括：
(1)：计算液体与固体间的粘滞作用力；之后：
(2)：计算固体的溶解与析出及其对流体运动造成的影响；
最后根据浅水波方程获得固体与液体两组分的速度场与深度场的变化率。
6.根据权利要求4或5所述的利用浅水波方程获得表面流流体运动的实时模拟方法，其特征在于，第3.2步所述的在每个时间步开始时采用半拉格朗日法对速度场与深度场进行移流计算包括：
第3.2.1、在考虑顶点i时，对于所有周围的顶点j，设其法向分别为ni、nj，将顶点j的3D坐标沿方向ni×nj旋转，使ni与旋转后的nj方向相同；
第3.2.2、此时两点i、j的切平面平行，计算两点局部切平面坐标系间的放射矩阵Mij；
第3.2.3、在半拉格朗日法进行速度插值时，即能够直接使用矩阵Mij和顶点j的速度uj插值获得顶点i的速度ui；
第3.2.4、进一步对深度场进行插值，方法是直接对原深度场进行插值。
7.根据权利要求4或5所述的利用浅水波方程获得表面流流体运动的实时模拟方法，其特征在于，第3.4步所述计算流体所受表面张力压强的方法是：
第3.4.1、判断网格点处于液体内部还是边缘；
第3.4.2、对于液体内部的格点，利用表面能的梯度获得表面张力，进而计算出表面张力压强；
第3.4.3、对于液体边缘的格点，利用表面能沿表面切平面方向的梯度，根据接触角数据，计算出边缘的表面张力，进而计算出表面张力压强。
8.根据权利要求5所述的利用浅水波方程获得表面流流体运动的实时模拟方法，其特征在于，所述计算固体的溶解与析出及其对流体运动造成影响的方法是：
A：根据用户指定的溶解度和溶解速率计算固体当前时间步的溶解量，或者当超过溶解度时的析出量；
B：根据步骤A中的数据计算液体因此损失或获得的动量，加入到液体速度变化率的计算中。
9.根据权利要求1所述的利用浅水波方程获得表面流流体运动的实时模拟方法，其特征在于，第4步所述与3D模拟器进行交互的方法是：
第4.1、在模拟开始前根据模型生成3D模拟器的粒子边界；
第4.2、同时将粒子3D模拟器进行一个时间步的计算；
第4.3、根据用户给定的阈值计算浅水波模拟器应向粒子3D模拟器传送的液体质量、位置和速度，并生成3D模拟器中的粒子；
第4.4、判断3D模拟器中运动轨迹与浅水波模拟器模型表面相交的粒子，根据粒子质量计算浅水波方程中的体积源项，影响下一时间步的浅水波方程模拟器计算，同时在3D模拟器中删除该粒子。
10.一种利用浅水波方程获得表面流流体运动的实时模拟系统，其特征在于该系统包括：
输入/输出模块：用于模拟过程中数据的输入输出，包括加载模拟过程中所需的物理参数与用户指定的参数，保存流体模拟器所获得的图片与数据结果；
预计算模块：用于预计算数据；
载入与更新模块：用于载入预计算数据，及在每时间步开始时需要对预计算数据重新计算情况下的重新计算；
移流计算模块：用于计算3D模型表面上的单组分或多组分流体的移流；
表面张力计算模块：用于计算单组分或多组分的表面流流体的表面张力；
流体速度与深度更新模块：用于计算单组分或多组分流体的表面流流体速度与深度的变化率；
与3D模拟器交互模块：用于处理与3D粒子模拟器的交互模拟。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710468199.4A CN107273617B (zh) | 2017-06-20 | 2017-06-20 | 一种利用浅水波方程获得表面流流体运动的实时模拟方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710468199.4A CN107273617B (zh) | 2017-06-20 | 2017-06-20 | 一种利用浅水波方程获得表面流流体运动的实时模拟方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107273617A CN107273617A (zh) | 2017-10-20 |
| CN107273617B true CN107273617B (zh) | 2019-07-26 |
Family
ID=60069517
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710468199.4A Active CN107273617B (zh) | 2017-06-20 | 2017-06-20 | 一种利用浅水波方程获得表面流流体运动的实时模拟方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107273617B (zh) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109840935A (zh) * | 2017-12-12 | 2019-06-04 | 中国科学院计算技术研究所 | 基于深度采集设备的波浪重建方法和系统 |
| CN108416838B (zh) * | 2018-03-02 | 2022-05-31 | 南开大学 | 一种利用相场理论进行二维与三维晶体生长模拟的方法及系统 |
| CN108229083A (zh) * | 2018-04-11 | 2018-06-29 | 南京航空航天大学 | 一种基于改进的有限差分格式的水流数值模拟方法 |
| CN111724225A (zh) * | 2019-03-21 | 2020-09-29 | 北京京东尚科信息技术有限公司 | 液体商品展示方法、终端设备、存储介质及电子设备 |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2003341315A (ja) * | 2002-05-28 | 2003-12-03 | Sumitomo Rubber Ind Ltd | タイヤの排水性能のシミュレーション方法、シミュレーション装置、及びシミュレーションプログラムを記録した記録媒体 |
| CN101930622B (zh) * | 2009-09-29 | 2012-03-28 | 北京航空航天大学 | 浅水水波的真实感建模与绘制 |
| CN103714575B (zh) * | 2013-12-30 | 2016-09-07 | 北京大学 | 一种sph与动态表面网格相结合的流体仿真方法 |
| CN104200015B (zh) * | 2014-08-20 | 2017-06-16 | 清华大学 | 一种流体模拟方法及装置 |
| CN106570305B (zh) * | 2015-10-09 | 2019-08-09 | 清华大学 | 基于亥姆霍兹自由能的多组分流体模拟方法及装置 |
-
2017
- 2017-06-20 CN CN201710468199.4A patent/CN107273617B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN107273617A (zh) | 2017-10-20 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107273617B (zh) | 一种利用浅水波方程获得表面流流体运动的实时模拟方法及系统 | |
| US7479963B2 (en) | Method and system for performing computer graphic simulation of a fluid using target-driven control | |
| US5802353A (en) | Haptic computer modeling system | |
| Shi et al. | Taming liquids for rapidly changing targets | |
| US7983884B2 (en) | Water particle manipulation | |
| US6515668B1 (en) | Computer graphics animation method and device | |
| Raveendran et al. | Controlling liquids using meshes | |
| US10467791B2 (en) | Motion edit method and apparatus for articulated object | |
| KR100568563B1 (ko) | 입자동역학 해석기법과 볼륨렌더링 기법을 이용한 실시간 유체유동 시뮬레이션 및 렌더링 방법 | |
| KR100568562B1 (ko) | 연속체 유체역학 해석기법과 볼륨렌더링 기법을 이용한 실시간 유체유동 시뮬레이션 및 렌더링 방법 | |
| Shi et al. | Inviscid and incompressible fluid simulation on triangle meshes | |
| Madill et al. | Target particle control of smoke simulation | |
| Kim et al. | Fire sprite animation using fire-flake texture and artificial motion blur | |
| Kim et al. | Efficient representation of detailed foam waves by incorporating projective space | |
| CN108416838B (zh) | 一种利用相场理论进行二维与三维晶体生长模拟的方法及系统 | |
| Thalmann et al. | The Making of the Xian terra-cotta Soldiers | |
| KR100705417B1 (ko) | 물체 표면의 젖음과 마름 표현 장치 및 방법 | |
| Di Fiore et al. | A framework for user control on stylized animation of gaseous phenomena | |
| Stiver et al. | Sketch based volumetric clouds | |
| Ai et al. | An improved shallow water equation model for water animation | |
| Pukki | Unreal Engine 4 for Automation | |
| CN104463933A (zh) | 一种基于三视图的2.5维卡通动画自动生成方法 | |
| Shi et al. | Object modeling and animation with smoke | |
| Wyvill et al. | Using soft objects in computer-generated character animation | |
| Chayanurak et al. | Ice Melting Simulation using SPH and Heat Transfer with Constant Ambient Temperature |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |