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

来源:互联网 发布:string数组去掉空格 编辑:程序博客网 时间:2024/04/26 20:41

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.
1 0
原创粉丝点击