SPH(光滑粒子流体动力学)流体模拟实现五:PCISPH
SPH(光滑粒子流体动力学)流体模拟实现五:PCISPH
我们知道真实的液体是不可压缩的,但我们在计算机中离散的计算流体运动,在一定的时间步长内,用标准的SPH方法求解,在在粒子聚集处容易发生挤压,造成压缩。有两种常用的方法模拟不可压缩性:1.在WCSPH(弱可压缩SPH)中,利用刚性状态方程(EOS)建模压力。2.通过求解压力泊松方程实现不可压缩性。但这两种方法都有很昂贵的计算费用。
文章“Predictive-Corrective Incompressible SPH”中,提出了一种预测矫正的方法,来使粒子达到不可压缩性,其性能上相比传统两种方法,更加的高效。
PCISPH模型
SPH概述
在拉格朗日粒子描述下, 控制流体运动的偏微分方程 Navier-Stokes 方程可表示为:
SPH方法的核心思想是以离散化粒子的形式来表征连续的场, 并对场量使用积分近似的方式进行计算。位置在xi粒子i的场量:
第i个粒子的密度计算公式为:
压力场直接从Navier-Stokes方程式推导而得:
PCISPH算法
在PCISPH方法中,速度和位置会及时更新,并估计新的粒子密度。然后,对于每个粒子,计算出参考密度的预测变化,并用于更新压力值,压力值又进入压力的重新计算。此过程一直迭代到收敛,即直到所有粒子密度波动均小于用户定义的阈值η(例如1%)。在完成矫正后,我们再更新速度和位置。详细算法流程图如下:
该算法总体思路如下:
1.计算每个粒子的邻居信息,并记录在邻接表内。
2.计算出了压力之外的所有其他力(黏力,重力)。
3.执行矫正循环:执行1).,2).,3).。
1).预测所有粒子新的的速度和位置。
2).预测所有粒子新的密度,以及计算新密度和旧密度之间的差值。
3).预测所有粒子的压力。
4.完成矫正后,利用矫正的压力计算粒子速度和位置。
压强导数
我们的算法是根据预测的密度计算新的压力值,因此我们需要找到它们之间变化的关系。我们的目的是找到一个压强p,该压强以这样一种方式改变粒子的位置,即预测的密度与参考密度相对应。
对于给定的内核平滑长度h,使用SPH密度求和公式计算时间点t + 1处的密度:
其中,。假定非常小,我们可以用一阶泰勒展开近似:
带入得:
我们令,因此密度增量公式为:
我们假设邻居具有相等的压强且密度为静止密度(根据不可压缩性条件),则结果是:
PCISPH算法在每次迭代校正时只矫正d流体粒子的压力,通过迭代计算出的压力让粒子之间不至于靠的过近(粒子之间距离太近可以理解为流体可压缩)。在只考虑压力的情况下,根据时间积分可以计算出粒子在该压力作用下的位移:
粒子i在获得粒子j的压力同时,也会对领居粒子j施加反作用力,因此:
同样,由于粒子i导致粒子j产生的位移为:
将位移增量带入密度增量公式:
因此为:
其中为:
公式表示了我们迭代过程中,需要不断变化的密度值,从而满足粒子不可压缩性。而改变的密度需要改变的压强。我们每次预测的粒子密度为:。预测密度和静止密度之间的误差为:。而我们的目标是让粒子密度为,这需要的密度改变量为。因此,我们需要施加的压强值为:
这个压强计算公式在邻居粒子数目不足的时候会导致计算错误,解决办法是进行一次预计算,即在流体粒子周围充满邻居粒子的情况下计算。我们可以直接在初始化流体时计算一次系数即可:
因此,我们的压强改变量为:
由于只要不满足不可压缩条件,我们就重复进行预测校正步骤,因此,我们需要在迭代中不断矫正压强:
算法实现
我们之前提到,在计算系数时需要在流体初始化时计算。因此我们使用函数_computeGradWValues()和_computeDensityErrorFactor()计算该系数:
void FluidSystem::_init(unsigned int maxPointCounts, const fBox3 &wallBox, const fBox3 &initFluidBox, const glm::vec3 &gravity){ m_pointBuffer.reset(maxPointCounts);m_sphWallBox=wallBox;m_gravityDir=gravity;//根据粒子间距计算粒子质量m_pointMass=m_restDensity*pow(m_pointDistance,3.0);//设定流体块_addFluidVolume(initFluidBox, m_pointDistance/m_unitScale);//MarchingCube算法属性m_mcMesh=rxMCMesh();m_gridContainer.init(wallBox, m_unitScale, m_smoothRadius*2.0, 1.0,m_rexSize,m_gridScale);//设置网格尺寸(2r)m_field=new float[(m_rexSize[0]+1)*(m_rexSize[1]+1)*(m_rexSize[2]+1)]();posData=std::vector<float>(3*m_pointBuffer.size(),0);m_hPosW.resize(3*m_pointBuffer.size(),0.0);m_gridContainer.insertParticles(&m_pointBuffer);//每帧刷新粒子位置m_neighborTable.reset(m_pointBuffer.size());//重置邻接表//计算W的梯度_computeGradWValues();//计算因子_computeDensityErrorFactor();}
其中_computeGradWValues()代码如下:
void FluidSystem::_computeGradWValues(){float h2=m_smoothRadius*m_smoothRadius;//h^2const int numP=m_pointBuffer.size();for (int i=0; i<numP; i++) {Point *pi=m_pointBuffer.get(i);pi->sum_grad_w=glm::vec3(0.0f);pi->sum_grad_w_dot=0.0f;}for (int i=0; i<numP; i++) {Point *pi=m_pointBuffer.get(i);pi->kernel_self=_computeNeighbor(i);m_neighborTable.point_commit();int neighborCounts=m_neighborTable.getNeighborCounts(i);//预测密度计算for (int j=0; j<neighborCounts; j++) {unsigned int neighborIndex;float r;m_neighborTable.getNeighborInfo(i, j, neighborIndex, r);Point* pj=m_pointBuffer.get(neighborIndex);//需要用预测的位置去重新计算距离和核值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;r=pow(r2,0.5f);float h=m_smoothRadius;glm::vec3 gradVec=(pi->pos-pj->pos)*m_kernelSpiky/r*(h-r)*(h-r);pi->sum_grad_w+=gradVec;pj->sum_grad_w-=gradVec;pi->sum_grad_w_dot+=glm::dot(gradVec,gradVec);pj->sum_grad_w_dot+=glm::dot(-1.0f*gradVec,-1.0f*gradVec);}}}}
我们在其中计算所有粒子的和。
之后我们在函数_computeDensityErrorFactor()中计算系数,代码如下:
void FluidSystem::_computeDensityErrorFactor(){float h2=m_smoothRadius*m_smoothRadius;int maxNeighborIndex = 0;int maxNeighs = 0;const int numParticles = m_pointBuffer.size();for (int i=0; i<numParticles; i++) {Point *pi=m_pointBuffer.get(i);int neighborCounts=m_neighborTable.getNeighborCounts(i);int numNeighbors=0;//预测密度计算for (int j=0; j<neighborCounts; j++) {unsigned int neighborIndex;float r;m_neighborTable.getNeighborInfo(i, j, neighborIndex, r);Point* pj=m_pointBuffer.get(neighborIndex);//需要用预测的位置去重新计算距离和核值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){numNeighbors++;}}//获取邻居最多的粒子,和邻居个数if (numNeighbors>maxNeighs) {maxNeighs=numNeighbors;maxNeighborIndex=i;}}//获取邻居最多的粒子Point* pmax=m_pointBuffer.get(maxNeighborIndex);//计算新的压力根据密度误差float restVol=m_pointMass/m_restDensity;float preFactor=2.0*restVol*restVol*m_deltaTime*m_deltaTime;float gradWTerm=glm::dot(pmax->sum_grad_w*(-1.0f),pmax->sum_grad_w)-pmax->sum_grad_w_dot;float divisor=preFactor*gradWTerm;m_densityErrorFactor=-1.0/divisor;const float factor = m_deltaTime / 0.0001f;float densityErrorFactorParameter = 0.05 * factor * factor;m_densityErrorFactor*=1.0*densityErrorFactorParameter;}
值得注意的是,我们还在系数之前乘上densityErrorFactorParameter影响系数,该系数受影响。如果不乘上该系数就会出错,论文里并没有明确提出。
计算完成预计算的系数后,我们的tick()函数如下:
void FluidSystem::tick(bool suf){//求解领居结构m_gridContainer.insertParticles(&m_pointBuffer);//每帧刷新粒子位置m_neighborTable.reset(m_pointBuffer.size());//重置邻接表//计算外力_computerExternalForces();//预测矫正_predictionCorrection();//更新速度和位置_updatePosAndVel();//用于构建表面glm::vec3 tem=m_sphWallBox.min-glm::vec3(1.0);if(suf){CalImplicitFieldDevice(m_rexSize, tem, glm::vec3((1.0/m_gridScale)/m_gridContainer.getDelta()), m_field);clearSuf();//清空表面数据m_mcMesh.CreateMeshV(m_field, tem, (1.0/m_gridScale)/m_gridContainer.getDelta(), m_rexSize, m_thre, m_vrts, m_nrms, m_face);}}
我们按照算法流程,首先求解领居结构,这在之前的章节里介绍过,这里就不多提了。
然后计算除了压力的所有外力,该函数_computerExternalForces()代码如下:
void FluidSystem::_computerExternalForces(){float h2=m_smoothRadius*m_smoothRadius;//h^2for(int i=0;i<m_pointBuffer.size();i++){Point* pi=m_pointBuffer.get(i);pi->forces=glm::vec3(0.0);//邻居粒子装载pi->kernel_self=_computeNeighbor(i);m_neighborTable.point_commit();//外力计算int neighborCounts=m_neighborTable.getNeighborCounts(i);const float restVolume=m_pointMass/m_restDensity;for(unsigned int j=0;j<neighborCounts;j++){unsigned int neighborIndex;float r;m_neighborTable.getNeighborInfo(i, j, neighborIndex, r);Point* pj=m_pointBuffer.get(neighborIndex);float h_r=m_smoothRadius-r;//F_Viscosity//m_kernelViscosity = 45.0f/(3.141592f * h^6);float vterm=m_kernelViscosity*m_viscosity*h_r*restVolume*restVolume;pi->forces-=(pi->velocity-pj->velocity)*vterm;}//F_gravitypi->forces+=m_gravityDir*m_pointMass;//F_boundarypi->forces+=_boundaryForce(pi)*m_pointMass;//初始化矫正因子pi->correction_pressure=0.0f;pi->correction_pressure_froce=glm::vec3(0.0);}}
计算完成外力后,进入矫正预测环节,该函数_predictionCorrection()代码为:
void FluidSystem::_predictionCorrection(){_density_error_too_large=true;int iteration=0;while ((iteration<minLoops)||((_density_error_too_large)&&(iteration<maxLoops))) {for(int i=0;i<m_pointBuffer.size();i++){Point* p=m_pointBuffer.get(i);_predictPositionAndVelocity(p);}//循环终止条件(所有粒子的最大预测密度)_max_predicted_density=0.0;for(int i=0;i<m_pointBuffer.size();i++){_computePredictedDensityAndPressure(i);}//循环终止条件float densityErrorInPercent=max(0.1f*_max_predicted_density-100.0f,0.0f);float maxDensityErrorAllowed=1.0f;//如果密度误差小于终止条件,则终止循环(小于密度误差阈值)if(densityErrorInPercent<maxDensityErrorAllowed)_density_error_too_large=false;for (int i=0; i<m_pointBuffer.size(); i++) {_computeCorrectivePressureForce(i);}iteration++;}
}
其中我们设置最小循环次数为3,最大循环次数为50。在矫正过程中,我们继续细分为3个步骤:
1).预测所有粒子的速度和位置,该函数_predictPositionAndVelocity(p)代码为:
void FluidSystem::_predictPositionAndVelocity(Point* p){// v' = v + delta_t * a// a = F / m// v' = v + delta_t * F * V / m// v' = v + delta_t * F * 1/density//计算预测速度和位置p->predicted_velocity=p->velocity+(p->forces+p->correction_pressure_froce)*m_deltaTime/m_pointMass;p->predicted_position=p->pos+p->predicted_velocity*m_deltaTime;//碰撞处理_collisionHandling(p->predicted_position,p->predicted_velocity);//初始化预测密度p->predicted_density=0.0;
}
值得注意的是,因为是预测的速度,所以需要好用矫正后的压力去计算。
2).计算预测的密度和压力,该函数_computePredictedDensityAndPressure(i)的代码为:
void FluidSystem::_computePredictedDensityAndPressure(int i){float h2=m_smoothRadius*m_smoothRadius;Point* pi=m_pointBuffer.get(i);int neighborCounts=m_neighborTable.getNeighborCounts(i);pi->predicted_density=pi->kernel_self*m_pointMass;//预测密度计算for (int j=0; j<neighborCounts; j++) {unsigned int neighborIndex;float r;m_neighborTable.getNeighborInfo(i, j, neighborIndex, r);Point* pj=m_pointBuffer.get(neighborIndex);//需要用预测的位置去重新计算距离和核值glm::vec3 pi_pj=(pi->predicted_position-pj->predicted_position)*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;pi->predicted_density+=m_kernelPoly6*pow(h2_r2, 3.0)*m_pointMass;}}auto sss=pi->predicted_density;//计算密度误差,仅考虑正向误差pi->density_error=max(pi->predicted_density-m_restDensity,0.0f);//更新压力,只矫正正压力,densityErrorFactor被预先计算并用作常数pi->correction_pressure+=max(pi->density_error*m_densityErrorFactor,0.0f);_max_predicted_density=max(_max_predicted_density,pi->predicted_density);}
值得注意的是,我们这里都只考虑正向的误差,即压缩的密度误差和压强误差,拉伸的我们不做考虑。
3).计算矫正后的压力。该函数_computeCorrectivePressureForce(i)代码为:
void FluidSystem::_computeCorrectivePressureForce(int i){float h2=m_smoothRadius*m_smoothRadius;Point* pi=m_pointBuffer.get(i);int neighborCounts=m_neighborTable.getNeighborCounts(i);pi->correction_pressure_froce=glm::vec3(0.0f);float densSq=m_restDensity*m_restDensity;float pi_pres=pi->correction_pressure;//更新压力for (int j=0; j<neighborCounts; j++) {unsigned int neighborIndex;float r;m_neighborTable.getNeighborInfo(i, j, neighborIndex, r);Point* pj=m_pointBuffer.get(neighborIndex);//需要用预测的位置去重新计算距离和核值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){r=pow(r2,0.5);float h=m_smoothRadius;float pj_pres=pj->correction_pressure;glm::vec3 gradientKer=pi_pj*m_kernelSpiky/r*(h-r)*(h-r);float grad=pi_pres/densSq+pj_pres/(m_restDensity*m_restDensity);pi->correction_pressure_froce-=m_pointMass*m_pointMass*grad*gradientKer;}}
}
这里,我们的预测矫正迭代代码就结束了,我们推出迭代后,利用矫正好的压力求解新的速度和位置,该函数_updatePosAndVel()代码如下:
void FluidSystem::_updatePosAndVel(){for(int i=0;i<m_pointBuffer.size();i++){Point* pi=m_pointBuffer.get(i);pi->forces+=pi->correction_pressure_froce;const float invMass=1.0/m_pointMass;// //Symplectic Euler Integration
// pi->velocity+=pi->forces*invMass*m_deltaTime;
// pi->pos+=pi->velocity*m_deltaTime/m_unitScale;//Leapfrog integrationglm::vec3 vnext = pi->velocity + pi->forces*invMass*m_deltaTime; // v(t+1/2) = v(t-1/2) + a(t) dtpi->velocity_eval = (pi->velocity + vnext)*0.5f; //v(t+1) = [v(t-1/2) + v(t+1/2)] * 0.5pi->velocity = vnext;pi->pos += vnext*m_deltaTime/m_unitScale; // p(t+1) = p(t) + v(t+1/2) dt//弹入位置数据posData[3*i]=pi->pos.x;posData[3*i+1]=pi->pos.y;posData[3*i+2]=pi->pos.z;}
}
完成结果图:
对比普通SPH,我们发现在挤压处,PCISPH有明显的矫正。普通SPH效果如下:
SPH(光滑粒子流体动力学)流体模拟实现五:PCISPH相关推荐
- 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(光滑粒子流体动力学)流体模拟实现二:SPH算法(4)-算法实现2
流体模拟(二) SPH算法实现2: 在前面一节我们完成了粒子缓存类,网格类和邻接表类.我们现在可以正式的整合在我们的流体系统类中了. 流体系统类 class FluidSystem{public:Fl ...
- SPH(光滑粒子流体动力学)流体模拟实现二:SPH算法(4)-算法实现1
流体模拟(二) SPH算法实现1: 由于我们计算每个粒子的状态时,都需要获得在它光滑核半径内(邻域内)的所有粒子信息.我们如果遍历每个粒子,计算欧式距离的话,那开销就过于庞大了.因此我们可以将我们的空 ...
- SPH(光滑粒子流体动力学)流体模拟实现二:SPH算法(3)-光滑核函数
流体模拟(二) 光滑核函数: sph中涉及的光滑核可以理解为:在一定的光滑核半径内,所受的力受距离权重的影响,距离越近所受影响越大.其表现形式如图所示. 这里我们便可以将流体看成一个个粒子的集合,每一 ...
- SPH(光滑粒子流体动力学)流体模拟实现四:各向异性(Anisotropic)表面光滑(1)
流体模拟四: 各向异性(Anisotropic)表面光滑(1) 表面的表示与定义 我们是以隐式表面来表示流体表面.定义一个足够大的参考域,它完全包含粒子表示的流体所在的空间D,即.流体表面隐含的定义为 ...
- SPH(光滑粒子流体动力学)流体模拟实现四:各向异性(Anisotropic)表面光滑(2)
流体模拟四: 各向异性(Anisotropic)表面光滑(2) 理论整理清楚了,我们看一下代码实现,我们首先在我们之前的FluidSystem类中添加几个成员函数和对象: class FluidSys ...
- SPH(光滑粒子流体动力学)流体模拟实现三:Marching Cube算法(2)
流体模拟(三) Marching Cube算法(2) 我们在之前的流体系统类里新加入一些函数和成员,用来引用我们的MC类,便可以获得生成有表面的流体模型了,效果如图: 添加简单的天空盒以及Phong氏 ...
最新文章
- Automation Test in Maya Plugin Development
- Spring-Data-MongoDB的Index注解的expireAfterSeconds参数不起作用?解决方案居然是这样的!...
- Spring Security OAuth2源码解析(一)
- linux调用odbc接口乱码,linux中pypyodbc读取GB编码mdb中文乱码解决办法
- 根据后序和中序求二叉树的先序
- chengxuyuan
- 如何让不懂信息化的甲方客户看懂需求文档,并确认签字?
- 未来函数在线检测_嵌入式实时操作系统任务栈溢出检测原理
- Nginx (1)---安装配置
- 西华大学c语言考试题,西华大学C语言程序设计复习题
- c 和java用cfb_一文彻底搞懂Java中的环境变量
- 关于VC6.0一些常见问题和解决方案
- 数据仓库与数据挖掘实践期末复习总结
- SQL Server 添加字段 修改字段 删除字段 语句
- windows下常见php集成环境安装包介绍
- 怎样使用Scanner(扫描仪),超级详细,不容错过!!!
- 在衣食住行上训练专注力
- Can‘t Update No tracked branch configured for branch
- ZInt支持中文例子
- matlab 保存.fig文件后无法保存的问题