第六講：流體模擬- 流體力學(xué)基礎(chǔ)

雖然綠燈沒怎么為我亮過仆葡，但我還是對(duì)生活充滿熱情习劫，這就是我理解的年輕 —— 汪曾祺《人間草木》

所謂的流體就是流動(dòng)的氣態(tài)涌庭、液態(tài)甚至固態(tài)的物體损敷，比如空氣、水屹逛、沙子等础废，這是經(jīng)典物理里面目前還依然沒能成功攻克的難題。

普通物體的形變可以分成三種：

拉伸Stretch
壓縮Compression
剪切形變Shear

而流體一個(gè)方面很難做到前面兩種形變罕模，另一個(gè)方面則十分容易發(fā)生第三種形變评腺，因?yàn)檫@個(gè)性質(zhì)的存在，導(dǎo)致計(jì)算機(jī)模擬變得很困難淑掌，因?yàn)椋?/p>

前者意味著當(dāng)受力之后由于抵抗拉伸與壓縮的能力很強(qiáng)蒿讥，這個(gè)受力狀態(tài)將會(huì)以擾動(dòng)的方式很快的傳播開來芋绸，也就是說任何一個(gè)擾動(dòng)將會(huì)導(dǎo)致整個(gè)流體中的所有位置都會(huì)受到影響（即使應(yīng)用了前面的Projection Dynamics，也很難計(jì)算出結(jié)果）
后者意味著擾動(dòng)會(huì)容易引起旋渦（湍流）摔敛，而旋渦的特性在于大旋渦會(huì)容易變成小旋渦，即旋渦會(huì)不斷分解苦酱，從而導(dǎo)致使用離散化方式計(jì)算的計(jì)算機(jī)模擬無法完全模擬出流體動(dòng)力學(xué)的細(xì)節(jié)。

流體問題可以從空間跟時(shí)間兩個(gè)角度來看待给猾，尺度有：

空間：從大到小敢伸，分子->平均自由程（分子相撞之前平均的移動(dòng)距離）->system scale（空氣流過飛機(jī)翅膀經(jīng)過的距離）
時(shí)間：t0(分子碰撞時(shí)間）-> t1（水分子布朗運(yùn)動(dòng)碰撞前的平均經(jīng)歷時(shí)間）->tc（空氣流過飛機(jī)翅膀需要的時(shí)間）

根據(jù)不同的尺度，有不同的模擬方式：

分子尺度：分子動(dòng)力學(xué)（合成钓丰、制藥）
自由程尺度：統(tǒng)計(jì)意義上的效果琢歇，Kinetic theory李茫，這是最近一段時(shí)間開始流行的模擬角度魄宏，Lattice Boltzmann Model
大尺度（system scale）：NS方程，這是計(jì)算機(jī)模擬在過去二三十年專注的區(qū)域

游戲模擬中常用的方式是后面兩種予跌，這兩種模擬方式各有千秋匕得，這里會(huì)優(yōu)先介紹大尺度下的模擬方式汁掠，后面有時(shí)間再介紹Kinetic theory方法考阱。

宏觀尺度下流體模擬的一個(gè)核心思想是Field View秽之，即場(chǎng)視角考榨，也就是說對(duì)于流體中任意一點(diǎn)的物理量（如速度河质、溫度、密度乐尊、壓強(qiáng)等）并不是固定的扔嵌，而是與對(duì)應(yīng)點(diǎn)位置相關(guān)的一個(gè)數(shù)值。

NS方程本質(zhì)上是牛頓第二定律：

$f = m \frac{d u}{d t}$

在流體中牺弄，速度是一個(gè)向量場(chǎng) $u(x, t)$ 势告，數(shù)值隨時(shí)間、位置發(fā)生變化回溺。

在流體中遗遵，我們對(duì)速度u進(jìn)行求導(dǎo)允粤，得到的不是加速度（回顧下加速度的定義类垫，單位時(shí)間內(nèi)速度的變化量），舉個(gè)例子购撼，某個(gè)小船在水面上行駛，其加速度可以按照下面方式進(jìn)行推導(dǎo)晃跺。

首先原始的速度為 $u_0 = u(t_0, x_0)$ ，經(jīng)過 $\Delta t$ 之后烹玉，其速度就變成了 $u(t_0 + \Delta t, x_0 + u_0 \Delta t)$ ，加速度就變成了：

$a = \lim_{\Delta t \rightarrow 0}\frac{u(t_0 + \Delta t, x_0 + u_0 \Delta t) - u_0}{\Delta t}$

用泰勒展開就變成了：

$a = \frac{D u}{D t} = u \nabla_x u + \frac{\partial u}{\partial t}$

這個(gè)加速度叫Material Derivative继效，這個(gè)公式中，前面一項(xiàng)表示的是時(shí)間不變凡简，空間上的速度變化帜乞，后面一項(xiàng)則是空間不變（即小船待在原地不動(dòng)）隨著時(shí)間流逝帶來的速度的變化挖函。

加速度有了，那么質(zhì)量怎么來定義呢必怜，流體沒有單個(gè)點(diǎn)的只能梳庆，只能用單位體積的質(zhì)量（也就是密度 $\rho$ ）來表示，那么我們有：

$\frac{D u}{D t} * \rho = \frac{f}{V}$

$u \nabla_x u * \rho + \frac{\partial u}{\partial t} * \rho = \frac{f}{V}$

上面的V是流體的體積（ $m = \rho * V$ )更米，而 $\frac{f}{V}$ 可以看成是力的密度：

