CN102831275A - 一种3d流体的仿真方法及系统 - Google Patents

Abstract

本发明提供了一种3D流体的仿真方法及系统，相应的方法包括：搜索每个流体粒子的邻居粒子，并通过核函数解出Navier-Stokes状态方程，获得流体粒子的粒子位置和粒子密度；计算每个流体粒子的象素浓度，再将流体粒子的象素浓度进行平滑处理后输出；将流体粒子的密度和位置转化成一个空间密度场，形成表面数据并得出表面的法向量；将表面数据和表面法向量处理成面中心和法向量格式，在光线的位置通过光子效应绘制屏幕坐标；根据象素流体粒子的象素浓度、屏幕坐标通过变换矩阵转换到世界坐标系下，并进行贴图和颜色合成，最后将合成的颜色输出。本发明无需依赖网格且计算速度较快、计算量和存储空间较小，能够获得较好的渲染效果。

Description

一种3D流体的仿真方法及系统

技术领域

本发明涉及一种3D流体的仿真方法及系统，属于图形仿真技术领域。

背景技术

SPH(Smoothed Particle Hydrodynamics，光滑粒子流体动力学)是现有技术中一种常用的流体仿真技术，其基本思想是将连续的流体(或固体)用相互作用的质点组来描述，各个物质点上承载各种物理量，包括质量、速度等，通过求解质点组的动力学方程和跟踪每个质点的运动轨道，求得整个系统的力学行为。

现有的SPH仿真技术主要包括Sphysics、Volume Rendering of SPH Data、OpenFOAM、Fluids v2.0和OpenTissue OpenTissue，其中的Sphysics是一个开源的SPH仿真平台，但由于其用Fortran语言编写且偏重于物理过程，导致计算速度相对较慢，且渲染效果较差；Volume Rendering of SPH Data是Stefan Auer模拟SPH流体使用的原型实例，但是基于体绘制的渲染，运算量和存储量巨大；OpenFOAM是一种基于网格的开源流体仿真系统，这与SPH的无网格理论相悖，网格方法计算准确但速度较慢，存储量大；Fluids v2.0是一个快速SPH计算系统，但该系统只给出了SPH的计算过程，缺乏有效的表面重绘算法和其它光影特效，绘制效果单一；OpenTissue OpenTissue仿真流体的几何边界较为固定，各个模块灵活性较差，不提供表面重建的模块，也没有光影特效。

发明内容

本发明为解决现有的SPH仿真技术中存在的必须依赖网格、计算速度较慢、渲染效果较差、运算量和存储量较大的问题，进而提供了一种3D流体的仿真方法及系统。为此，本发明提供了如下的技术方案：

一种3D流体的仿真方法，包括：

通过空间Hash方法搜索每个流体粒子的邻居粒子，并通过核函数解出Navier-Stokes状态方程，获得流体粒子的粒子位置和粒子密度，并输出；

通过OpenGL颜色混合模式和Billboard方法并根据每个流体粒子的粒子位置，计算出每个流体粒子的象素密度，再通过Cuda将所述流体粒子的象素浓度进行平滑处理后输出；

将所述流体粒子的密度和位置转化成一个空间密度场，通过设定一个阈值将空间密度场转换成Marching Cube数据构造，形成表面数据并得出表面的法向量；

将所述表面数据和表面法向量处理成面中心和法向量格式，根据Snell原理，计算出光线的位置，在所述光线的位置通过光子效应绘制屏幕坐标；

根据所述象素流体粒子的象素浓度、屏幕坐标通过变换矩阵转换到世界坐标系下，并根据转换到世界坐标系下的所述象素流体粒子的象素浓度、屏幕坐标进行贴图和颜色合成，最后将合成的颜色输出。

一种3D流体的仿真系统，包括：

粒子状态计算模块，用于通过空间Hash方法搜索每个流体粒子的邻居粒子，并通过核函数解出Navier-Stokes状态方程，获得流体粒子的粒子位置和粒子密度，并输出；

粒子浓度计算模块，用于通过OpenGL颜色混合模式和Billboard方法并根据每个流体粒子的粒子位置，计算出每个流体粒子的象素浓度，再通过Cuda将所述流体粒子的象素浓度进行平滑处理后输出；

Marching Cube表面重建模块，用于将所述流体粒子的密度和位置转化成一个空间密度场，通过设定一个阈值将空间密度场转换成Marching Cube数据构造，形成表面数据并得出表面的法向量；

