納維-斯托克斯方程式

納維－斯托克斯方程式（Navier-Stokes equations），是一組描述液體、空氣等流體的運動的偏微分方程式。該方程式以法國物理學家克勞德-路易·納維和愛爾蘭物理學家喬治·斯托克斯的名字命名。

納維爾－斯托克斯方程式描述了牛頓流體在運動過程中需要滿足的動量和質量守恆性質。在納維－斯托克斯方程式的一些應用中，會將該方程式與狀態方程式一同列出，以說明流體壓強、溫度和密度三者之間的關係。

納維－斯托克斯方程式是將牛頓第二運動定律應用到流體動力學之後得到的結果。與歐拉方程式不同，納維－斯托克斯方程式考慮了流體的粘性。它假設流體在運動過程中，其所受的應力是正比於速度梯度的擴散黏性力和壓力的總和。因此，納維－斯托克斯方程式可以描述有黏性流體的運動；而歐拉方程式僅能描述無粘性流的運動。在流體黏性為零時，納維-斯托克斯方程式退化為歐拉方程式。

納維－斯托克斯方程式的應用十分廣泛。它可以用於模擬天氣、洋流、管道中的水流、星系中恆星的運動、翼型周圍的氣流，也可以用於設計飛行器和車輛、研究血液循環、設計電站、分析污染效應等等。如果將納維爾－斯托克斯方程式與麥克斯韋方程組聯立，還可以研究磁流體力學。

納維－斯托克斯方程式難以求解。實用上，也只有最簡單的情況才能用這種方法獲得已知解。這些情況通常涉及穩定態（流場不隨時間變化）的非紊流，其中流體的粘滯係數很大或者其速度很小（低雷諾數）。對於更複雜的情形，例如厄爾尼諾這樣的全球性氣象系統或機翼的升力，現時僅能藉助計算機求出納維-斯托克斯方程式的數值解。這個科學領域稱為計算流體力學。

雖然紊流是日常經驗中就可以遇到的，但這類非線性問題在理論上極難求解，仍未能證明三維空間中是否總存在光滑解，甚至有界解。此問題稱為納維-斯托克斯存在性與光滑性。克雷數學學院於2000年5月21日列入七大未解難題，懸賞一百萬美元，獎勵找到數學證明或反例的任何人。

性質

非線性

納維－斯托克斯方程式的一般形式是非線性的偏微分方程式，所以在大多數實際情況下仍是如此。在特定情況，如一維流和斯托克斯流（又稱蠕動流）下，方程式可以簡化為線性方程組。非線性項的出現，令大部分問題變得很難，甚至無法求解。另一方面，方程組能描述亂流，而非線性正是亂流出現的重要因素。

方程式中的非線性項是對流加速（與點速度變化相關聯的加速度），因此，任何對流，無論亂流與否，都會涉及非線性。有對流但仍為層流（而非亂流）的例子如下：考慮黏性液體（如油），流經一個細小並逐漸收窄的噴嘴。此種流，不論能否確切解出，通常都能透徹研究、理解。

亂流

參見：亂流

亂流是時變的混沌行為，這種行為常見於許多流體流動中。人們普遍認為，亂流的成因，是整個流體的慣性：時變加速度與對流加速度疊加，以產生亂流告終。因此慣性影響很小的流體，往往是層流（雷諾數量化了流所受慣性的大小）。雖然不完全確定，但一般相信納維－斯托克斯方程式能夠合理地描述亂流。

納維-斯托克斯方程式關於亂流的數值解是非常難得到的，而且由於亂流之中，有多個顯著不同的混合長度尺度，若要得到穩定解，所需要的解像度要極度精細，於是計算或直接數值模擬的時間長得不可行。若試圖用解層流的方法來解決亂流問題，通常會得到時間不穩定的解，而不能適當收斂。為了解決這個問題，計算流體力學中，實際模擬亂流的程序，多採用雷諾平均納維－斯托克斯方程式（RANS）等時間平均方程式，再輔以各亂流模型，如Spalart-Allmaras、k–ω、k–ε、SST，以添加另外的方程式。另一種數值解法是大渦模擬（LES）。這種方法比RANS方法，佔用更多計算時間和記憶體空間，但效果較好，因為LES明確解出較大的亂流尺度。