對(duì)于重力而言，其密度就是單位質(zhì)量上的力栏笆，也就是 $\rho g = \frac{m g}{V}$
壓力，以2D為例七婴，在流體中任意框選出一個(gè)box中，四邊都有壓強(qiáng)户盯，且左右上下壓強(qiáng)不一定相同吗伤，從而形成一個(gè)力足淆，這個(gè)力就是壓強(qiáng)的梯度： $- \nabla P$
viscosity（抵抗剪切形變的力），這一項(xiàng)就直接給出公式為 $\mu \nabla \cdot \nabla u$ 丹鸿，被稱為viscosity term，其中的 $\mu$ 稱為dynamic viscosity门怪。

在上面這幾項(xiàng)力的作用下，我們就有最終的NS方程：

$\frac{\partial u}{\partial t} * \rho+ \rho(u \nabla_x) u = \rho g - \nabla P + \mu \nabla \cdot \nabla u$

在很多流體模擬中， $\mu$ 很小趨近于0嚼锄，所以通常直接省略最后面一項(xiàng)。

這個(gè)方程中沧侥，u是不知道的，壓強(qiáng)P也是不知道的旺罢，密度也是不知道的正卧，但是只有一個(gè)方程，因此是不足以求解的窘行，要想求解抽高，還需要補(bǔ)上兩個(gè)方程。

根據(jù)質(zhì)量守恒方程，我們有（難受莹桅，這里沒聽到诈泼，看了半天也沒看出來是如何推導(dǎo)的，后面有機(jī)會(huì)再補(bǔ)上對(duì)應(yīng)的說明）：

$\fracv3mex31{dt}\int_{V_0}\rho(x) dx = - \oint \rho u dA$

這個(gè)公式可以變化成（也錯(cuò)過了這個(gè)變化的推導(dǎo)，-_-||）：

$\int_{V_0} \fracc2isk7z{dt}\rho(x) dx = - \int_{V_0} \nabla(\rho u) dx$

最后我們將兩者合并，得到：

$\int_{V_0} [\frac2y8mits{dt}\rho(x) + \nabla(\rho u)] dx = 0$

由于質(zhì)量不能為負(fù)偏灿，所以我們有：

$\fracnjmel7t{dt}\rho(x) + \nabla(\rho u) = 0$

也就是：

$\fractphhdhv{dt}\rho(x) + u \nabla \rho + \rho \nabla u = 0$

而等式左邊前面兩部分（為什么忿墅？參考前面提到的Material Derivative）又可以表示成： $\frac{D \rho}{D t}$ ，所以我們有：

$\frac{D \rho}{D t} + \rho \nabla u = 0$

而基于流體的不可壓縮性原理，我們可以得到速度的梯度為0：

$\nabla u = 0$

好了呼畸，有了上面這些條件，我們就可以著手來求解NS方程了儒陨，這里介紹一種經(jīng)典的解法（也就是Houndini中的解法），即通過time integrator的方式來求解笛园，對(duì)NS方程進(jìn)行分析，我們可以知道棵红，整個(gè)方程中只有第一項(xiàng) $\frac{\partial u}{\partial t}$ 是跟時(shí)間有關(guān)系的哟冬，也就是說可岂，我們可以將這一項(xiàng)按照時(shí)間寫成如下形式：

$\frac{\partial u}{\partial t} \approx \frac{u^{n+1} - u^n}{\Delta t}$

接下來只需要將NS方程中剩下的幾項(xiàng)估計(jì)出來稚茅，就能求得 $u^{n+1}$ ，但是由于這里有多項(xiàng)需要估計(jì)欺税，且壓強(qiáng)是不可知的，我們要如何進(jìn)行估計(jì)與求解呢歼秽？

這里常用的一種方法叫做Operator Splitting，簡(jiǎn)單來說，如果我們有一個(gè)簡(jiǎn)單的微分方程：

$\frac{dx}{dt} = f(x) + g(x)$

根據(jù)顯示歐拉方程做院，我們轉(zhuǎn)換一下就是：

$x^{n+1} = x^n + f(x^n) \Delta t+ g(x^n) \Delta t$

這個(gè)公式可以將之拆分成兩步來求得：

$x^* = x^n + f(x^n) \Delta t$

$x^{n+1} = x^* + g(x^n) \Delta t$

我們知道柑营，顯式歐拉跟隱式歐拉之間的區(qū)別在于右邊的函數(shù)使用的x是 $x^n$ 還是 $x^{n+1}$ 酒奶，其實(shí)都是一種近似，而在這里我們也可以用其他的x來近似另伍，所以這里我們可以將上面的第二個(gè)式子改寫成：

$x^{n+1} \approx x^* + g(x^*) \Delta t$

所以問題就變成了温艇，我們需要求解兩個(gè)簡(jiǎn)單的微分方程：

$\frac{dx}{dt} = f(x)$

$\frac{dx}{dt} = g(x)$

回到NS方程來說，我們就可以采用同樣的思想將之拆解成四個(gè)步驟的簡(jiǎn)單微分方程的求解（實(shí)際上，最后一項(xiàng)粘稠性非常小绣否，可以先暫時(shí)忽略）。

首先，我們先看左邊的一項(xiàng)：

$\frac{\partial u}{\partial t} \rho = -u \nabla u$

右邊的 $\nabla u$ 是u的梯度苹支，也就是空間中的微分究反，一種最簡(jiǎn)單的方式（經(jīng)典軟件使用的方法會(huì)更復(fù)雜）就是將一塊區(qū)域想象成由格子來定義的場(chǎng)精耐，每個(gè)采樣點(diǎn)都落在格子的頂點(diǎn)上向胡，那么通過 $\frac{u_{n+1} - u_n}{\Delta s}$ 就能求解出 $\nabla u$

