CN110717269A - 一种基于网格和粒子耦合的流体表面细节保护方法 - Google Patents
一种基于网格和粒子耦合的流体表面细节保护方法 Download PDFInfo
- Publication number
- CN110717269A CN110717269A CN201910956536.3A CN201910956536A CN110717269A CN 110717269 A CN110717269 A CN 110717269A CN 201910956536 A CN201910956536 A CN 201910956536A CN 110717269 A CN110717269 A CN 110717269A
- Authority
- CN
- China
- Prior art keywords
- grid
- particle
- fluid
- particles
- mesh
- 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.)
- Granted
Links
- 239000002245 particle Substances 0.000 title claims abstract description 181
- 239000012530 fluid Substances 0.000 title claims abstract description 94
- 238000000034 method Methods 0.000 title claims abstract description 82
- 238000010168 coupling process Methods 0.000 title claims abstract description 37
- 230000008878 coupling Effects 0.000 title claims abstract description 33
- 238000005859 coupling reaction Methods 0.000 title claims abstract description 33
- 230000002159 abnormal effect Effects 0.000 claims abstract description 35
- 238000004088 simulation Methods 0.000 claims abstract description 34
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 23
- 238000009877 rendering Methods 0.000 claims abstract description 14
- 210000001124 body fluid Anatomy 0.000 claims abstract description 8
- 239000010839 body fluid Substances 0.000 claims abstract description 8
- 238000001914 filtration Methods 0.000 claims abstract description 4
- 238000005286 illumination Methods 0.000 claims abstract description 3
- 238000004364 calculation method Methods 0.000 claims description 27
- 238000005315 distribution function Methods 0.000 claims description 20
- 238000012545 processing Methods 0.000 claims description 9
- 230000005540 biological transmission Effects 0.000 claims 1
- 238000003672 processing method Methods 0.000 claims 1
- 238000012546 transfer Methods 0.000 abstract description 2
- 230000005012 migration Effects 0.000 description 11
- 238000013508 migration Methods 0.000 description 11
- 230000033001 locomotion Effects 0.000 description 9
- 210000004027 cell Anatomy 0.000 description 8
- 230000000694 effects Effects 0.000 description 7
- 238000010586 diagram Methods 0.000 description 5
- 230000008569 process Effects 0.000 description 5
- 230000000007 visual effect Effects 0.000 description 5
- 238000013213 extrapolation Methods 0.000 description 4
- 230000008901 benefit Effects 0.000 description 3
- 230000007547 defect Effects 0.000 description 3
- 238000013461 design Methods 0.000 description 3
- 210000000712 G cell Anatomy 0.000 description 2
- 230000002457 bidirectional effect Effects 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 230000000704 physical effect Effects 0.000 description 2
- 239000007787 solid Substances 0.000 description 2
- 230000001133 acceleration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 239000007788 liquid Substances 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000005293 physical law Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/15—Correlation function computation including computation of convolution operations
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Computing Systems (AREA)
- Operations Research (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明涉及一种基于网格和粒子耦合的流体表面细节保护方法：(1)求解纳维‑斯特克斯方程(Navier‑Stokes，N‑S)采用网格方法，构成主体流体的仿真模型；(2)改进LBM‑VOF方法，并利用改进的VOF‑LBM耦合算法追踪流体的表面；(3)在异常表面网格位置生成粒子，然后通过粒子方法对粒子进行演化；(4)设计网格和粒子的耦合算法，将网格流体和粒子流体集成到同一个场景中，以此保证整个流场的物理守恒以及网格和粒子之间进行合理的物理信息传递；(5)最后利用屏幕空间方法进行真实感渲染，依次通过绘制球体，并计算每个像素点的深度值，进行深度滤波，根据深度值和像素点的位置信息求法向量，以及光照渲染，在GPU上实现一个逼真和实时的流体渲染。
Description
技术领域
本发明涉及一种基于网格和粒子耦合的流体表面细节保护方法，具体地说是一种基于物 理的图形学动画模拟方法，其部分技术涉及到构建高精度的主体流体仿真模型、构建与LBM 耦合的流体自由表面追踪模型、粒子的生成与演化、网格与粒子的耦合算法以及力学相关理 论。主要应用于各种真实物体三维逼真模拟，特别是应用于娱乐游戏、电影特效等领域。
背景技术
近些年来，计算机图形学领域研究的热点与难点之一就是真实感的场景、动画和现象的 模拟，它在制造工程、影视特效、三维动画、电脑游戏与虚拟现实等领域应用广泛。为了模 拟出和现实世界相差无几的真实感图形效果，可以基于物理原理和实验数据构建能准确描述 物理现象的物理模型。基于物理模型的自然现象的绘制可以模拟出更加逼真的效果，但计算 过于复杂，实时性得不到保证。因此，在基于物理的流体仿真中既要减少计算代价，又需要 尽可能的保持逼真的视觉效果。
当前计算流体力学中描述流体运动的基本方法有三种：基于粒子的拉格朗日方法、基于 网格的欧拉方法以及LBM方法。拉格朗日方法以SPH为代表，是一种基于粒子的方法，着眼于 流体质点的运动追踪上，考虑所有质点的运动，进而形成整个流体运动，其优点在于简单直 观，容易理解，可模拟流体的运动细节，但是计算量大，难以选取合适的光滑核函数以及在 特大变形的情况下容易发生畸变；欧拉方法基于网格，它不具体跟踪某一质点的具体运动过 程，而是将流场作为对象。欧拉法有利于分析物理场属性，容易构造液体运动表面，能够处 理变形与扭曲问题，但难以表现流体的表面细节；LBM方法是一种基于欧拉的“介观”的流 体仿真方法，它忽略了大量分子运动的细节，只保留了分子运动速度的部分统计特征，该方 法便于处理复杂边界，算法简单，具有很高的并行度，但依旧是基于网格的描述方法，同时 在细节仿真上依然存在不足。
文献1-Koerner C,Thies M,Hofmann T,et al.Lattice Boltzmann Model forFree Surface Flow for Modeling Foaming[J].Journal of Statistical Physics,2005, 121(1-2):179-196，针对在表面网格重构分布函数的问题，提出了一种基于动量交换法的表 面网格重构方法，但是此方法数值精度低，平衡状态下表面网格在视觉效果上还会发生抖动 现象，同时文献1还针对异常表面网格提出了一套特殊的处理方法，一定程度上缓解了该问 题，但是难以保证物理量的守恒。文献2-郝爱民,李帅,高阳.一种基于欧拉-拉格朗日耦 合方法的流体仿真方法，该方法虽然也将LBM与SPH耦合在一起，但是与本发明有着本质的不 同，首先是目标不同，文献2着重于增加流体细节，从而达到丰富的视觉效果，而本发明着 重于保护流体仿真的数值精度，从而达到保护流体细节的目的；然后是粒子的生成不同，文 献2为了达到增加细节的目的，在表面外生成粒子，而本发明主要是保护仿真的精度，用粒 子解决异常网格问题，因此只在异常网格处生成粒子；最后是耦合方式不同，文献2采用PLSM 进行表面追踪，并以此为媒介耦合LBM和SPH，本发明采用改进的VOF方法进行表面追踪，数 值精度更高。简而言之，文献2为了增加流体细节，可能会破坏原流体的演化规律，而本发 明是为了保护原流体演化规律的正常演化。
为了获得更好的真实感模拟，本发明提出了一种基于网格和粒子耦合的流体表面细节保 护方法，该方法主要由五部分组成：构建基于LBM的主体流体仿真模型、VOF追踪流体表面， 并标记异常界面网格、生成粒子并基于SPH方法进行演化、设计耦合算法耦合网格流体与粒 子流体以及真实感渲染，最后实现一个实时、多细节、逼真和稳定的流体仿真。
发明内容
本发明的技术解决问题：克服现有技术的不足，提供一种基于网格和粒子耦合的流体表 面细节保护方法，既可以摆脱传统网格法存在的诸多缺陷，又可以很好地保护流体的细节， 同时满足了流体仿真对实时性的要求。
本发明的技术解决方案：一种基于网格和粒子耦合的流体表面细节保护方法，步骤包括：
(1)求解纳维-斯特克斯方程(Navier-Stokes，N-S)采用格子莫尔兹曼方程(Lattice Boltzmann Method，LBM)，LBM的碰撞模型采用线性Bhatnagar-Gross-Krook(BGK)碰撞模 型，同时采用n维离散空间的m个速度多维离散网格模型(DnQm模型)，构成主体流体的仿 真模型，然后经过演化得到每个网格的物理信息；
(2)根据步骤(1)得到的每个网格的物理信息，利用改进的LBM-VOF耦合算法追踪流 体的表面，追踪流体表面，重构表面网格分布函数，该步骤得到了主体流体的表面以及表面 新的分布函数，并标记流体表面网格中异常表面网格的位置；
(3)在步骤(2)获得异常表面网格位置后，首先对异常表面网格进行处理，即在该位 置生成粒子代替原异常表面网格，同时粒子具备原异常表面网格的物理信息，然后通过Smooth Particle Hydrodynamics(SPH)方法对粒子进行演化，得到粒子新的物理信息；
(4)根据步骤(1)所得网格的物理信息和步骤(3)所得粒子的物理信息设计耦合算法，将网格流体和粒子流体集成到同一个场景中，并利用耦合算法将粒子流体转换为网格流 体，并重新计算网格的物理信息，包括质量、速度、位置，以此保证整个流场的物理守恒以 及网格和粒子之间进行合理的物理信息传递，并得到网格新的物理信息；
(5)将步骤(3)和步骤(4)所得到的粒子和网格的物理信息，利用屏幕空间方法，依次通过绘制球体，并计算每个像素点的深度值，深度滤波，根据深度值和像素点的位置求法向量，以及光照渲染，在GPU上实现一个逼真和实时的流体渲染。
本发明的原理：
(1)传统的LBM方法便于处理复杂边界、算法简单易实现以及并行性高等多个优势构 建对主体流体仿真的模型，LBM方法可以分为碰撞和迁移两部分，公式分别为：
fi(x+ciδt，t+δt)＝fi′(x，t)
其中，i为离散方向的下标，δt为单位时间步长，ci为离散速度的方向向量，ciδt为单位时间 的迁移步长，τ为无量纲松弛时间。
(2)传统的LBM-VOF耦合算法通过计算每个流体格子的质量变化来追踪流体的自由表 面，将模拟区域大致划分为三类网格：流体网格(记作”F”)，表面网格(记作”I”)和 没有流体的网格(记作”G”)，严格意义上I网格必须与F网格相邻，F网格必须被F网格 或I网格包围(即F网格不与G网格相邻)，同时F网格与G网格之间的转换必须先转换为I 网格。但在实际情况中，会出现I网格不与F网格相邻的情况(本发明将此类I网格称为异 常表面网格)，传统解决方案是特殊处理I网格的质量交换以及降低网格类型的转换条件， 但是会导致物理量不守恒和数值解的不稳定；针对LBM-VOF方法产生的异常表面网格，本发 明不再人工干预这些异常表面网格的演化，本发明采用粒子的方法处理这些异常表面网格， 利用SPH方法利于刻画流体细节的优势进行粒子的演化，保证了物理守恒和数值解的稳定性，同时丰富了流体表面的细节。
(3)在主体流体仿真过程中，在I网格处需要考虑G网格对I网格的作用力，而G网格的物理量是未知的，传统的LBM-VOF耦合算法通过动量交换法估计G网格对I网格的作用力，但是误差较大。本发明借鉴固体边界问题的压力非平衡态外推格式方法，将其耦合到LBM-VOF耦合算法中，此方法具有二阶精度，其具体公式为：
fi(I)＝fi(G)
其中，表面网格“I”在空网格“G”的i方向，i与
互为反向，fi(I)为表面网格重构后的分布函数在
方向上的分量，fi(G)为G网格对I网格作用力在i方向上的分量，
为G网格 i方向平衡态分量，
为G网格i方向非平衡态分量，ρG为G网格密度，一般为常数，u 为I网格速度，ρI为I网格密度，
为I网格
方向上的分量，其中
的计算不同于压力 非平衡态外推格式方法，其仿真结果视觉效果更逼真，流体细节更多。
(4)对异常表面网格进行处理，即在该异常表面网格位置生成粒子代替原异常表面网 格，同时粒子具备原异常表面网格的物理信息，然后通过SPH方法对粒子进行演化，得到粒 子新的物理信息，具体步骤如下：
(41)计算网格转换为粒子后，粒子的位置信息，具体公式如下：
posparticle＝posgrid-(1-ε)·n
其中，posgrid为表面网格的位置，posparticle为粒子的位置，ε为I网格的体积积分，n为 表面网格的单位法向量。
(42)物理量的传递，具体公式如下：
massparticle＝massgrid
vparticle＝vgrid
ρparticle＝ε·mgrid
其中，massgrid和vgrid分别为表面网格的质量与速度，massparticle和vparticle为所生成 粒子的质量与速度，ε为表面网格的体积积分，ρparticle为粒子的密度。
(43)粒子位置信息和物理信息计算完成后，删除转化为粒子的网格，并通过SPH方法 对粒子进行演化，得到粒子新的物理信息；
(5)网格流体与粒子流体耦合算法的设计，主要是为了保证模拟区域内流体在物理规 律的约束下演化，具体步骤如下：
(51)计算粒子所在的网格位置，粒子所在的网格位置具体计算公式如下：
其中，posparticle为粒子演化后的位置信息，posgrid为粒子所在网格的位置。
(52)如果步骤(51)所求的posgrid所在位置的类型为表面网格，继续执行步骤(53)。
(53)为保证流体仿真的质量守恒，需要对表面网格的质量进行重新计算，计算公式如 下：
m′grid＝mparticle+mgrid
其中，m′grid为表面网格重新计算后的质量，mparticle为粒子的质量，mgrid为表面网格 重新计算前的质量。
(54)为保证流体仿真的动量守恒，根据步骤(53)计算的m′grid，重新计算表面网格的速度，新的速度v′grid要满足以下公式：
m′gridv′grid＝mpraticle vparticle+mgridvgrid
其中，v′grid为表面网格重新计算后速度，m′grid为步骤(53)重新计算后的表面网格质 量，mpraticle，vparticle分别为粒子的质量和速度，mgrid，vgrid分别为表面网格重新计算前的质量和速度。
(55)根据步骤(54)所计算的v′grid，重新计算表面网格的分布函数，计算公式如下：
其中，i为离散方向的下标，τ为无量纲松弛时间，fi′为表面网格计算后的在i方向上分布 函数分量，fi为表面网格计算前在i方向上分布函数分量，v′grid为重新计算后速度，fi (eq)(v′grid) 为根据v′grid计算的在i方向上的平衡分布函数。
(56)删除所在位置在表面网格的粒子，并由步骤(51)、(52)、(53)、(54)、(55)更新了表面网格的物理信息，得到新的表面网格；
本发明与现有技术相比的优点在于：
(1)本发明提出一种新的重建表面网格分布函数的方法，从而改进了LBM-VOF方法， 使LBM-VOF算法具有更高的鲁棒性以及提高了流体表面追踪的精度。
(2)本发明针对异常表面网格问题，设计了一种新的网格-粒子双向耦合方法，在保证 流体仿真的物理守恒的前提下，解决了异常表面网格的问题，保护了流体表面的细节，还通 过粒子方法一定程度的丰富了表面细节，最终有效提高了视觉效果的逼真度。
(3)本发明提出的网格-粒子双向耦合方法，其中LBM的空间网格可以作为SPH领域搜 索的空间网格，保证了LBM方法和SPH方法的契合度，减少了不必要的计算；同时继承了LBM 具有很高并行性的特点，可以通过GPU并行技术加速算法。
附图说明
图1为本发明方法实现流程图；
图2为I网格迁移示意图；
图3为正常网格分布示意图；
图4为异常表面网格处理示意图；
图5为粒子与网格耦合示意图；
图6为I网格绘制位置计算示意图；
图7为本发明方法效果对比图；
图8为最终真实感渲染图。
具体实施方式
下面结合附图实施例对本发明进行具体描述。
本发明一种基于网格和粒子耦合的流体表面细节保护方法，如图1所示，具体实现过程 如下：
(1)构建基于LBM的主体流体仿真模型
LBM方法是一种基于欧拉的“介观”的流体仿真方法，它忽略了大量分子运动的细节， 只保留了分子运动速度的部分统计特征，该方法便于处理复杂边界，算法简单，具有很高的 并行度。Boltzmann方程可以表示为：
其中，f为分子的分布函数，ξ为分子速度，a为分子加速度，Q(f，f)为碰撞项。
在boltzmann方程的基础上用分子分布函数代替分子本身进行演化，其演化方程直接采 用格子boltzmann方程，并根据分布函数直接计算流体的密度和速度，其方程被分为碰撞和 迁移两个步骤，具体方程分别为：
fi′(x，t)＝fi(x，t)+Q(f，f)
fi(x+ciδt，t+δt)＝fi′(x，t)
其中，i为离散方向的下标，δt为单位时间步长，ci为离散速度的方向向量，ciδt为单位 时间的迁移步长。
同时通过BGK碰撞模型将碰撞项线性化，新的碰撞步计算方程如下：
其中，τ为无量纲松弛时间。
根据Boltzmann H定理和DnQm网格离散模型，fi (eq)(x，t)的计算公式可以近似为：
其中，
ωi为i方向分布函数的权重，ei为i方向的方向向量，ρ为宏观密度，u为 宏观速度。
(2)构建与LBM耦合的流体自由表面追踪
表面追踪是指提取流体的表面信息，表面的精度对仿真的精度有着重要的影响。本发明 的流体自由表面追踪可以分为以下三个步骤：
(21)在迁移过程中，表面网格需要估计G网格对I网格的作用力fi(G)，如图2所示，由于G网格未参与演化，其分布函数是未知的，因此需要估计G网格对I网格的作用力，本 发明借鉴固体边界问题的压力非平衡态外推格式方法，将其耦合到LBM-VOF耦合算法中，此方法具有二阶精度，其具体公式为：
fi(I)＝fi(G)
其中，表面网格“I”在空网格“G”的i方向，i与
互为反向，fi(I)为表面网格重构后的分布函数在
方向上的分量，fi(G)为G网格对I网格作用力在i方向上的分量，为G网格 i方向平衡态分量，为G网格i方向非平衡态分量，ρG为G网格密度，一般为常数，u 为I网格速度，ρI为I网格密度，
为I网格
方向上的分量，其中
的计算不同于压力 非平衡态外推格式方法，其仿真结果视觉效果更逼真，流体细节更多。
(22)然后在LBM方法迁移步计算I网格与I网格、F网格与F网格和I网格与F网格之间的 质量迁移，其中F网格与F网格和I网格与F网格之间的i方向质量迁移计算公式为：Δmi＝fin- fout，其中fin为i方向进来的量，fout为从i方向出去的量，I网格与I网格的质量迁移计算公式 为：Δmi＝ε·(fin-fout)，
其中ε2，ε1为两网格的体积积分即质量与密度之比。
(23)根据迁移后的质量，重新计算体积积分，再根据体积积分更新表面网格的类型， 具体判断如下：ε＞1，I网格被充满，变为F网格；ε＜0，I网格被置空，变为G网格；同时要从新分配那些被充满或被置空网格的质量，以及保证F网格不与G网格相邻，正常流体的网 格分布如图3所示，F网格被F网格或I网格包围，并且不与G网格相邻，I网格与F网格相邻。 然而实际中，会出现不与F网格相邻的I网格，本发明称为异常表面网格，如图4所示，I网格 就为异常表面网格。本发明提出了一种全新的方法处理异常表面网格，从粒子的角度出发， 用粒子代替异常表面网格进行演化，因此本步骤要标记这些异常表面网格，为下一步做准备。
(3)粒子的生成与演化
(31)粒子的生成：如图4所示，对异常表面网格进行处理，即在该异常表面网格位置 生成粒子代替原异常表面网格，同时粒子具备原异常表面网格的物理信息，具体步骤如下：
1)计算网格转换为粒子后，粒子的位置信息，具体公式如下：
posparticle＝posgrid-(1-ε)·n
其中，posgrid为表面网格的位置，posparticle为粒子的位置，ε为I网格的体积积分，n为 表面网格的法向量。
2)物理量的传递，具体公式如下：
massparticle＝massgrid
vparticle＝vgrid
ρparticle＝ε·mgrid
其中，massgrid和vgrid分别为表面网格的质量与速度，massparticle和vparticle为所生成 粒子的质量与速度，ε为表面网格的体积积分，ρparticle为粒子的密度。
3)粒子位置信息和物理信息计算完成后，在posparticle位置生成新的粒子，并且粒子质 量、速度和密度分别为massparticle、vparticle和ρparticle。
(32)粒子的演化：本发明采用光滑粒子流体动力学方法(SPH)，SPH方法是一种利用 粒子插值法计算流体力学变量的拉格朗日无网格方法。在SPH方法中，整个流场被离散成一 系列粒子，每个粒子都携带了一些物理信息，如质量、速度、位置等。SPH方程的构造通常 按照两个关键步进行。第一步为积分表示法，又称“核函数逼近”；第二步为粒子近似法，称“粒子逼近”。从而将积分形式转化成粒子求和的级数形式，从而场内的每个粒子的物理属性可以通过支持域内所有粒子共同插值得到。因此，得到如下表达形式：
其中，范围j为xi粒子支持域内的粒子，mj为粒子j的质量，ρi是指的所求粒子j的密度。 f(x)代表了粒子的任意一种物理属性，W是光滑核函数，h是光滑核函数的影响半径。粒子运 动由经典的Navier-Stokes方程(N-S方程)控制，SPH需要对不可压缩流体的N-S方程进行离 散化：
其中，若μ是动力粘性系数，v是μ/ρ为运动粘性系数，ρ是密度，p是压强，
是速度，F 是重力等外力。方程等号右边依次为粘滞力、压力和外力，每个时间步进行各种力的求解， 获得粒子的速度，更新位置，驱动粒子运动。
(4)耦合网格流体与粒子流体
网格流体与粒子流体耦合算法的设计，主要是为了保证模拟区域内流体在物理规律的约 束下演化，如图5所示，粒子演化后所在网格类型是I网格，符合耦合条件，耦合算法具体 步骤如下：
(41)计算粒子所在的网格位置，粒子所在的网格位置具体计算公式如下：
其中，posparticle为粒子演化后的位置信息，posgrid为粒子所在网格的位置。
(42)如果步骤(41)所求的posgrid所在位置的类型为表面网格，继续执行步骤(43)， 否则执行步骤(47)。
(43)为保证系统的质量守恒，需要对表面网格的质量进行重新计算，计算公式如下：
m′grid＝mparticle+mgrid
其中，m′grid为表面网格重新计算后的质量，mparticle为粒子的质量，mgrid为表面网格 重新计算前的质量。
(44)为保证系统的动量守恒，根据步骤(43)计算的m′grid，重新计算表面网格的速度，新的速度v′grid要满足以下公式：
m′gridv′grid＝mpraticlevparticle+mgridvgrid
其中，v′grid为表面网格重新计算后速度，m′grid为步骤(43)重新计算后的表面网格质 量，mpraticle，vparticle分别为粒子的质量和速度，mgrid，vgrid分别为表面网格重新计算前的质量和速度。
(45)根据步骤(44)所计算的v′grid，重新计算表面网格的分布函数，计算公式如下：
其中，i为离散方向的下标，τ为无量纲松弛时间，fi′为表面网格计算后的在i方向上分布 函数分量，fi为表面网格计算前在i方向上分布函数分量，v′grid为重新计算后速度，fi (eq)(v′grid) 为根据v′grid计算的在i方向上的平衡分布函数。
(46)删除所在位置在表面网格的粒子。
(47)剩下的网格和粒子继续演化，然后重复步骤(41)、(42)、(43)、(44)、(45)、(46)、 (47)。
(5)真实感渲染
本发明采用屏幕空间法渲染流体。将步骤(4)所得到的流体信息通过绘制精灵点，深 度滤波，根据深度值和UV值求法向量等一系列步骤在GPU上实现一个逼真和实时的流体渲 染。其中精灵点的半径与体积积分相关联，即：
其中k与空间的维度有关系，二 维空间为
三维空间为
如图6所示，为了避免仿真结果在视觉上流体的不连续感，绘制I网格的精灵点时，需要重新计算位置，使流体在空间上和视觉上更连续，具体计算公式如下：
pos′new＝pos+(1-ε)·ei
其中，pos为I网格的原位置，pos′new为，ε为I网格的体积积分，ei表示F网格在I网格的i方向的方向向量。
图7展示了传统LBM-VOF方法和本发明的方法在经过相同时间演化后的效果对比图，左 图为传统LBM-VOF方法，可以看出相较与本发明丢失了很多细节，右图为本发明的方法，其 中画圈内(不仅限于圈内)即为使用本发明所产生的细节，可以看出本发明通过网格-粒子 耦合方法不仅保证了网格的仿真精度，保护了流体表面的细节，还通过粒子方法增加了表面 的细节；图8为利用本发明方法进行流体仿真，再利用屏幕空间法真实感绘制的效果图，由 于算法核心计算均在GPU上完成，因此帧率也达到105FPS(网格大小：40X40X40，GPU：NVIDIA Quadro K5200)，完全满足流体仿真对实时性的要求。
以上所述，仅为本发明部分具体实施方式，但本发明的保护范围并不局限于此，任何熟 悉本领域的人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明 的保护范围之内。
Claims (4)
1.一种基于网格和粒子耦合的流体表面细节保护方法，其特征在于，步骤包括：
(1)采用格子莫尔兹曼方程LBM(Lattice Boltzmann Method，LBM)求解纳维-斯特克斯方程(Navier-Stokes，N-S)，LBM的碰撞模型采用线性Bhatnagar-Gross-Krook(BGK)碰撞模型，同时采用n维离散空间的m个速度多维离散网格模型(DnQm模型)，构成主体流体的仿真模型，然后经过演化得到每个网格的物理信息；
(2)根据步骤(1)得到的每个网格的物理信息，利用改进的LBM-VOF耦合算法追踪流体表面，并重构流体表面网格的分布函数，得到主体流体的表面以及表面新的分布函数，并标记流体表面网格中异常表面网格位置；
(3)在步骤(2)获得异常表面网格位置后，首先对异常表面网格进行处理，即在该异常表面网格位置生成粒子代替原异常表面网格，同时粒子具备原异常表面网格的物理信息，然后通过Smooth Particle Hydrodynamics(SPH)方法对粒子进行演化，得到粒子新的物理信息；
(4)根据步骤(1)所得网格的物理信息和步骤(3)所得粒子的物理信息设计耦合算法，将网格流体和粒子流体集成到同一个场景中，并利用耦合算法将粒子流体转换为网格流体，并重新计算网格的物理信息，物理信息包括质量、速度、位置，保证整个流场的物理守恒以及网格和粒子之间进行合理的物理信息传递，并得到网格新的物理信息；
(5)将步骤(3)得到的粒子新的物理信息和步骤(4)网格新的物理信息，利用屏幕空间方法，依次通过绘制球体，并计算每个像素点的深度值，进行深度滤波，根据深度值和像素点的位置信息求法向量，以及光照渲染，在GPU上实现一个逼真和实时的流体渲染。
2.根据权利要求1所述的一种基于网格和粒子耦合的流体表面细节保护方法，其特征在于：所述步骤(2)中，重构流体表面网格的分布函数fi(I)，具体公式如下：
fi(I)＝fi(G)
其中，表面网格“I”在空网格“G”的i方向，i与
互为反向,fi(I)为表面网格重构后的分布函数在
方向上的分量，fi(G)为G网格对I网格作用力在i方向上的分量，
为G网格i方向平衡态分量，
为G网格i方向非平衡态分量，ρG为G网格密度，为常数，u为I网格速度，ρI为I网格密度，
为I网格
方向上的分量。
3.根据权利要求1所述的一种基于网格和粒子耦合的流体表面细节保护方法，其特征在于：所述步骤(3)中异常表面网格的处理方法是用粒子代替异常表面网格，网格转换为粒子的具体公式如下：
(31)计算网格转换为粒子后，粒子的位置信息，具体公式如下：
posparticle＝posgrid-(1-ε)·n
其中，posgrid为表面网格的位置，posparticle为粒子的位置，ε为I网格的体积积分，n为表面网格的法向量；
(32)物理量的传递，具体公式如下：
massparticle＝massgrid
vparticle＝vgrid
ρparticle＝ε·mgrid
其中，massgrid和vgrid分别为表面网格的质量与速度，massparticle和vparticle为所生成粒子的质量与速度，ε为表面网格的体积积分，ρparticle为粒子的密度；
(33)粒子位置信息和物理信息计算完成后，删除转化为粒子的网格，并通过SPH方法对粒子进行演化，得到粒子新的物理信息。
4.根据权利要求1所述的一种基于网格和粒子耦合的流体表面细节保护方法，其特征在于：所述步骤(4)中将粒子流体转换为网格流体的具体步骤如下：
(41)计算粒子所在的网格位置，粒子所在的网格位置具体计算公式如下：
其中，posparticle为粒子演化后的位置信息，posgrid为粒子所在网格的位置；
(42)如果步骤(41)所求的posgrid所在位置的类型为表面网格，继续执行步骤(43)，否则执行步骤(47)；
(43)为保证流体仿真的质量守恒，对表面网格的质量进行重新计算，计算公式如下：
m′grid＝mparticle+mgrid
其中，m′grid为表面网格重新计算后的质量，mparticle为粒子的质量，mgrid为表面网格重新计算前的质量；
(44)为保证流体仿真的动量守恒，根据步骤(43)计算的m′grid，重新计算表面网格的速度，新的速度v′grid满足以下公式：
m′gridv′grid＝mpraticlevparticle+mgridvgrid
其中，v′grid为表面网格重新计算后速度，m′grid为步骤(43)重新计算后的表面网格质量，mpraticle，vparticle分别为粒子的质量和速度，mgrid，vgrid分别为表面网格重新计算前的质量和速度；
(45)根据步骤(44)所计算的v′grid，重新计算表面网格的分布函数，计算公式如下：
其中，i为离散方向的下标，τ为无量纲松弛时间，fi′为表面网格计算后的在i方向上分布函数分量，fi为表面网格计算前在i方向上分布函数分量，v′grid为重新计算后速度，fi (eq)(v′grid)为根据v′grid计算的在i方向上的平衡分布函数；
(46)删除所在位置在表面网格的粒子，并由步骤(41)、(42)、(43)、(44)、(45)更新了表面网格的物理信息，得到新的表面网格。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910956536.3A CN110717269B (zh) | 2019-10-10 | 2019-10-10 | 一种基于网格和粒子耦合的流体表面细节保护方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910956536.3A CN110717269B (zh) | 2019-10-10 | 2019-10-10 | 一种基于网格和粒子耦合的流体表面细节保护方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110717269A true CN110717269A (zh) | 2020-01-21 |
| CN110717269B CN110717269B (zh) | 2023-07-25 |
Family
ID=69212333
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910956536.3A Active CN110717269B (zh) | 2019-10-10 | 2019-10-10 | 一种基于网格和粒子耦合的流体表面细节保护方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110717269B (zh) |
Cited By (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111368487A (zh) * | 2020-03-17 | 2020-07-03 | 广西师范大学 | 一种基于晶格Boltzmann模型模拟颗粒周期性运动的流场处理方法 |
| CN112347708A (zh) * | 2020-09-28 | 2021-02-09 | 西安电子科技大学 | 一种基于物理过程的实时光学烟幕的渲染方法 |
| CN112862942A (zh) * | 2021-02-08 | 2021-05-28 | 腾讯科技（深圳）有限公司 | 物理特效模拟方法、装置、电子设备和存储介质 |
| CN113935254A (zh) * | 2020-06-29 | 2022-01-14 | 达索系统西姆利亚公司 | 使用表面算法模拟物理过程的计算机系统 |
| CN115461159A (zh) * | 2020-05-08 | 2022-12-09 | 发那科株式会社 | 模拟装置 |
| CN115618663A (zh) * | 2021-07-16 | 2023-01-17 | 合肥本源量子计算科技有限责任公司 | 一种网格方程与物理方程耦合的量子求解方法及装置 |
| CN117475040A (zh) * | 2023-10-09 | 2024-01-30 | 北京航空航天大学 | 基于隐式粒子插值的格子玻尔兹曼流体动画仿真方法 |
| CN117593486A (zh) * | 2024-01-19 | 2024-02-23 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于空间粒子的网格重构方法和装置 |
| WO2024099206A1 (zh) * | 2022-11-11 | 2024-05-16 | 华为技术有限公司 | 一种图形界面处理方法以及装置 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20090267951A1 (en) * | 2008-04-28 | 2009-10-29 | Institute For Information Industry | Method for rendering fluid |
| CN103679802A (zh) * | 2013-12-01 | 2014-03-26 | 北京航空航天大学 | 基于屏幕空间的sph流体表面实时绘制方法 |
| CN104268943A (zh) * | 2014-09-28 | 2015-01-07 | 北京航空航天大学 | 一种基于欧拉-拉格朗日耦合方法的流体仿真方法 |
| CN108717722A (zh) * | 2018-04-10 | 2018-10-30 | 天津大学 | 基于深度学习和sph框架的流体动画生成方法及装置 |
-
2019
- 2019-10-10 CN CN201910956536.3A patent/CN110717269B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20090267951A1 (en) * | 2008-04-28 | 2009-10-29 | Institute For Information Industry | Method for rendering fluid |
| CN103679802A (zh) * | 2013-12-01 | 2014-03-26 | 北京航空航天大学 | 基于屏幕空间的sph流体表面实时绘制方法 |
| CN104268943A (zh) * | 2014-09-28 | 2015-01-07 | 北京航空航天大学 | 一种基于欧拉-拉格朗日耦合方法的流体仿真方法 |
| CN108717722A (zh) * | 2018-04-10 | 2018-10-30 | 天津大学 | 基于深度学习和sph框架的流体动画生成方法及装置 |
Non-Patent Citations (1)
| Title |
|---|
| 张凤全等: "细节保护的流体表面绘制方法" * |
Cited By (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111368487A (zh) * | 2020-03-17 | 2020-07-03 | 广西师范大学 | 一种基于晶格Boltzmann模型模拟颗粒周期性运动的流场处理方法 |
| CN111368487B (zh) * | 2020-03-17 | 2023-07-18 | 广西师范大学 | 一种基于晶格Boltzmann模型模拟颗粒周期性运动的流场处理方法 |
| CN115461159A (zh) * | 2020-05-08 | 2022-12-09 | 发那科株式会社 | 模拟装置 |
| CN113935254A (zh) * | 2020-06-29 | 2022-01-14 | 达索系统西姆利亚公司 | 使用表面算法模拟物理过程的计算机系统 |
| CN112347708A (zh) * | 2020-09-28 | 2021-02-09 | 西安电子科技大学 | 一种基于物理过程的实时光学烟幕的渲染方法 |
| CN112862942A (zh) * | 2021-02-08 | 2021-05-28 | 腾讯科技（深圳）有限公司 | 物理特效模拟方法、装置、电子设备和存储介质 |
| CN112862942B (zh) * | 2021-02-08 | 2024-04-09 | 腾讯科技（深圳）有限公司 | 物理特效模拟方法、装置、电子设备和存储介质 |
| CN115618663A (zh) * | 2021-07-16 | 2023-01-17 | 合肥本源量子计算科技有限责任公司 | 一种网格方程与物理方程耦合的量子求解方法及装置 |
| WO2024099206A1 (zh) * | 2022-11-11 | 2024-05-16 | 华为技术有限公司 | 一种图形界面处理方法以及装置 |
| CN117475040A (zh) * | 2023-10-09 | 2024-01-30 | 北京航空航天大学 | 基于隐式粒子插值的格子玻尔兹曼流体动画仿真方法 |
| CN117593486A (zh) * | 2024-01-19 | 2024-02-23 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于空间粒子的网格重构方法和装置 |
| CN117593486B (zh) * | 2024-01-19 | 2024-04-09 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于空间粒子的网格重构方法和装置 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110717269B (zh) | 2023-07-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110717269A (zh) | 一种基于网格和粒子耦合的流体表面细节保护方法 | |
| Stam | Real-time fluid dynamics for games | |
| Feldman et al. | Animating gases with hybrid meshes | |
| He et al. | Local poisson SPH for viscous incompressible fluids | |
| US7565276B2 (en) | Method of simulating detailed movements of fluids using derivative particles | |
| Song et al. | Stable but nondissipative water | |
| Shao et al. | Stable and fast fluid–solid coupling for incompressible SPH | |
| Bargteil et al. | Animation of deformable bodies with quadratic Bézier finite elements | |
| Wang et al. | Solving general shallow wave equations on surfaces | |
| CN108717722A (zh) | 基于深度学习和sph框架的流体动画生成方法及装置 | |
| Thürey et al. | Interactive free surface fluids with the lattice Boltzmann method | |
| Liu et al. | Turbulent details simulation for SPH fluids via vorticity refinement | |
| Shao et al. | Realistic and stable simulation of turbulent details behind objects in smoothed‐particle hydrodynamics fluids | |
| CN109215100A (zh) | 一种混合流体相变动画生成方法及装置 | |
| Fan et al. | Adapted unstructured LBM for flow simulation on curved surfaces | |
| Zhang et al. | A SPH-based method for interactive fluids simulation on the multi-GPU | |
| Gao et al. | Accelerating liquid simulation with an improved data‐driven method | |
| Im et al. | Visual simulation of rapidly freezing water based on crystallization | |
| Wendt et al. | Finite volume flow simulations on arbitrary domains | |
| Cline et al. | Fluid flow for the rest of us: Tutorial of the marker and cell method in computer graphics | |
| Wu et al. | Improved divergence‐free smoothed particle hydrodynamics via priority of divergence‐free solver and SOR | |
| Fujisawa et al. | Particle-based shallow water simulation with splashes and breaking waves | |
| Yang et al. | Interactive coupling between a tree and raindrops | |
| Alduán et al. | Efficient and Robust Position-Based Fluids for VFX. | |
| Abdelnaim et al. | Fluid-structure interactions simulation and visualization using ISPH approach |
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 |