適用性

連同其他方程式（如質量守恆定律）和良好的邊界條件一併考慮時，納維－斯托克斯方程式似乎是流體運動的精確模型；甚至亂流（平均而言）也符合實際觀察結果。

納維－斯托克斯方程式假定所研究的流體連續（無限可分，而不是由粒子組成），且不具相對論流速。在非常小的尺度或極端條件下，由離散分子組成的真實流體，與由納維－斯托克斯方程式描繪的連續流體，將產生不同的結果。例如，大梯度流的流體內層，有毛細現象。對於大克努森數的問題，用統計力學的波茲曼方程式可能更適合 ，甚至要用分子動力學或其他混合方法。

另一個限制是方程式的複雜性。要刻劃一般常見的流體類，有行之有效的公式，但對於較罕見的類別，應用納維－斯托克斯方程式時，往往會得到非常複雜的描述，甚至是未解難題。出於這個原因，這些方程式通常用於描述牛頓流體。研究這種液體是「簡單」的，因為粘度模型最終被線性化；真正描述其他類型流體（如血液）的普遍模型，截至2012年還不存在。

基本假設

在解釋納維－斯托克斯方程式的細節之前，首先對流體的性質作幾個假設。第一個假設是流體是連續的。這強調其內部沒有空隙，例如，汽水中有氣泡，便不屬此例；同樣，包含霧狀粒子的氣體亦不屬此例。另一個必要的假設是，所有涉及到的場，即壓強、速度、密度、溫度等，全部都可微。

該方程式從質量守恆、動量守恆、能量守恆等基本原理導出。推導過程中，經常考慮一個有限的任意體積，稱為控制體積，並對其應用這些原理。該有限體積記為$\Omega$，而其表面記為$\partial \Omega$。該控制體積可以在空間中固定，也可能隨着流體運動。這會導致一些特殊的結果，見下文。

隨質導數

運動流體的屬性變化，譬如大氣風速的變化，可以有兩種不同的方法來測量：風速儀可以固定在氣象站上，也可以隨氣象氣球飄動。第一種情況下，風速儀測量的速度，是所有運動的粒子經過一個固定點的速度；而第二種情況下，儀器測量的數值，是其隨着流體運動時，速度的變化。同樣的論證對於密度、溫度等物理量的測量也是成立的。因此，當作微分時必須區分兩種情況。第一種情況稱為空間導數或者歐拉導數。第二種情況稱為實質導數或拉格朗日導數。

隨質導數定義為算子：

${\frac {\mathrm {D} }{\mathrm {D} t}}(\star )={\frac {\partial (\star )}{\partial t}}+(\mathbf {v} \cdot \nabla )(\star )$

其中$\mathbf {v}$是流體的速度。方程式右邊的第一項是普通的歐拉導數（也就是在靜止參照系中的導數），而第二項表示由於流體的運動帶來的變化。這個效應稱為移流。

L在一個控制體積$\Omega$上守恆，可用以下積分式表示：

${\frac {\mathrm {d} }{\mathrm {d} t}}\int _{\Omega }\mathbf {L} \,\mathrm {d} \Omega =0.$

設Ω共動（隨流體移動），則它隨着時間而改變，所以不能將時間導數和積分簡單的交換，但有：

$0={\frac {\mathrm {d} }{\mathrm {d} t}}\int _{\Omega }\mathbf {L} \,\mathrm {d} \Omega =\int _{\Omega }{\frac {\partial }{\partial t}}\mathbf {L} \,\mathrm {d} \Omega +\int _{\partial \Omega }\mathbf {L} \left(\mathbf {v} \cdot \mathbf {n} \,\mathrm {d} \partial \Omega \right)=\int _{\Omega }\left[{\frac {\partial }{\partial t}}\mathbf {L} +\nabla \cdot \left(\mathbf {L} \mathbf {v} \right)\right]\mathrm {d} \Omega .$

因為這個表達式對於所有$\Omega$成立，它可以簡化為：

${\frac {\mathrm {D} }{\mathrm {D} t}}\mathbf {L} +\left(\nabla \cdot \mathbf {v} \right)\mathbf {L} ={\frac {\partial }{\partial t}}\mathbf {L} +\nabla \cdot \left(\mathbf {v} \mathbf {L} \right)=0$

（第一個等號用到隨質導數的定義，以及對$\nabla \cdot \left(\mathbf {v} \mathbf {L} \right)$用到導數的積法則）對於不是密度的量（因而它不必在空間中積分），${\frac {\mathrm {D} }{\mathrm {D} t}}$給出了正確的共動時間導數。

守恆定律

主條目：守恆定律

在上式中，依次取$\mathbf {L}$為下列各守恆量：

• 質量
• 能量
• 動量
• 角動量

並且用狀態定律配合一些假設（比如假設流體為理想氣體等），將各物理量寫成由各守恆量表示的形式，就能推導出NS方程式。

連續性方程式

質量的守恆寫作：

${\frac {\partial \rho }{\partial t}}+\nabla \cdot (\rho \mathbf {v} )=0$

其中$\rho$是流體的密度。

在不可壓縮流體的情況，$\rho$是常數，不取決於時間或位置，於是方程式簡化為：

$\nabla \cdot \mathbf {v} =0$

動量守恆

動量守恆寫作：

${\frac {\partial }{\partial t}}\left(\rho \mathbf {v} \right)+\nabla \cdot (\rho \mathbf {v} \otimes \mathbf {v} )=\sum \rho \mathbf {f} ,$

注意$\mathbf {v} \otimes \mathbf {v}$是一個張量，$\otimes$代表張量積。

可以進一步簡化，利用連續性方程式，上式變成：

$\rho {\frac {D\mathbf {v} }{Dt}}=\sum \rho \mathbf {f} ,$

可以認出這就是通常牛頓第二定律$m\mathbf {a} =\mathbf {F}$。

方程組

一般形式

方程組的形式

納維-斯托克斯方程式的一般形式是：

$\rho {\frac {\mathrm {D} \mathbf {v} }{\mathrm {D} t}}=\nabla \cdot \mathbb {P} +\rho \mathbf {f}$

關於動量守恆。張量$\mathbb {P}$代表施加在一個流體粒子上的表面力（應力張量）。除非流體是由象旋渦這樣的旋轉自由度組成，$\mathbb {P}$是一個對稱張量。一般來講，我們有如下形式：

$\mathbb {P} ={\begin{pmatrix}\sigma _{xx}&\tau _{xy}&\tau _{xz}\\\tau _{yx}&\sigma _{yy}&\tau _{yz}\\\tau _{zx}&\tau _{zy}&\sigma _{zz}\end{pmatrix}}=-{\begin{pmatrix}p&0&0\\0&p&0\\0&0&p\end{pmatrix}}+{\begin{pmatrix}\sigma _{xx}+p&\tau _{xy}&\tau _{xz}\\\tau _{yx}&\sigma _{yy}+p&\tau _{yz}\\\tau _{zx}&\tau _{zy}&\sigma _{zz}+p\end{pmatrix}}$

其中$\sigma$是法向約束，而$\tau$是切向約束。

跡$\sigma _{xx}+\sigma _{yy}+\sigma _{zz}$在流體處於平衡態時為0。這等價於流體粒子上的法向力的積分為0。

我們再加上連續性方程式：

${\frac {\mathrm {D} \rho }{\mathrm {D} t}}+\rho \nabla \cdot \mathbf {v} =0$

對於處於平衡的液體，$\mathbb {P}$的跡是3p。

其中

p是壓強

最後，我們得到：

$\rho {\frac {\mathrm {D} \mathbf {v} }{\mathrm {D} t}}=-\nabla p+\nabla \cdot \mathbb {T} +\rho \mathbf {f}$

其中$\mathbb {T}$是$\mathbb {P}$的非對角線部分。

閉合問題

這些方程式是不完整的，對具有不同物理性質的流體的可以有不同的解。要對它們進行完備化，必須對$\mathbb {P}$的形式作一些假設。例如在理想流體的情況$\tau$分量為0。用於完備方程組的方程式是狀態方程式。

再如，壓強可以主要是密度和溫度的函數。

要求解的變量是速度的各個分量，流體密度，靜壓力，和溫度。流場假定為可微並連續，使得這些平衡得以用偏微分方程式表達。這些方程式可以轉化為渦度和流函數這些次變量的威爾金森方程組。解依賴於流體的性質（例如粘滯度、比熱、和熱導率），並且依賴於所研究的區域的邊界條件。

$\mathbb {P}$的分量是流體的一個無窮小元上面的約束。它們代表垂直和剪切約束。$\mathbb {P}$是對稱的，除非存在非零的自旋密度。

所謂非牛頓流體就是其中該張量沒有特殊性質使得方程式的特殊解出現的流體

特殊形式

這些是問題的特定的常見簡化，有時解是已知的。

牛頓流體

主條目：牛頓流體

在牛頓流體中，如下假設成立：

$\tau _{ij}=\mu \left({\frac {\partial v_{i}}{\partial x_{j}}}+{\frac {\partial v_{j}}{\partial x_{i}}}-{\frac {2}{3}}\delta _{ij}\nabla \cdot \mathbf {v} \right)$

其中

$\mu$是液體的粘滯度。

$\rho \left({\frac {\partial \mathbf {v} }{\partial t}}+\nabla _{\mathbf {v} }\mathbf {v} \right)=\rho \mathbf {f} -\nabla p+\mu \left(\Delta \mathbf {v} +{\frac {1}{3}}\nabla \left(\nabla \cdot \mathbf {v} \right)\right)$

$\rho \left({\frac {\partial v_{i}}{\partial t}}+v_{j}{\frac {\partial v_{i}}{\partial x_{j}}}\right)=\rho f_{i}-{\frac {\partial p}{\partial x_{i}}}+\mu \left({\frac {\partial ^{2}v_{i}}{\partial x_{j}\partial x_{j}}}+{\frac {1}{3}}{\frac {\partial ^{2}v_{j}}{\partial x_{i}\partial x_{j}}}\right)$

其中為簡化書寫，對腳標使用了愛因斯坦求和約定。

不採用簡化書寫的完整形式非常繁瑣，分別為：

動量守恆：

$\rho \cdot \left({\partial u \over \partial t}+u{\partial u \over \partial x}+v{\partial u \over \partial y}+w{\partial u \over \partial z}\right)=k_{x}-{\partial p \over \partial x}+{\partial \over \partial x}\left[\mu \cdot \left(2\cdot {\partial u \over \partial x}-{\frac {2}{3}}\cdot (\nabla \cdot {\vec {v}})\right)\right]+{\partial \over \partial y}\left[\mu \cdot \left({\partial u \over \partial y}+{\partial v \over \partial x}\right)\right]+{\partial \over \partial z}\left[\mu \cdot \left({\partial w \over \partial x}+{\partial u \over \partial z}\right)\right]$

$\rho \cdot \left({\partial v \over \partial t}+u{\partial v \over \partial x}+v{\partial v \over \partial y}+w{\partial v \over \partial z}\right)=k_{y}-{\partial p \over \partial y}+{\partial \over \partial y}\left[\mu \cdot \left(2\cdot {\partial v \over \partial y}-{\frac {2}{3}}\cdot (\nabla \cdot {\vec {v}})\right)\right]+{\partial \over \partial z}\left[\mu \cdot \left({\partial v \over \partial z}+{\partial w \over \partial y}\right)\right]+{\partial \over \partial x}\left[\mu \cdot \left({\partial u \over \partial y}+{\partial v \over \partial x}\right)\right]$

$\rho \cdot \left({\partial w \over \partial t}+u{\partial w \over \partial x}+v{\partial w \over \partial y}+w{\partial w \over \partial z}\right)=k_{z}-{\partial p \over \partial z}+{\partial \over \partial z}\left[\mu \cdot \left(2\cdot {\partial w \over \partial z}-{\frac {2}{3}}\cdot (\nabla \cdot {\vec {v}})\right)\right]+{\partial \over \partial x}\left[\mu \cdot \left({\partial w \over \partial x}+{\partial u \over \partial z}\right)\right]+{\partial \over \partial y}\left[\mu \cdot \left({\partial v \over \partial z}+{\partial w \over \partial y}\right)\right]$

質量守恆：

${\partial \rho \over \partial t}+{\partial (\rho \cdot u) \over \partial x}+{\partial (\rho \cdot v) \over \partial y}+{\partial (\rho \cdot w) \over \partial z}=0$

因為密度是一個未知數，我們需要另一個方程式。

能量守恆：

$\rho \left({\partial e \over \partial t}+u{\partial e \over \partial x}+v{\partial e \over \partial y}+w{\partial e \over \partial z}\right)=\left({\partial \over \partial x}\left(\lambda \cdot {\partial T \over \partial x}\right)+{\partial \over \partial y}\left(\lambda \cdot {\partial T \over \partial y}\right)+{\partial \over \partial z}\left(\lambda \cdot {\partial T \over \partial z}\right)\right)-p\cdot \left(\nabla \cdot {\vec {v}}\right)+{\vec {k}}\cdot {\vec {v}}+\rho \cdot {\dot {q}}_{s}+\mu \cdot \Phi$

其中：

$\Phi =2\cdot \left[\left({\partial u \over \partial x}\right)^{2}+\left({\partial v \over \partial y}\right)^{2}+\left({\partial w \over \partial z}\right)^{2}\right]+\left({\partial v \over \partial x}+{\partial u \over \partial y}\right)^{2}+\left({\partial w \over \partial y}+{\partial v \over \partial z}\right)^{2}+\left({\partial u \over \partial z}+{\partial w \over \partial x}\right)^{2}-{\frac {2}{3}}\cdot \left({\partial u \over \partial x}+{\partial v \over \partial y}+{\partial w \over \partial z}\right)^{2}$

假設一個理想氣體：

$e=c_{p}\cdot T-{\frac {p}{\rho }}$

上面是一共6個方程式6個未知數的系統。（u，v，w，T，e以及　$\rho$）。

賓漢（Bingham）流體

主條目：賓漢流體

在賓漢流體中，我們有稍微不同的假設：

$\tau _{ij}=\tau _{0}+\mu {\frac {\partial v_{i}}{\partial x_{j}}},\;{\frac {\partial v_{i}}{\partial x_{j}}}>0$

那些流體在開始流動之前能夠承受一定的剪切。牙膏是一個例子。

冪律流體

主條目：冪律流體

這是一種理想化的流體，其剪切應力，$\tau$，由下式給出

$\tau =K\left({\frac {\partial u}{\partial y}}\right)^{n}$

該形式對於模擬各種一般流體有用。

不可壓縮流體

主條目：不可壓縮流體

其納維－斯托克斯方程式（Navier-Stoke equation）分為動量守恆公式

$\quad \overbrace {\rho {\Big (}\underbrace {\frac {\partial \mathbf {v} }{\partial t}} _{\begin{smallmatrix}{\text{非 穩 態}}\\{\text{加 速}}\end{smallmatrix}}+\underbrace {(\mathbf {v} \cdot \nabla )\mathbf {v} } _{\begin{smallmatrix}{\text{對 流}}\\{\text{加 速}}\end{smallmatrix}}{\Big )}} =\underbrace {-\nabla p} _{\begin{smallmatrix}{\text{壓 強}}\\{\text{梯 度}}\end{smallmatrix}}+\underbrace {\mu \nabla ^{2}\mathbf {v} } _{\text{黏 滯 力}}+\underbrace {\mathbf {f} } _{\begin{smallmatrix}{\text{其 他 力}}\end{smallmatrix}}$

