这种流体模拟基于Stam 的一篇论文。在第 7 页,他描述了平流背后的基本思想:
从两个网格开始:一个包含上一个时间步的密度值,另一个包含新值。对于后者的每个网格单元,我们通过速度场向后跟踪单元的中心位置。然后,我们从先前密度值的网格中线性插值,并将该值分配给当前网格单元。
建议代码。两个密度网格分别是d
和d0
,u
和v
是速度分量,dt
是时间步长,N
(全局)是网格大小,b
可以忽略:
void advect(int b, vfloat &d, const vfloat &d0, const vfloat &u, const vfloat &v, float dt, std::vector<bool> &bound)
{
float dt0 = dt*N;
for (int i=1; i<=N; i++)
{
for (int j=1; j<=N; j++)
{
float x = i - dt0*u[IX(i,j)];
float y = j - dt0*v[IX(i,j)];
if (x<0.5) x=0.5; if (x>N+0.5) x=N+0.5;
int i0=(int)x; int i1=i0+1;
if (y<0.5) y=0.5; if (y>N+0.5) y=N+0.5;
int j0=(int)y; int j1=j0+1;
float s1 = x-i0; float s0 = 1-s1; float t1 = y-j0; float t0 = 1-t1;
d[IX(i,j)] = s0*(t0*d0[IX(i0,j0)] + t1*d0[IX(i0,j1)]) +
s1*(t0*d0[IX(i1,j0)] + t1*d0[IX(i1,j1)]);
}
}
set_bnd(b, d, bound);
}
这种方法简洁且效果很好,但是实现对象边界对我来说很棘手,因为值是向后跟踪和插值的。我目前的解决方案是,如果旁边有一个空白空间(或多个空间),则简单地将密度推出边界,但这是不准确的,会导致密度增加,尤其是在角落和具有对角线速度的区域。仅在视觉上准确。我现在正在寻找“正确性”。
我的边界代码的相关部分:
void set_bnd(const int b, vfloat &x, std::vector<bool> &bound)
{
//...
for (int i=1; i<=N; i++)
{
for (int j=1; j<=N; j++)
{
if (bound[IX(i,j)])
{
//...
else if (b==0)
{
// Distribute density from bound to surrounding cells
int nearby_count = !bound[IX(i+1,j)] + !bound[IX(i-1,j)] + !bound[IX(i,j+1)] + !bound[IX(i,j-1)];
if (!nearby_count) x[IX(i,j)] = 0;
else
x[IX(i,j)] = ((bound[IX(i+1,j)] ? 0 : x[IX(i+1,j)]) +
(bound[IX(i-1,j)] ? 0 : x[IX(i-1,j)]) +
(bound[IX(i,j+1)] ? 0 : x[IX(i,j+1)]) +
(bound[IX(i,j-1)] ? 0 : x[IX(i,j-1)])) / surround;
}
}
}
}
}
bound
是具有行和列的布尔0
向量N+1
。bound
通过将单元格坐标设置为,在主循环之前设置边界对象1
。
该论文含糊地说“然后我们只需在set_bnd()
例程中添加一些代码,以根据其直接邻居的值填充被占用单元格的值”,这就是我正在做的事情。我正在寻找一种更准确地实现边界的方法,即具有非流体固体边界,并可能最终支持多种流体的边界。视觉质量比物理正确性重要得多。