之后就可以用顯式歐拉的方法來求解出 $u^{n+1}$ 了专执，不過由于顯式歐拉本身不穩(wěn)定攀痊，會(huì)不斷震蕩跑飛，而如果用隱式歐拉棘街，則需要求解非線性方程，加大了復(fù)雜度承边。

而上面這個(gè)方程其實(shí)可以表示成：

$\frac{Du}{Dt} = 0$

可以想象有個(gè)particle在流體中流動(dòng)遭殉，如上圖中所示，從格子上頂點(diǎn)流動(dòng)到右上角的紅點(diǎn)上博助，也就是說险污，上面公式的物理意義在于富岳，在經(jīng)過這個(gè)流動(dòng)之后蛔糯，兩者的速度并不會(huì)發(fā)生變化，也就是說我們可以采用一種叫做Semi-Lagrangian的方法來求解窖式。

以上圖進(jìn)行說明蚁飒，假設(shè)我們想要知道上圖中綠色方框中的一點(diǎn)的速度，我們假設(shè)這一點(diǎn)是來自于圖上五角星位置萝喘，那么我們就可以知道上一個(gè)timestep上五角星位置跟當(dāng)前timestep上綠色box位置的速度是相同的淮逻，這里有兩個(gè)問題：

如何從綠色方框位置反推出五角星位置，這就是semi的由來蜒灰，也就是根據(jù)當(dāng)前綠色方框的位置朝著時(shí)間相反的方向以當(dāng)前的速度移動(dòng)一段距離弦蹂，就知道五角星位置
找到五角星位置后，我們?cè)趺粗牢褰切俏恢玫乃俣惹拷眩挥捎谖覀冎理旤c(diǎn)上的速度凸椿，我們就可以通過插值得到五角星位置的速度了。

這種算法的好處在于：

整個(gè)系統(tǒng)是穩(wěn)定的翅溺，因?yàn)槲覀冎皇菍⒛骋稽c(diǎn)的速度拷貝到另一點(diǎn)脑漫，而不會(huì)出現(xiàn)速度的增加，所以不會(huì)爆掉
求解方便咙崎，不需要求解非線性方程
Energy Disappearance會(huì)小一些优幸，隱式歐拉會(huì)隨著時(shí)間的推動(dòng)，能量會(huì)不斷磨損褪猛，而如果這個(gè)過程過于迅速的話网杆，流體的效果就一閃而逝了，看起來就不真實(shí)，而這里的Semi-Lagrangian方法就可以稍微降低消失的速度碳却，使得效果更好一點(diǎn)队秩。

經(jīng)過這一步計(jì)算之后，我們就根據(jù) $u^n$ 得到了 $u^*$ 昼浦，接下來我們考慮右邊的第一個(gè)式子：

$\frac{\partial u}{\partial t} = g$

由于右邊是一個(gè)常量馍资，因此這個(gè)公式可以直接使用顯式歐拉求解得到，從而得到了 $\tilde u$ 关噪，接著再看下一個(gè)式子：

$\frac{\partial u}{\partial t} = -\frac{1}{\rho} \nabla P$

由于我們不知道壓強(qiáng)P鸟蟹，所以需要將流體的不可壓縮特性用起來，也就是說使兔，用近似公式建钥，我們有：

$u^{n+1} = \tilde u -\frac{1}{\rho} \nabla P \Delta t$

而由于 $\nabla u = 0$ （注意，這里的梯度是空間上的差分或者微分）虐沥，也就是：

$\nabla u^{n+1} = \nabla \tilde u -\frac{\Delta t}{\rho} \nabla \cdot \nabla P = 0$

根據(jù)這個(gè)式子我們有：
$\frac{\rho \cdot \nabla \tilde u}{\Delta t} = \nabla \cdot \nabla P$

這個(gè)公式叫做Pressure Projection锦针，也就是說，我們有一個(gè)壓強(qiáng)P置蜀，這個(gè)參數(shù)存在于上圖中的格子的頂點(diǎn)上奈搜，而這個(gè)數(shù)據(jù)我們就可以通過有限差分公式估計(jì)出其二階導(dǎo)數(shù)，而由于這個(gè)方程中我們不知道壓強(qiáng)P盯荤，而密度 $\rho$ 是常量（不可壓縮）馋吗， $\Delta t$ 也是已知常量，而 $\tilde u$ 是上一步計(jì)算得到的結(jié)果秋秤，所以我們就可以求解得到P（線性方程組求解）宏粤，將這個(gè)代入前面的公式，我們就能夠得到 $u^{n+1}$ 灼卢。

$u^{n+1} = \tilde u -\frac{1}{\rho} \nabla P \Delta t$

這了需要注意的是绍哎，我們用P來估計(jì) $u^{n+1}$ ，采用的是 $\nabla P$ 鞋真，而這個(gè)有限差分估計(jì)出來的結(jié)果實(shí)際上是上圖中頂點(diǎn)中心處的數(shù)值崇堰，用這個(gè)數(shù)值估算出來的 $u^{n+1}$ 也就是中心位置的速度，這就是為什么所有的流體算法將上圖中的格子叫做Staggered Grid涩咖，也就是速度是存放在格子的每條邊的中點(diǎn)位置上的海诲。

求解就是這個(gè)過程，這里需要說幾點(diǎn)：

