【CG物理模拟系列】流体模拟--粒子法之SPH(理论)

原创 2016年05月31日 13:46:30

SPH法介绍


SPH(Smoothed Particle Hydrodynamics)是早年银河系碰撞,天体形成等宇宙物理学模拟所使用的方法[Lucy1977],近年来被应用到流体,热等其他现象处理中。但是,由于是曾经宇宙物理学中常产生的压缩性流体问题的处理方法,并不适用于像水和流速较慢的空气(流速<<音速)这样的非压缩性流体 (强行转换非压缩性的研究近年来也有很多[Becker2007,Solenthaler2009])。

[Desbrun1996]将SPH法第一次应用于CG领域, 此后[Muller2003]将其改进,使其也能适用于3次元的粘性流体。 以下,主要对[Muller2003]的方法进行详细介绍。
SPH法的离散化方程

SPH法中的物理量sph_eq_phi.gif的离散化方程如下所示。
sph_eq_01.gif

这里,N是临近粒子的集合,m是粒子质量,sph_eq_rho.gif是粒子密度, W是kernel方程。内核方程的有效半径是h,有效半径外值为0(Compact Support )。 同时,它的积分设定为1(sph_eq_02.gif)

根据上述方程,通过周围粒子的权重相加和来近似计算物理量。 另外,物理量的梯度sph_eq_nabla_phi.gif可以使用内核函数的导数来表示。 物理量sph_eq_phi.gif的梯度公式为,
sph_eq_05.gif

物理量sph_eq_phi.gif的laplacian算子为,
sph_eq_06.gif


使用上述方程,流体的密度方程如下所示。
sph_eq_07.gif

SPH法中NS方程(Navier-Stokes equations)的解法

支配方程(非压缩性的navier-stokes方程)如下所示
sph_eq_08.gif
sph_eq_09.gif

这里,sph_eq_Vu.gif是流体速度,sph_eq_10.gif是动粘性系数,sph_eq_rho.gif是流体的密度,p为圧力,sph_eq_Vf_ext.gif是外力,包含重力,表面张力等。 式一是质量守恒方程,式二为动量守恒方程。

先将支配方程用粒子离散化后,通过SPH法求解。为了将粒子表示为液体,只要粒子质量保持不变就可以保证质量的守恒性,因此,不必求解质量守恒方程(这里仅质量守恒,体积不守恒)。质量保存式パーティクル質量が変化しないかぎり質量保存性が常に保持され,質量保存式を解く必要がない(あくまで質量保存であり,体積保存ではないことに注意).同时,粒子法不必像网格法一样计算平流项,只需移动粒子即可。为此,不必离散化平流项sph_eq_11.gif。所以,各粒子i的速度sph_eq_Vu_i.gif的计算公式更新如下。
sph_eq_12.gif
sph_eq_13.gif

sph_eq_Vf_i.gif的各项随后叙述。

粘性扩散项 


粘性项根据流体的速度构成。速度sph_eq_Vu.giflaplacian算子离散化后( 把离散化方程sph_eq_phi.gifsph_eq_Vu.gif替换),由于粒子i到j 所施加的力与j到i所施加的力不同(asymmetric) 。 不做处理的话,速度快的粒子对速度慢的粒子有影响,反之,影响力小的粒子的速度则会发散。

粘性力便是由流体的速度差来决定的力,通过下式即可实现symmetric的作用。
sph_eq_14.gif

圧力项


圧力项是由压力值的梯度组成。通过使用梯度值用的离散化方程,虽然能离散化压力项,但是也会发生同粘性扩散项相同的asymmetric问题,根据[Muller2003],
sph_eq_15.gif

可以使用压力的平均值来解决。

sph_eq_p_i.gif可以通过表示压力和密度关系的状态方程计算。理想气体的状态方程为,
sph_eq_16a.gif

这里,k是气体常数。大多数的SPH模拟,使用的是[Desbrun1996]提出的公式。
sph_eq_16.gif

这里,sph_eq_rho_0.gif是模拟对象的流体密度。例如,水的密度为sph_eq_17.gif,编程时,初期设定使用最大密度即可。这个公式简单易懂,但会导致流体压缩。非压缩性强制转换中最常用的方法是 求解Poisson 方程sph_eq_18.gif,但为此也要求解一个庞大的线形系统,非常耗时。

[Becker2007]使用Tait方程,可以确保一定程度的非压缩性。Tait方程如下.
sph_eq_19.gif

这里,sph_eq_20.gif


B是圧力定数.
sph_eq_21.gif

这里,sph_eq_c_s.gif是音速。通过 sph_eq_22.gif,从sph_eq_23.gif的关系中计算得出。 密度的変化率sph_eq_eta.gif常用0.01表示。

