計算流體模擬理論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 壓力梯度方程離散:

這個方程我以爲纔是套路中的套路。有了這個壓力，直接帶入這個公式 ，就能夠求出無散的流體。