$u^*$ 跟 $\tilde u$ 的求解是解耦的檩互，因此可以并行完成特幔，且可以放在GPU上來實(shí)現(xiàn)（每個(gè)格點(diǎn)是獨(dú)立的）
最后P的求解是需要求解一個(gè)線性方程組，需要將整個(gè)Grid上的每個(gè)格點(diǎn)都耦合在一起闸昨，因?yàn)榱黧w的不可壓縮性蚯斯，因此在任意位置給流體一個(gè)壓強(qiáng)或者擾動(dòng)薄风，就會(huì)在無限短的時(shí)間內(nèi)導(dǎo)致全盤的影響，因此需要求解線性方程組拍嵌，這個(gè)計(jì)算會(huì)很麻煩村刨；目前經(jīng)典解算方法都是采用這種方式，比如Houndini撰茎，雖然可以通過GPU來甲酸，但是因?yàn)椴惶貌⑿写蛲荩堑倪^程龄糊，在這一步就會(huì)很慢。
彈性體的時(shí)候募疮，我們說過炫惩，一個(gè)約束條件下的求解，可以看成是一個(gè)優(yōu)化問題阿浓；這里我們可以采用同樣的方式來考慮這個(gè)問題他嚷，假設(shè)我們沒有流體的不可壓縮條件，那么我們就可以忽視NS方程中的P那一項(xiàng)芭毙，也就是說我們只需要經(jīng)過前面兩步筋蓖，不考慮最后P的求解即可，但實(shí)際上不可壓縮性是實(shí)際存在的退敦，因此我們需要求解如下優(yōu)化方程：

$min \frac{1}{2}(u^{n+1} - \tilde u)\rho(u^{n+1} - \tilde u)$

也就是說 $(u^{n+1} - \tilde u)$ 最小的時(shí)候我們得到最優(yōu)解粘咖，但是由于不可壓縮性，我們需要滿足：

$\nabla u^{n+1} = 0$

而這種條件我們就需要使用Lagrangian Multiplier方法來求解侈百，其求解得到的結(jié)果跟前面Operator Splitting得到的結(jié)論是一樣的瓮下，從這個(gè)角度來看，壓強(qiáng)P不是一個(gè)真正的物理量钝域，而是一個(gè)力讽坏，其作用是使得上面的流體具有不可壓縮性，也就是說用來反抗流體中的每一點(diǎn)的速度變化的（慣性力）例证。

另外需要介紹的是流體在一個(gè)容器內(nèi)流動(dòng)的邊界條件：

邊界上每一點(diǎn)的法線方向上的速度為0（從而避免流體穿過容器）
如果流體中有一個(gè)擾動(dòng)物路呜，比如一個(gè)葉片，那么葉片上的速度與葉片上流體的速度是一樣的

下面來介紹一下游戲中常用的SPH（Smoothed Particle Hydrodynamics织咧，最開始來自于一些核試驗(yàn)?zāi)M領(lǐng)域）算法拣宰，因?yàn)镹S中的線性方程組的求解很難做到實(shí)時(shí)，而很多電影特效中常用的就是NS Solver方法（為什么烦感，因?yàn)榫雀邌嵫采纾浚?/p>

SPH（Smoothed Particle Hydrodynamics）

SPH的思想是，我們將流體看成是一堆particle手趣，想象particle可以在其中移動(dòng)晌该，我們可以為每個(gè)particle賦予一套物理量肥荔，假設(shè)我們稱之為 $f_i$

假設(shè)我們知道了每個(gè)粒子上的物理量，那么對(duì)于任意一個(gè)位置x朝群，我們就需要能計(jì)算出其物理量 $f_x$ 燕耿，從而計(jì)算得到其gradient，甚至二階gradient（Laplacian）姜胖，我們就能將這些數(shù)據(jù)代入NS方程中誉帅，完成粒子上速度的更新，從而用這個(gè)速度對(duì)粒子進(jìn)行移動(dòng)右莱，就能得到下一個(gè)timestep的粒子蚜锨。

也就是說，其中心思想是根據(jù)粒子的物理量慢蜓，求解出整個(gè)場(chǎng)域的物理量亚再。

這個(gè)算法中最核心的東西就是一個(gè) $\delta(x)$ ，這不是一個(gè)函數(shù)晨抡，而是一個(gè)分布：

在x=0處有一個(gè)非零值氛悬，其他位置結(jié)果為0，定義可以表示為：

$f(x) = \int f(y) \delta(y-x) dx$

實(shí)際上這個(gè)分布是一個(gè)理想的分布耘柱，現(xiàn)實(shí)中是沒有的如捅，也不好給出表達(dá)式，而我們可以近似一下调煎，給出如下的近似分布：

$f(x) \approx \int f(y) W(|y-x|,h) dx$

這里的W稱之為kernel函數(shù)伪朽，h用來控制上面近似分布的寬度，h越大汛蝙，分布越寬烈涮。

有了這個(gè)近似的分布函數(shù)，我們就可以表達(dá)f(x)窖剑，也就是對(duì)于任何的x坚洽，我們可以將f(x)表達(dá)為一個(gè)求和項(xiàng)：

$f(x) = \sum_i f_i W(|x-x_i|, h) w_i$

最后的 $w_i$ 是權(quán)重，而W是kernel函數(shù)西土，我們是知道其表達(dá)式的讶舰，而這個(gè)函數(shù)在|y-x|大于某個(gè)范圍的時(shí)候，這個(gè)函數(shù)返回為0需了。

這里我們先不管f表示的是什么跳昼，就先看成是任意的物理量，我們推導(dǎo)出任意位置的f(x)肋乍， $\nabla f$ 跟 $\nabla \cdot \nabla f$ 鹅颊，就可以將實(shí)際的物理量（比如速度、壓強(qiáng)）代入進(jìn)去墓造，并放在NS方程中堪伍，就能得到 $u^{n+1}$ 锚烦，并根據(jù)這個(gè)速度對(duì)粒子進(jìn)行移動(dòng)，就能得到下一個(gè)timestep的初始數(shù)據(jù)帝雇，完成了一輪迭代涮俄。

