引
“咚咚咚。”
“門外是誰？”
“原來白羊在家啊！叫了半天也沒人答應，我還以為白羊不在家呢。”
忽然想到什麼，齊夏渾身一震。一陣細小的滴答聲在耳邊陡然響起，整個人汗毛豎立。眼前景象如雨滴落入水面般蕩起漣漪，變得支離破碎。
……
“門外不是我”
序
"There are two ways of doing calculations in theoretical physics.
One way, and this is the way I prefer, is to have a clear physical picture of the process you are calculating.
The other way is to have a precise and self-consistent mathematical formalism.
...
You have neither."
——Enrico Fermi
1 概述
水體模擬一直是計算機圖形學的熱點研究領域之一，也是遊戲開發、影視製作等領域必不可少的技術手段。如果VR技術進一步發展想要實現元宇宙的夢想，如何在虛擬世界呈現出像現實世界一般逼真的流體效果更是一個繞不過去的話題。從半衰期(Half-Life, Valve 1998)的純動畫紋理貼圖，到俠盜獵車手5(Grand Theft Auto V, Rockstar 2013)的紋理動畫加粒子特效，到死亡擱淺(Death Stranding, Kojima Project 2019)的弱交互水面，再到最後生還者：第一部(The Last of Us Part I, Naughty Dog 2022)的強物理交互水體，可以看到水體模擬技術不斷進步的一個演化過程。
流體力學研究現實中流體作用和運動規律，計算流體力學希望準確快速地再現實際物理流動。看似與計算機圖形學的水體模擬目的相同，但計算流體力學並不涉及圖像的渲染和呈現過程，而水體模擬也不需要完全遵循物理規律：對它而言更重要的是讓觀察者“認為”水足夠擬真。二者的交集為流體力學的基礎研究方法，即對物理對象或過程進行建模並用數學語言描述。由於不懂計算機圖形學，以下主要內容均為交集部分。
淺水方程(Shallow Water Equation)是從流體力學基本方程（組），納維-斯托克斯方程(Navier-Stokes Equations，下稱NS方程)基於若干假設簡化得出。聖維南伯爵於1871年提出的聖維南方程就是淺水方程的1維形式，而淺水方程至今仍有重要地位，被廣泛應用在海洋河流流動、潮汐、潰壩等眾多領域。20世紀80年代左右，人們開始研究使用計算機對自然現象進行建模渲染。到90年代前後，將流體視為大量粒子系統並追蹤粒子的拉格朗日法，和將流體視為空間流場的歐拉法相繼被引入水體模擬渲染。1990年，Kass與Miller首次提出用淺水方程來模擬水面，在當時獲得了不錯的效果。
本文主要介紹淺水方程的物理過程（雖然以上舉例的遊戲中未必在使用淺水方程），而不涉及計算機圖形學水體模擬部分的光照計算、光柵化等渲染過程。先介紹流體靜力平衡假設，再引入邊界條件，最後分別用微分法和積分法推導2維淺水方程。此前，先給出NS方程的一般形式用於後續計算，眾多流體力學、計算流體力學參考書給出了NS方程本身的具體推導過程，不再贅述。
由質量守恆得出的連續性方程：
由牛頓第二定律得出的三方向動量方程：
其中考慮了地球自轉產生的科氏力的影響。未用矢量或張量表示便於後續對非定常項、對流項等進行分析和變形。動量方程中的科氏力項為
其中θ為當地緯度，Ω為星球自轉角速度。當考慮與星球半徑一個數量級的大尺度洋流運動或風暴、氣旋運動等，科氏力項通常不可忽略。本文中不考慮這種極大尺度。如質量效應2(Mass Effect 2, BioWare 2010)中在不同星球上採礦時，如果需要模擬海洋或氣旋流動（當然可能直接做個動畫貼圖更省事還節省性能開銷）。看評價褒貶不一還沒入星空(Starfield, Bethesda 2023)，不知道里面是否需要模擬行星級別的流動。
粘性力中，切應力項由廣義牛頓內摩擦定律給出：
正應力項由斯托克斯假設給出（ii非啞標）
2 基本假設
2.1 流體密度不變
與氣體相比，固體和液體的可壓縮性要小得多，在溫度沒有劇烈變化的情況下通常可視為恆密度。對于海洋等研究對象，其含鹽量、溫度、含沙量等會對密度造成一定影響，但暫時不考慮，也認為具有恆定的密度。在2.6.3節中有布辛尼斯克(Boussinesq)假設供修正。在密度為常數的情況下，(1.1)~(1.4)式的密度可以提到括號外。(1.1)式退化成
2.2 切平面近似
在流體尺度與地球尺度相比較小的情況下，流體在地球表面上的運動可看作是在與地球接觸處的切平面上的運動。此時對切平面上物體的運動而言，起主導作用的是垂直於運動平面方向的Ωz，而忽略Ωy。
2.3 淺水近似
淺水方程要求流體深度(z)方向尺度遠小於長度和寬度(x-y)方向尺度，具體原因見2.7節。
2.4 平面高度函數假設
認為水面與水底的高度可表示為x y t的函數。見(3.1)和(3.2)式。
2.5 拉格朗日假設
拉格朗日假設流體微元一旦進入水面或水底，將永遠保持在該曲面上。
2.6 流體靜力平衡近似
流體靜力平衡近似認為流體在豎直方向的加速度可以忽略不計。從最簡單的角度看，在(1.4)式中忽略科氏力、粘性力，等式左邊為z方向的加速度Dw/Dt，右邊只剩壓力項和浮力項。當豎直方向的加速度與重力相比很小時，Dw/Dt≈0，所以在靜力平衡狀態下流體的壓力梯度與浮力項平衡。定量分析如下：
2.6.1 量綱分析
考慮z方向動量方程(1.4)式的各項相對大小。有如下特徵尺度：x、y方向的長度尺度L，速度尺度U；z方向深度尺度H<<L，速度未知；地球自傳角速度尺度Ω，則根據(2.1)式可知w方向速度尺度為
對於(1.4)式中各項，統一除以密度再分析有
如果根據特徵速度和特徵長度構造特徵時間尺度，則
將(2.2)和(2.4)帶入(2.3)，得到非定常項和對流項相同，正應力和切應力相同。U~L二者一個量級，W很小，非定常項和對流項為小量。地球自轉角速度很慢而U不大，科氏力項也非常小。粘性項中由於運動粘度ν在10-3量級遠小於H，故正應力項也是小量。所以，只有壓力梯度項能夠跟浮力項平衡，其他項都可忽略。
如果引入3個無量綱數，分別為慣性力與粘性力之比、與科氏力之比、與壓力梯度項之比：
其中f=2Ωsinθ，再用淺水波內的“聲速”，即特徵信號傳遞時間（見2.6.4節）當作時間尺度，
對(2.3)各項再除以g，可化為
其中Re數在100~1000量級。Ro數在河流、海岸附近為1左右，研究地球大尺度洋流時Ro在10-3量級（此時不可忽略，文中不考慮）。F數很難超過0.2。各值帶入(2.9)也能發現只有壓力梯度項能與浮力項相當。
綜上，在流體靜力平衡的假設下，z方向動量方程(1.4)退化為
式(2.10)也稱靜水壓平衡方程或歐拉靜平衡方程。如果將水壓分為兩部分
其中p'為修正壓力，而ph有
帶入無粘無科氏力的動量方程得
因此對於常密度流體而言，重力加速度等效於對壓力分佈進行修正，使原動量方程的浮力項消失，而重力本身不會產生動力學效應。
如果認為水面壓力ps為常量並取為參考點，ps=0，Zs為水面高度，是x、y、t的函數。積分(2.10)得
2.6.2 水平速度與深度無關
式(2.14)對x求偏導得
可見(2.15)式與z無關。帶入x方向動量方程(1.2)，可見x方向加速度與z無關（z方向粘性項通過量綱分析可忽略）。y方向同理。因此u、v與z無關。
2.6.3 布辛尼斯克假設
流體靜力近似在很多情況下較為準確，但誤差隨海水密度變化而增加，尤其是海洋含鹽量等有一定變化時。由於整體流體的密度變化量一般要小於百分之幾的量級，這樣的改變量對於慣性項和粘性項而言幾乎沒有影響，但對浮力項有較大影響。布辛尼斯克(Boussinesq)提出，在動量方程中，可以認為除了浮力項之外的密度為常數，而浮力項本身的密度會受到其他因素影響。此時密度視為一個常量加變化量
其中ρ0為常量，而變化量δρ<<ρ0。帶入z方向動量方程(1.4)
根據假設，左側密度仍視為ρ0，而右側不變，再忽略科氏力和粘性項，得到
對p也進行拆分，
其中p0對應靜力平衡、密度為ρ0時的壓力梯度，δp為密度波動導致的壓力波動
(2.19)和(2.20)帶入(2.18)得
其中
在該假設下，x方向動量方程為
y方向同理。前有δρ<<ρ0，而式(2.22)中不忽略b的原因是g相對來說比較大（從量綱分析中也可看出），因此b未必很小，不可忽略。而b就是在布辛尼斯克假設下的浮力項，修正了純靜力平衡時浮力項缺失的問題。
浮力項密度改變的情況下，對(2.10)積分，水面壓強仍認為相等且為0，得
帶入(2.16)得
流體密度的變化量可用有效示蹤劑的濃度來表示，如含鹽量、含沙量等，溫度也可包含在內，有
其中β為膨脹係數，為恆壓下密度變化率與原密度之比
對(2.27)線性近似得
有多個量同時變化時累加即為(2.26)。
後續不考慮在豎直方向上流體密度的變化，因為一旦這樣，會在靜力平衡方程中引入z，使得x、y方向速度與z有關。因此認為流體在豎直方向沒有密度變化，而允許在水平方向有密度變化，且仍可用膨脹係數來描述。
2.6.4 氣體等效
流體密度不變時，引入單位面積密度和單位長度壓力
Zs(x, y, t)為海面高度，Zb(x, y)為海底高度，不考慮海底移動的情況。(2.29)對時間的導數為
帶入(2.14)得
(2.32)滿足絕熱關係，且絕熱指數γ=2。根據氣體動力學聲速公式可得
h為水深，h=Zs-Zb。
式(2.33)也可從淺水方程得到的波動方程線性化後求解特徵值得出。
2.7 再看淺水近似
在(2.23)中由量綱分析有（T取L/U）
帶入(2.21)得
因此若要忽略z方向的加速度，滿足流體靜力平衡假設，就需要H<<L，這也是叫淺水方程的一個原因。另外z方向密度不變也是後續積分的前提，淺水近似也能使該假設更合理。
2.7.1 z方向密度不變
為能得到深度平均的2維淺水方程，密度在z方向的變化會破壞深度平均的簡潔。因此不考慮布辛尼斯克對於豎直z方向密度的變化相關假設，而允許密度在水平xy方向有改變。此時對於式(2.10)積分得
其中ps為水面大氣壓。對上式求x方向偏導得
y方向同理。可見其中出現了密度在x方向的導數項，如有必要可用(2.26)式帶入膨脹係數具體計算。
2.7.2 z方向加速度和速度很小
由於靜力平衡假設，流體在z方向認為無加速度，但這並不意味流體在z方向無速度，而是z方向速度變化很慢；更嚴謹應該是z方向流體的加速度與g相比非常小，即使認為z方向加速度為0，也能滿足靜力平衡近似得到的結果。
另外還有淺水近似、平面高度函數假設和拉格朗日假設，流體在z方向的速度必須非常小才能不破壞表面。而u、v與z無關，水深h又很小，通常可認為u、v在z方向是均勻的，但不均勻也可以，後續具體討論。
2.8 3維淺水方程
綜上，3維密度xy可變下靜力平衡切平面近似的x方向動量方程為
其中ps為水面大氣壓力，這裡認為其在xy方向也有變化。密度由前述也允許在xy方向有變化。y方向動量方程同理得
z方向動量方程已退化進xy方向。連續方程仍為(2.1)。
3 邊界條件
3.1 數學
淺水方程在數學上要求，水面ζs和水底ζb光滑且沒有洞，它們的運動是一個與時間有關的連續拓撲變形。還認為，水面和水底的高度能夠看作x y t的函數，即
其中Zs與Zb為水面與水底的高度函數。由(3.1)可得
即為水面的邊界條件。
若流動無旋，有速度勢函數ψ滿足
在水底速度勢函數在底面法向方向的偏導應為0，即
因此水底的邊界條件為
若水底可能移動，則水底邊界條件與水面相同，將(3.3)式中的Zs換Zb即可。
3.2 物理
拉格朗日假設即對應3.1節數學上對水面與水底的要求。
從物理約束上看，對於不可穿透的水底，應該有流體垂直於底面方向的速度分量為0，即
其中V為流體速度矢量，而空間曲面ζ的法向量為
帶入(3.7)即可得到(3.6)的結果。對於水底還有無滑移條件
對於水面，為保證平面高度函數假設，需要使水面流體微元相對於水面的法向速度為0，即為(3.3)式。
水面的物理約束不夠清晰，提供兩種不夠嚴謹（甚至錯誤）的“民科”思路：
3.2.1 時空速度？
空間中運動的曲面ζ(x, y, z, t)上面一點的座標為
該點在時空中的速度
運動曲面的時空法向量
水面的相對不可穿透性有
得到(3.3)的約束式。
3.2.2 表面守恆？
如果將表面流體微元的高度h=(Zs-Zb)進行積分，得到的是整個流場內流體的體積。又因為不可壓縮性，質量守恆就是體積保持不變，因此該積分為定值，有
但是因為積分區域Zs在隨時間變化，導數不能直接與積分交換順序。考慮三維萊布尼茨積分法則：
其中Vs為運動表面的速度矢量由於表面並不是整個在空間平動，而是各部分有各自的速度還帶彎曲，無法確定Vs。而且(3.15)式內被積分量為矢量，但h是標量，也許可以指定h方向豎直向上構造出h矢量，但後面也不知道能不能算。
感覺這種方法應該是可以，對應於控制體的積分法，只是數學不好算不來。
3.2.3 水面剪切外力？
如果有風，對錶面還有剪切力約束
其中τsx為風產生的剪切力，可認為已知。（圖中a為水深h，h為表面高度Zs）
沒有發現(3.16)是如何得出的。能看出右側為微元體x平面的切應力矢量，點乘表面法向量，但是為什麼它就等於風吹產生的剪切外力？
4 2維淺水方程：微分法
4.1 連續性方程
取水中任一微元控制體分析：
A為底面面積，S為側面面積，n為對應面的法向量。流入為負流出為正。
4.1.1 uv在z方向均勻
單位時間內淨流入控制體的質量為
其中V2=(u,v)為流體在xy平面的速度，h=Zs-Zb為水深，dl為底邊線元。由格林公式
而控制體內質量變化率為
由(4.1)~(4.3)得
所以
對於不可壓流體，有
4.1.2 uv在z方向不均勻
若u和v在z方向不均勻，記深度平均速度為
將u和v可分解為平均量和脈動量
帶入(4.2)，應該還有一個假設，速度脈動量的積分為0，才能在積分中消去速度脈動，只剩平均值。步驟同上，不可壓的結果為
4.2 動量方程，uv均勻
x方向，控制體內動量的時間變化率為
x方向受力分析如下（水面和水底仍用Zs和Zb表示，而非ζ和h；不考慮粘性力和水面力）：
4.2.1 x方向側壓力
截面上單位寬度壓力為豎直方向壓力的積分，取水面壓力為0，有
控制體x方向側壓力為
4.2.2 x方向海底力
海底對控制體的力為水壓的反作用力，有
其中nb為Zb的法向量，見(3.8)式。
4.2.3 x方向流入量
x方向淨流入控制體的動量為
4.2.4 x方向動量方程
綜上各式去掉積分號，得
不可壓時有
此即為無粘、無科氏力、無水錶面力的x方向動量方程。y方向同理。
uv不均勻時同連續性方程，將平均量和脈動量替換原速度，以及脈動量積分為0的假設，應該能得出(4.16)式同樣的形式，其中各速度均換成平均量。
5 2維淺水方程：積分法
5.1 連續性方程，uv均勻
積分法通過將3維淺水方程中的連續方程在z方向積分得到。對不可壓縮的連續性方程(2.1)在z方向積分，有
其中
使用了萊布尼茨積分法則。第二項同理得
第三項無需使用萊布尼茲積分法則，直接積分得
將(5.2)~(5.4)代入(5.1)，結合水面和水底的邊界條件(3.3)式，可得到(4.6)式。
如果uv分解為平均量和脈動量，同樣假設脈動量的積分為0即可得到(4.9)式。
5.2 動量方程，uv不均勻
對3維淺水方程中的x方向動量方程(2.38)在z方向積分，有
5.2.1 非定常項
5.2.2 對流項
對流項第二項積分為
而其中
等號右側第二項為速度空間分佈不均引起的額外附加項，稱微分對流項或色散項，類似雷諾應力。
對流項中第一項同理，
而第三項為
對流項三項之和結合水面和底面的邊界條件(3.3)又可進一步化簡，同(5.2)~(5.4)。
5.2.3 科氏力項
非大尺度流動情況下，緯度幾乎不變，科氏力可視為常量
5.2.4 浮力項與大氣壓力項
與z無關
5.2.5 密度項
5.2.6 粘性項
密度倒數為常數，積分中先不考慮
與(3.16)式類似，令
ns與nb分別為Zs與Zb的法向量，τxs與τxb為Zs與Zb面x方向的切應力，則(5.14)式可化為
所以(3.16)式的邊界條件確實是x方向的剪切力
但是這是從積分結果得出的結論，而如何通過對微元體受力（有風吹提供的外部剪切力）的角度出發得出這一邊界條件？
5.2.7 源項
由於認為z方向無外力，故如果有外力，任何xy方向的外力均可直接積分
5.2.8 綜合
將以上積分結果對應原動量方程位置替換，即可得到動量方程。
5.3 2維淺水方程
連續性方程由(4.6)給出
動量方程為
其中
而(5.21)式中的τij又可由(1.8)和(1.9)給出。(5.21)式中還包含雷諾應力項，可認為是在(1.2)~(1.4)式進行了時間平均後進行的後續淺水方程推導。整體流程不受影響，只是粘性應力項的公式(1.8)和(1.9)中需要額外增加雷諾應力項。而在最後的(5.19)和(5.20)式中，雷諾應力項被放入了Tij裡。
記
淺水方程介紹完畢，現在能把“引”中的文字渲染成畫面嗎？
淺水方程不覺淺，曾是驚鴻照影來。驚鴻不僅是遊戲、影視作品中水面映出的驚豔角色，也是在探索流體性質乃至世界奧秘道路上無數的前輩。他們的成果並未消散，而是化作天上的星星璀璨世間，也算飛昇星域了吧。而遊戲，一定是科學與藝術、理性與感性最前沿、最成熟的曼妙融合。它的價值還遠未被髮掘完全。
驚鴻一面
做水體模擬真需要搞這些物理公式嗎？顯然並不需要。雖然行星級流動可能需要含科氏力的NS方程，但也如概述中所說，流體模擬與計算流體力學最根本的差別在於目的不同。算力時間有限，對遊戲和影視用極高精度的仿真是一種不必要的浪費。如文明帝國6(Civilization VI, Firaxis Games 2016)的海洋和河流模擬大概率是非物理的，這絲毫不影響遊戲內容。而一些簡化到不能再簡化甚至從本質上說是錯誤的模型，卻可能在水體模擬領域內有奇效。
Fournier和Reeves對此一針見血地評論道：
"We do not expect any 'physical' answer from our model. In fact, we do not ask any question except 'does it look like the real thing'. So a model which might be useless to an oceanographer might turn out to be good for our purposes. By the same token, we will feel free to adjust coefficients to achieve the needed effect, even if we are not quite sure if doing so has a physical meaning. This is of course extremely bad form when building physical models."
自己直言“糟糕的物理模型”恰好是對費米批評戴森直截了當的回應。而“隨心所欲調整係數”這一點又精準踩在費米雷點上——序中費米與戴森的談話還有後續：
Fermi asked him(Dyson) how many free parameters he had used to obtain the fit.
Smiling after being told "Four," Fermi remarked,
"I remember my old friend Johnny von Neumann used to say, with four parameters I can fit an elephant, and with five I can make him wiggle his trunk."
用四個參數馮·諾依曼都能擬合一個大象，再多給1個參數能讓大象鼻子扭起來，你們水體模擬又用了多少個參數呢？當然還是因為兩者目的不同，無關對錯。
但以上“複雜”的淺水方程真的正確嗎？注意在淺水方程的簡化裡有一個關鍵假設，是海面高度能被視為x、y、t的函數。這意味在一個確定的x、y、t下，海面高度唯一確定。但顯然海浪捲起的情況就不是這樣，Fournier和Reeves還使用神奈川衝浪裡來形象地說明這一點（想必FGO玩家對於這張圖不會陌生）。而這不僅僅是淺水方程的假設，很多其他方法也採用了這一假設。所以是複雜的模型也不一定正確，還是以上的模型遠不夠複雜？
就淺水方程而言以上也只是入門，還有很多內容如旋度勢守恆、能量方程、線性非線性淺水波等根本看不懂，非專業方向，上述內容也幾乎一定有錯誤的地方。為了學習這些入門內容已經好幾天沒開過steam。有興趣的玩家可以繼續深入探究。
打遊戲去了。
尾
“夏，你知道嗎？
這世上的路有許多條，
每個人都有屬於自己的那條。”
參考
[1] 《十日終焉》, 殺蟲隊隊員
[2] Jeff Bezos, 偏微分方程：澆滅了傑夫·貝索斯的物理夢
[3] Segre G., Hoerlin B. The Pope of Physics: Enrico Fermi and the Birth of the Atomic Age.
[4] Mehra J., Rechenberg H. The historical development of Quantum Theory
[5] 王洪偉, 我所理解的流體力學
[6] 陳懋章, 粘性流體力學基礎
[7] 劉沛清, 空氣動力學
[8] Wilcox D C .Turbulence Modeling for CFD
[9] Lamb S H .Hydrodynamics
[10] De-Leon Y, Linear waves in midlatitudes on the rotating spherical earth.
[11] J.J. Stoker. Water waves: the Mathematical Theory with Applications
[12] Arakawa, A. Computational design for long-term numerical integration of the equations for fluid motion: two-dimensional incompressible flow
[13] Heaps, N.S. Three-dimensional coastal ocean models
[14] Wubs, F.W. Numerical solution of the shallow-water equations
[15] Peachey D R .Modeling waves and surf
[16] Fedotova Z I , Khakimzyanov G S .Nonlinear dispersive equations of the shallow water on the rotating sphere
[17] Paldor N. Shallow Water Waves on the Rotating Earth
[18] Vreugdenhil C. B. Numerical Methods for Shallow-Water Flow
[19] Bates P. D. A New Method for Moving-Boundary Hydrodynamic Problems in Shallow Water.
[20] Li B., Fleming C., Three-dimensional Model of Navier-Stokes Equations for Water Waves.
[21] Vallis G K .Essentials of Atmospheric and Oceanic Dynamics
[22] Khakimzyanov G S Dispersive Shallow Water Waves: Theory, Modeling, and Numerical Methods
[23] Kubatko E .Development, implementation, and verification of hp discontinuous Galerkin models for shallow water hydrodynamics and transport.
[24] Fournier A , Reeves W T .A Simple Model of Ocean Waves.
[25] D.R. Lynch, A wave equation model for finite element tidal computations
[26] A. Oliveiria, Mass balance in Eulerian-Lagrangian transport simulations in estuaries
[27] University O .Waves, Tides & Shallow-Water Processes
[28] Deutsch.The Shallow Water Wave Equations: Formulation, Analysis and Application
[29] Roe, P. L. , Upwind Schemes Using Various Formulations of the Euler Equations, Numerical Methods for the Euler Equations of Fluid Dynamics
[30] Ouazar, D. , Computer Methods and Water Resources