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

最新推荐文章于 2022-11-23 20:59:43 发布
最新推荐文章于 2022-11-23 20:59:43 发布
阅读量2.3k 收藏 22
点赞数 7
分类专栏： 计算机图形学 # 流体力学 文章标签： 图形学
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_43940314/article/details/120922210
版权

简介

参考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求相邻粒子。

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

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

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

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

具体实现2

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

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

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

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

key (particleID)value(cellID)
1260
2550
7710
1400
3311
2331
4391
5912
722
……

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

  1. 根据粒子编号，直接索引到网格编号
  2. 根据网格编号，在该哈希表内向上向下搜索，直到cellID!=当前粒子的cellID。同样地，对相邻的四个网格都这么处理。
  3. 在这四个网格内计算和当前粒子距离，假如在核半径内，就认为它在邻域内，存入邻域链表（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-1
        neiborCell2=thisCellID-parNumPerRow-1
        neiborCell3=thisCellID-parNumPerRow
    elif fracPartX>0.5 and fracPartY<=0.5:
        neiborCell1=thisCellID+1
        neiborCell2=thisCellID+parNumPerRow
        neiborCell3=thisCellID+parNumPerRow+1
    elif fracPartX>0.5 and fracPartY>0.5:
        neiborCell1=thisCellID-parNumPerRow
        neiborCell2=thisCellID-parNumPerRow+1
        neiborCell3=thisCellID+1
    elif fracPartX<=0.5 and fracPartY>0.5:
        neiborCell1=thisCellID-parNumPerRow-1
        neiborCell2=thisCellID-parNumPerRow
        neiborCell3=thisCellID-1