令r = y - x，我們可以給出W的一個(gè)示例函數(shù)：

$W(r,h)=\sigma \begin{cases} 6(q^3-q^2)+1 & {0 \leq q \leq \frac{1}{2}}\\ 2(1-q)^3 & { \frac{1}{2} < q \leq 1}\\ 0 & {otherwise} \end{cases}$

其中 $q=\frac{|r|}{h}$ 尸闸。

上面已經(jīng)給出了f(x)的求和項(xiàng)彻亲，那么其gradient就可以在前面公式兩邊加一個(gè)gradient就行，那么由于 $f_i$ 是常量吮廉， $w_i$ 也是不變的苞尝，那么gradient就施加在W上，但是這樣的做法雖然形式上沒問題茧痕，實(shí)際上我們?cè)趫?zhí)行中就會(huì)發(fā)現(xiàn)運(yùn)行出來的結(jié)果就會(huì)爆掉（表現(xiàn)跟顯式歐拉一樣），這是因?yàn)榱Ｗ釉谝苿?dòng)過程中恼除，有些地方粒子會(huì)密集一點(diǎn)踪旷，有些會(huì)稀疏。

在密集的位置豁辉，這個(gè)積分變成求和的近似是相對(duì)準(zhǔn)確的令野，但是稀疏區(qū)域的近似就不準(zhǔn)確，就會(huì)導(dǎo)致數(shù)值上的震蕩徽级，即誤差會(huì)不斷累計(jì)气破，導(dǎo)致爆掉。

所以餐抢，這個(gè)給我們的啟發(fā)就是在密集區(qū)域现使，可以按照這種方式來求取，但是在稀疏的區(qū)域就需要另外處理旷痕，也就是首先我們得知道各個(gè)位置的粒子密度碳锈，我們需要先計(jì)算 $\nabla (\frac{f}{\rho})$ ，也就是用粒子的密度來對(duì)物理量進(jìn)行歸一化欺抗。

假設(shè)我們定義粒子的密度為：

$\rho_i = \sum_j m_j W_{ij}$

其中 $m_j$ 表示的是j點(diǎn)的質(zhì)量售碳， $W_{ij} = W(|x_i - x_j|,h)$

那么用鏈?zhǔn)椒▌t（復(fù)合函數(shù)的微分）展開就可以得到：

$\nabla (\frac{f}{\rho}) = \frac{\nabla f}{\rho} - \frac{f}{\rho^2}\nabla \rho$

變換一下就得到了：

$\nabla f = \rho \nabla (\frac{f}{\rho}) + \frac{f}{\rho}\nabla \rho = \rho \cdot (\nabla (\frac{f}{\rho}) + \frac{f}{\rho^2}\nabla \rho)$

右邊有兩個(gè)gradient，我們就可以用前面的求和（可能導(dǎo)致不穩(wěn)定的）公式分別來近似這兩個(gè)式子

$\nabla (\frac{f}{\rho}) = \sum_j \frac{f_j}{\rho_j^2} \nabla W_{j} w_j \rho_j = \sum_j \frac{f_j}{\rho_j^2} \nabla W_{j} m_j$

$w_j \rho_j = m_j$ 表示的是某個(gè)位置的質(zhì)量绞呈，可以用這個(gè)位置的密度乘上其權(quán)重算得贸人。

$\frac{f}{\rho^2}\nabla \rho = \sum_j m_j \frac{f}{\rho^2} \nabla W_{j}$

最后得到一個(gè)表達(dá)式(定義在任意一個(gè)粒子上的位置，而非任意位置）：

$\nabla f(x_i) = \nabla f_i = \rho_i \sum_j m_j(\frac{f_j}{\rho_j^2} + \frac{f_i}{\rho_i^2})\nabla_i W_{ij}$

而這個(gè)式子則能夠避免簡(jiǎn)單gradient導(dǎo)致的不穩(wěn)定問題佃声，因?yàn)橛眠@個(gè)公式計(jì)算艺智，可以保證動(dòng)量是守恒的。對(duì)于上面的式子圾亏，每個(gè)粒子的動(dòng)量之和有：

$\sum_i m_i u_i = \sum_i m_i(u_i -\frac{1}{\rho_i} \nabla P_i)$

上面式子中的壓強(qiáng)來自于NS方程（除以了等號(hào)左邊的密度）力惯， $\frac{1}{\rho_i} \nabla P_i$ 表示的是由于壓強(qiáng)（內(nèi)力）導(dǎo)致的速度變化碗誉。因?yàn)槿绻黧w的速度的變化是由于內(nèi)力（即壓強(qiáng)，準(zhǔn)確來說是壓強(qiáng)的梯度）導(dǎo)致的父晶，那么就可以保證動(dòng)量是守恒的哮缺。

移除同類項(xiàng)：

$\sum_i m_i \frac{1}{\rho_i} \nabla P_i = 0$

將前面 $\nabla f_i$ 的求和式子代入 $\nabla P_i$ 就有：

$\sum_i m_i \sum_j m_j (\frac{f_j}{\rho_j^2} + \frac{f_i}{\rho_i^2})\nabla_i W_{ij} = 0$

那么這個(gè)等式為什么成立呢？

同樣甲喝，我們就可以求得Laplacian式子尝苇，這個(gè)有很多個(gè)版本：