和質量守恆公式

$\nabla \cdot \mathbf {v} =0$。

其中，對不可壓縮牛頓流體來說，只有對流項（convective terms）為非線性形式。對流加速度（convective acceleration）來自於流體流動隨空間之變化所產生的速度改變，例如：當流體通過一個漸縮噴嘴（convergent nozzle）時，流體產生加速之情況。由於此項的存在，對於暫態運動中的流體來說，其流場速度變化不再單是時間的函數，亦與空間有關。

另外一個重要的觀察重點，在於黏滯力（viscosity）在流場中的以流體速度作拉普拉斯運算來表現。這暗示了在牛頓流體中，黏滯力為動量擴散（diffusion of momentum），與熱擴散方程式非常類似。

$e_{ij}={\frac {1}{2}}\left({\frac {\partial u_{i}}{\partial x_{j}}}+{\frac {\partial u_{j}}{\partial x_{i}}}\right)$;

$\Delta =e_{ii}$是散度，

$\delta _{ij}$是克羅內克記號。

若$\mu$在整個流體上均勻，動量方程式簡化為

$\rho {\frac {Du_{i}}{Dt}}=\rho f_{i}-{\frac {\partial p}{\partial x_{i}}}+\mu \left({\frac {\partial ^{2}u_{i}}{\partial x_{j}\partial x_{j}}}+{\frac {1}{3}}{\frac {\partial \Delta }{\partial x_{i}}}\right)$

（若$\mu =0$這個方程式稱為歐拉方程式；那裏的重點是可壓縮流和衝擊波）。

如果現在再有$\rho$為常數，我們得到如下系統：

$\rho \left({\partial v_{x} \over \partial t}+v_{x}{\partial v_{x} \over \partial x}+v_{y}{\partial v_{x} \over \partial y}+v_{z}{\partial v_{x} \over \partial z}\right)=\mu \left[{\partial ^{2}v_{x} \over \partial x^{2}}+{\partial ^{2}v_{x} \over \partial y^{2}}+{\partial ^{2}v_{x} \over \partial z^{2}}\right]-{\partial p \over \partial x}+\rho g_{x}$

$\rho \left({\partial v_{y} \over \partial t}+v_{x}{\partial v_{y} \over \partial x}+v_{y}{\partial v_{y} \over \partial y}+v_{z}{\partial v_{y} \over \partial z}\right)=\mu \left[{\partial ^{2}v_{y} \over \partial x^{2}}+{\partial ^{2}v_{y} \over \partial y^{2}}+{\partial ^{2}v_{y} \over \partial z^{2}}\right]-{\partial p \over \partial y}+\rho g_{y}$

$\rho \left({\partial v_{z} \over \partial t}+v_{x}{\partial v_{z} \over \partial x}+v_{y}{\partial v_{z} \over \partial y}+v_{z}{\partial v_{z} \over \partial z}\right)=\mu \left[{\partial ^{2}v_{z} \over \partial x^{2}}+{\partial ^{2}v_{z} \over \partial y^{2}}+{\partial ^{2}v_{z} \over \partial z^{2}}\right]-{\partial p \over \partial z}+\rho g_{z}$

連續性方程式（假設不可壓縮性）：

${\partial v_{x} \over \partial x}+{\partial v_{y} \over \partial y}+{\partial v_{z} \over \partial z}=0$

N-S方程式的簡化版本。採用《不可壓縮流》，Ronald Panton所著第二版

注意納維－斯托克斯方程式僅可近似描述液體流，而且在非常小的尺度或極端條件下，由離散的分子和其他物質（例如懸浮粒子和溶解的氣體）的混合體組成的真實流體，會產生和納維－斯托克斯方程式所描述的連續並且齊性的液體不同的結果。依賴於問題的納森數，統計力學可能是一個更合適的方法。但是，納維－斯托克斯方程式對於很大範圍的實際問題是有效的，只要記住他們的缺陷是天生的就可以了。

程式模擬

