计算流体模拟理论2

二维的NS-方程:

 

 

这个方程必定要拆分红部分才能解出来。

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

 

<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>推出压力方程:

梯度压力方程:

离散，下面有具体的压力梯度离散过程。

 

这个为压力梯度方程。只要求出压力p就能够获得无散速度场。

接下来不可压缩流体条件:

 

使用中心差分法离散。这样没有null-space,由于速度场特殊的储存方式

 

接下来推出怎么样获得下一个时间步上的流体速度场是无散度,这里有个技巧，咱们不能直接使用5.4式。可是:

首先把速度散度公式写成下一步时间的离散形式:

 

把梯度压力方程带入能够获得:

 

观察此方程，右边就是不可压缩 负的速度梯度, 这部分叫作rhs,右手方程,这个是已知的。

因此方程中只有压力p未知。求出压力p便可。

 

 

 

 

<2>求压力右手方程。一般就叫作RHS,就是负的速度场的散度。-divergence(u)

    这个是简单的。

 

 

<3>求解压力.

 

 

    这个有好多方法求，因为是一个AX=b 的形式，大部分文章是以PCG/GAUSS-SEIDEI/GAUSS-SEIDEI SOR/JACOBI方式

    因为在numpy方便，我直接把方程用切片法作了。

 

 

 <4>把求的pressure field 带入压力梯度更新方程，求出无散度的速度。

       Pressure gradient 压力梯度方程离散:

这个方程我以为才是套路中的套路。有了这个压力，直接带入这个公式 ，就能够求出无散的流体。