$\nabla^2 f(x_i) = \sum_j \frac{m_j}{\rho_j} f_j \nabla^2 W_{ij}$

另一個(gè)常用的版本為：

$\nabla^2 f(x_i) \approx 2(d+2) \sum_j \frac{m_j}{\rho_j} \frac{A_{ij} \cdot r_{ij}}{|r_{ij}|^2} \nabla_i W_{ij}$

其中d表示的是維度，2維則取值為2埠胖，另：

$A_{ij} = \frac{f_j}{\rho_j^2} + \frac{f_i}{\rho_i^2}$

有了一階二階梯度公式之后糠溜，我們?cè)賮砬蠼釴S方程：

$\frac{\partial u}{\partial t} * \rho+ \rho(u \nabla_x) u = \rho g - \nabla P + \mu \nabla \cdot \nabla u$

根據(jù)Operator Splitting的思想，我們分步進(jìn)行處理直撤，首先：

$\frac{Du}{Dt} = \frac{\partial u}{\partial t}+ (u \nabla_x) u = \frac{\mu}{\rho} \nabla \cdot \nabla u + g$

根據(jù)前面的推導(dǎo)非竿，右邊的式子中我們都是知道的（ $\mu$ 是dynamic viscosity），我們就能夠通過顯式歐拉的方法完成速度的求解谋竖。此外红柱，這個(gè)過程對(duì)于每個(gè)粒子而言是完全獨(dú)立的，因此可以放在GPU上來計(jì)算蓖乘。

之后根據(jù) $\frac{D \rho}{Dt}=0$ 計(jì)算出 $\nabla P$ 锤悄，壓強(qiáng)跟密度有很多種映射表達(dá)方式，這里選取最簡(jiǎn)單的一種：

$P = k (\rho - \rho_0)$

其中 $\rho_0$ 表示的是期望的密度嘉抒，選取這種方式的好處在于各個(gè)粒子之間的計(jì)算是相互獨(dú)立的方便GPU并行零聚。

另外，這個(gè)公式也就說明了雖然我們一直在說流體的不可壓縮性些侍，但實(shí)際上如氣體在空間中的流動(dòng)隶症，其密度是會(huì)發(fā)生變化的，而這個(gè)公式實(shí)際上相當(dāng)于添加了一個(gè)懲罰項(xiàng)岗宣，目的是維持流體的不可壓縮性沿腰。

當(dāng)然，學(xué)術(shù)界也有其他的表達(dá)式是能確保流體的不可壓縮性的狈定，但是這類的方法會(huì)導(dǎo)致需要進(jìn)行多輪迭代才能輸出結(jié)果颂龙，從而打破了GPU的并行性，這也是SPH方法尷尬的地方纽什，雖然是想要設(shè)計(jì)成GPU并行的措嵌，但是卻沒有辦法真正并行。

再之后就求解：

$\frac{D u}{Dt} = -\frac{\nabla P}{\rho}$

并根據(jù)這個(gè)式子完成速度u的更新芦缰。

最后根據(jù)：

$\frac{D x}{Dt} = u$

求得更新后的位置企巢。

下面來介紹一下Lattice Boltzman Method。

Lattice Boltzman Method

這個(gè)方法是上個(gè)世紀(jì)七十年代就已經(jīng)提出來的方法让蕾，這個(gè)方法來自于Kinetic Theory（統(tǒng)計(jì)物理學(xué)）浪规，之所以這么久以來一直不溫不火是因?yàn)檫@個(gè)方法不夠精確或听，但是2006年左右出現(xiàn)了一系列的方法提升了這個(gè)方法的精確性。

而且這個(gè)方法本身沒有迭代的過程笋婿，十分適合GPU并行計(jì)算誉裆。

統(tǒng)計(jì)學(xué)中的一個(gè)概念叫做Distribution Function（分布函數(shù)） $f(x, \xi, t)$ ，這個(gè)函數(shù)的三個(gè)參數(shù)分別對(duì)應(yīng)位置缸濒、速度與時(shí)間足丢，表示的是隨著三個(gè)參數(shù)的不同，某種現(xiàn)象發(fā)生的概率庇配，這里我們考慮的是流體中微觀層面（如布朗運(yùn)動(dòng)）現(xiàn)象的概率斩跌，比如在任何時(shí)間t在位置x處觀察到的以速度 $\xi$ 移動(dòng)的空氣分子的概率（或者密度）。

對(duì)這個(gè)函數(shù)沿著速度（三維）進(jìn)行積分捞慌，就能得到當(dāng)前時(shí)刻對(duì)應(yīng)位置處的空氣分子的密度（這是一個(gè)標(biāo)量）：

$\int f(x, \xi, t) d^3 \xi = \rho(x,t)$

而這個(gè)公式可以稱之為0階動(dòng)量（0-th momentum)耀鸦，而一階動(dòng)量則可以表示為：

$\int \xi f(x, \xi, t) d^3 \xi = \rho u$

這里的u是宏觀層面的速度（即NS方程中的速度），其實(shí)左邊項(xiàng)可以看成是微觀層面的速度乘以微觀的速度（密度）啸澡，積分之后就是宏觀層面的動(dòng)量了袖订，這個(gè)結(jié)果是一個(gè)向量。

從上面兩個(gè)公式也可以窺見锻霎，LBM方法的要點(diǎn)就是拋開宏觀層面的計(jì)算著角，我們關(guān)心微觀層面的概率分布函數(shù)是什么樣的揪漩。

我們先來看下旋恼，概率分布函數(shù)的微分形式應(yīng)該是什么樣的：

