简介

参考doyub Kim那本《Fluid Engine Development》写一个最简单的弱可压SPH。
目前有BUG, 粒子太散了

效果展示

CSDN有图片大小限制，大概就这样

代码

https://github.com/chunleili/WCSPHTaichiHW

正文结束

以下是旧的代码文档，新文档请看github链接里面的word

1程序概述

二维SPH小程序

编程语言：taichi

2程序主要模块说明

主要步骤包含：
1.定义核函数
2.定义邻域搜索函数
3.计算密度
4.计算压力
5.利用核函数计算三种力并得到加速度
6.时间推进
7.处理边界碰撞
8.初始化（只需一次）

附录：
A.邻域搜索具体实现
B.时间步的选定
C.给定密度的选定

2.1核函数

函数/模块名称：kernel()
输入：double r，为从当前粒子中心指向相邻粒子j中心的距离
作用：利用如下核函数进行插值离散

输出：double interplolatedValue, 即插值之后的值。

2.2 邻域搜索

函数/模块名称：neighborSearch()
输入：double currentParticleID 当前粒子的ID
作用：根据当前粒子找到其周围粒子。
输出：周围粒子的ID（多个）。
注：因粒子位置存储在全局变量position中，故只要给出ID即可找到对应位置，即一个哈希表。具体邻域搜索算法需要采用gridBased搜索建立背景网格。

2.3 计算密度

函数/模块名称：computeDensity()
输入： mass, 粒子质量(全局常量)
作用：计算密度场，计算公式如下：

输出：密度场Density（全局场）。

2.4 计算压力

函数/模块名称：computePressure()
输入：密度场density；给定密度density0；常量参数EOSgamma；常量参数EOSkappa
作用：根据理想气体状态方程（EOS）计算压力场。计算公式如下：

输出：压力场（全局场）。

2.5 计算力和加速度

分为三个子模块：压力梯度力，粘性力和加速度

2.5.1 计算压力梯度力

函数/模块名称：computePressureGradientForce ()
输入：pressure，压力场（需辅助函数额外计算，全局场）
作用：计算压力梯度力，公式为：

输出：压力梯度力场pressureGradientForce（全局场）

2.5.2 计算粘性力

函数/模块名称：computeViscosityForce ()
输入：velocity，速度场（全局场）；viscosity, 粘度（全局常量）
作用：计算粘性力，公式为：

输出：粘性力场viscosityForce（全局场）

2.5.3 叠加力并计算加速度

函数/模块名称：computeAcceleration()
输入：已计算出的两个力场：pressureGradientForce; viscosityForce
重力加速度gravityAcceleration
旧加速度场acceleration
作用：计算加速度，公式为：

输出：加速度（全局场）

2.6 时间推进

函数/模块名称：advanceTime ()
输入：旧位置场position；旧速度场velocity；计算得到的加速度acceleration
作用：向前推进时间步。利用欧拉显式时间积分更新粒子位置和速度，公式为：

输出：新速度场velocity和加速度场position（全局场）

2.7 处理边界碰撞

函数/模块名称：boundaryCollision()
输入：粒子当前位置position，边界位置，粒子速度velcotiy，弹性恢复系数restitutionCoeff; 摩擦系数frictionCoeff
作用：判断粒子是否穿过壁面边界，如果穿过，则将粒子位置挪到最近的边界点处，速度修正为：

输出：粒子新速度场velocity和新位置场position

2.8 初始化

函数/模块名称：initialization()
输入：无
作用：按照初始条件初始化如下场：
位置场position：按照粒子初始排列位置给值，为图示方块。
速度场velocity：全0
密度场density：无需给定，由computeDensity根据m计算
压力场pressure：无需给定，按照EOS自动计算。
其余场如压力梯度力和加速度等均自动初始化为0。

设置水方块大小为0.20.2 m
计算域大小为1.01.0 m

输出：粒子新速度场velocity和新位置场position

3 变量说明

3.1 全局常量

粒子数目 particleNumber
粒子给定密度 density0
核半径 kernelRadius
粒子质量 mass =pikernelRadius**2density0
整个计算区域大小（方形） wholeBoundarySize
水区域大小（方形） waterBoundarySize
弹性恢复系数 restitutionCoeff
摩擦系数 frictionCoeff
EOS参数(指数) EOSgamma
EOS参数(乘数） EOSkappa
粘度 viscosity
时间步 dt

3.2全局场

场大小均为 particleNumber
位置场 position 2D
速度场 velocity 2D
密度场 density 1D
压力场 pressure 1D
加速度场 acceleration 2D
压力梯度力场pressureGradientForce 2D
粘性力场 viscosityForce 2D
重力加速度 gravityAcceleration 2D

注：重力加速度场定义为场而非常量是为了方便拓展，其他任何附加的体积力都可以加到重力加速度场上。

附录A 邻域搜索

算法为基于网格的邻域搜索

注：这里的网格只是一种邻域搜索时使用辅助的背景网格，不是真正的网格。不要和网格法或者混合法里的网格搞混。

原理大概就是：

我们想要求解的问题：知道粒子的编号（particleID），想求和它在相邻网格的粒子的编号（neighbors）。
因为只要知道相邻网格内的粒子，就只需要搜索相邻网格而不用搜索所有粒子了。

如图：

图中x代表当前粒子，虚线圆圈代表其邻域（也就是核半径内区域），黑点代表其他粒子。
可见，只有相邻四个网格内的粒子有可能是在邻域内的，这样就不必搜索其他粒子，大大节约了计算量。

具体实现1

设置两个数组（两个哈希表）：

数组1（particleID2CellID[n]）：particleID作为下标（哈希表的键），cellID作为存储的值（哈希表的值）。
数组2（cellID2particleID[m]）：cellID作为下标（哈希表的键），particleID作为存储的值（哈希表的值，有多个）。这是个二维数组，其第二维长度是不定的，其形状示意图如下：
n代表粒子总数，m代表网格总数。
其第一维代表的是网格，按网格编号排列，第二位是该网格内的粒子。

其中cellID不用真正建立网格，只要根据Cartisian 坐标得到即可
假设计算域大小为boundX, boundY
于是每行排列的网格数目为：parNumPerRow=boundX/cellSize
其中cellSize可以设置为2倍核半径：cellSize=2*kernelRadius
根据当前粒子位置可以计算cellID：
cellID=floor(position[particleID].x/cellSize)+floor(position[particleID].y/cellSize)*parNumPerRow
floor代表舍弃小数点后，只保留整数
即为一行一行排列，每行满了之后再往上排列。

每次时间步开始的时候，先更新这两个数组。
更新方法：particleID2CellID即根据上面计算的公式得到cellID，与此同时，把相应的粒子编号放入cellD2ParticleID。

于是回到原问题：已知粒子ID求相邻粒子。

根据particleID找到所在网格cellID
根据网格cellID找到其相邻的四个网格，cellID, cellID-1, cellID-parNumPerRow, cellID-parNumPerRow-1（这里是为了和图示相符，具体要根据位置小数点后面的数字再做判断, 附录A.1详细分类讨论）
在这四个网格内计算和当前粒子距离，假如在核半径内，就认为它在邻域内，可以存入一个邻域链表（neiborList）。

注：维持邻域链表是为了以后访问邻域信息方便。

邻域链表也是这样的形状，只不过现在第一维度代表的是粒子，之前第一维代表的是网格。

这个算法的计算量都在于维持这两个哈希表了，所以每个时间步的算法复杂度为2n

具体实现2

可以不用维持两个哈希表，只需要一个即可。但是每个时间步需要进行排序。

由于需要排序，又要能随机访问，可以用类似于c++的map来实现。

该哈希表的键为粒子ID。值为cellID

每个时间步开始的时候，按照cellID排序
排序之后的表:

key (particleID)	value(cellID)
126	0
255	0
771	0
140	0
331	1
233	1
439	1
591	2
72	2
…	…

当已知粒子编号需要检索邻域粒子的时候：

根据粒子编号，直接索引到网格编号
根据网格编号，在该哈希表内向上向下搜索，直到cellID!=当前粒子的cellID。同样地，对相邻的四个网格都这么处理。
在这四个网格内计算和当前粒子距离，假如在核半径内，就认为它在邻域内，存入邻域链表（neiborList）。

和上面的算法相比，每个时间步需要进行一次搜索，同时还需要进行上下双向搜索。我们就姑且假设双向搜索只需要搜索10次吧（因为一个网格超过10个粒子就太密了）。算法复杂度为nlogn+10m。而恰好我们假设的是一个网格内有10个粒子，所以其实这里可以大致认为10m=n。所以可以认为算法复杂度为nlogn+n

附录 A.1 相邻的四个网格

按照当前粒子在所在网格的位置，可分为四种情况，如图所示

橘色代表当前粒子落在的位置，标号123表示相邻的3个网格
那么就按照粒子位置的小数部分来分类

    fracPartX=position[thisID].x-int(position[thisID].x)fracPartY=position[thisID].y-int(position[thisID].y)if fracPartX<=0.5 and fracPartY<=0.5:neiborCell1=thisCellID-1neiborCell2=thisCellID-parNumPerRow-1neiborCell3=thisCellID-parNumPerRowelif fracPartX>0.5 and fracPartY<=0.5:neiborCell1=thisCellID+1neiborCell2=thisCellID+parNumPerRowneiborCell3=thisCellID+parNumPerRow+1elif fracPartX>0.5 and fracPartY>0.5:neiborCell1=thisCellID-parNumPerRowneiborCell2=thisCellID-parNumPerRow+1neiborCell3=thisCellID+1elif fracPartX<=0.5 and fracPartY>0.5:neiborCell1=thisCellID-parNumPerRow-1neiborCell2=thisCellID-parNumPerRowneiborCell3=thisCellID-1

利用taichi写一个最简单的SPH(光滑粒子动力学）相关推荐

利用python写一个简单的双色球彩票系统
利用python写一个简单的双色球彩票系统 1.设置每次买的号码一样写一个双色球彩票系统,系统可以随机产生一组数据,一组彩票数据有六位数,这六位数的的取值范围是0和1. 一张彩票是两块钱,用户可以选 ...
学了C语言，如何利用CURL写一个下载程序？—用nmake编译CURL并安装
在这一系列的前一篇文章学了C语言,如何为下载狂人写一个磁盘剩余容量监控程序?中,我们为下载狂人写了一个程序来监视磁盘的剩余容量,防止下载的东西撑爆了硬盘.可是,这两天,他又抱怨他的下载程序不好用,让我 ...
python抽奖游戏_利用Python写一个抽奖程序，解密游戏内抽奖的秘密
原标题:利用Python写一个抽奖程序,解密游戏内抽奖的秘密前言本文的文字及图片来源于网络,仅供学习.交流使用,不具有任何商业用途,版权归原作者所有,如有问题请及时联系我们以作处理. 作者: 极客 ...
利用Flutter写一个跨平台的果核APP（4）——数据存储
前言目前我们已经实现了几个界面,今天这篇文章开始着手进行登录页的制作,主要流程就是获取输入框中的内容,发送给后台进行验证,如果成功将返回信息保存在本地并跳转至首页,如果失败就提示用户重新输入. 在这 ...
利用SpringCloud搭建一个最简单的微服务框架
利用SpringCloud搭建一个最简单的微服务框架 https://blog.csdn.net/caicongyang/article/details/52974406 1.微服务微服务主要包含服 ...
Java入门知识，写一个最简单java程序
本文目录一.Java语言的简介二.写一个最简单Java程序 1.Notepad配置 2.最简单的Java程序 3.代码分析 4.编译与编译常见错误 5.执行java程序一.Java语言的简介 0 ...
利用Cocos2d-x写一个程序读取传奇wzl文件
Cocos2d-x是一个用于游戏开发的开源框架,它提供了用于制作2D游戏的工具和功能.若要利用Cocos2d-x读取传奇wzl文件,需要对wzl文件的格式进行分析,并使用Cocos2d-x提供的读取文 ...
用JAVA写一个最简单的飞翔的小鸟
如果你想写一个最简单的飞翔的小鸟的 Java 程序,可以先了解 Java 的图形绘制功能.Java 提供了一个叫做 Graphics 的图形绘制类,可以用来绘制图形.填充颜色.画线等. 你可以通过创建 ...
利用itchat写一个聊天机器人
利用itchat写一个聊天机器人聊天机器人图灵机器人需要的库 **自动回复私聊消息** **自动回复群聊消息** 结语: 聊天机器人偶然在CSDN上看到大佬用20行教你写一个聊天机器人,觉得甚 ...

利用taichi写一个最简单的SPH(光滑粒子动力学）

简介