这个方法保有很高的非压缩性(1%以下),为了使计算更安定,timestep的值要设定的很小才可以。timestep的大小可以通过CFL(Courant-Friedrichs-Lewy)计算
sph_eq_24.gif

这里,sph_eq_Vf_a.gif是外力,sph_eq_alpha.gif是粘性定数,根据[Becker2007]取值范围为0.08~0.5。

参考文献 

[Becker2007] M. Becker and M. Teschner, Weakly Compressible SPH for Free Furface Flows, In Proc. SCA2007, pp.209-217, 2007.

[Desbrun1996] M. Desbrun and M.-P. Cani, Smoothed Particles: A New Paradigm for Animating Highly Deformable Bodies, Eurographics Workshop on Computer Animation and Simulation (EGCAS), pp.61-76. 1996.

[Lucy1977] L. B. Lucy, A Numerical Approach to the Testing of the Fission Hypothesis, The Astronomical Journal, Vol.82, No.12, pp.1013-1024, 1977.

[Muller2003] M. Muller, D. Charypar and M. Gross, Particle-based Fluid Simulation for Interactive Applications, In Proc. SCA2003, pp.154-159, 2003.

[Solenthaler2009] B. Solenthaler and R. Pajarola, Predictive-Corrective Incompressible SPH, In Proc. SIGGRAPH2009, pp.1-6, 2009.
版权声明:本文为博主原创文章,未经博主允许不得转载。

HTML5实验:JavaScript模拟流体效果

把现实世界当中的物体模拟到计算机当中,一些简单的物理实验、碰撞旋转等等难度还是不算很大,难度较大的应当算流体模拟。   本文将在Canvas当中模拟出一个2D平面内的水珠,涉及的知识点和技巧包括:J...

【CG物理模拟系列】粒子法--表面生成手法(上)

使用粒子法模拟水流时,水流表面生成方法主要分为以下两类。 表面隐函数 Muller的色彩函数 [1]Zhu and Bridson的方法 [2]Adams等的方法 [3]Anisotropic ...

【CG物理模拟系列】流体模拟--粒子法之SPH(代码讲解)

WCSPH,PCISPH,IISPH等研究方法,其本质都是以非压缩性为目标,求解Navier-Stokes方程。 本文以WCSPH为例,讲解下SPH方法代码的实现。 代码讲解 sph_type...

【CG物理模拟系列】流体模拟--粒子法之SPH(实现)

邻域搜索的效率化 SPH等粒子法,由于需要考虑到邻域粒子带来的影响,通常邻域搜索都会消耗大量时间。如果我们只是单纯的计算所有粒子组合的欧氏距离的话,计算时间只会呈指数增加。 而空间分割法的...

【CG物理模拟系列】流体模拟--粒子法之SPH法的加权函数计算

SPH法的加权函数 Poly6 kernel[1] 用于密度计算等。由于r只存在于2次项中,可以省去平方根的计算。 梯度为, 拉普拉斯算子为, 如图所示, ...

【CG物理模拟系列】流体模拟--粒子法之MPS法(理论)

MPS法  前面的文章里我们讲过SPH曾是为了处理压缩性流体问题而提出的方法,与之相对,这一篇来说说用粒子法处理非压缩性流体的研究方法--Moving Particle Semi-implicit...

【CG物理模拟系列】粒子法--表面生成手法(下)

这一篇来说说 网格生成方法 中的 Screen Space Mesh 法。 Screen Space Mesh  一般情况下,从隐函数曲面中提取出记载着表面数据的三角形网格面时,我们常用像Ma...

基于SPH的流体模拟实践和一些技巧总结

SPH的流体模拟是目前大多数游戏所采用的模拟流体方法,特点是简单,十分容易实现,相比与基于Grid的Eulerian方法更加简单和高速,本文主要介绍一下使用SPH的流体模拟中一些常用的技巧和数据结构。...

【CG物理模拟系列】弹性体模拟--Position-based法之Shape Matching

Position-based法 Position-based法与传统的力学基础方法不同,根据构成物体的顶点等元素间的约束条件(Constraint),直接变更顶点的位置坐标的方法。        ...

【CG物理模拟系列】开篇:介绍(下)

上一篇介绍了CG物体模拟的定义,流程及种类,这一篇讲下物理模拟常用手法,物理模拟引擎,从物理模拟+3DCG程序的编写・到导出结果动画的处理顺序。 物理模拟常用手法 粒子法(Particle Me...
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:【CG物理模拟系列】流体模拟--粒子法之SPH(理论)
举报原因:
原因补充:

(最多只允许输入30个字)