计算流体模拟理论2

二维的NS-方程:

这个方程一定要拆分成部分才能解出来。

这里面我感觉只要把泊松方程解法搞定，基本快出山写最简单的 "计算流体" 完全没问题

以下是做了个初始的source field,用python numpy 先快速撸了一遍算法。

并且重新实现3d版本在Houdini中，houdini有更好的可视化.

velocity advection 是RK2,

density advection 一阶的,

然后发现density field确实跟不上velocity field

从零蛋编写无粘流体最简单的路线:

<0>为什么要使用MAC-GRID:

按照作者的说法，使用这种grid,最大的优势是解决了中心差分法零空间问题。

速度场一定要存在face center.

所以在Houdini看vel field 是一定存在grid face center上的,而不是grid cell center!

而density,temperature....这类场是存在grid cell上的。

总结下一个inviscid fluid(欧拉无粘流体)实现过程:

<1>对流。

最重要的是理解材质导数.为什么这个玩意能等于0? 详见欧拉与拉格朗日的观点。

这个差分法是万万不可取，因为会出现null-space,试着想想如果有3个点，中间的点高一点，两边一样高，那么用中心差分法算出来中间点的斜率将是0.

所以使用无条件稳定的半拉格朗日,这个其实完全是解ODE问题了。而不是PDE了

比如要得到xp的值，只需用用向前euler 向后追踪法。当然这个方法比较废物。作者推荐RK家族的算法比如:RK2:

def RK2(x, y, dt, ugrid, vgrid):
    nu = neibours_value(ugrid, x, y, "u")
    nv = neibours_value(vgrid, x, y, "v")

    xmid = x - 1.0 / 2.0 * dt * nu
    ymid = y - 1.0 / 2.0 * dt * nv

    umid = neibours_value(ugrid, xmid, ymid)
    vmid = neibours_value(vgrid, xmid, ymid)

    x -= dt * umid
    y -= dt * vmid

    return (x, y)

(你去看houdini上面的gas advect 还有RK5)

得到位置直接用bilerp()插值法求出xp的量。

def bilerp(f00, f10, f01, f11, tx, ty):
    """ FIGURE : first lerp in top x,then bottom x, then along y axis
    f00*----------.tx-----------*f10
    |             |             |
    |             |             |
    |             .ty           |
    |             |             |
    |             |             |
    f01*----------.tx-----------*f11
    """
    return lerp(lerp(f00, f10, tx), lerp(f01, f11, tx), ty)

那么这个量就是作为下一个时间步的量。

测试这个最简单的就是创建一个简单的恒定区域网格速度，然后让自己的初始的density是不是根据semi-lagrangian方法能advection.

<1.1>推出压力方程:

梯度压力方程: