【USparkle專欄】如果你深懷絕技，愛“搞點研究”，樂于分享也博采眾長，我們期待你的加入，讓智慧的火花碰撞交織，讓知識的傳遞生生不息！

這是侑虎科技第1787篇文章，感謝作者楊超wantnon供稿。歡迎轉發(fā)分享，未經(jīng)作者授權請勿轉載。如果您有任何獨到的見解或者發(fā)現(xiàn)也歡迎聯(lián)系我們，一起探討。（QQ群：793972859）

作者主頁：

https://www.zhihu.com/people/wantnon

一、原理

參考這個論文：

《Real-time Simulation of Large Bodies of Water with Small Scale Details》

https://matthias-research.github.io/pages/publications/hfFluid.pdf

核心是這兩個公式：

我在這篇《波動方程與淺水方程》中作了說明，轉貼如下：

淺水方程：

（1）動量方程

（2） 質(zhì)量守恒

其實淺水方程的動量方程，就是在NS方程的動量方程

中令p=ρgη（η為海拔）得來，物理含義是表面隆起（或凹陷）產(chǎn)生壓強（或負壓強）。

注意η和h的區(qū)別，h是水深，η是海拔，如果設河床高度為H，則有η=H+h。H是x，y的函數(shù)，h是x，y，t的函數(shù)。

對于普通水可以忽略粘性，進一步簡化為：

論文中所用公式為：

這其實就是淺水方程。說明如下：

D為物質(zhì)導數(shù)，定義為

將（2）式中物質(zhì)導數(shù)展開得

可見就是淺水方程的動量方程忽略掉粘性項。再看（1）式，將其中物質(zhì)導數(shù)展開得

又因為h是標量v是向量且均為空間的函數(shù)，根據(jù)散度的運算性質(zhì)，有

代入得

可見就是淺水方程的質(zhì)量守恒方程。

注：

1. 淺水方程的動量方程（2）用的是海拔η ，質(zhì)量守恒方程（1）用的是水深h，這正是它可以適用于凹凸不平地形，模擬洪泛的原因。假如我們強行讓兩方程都用h，那就相當于是假定海拔等于水深，即水底是平面，則不再具備模擬洪泛的能力，但仍可模擬水波。

質(zhì)量守恒方程：

求解，對應論文中：

動量方程：

求解，其中 求解對流項：

對應論文中：

由于沒給公式，我暫時忽略了這項。

求解壓力和外力項：

對應論文中：

二、MAC網(wǎng)格

論文中的離散化形式，含有：

這種一半的項，是因為引入了稱作MAC網(wǎng)格的技巧，即將流入/流出速度看作是在cell的邊上，在處理散度相關問題時，經(jīng)常用到。

下圖很直觀：

Eulerian Fluid Simulation Notes [Part 1]

https://quangduong.me/notes/eulerian_fluid_sim_p1/

Bridson's notes：nov17-cs533d-slides.ppt

https://www.cs.ubc.ca/~rbridson/courses/533d-fall-2005/nov17-cs533d-slides.pdf

三、實現(xiàn)

最終是要用Shader實現(xiàn)，但初期驗證算法出于Debug方便考慮，還是先用Houdini實現(xiàn)。

當然，我們可以真的在Houdini里建一個非?！熬呦蟮摹盡AC網(wǎng)格，就是每個格子表示一個cell，cell的邊界上可以附加速度向量。

但為了將來轉Shader更直接，我用的更接近二維數(shù)組的結構，我用一個二維點陣代表Grid，每個Point代表一個cell，Point有如下屬性：

f@h=0; f@u_r=0;//即u[i+1/2,j] f@w_d=0;//即w[i,j+1/2] f@h_r=0;//即h[i+1/2,j] f@h_d=0;//即h[i,j+1/2]

至于為啥沒有u_l（u[i-1/2,j]），w_u（w[i,j-1/2]），h_l（h[i-1/2,j]），h_u（h[i,j-1/2]），是因為當前Point的u_l其實就是左邊格子的u_r，即它已經(jīng)在左邊格子中存過的，就沒有必要在此格中再存一遍了。如果我們想獲得u_l，只需如下：

int ID=@ptnum;

int i=floor(ID/N);

int j=ID%N;

int ID_L=j*N+i-1;

float u_l=point(0,"u_r",ID_L)

w_u，h_l，h_u同理。

