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

WCSPH,PCISPH,IISPH等研究方法，其本质都是以非压缩性为目标，求解Navier-Stokes方程。

本文以WCSPH为例，讲解下SPH方法代码的实现。

代码讲解

sph_type.h里定义几个vector函数类型

#ifndef _SPHTYPE_H#define _SPHTYPE_Htypedef unsigned  int uint;struct Vec3_f{float x;float y;float z;};struct Vec3_i{int x;int y;int z;};struct uint3{uint x;uint y;uint z;};#endif  //_SPHTYPE_H

sph_timer.h是一个时间循环函数，实现如下

#ifndef _SPHTIMER_H#define _SPHTIMER_H#include <windows.h>class Timer{private:int frames;int update_time;int last_time;double FPS;public:Timer();void update();double get_fps();};#endif//_SPHTIMER_H

这是sph_timer.cpp的具体内容

#include "sph_timer.h"Timer::Timer(){frames=0;update_time=1000;last_time=0;FPS=0;}void Timer::update(){frames++;if(GetTickCount()-last_time > update_time){FPS=((double)frames/(double)(GetTickCount()-last_time))*1000.0;last_time=GetTickCount();frames=0;}}double Timer::get_fps(){return FPS;}

sph_header.h包含了一些要用到的头文件

#ifndef _SPHHEADER_H#define _SPHHEADER_H#include <stdio.h>#include <stdlib.h>#include <math.h>#include <time.h>//#include "rx_utility.h"#define PI 3.141592f#define INF 1E-12f#define BOUNDARY 0.0001f#endif

sph_data.h包含了一些关于窗口的定义

#ifndef _SPHDATA_H#define _SPHDATA_H#include "sph_header.h"#include "sph_type.h"float window_width=1000;float window_height=1000;float xRot = 15.0f;float yRot = 0.0f;float xTrans = 0.0f;float yTrans = 0.0f;float zTrans = -35.0;int ox;int oy;int buttonState;float xRotLength = 0.0f;float yRotLength = 0.0f;Vec3_f real_world_origin;Vec3_f real_world_side;Vec3_f sim_ratio;//size of the worldfloat world_witdth;float world_height;float world_length;#endif //_SPHDATA_H

以上是这个程序要用的的一些窗口及数据设置，如果大家有自己的工具库和函数框架，直接修改下文相关变量。

----------------------------------------------------------------------------------------------------------

----------------------------------------------------------------------------------------------------------

本番来了！！！

sph_system.h中定义了一个particle()类，包含粒子的id，速度，加速度，密度，压力等信息。SPHSystem()类中包含了一些函数定义，及组成SPH算法的共有私有函数类，作用见注释。

#ifndef __SPHSYSTEM_H__#define __SPHSYSTEM_H__#include "sph_type.h"//粒子类，包含粒子的各项属性class Particle{public:uint id;float3 pos;//位置float3 vel;//速度float3 acc;//加速度float3 ev;float dens;float pres;float surf_norm;Particle *next;};class SPHSystem{public:uint max_particle;uint num_particle;float kernel;float mass;float3 world_size;float cell_size;//单位空间大小uint3 grid_size;//网格大小uint tot_cell;//单元总数float3 gravity;float wall_damping;float rest_density;float gas_constant;//气体定数float viscosity;//粘性float time_step;float surf_norm;float surf_coe;//各项系数float poly6_value;float spiky_value;float visco_value;float grad_poly6;float lplc_poly6;float kernel_2;float self_dens;float self_lplc_color;Particle *mem;Particle **cell;uint sys_running;public:SPHSystem();~SPHSystem();void update();void init_system();//初始化void add_particle(float3 pos, float3 vel);//添加粒子private:void build_table();//创建表格void comp_dens_pres();//计算密度，压力void comp_force_adv();//计算压力项和粘性项void advection();//计算移流private:int3 calc_cell_pos(float3 p);//计算粒子p所属的单元uint calc_cell_hash(int3 cell_pos);//计算单元的hash值};#endif

sph_system.cpp是其具体实现，也是这个工程的主函数。