參考MIT18086_NAVIERSTOKES

function mit18086_navierstokes
%MIT18086_NAVIERSTOKES
%    Solves the incompressible Navier-Stokes equations in a
%    rectangular domain with prescribed velocities along the
%    boundary. The solution method is finite differencing on
%    a staggered grid with implicit diffusion and a Chorin
%    projection method for the pressure.
%    Visualization is done by a colormap-isoline plot for
%    pressure and normalized quiver and streamline plot for
%    the velocity field.
%    The standard setup solves a lid driven cavity problem.

% 07/2007 by Benjamin Seibold
%            http://www-math.mit.edu/~seibold/
% Feel free to modify for teaching and learning.
%-----------------------------------------------------------------------
Re = 1e2;     % Reynolds number
dt = 1e-2;    % time step
tf = 4e-0;    % final time
lx = 1;       % width of box
ly = 1;       % height of box
nx = 90;      % number of x-gridpoints
ny = 90;      % number of y-gridpoints
nsteps = 10;  % number of steps with graphic output
%-----------------------------------------------------------------------
nt = ceil(tf/dt); dt = tf/nt;
x = linspace(0,lx,nx+1); hx = lx/nx;
y = linspace(0,ly,ny+1); hy = ly/ny;
[X,Y] = meshgrid(y,x);
%-----------------------------------------------------------------------
% initial conditions
U = zeros(nx-1,ny); V = zeros(nx,ny-1);
% boundary conditions
uN = x*0+1;    vN = avg(x)*0;
uS = x*0;      vS = avg(x)*0;
uW = avg(y)*0; vW = y*0;
uE = avg(y)*0; vE = y*0;
%-----------------------------------------------------------------------
Ubc = dt/Re*([2*uS(2:end-1)' zeros(nx-1,ny-2) 2*uN(2:end-1)']/hx^2+...
      [uW;zeros(nx-3,ny);uE]/hy^2);
Vbc = dt/Re*([vS' zeros(nx,ny-3) vN']/hx^2+...
      [2*vW(2:end-1);zeros(nx-2,ny-1);2*vE(2:end-1)]/hy^2);

fprintf('initialization')
Lp = kron(speye(ny),K1(nx,hx,1))+kron(K1(ny,hy,1),speye(nx));
Lp(1,1) = 3/2*Lp(1,1);
perp = symamd(Lp); Rp = chol(Lp(perp,perp)); Rpt = Rp';
Lu = speye((nx-1)*ny)+dt/Re*(kron(speye(ny),K1(nx-1,hx,2))+...
     kron(K1(ny,hy,3),speye(nx-1)));
peru = symamd(Lu); Ru = chol(Lu(peru,peru)); Rut = Ru';
Lv = speye(nx*(ny-1))+dt/Re*(kron(speye(ny-1),K1(nx,hx,3))+...
     kron(K1(ny-1,hy,2),speye(nx)));
perv = symamd(Lv); Rv = chol(Lv(perv,perv)); Rvt = Rv';
Lq = kron(speye(ny-1),K1(nx-1,hx,2))+kron(K1(ny-1,hy,2),speye(nx-1));
perq = symamd(Lq); Rq = chol(Lq(perq,perq)); Rqt = Rq';