Caustics效应计算模块，用于将所述表面数据和表面法向量处理成面中心和法向量格式，根据Snell原理，计算出光线的位置，在所述光线的位置通过光子效应绘制屏幕坐标；

CGShader着色模块，用于根据所述象素流体粒子的象素浓度、屏幕坐标通过变换矩阵转换到世界坐标系下，并根据转换到世界坐标系下的所述象素流体粒子的象素浓度、屏幕坐标进行贴图和颜色合成，最后将合成的颜色输出。

本发明实施方式提供的技术方案通过Navier-Stokes状态方程获得流体粒子的粒子位置和粒子密度，无需依赖网格且计算速度较快、计算量和存储空间较小；借助于Cuda的GPU平台，进一步提升了运算速度；使用OpenGL混合和Billboard技术，使OpenGL可以自动生成浓度信息，避免了不必要的计算，同时加快了计算速度；在Marching Cube算法的基础上，利用光线的折射原理，计算面片的折射位置，通过OpenGL进行混合，能够获得较好的渲染效果。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1是本发明的具体实施方式提供的3D流体的仿真方法的流程示意图；

图2是本发明的具体实施方式提供的核函数一维曲线图；

图3是本发明的具体实施方式提供的采用蛙跳式更新的粒子的运动状态示意图；

图4是本发明的具体实施方式提供的空间Hash方法查找邻居粒子的示意图；

图5是本发明的具体实施方式提供的采用固体碰撞模拟固壁反弹作用的示意图；

图6是本发明的具体实施方式提供的Caustics效应计算示意图。

图7是本发明的具体实施方式提供的3D流体的仿真系统的结构示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本具体实施方式提供的是一种基于SPH的拉格朗日方法，该方法是一种流体仿真引擎，通过跟踪粒子状态，从而模拟整个流体状态。该流体仿真引擎的设计是基于OpenGL图形接口之上，借助于NVIDIA GPU的Cuda的强大并行计算平台，从而实现3D流体的仿真过程和渲染过程。为了解决相应的技术问题，本发明的具体实施方式提供了一种3D流体的仿真方法，如图1所示，该方法具体可以包括：

步骤11，通过空间Hash方法搜索每个流体粒子的邻居粒子，并通过核函数解出Navier-Stokes状态方程，获得流体粒子的粒子位置和粒子密度，并输出。

具体的，Navier-Stokes方程的代数形式可以表示成：

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

$\▿\·u=0$

其中ρ表示密度，u表示速度，p表示压力，μ表示粘度系数，f表示外力。

SPH方法是基于插值的方法，插值的函数称之为核函数，其代数形式如下：

A_I(r)＝∫_ΩA(r`)W(r-r`,h)dr`

离散形式如下：

其中，Ω是差值区域，r`是差值区域内的积分变量。

$A_S\left(r\right)=\underset{j}{\mathrm{\Σ}}{A}_j{V}_{j}W(r-r_j,h)$

按照体积与密度的关系：

$V=\frac{m}{\mathrm{\ρ}}$

所以离散化形式可以变换成：

$A_S\left(r\right)=\underset{j}{\mathrm{\Σ}}{A}_j\frac{m_j}{{\mathrm{\ρ}}_j}W(r-r_j,h)$

对等式两边求偏导：

$\frac{\∂}{\∂x}{A}_S\left(r_i\right)=\frac{\∂}{\∂x}\left(A_j\frac{m_j}{{\mathrm{\ρ}}_j}W(r-r_j,h)\right)$

化简得：

$\frac{\∂}{\∂x}{A}_S\left(r_i\right)=\frac{\∂}{\∂x}\left(A_j\frac{m_j}{{\mathrm{\ρ}}_j}\right)W(r-r_j,h)+\left(A_j\frac{m_j}{{\mathrm{\ρ}}_j}\right)\frac{\∂}{\∂x}W(r-r_j,h)$

$=0\·W(r-r_j,h)+A_j\frac{m_j}{{\mathrm{\ρ}}_j}\frac{\∂}{\∂x}W(r-r_j,h)$

$=A_j\frac{m_j}{{\mathrm{\ρ}}_j}\frac{\∂}{\∂x}W(r-r_j,h)$