$\frac{df}{dt} = \frac{\partial f}{\partial t} + \nabla_x f \cdot \frac{\partial x}{\partial t} + \nabla_{\xi} f \cdot \frac{\partial \xi}{\partial t} = \frac{\partial f}{\partial t} + \nabla_x f \cdot \xi + \nabla_{\xi} f \cdot F = \Omega(f)$

上面式子中的F是外力（其實(shí)是加速度，這里假設(shè)質(zhì)量為1）奄容， $\Omega$ 叫做Collision Operator冰更，如果這個(gè)函數(shù)是知道的，那么后面這個(gè)等式被叫做Boltzman Equation昂勒，不過歷史上大家是不知道這個(gè)公式是什么樣的蜀细，因此大家會(huì)通過一些條件來推測(cè)一些經(jīng)驗(yàn)公式。

質(zhì)量守恒
$\int \Omega(f) d^3 \xi = 0$
動(dòng)量守恒
$\int \xi \Omega(f) d^3 \xi = 0$
能量守恒
$\int |\xi|^2 \Omega(f) d^3 \xi = 0$

根據(jù)這三個(gè)條件戈盈，大家湊出一個(gè)模型奠衔，叫做BGK Model：
$\Omega(f) = -\frac{1}{\tau}(f - f^{eq}) = \frac{df}{dt}$

這個(gè)公式中 $\tau$ 是時(shí)間的常量，而 $f^{eq}$ 則是將氣體在靜止很長(zhǎng)一段時(shí)間之后塘娶，溫度達(dá)到平衡之后归斤，對(duì)應(yīng)的distribution function。

右邊的等式就是一個(gè)常微分方程刁岸，手算之后其結(jié)果輸出為：
$f(x, \xi, t) = f^{eq}+[f(x, \xi, 0)-f^{eq}]e^{-\frac{t}{\tau}}$

這個(gè)公式的意思是脏里，隨著時(shí)間的推移，概率將會(huì)以指數(shù)的形式遞減虹曙，最終達(dá)到 $f^{eq}$ 迫横，流體密度越大番舆， $\tau$ 越小聂儒，越快達(dá)到平衡患蹂，也就是說响迂， $\tau$ 控制的是流體的粘性削锰。

BGK方法雖然簡(jiǎn)單护侮，但是有一個(gè)比較顯著的問題在于準(zhǔn)確性與真實(shí)情況相去甚遠(yuǎn)拳氢，畢竟是拼湊出來的公式践剂，后面大家觀察到弓摘，這個(gè)不準(zhǔn)確性來自于圣蝎，我們這里將f看成是一個(gè)隨著時(shí)間均勻趨向平衡的函數(shù)刃宵，但是實(shí)際上f代表的各個(gè)特性，如質(zhì)量徘公、動(dòng)量以及能量（也就是不同階數(shù)的梯度）其趨向平衡的速度是不同的牲证。

基于上述發(fā)現(xiàn)，2006年之后出現(xiàn)了一系列新的方法大大提升了整個(gè)方案的精確性（當(dāng)然也會(huì)變得更為復(fù)雜）关面，時(shí)間關(guān)系坦袍，這里就不做擴(kuò)展。

不管怎么說等太，我們都有了一個(gè)Collision Operator捂齐，接下來要做的事情就是根據(jù)timestep對(duì)其進(jìn)行積分以求得f。

先來看一個(gè)叫做Hermite Interpolation的方法缩抡，以2D Domain為例奠宜，我們將流體場(chǎng)域分成Grid，每個(gè)頂點(diǎn)上對(duì)應(yīng)的就是位置x瞻想，對(duì)于每個(gè)頂點(diǎn)而言压真，我們存儲(chǔ)的不是一個(gè)數(shù)值，而是多個(gè)數(shù)值蘑险，每個(gè)數(shù)值對(duì)應(yīng)的f在同一時(shí)刻下不同方向下的數(shù)值 $f(x, \xi_i, t)$ （即不同速度的概率滴肿？），這些方向分別對(duì)應(yīng)于當(dāng)前頂點(diǎn)指向相鄰頂點(diǎn)（加上自身）佃迄，總共9個(gè)方向泼差，因此總結(jié)為D2Q9，在3D情況下呵俏，就是D3Q27等等堆缘。

下面我們看下，如何根據(jù)Collision Operator求得f柴信。

$\frac{\partial f}{\partial t} + \nabla_x f \cdot \xi + \nabla_{\xi} f \cdot F = \Omega(f)$

同樣套啤，這里采用Splitting Operator。

先考慮前兩項(xiàng)：

$\frac{\partial f}{\partial t} + \nabla_x f \cdot \xi = 0$

這個(gè)就是advection equation，在NS方程求解的時(shí)候潜沦，advection是速度場(chǎng)萄涯，當(dāng)時(shí)是使用Semi-Lagrangian方法求解的，而這里更為簡(jiǎn)單唆鸡，考慮剛才的格子D2Q9涝影，如果我們只考慮一個(gè)方向的分量，按照Semi-Lagrangian方法争占，由于頂點(diǎn)上存儲(chǔ)的指向各個(gè)相鄰頂點(diǎn)的速度其長(zhǎng)度正好在一個(gè)timestep上走到了相鄰頂點(diǎn)上燃逻，因此對(duì)于相鄰頂點(diǎn)在當(dāng)前timestep的物理量就根本不用求解，只需要將此頂點(diǎn)上的物理量拷貝過去即可臂痕，從而就能完成這個(gè)方程的求解（文獻(xiàn)中稱之為Streaming）伯襟。