fprintf(', time loop\n--20%%--40%%--60%%--80%%-100%%\n')
for k = 1:nt
   % treat nonlinear terms
   gamma = min(1.2*dt*max(max(max(abs(U)))/hx,max(max(abs(V)))/hy),1);
   Ue = [uW;U;uE]; Ue = [2*uS'-Ue(:,1) Ue 2*uN'-Ue(:,end)];
   Ve = [vS' V vN']; Ve = [2*vW-Ve(1,:);Ve;2*vE-Ve(end,:)];
   Ua = avg(Ue')'; Ud = diff(Ue')'/2;
   Va = avg(Ve);   Vd = diff(Ve)/2;
   UVx = diff(Ua.*Va-gamma*abs(Ua).*Vd)/hx;
   UVy = diff((Ua.*Va-gamma*Ud.*abs(Va))')'/hy;
   Ua = avg(Ue(:,2:end-1));   Ud = diff(Ue(:,2:end-1))/2;
   Va = avg(Ve(2:end-1,:)')'; Vd = diff(Ve(2:end-1,:)')'/2;
   U2x = diff(Ua.^2-gamma*abs(Ua).*Ud)/hx;
   V2y = diff((Va.^2-gamma*abs(Va).*Vd)')'/hy;
   U = U-dt*(UVy(2:end-1,:)+U2x);
   V = V-dt*(UVx(:,2:end-1)+V2y);
   
   % implicit viscosity
   rhs = reshape(U+Ubc,[],1);
   u(peru) = Ru\(Rut\rhs(peru));
   U = reshape(u,nx-1,ny);
   rhs = reshape(V+Vbc,[],1);
   v(perv) = Rv\(Rvt\rhs(perv));
   V = reshape(v,nx,ny-1);
   
   % pressure correction
   rhs = reshape(diff([uW;U;uE])/hx+diff([vS' V vN']')'/hy,[],1);
   p(perp) = -Rp\(Rpt\rhs(perp));
   P = reshape(p,nx,ny);
   U = U-diff(P)/hx;
   V = V-diff(P')'/hy;
   
   % visualization
   if floor(25*k/nt)>floor(25*(k-1)/nt), fprintf('.'), end
   if k==1|floor(nsteps*k/nt)>floor(nsteps*(k-1)/nt)
      % stream function
      rhs = reshape(diff(U')'/hy-diff(V)/hx,[],1);
      q(perq) = Rq\(Rqt\rhs(perq));
      Q = zeros(nx+1,ny+1);
      Q(2:end-1,2:end-1) = reshape(q,nx-1,ny-1);
      clf, contourf(avg(x),avg(y),P',20,'w-'), hold on
      contour(x,y,Q',20,'k-');
      Ue = [uS' avg([uW;U;uE]')' uN'];
      Ve = [vW;avg([vS' V vN']);vE];
      Len = sqrt(Ue.^2+Ve.^2+eps);
      quiver(x,y,(Ue./Len)',(Ve./Len)',.4,'k-')
      hold off, axis equal, axis([0 lx 0 ly])
      p = sort(p); caxis(p([8 end-7]))
      title(sprintf('Re = %0.1g   t = %0.2g',Re,k*dt))
      drawnow
   end
end
fprintf('\n')

%=======================================================================

function B = avg(A,k)
if nargin<2, k = 1; end
if size(A,1)==1, A = A'; end
if k<2, B = (A(2:end,:)+A(1:end-1,:))/2; else, B = avg(A,k-1); end
if size(A,2)==1, B = B'; end

function A = K1(n,h,a11)
% a11: Neumann=1, Dirichlet=2, Dirichlet mid=3;
A = spdiags([-1 a11 0;ones(n-2,1)*[-1 2 -1];0 a11 -1],-1:1,n,n)'/h^2;

模擬結果

執行MATLAB 即可觀察4秒的模擬結果

參看

• 雷諾數
• 馬赫數
• 雷諾平均納維－斯托克斯方程式
• 納維-斯托克斯存在性與光滑性

參考文獻

• Inge L. Rhyming Dynamique des fluides, 1991 PPUR.
• Polyanin A.D., Kutepov A.M., Vyazmin A.V., Kazenin D.A., Hydrodynamics, Mass and Heat Transfer in Chemical Engineering, Taylor & Francis, London, 2002. ISBN 0-415-27237-8.

外部連結

• 克雷數學研究院納維-斯托克斯方程式大獎
  • 該問題的正式命題
• 納維-斯托克斯方程式的一個推導
• 納維-斯托克斯方程式的推導 （頁面存檔備份，存於互聯網檔案館）
• NASA關於納維-斯托克斯方程式的網頁
• 納維-斯托克斯方程式（一些精確解） （頁面存檔備份，存於互聯網檔案館），位於EqWorld:數學方程式的世界