Simulation
Simulation
模擬。真實現象化作數值,以數值模仿真實現象。
Mechanics
力學。距離、時間、質量,討論三者之間的關係。
mechanics - statics 靜力學:有力,呈穩態(處處總和為零) - kinematics 運動學:無力 - dynamics = kinects 動力學:有力
System
系統。大量物體,相互牽連。
system - structure 結構系統:不動 - mechanism 機械系統:會動
UVa 11574
Particle Simulation
Particle Simulation
牛頓發現運動定律F=ma。歐拉發現物體運動可以簡化成質心運動,物體軌跡等於質心軌跡。從質點開始介紹,事情比較簡單。
質點移動
質點。擁有位置、速度、加速度。
移動。加速度改變了速度,速度改變了位置。
兩種方式:數值解(數值模擬)、符號解(公式解)。
兩者優缺點互補,通常採用數值模擬。優點:簡單明快。缺點:誤差逐步累積。
步法。數值模擬,Δt並非無限微小,導致軌跡逐漸偏斜。大家發明各種步法,盡量撥正軌跡。
兩種步法:Euler method、Verlet method。
兩者優缺點互補,通常採用Euler method。優點:簡單明快。缺點:軌跡較歪。
時間。以無窮迴圈,更新時間、更新位置、重新繪圖。
兩種更新時間方式:定量增加(時間差Δt是定值)、浮動增加(時間t是當前時刻,時間差Δt是前後兩幀時刻差距)。
兩者優缺點互補,建議採用定量增加。優點:步伐穩定,軌跡端正。缺點:動畫播放速度不符合真實時間流逝。
繪圖。質點畫成什麼形狀都行。質點位置視作像素座標。
質點施力
推進力。F = ma。
摩擦力。速度乘以摩擦係數等於摩擦力。
彈簧
彈力。可以類比為推進力。
虎克定律:位置差異乘以彈性係數等於彈力。
阻尼力。可以類比為摩擦力。
速度乘以阻尼係數等於阻尼力。
當力不在水平垂直方向,必須求得水平分量、垂直分量。
阻尼力,當速度不在彈簧方向,必須投影至彈簧方向。
位移、旋轉、縮放
幾何變換有三個獨立因素:位移、旋轉、縮放。
先前只介紹位移。此處同時使用位移、旋轉、縮放。
活用三種獨立因素,調整彈簧參數,製作各種動畫特效:
參數設定
參數。質量、長度、時間、彈性係數、阻尼係數,設定成真實世界的數值。讓螢幕座標、真實世界座標兩者一致。我沒有研究。
兩個需要考慮的項目:像素寬度與真實長度的比例(dpi)、重新繪圖的時間間隔(frame rate)。
Rigid Body Simulation
Rigid Body Simulation
資源
Physics engine list, books, PhD thesis, on-line resources https://pybullet.org/Bullet/phpBB3/viewtopic.php?t=63 pyBullet tutorial https://pybullet.org/wordpress/index.php/forum-2/ Box2D tutorial https://box2d.org/publications/ Simbody paper https://simtk.org/docman/?group_id=47 bepuphysics tutorial https://www.bepuentertainment.com/
(Equation)(Under Construction!)
Body / Motion
物體。質點連成一串、連成一片、連成一團。
運動。質點位置改變。
歐拉發現物體運動可以視作質心運動。質心是質點座標平均數。
Inertia / Force
慣量。改變物體運動狀態的難易度。
力。改變物體運動狀態。
牛頓發現物體速度恆定不變。除非施力,才能改變速度。
科普教育當中,物體是主角,質量是關鍵參數。質量意義是「物體與生俱來的量」,不增不減。
靜力學當中,重量是主角,質量是係數。質量意義是「重量除以重力的比值」。質量越大,比值越大。
動力學當中,位置是主角,慣量是係數。慣量意義是「阻礙位置變化的程度」。慣量越大,位置變化越少。
當位置變化來自位移運動,慣量恰是質量。當位置變化來自旋轉運動,慣量稱作「轉動慣量」或者「慣性矩」,需要特別計算。
旋轉慣量取決於物體造型。每個質點分別計算慣量,通通加總。維基百科整理了一份列表。
Translational Motion / Rotational Motion
數學幾何變換有三個因素:位移、旋轉、縮放。
物理運動有兩個因素:位移運動(移動)、旋轉運動(轉動)。
(縮放運動沒有特地命名。請見柔體運動章節。)
質點運動有一種幾何變換:位移。 剛體運動有兩種幾何變換:位移、旋轉。 柔體運動有三種幾何變換:位移、旋轉、縮放。
兩種運動各有一套物理量。內容對應,名稱有變。
物理量有三類:位置變化、運動量、運動量變化。
位移運動(移動)的物理量:
name | notation 丨中文翻譯丨意義 ---------------------|--------------丨一一一一丨一一一一 position | x 丨位置 丨 velocity | v = dx/dt 丨速度 丨位置變化 acceleration | a = dv/dt 丨加速度 丨速度變化 ---------------------|--------------丨一一一一丨一一一一 mass | m 丨質量 丨位移慣量 momentum | p = mv 丨動量 丨 energy | E = ½mv² 丨能量 丨 ---------------------|--------------丨一一一一丨一一一一 mass flow rate | ṁ = dm/dt 丨質量流量丨質量變化 force | F = dp/dt 丨力 丨動量變化 power | P = dE/dt 丨功率 丨能量變化
旋轉運動(轉動)的物理量:
angular position | θ 丨角位置 丨 angular velocity | ω = dθ/dt 丨角速度 丨角位置變化 angular acceleration | α = dω/dt 丨角加速度丨角速度變化 ---------------------|--------------丨一一一一丨一一一一一一 rotational inertia | I 丨轉動慣量丨 angular momentum | L = Iω = r×p 丨角動量 丨 rotational energy | E = ½Iω² 丨轉動動能丨 ---------------------|--------------丨一一一一丨一一一一一一 ? | İ = dI/dt 丨? 丨轉動慣量變化 torque | τ = Iα = r×F 丨力矩 丨角動量變化 rotational power | P = τ∙ω 丨轉動功率丨轉動動能變化
位移運動(移動)的物理量,補充內容:
distance | s = ‖x' - x‖ 丨距離 丨位置差異(純量) displacement | r = x' - x 丨移位 丨位置差異(向量) ---------------------|--------------丨一一一一丨一一一一一一一一 speed | v = ‖v‖ 丨速率 丨位置變化(純量) velocity | v = dx/dt 丨速度 丨位置變化(向量) ---------------------|--------------丨一一一一丨一一一一一一一一 impulse | J = ∫ F dt 丨衝量 丨動量差異(向量) work | W = ∫ F∙dr 丨功 丨能量差異(純量)
旋轉運動(轉動)有三個軸,總共三組物理量。
Linear Motion / Circular Motion
https://en.wikipedia.org/wiki/Circular_motion https://en.wikipedia.org/wiki/Acceleration#Tangential_and_centripetal_acceleration
位移運動分為直線運動和圓周運動,產生直線軌跡和曲線軌跡。
旋轉運動分為直線運動和曲線運動,在球面上產生直線軌跡與曲線軌跡。
數學幾何變換,只考慮起點終點,不考慮軌跡。
物理運動,需要考慮慣量。
Simple Motion / Complex Motion
硬幣在地上滾動,包含兩種運動。轉軸轉動,轉軸移動。當硬幣往旁邊偏,那麼移動軌跡從直線變成曲線。
陀螺有兩種運動:轉軸做轉動spin,轉軸做移動之圓周運動precession。
當物體形如球拍,物體每隔一段時間會翻面,稱作Dzhanibekov effect。當物體形如陀螺,物體會側滾,稱作Gyroscopic precession。
滑動Slipping、滾動Rolling
static friction | cs 丨靜摩擦力丨 kinetic friction | ck 丨動摩擦力丨
質點的圓周運動,可以用來計算物體的旋轉運動。
radius | r 丨半徑(不是位移) -------------------------|------------------丨一一一一一一一一 angular position | θ = s /r 丨角位置 angular velocity | ω = vt/r 丨角速度 angular acceleration | α = at/r 丨角加速度 -------------------------|------------------丨一一一一一一一一 tangential velocity | vt = ωr 丨切線速度 tangential acceleration | at = αr 丨法線加速度 centripetal acceleration | ac = ω²r = vt²/r 丨向心加速度 centripetal force | fc 丨向心力
Collision / Contact
兩個物體運動,質點相互緊靠,運動量相互轉移。
彈性碰撞、非彈性碰撞。
一、偵測碰撞。
二、接觸位置。
三、動量轉移。
Action Force / Reaction Force
兩個物體相互施力。一個物體當主角,一個物體當配角。主角施予配角的力,稱作作用力。配角施予主角的力,稱作反作用力。
作用力Action Force。
gravity | g = Gm₁m₂/r² 丨引力 丨朝向質心的作用力 lift | 丨升力 丨流體介質,速度垂直方向的作用力 centripetal force | 丨向心力丨圓周運動,速度垂直方向的作用力
反作用力Reaction force:與作用力同時發生,方向相反。
drag | d = k₁v + k₂v² 丨阻力 丨速度方向的反作用力 friction | 丨摩擦力丨固體介面,速度方向的反作用力 damping | 丨阻尼 丨振盪,速度方向的反作用力 buoyancy | f = ρVg 丨浮力 丨流體介質,引力方向的反作用力 tension | 丨張力 丨繩索,施力方向的反作用力 restoring force | F = -ks 丨恢復力丨彈簧,施力方向的反作用力 tensile force | 丨抗拉力丨材質,施力方向的反作用力
Fracture
破碎。一個物體變成多個物體。
(Rotation)(Under Construction!)
Rotation
旋轉。既是物理學詞彙,亦是數學詞彙。
物理學當中,旋轉區分為二維旋轉和三維旋轉。
數學當中,旋轉可以寫成矩陣。正規正交矩陣,容積+1。二維旋轉是2x2矩陣。三維旋轉是3×3矩陣。
4個數字、9個數字。難以記憶、難以控制。必須簡化。
二維旋轉的簡易表示法:一個轉角θ。
設定兩個向量,從一個向量旋轉至另一個向量。兩個向量的夾角,作為轉角。運用三角函數,求得矩陣數值。
[ cosθ -sinθ ] [ sinθ cosθ ]
三維旋轉的簡易表示法:一個轉軸e、一個轉角θ。
設定兩個向量,從一個向量旋轉至另一個向量。兩個向量形成一個平面。兩個向量的叉積,得到平面法向量,作為轉軸。兩個向量的夾角,作為轉角。旋轉化作自旋。
Spin
自旋:一個轉軸(向量)、一個轉角(純量),沿著轉軸旋轉。
自旋有三種計算方式:向量、矩陣、四元數。
v' = v cosθ + (e × v) sinθ + e (e ∙ v) (1 - cosθ)
矩陣:進一步使用Rodrigues' formula求得變換矩陣。
v' = R v R = I + sinθ E + (1 - cosθ) E² [ 0 -ez ey ] E = [ ez 0 -ex ] [ -ey ex 0 ] (Ev = e × v)
當轉軸是標準座標軸,變換矩陣是Givens matrix。
[ 1 0 0 ] Rx(θ) = [ 0 cosθ -sinθ ] roll [ 0 sinθ cosθ ] (spin around x-axis) [ cosθ 0 -sinθ ] Ry(θ) = [ 0 1 0 ] pitch [ sinθ 0 cosθ ] (spin around y-axis) [ cosθ -sinθ 0 ] Rz(θ) = [ sinθ cosθ 0 ] yaw [ 0 0 1 ] (spin around z-axis)
四元數:我不知道公式名稱。
a + bi + cj + dk。a b c d是實數。i j k是虛數。
虛數乘法規則如下所示。
| 1 i j k --+--------------- 1 | 1 i j k i² = j² = k² = ijk = -1 i | i −1 k −j ij = k , jk = i , ki = j j | j −k −1 i ji = -k , kj = -i , ik = -j k | k j −i −1
自旋是四元數的共軛乘法q p q⁻¹。
向量 v = (vx,vy,vz) <=> v = 0 + vxi + vyj + vzk 轉軸 e = (ex,ey,ez) <=> e = 0 + exi + eyj + ezk 轉角 θ 自旋 v' = q v q⁻¹ q = exp(e ½θ) = cos(½θ) + exsin(½θ)i + eysin(½θ)j + ezsin(½θ)k
多個四元數可以併成一個四元數,如同矩陣。
v' = q₃ q₂ q₁ v q₁⁻¹ q₂⁻¹ q₃⁻¹ = q v q⁻¹ where q = q₃ q₂ q₁
進一步使用Euler–Rodrigues formula求得變換矩陣。
[ a²+b²-c²-d² 2(bc-ad) 2(bd+ac) ] v' = [ 2(bc+ad) a²+c²-b²-d² 2(cd-ab) ] v [ 2(bd-ac) 2(cd+ab) a²+d²-b²-c² ] a = cos(½θ) b = ex sin(½θ) c = ey sin(½θ) d = ez sin(½θ) (a²+b²+c²+d² = 1)
Euler Angles
根據Euler rotation theorem:一次旋轉等價於一次自旋。
然而目前沒有機械裝置能夠指定轉軸方向(我沒聽過)。
因此大家將一次旋轉化作三次自旋。以標準座標軸作為轉軸。
然而卻產生了gimbal lock,存在無法抵達的旋轉。稍後介紹。
三次自旋有兩種方式:
Euler angles: z-x-z 第一轉軸、第三轉軸,兩者相同。 Tait–Bryan angles: z-y-x 標準座標軸依序作為轉軸。
實務上使用Tait–Bryan angles。三次自旋是三個轉軸與三個轉角。轉軸採用標準座標軸,互相垂直。轉角稱作偏航角yaw、俯仰角pitch、橫滾角roll。
轉軸動向有兩種方式:
extrinsic rotation:沿著世界座標軸自旋(所有轉軸固定不動)。 intrinsic rotation:沿著物體座標軸自旋(其餘轉軸隨之旋轉)。
矩陣連乘順序、四元數乘法順序,左右顛倒,extrinsic rotation就變成intrinsic rotation。
Gimbal Lock
環架:三個環,對應三個轉軸。當一個環自旋,所有外環轉軸固定不動,相當方便。紀錄三個環的三個旋角,完成旋轉。
環架的缺點:當其中兩個環對齊,則存在無法抵達的旋轉。這個現象稱作gimbal lock。
歐拉角和四元數也有相同缺點。無論轉軸是世界座標軸和物體座標軸,都無法避免此缺點。
簡易理解方式:寫下三個自旋矩陣,然後設定適當的轉角,使得其中一個自旋矩陣形成恆等矩陣。剩下兩個自旋矩陣相乘結果,恰好退化成一個自旋矩陣。自旋方位短少一個。自由度從3變2。
let β = π/2 [ 1 0 0 ] [ cosβ 0 -sinβ ] [ cosγ -sinγ 0 ] R = [ 0 cosα -sinα ] [ 0 1 0 ] [ sinγ cosγ 0 ] [ 0 sinα cosα ] [ sinβ 0 cosβ ] [ 0 0 1 ] roll pitch yaw [ 1 0 0 ] [ 1 0 0 ] [ cosγ -sinγ 0 ] = [ 0 cosα -sinα ] [ 0 1 0 ] [ sinγ cosγ 0 ] [ 0 sinα cosα ] [ 0 0 1 ] [ 0 0 1 ] [ 0 0 1 ] = [ sin(α+γ) cos(α+γ) 0 ] [ -cos(α+γ) sin(α+γ) 0 ]
gimbal lock引發許多問題,例如登月太空船Apollo 11因而難以瞄準星星進行定位。太空人當場表示:給我第四個環當作聖誕禮物好不好。
想要避免gimbal lock,一種方式是追加一個自旋,另一種方式是混用extrinsic rotation和intrinsic rotation。【尚待確認】
(Collision)(Under Construction!)
Collision / Contact / Intersection
名詞collision。動詞collide。物理學詞彙。力學詞彙。
碰撞。兩個形狀,相互接觸。
名詞contact。動詞contact。物理學詞彙。力學詞彙。
碰觸。兩個形狀,接觸地點。
名詞intersection。數學詞彙。幾何學詞彙。
交集。兩個形狀,重疊區域。
動詞intersect。數學詞彙。幾何學詞彙。
相交。兩個形狀,相互交錯、相互重疊。
「偵測碰撞」等同「偵測交集是否為空集合」。
碰撞止於形狀邊界。相交深入形狀內部。電腦做計算,一切皆數學。儘管我們談論物理碰撞模擬,但是此處介紹數學相交演算法。
Narrow-phase Collision
兩個物體碰撞。根據物體形狀,有著各種演算法。
凸多面體有著高速演算法。大家把一個多面體切割成多個凸多面體,以便降低計算時間。
是否碰撞、侵入距離。
Convex Polygon Collision
二維凸多邊形碰撞。O(N)。
樸素算法:二維碰撞即是某兩條線段碰撞。枚舉各種線段配對,偵測是否碰撞。O(N²)。
不知名算法:兩個凸多邊形的交集,恰是凸多邊形。兩個凸多邊形的所有邊,涵蓋了交集的所有邊。保留交集的邊,捨棄其餘的邊,以此設計演算法。
一個凸多邊形的所有邊,恰好依照角度排序。實施merge sort,合併兩個凸多邊形的所有邊,依照角度排序。
依序掃描所有邊,每當凸多邊形更替,就偵測前邊與後邊是否碰撞。若碰撞,後邊是交集;若無碰撞,後邊不是交集。
上述過程是multi-pass。亦可改成one-pass,一邊合併、一邊偵測碰撞。
Convex Polyhedron Collision
三維凸多面體碰撞,有SAT和GJK兩種演算法,都是O(N²)。
樸素算法:三維碰撞即是某兩個三角形碰撞。枚舉各種三角形配對,偵測是否碰撞。O(N²)。
separating axis theorem:三角形碰撞改成影子碰撞。叉積改成點積,節省時間。
一、首先介紹影子。
三角形視作牆壁。另一個三角形視作飛鏢。碰撞視作飛鏢穿牆。
三角形頂點,投影到牆壁法向量,偵測影子碰撞。
二、三角形碰撞、影子碰撞,兩者不等價。
三角形為主角: 一、三角形碰撞,保證影子碰撞。 二、三角形無碰撞,不保證影子有無碰撞。 影子為主角: 三、影子碰撞,不保證三角形有無碰撞。 四、影子無碰撞,保證三角形無碰撞。
此演算法以影子為主角。因此採用第四個性質:影子無碰撞。
三、設計演算法。
凸多面體A,每一個面分別作為牆壁。
凸多面體A與B,每一個點皆投影到牆壁法向量。
偵測影子無碰撞。法向量朝外視作正向。投影處越外視作越高。A中最高點低於B中最低點,或者,B中最高點低於A中最低點。
四、援引hyperplane separation theorem。
兩個凸的形狀,能以任意超平面隔開,即無碰撞。
一旦影子無碰撞,立即結束演算法,節省時間。
Gilbert–Johnson–Keerthi algorithm:枚舉順序,朝向交集。以此減少枚舉數量,節省時間。
一、首先介紹Minkowski difference。
形狀A、形狀B,兩者的Minkowski difference,是形狀C。
A中一點,B中一點,兩者座標相減,得到C中一點。
A中一點,B中一點,枚舉各種配對,得到整個形狀。
二、如何利用Minkowski difference偵測碰撞呢?
兩個凸多面體的Minkowski difference,恰是一個凸多面體。
證明很簡單。凸的定義,寫成數學式子,相減一下就行了。
兩個凸多面體碰撞,那麼Minkowski difference包含原點。沒有碰撞,那麼不含原點。
證明很簡單。交點,座標相減是0。無交點,座標相減非0。
最後得到一個偵測碰撞演算法:先計算Minkowski difference,再判斷是否包含原點。
三、如何找出Minkowski difference,並且偵測碰撞呢?
樸素算法:求出整個Minkowski difference。
取巧算法:求出部分Minkowski difference,節省時間。
一邊建立四面體剖分,一邊移動至相鄰四面體。逐步靠近原點。
首先求出任意四個頂點,形成四面體。每回合找到一個新頂點,並且移動至相鄰四面體。逐步靠近原點。
Minkowski difference恰是凸的,四面體剖分恰是連通的,得以使用貪心法。只要逐步靠近原點,終會抵達最近之處。
補充說明一下。二維的情況下,移動過程可以視作Minkowski difference所有頂點實施binary search,找到最靠近原點的頂點。
所有頂點形成一條鏈。每回合找到一個新頂點,將鏈截成兩段,並且保留靠近原點的那一段。
四、如何找出Minkowski difference的新頂點呢?
找出一個朝向原點的方向向量。A中最高點、B中最低點,兩者相減,得到Minkowski difference的新頂點。
這宛如SAT。本來枚舉A的面,現在枚舉A與B的頂點。本來使用法向量,現在使用朝向交集的方向向量。由於走向原點,只需枚舉部分頂點,不需枚舉所有面,節省時間。
還可以順手偵測影子無碰撞,立即結束演算法。由於走向一致,只需檢查A中最高點低於B中最低點,不需檢查B中最高點低於A中最低點,又節省一半時間。
五、如何找出朝向原點的方向向量,進而求得新頂點呢?
求第一點:使用任意方向向量。例如X軸方向。
求第二點:使用反方向。
求第三點:前兩點形成線段。再與原點形成三角形。三角形法向量、線段向量,兩者叉積,得到朝向原點的向量。
求第四點:前三點形成三角形。三角形法向量,正向與反向,找出朝向原點者。
求後續各點:前四點形成四面體。三個新三角形、一個舊三角形。三個新三角形法向量,找出朝向原點者。任一皆可。
六、補充一下。GJK也能計算最短距離。
兩個凸多面體的最短距離,即是Minkowski difference與原點的最短距離。
四面體的移動過程當中,隨時計算四面體與原點的最短距離。
Farthest Point in a Direction
SAT和GJK都用到了「給定一個方向向量,找到凸多面體中的最高點」。
樸素算法:凸多面體的所有點,投影到方向向量,得到投影量,取最大值。O(N)。
凸多邊形:預先切成上鏈與下鏈。兩鏈都是單峰函數,無論方向向量為何。兩鏈分別實施bisection method求得極值,其中一個極值是最大值。O(logN)。
凸多面體:預先切成上殼與下殼。兩殼都是單峰函數,無論方向向量為何。使用quadtree儲存頂點,方便實施bisection method。O(N)。
Simple Polygon Collision
二維簡單多邊形碰撞。O(NlogN)。
Bentley–Ottmann algorithm:直接視作大量線段相交問題,直接使用其演算法。平移的掃描線,輔以二元搜尋樹。
Simple Polyhedron Collision
三維簡單多面體碰撞,好像是O(N²)。
Lin–Canny Algorithm:我沒有研究。
V-Clip:我沒有研究。
Bounding Volume
每個物體不分青紅皂白,實施SAT或CJK偵測碰撞,O(N)。
替物體設計簡易外框,得以迅速偵測碰撞,O(1),節省時間。外框沒有碰撞,保證物體沒有碰撞,大可不必實施SAT或CJK。外框有碰撞,才實施SAT或CJK,偵測物體是否碰撞。
外框偵測碰撞,些微增加時間;物體偵測碰撞,大幅減少時間。剛體運動當中,少數時候有碰撞,多數時候沒有碰撞。利大於弊。
外框造型有三種:
bounding circle / bounding sphere:體積最小的球。紀錄中心座標、半徑。平均O(N)。
bounding box / bounding volume:體積最小的長方體,細分為AABB和OBB。
axis-aligned bounding box (AABB):必須對齊世界座標軸。紀錄中心座標、長寬高。O(N)。
oriented bounding box (OBB):不必對齊世界座標軸。記錄中心座標、長寬高、轉角。二維O(NlogN)。三維O(N³)。
剛體運動,適合採用OBB。OBB外框最貼切,誤判最少,節省最多時間。儘管OBB時間複雜度極高,但是剛體不會變形,只需預先計算一次。
柔體運動,不適合採用OBB。柔體隨時變形,重算太花時間。
Circle Collision / Sphere Collision
二維圓形碰撞、三維球體碰撞。O(1)。
樸素算法:圓心連線。若碰撞,圓心距離小於等於半徑相加。
Axis-aligned Box Collision / Axis-aligned Cuboid Collision
二維擺正長方形碰撞、三維擺正長方體碰撞。O(1)。
樸素算法:SAT。「A的X座標最大值」小於「B的X座標最小值」則無碰撞。AB對調,如法炮製。Y座標、Z座標,如法炮製。
Box Collision / Rectangular Cuboid Collision
二維長方形碰撞、三維長方體碰撞。O(1)。
樸素算法:SAT。A轉正,B隨之旋轉,檢查座標極值大小。B轉正,如法炮製。
Broad-phase Collision
大量物體碰撞。使用區域資料結構,儲存所有物體,節省時間。
樸素算法:枚舉所有物體配對。建立O(1),偵測O(N²)。
uniform grid:找出每個物體的bounding sphere的直徑的最大值。建立網狀方格,格子寬度是直徑。每當發生碰撞,物體一定位在相同格子或相鄰格子。實施bucket sort,所有物體存入對應格子。建立O(N+Kᴰ),偵測平均O(N+3ᴰ)。
quadtree / octree:當相同區域有多個物體,就遞迴分割區域,直到每個區域只有一個物體。首先(由葉往根)逐步擴大區域寬度,直到超過前述直徑。然後(由根往葉)檢查轄下所有區域所有物體。建立O(N+K2ᴰ),偵測平均O(N+K)。
bounding region hierarchy / bounding volume hierarchy:物體計算AABB。外框依照物體位置排序。相鄰兩個外框不斷併成一個外框,形成二元樹。(由根往葉)先檢查大外框碰撞,再檢查小外框碰撞。建立O(DNlogN),偵測平均O(N)。
當物體座標隨時改變,要嘛每次重算,要嘛使用動態資料結構。我沒有研究。Box2D作者正在實作,預計出現於下一版本。我到時候再來研究。
(Particle)(Under Construction!)
Collision Resolution
朝著法向量,將接觸點挪到表面。
讓侵略者盡快離開,給予額外動量。侵入距離越深,動量越大,宛如彈簧。
Tunneling / Speculative Collision
一個物體,速度太快。另一個物體,形狀太薄。導致一個物體穿越另一個物體,偵測不到碰撞。
物體移動範圍,形成外框。以此外框實施碰撞偵測。
剛體運動當中,只有移動和轉動。移動是平行四邊形,轉動是扇形。
Linear Complementarity Problem
Soft Body Simulation(Under Construction!)
Soft Body Simulation(Deformable Body Simulation)
備忘
n-body https://en.wikipedia.org/wiki/Barnes–Hut_simulation cloth http://www.cs.cornell.edu/projects/YarnCloth/ hair http://www.andyselle.com/papers/9/ soft object http://www.alecrivers.com/fastlsm/ water wave http://www.cemyuksel.com/research/waveparticles/ matter http://chemists.princeton.edu/torquato/
資源
https://en.wikipedia.org/wiki/Soft_body_dynamics https://www.seas.upenn.edu/~cffjiang/
(Equation)(Under Construction!)
elasticity / plasticity
彈性。物體會恢復形狀。
塑性。物體不會恢復形狀。
deformation
物體視作質點相連。運動視作質點連動。縮放運動水平於質點連接方向。旋轉運動垂直於質點連接方向。
縮放運動、旋轉運動,通通視作位移運動。
旋轉運動化作位移運動之圓周運動。一個質點作為圓心,另一個質點作為圓周。
縮放運動化作位移運動之彈簧運動。質點距離視作彈簧長度。
柔體運動視作質點位移運動,稱作deformation。
物理學沒有「縮放運動」這個詞彙。我們懷念他。
剛體運動和柔體運動採用了完全不同的策略。剛體運動,細分為移動和轉動,形成兩套相似的物理量。柔體運動,一律視作移動,只有一套物理量。
http://ddg.cs.columbia.edu/SIGGRAPH05/Shells.pdf deformation mapping: x̕ = f(x) transformation of point deformation gradient: t̕ = ∇f(t) transformation of tangent vector strain energy: ‖t̕‖² - ‖t‖² difference of squared length ‖t̕‖² - ‖t‖² = ‖∇f(t)‖² - ‖t‖² = tᵀ ∇fᵀ∇f t - tᵀt = tᵀ (∇fᵀ∇f - I) t
force: F displacement: ΔL stress: σ = F/A ---> df/dx² area: A strain: ε = ΔL/L ---> du/dx = ∇ₓu length: L stress-strain curve: hooke: σ = Eε young modulus / stiffness tensor: E hooke: f = kΔL strain energy: U = 1/2Vσε = 1/2fΔL elastic energy: U = 1/2kΔL²
deformation: F(x) deformation gradient (transformation matrix): ∇F(x)ᵀ = J ---> F displacement gradient: ∇u = Fᵀ-I Cauchy infinitesimal strain tensor (even part): E = [(F-I)+(F-I)ᵀ]/2 = (F+Fᵀ)/2-I = (∇u + ∇uᵀ)/2 Cauchy infinitesimal rotation tensor (odd part): R = [(F-I)-(F-I)ᵀ]/2 = (F-Fᵀ)/2-I Green finite deformation tensor: C = FᵀF Green finite strain tensor: E = (C-I)/2 = (FᵀF-I)/2 = (∇u + ∇uᵀ + ∇uᵀ∇u)/2 strain->stress: Eᵀ = E
positive eigenvalue in strain tensor -> lengthen negative eigenvalue in strain tensor -> shorten reflection in rotation tensor -> shorten?
Poisson's ratio: dεy / dεₓ Lamé's coefficients:
Navier–Cauchy Equation 1. strain-displacement: E = 1/2(∇u + ∇uᵀ) 2. equilibrium: div(S) = 0 and S = Sᵀ 3. stress-strain (hooke): S = C[E] stiffness tensor (order-4): C
Stretch / Bend / Twist
(Particle)(Under Construction!)
全域、局部
全域求解:物理限制就是物理等式。所有等式變成超大矩陣。
局部XPBD:位置差異、容積差異,作為彈力。以參數調整力的大小。
https://matthias-research.github.io/pages/tenMinutePhysics/ 更新位置 x' = x + v dt 追加限制:修正位置 x" = x' + solve(C, x') 新位置與舊位置差異,反推新速度 v' = (x" - x) / dt
最佳化與物理現象
物理 丨數學 一一一一十一一一一一一一一一一 力 丨向量 合力 丨線性組合(加權總和) 如果有很多力,可以直接乘上倍率再相加。 力>位能丨積分 力對位置積分=位能。 位能>力丨微分 位能對位置微分=力。 位能 丨二次型 拋物面函數、凸函數,有唯一最小值。
物理 丨數學 丨位能 力 一一一一十一一一一一一十一一一一一一一一一一一一一一一一一一一一 彈力 丨微分(梯度)丨½ k (x - l)² 丨k (x - l) 約束力 丨微分(梯度)丨½ k c(x)² 丨k c(x) c′(x) 慣性運動丨梯度等於零 丨min 位能 丨力 = 0 靜力平衡丨梯度等於零 丨min 位能總和 丨合力 = 0 施加外力丨拉格朗日乘數丨min 位能 st 約束能丨力 + λ ⋅ 約束力 = 0
例如彈簧系統。統計所有彈簧的彈力位能,求最小值,求最小值所在位置(每個彈簧的端點座標)。梯度等於零,就是靜力平衡,就是穩定狀態。實際應用諸如建築結構分析、化學結構分析。
例如軌道運動。物體位置限制在軌道上面,動能與位能的總和是常數。根據拉格朗日乘數,限制能量大小,等同施加外力;此例當中就是施加向心力。實際應用諸如太空衛星。
相關應用:
Structural Optimization: Topology Optimization。調整模型結構。 Molecular Modeling: Energy Minimization。調整化學結構。
位置與速度盡量相同,能量必須固定
min sum ‖xi' - xi‖² + ‖vi' - vi‖² subject to c(x) = 0 x x'是先後位置,v v'是先後速度。c(x)是能量。
二次函數有著高速演算法:Newton–Lagrange Method與Projected Gradient Method兩者合體。
乘子法-->微分方程組--->手動微分--->係數和乘數取出來變向量--->解大矩陣
qk+1 = argmin 1/2 ‖q - qₖ‖²D subject to c(qₖ) + ∇c(qₖ)ᵀ(q - qₖ) = 0 q solve λk+1: [∇c(qₖ)ᵀ D⁻¹ ∇c(qₖ)] λk+1 = c(qₖ) solve qk+1: qk+1 = qk - D⁻¹ ∇c(qk) λk+1
《FEPR: Fast Energy Projection for Real-Time Simulation of Deformable Objects》
Cloth
Water Wave
Strand
Stretching
Fluid Simulation
Fluid Simulation
演算法一覽
流體模擬有兩種策略:
粒:紀錄每一個水滴的位置、速度、加速度(力)。
格:紀錄每一個位置的水滴總數量、總速度、總加速度(力)。
物理學術語分別是拉格朗日描述、歐拉描述。
粒格衍生許多演算法:粒連續SPH、粒離散SEC、粒方程PBD、粒化格PIC、格八方LBM、格中粒MAC、格中格MG。
近期的演算法:PCISPH、PIC/FLIP、Affine PIC、……。
關鍵字一覽
本文並未涉及這些內容。
粒子間作用力 Lennard–Jones potential、intermolecular force 內聚力、附著力 cohesion、adhesion 浮力、摩擦力 buoyancy、friction 障礙物 immersed boundary method、Werner–Wengle wall function 多相流 multiphase flow 交融 Rayleigh–Taylor instability、Kelvin–Helmholtz instability 層流、亂流 laminar、turbulence (mean + fluctuation) 無因次量 Reynolds number、Froude number 溫度、能量 temperature、energy 角動量 angular monentum、vorticity、helicity 流線 streamline、pathline、streakline、timeline 底形 bedform 氣泡流、氣塊流 bubble flow、slug flow 顆粒流 granular flow
(Grid)
格
格子。建立40×40個格子。
只紀錄速度。不紀錄數量:假設粒子密度均勻不變,假設數量是常數。不紀錄加速度:加速度乘上時間,直接添加到速度。
配置。兩種方式。
colocated grid:速度在格中央。 staggered grid:X速度在格左右,Y速度在格上下。
此處採用staggered grid。兩點優勢:一、水衝撞牆,會往兩側散開,不會只往其中一側流動。二、計算散度,會更加精確,不會省略中央速度。
上下添加一排X速度,以便內插上下牆壁附近的X速度。左右添加一排Y速度,以便內插左右牆壁附近的Y速度。
牆壁。四種方式。數學術語是邊界條件。
no-slip condition:牆壁沾黏,切線速度是特定值。 free-slip condition:牆壁光滑,切線速度等於鄰水速度。 inflow condition:牆壁插管,法線速度是特定值。 outflow condition:牆壁虛設,法線速度等於鄰水速度。
此處採用no-slip。水被牆壁帶動,牆壁切線速度是自訂數值;水不可穿牆,牆壁法線速度是零。
牆內速度、牆外速度,內插得到牆壁速度。以此反推牆外速度。
運動。考慮兩種現象。
diffuse:擴散。格中粒子均勻四散,其速度也隨之均勻四散。 ∂ ∂ ∂ ∂ ∂ ∂ ―― u = ν ( ――― u + ――― u ) ―― v = ν ( ――― v + ――― v ) ∂t ∂x² ∂y² , ∂t ∂x² ∂y² ^^^ ^^^ 這兩個是擴散速率。 小寫希臘字母ν,唸作/njuː/,寫作nu。不是小寫英文字母v。 flow:流動。格中粒子依照自身速度集體移動,其速度也隨之集體移動。 ∂ ∂ ∂ ∂ ∂ ∂ ―― u = - ( ―― uu + ―― vu ) ―― v = - ( ―― uv + ―― vv ) ∂t ∂x ∂y , ∂t ∂x ∂y
數學術語是Navier–Stokes Equation與Euler Method。
步幅。移動不得穿越格子,否則速度可能高速逸散。保守策略是移動距離小於格子寬度。
數學術語是穩定條件CFL condition。格子邊長Δx越窄,時間差距Δt越小。格子越細膩,步幅越短小,動畫越緩慢。
不可壓縮。格中粒子速度,出入總量相等,流動順暢。
incompressibility:不可壓縮。 ∂ ∂ ―― u + ―― v = 0 ∂x ∂y
按理來說,擴散、流動、不可壓縮必須同時聯立,一口氣求得新速度。不幸的是,大家尚未發現漂亮的解法。別無選擇,只好逐個計算。首先擴散與流動,接著校正速度場,使之不可壓縮。
數學術語是散旋諧分解。依序計算散度、散勢(壓力)、散場。速度場減去散場,使得速度場的散度變成零,成為無散場。
數學術語是投影。給定速度場,找到平方誤差最小的無散場。
鄰居。投影時,牆壁不計入。牆壁壓力設為零,作為無效的分子;事先算好鄰居有幾格,作為分母。
求解。一次方程組求解,兩種演算法。
SOR:稀疏矩陣,適合鬆弛法。容易實作,但是收斂速度慢。 PCG:半正定矩陣,適合共軛梯度法。請見這篇講義。
數學術語是Discrete Poisson Equation。
繪圖。速度場通常畫成軌跡圖、箭矢圖。亦得入境隨俗:
粒:散布大量粒子,沿著速度場跑。 格:建立粒子數量場(粒子濃度場),沿著速度場跑。
速度場繪圖方式:軌跡圖
速度場繪圖方式:箭矢圖
速度場繪圖方式:粒
粒子。隨機散布大量粒子。
繪圖。粒子畫成什麼形狀都行。
流動。粒子沿著速度場跑。總共需要考慮:內插、步幅、步法、牆壁。
內插。取得更精確的速度。速度位於格線,加0.5微調位置。
步幅。採用更大的步幅,以便製作即時動畫。
步法。三種步法。以邪門歪道撥正軌跡。大家習慣採用one-step method,簡單走一步。數學術語是積分方法。
one-step method 1. 抓本地速度,走一步。 predictor-corrector method 1. 抓本地速度,走一步,到終點。 2. 抓本地與終點的平均速度,從頭重新走一步。 midpoint method 1. 抓本地速度,走半步,到中點。 2. 抓中點速度,從頭重新走一步。
牆壁。Δt畢竟不是無限微小,靠牆的水容易一步穿牆。穿牆粒子,放回牆內。
速度場繪圖方式:格
格子。建立粒子數量場。
繪圖。每一格根據粒子數量(粒子濃度),決定顏色深淺。
流動。粒子數量沿著速度場跑。總共需要考慮:內插、步幅、步向、步法、牆壁、內插。
內插。取得更精確的速度。速度位於格線,加0.5微調位置。
步幅。採用更大的步幅,以便製作即時動畫。
步向。兩種方向。
forward 前進,往後補值。 優點是動量守恆。 缺點是空格與疊格。解法是縮小步幅、減少Δt。 backward 倒退,往前取值。 優點是流向精確。 缺點是零速度處無法流入、零速度處無限流出。解法是強制擴散。
步法。大家習慣採用backward步向、one-step method步法,兩者合稱semi-Lagrangian method。
牆壁。Δt畢竟不是無限微小,靠牆的水容易一步穿牆。一旦穿牆,立即挪回牆內。
內插。更精確的粒子數量。副作用是擴散。兩種方式。
前進與分配:根據面積大小分攤。 倒退與內插:一次內插、或者Catmull–Rom interpolation。 牆壁事先設定粒子數量,方便內插。
補充說明:步幅增長
先前章節,計算速度場,避免增長步幅,以滿足CFL condition;繪製速度場,作弊增長步幅,以便讓動畫順暢。然而此舉導致繪圖結果不符合真實流體。
事實上,計算速度場,可以增長步幅,但是必須改變擴散與流動的演算法。請見這篇文章。
步幅。計算速度場、繪製速度場,步幅一致,不必作弊。
擴散。速度向四周擴散。
數學術語是implicit Euler method。
流動。粒子數量場沿著速度場跑,速度場也跟著一起跑。
數學術語是semi-Lagrangian method。
守恆。Δt變大,流動之後誤差變大,各個物理量不再守恆。必須校正速度場,令物理量守恆。兩種方法。
質量守恆:流動前後,粒子數量總和不變。 動量守恆:流動前後,粒子動量總和不變。 能量守恆:流動前後,粒子動能總和不變。
projection method 1. flow,走一步。V1 = V0 + V0 ⋅ Δt 2. project,減掉散場。V1p = V1 - div_part(V1) reflection method [Zehnder et al 2018] 1. flow,走半步。V½ = V0 + V0 ⋅ ½Δt 2. project,減掉散場。V½p = V½ - div_part(V½) 3. reflect,再減掉一次散場。V½r = V½ - div_part(V½) - div_part(V½) 4. flow,走半步。V1 = V½r + V½p ⋅ ½Δt 5. project,減掉散場。V1p = V1 - div_part(V1)
校正。Δt變大,流動之後誤差變大。有人校正速度場。
modified MacCormack method 1. 更新速度場。V0 -> V1 2. 速度變號,再更新速度場。V1 -> V2 理論上回復原狀:V0 = V2 實際上產生誤差:V0 ≠ V2 3. 來回更新兩次,因此更新一次的誤差是一半:(V2 - V0) / 2 4. 正確答案是更新一次,再扣掉更新一次的誤差:V1 - (V2 - V0) / 2 5. 裁切答案,不超過V1四周的最大值和最小值。
格:空氣與水
水域。以粒或格,紀錄水域資訊,沿著速度場跑。四種解法。
marker particles 水域建立排列整齊的大量粒子,沿著速度場跑。 檢查每一格,若格中駐留足夠粒子,則視作水域。 粒子流動,間距改變,需要重置粒子。 front tracking 水面建立網格,沿著速度場跑。 網格流動,間距改變,需要修正網格。 間距太近則刪除一點,間距太遠則插入新點。 volume of fluid 建立水量場,沿著速度場跑。 格子流動,水量分配給目的地周遭四格,根據面積比例。 level set 水面到格子的距離,建立距離場,沿著速度場跑。 距離場流動,距離必錯,需要修正距離場。 重新估計水面位置,重新計算水面到格子的距離。
此處採用level set。建立距離場,水中為負值,水外為正值,水面為零。
我的實作,牆壁的距離值設為半格寬度。我沒查到相關資料。
重算距離場。兩個步驟。
水面:距離異號的相鄰兩格,做一次內插,找到距離為零之處。 水內水外:從水面開始BFS遍歷每一格,並且計算距離。 八鄰居取距離最小者,只取已有距離的格子。
我的實作,牆壁與水沒有做一次內插。我沒查到相關資料。
外插。如果採用倒退步法,則空氣必須有速度,而水才能流入空氣。故以水域速度設定空氣速度。兩種解法。
最短距離:找到距離最近的水面格子,速度設成相同。 逐層擴散:從水面開始,一層層設定速度。 八鄰居的速度平均數,只取已有速度的格子。
我的實作,距離場的計算順序不正確。既不是最短路徑Dijkstra,也不是逐層設定BFS。導致速度外插有時候優先遭遇距離較小、但是離水面較遠的孤立格子。目前的應急措施是塞回queue。
不可壓縮。水域才投影,水外不必投影。
求解。水域才求解,水外不必求解。空氣壓力為零,須採計。
牆壁速度:方便起見設為鄰水速度外插。 理論上是零。 空氣速度:方便起見設為水面速度外插。 理論上是零。 牆壁壓力:不計,方便起見設為零。 理論上是水壓:水深乘以密度。 空氣壓力:採計,方便起見設為零。 理論上是一大氣壓。
重力。垂直速度,添加常數。
補充說明:redistance
redistance可以跟project並排,視作距離場校正。redistance也可以跟air並排,視作距離場外插,視作flow的事前準備工作。我覺得後者是比較好的。以程式碼結構來說,速度場外插和距離場外插是可以合併一起做的。
補充說明:multigrid level set
level set是舊方法了。有個新方法是multigrid level set。請參考《Real-Time Eulerian Water Simulation Using a Restricted Tall Cell Grid》。然而這是NVIDIA的專利。
補充說明:advect命名錯誤
許多流體模擬文獻,flow誤植為advect,流動誤認為平流。
flow是流動(連續方程式)。project讓流動變得不可壓縮。advect是平流(平流方程式)。flow和project合體才是advect。
(Equation)
微觀與宏觀
粒:壓縮力、黏滯力。格:擴散、流動、不可壓縮。
粒微觀,格宏觀。有趣的是,迄今無人發現微觀與宏觀之間的數學關係式。
擴散包含了壓縮力和黏滯力。然而大家習慣將擴散視作黏滯,將雷諾數視作黏滯程度。雖然不完全正確,但是容易理解。
Continuity Equation
連續方程式。描述流動守恆的方程式。一個地點的物理量之所以增加減少,是因為跟隔壁跑來跑去。不會無中生有,不會無端消失。
時間變化。這個時刻、上個時刻的物理量差距。
物理量可以是質量、體積、力、………。
∂ ―― q ---> q̇ ∂t
空間變化。這個地點、隔壁地點的物理量差距。考慮流動速度。
物理量必須適合各維度加總,例如質量、體積、………。
∂ ∂ ―― qu + ―― qv ---> div(qu) ∂x ∂y
連續方程式。時間變化=空間變化。
∂ ∂ ∂ ―― q + ( ―― qu + ―― qv ) = 0 ---> q̇ + div(qu) = 0 ∂t ∂x ∂y ^^^^ ^^^^^^^^^^^^^^^^^ time space difference
Advection Equation
平流方程式。水位差距與流量成正比。水位差距會移動。
連續方程式,拆成平流與壓縮。微分的乘法公式:前微後不微、前不微後微。
∂ ∂ ∂ ∂ ∂ ―― q + ( u ―― q + v ―― q ) + q ( ―― u + ―― v ) = 0 ∂t ∂x ∂y ∂x ∂y ^^^^^^^^^^^^^^^^^^^ ^^^^^^^^^^^^^^^ advection compression
平流。物理量往隔壁流動,差距越大,流量越大。
∂ ∂ ( u ―― q + v ―― q ) ---> u dot grad(q) ∂x ∂y
壓縮。物理量聚集或者散開。
∂ ∂ q ( ―― u + ―― v ) ---> q div(u) ∂x ∂y
不可壓縮。壓縮項為零,速度的散度為零。速度出入總量相等,不生不滅;物質流動順暢,不推不擠。
∂ ∂ q ( ―― u + ―― v ) = 0 ---> q div(u) = 0 ∂x ∂y
平流方程式。時間變化=平流。
∂ ∂ ∂ ―― q + ( u ―― q + v ―― q ) = 0 ---> q̇ + (u dot grad(q)) = 0 ∂t ∂x ∂y ^^^^ ^^^^^^^^^^^^^^^^^^^ time advection
有人簡寫成大寫D,稱為「質點導數」。
D ―― q = 0 Dt
Navier–Stokes Equation
一、密度:連續方程式。
時間變化=空間變化。
∂ ∂ ∂ ―― ρ + ( ―― ρu + ―― ρv ) = 0 ∂t ∂x ∂y
二、動量密度:連續方程式,再追加兩種力。
時間變化+空間變化=自身位勢(壓力)+相互動勢(應力)。
{ ∂ ∂ ∂ ∂ ∂ ∂ { ―― ρu + ( ―― ρuu + ―― ρvu ) = - ―― p + ( ―― τ + ―― τ ) { ∂t ∂x ∂y ∂x ∂x ˣˣ ∂y ˣʸ { { ∂ ∂ ∂ ∂ ∂ ∂ { ―― ρv + ( ―― ρuv + ―― ρvv ) = - ―― p + ( ―― τ + ―― τ ) { ∂t ∂x ∂y ∂x ∂x ʸˣ ∂y ʸʸ ^^^^^ ^^^^^^^^^^^^^^^^^^^ ^^^^ ^^^^^^^^^^^^^^^^^^^ time space intra inter difference difference force force
三、能量密度:我不熟悉,就不講了。
Navier–Stokes Equation for Water
水是流體的特例。因此方程組可以簡化。
一、水滿足牛頓黏滯定律:應力與速度梯度成正比。
因此應力簡化成擴散。
∂ ∂ ―― τ + ―― τ ∂x ˣˣ ∂y ˣʸ ∂ ∂ ∂ ∂ = ―― ( μ ―― u ) + ―― ( μ ―― u ) ∂x ∂x ∂y ∂y ∂ ∂ = μ ( ――― u + ――― u ) ---> μ laplace(u) ∂x² ∂y²
二、水幾乎不可壓縮:速度散度為零。
因此空間變化簡化成平流。
∂ ∂ ―― ρuu + ―― ρvu ∂x ∂y ^^^^^^^^^^^^^^^ space difference ∂ ∂ ∂ ∂ = ( u ―― ρu + v ―― ρu ) + ρu ( ―― u + ―― v ) ∂x ∂y ∂x ∂y ^^^^^^^^^^^^^^^^^^^^^ ^^^^^^^^^^^^^^^ advection compression = 0 ∂ ∂ = ( u ―― ρu + v ―― ρu ) ---> u dot grad(ρu) ∂x ∂y
三、水幾乎密度均勻不變:密度是常數。
因此時空變化可以一併提出密度。
∂ ∂ ∂ ―― ρu + ( u ―― ρu + v ―― ρu ) ∂t ∂x ∂y ∂ ∂ ∂ = ρ [ ―― u + ( u ―― u + v ―― u ) ] ---> ρ [ u̇ + (u dot grad(u)) ] ∂t ∂x ∂y
水的Navier–Stokes Equation。教科書通常只介紹它。
{ ρ [ u̇ + (u dot grad(u)) ] = - grad(p) + μ laplace(u) { div(u) = 0
{ ρ ( u̇ + u·∇u ) = -∇p + μ ∆u { ∇·u = 0
{ u̇ + u·∇u = -∇p/ρ + ν ∆u (let μ/ρ = ν) { ∇·u = 0
Navier–Stokes Equation for 2D Water
水的Navier–Stokes Equation,二維的情況下,有兩種形式。
一、壓力與速度(原式):可精確模擬層流。
pressure-velocity formulation { u̇ + (u dot grad(u)) = -(grad(p) / ρ) + ν laplace(u) { div(u) = 0
二、流函數與渦旋:可精確模擬亂流。例如自然現象leapfrogging vortex ring、microburst。
streamfunction-vorticity formulation { ̇ω - (grad(ψ) cross grad(ω)) = ν laplace(ω) { laplace(ψ) = -ω
推導過程請見這篇講義。原式套用旋度,速度改成位勢。
Schrödinger Equation for Water
計算幾何有個手法,讓實數變複數:二維座標(x,y)改成複數x+yi,二維座標位移、縮放、旋轉改成複數加法、倍率、乘法。
流體力學也有個手法,讓實數變複數:Navier–Stokes equation改成Schrödinger equation,三維實數改成二維複數,質量(密度)改成複數長度平方,速度改成複數角度。請見這篇論文。