最后考慮剩下的部分，在剛才計(jì)算的結(jié)果上加上Collision Operator（這個(gè)是已知的握童，前面說過的經(jīng)驗(yàn)?zāi)Ｐ停┮约巴饬Φ淖饔媚饭郑偷玫搅讼旅娴氖阶?/p>

$\frac{\partial f}{\partial t} = \Omega(f) - \nabla_{\xi} f \cdot F$

這個(gè)方法由于在每個(gè)頂點(diǎn)上進(jìn)行的計(jì)算都是完全獨(dú)立的，因此適合在GPU上計(jì)算澡绩，適合用在實(shí)時(shí)場(chǎng)景稽揭。

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者

人面猴
序言：七十年代末，一起剝皮案震驚了整個(gè)濱河市肥卡，隨后出現(xiàn)的幾起案子溪掀，更是在濱河造成了極大的恐慌，老刑警劉巖步鉴，帶你破解...
沈念sama閱讀 212,816評(píng)論 6贊 492
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件揪胃，死亡現(xiàn)場(chǎng)離奇詭異，居然都是意外死亡唠叛，警方通過查閱死者的電腦和手機(jī)只嚣，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 90,729評(píng)論 3贊 385
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進(jìn)店門沮稚，熙熙樓的掌柜王于貴愁眉苦臉地迎上來艺沼，“玉大人，你說我怎么就攤上這事蕴掏≌习悖” “怎么了？”我有些...
開封第一講書人閱讀 158,300評(píng)論 0贊 348
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵盛杰，是天一觀的道長(zhǎng)挽荡。經(jīng)常有香客問我，道長(zhǎng)即供，這世上最難降的妖魔是什么定拟？我笑而不...
開封第一講書人閱讀 56,780評(píng)論 1贊 285
?港島之戀（遺憾婚禮）
正文為了忘掉前任，我火速辦了婚禮逗嫡，結(jié)果婚禮上青自，老公的妹妹穿的比我還像新娘株依。我一直安慰自己，他們只是感情好延窜，可當(dāng)我...
茶點(diǎn)故事閱讀 65,890評(píng)論 6贊 385
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布恋腕。她就那樣靜靜地躺著，像睡著了一般逆瑞。火紅的嫁衣襯著肌膚如雪荠藤。梳的紋絲不亂的頭發(fā)上，一...
開封第一講書人閱讀 50,084評(píng)論 1贊 291
城市分裂傳說
那天获高，我揣著相機(jī)與錄音哈肖，去河邊找鬼。笑死念秧，一個(gè)胖子當(dāng)著我的面吹牛牡彻，可吹牛的內(nèi)容都是我干的。我是一名探鬼主播出爹，決...
沈念sama閱讀 39,151評(píng)論 3贊 410
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼庄吼，長(zhǎng)吁一口氣：“原來是場(chǎng)噩夢(mèng)啊……” “哼！你這毒婦竟也來了严就？” 一聲冷哼從身側(cè)響起总寻，我...
開封第一講書人閱讀 37,912評(píng)論 0贊 268
萬榮殺人案實(shí)錄
序言：老撾萬榮一對(duì)情侶失蹤，失蹤者是張志新（化名）和其女友劉穎梢为，沒想到半個(gè)月后渐行，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體，經(jīng)...
沈念sama閱讀 44,355評(píng)論 1贊 303
?護(hù)林員之死
正文獨(dú)居荒郊野嶺守林人離奇死亡铸董，尸身上長(zhǎng)有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點(diǎn)故事閱讀 36,666評(píng)論 2贊 327
?白月光啟示錄
正文我和宋清朗相戀三年祟印，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了。大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片粟害。...
茶點(diǎn)故事閱讀 38,809評(píng)論 1贊 341
活死人
序言：一個(gè)原本活蹦亂跳的男人離奇死亡蕴忆，死狀恐怖，靈堂內(nèi)的尸體忽然破棺而出悲幅，到底是詐尸還是另有隱情套鹅，我是刑警寧澤，帶...
沈念sama閱讀 34,504評(píng)論 4贊 334
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布汰具，位于F島的核電站卓鹿，受9級(jí)特大地震影響，放射性物質(zhì)發(fā)生泄漏留荔。R本人自食惡果不足惜吟孙，卻給世界環(huán)境...
茶點(diǎn)故事閱讀 40,150評(píng)論 3贊 317
男人毒藥：我在死后第九天來索命
文/蒙蒙一、第九天我趴在偏房一處隱蔽的房頂上張望。院中可真熱鬧杰妓，春花似錦肥隆、人聲如沸。這莊子的主人今日做“春日...
開封第一講書人閱讀 30,882評(píng)論 0贊 21
一樁弒父案栋艳，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽。三九已至句各，卻和暖如春吸占，著一層夾襖步出監(jiān)牢的瞬間，已是汗流浹背凿宾。一陣腳步聲響...
開封第一講書人閱讀 32,121評(píng)論 1贊 267
情欲美人皮
我被黑心中介騙來泰國(guó)打工矾屯，沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留，地道東北人初厚。一個(gè)月前我還...
沈念sama閱讀 46,628評(píng)論 2贊 362
代替公主和親
正文我出身青樓件蚕，卻偏偏與公主長(zhǎng)得像，于是被迫代替她去往敵國(guó)和親产禾。傳聞我的和親對(duì)象是個(gè)殘疾皇子排作，可洞房花燭夜當(dāng)晚...
茶點(diǎn)故事閱讀 43,724評(píng)論 2贊 351

第六講：流體模擬- 流體力學(xué)基礎(chǔ)

SPH（Smoothed Particle Hydrodynamics）

Lattice Boltzman Method

推薦閱讀更多精彩內(nèi)容