#include "sph_system.h"#include "sph_header.h"SPHSystem::SPHSystem(){max_particle=30000;num_particle=0;kernel=0.04f;mass=0.02f;//初始化空间大小，并计算网格sizeworld_size.x=0.64f;world_size.y=0.64f;world_size.z=0.64f;cell_size=kernel;grid_size.x=(uint)ceil(world_size.x/cell_size);grid_size.y=(uint)ceil(world_size.y/cell_size);grid_size.z=(uint)ceil(world_size.z/cell_size);tot_cell=grid_size.x*grid_size.y*grid_size.z;gravity.x=0.0f; gravity.y=-6.8f;gravity.z=0.0f;wall_damping=-0.5f;rest_density=1000.0f;//初始密度gas_constant=1.0f;//气体系数viscosity=6.5f;time_step=0.003f;surf_norm=6.0f;surf_coe=0.1f;poly6_value=315.0f/(64.0f * PI * pow(kernel, 9));;spiky_value=-45.0f/(PI * pow(kernel, 6));visco_value=45.0f/(PI * pow(kernel, 6));grad_poly6=-945/(32 * PI * pow(kernel, 9));lplc_poly6=-945/(8 * PI * pow(kernel, 9));kernel_2=kernel*kernel;self_dens=mass*poly6_value*pow(kernel, 6);self_lplc_color=lplc_poly6*mass*kernel_2*(0-3/4*kernel_2);mem=(Particle *)malloc(sizeof(Particle)*max_particle);cell=(Particle **)malloc(sizeof(Particle *)*tot_cell);sys_running=0;}SPHSystem::~SPHSystem(){free(mem);free(cell);}void SPHSystem::update(){if(sys_running == 0){return;}build_table();comp_dens_pres();comp_force_adv();advection();}void SPHSystem::init_system(){float3 pos;float3 vel;vel.x=0.0f;vel.y=0.0f;vel.z=0.0f;for(pos.x=world_size.x*0.2f; pos.x<world_size.x*0.8f; pos.x+=(kernel*0.5f)){for(pos.y=world_size.y*0.5f; pos.y<world_size.y*0.9f; pos.y+=(kernel*0.5f)){for(pos.z=world_size.z*0.2f; pos.z<world_size.z*0.8f; pos.z+=(kernel*0.5f)){add_particle(pos, vel);}}}printf("Init Particle: %u\n", num_particle);}//添加粒子void SPHSystem::add_particle(float3 pos, float3 vel){Particle *p=&(mem[num_particle]);p->id=num_particle;p->pos=pos;p->vel=vel;p->acc.x=0.0f;p->acc.y=0.0f;p->acc.z=0.0f;p->ev.x=0.0f;p->ev.y=0.0f;p->ev.z=0.0f;p->dens=rest_density;p->pres=0.0f;p->next=NULL;num_particle++;}void SPHSystem::build_table(){Particle *p;uint hash;for(uint i=0; i<tot_cell; i++){cell[i]=NULL;}//创建hash表for(uint i=0; i<num_particle; i++){p=&(mem[i]);hash=calc_cell_hash(calc_cell_pos(p->pos));if(cell[hash] == NULL){p->next=NULL;cell[hash]=p;}else{p->next=cell[hash];cell[hash]=p;}}}void SPHSystem::comp_dens_pres(){Particle *p;//当前粒子Particle *np;//周围粒子int3 cell_pos;//粒子所属cell坐标int3 near_pos;//周围粒子坐标uint hash;float3 rel_pos;float r2;for(uint i=0; i<num_particle; i++){p=&(mem[i]); cell_pos=calc_cell_pos(p->pos);p->dens=0.0f;p->pres=0.0f;for(int x=-1; x<=1; x++){for(int y=-1; y<=1; y++){for(int z=-1; z<=1; z++){near_pos.x=cell_pos.x+x;near_pos.y=cell_pos.y+y;near_pos.z=cell_pos.z+z;hash=calc_cell_hash(near_pos);if(hash == 0xffffffff){continue;}np=cell[hash];//扫描kernel_2范围内的粒子while(np != NULL){rel_pos.x=np->pos.x-p->pos.x;rel_pos.y=np->pos.y-p->pos.y;rel_pos.z=np->pos.z-p->pos.z;r2=rel_pos.x*rel_pos.x+rel_pos.y*rel_pos.y+rel_pos.z*rel_pos.z;if(r2<INF || r2>=kernel_2){np=np->next;continue;}p->dens=p->dens + mass * poly6_value * pow(kernel_2-r2, 3);np=np->next;}}}}//密度及压力p->dens=p->dens+self_dens;p->pres=(pow(p->dens / rest_density, 7) - 1) *gas_constant;//这里使用了Tait方程，详见理论那篇文章}}void SPHSystem::comp_force_adv(){Particle *p;//当前粒子Particle *np;//周围粒子int3 cell_pos;//所属cell位置int3 near_pos;//检测周围粒子位置uint hash;float3 rel_pos;float3 rel_vel;float r2;float r;float kernel_r;float V;float pres_kernel;float visc_kernel;float temp_force;float3 grad_color;float lplc_color;for(uint i=0; i<num_particle; i++){p=&(mem[i]); cell_pos=calc_cell_pos(p->pos);p->acc.x=0.0f;p->acc.y=0.0f;p->acc.z=0.0f;grad_color.x=0.0f;grad_color.y=0.0f;grad_color.z=0.0f;lplc_color=0.0f;for(int x=-1; x<=1; x++){for(int y=-1; y<=1; y++){for(int z=-1; z<=1; z++){near_pos.x=cell_pos.x+x;near_pos.y=cell_pos.y+y;near_pos.z=cell_pos.z+z;hash=calc_cell_hash(near_pos);if(hash == 0xffffffff){continue;}np=cell[hash];while(np != NULL){rel_pos.x=p->pos.x-np->pos.x;rel_pos.y=p->pos.y-np->pos.y;rel_pos.z=p->pos.z-np->pos.z;r2=rel_pos.x*rel_pos.x+rel_pos.y*rel_pos.y+rel_pos.z*rel_pos.z;if(r2 < kernel_2 && r2 > INF){r=sqrt(r2);V=mass/np->dens/2;kernel_r=kernel-r;//计算压力项pres_kernel=spiky_value * kernel_r * kernel_r;temp_force=V * (p->pres+np->pres) * pres_kernel;p->acc.x=p->acc.x-rel_pos.x*temp_force/r;p->acc.y=p->acc.y-rel_pos.y*temp_force/r;p->acc.z=p->acc.z-rel_pos.z*temp_force/r;rel_vel.x=np->ev.x-p->ev.x;rel_vel.y=np->ev.y-p->ev.y;rel_vel.z=np->ev.z-p->ev.z;//计算粘性项visc_kernel=visco_value*(kernel-r);temp_force=V * viscosity * visc_kernel;p->acc.x=p->acc.x + rel_vel.x*temp_force; p->acc.y=p->acc.y + rel_vel.y*temp_force; p->acc.z=p->acc.z + rel_vel.z*temp_force; float temp=(-1) * grad_poly6 * V * pow(kernel_2-r2, 2);grad_color.x += temp * rel_pos.x;grad_color.y += temp * rel_pos.y;grad_color.z += temp * rel_pos.z;lplc_color += lplc_poly6 * V * (kernel_2-r2) * (r2-3/4*(kernel_2-r2));}np=np->next;}}}}lplc_color+=self_lplc_color/p->dens;p->surf_norm=sqrt(grad_color.x*grad_color.x+grad_color.y*grad_color.y+grad_color.z*grad_color.z);if(p->surf_norm > surf_norm){p->acc.x+=surf_coe * lplc_color * grad_color.x / p->surf_norm;p->acc.y+=surf_coe * lplc_color * grad_color.y / p->surf_norm;p->acc.z+=surf_coe * lplc_color * grad_color.z / p->surf_norm;}}}void SPHSystem::advection(){Particle *p;for(uint i=0; i<num_particle; i++){p=&(mem[i]);//速度及位置更新p->vel.x=p->vel.x+p->acc.x*time_step/p->dens+gravity.x*time_step;p->vel.y=p->vel.y+p->acc.y*time_step/p->dens+gravity.y*time_step;p->vel.z=p->vel.z+p->acc.z*time_step/p->dens+gravity.z*time_step;p->pos.x=p->pos.x+p->vel.x*time_step;p->pos.y=p->pos.y+p->vel.y*time_step;p->pos.z=p->pos.z+p->vel.z*time_step;//墙壁if(p->pos.x >= world_size.x-BOUNDARY){p->vel.x=p->vel.x*wall_damping;p->pos.x=world_size.x-BOUNDARY;}if(p->pos.x < 0.0f){p->vel.x=p->vel.x*wall_damping;p->pos.x=0.0f;}if(p->pos.y >= world_size.y-BOUNDARY){p->vel.y=p->vel.y*wall_damping;p->pos.y=world_size.y-BOUNDARY;}if(p->pos.y < 0.0f){p->vel.y=p->vel.y*wall_damping;p->pos.y=0.0f;}if(p->pos.z >= world_size.z-BOUNDARY){p->vel.z=p->vel.z*wall_damping;p->pos.z=world_size.z-BOUNDARY;}if(p->pos.z < 0.0f){p->vel.z=p->vel.z*wall_damping;p->pos.z=0.0f;}p->ev.x=(p->ev.x+p->vel.x)/2;p->ev.y=(p->ev.y+p->vel.y)/2;p->ev.z=(p->ev.z+p->vel.z)/2;}}//计算粒子p所属的cell位置int3 SPHSystem::calc_cell_pos(float3 p){int3 cell_pos;cell_pos.x = int(floor((p.x) / cell_size));cell_pos.y = int(floor((p.y) / cell_size));cell_pos.z = int(floor((p.z) / cell_size));    return cell_pos;}//cell hash值计算uint SPHSystem::calc_cell_hash(int3 cell_pos){if(cell_pos.x<0 || cell_pos.x>=(int)grid_size.x || cell_pos.y<0 || cell_pos.y>=(int)grid_size.y || cell_pos.z<0 || cell_pos.z>=(int)grid_size.z){return (uint)0xffffffff;}cell_pos.x = cell_pos.x & (grid_size.x-1);      cell_pos.y = cell_pos.y & (grid_size.y-1);  cell_pos.z = cell_pos.z & (grid_size.z-1);  return ((uint)(cell_pos.z))*grid_size.y*grid_size.x + ((uint)(cell_pos.y))*grid_size.x + (uint)(cell_pos.x);}

没有绘制水面，只是实现了最基本的SPH算法原理，以供参考。
以上~~

参考文献 

Matthias Müller, David Charypar, and Markus Gross,  Particle-Based Fluid Simulation for Interactive Applications, SCA 2003.