這也正是論文中這段：

只給出

表達式，而沒有寫

表達式的原因。因為當前cell的

正是左邊cell和上邊cell（假設我的Grid是y軸向下）的

只要左邊cell和上邊cell已經(jīng)計算過，就可以直接訪問而無需在本cell中再重新計算一次。

同理，下面這段只給

也是同樣原因：

然后，就只是抄公式了。

四、穩(wěn)定性

由于是離散模擬，而且用的是最簡單的顯式積分，所以會有數(shù)值穩(wěn)定性問題，論文中列了幾個措施可以提高穩(wěn)定性：

Overshooting Reduction：

在論文的2.2節(jié)，提到有水與沒水交界處，h的振幅會因數(shù)值問題而被放大，導致邊緣出現(xiàn)尖刺甚至撕裂，用如下方法緩解：

Stability Enhancements：

在論文的2.1.5節(jié)，提到如下約束可提高穩(wěn)定性：

五、參數(shù)對效果的影響

g，α，Δx，Δt對效果均有影響。

例如，不同α值對比：

α=0.5

α=0.6

可見，α =0.5和α=0.6相比，由于速度限制更狠，顯得更粘稠。

以上我們實現(xiàn)了凹凸不平地形上的淺水方程模擬，但忽略了對流項。

以下作了些改進，主要添加了對流項、渦度約束和邊界條件。

效果：

六、添加對流項

對流項：

有兩種實現(xiàn)方式：一種是Gather方法，一種是直接離散化方法。

（一）Gather方法

實際上在Gather方法之前，還有一個Scatter方法，這兩種方法都是拋開對流的數(shù)學公式，而直接從公式背后的物理含義出發(fā)的模擬方法。Scatter方法就是把當前cell的屬性傳輸出去，Gather方法就是把上游點的屬性傳輸給當前cell。在《NS方程與2d流體模擬》中說過，Gather方法由于對GPU友好，更常用。

對流的Gather方法實現(xiàn)，可以參考《游戲中的流體模擬（Fluid Simulation），技術美術教程 》。

由于我要在之前的基礎上加對流，而那個模擬過程中只使用到邊速度，但上面教程中的實現(xiàn)用的是cell速度，所以需要在邊速度和cell速度之間互轉。

過程如下：

Pass1：由邊速度計算cell速度。

Pass2：計算對流Advect。

1. 用9個格加權平均估算當前cell平滑速度。

2. 將當前cell平滑速度乘以dt再反向，以此作為偏移量去采樣上游目標點的速度（用雙線性插值）。

3. 用采到的速度替換當前cell的速度。

Pass3：由cell速度計算邊速度。

其中需要注意的點有：

1. 邊速度與cell速度之間的互轉。

我問的ChatGPT：

由于兩個問題我是連續(xù)問的，所以它給出的公式也恰好是互逆的。我驗證了一下，如果把P ass2屏蔽，只剩Pass1和Pass3，結果確實不受影響。

2. 無水處的速度

用9個cell的速度加權平均計算當前cell平滑速度時，不用管各cell是否有水，照加不誤。采樣上游點速度覆蓋當前cell速度時，也不用管上游點處是否有水，照覆蓋不誤。原因是，只要能保證無水處速度是0，上面的計算就都是合理的。至于如何保證無水處速度為0，見后面"邊界條件"一節(jié)。

3. 數(shù)值穩(wěn)定性及對流快慢和強度控制

以下兩項措施可以提高對流的數(shù)值穩(wěn)定性：

（1）縮短步長：將偏移-V*dt采樣改為偏移-V*dt*0.5采樣，這樣雖然會導致對流變慢（如果滴入墨水的話，體現(xiàn)為墨水擴散變慢），但有助于數(shù)值穩(wěn)定性。如果將步長減半的同時又不想對流變慢，可以在一幀內(nèi)將對流迭代兩次。

（2）減弱影響：將“用目標點速度覆蓋當前點速度”這一步，改為不完全覆蓋，而是按一定比例（如0.5）進行混合。這樣雖然會使對流效果減弱，但也有助于數(shù)值穩(wěn)定性。

是否縮短步長和減弱影響，也未必總是出于數(shù)值穩(wěn)定性考慮，有時候為了效果，也需要對對流的快速和強度進行控制。

4. 邊框處理