beidou111
关注
  • 7
    点赞
  • 踩
  • 22
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
经典流体模拟—SPH源代码
08-01
自己最近用SPH方法做的一个经典流体模拟实例，花了很长时间才圆满完成，已实测，可行！
taichi生成sph
最新发布
apple_64921913的博客
02-28 364
ti.kernel@ti.kernelpass@ti.kernel。
【taichi】关于SPH_Taichi的探索与尝试
beidou111的博客
11-23 757
这是太极实习生张铭睿同学做的SPH太极库。请看是作者的讲解。我们来探索一下这个库。这只是个笔记，写得会比较细碎。
wcsph:使用taichi创建sph模拟
05-25
弱压缩liqiud模拟 该代码主要是从移植的代码 该代码使用编程语言 当前执行器算法 [SESPH] [PCISPH] [IISPH] [DFSPH] 怎么跑 首先配置您的anaconda工作区，然后打开anaconda提示符 其次，您需要安装语言， pip安装taichi 最后，您输入ti dfsph.py ，仅此而已 该项目产生的一些图像 位置盘样本（完全并行），用于边界处理 隐式粘度求解器 使用前提条件cg解决粘度： 表面张力（不带和带）： 各向异性网格重建： 行军魔方 使用各向异性内核建立体积场 渲染图像 --- 其他算法： 哈希网格 cfl时间步 推荐人 M Weiler 2018：适用于SPH流体的物理一致性隐式粘度求解器 Nadir Akinci 2013：SPH流体的通用表面张力和附着力 JIHUN YU 2013：重构基于粒子的流体的表面。
SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（3）-光滑核函数
小龙虾的博客
10-17 4290
流体模拟（二） 光滑核函数： sph中涉及的光滑核可以理解为：在一定的光滑核半径内，所受的力受距离权重的影响，距离越近所受影响越大。其表现形式如图所示。 这里我们便可以将流体看成一个个粒子的集合，每一个粒子都受周围一定范围内的其他粒子影响。该粒子的最终属性都由周围所有粒子的加权值来决定。如下图所示： 我们在流体中随机选择一个粒子,在光滑核半径h范围内有多个粒子，位置分别是，…，则该...
SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（4）-算法实现1
小龙虾的博客
10-21 7735
流体模拟（二） SPH算法实现1： 由于我们计算每个粒子的状态时，都需要获得在它光滑核半径内（邻域内）的所有粒子信息。我们如果遍历每个粒子，计算欧式距离的话，那开销就过于庞大了。因此我们可以将我们的空间划分成多个网格。 如上图所示，我们可以把空间划分成均匀网格，从而使我们在获取每个粒子的邻域粒子时，只需要访问该粒子的周围网格获取邻域粒子。这时候我们可以使用哈希表的方式，给每个空间网格赋予...
zygisk-taichi-v8.1.0.zip
03-14
【zygisk-taichi-v8.1.0.zip】是一个特定版本的模块化系统增强工具，主要用于Android设备。Zygisk是Xposed框架的一种新型实现，由知名的开发者"Taichi"开发。Xposed框架是一个在Android系统底层运行的插件框架，允许...
LBM_Taichi:taichi编程语言实现的基于Lattice Boltzmann方法的流体求解器
03-08
用法为了用数值方法解决流体动力学问题，应给出畴尺寸，流体性质，边界条件和初始条件。 在此代码中，可以通过实例化求解器来指定这些参数： lbm = lbm_solver(nx, ny, niu, bc_type, bc_value)每个参数的含义是： ...
docs.taichi.graphics:Taichi 文档站点的主页
08-05
如果您想更新任何 Taichi 文档，请转到，请不要在此工作！ 此 repo 仅托管文档的镜像版本，以便构建和部署文档网站。 如果您想为任何 Taichi 文档的翻译做出贡献，请转至。 我们的页面包含有关 i18n 工作流程的更...
taichi_7.0.0.release.apk
05-09
太极是一个能够运行 Xposed 模块的框架，模块能通过它改变系统和应用的行为。太极既能以传统的 Root/...简单来说，太极就是一个 类 Xposed框架，它能够加载 Xposed 模块、修改系统和 APP、拦截方法，执行 hook 逻辑等。
光滑粒子流体动力学（SPH）源程序
10-08
Dr. Liu 在《光滑粒子流体动力学--一种无网格方法》一书中提供的光滑粒子流体动力学（SPH）源程序，采用Fortran语言编写，用于计算可压缩流体力学问题。
SPH无网格程序
03-31
SPH无网格程序，可以调试使用的工程文件，C++和Fortran语音混合编程！
用sph方法模拟流动，亲测可用，适合初学者学习
07-20
自己最近用sph方法模拟的流动，亲测可用，适合初学者学习，经典应与人分享，拿去不谢！
sph代码(美国)
04-03
美国人编制的溃坝的sph源程序，非常经典，适合初学者入门使用！
SPH 流体动画源代码
11-24
光滑流体粒子动力学（SPH）方法简单实例代码，包括可运行程序。 其中，图形绘制部分采用 DirectX 引擎。 具体分析过程，见博客：http://blog.csdn.net/ironyoung/article/details/8068929
【CG物理模拟系列】流体模拟--粒子法之SPH（代码讲解）
【计算机图形】Physics-Based Animation-小8君的专栏
10-25 5025
WCSPH,PCISPH,IISPH等研究方法，其本质都是以非压缩性为目标，求解Navier-Stokes方程。 本文以WCSPH为例，讲解下SPH方法代码的实现。 代码讲解 sph_type.h里定义几个vector函数类型 #ifndef _SPHTYPE_H #define _SPHTYPE_H typedef unsigned int uint;
SPH算法简介（四）:Hello,SPH
狼牙
11-06 5531
SPH算法简介（四）:Hello,SPH 2011年04月2日 |本网站遵守CC版权协议 转载请注明出自www.thecodeway.com     上几节，我们推导出一大推复杂无比的公式，似乎有点纸上谈兵，这节来点真的，写一个可以运行的SPH系统，下面就是SPH基本的运算流程 初始化粒子，为每个粒子赋上初始位置根据公式3.7计算每个粒子的密度根据公式3.10计算每个粒子的压强根据公
SPH（光滑粒子流体动力学）流体模拟实现四：各向异性(Anisotropic)表面光滑(2)
小龙虾的博客
12-01 1623
流体模拟四： 各向异性(Anisotropic)表面光滑(2) 理论整理清楚了，我们看一下代码实现,我们首先在我们之前的FluidSystem类中添加几个成员函数和对象： class FluidSystem{ ......//之前的代码 void CalAnisotropicKernel();//各项异性核计算 ...
SPH法理论(1)
RMG的专栏
04-22 3030
SPH法理论与实践 前言 SPH(Smoothed Particle Hydrodynamics)法，中文称其为光滑粒子动力学。最早是为了模拟银河系中天体之间的碰撞与形成等等宇宙物理学的现象而提出的算法[Lucy1977]1。虽然如此，近些年来SPH法在流体模拟热扩散等等的领域内都被广泛应用了起来。而CG(Computer Graphics)相关的研究者[Desbrun1996]2最早在1996年...
用taichi编程语言写一个欧拉视角的流体仿真程序，带注释
05-26
好的，以下是一个简单的欧拉视角的流体仿真程序，使用了Taichi编程语言。代码中包含了注释，希望能够帮助您理解。 ```python import taichi as ti ti.init(arch=ti.gpu) # 初始化Taichi，指定使用GPU # 定义仿真参数 dim = 2 # 维度 res = 128 # 分辨率 dt = 1e-4 # 时间步长 viscosity = 0.01 # 粘性系数 gravity = ti.Vector([0.0, -9.8]) # 重力 # 声明变量 vel = ti.Vector(dim, dt=ti.f32, shape=(res, res)) # 速度场 vel_new = ti.Vector(dim, dt=ti.f32, shape=(res, res)) # 新的速度场 density = ti.var(dt=ti.f32, shape=(res, res)) # 密度场 density_new = ti.var(dt=ti.f32, shape=(res, res)) # 新的密度场 # 定义函数，用于计算速度场的散度 @ti.kernel def divergence(): for i, j in vel: vel_new[i, j] = vel[i, j] + dt * (-gravity + viscosity * ti.Matrix([[0,-1],[1,0]]) @ vel.grad[i, j]) if i > 0 and j > 0 and i < res-1 and j < res-1: density_new[i, j] = density[i, j] + dt * (-vel_new[i, j].dot(ti.Vector([density.grad[i, j][0], density.grad[i, j][1]]))) # 交换速度场和新的速度场，密度场和新的密度场 vel, vel_new = vel_new, vel density, density_new = density_new, density # 主循环 for step in range(1000): divergence() # 计算速度场的散度 if step % 50 == 0: ti.imshow(density.to_numpy(), cmap='gray') # 显示密度场 ``` 上述代码实现的是一个简单的二维流体仿真，使用了欧拉视角的方法。程序中主要包含了以下几个部分： 1. 初始化Taichi，并设置使用GPU加速。 2. 定义了一些仿真参数，包括维度、分辨率、时间步长、粘性系数和重力。 3. 声明了一些用于存储速度场、密度场的变量，以及用于更新速度场、密度场的变量。 4. 定义了一个计算速度场的散度的函数，其中使用了Taichi的自动求导功能来计算速度场的梯度，并使用旋度算子计算速度场的散度。 5. 在主循环中调用计算速度场的散度的函数，并在每50个时间步长之后显示密度场。 需要注意的是，上述代码只是一个非常简单的流体仿真程序，还有很多方面需要进一步完善，如边界条件、数值稳定性等。但是这个程序可以作为入门级别的流体仿真教程，帮助初学者理解欧拉视角的流体仿真原理。
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值