流体模拟（二）

SPH算法实现2：

在前面一节我们完成了粒子缓存类，网格类和邻接表类。我们现在可以正式的整合在我们的流体系统类中了。

流体系统类

    class FluidSystem{public:FluidSystem();void init(unsigned short maxPointCounts,const glm::vec3& wallBox_min,const glm::vec3& wallBox_max,const glm::vec3& initFluidBox_min,const glm::vec3& initFluidBox_max,const glm::vec3& gravity){_init(maxPointCounts, fBox3(wallBox_min,wallBox_max),fBox3(initFluidBox_min,initFluidBox_max), gravity);}//获取点的尺寸（字节）unsigned int getPointStride()const {return sizeof(Point);}//获取点的数量unsigned int getPointCounts()const{return m_pointBuffer.size();}//获取流体点缓存const Point* getPointBuf()const{return m_pointBuffer.get(0);}//逻辑桢void tick();~FluidSystem();private://初始化系统void _init(unsigned short maxPointCounts,const fBox3& wallBox,const fBox3& initFluidBox, const glm::vec3& gravity);//计算密度，压力以及相邻关系void _computerPressure();//计算加速度void _computerForce();//移动粒子void _advance();//创建初始液体块void _addFluidVolume(const fBox3& fluidBox,float spacing);//数据成员PointBuffer m_pointBuffer;GridContainer m_gridContainer;NeighborTable m_neighborTable;//点位置缓存数据(x,y,z)std::vector<float>posData;//SPH光滑核float m_kernelPoly6;float m_kernelSpiky;float m_kernelViscosity;//其他参数float m_pointDistance;//半径float m_unitScale;//尺寸单位float m_viscosity;//粘性float m_restDensity;//静态密度float m_pointMass;//质量float m_smoothRadius;//光滑核半径float m_gasConstantK;//气体常量kfloat m_boundaryStiffness;//边界刚性float m_boundaryDampening;//边界阻尼float m_speedLimiting;//速度限制glm::vec3 m_gravityDir;//重力矢量int m_rexSize[3];//网格尺寸fBox3 m_sphWallBox;};

目前流体系统还是比较简单的，设定我们SPH需要的各种属性，并建立我们上节创建的三种类成员。算法主要在初始化函数_init,密度和压力计算函数_computerPressure,加速度计算函数_computerForce和移动粒子函数_advance。我们的动画在tick帧函数中进行。最后我们计算后的粒子位置在每帧计算后更新posData缓存数据。

然后是我们的实现代码：

 FluidSystem::FluidSystem(){m_unitScale            = 0.004f;            // 尺寸单位m_viscosity            = 1.0f;                // 粘度m_restDensity        = 1000.f;            // 密度m_pointMass            = 0.0006f;            // 粒子质量m_gasConstantK        = 1.0f;                // 理想气体方程常量m_smoothRadius        = 0.01f;            // 光滑核半径m_pointDistance       =0.0;m_rexSize[0]=0;m_rexSize[1]=0;m_rexSize[2]=0;m_boundaryStiffness = 10000.f;m_boundaryDampening = 256;m_speedLimiting        =    200;//Poly6 Kernelm_kernelPoly6 = 315.0f/(64.0f * 3.141592f * pow(m_smoothRadius, 9));//Spiky Kernelm_kernelSpiky = -45.0f/(3.141592f * pow(m_smoothRadius, 6));//Viscosity Kernelm_kernelViscosity = 45.0f/(3.141592f * pow(m_smoothRadius, 6));}FluidSystem::~FluidSystem(){}//构造流体中的点void FluidSystem::_addFluidVolume(const fBox3 &fluidBox, float spacing){for(float z=fluidBox.max.z; z>=fluidBox.min.z; z-=spacing){for(float y=fluidBox.min.y; y<=fluidBox.max.y; y+=spacing){for(float x=fluidBox.min.x; x<=fluidBox.max.x; x+=spacing){Point* p = m_pointBuffer.addPointReuse();p->pos=glm::vec3(x,y,z);}}}}void FluidSystem::_init(unsigned short maxPointCounts, const fBox3 &wallBox, const fBox3 &initFluidBox, const glm::vec3 &gravity){m_pointBuffer.reset(maxPointCounts);m_sphWallBox=wallBox;m_gravityDir=gravity;m_pointDistance=pow(m_pointMass/m_restDensity, 1.0/3.0);//计算粒子间距_addFluidVolume(initFluidBox, m_pointDistance/m_unitScale);m_gridContainer.init(wallBox, m_unitScale, m_smoothRadius*2.0, 1.0,m_rexSize);//设置网格尺寸(2r)posData=std::vector<float>(3*m_pointBuffer.size(),0);}

在构造函数里，我们给属性赋初值，并根据公式3.6:

公式3.12:

公式3.17：

计算该三者的光滑核系数。_addFluidVolume函数，我们用box3创建流体粒子块，并计算粒子位置然，后都加入粒子缓存。在_init函数里，我们初始化流体演示的容器，重力。然后根据粒子质量和密度计算出半径（即粒子间距），然后加入流体块，初始化空间网格。最后初始化posData的粒子位置缓存。

然后我们看计算压力，密度以及相邻关系代码：

void FluidSystem::_computerPressure(){float h2=m_smoothRadius*m_smoothRadius;//h^2m_neighborTable.reset(m_pointBuffer.size());//重置邻接表for(unsigned int i=0;i<m_pointBuffer.size();i++){Point* pi=m_pointBuffer.get(i);float sum=0.0;m_neighborTable.point_prepare(i);int gridCell[8];m_gridContainer.findCells(pi->pos, m_smoothRadius/m_unitScale, gridCell);for(int cell=0;cell<8;cell++){if (gridCell[cell]==-1) continue;int pndx=m_gridContainer.getGridData(gridCell[cell]);while (pndx!=-1) {Point* pj=m_pointBuffer.get(pndx);if(pj==pi)sum+=pow(h2,3.0);else{glm::vec3 pi_pj=(pi->pos-pj->pos)*m_unitScale;float r2=pi_pj.x*pi_pj.x+pi_pj.y*pi_pj.y+pi_pj.z*pi_pj.z;if(h2>r2){float h2_r2=h2-r2;sum+=pow(h2_r2, 3.0);if(!m_neighborTable.point_add_neighbor(pndx, glm::sqrt(r2)))goto NEIGHBOR_FULL;}}pndx=pj->next;}}NEIGHBOR_FULL:m_neighborTable.point_commit();pi->density=m_kernelPoly6*m_pointMass*sum;pi->pressure=(pi->density-m_restDensity)*m_gasConstantK;}}

在计算邻接关系之前我们先重置邻接表。然后遍历粒子缓存中所有的粒子。每次遍历时首先获取粒子pi，并做初始化处理，之后我们设置一个数组gridCell去获取当前粒子的周围网格，用于我们处理邻域粒子。在每个网格内，我们通过next遍历该网格内的粒子直到next为-1。我们获取距离在光滑核半径内的领域粒子，通过公式3.7:

计算最终pi粒子的密度。然后通过公式3.10:

计算最终pi粒子受到的压力。最后将邻接粒子信息存入邻接表。

我们再来看计算加速度的代码实现：

    void FluidSystem::_computerForce(){float h2=m_smoothRadius*m_smoothRadius;for(unsigned int i=0;i<m_pointBuffer.size();i++){Point* pi=m_pointBuffer.get(i);glm::vec3 accel_sum=glm::vec3(0.0);int neighborCounts=m_neighborTable.getNeighborCounts(i);for(unsigned int j=0;j<neighborCounts;j++){unsigned short neighborIndex;float r;m_neighborTable.getNeighborInfo(i, j, neighborIndex, r);Point* pj=m_pointBuffer.get(neighborIndex);glm::vec3 ri_rj=(pi->pos-pj->pos)*m_unitScale;float h_r=m_smoothRadius-r;float pterm=-m_pointMass*m_kernelSpiky*h_r*h_r*(pi->pressure+pj->pressure)/(2.f * pi->density * pj->density);accel_sum+=ri_rj*pterm/r;float vterm=m_kernelViscosity*m_viscosity*h_r*m_pointMass/(pi->density * pj->density);accel_sum+=(pj->velocity_eval - pi->velocity_eval)*vterm;}pi->accel=accel_sum;}}

在计算加速度时，我们同样遍历粒子缓存中的每个粒子，并通过该粒子索引访问邻接表内的邻域粒子，通过对每个邻域粒子计算公式3.19:

这里的重力g我们在移动粒子中加入。这样遍历完所有的邻域粒子，我们的最终加速度也就计算完成了。

接下来便可以看粒子移动的代码：

void FluidSystem::_advance(){float deltaTime=0.003;float SL2=m_speedLimiting*m_speedLimiting;//posData.clear();for (unsigned int i=0; i<m_pointBuffer.size(); i++) {Point* p=m_pointBuffer.get(i);glm::vec3 accel=p->accel;//获取p的加速度float accel2=accel.x*accel.x+accel.y*accel.y+accel.z*accel.z;if(accel2>SL2)//速度限制accel*=m_speedLimiting/glm::sqrt(accel2);float diff;//边界情况// Z方向边界diff = 1 * m_unitScale - (p->pos.z - m_sphWallBox.min.z)*m_unitScale;if (diff > 0.0f ){glm::vec3 norm(0, 0, 1.0);float adj = m_boundaryStiffness * diff - m_boundaryDampening * glm::dot ( norm,p->velocity_eval );accel.x += adj * norm.x;accel.y += adj * norm.y;accel.z += adj * norm.z;}diff = 1 * m_unitScale - (m_sphWallBox.max.z - p->pos.z)*m_unitScale;if (diff > 0.0f){glm::vec3 norm( 0, 0, -1.0);float adj = m_boundaryStiffness * diff - m_boundaryDampening *glm::dot ( norm,p->velocity_eval );accel.x += adj * norm.x;accel.y += adj * norm.y;accel.z += adj * norm.z;}//X方向边界diff = 1 * m_unitScale - (p->pos.x - m_sphWallBox.min.x)*m_unitScale;if (diff > 0.0f ){glm::vec3 norm(1.0, 0, 0);float adj = m_boundaryStiffness * diff - m_boundaryDampening * glm::dot ( norm,p->velocity_eval ) ;accel.x += adj * norm.x;accel.y += adj * norm.y;accel.z += adj * norm.z;}diff = 1 * m_unitScale - (m_sphWallBox.max.x - p->pos.x)*m_unitScale;if (diff > 0.0f){glm::vec3 norm(-1.0, 0, 0);float adj = m_boundaryStiffness * diff - m_boundaryDampening * glm::dot ( norm,p->velocity_eval );accel.x += adj * norm.x;accel.y += adj * norm.y;accel.z += adj * norm.z;}//Y方向边界diff = 1 * m_unitScale - ( p->pos.y - m_sphWallBox.min.y )*m_unitScale;if (diff > 0.0f){glm::vec3 norm(0, 1.0, 0);float adj = m_boundaryStiffness * diff - m_boundaryDampening * glm::dot ( norm,p->velocity_eval );accel.x += adj * norm.x;accel.y += adj * norm.y;accel.z += adj * norm.z;}diff = 1 * m_unitScale - ( m_sphWallBox.max.y - p->pos.y )*m_unitScale;if (diff > 0.0f){glm::vec3 norm(0, -1.0, 0);float adj = m_boundaryStiffness * diff - m_boundaryDampening * glm::dot ( norm,p->velocity_eval );accel.x += adj * norm.x;accel.y += adj * norm.y;accel.z += adj * norm.z;}// 重力作用accel += m_gravityDir;// 位置计算----------------------------glm::vec3 vnext = p->velocity + accel*deltaTime;            // v(t+1/2) = v(t-1/2) + a(t) dtp->velocity_eval = (p->velocity + vnext)*0.5f;                // v(t+1) = [v(t-1/2) + v(t+1/2)] * 0.5        used to compute forces laterp->velocity = vnext;p->pos += vnext*deltaTime/m_unitScale;        // p(t+1) = p(t) + v(t+1/2) dt//弹入位置数据posData[3*i]=p->pos.x;posData[3*i+1]=p->pos.y;posData[3*i+2]=p->pos.z;}}

我们每隔0.03秒计算位置，防止速度太快我们还设置了速度限制。然后同样遍历粒子缓存中的每个粒子，在获取每个粒子时对我们的流体容器做简单的碰撞计算，然后考虑重力作用。最后根据加速度和速度计算最后的粒子位置，最后将该点的位置储存在位置缓存中。

然后就是我们的帧动画函数：

    void FluidSystem::tick(){m_gridContainer.insertParticles(&m_pointBuffer);//每帧刷新粒子位置_computerPressure();_computerForce();_advance();}

每次线刷新粒子位置，然后对其计算压力等属性，然后计算加速度，最后计算最终位置，然后重复该循环便可以获得我们的流体动画了。

我们通过opengl来绘制该效果，每次只要VBO里保存我们更新的粒子位置，便能直接绘出效果，效果图如下：

到这里我们的SPH的实现也就结束了，可是光看到这些粒子，怎么把它变成流体呢，我们下一节会探讨一下如何给其加表面。

SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（4）-算法实现2相关推荐

SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（4）-算法实现1
流体模拟(二) SPH算法实现1: 由于我们计算每个粒子的状态时,都需要获得在它光滑核半径内(邻域内)的所有粒子信息.我们如果遍历每个粒子,计算欧式距离的话,那开销就过于庞大了.因此我们可以将我们的空 ...
SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（3）-光滑核函数
流体模拟(二) 光滑核函数: sph中涉及的光滑核可以理解为:在一定的光滑核半径内,所受的力受距离权重的影响,距离越近所受影响越大.其表现形式如图所示. 这里我们便可以将流体看成一个个粒子的集合,每一 ...
SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（2）-粒子受力分析
流体模拟(二) SPH算法的粒子受力分析: SPH算法的基本设想,就是将连续的流体想象成一个个相互作用的微粒,这些粒子相互影响,共同形成了复杂的流体运动.其实现的原理则是我们在初始空间里创建多个粒子, ...
SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（1）-数学原理
流体模拟(二) SPH算法的数学原理: 标量场和矢量场如果空间区域内任意一点P,都有一个确定的数f(P),则称这个空间区域内确定了一个标量场,如果空间区域内任意一点P,都有一个确定的向量vF(P), ...
SPH（光滑粒子流体动力学）流体模拟实现六：Position Based Fluid（PBF）
SPH(光滑粒子流体动力学)流体模拟实现六:Position Based Fluid(PBF) PBF方法和前篇提到的PCISPH方法类似,都属于迭代矫正法.PCISPH是通过迭代预测压力,通过压力变 ...
SPH（光滑粒子流体动力学）流体模拟实现七：屏幕空间流体渲染（SSF）
SPH(光滑粒子流体动力学)流体模拟实现七:屏幕空间流体渲染(SSF) 之前都是用的Marching Cube重建流体表面的方法.近来为了做对比实验,尝试了屏幕空间流体渲染的方法.发现屏幕空间的方法不 ...
SPH（光滑粒子流体动力学）流体模拟实现：算法总览
流体模拟(一) 流体模拟算法总体流程: 流体现象广泛存在于自然界.日常生活以及工业生产中,对流体的模拟即流体动画, 一直是基于物理的动画以及计算机图形学的重要研究内容.目前, 基于物理模拟的流体动画 ...
SPH（光滑粒子流体动力学）流体模拟实现三：Marching Cube算法（1）
流体模拟(三) Marching Cube算法(1) 我们在实现流体表面重建时,需要事先在空间中划分网格,我们的流体系统正好已经完成了此项工作.其次利用Marching Cube算法计算出构成表面的 ...
SPH（光滑粒子流体动力学）流体模拟实现三：Marching Cube算法（2）
流体模拟(三) Marching Cube算法(2) 我们在之前的流体系统类里新加入一些函数和成员,用来引用我们的MC类,便可以获得生成有表面的流体模型了,效果如图: 添加简单的天空盒以及Phong氏 ...

SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（4）-算法实现2

流体模拟（二）

SPH算法实现2：

流体系统类

SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（4）-算法实现2相关推荐

最新文章

热门文章