計算當前cell平滑速度時，9個cell中超出邊框者要去掉。采樣上游點時，如果上游點超出邊框，則放棄當前cell的對流計算。

（二）直接離散化方法

就是將

展開為（設v=(u,w)）：

其中右邊各項離散化如下：

這里用的是 迎風格式，參考ChatGPT：

同樣可以用縮短步長和減弱影響對對流快慢和強度進行控制，以及提高數(shù)值穩(wěn)定性。

七、添加墨水

要想讓對流效果明顯地顯示出來，需要添加墨水。

方法如下：

1. 標記一些cell為inkSource。

2. 每幀對標為inkSource的cell追加ink值：ink+=inkAmount*dt

3. 在解算對流時：

如果用的是Gather方法，則除了用上游點的速度覆蓋當前點速度外，還要用上游點ink覆蓋當前點ink。如果用的是直接離散化方法，則除了計算du、dw，(u,w)+=(du,dw)，還要計算dink，ink+=dink。

4. 最后將ink可視化。

八、兩種對流方法效果對比

下面視頻中左邊為Gather，右邊為直接離散化：

可見，大致效果接近，由于兩種方法很不同，所以基本可佐證實現(xiàn)是正確的。

九、添加渦度約束

方法在《NS方程與2d流體模擬》中已經(jīng)寫過。

需要注意的一點是，如果GetCurl函數(shù)中沒有除以Δx，則要在外面補上。另外為了讓效果與時間解耦，最終速度調(diào)整量還應乘以Δt。

有無渦度約束對比：

左： 無渦 度約束，右： 有渦度約束

可見，區(qū)別還是很明顯的。

十、添加邊界條件

1. 邊框處理

在之前的文章中，我沒有限制水不能流出邊框，所以水一直沿著河道流，永遠不會漫上來。

這次，我加了水不能流出邊框的限制，只需將邊框視作障礙物，每幀末尾將邊框cell的邊速度和cell速度均標為0即可。這在論文《hfFluid.pdf》中有寫：

2. 使無水處速度為0

淺水方程本身是無法保證無水處速度為0的，因為從

來看，當無水，即h=0時，變?yōu)?/p>

即仍會 由于地形本身有梯度而產(chǎn)生速度。 所以，要想讓無水處速度為0，只能通過人為添加約束實現(xiàn)。

其實，上面Boundary Conditions這段已經(jīng)給出了方法，就是通過

判斷右側 邊（i+1/2,j）是否為“岸”。 類似方法也可以判斷出左、上、下各邊是否為岸。

而且分析一下會發(fā)現(xiàn)，這種判斷方法下，不僅岸與非岸的交界邊會被判為岸，岸上的邊也都會被判為岸，這正是我們想要的。于是對于判斷為岸的邊，我們在每幀末尾將其邊速度強制設為0，則可實現(xiàn)無水處速度為0的目標。至于cell速度，由于它是由邊速度算出的，所以只要邊速度標對了，cell速度自然也就對了。

另外，由于我們是在每幀末尾將無水處速度標為0的，所以我將對流的計算放在最前面，以保證對流計算是在無水cell速度為0的情況下進行的。

下圖是不加岸邊速度置0（左）與加岸邊速度置0（右），速度場及最終效果對比：

左： 不 加岸邊速度置0，右： 加岸邊速度置0

看起來似乎差別不大，但有一點，如果近看的話，會發(fā)現(xiàn)岸邊速度不置0的模擬結果，水會在岸邊隆起，不太合理，而加了岸邊速度置0的模擬結果，水與岸邊則是比較平滑的過度，如圖所示：

左： 不加岸 邊速度置0，右： 加岸邊速度置0

參考文章：

《波動方程與淺水方程》

https://zhuanlan.zhihu.com/p/23947310903

《NS方程與2d流體模擬》

https://zhuanlan.zhihu.com/p/23421572373

《游戲中的流體模擬（Fluid Simulation），技術美術教程》

https://zhuanlan.zhihu.com/p/627468747

《hfFluid.pdf 》

https://matthias-research.github.io/pages/publications/hfFluid.pdf

文末，再次感謝 楊超wantnon 的分享， 作者主頁： https://www.zhihu.com/people/wantnon， 如果您有任何獨到的見解或者發(fā)現(xiàn)也歡迎聯(lián)系我們，一起探討。（QQ群： 793972859 ）。

近期精彩回顧