模拟流体粒子运动
源码：https://files.cnblogs.com/flash3d/fluid.rar
拖动框框，可以让流体出现惯性现象（是这说法么？+_+）。
今儿在一个小日本网站上看到人家用as模拟流体运动，炫目无比，于是模仿做了个类似的效果。
下面主要讲其中涉及到的数学物理知识。
关于流体运动，上网查了下，捞到一些公式。由于用用在程序用模拟粒子运动，而不是进行精确的科学计算，所以，本人特意将公式化为最简，将一些能默认的系数默认，能忽悠的函数忽悠，得出效率和效果比较平衡的算法（其中两个光滑核函数参考小日本的）。
所谓模拟流体运动，并不是说每个粒子代表一个分子，而是比如你将一定量的羽毛放在风中，让其飞扬，羽毛的运动就能近似模拟风的运动。在下在构思的时候，曾试图通过范德华力这些分子层面的东西来模拟，走了一些弯路。
作为一个放入风中的羽毛，它会受到什么力呢？
羽毛会受到重力，来自其他羽毛带来的压力（并非直接接触，而是通过驱动风），风的阻力。
将这些力合成，就是羽毛所受合力。
知道所受合力，根据牛顿定律，就能求出速度和位移。
问题I在于如何求这三个力？
重力好求的，压力和阻力怎么求呢？
网上找了一些公式。
一个羽毛受到所有影响范围内的羽毛的作用，我们认为被作用的那个羽毛，它的各个属性被影响后的值，一般遵循以下规律（科学家研究的）：
假设被作用的那个羽毛叫F，所有影响那个羽毛分别为F0 F1 F2...FN。
那么羽毛FJ对羽毛F的A属性的影响是这么计算的：F.AJ=FJ.A*FJ.M/FJ.RO*W(DIST,RANGE)
其中M属性为质量，RO属性为密度，DIST为两个羽毛之间的距离，RANGE为影响范围的半径，W函数为光滑核函数。
光滑核函数是一类满足特定性质的函数，其通式为W(H,R)，其函数值随H的绝对值增大而衰减，直到H的绝对值大于等于R的时候，函数值就为0啦。光滑核函数还有一些性质，与我们的程序无关，比如他是一个偶函数，他的积分为1。要注意下，这个通式，他的自变量只有H一个，R是一个常数。
光滑核函数的意义在于，它可以描述一个粒子对周围影响的强烈程度和周围事物与该粒子的距离的关系。
于是，F.A=F.A0+F.A1+F.A2...+F.AJ+...+F.AN
首先计算粒子所受压力F.FP。
一个粒子所受压力是周围粒子固有的压力的累加。那么粒子固有压力是什么呢？就是压强啦。
高中学过气体压强公式，P=K(RO-RO')，RO-RO'是内外密度差，K是某热力学常数，P就是压强。
所以，要求所受压力，还必须求出每个粒子所在位置的密度。动用刚才那个属性叠加公式。
F.ROJ=FJ.RO*FJ.M/FJ.RO*W(DIST,RANGE)
密度能约掉，得到F.ROJ=FJ.M*W(DIST,RANGE)
其中，每个粒子的质量默认为1，于是F.ROJ=W(DIST,RANGE)
这里，我再讲一下光滑核函数，其是一类函数，满足其给定性质的均被称为光滑核函数。所以每个属性的影响叠加所用的光滑核函数也是不一样的，具体取哪个函数，是要看要叠加的属性的性质，或者看你需要的效果了。
这里我参考了小日本程序里面的光滑核函数：W(DIST,RANGE)=(1-DIST/RANGE)^2 自变量范围0到RANGE
这样 我们就能讲两个粒子的距离求出来，从而算出粒子间相互影响的量，并累加到粒子的密度中去。
代码实现
for (i = 0; i < numParticles; i++)
{
pi = particles[i];//取外循环要操作的粒子
pi.numNeighbors = 0;//将粒子邻居至空
pi.density = 0;//将粒子密度属性初始化
/*内循环，用握手规则减少循环次数*/
for (j = 0; j < i; j++)
{
pj = particles[j];//取内循环粒子
dist = Math.pow(pi.x - pj.x, 2) + Math.pow(pi.y - pj.y, 2);//计算两个粒子的距离平方
/*如果两个粒子距离在范围内进入判断*/
if (dist < RANGE2)
{
dist = Math.sqrt(dist);//算出距离
weight = Math.pow(1 - dist / RANGE, 2);//算出光滑核函数的值
//将值累加到两个粒子的密度里
pi.density += weight;
pj.density += weight;
//互加为邻居
pi.neighbors[pi.numNeighbors++] = pj;
pj.neighbors[pj.numNeighbors++] = pi;
}
}
有了密度，我们就能算压强啦，P=K(RO-RO')
RO'是流体在静态的时候全局密度，在代码中对应为DENSITY
代码实现
/*算出粒子的压力*/
for (i = 0; i < numParticles; i++)
{
pi = particles[i];
if (pi.density < DENSITY) pi.density = DENSITY;
pi.pressure = pi.density - DENSITY;//粒子压力等于K(粒子密度-静态密度)，K默认为1
}
有了压强，终于能计算粒子对其他粒子的压力了。
这里聪明的孩子马上意识到，压强是标量，压力是矢量，标量累加怎么就成了矢量了呢？
公式上的做法是将压强做一介微分和，也就是哈密顿算子（这个可以去百度），不过，我们程序上的做法是先算标量，再处理矢量。因为力的方向，其实我们早就知道了，从发力粒子指向受力粒子。
所以得到的公式就是F.PJ=FJ.P*FJ.M/FJ.RO*W(DIST,RANGE)
不过不幸的是，这个公式是“不平衡”的，也就是说，位于不同压强区的两个粒子之间的作用力不等，所以计算中一般使用双方粒子压强的算术平均值代替单个粒子的压力
F.PJ=(FJ.P+F.P)/2*FJ.M/FJ.RO*W(DIST,RANGE)
其中光滑核函数W(DIST,RANGE)=1 - DIST / RANGE
质量为1
F.PJ=(FJ.P+F.P)/2/FJ.RO*(1 - DIST / RANGE)
求出作用力之后，将其分成x和y两个分量分别累加
其中还要乘上压力系数也就是表面积。
F.PJ=F.PJ*PRESSURE
代码实现
dx = pi.x - pj.x;
dy = pi.y - pj.y;
dist = Math.sqrt(Math.pow(dx, 2) + Math.pow(dy, 2));//计算两个粒子的距离
weight = 1 - dist / RANGE;//光滑核函数
pressure = weight * (pi.pressure + pj.pressure) / (2 * pj.density) * PRESSURE;//计算粒子对粒子的压力标量
/*将压力分成向量的两个值累计到受力中*/
dist = 1 / dist;
dx *= dist;
dy *= dist;
pi.fx += dx * pressure;
pi.fy += dy * pressure
然后是计算阻力，阻力是粘度系数乘上运动速度差。所以，我们对运动速度差进行属性叠加。
F.VJ=FJ.V*FJ.M/FJ.RO*W(DIST,RANGE)
F.VJ=(F.v-FJ.v)*FJ.M/FJ.RO*W(DIST,RANGE)
其中光滑核函数W(DIST,RANGE)=1 - DIST / RANGE
于是化简得到F.VJ=(F.v-FJ.v)/FJ.RO*(1 - DIST / RANGE)
最后得到阻力F.VJ=F.VJ*VISCOSITY
在代码中，由于速度差是矢量，所以，乘上速度的时候，分别在x和y两个方向处理。
代码实现
viscosity = weight / pj.density * VISCOSITY;//计算除速度外的阻力系数
//计算粒子相对速度的向量
dx = pi.vx - pj.vx;
dy = pi.vy - pj.vy;
//乘入速度因子得出粒子对粒子的压力累加入受力
pi.fx -= dx * viscosity;
pi.fy -= dy * viscosity
最后剩下其他外力作用了
其中有重力，瓶子被移动引擎的加速度，还有碰到墙壁反弹的力
重力直接累加到受力内
瓶子的加速度反向累加到受力内
墙壁的弹力，这里将墙壁想象成弹簧，弹性系数K，那么如果粒子刚接触到墙壁被视为开始接触弹簧，那么粒子超出墙壁的范围，被视为弹簧的压缩量L，那么产生的力F=-K*L，这个力是被即时累加到速度里面的。
重力和瓶子加速度累加的代码
for (i = 0; i < numParticles; i++)
{
pi = particles[i];
//考虑外部瓶子受力，将其反向累计
pi.fx -= bottle_fx;
pi.fy -= bottle_fy;
pi.fy += Main.GRAVITY;
pi.move();//将受力作用成速度和位移
}
撞墙后反馈力的代码（弹性系数为1）
//如果撞墙，速度为0，如果穿过墙，则受到和超过部分相同的力的反馈
if (x < 4) vx += 4 - x;
if (y < 4) vy += 4 - y;
if (x > 461) vx += 461 - x;
if (y > 461) vy += 461 - y;
至于其他部分，如力怎么累加到速度，速度怎么累加到位移之类的问题，我就先不讲了。祝大家周末愉快！！
感谢四楼提醒，本文部分公式参考自http://www.thecodeway.com/blog/?p=1557