其中，V是积分元体积，m是积分元质量，d(r`)＝v_j＝m_j/ρ_j。

上面等式说明对As的导数形式和对核函数的导数形式相等。

核函数一般要满足4个条件，首先是归一化特性，其数学形式如下：

∫_ΩW(r,h)dr＝1

非负性，即核函数在每一点都应该非负：

W(r,h)≥0

对称性，即在对称位置上的值相等：

W(r,h)＝W(-r,h)

还应该满足极值性：

$\underset{h\→0}{\lim }W(r,h)=\mathrm{\δ}\left(r\right)$

其中

本具体实施方式采用的多项式核函数，其形式如下：

$W_{\mathrm{default}}(r,t)=\frac{315}{64\mathrm{\π}{h}^9}\left\{\begin{array}{c}{(r^2-{\left|\right|r\left|\right|}^2)}^3,0\≤\left|\right|r\left|\right|\≤h\\ 0,\left|\right|r\left|\right|>h\end{array}\right.$

其一维曲线图如图2所示。

流体粒子的象素密度计算公式为：

${\mathrm{\ρ}}_i=\mathrm{\ρ}\left(r_i\right)=\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j\frac{m_j}{{\mathrm{\ρ}}_j}W(r_i-r_j,h)=\underset{j}{\mathrm{\Σ}}{m}_{j}W(r_i-r_j,h)$

流体粒子的象素压力计算公式为：

$f_i^{\mathrm{pressure}}=-\underset{j\≠i}{\mathrm{\Σ}}\frac{p_i+p_j}{2}\frac{m_j}{{\mathrm{\ρ}}_j}\▿W(r_i-r_j,h)$

其中压力的计算使用拟压缩方法，即将流体看成可压缩的状态，其计算方法如下：

流体粒子的象素粘力计算公式为：

$\mathrm{\θ}={\left(\frac{\mathrm{\ρ}}{{\mathrm{\ρ}}_0}\right)}^7$

$\mathrm{\&epsiv;}=\left\{\begin{array}{c}\mathrm{\&theta;},a<\mathrm{\&theta;}<b\\ a,\mathrm{\&theta;}<a\\ b,\mathrm{\&theta;}>b\end{array}\right.$

p＝k`ρ₀(ε-1.0)

$f_i^{\mathrm{vis}\cos \mathrm{ity}}=\mathrm{\&mu;}\underset{j\&NotEqual;i}{\mathrm{\&Sigma;}}(u_j-u_i)\frac{m_j}{{\mathrm{\&rho;}}_j}{\&dtri;}^{2}W(r_i-r_j,h)$

流体粒子的象素粒子力组成包括粒子受到的力包括内力和外力。

F＝f^int ernal+f^external

内力包括上述的压力、粘力等，而最典型的外力就是重力。

$f_i^{\mathrm{gravity}}={\mathrm{\&rho;}}_{i}g$

当粒子物理状态计算完成后，要计算粒子的运动状态。

粒子的运动状态采用蛙跳式更新，如图3所示。

$u_t\&ap;\frac{u_{t-0.5\&CenterDot;\mathrm{\&Delta;t}}+u_{t+0.5\&CenterDot;\mathrm{\&Delta;t}}}{2}$

u_t+0.5·Δt  ＝u_t-0.5·Δt+Δta_t

r_t+Δt＝r_t+Δtu_t+0.5·t

Navier-Stokes方程计算的基础是利用核函数进行插值，而核函数计算的前提就是要知道周围有哪些粒子。本系统使用空间Hash方法，该方法能快速将自己定位到网格，继而可以搜索到周围网格的粒子，其结果如图4所示。

对于SPH边界处理方案，目前尚无一个好的公认的解决方法，大都通过模拟边界力和固体碰撞，来实现边界作用，本具体实施方式采用固体碰撞模型模拟固壁反弹作用，其结果如图5所示。

r_i＝r_i+dn

r_i＝cp

u_i＝u_i-(1+c_r)(u_i·n)n

其中，r_i是粒子i的位置，n是边界的法向，cp是碰撞点，c_r是恢复系数。

当粒子向边界碰撞时，边界对粒子的作用垂直于边界，作用效果不改变粒子的物理状态，而是改变粒子的运动状态，即粒子的位置和速度。

步骤12，通过OpenGL颜色混合模式和Billboard方法并根据每个流体粒子的粒子位置，计算出每个流体粒子的象素浓度，再通过Cuda将流体粒子的象素浓度进行平滑处理后输出。

具体的，浓度可以使流体看上去有深浅区别，本具体实施方式采用基于BillBoard颜色叠加的方法，将每个粒子的密度场看成一个逐渐减弱的亮度，密度场的叠加转换成的颜色的叠加，该叠加过程可以通过OpenGL颜色混合完成。具体为：在每个粒子位置使用OpenGL绘制一个矩形，并带有纹理贴图，该矩形始终朝向观察者，绘制结果会缓冲到FrameBuffer中，Cuda从FrameBuffer中取出结果进行光滑处理，将处理结果送CGShander模块进行合成。

步骤13，将流体粒子的密度和位置转化成一个空间密度场，通过设定一个阈值将空间密度场转换成Marching Cube数据构造，形成表面数据并得出表面的法向量。

具体的，本步骤是对流体粒子表面重建，包括：使用SPH方法模拟流体的结果是得到一系列的粒子位置和这些粒子的物理状态，表面重建就是将离散的粒子点，渲染成连续的流体，对粒子所表示的流体可以看成粒子表示的密度场，每个粒子对周围空间的密度的影响，从而形成了一个密度场，或者可以简单理解为浓度场。

在规则体数据场中，由八个数据点围成的小立方体称为体素，而这八个数据点就是该体素的角点。体素内的点，其属性值可以由八个角点的属性值进行三线性插值得到。八个数据点的属性决定了15种面的构造方式，这15种方式通过旋转等方式可以表示256个不同的构造。

步骤14，将表面数据和表面法向量处理成面中心和法向量格式，根据Snell原理，计算出光线的位置，在光线的位置通过光子效应绘制屏幕坐标。

具体的，Caustics效应是由于光的折射和反射造成的光的聚集现象，Caustics的计算可以基于Marching Cube表面结果进行，Marching Cube表面输出的数据组织结构是三角形的顶点和顶点处的法向量。使用一个三角形中心点和中心点的法向量就可以确定该面片的属性。如图6所示，入射光线I，经过p点和N所决定的平面发生折射和反射，而折射光线与接收平面（法向量是Ng）相交于P’，当大多数表面片元结果折射后汇聚于接收平面时，形成Caustics效应。在P’处同样使用OpenGL混合模式合成，将合成结果给Cuda做平滑处理，获得相应的屏幕坐标。

步骤15，根据象素流体粒子的象素浓度、屏幕坐标通过变换矩阵转换到世界坐标系下，并根据转换到世界坐标系下的象素流体粒子的象素浓度、屏幕坐标进行贴图和颜色合成，最后将合成的颜色输出。

具体的，本步骤采用的CG语言是一种可编程着色语言，功能强大，跨平台可以用来实现着色需求。在顶点编程中，将象素流体粒子的象素浓度、屏幕坐标通过变换矩阵转换到世界坐标系下，再根据转换到世界坐标系下的象素流体粒子的象素浓度、屏幕坐标进行贴图和颜色合成，最后将合成的颜色输出。

采用本具体实施方式提供的技术方案，通过Navier-Stokes状态方程获得流体粒子的粒子位置和粒子密度，无需依赖网格且计算速度较快、计算量和存储空间较小；借助于CudaGPU平台，进一步提升了运算速度；使用OpenGL混合和Billboard技术，使OpenGL可以自动生成浓度信息，避免了不必要的计算，同时加快了计算速度；在Marching Cube算法的基础上，利用光线的折射原理，计算面片的折射位置，通过OpenGL进行混合，能够获得较好的渲染效果。

需要说明的是，本领域普通技术人员可以理解实现上述各方法实施例中的全部或部分步骤是可以通过程序来指令相关的硬件完成，相应的程序可以存储于一种计算机存储介质中。

本发明的具体实施方式还提供了一种3D流体的仿真系统，如图7所示，包括：

粒子状态计算模块71，用于通过空间Hash方法搜索每个流体粒子的邻居粒子，并通过核函数解出Navier-Stokes状态方程，获得流体粒子的粒子位置和粒子密度，并输出；

粒子浓度计算模块72，用于通过OpenGL颜色混合模式和Billboard方法并根据每个流体粒子的粒子位置，计算出每个流体粒子的象素浓度，再通过Cuda将流体粒子的象素浓度进行平滑处理后输出；

Marching Cube表面重建模块73，用于将流体粒子的密度和位置转化成一个空间密度场，通过设定一个阈值将空间密度场转换成Marching Cube数据构造，形成表面数据并得出表面的法向量；

Caustics效应计算模块74，用于将表面数据和表面法向量处理成面中心和法向量格式，根据Snell原理，计算出光线的位置，在光线的位置通过光子效应绘制屏幕坐标；

CGShader着色模块75，用于根据象素流体粒子的象素浓度、屏幕坐标通过变换矩阵转换到世界坐标系下，并根据变换矩阵转换到世界坐标系下的象素流体粒子的象素浓度、屏幕坐标进行贴图和颜色合成，最后将合成的颜色输出。

可选的，在粒子状态计算模块71中可以包括：边界作用子模块，用于在通过核函数解出Navier-Stokes状态方程后，对流体粒子施加边界作用。

上述3D流体的仿真系统中包含的各模块的处理功能的具体实现方式在之前的方法实施例中已经描述，在此不再重复描述。

采用本具体实施方式提供的技术方案，通过Navier-Stokes状态方程获得流体粒子的粒子位置和粒子密度，无需依赖网格且计算速度较快、计算量和存储空间较小；借助于Cuda的GPU平台，进一步提升了运算速度；使用OpenGL混合和Bil lboard技术，使OpenGL可以自动生成浓度信息，避免了不必要的计算，同时加快了计算速度；在Marching Cube算法的基础上，利用光线的折射原理，计算面片的折射位置，通过OpenGL进行混合，能够获得较好的渲染效果。

需要注意的是，上述系统实施例中，所包括的各个模块只是按照功能逻辑进行划分的，但并不局限于上述的划分，只要能够实现相应的功能即可；另外，各功能模块的具体名称也只是为了便于相互区分，并不用于限制本发明的保护范围。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明实施例揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (4)

1.一种3D流体的仿真方法，其特征在于，包括：

通过空间Hash方法搜索每个流体粒子的邻居粒子，并通过核函数解出Navier-Stokes状态方程，获得流体粒子的粒子位置和粒子密度，并输出；

通过OpenGL颜色混合模式和Billboard方法并根据每个流体粒子的粒子位置，计算出每个流体粒子的象素浓度，再通过Cuda将所述流体粒子的象素浓度进行平滑处理后输出；

将所述流体粒子的密度和位置转化成一个空间密度场，通过设定一个阈值将空间密度场转换成Marching Cube数据构造，形成表面数据并得出表面的法向量；

将所述表面数据和表面法向量处理成面中心和法向量格式，根据Snell原理，计算出光线的位置，在所述光线的位置通过光子效应绘制屏幕坐标；

根据所述象素流体粒子的象素浓度、屏幕坐标通过变换矩阵转换到世界坐标系下，并根据变换矩阵转换到世界坐标系下的所述象素流体粒子的象素浓度、屏幕坐标进行贴图和颜色合成，最后将合成的颜色输出。

2.根据权利要求1所述的3D流体的仿真方法，其特征在于，在通过核函数解出Navier-Stokes状态方程后还包括：对所述流体粒子施加边界作用。

3.一种3D流体的仿真系统，其特征在于，包括：

粒子状态计算模块，用于通过空间Hash方法搜索每个流体粒子的邻居粒子，并通过核函数解出Navier-Stokes状态方程，获得流体粒子的粒子位置和粒子密度，并输出；

粒子浓度计算模块，用于通过OpenGL颜色混合模式和Billboard方法并根据每个流体粒子的粒子位置，计算出每个流体粒子的象素浓度，再通过Cuda将所述流体粒子的象素浓度进行平滑处理后输出；

Marching Cube表面重建模块，用于将所述流体粒子的密度和位置转化成一个空间密度场，通过设定一个阈值将空间密度场转换成Marching Cube数据构造，形成表面数据并得出表面的法向量；

Caustics效应计算模块，用于将所述表面数据和表面法向量处理成面中心和法向量格式，根据Snell原理，计算出光线的位置，在所述光线的位置通过光子效应绘制屏幕坐标；

CGShader着色模块，用于根据所述象素流体粒子的象素浓度、屏幕坐标通过变换矩阵转换到世界坐标系下，并根据变换矩阵转换到世界坐标系下的所述象素流体粒子的象素浓度、屏幕坐标进行贴图和颜色合成，最后将合成的颜色输出。

4.根据权利要求3所述的3D流体的仿真系统，其特征在于，在粒子状态计算模块中包括：

边界作用子模块，用于在通过核函数解出Navier-Stokes状态方程后，对所述流体粒子施加边界作用。

Images