P5.js——基于SPH并融入动画技术的交互
本次实验是基于P5.js和SPH做的一个类似流体的一个交互艺术程序。
背景与设想
流体流动存在两种运动状态：层流和湍流。倘流速很慢，流体会分层流动，互不混合，此乃层流。倘流速增加，越来越快，流体开始出波动性摆动，此情况称之为过渡流。当流速继续增加，达到流线不能清楚分辨，会出现很多漩涡，这便是湍流，又称作乱流、扰流或紊流。
当流体中发生了层与层之间的相对运动时，速度快的层对速度慢的层产生了一个拖动力使它加速，而速度慢的流体层对速度快的就有一个阻止它向前运动的阻力，拖动力和阻力是大小相等方向相反的一对力，分别作用在两个紧挨着但速度不同的流体层上，这就是流体粘性的表现，称为内摩擦力或叫粘滞力。为了维持流体的运动就必须消耗能量来克服由于内摩擦力产生的能量损失，这就是流体运动时会造成能量损失的原因。实际上，粘性是流体阻止发生剪切变形和角变形的一种特性。这是由于内聚力的存在和流体层间的动量交换而造成的。内摩擦力就是这种特性的表现形式。当流体处于静止或各部分之间相对速度为零时，流体的粘性就表现不出来，内摩擦阻力也就等于零。
对于自然界中流体的运动做的一个模拟，可以用鼠标操控粒子的流向。自然与技术的融合可以产生与众不同的美与创意。
简述SPH
SPH (Smoothed Particle Hydrodynamics)是光滑粒子流体动力学方法的缩写，是在近20多年来逐步发展起来的一种无网格方法。该方法的基本思想是将连续的流体用相互作用的粒子来描述，各个粒子上承载各种物理量，包括质量、速度等，通过求解粒子的动力学方程和跟踪每个粒子的运动轨道，求得整个系统的力学行为。
SPH方法
1、对每个粒子，通过它附近周围的所有和它的距离半径小于r的所有粒子来计算出该粒子点的密度r被称为Smooth Length，通俗的说就是每个粒子的最大影响半径。
2、对于每一个粒子，通过理想气体状态方程根据上一步计算出的密度从而算出它的压强。
3、对于每一个粒子，根据它附近周围距离半径小于r的所有粒子的压强（上一步已经计算出）来计算出该粒子由于附近压强差而导致受到的压力差。
4、对于每一个粒子根据它附近周围距离半径小于r的所有粒子的速度差来计算出该粒子由于附近速度差而导致受到的粘滞力。
粒子的受力分析
作用在粒子上的力由三部分组成：
分别是外力、压力、粘滞力。
1.外力一般是重力。（本次实验中附加了一个鼠标拖动的外力）
(PS：这里的力F的量纲发生了变化，正常情况下，F=m*g。但在SPH算法里，流体的质量是由流体单元的密度决定的，所以一般用密度代替质量，后面的分析都是用这个量纲的“作用力”)
2.压力是由流体内部的压力差产生的作用力，试想一下在水管中流动的液体，进水口区域的压力一定会比出水口区域大，所以液体才会源源不断的流动，数值上，它等于压力场的梯度，方向由压力高的区域指向压力低的区域。
3.粘滞力是由粒子之间的速度差引起的，设想在流动的液体内部，快速流动的部分会施加类似于剪切力的作用力到速度慢的部分，这个力的大小跟流体的粘度系数μ以及速度差有关。
最后得到粒子的受力公式 ：
加速度形式：
SPH算法同样涉及到“光滑核”的概念，可以这样理解这个概念，粒子的属性都会“扩散”到周围，并且随着距离的增加影响逐渐变小，这种随着距离而衰减的函数被称为“光滑核”函数，最大影响半径为“光滑核半径”。
函数形式
1.计算使用的光滑核函数称为Poly6函数，具体形式为:
其中KPoly6 是一个固定的系数，根据光滑核的规整属性，通过积分计算出这个系数的具体值，在极坐标中计算积分：
2.压力计算中使用的光滑核函数称为Spiky函数
3.粘度的光滑核函数形式如下：
SPH基本的运算流程
(PS：后续代码是基于上述理论公式进行调整改编，最终达到个人满意的效果)
代码展示
先进行初始化，在setup()函数中调用 initialize()
function initialize() {
//产生粒子
makeParticles();
//初始化粒子数组
initializeArrays();
//初始化参数
initializeParams();
}
产生粒子的函数代码如下（其中有些参数的含义自行查看最上面定义时的注释）
function makeParticles() {
posX = [];
posY = [];
N = 0;
var half = PARTICLE_SPACING / 2.0;
for (var x = half; x < INITIAL_FLUID_WIDTH; x += PARTICLE_SPACING) {
for (var y = half; y < INITIAL_FLUID_HEIGHT; y += PARTICLE_SPACING) {
posX.push(x + random(-0.5, 0.5) * half);
posY.push(y + random(-0.5, 0.5) * half);
N++;
}
}
}
初始化上面产生粒子的数组
function initializeArrays() {
velX = new Array(N);
velY = new Array(N);
accX = new Array(N);
accY = new Array(N);
dens = new Array(N);
press = new Array(N);
neighbor = new Array(N);
for (var i = 0; i < N; i++) {
velX[i] = 0.0;
velY[i] = 0.0;
accX[i] = 0.0;
accY[i] = 0.0;
dens[i] = 0.0;
press[i] = 0.0;
neighbor[i] = [];
}
}
参数初始化函数
function initializeParams() {
KERNEL_POLY6_COEF = 4.0 / (PI * sq(KERNEL_RADIUS));
KERNEL_SPIKY_GRAD_COEF = -30.0 / (PI * pow(KERNEL_RADIUS, 3.0));
KERNEL_VISCO_LAP_COEF = 20.0 / (3.0 * PI * pow(KERNEL_RADIUS, 4.0));
calcParticleMass();
}
draw()中的函数，绘制各个粒子并更新
function draw() {
//console.log('frameCount: ' + frameCount)
render();
update();
}
function render() {
background(0);
var xScale = width / POOL_WIDTH;
var yScale = height / POOL_HEIGHT;
var kRadX = KERNEL_RADIUS * xScale;
var kRadY = KERNEL_RADIUS * yScale;
noStroke();
for (var i = 0; i < N; i++) {
var x = posX[i] * xScale;
var y = height - posY[i] * yScale;
//fill(100, 100, 255, 50);
fill(150,10,i/5,60);
ellipse(x, y, kRadX, kRadY);
}
for (var i = 0; i < N; i++) {
var x = posX[i] * xScale;
var y = height - posY[i] * yScale;
fill(255,255,i/5,90);
ellipse(x, y, 5, 5);
}
}
更新粒子的位置，重新计算各个参数
//更新粒子位置并计算各个属性
function update() {
calcNeighbors();//邻域
calcDensity();//密度
calcExternalForce();//外力
calcGravity();//重力
calcViscosity();//粘滞力
calcPressure();//压力
calcPressureGradient();//梯度压力
updateParticles();//更新粒子
calcCollision();//碰撞
}
其中外力是根据鼠标拖动粒子来计算的
function calcExternalForce() {
if (!mouseIsPressed) return;
var mx = mouseX / width * POOL_WIDTH;
var my = (1.0 - mouseY / height) * POOL_HEIGHT;
var pmx = pmouseX / width * POOL_WIDTH;
var pmy = (1.0 - pmouseY / height) * POOL_HEIGHT;
var mlength = sqrt(sq(mx - pmx) + sq(my - pmy));
if (mlength < 0.0000001) return;
// var dirX = (mx - pmx);
// var dirY = (my - pmy);
var dirX = (mx - pmx) / mlength;
var dirY = (my - pmy) / mlength;
for (var i = 0; i < N; i++) {
var d = sqrt(sq(posX[i] - mx) + sq(posY[i] - my));
if (d > EXTERNAL_FORCE_RADIUS) continue;
var decay = 1.0 - d / EXTERNAL_FORCE_RADIUS;
accX[i] += PARTICLE_MASS * dirX * decay * EXTERNAL_FORCE_SCALE;
accY[i] += PARTICLE_MASS * dirY * decay * EXTERNAL_FORCE_SCALE;
}
}
碰撞相关代码
function calcCollision() {
for (var i = 0; i < N; i++) {
if (posX[i] < 0.0) {
velX[i] *= velX[i] < 0.0 ? -1.0 : 1.0;
velX[i] *= RESTITUTION_SCALE;
velY[i] *= FRICTION_SCALE;
posX[i] += velX[i] * TIME_STEP;
posX[i] = posX[i] < 0.0 ? 0.0 : posX[i];
}
if (posX[i] > POOL_WIDTH) {
velX[i] *= velX[i] > 0.0 ? -1.0 : 1.0;
velX[i] *= RESTITUTION_SCALE;
velY[i] *= FRICTION_SCALE;
posX[i] += velX[i] * TIME_STEP;
posX[i] = posX[i] > POOL_WIDTH ? POOL_WIDTH : posX[i];
}
if (posY[i] < 0.0) {
velY[i] *= velY[i] < 0.0 ? -1.0 : 1.0;
velY[i] *= RESTITUTION_SCALE;
velX[i] *= FRICTION_SCALE;
posY[i] += velY[i] * TIME_STEP;
posY[i] = posY[i] < 0.0 ? 0.0 : posY[i];
}
if (posY[i] > POOL_HEIGHT) {
velY[i] *= velY[i] > 0.0 ? -1.0 : 1.0;
velY[i] *= RESTITUTION_SCALE;
velX[i] *= FRICTION_SCALE;
posY[i] += velY[i] * TIME_STEP;
posY[i] = posY[i] > POOL_HEIGHT ? POOL_HEIGHT : posY[i];
}
}
}
效果截图
参考链接
[1]https://thecodeway.com/blog/?p=83
[2]https://en.wikipedia.org/wiki/Smoothed-particle_hydrodynamics
[3]《代码本色》Daniel Shiffman(第1、2、3、4、5章)
[4]https://blog.****.net/kuaxianpan2004/article/details/81353786