纳维-斯托克斯方程

非线性

纳维－斯托克斯方程的一般形式是非线性的偏微分方程，所以在大多数实际情况下仍是如此。^[4][5]在特定情况，如一维流和斯托克斯流（又称蠕动流）下，方程可以简化为线性方程组。非线性项的出现，令大部分问题变得很难，甚至无法求解。另一方面，方程组能描述湍流，而非线性正是湍流出现的重要因素。

方程中的非线性项是对流加速（与点速度变化相关联的加速度），因此，任何对流，无论湍流与否，都会涉及非线性。有对流但仍为层流（而非湍流）的例子如下：考虑黏性液体（如油），流经一个细小并逐渐收窄的喷嘴（英语：nozzle）。此种流，不论能否确切解出，通常都能透彻研究、理解。^[6]

湍流

湍流是时变的混沌行为，这种行为常见于许多流体流动中。人们普遍认为，湍流的成因，是整个流体的惯性：时变加速度与对流加速度叠加，以产生湍流告终。因此惯性影响很小的流体，往往是层流（雷诺数量化了流所受惯性的大小）。虽然不完全确定，但一般相信纳维－斯托克斯方程能够合理地描述湍流。^[7]

纳维-斯托克斯方程关于湍流的数值解是非常难得到的，而且由于湍流之中，有多个显著不同的混合长度（英语：Turbulence modeling）尺度，若要得到稳定解，所需要的分辨率要极度精细，于是计算或直接数值模拟（英语：Direct numerical simulation）的时间长得不可行。若试图用解层流的方法来解决湍流问题，通常会得到时间不稳定的解，而不能适当收敛。为了解决这个问题，计算流体力学中，实际模拟湍流的程序，多采用雷诺平均纳维－斯托克斯方程（RANS）等时间平均方程，再辅以各湍流模型，如Spalart-Allmaras（英语：Spalart–Allmaras turbulence model）、k–ω（英语：k–omega turbulence model）、k–ε（英语：turbulence kinetic energy）、SST（英语：Menter's Shear Stress Transport），以添加另外的方程。另一种数值解法是大涡模拟（LES）。这种方法比RANS方法，占用更多计算时间和内存空间，但效果较好，因为LES明确解出较大的湍流尺度。

适用性

连同其他方程（如质量守恒定律）和良好的边界条件一并考虑时，纳维－斯托克斯方程似乎是流体运动的精确模型；甚至湍流（平均而言）也符合实际观察结果。

纳维－斯托克斯方程假定所研究的流体连续（无限可分，而不是由粒子组成），且不具相对论流速。在非常小的尺度或极端条件下，由离散分子组成的真实流体，与由纳维－斯托克斯方程描绘的连续流体，将产生不同的结果。例如，大梯度流的流体内层，有毛细现象。^[8]对于大克努森数的问题，用统计力学的玻尔兹曼方程可能更适合^[9] ，甚至要用分子动力学或其他混合方法^[10]。

另一个限制是方程的复杂性。要刻划一般常见的流体类，有行之有效的公式，但对于较罕见的类别，应用纳维－斯托克斯方程时，往往会得到非常复杂的描述，甚至是未解难题。出于这个原因，这些方程通常用于描述牛顿流体。研究这种液体是“简单”的，因为粘度模型最终被线性化；真正描述其他类型流体（如血液）的普遍模型，截至2012年还不存在。^[11]

在解释纳维－斯托克斯方程的细节之前，首先对流体的性质作几个假设。第一个假设是流体是连续的。这强调其内部没有空隙，例如，汽水中有气泡，便不属此例；同样，包含雾状粒子的气体亦不属此例。另一个必要的假设是，所有涉及到的场，即压强、速度、密度、温度等，全部都可微。

该方程从质量守恒、动量守恒、能量守恒等基本原理导出。推导过程中，经常考虑一个有限的任意体积，称为控制体积，并对其应用这些原理。该有限体积记为${\displaystyle \Omega }$ ，而其表面记为${\displaystyle \partial \Omega }$ 。该控制体积可以在空间中固定，也可能随着流体运动。这会导致一些特殊的结果，见下文。

随质导数

运动流体的属性变化，譬如大气风速的变化，可以有两种不同的方法来测量：风速仪可以固定在气象站上，也可以随气象气球飘动。第一种情况下，风速仪测量的速度，是所有运动的粒子经过一个固定点的速度；而第二种情况下，仪器测量的数值，是其随着流体运动时，速度的变化。同样的论证对于密度、温度等物理量的测量也是成立的。因此，当作微分时必须区分两种情况。第一种情况称为空间导数或者欧拉导数。第二种情况称为实质导数或拉格朗日导数。

随质导数定义为算子：

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

其中${\displaystyle \mathbf {v} }$ 是流体的速度。方程右边的第一项是普通的欧拉导数（也就是在静止参照系中的导数），而第二项表示由于流体的运动带来的变化。这个效应称为移流（英语：advection）。

L在一个控制体积${\displaystyle \Omega }$ 上守恒，可用以下积分式表示：

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

设Ω共动（随流体移动），则它随着时间而改变，所以不能将时间导数和积分简单的交换，但有：

${\displaystyle 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 .}$ 

因为这个表达式对于所有${\displaystyle \Omega }$ 成立，它可以简化为：

${\displaystyle {\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}$ 

（第一个等号用到随质导数的定义，以及对${\displaystyle \nabla \cdot \left(\mathbf {v} \mathbf {L} \right)}$ 用到导数的积法则）对于不是密度的量（因而它不必在空间中积分），${\displaystyle {\frac {\mathrm {D} }{\mathrm {D} t}}}$ 给出了正确的共动时间导数。

守恒定律

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

• 质量
• 能量
• 动量
• 角动量

并且用状态定律来闭合^{[需要解释]}，就能推导出NS方程。

连续性方程

质量的守恒写作：

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

其中${\displaystyle \rho }$ 是流体的密度。

在不可压缩流体的情况，${\displaystyle \rho }$ 是常数，不取决于时间或位置，于是方程简化为：

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

动量守恒

动量守恒写作：

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

注意${\displaystyle \mathbf {v} \otimes \mathbf {v} }$ 是一个张量，${\displaystyle \otimes }$ 代表张量积。

可以进一步简化，利用连续性方程，上式变成：

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

可以认出这就是通常牛顿第二定律${\displaystyle m\mathbf {a} =\mathbf {F} }$ 。

一般形式

方程组的形式

纳维-斯托克斯方程的一般形式是：

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

关于动量守恒。张量${\displaystyle \mathbb {P} }$ 代表施加在一个流体粒子上的表面力（应力张量）。除非流体是由象旋涡这样的旋转自由度组成，${\displaystyle \mathbb {P} }$ 是一个对称张量。一般来讲，我们有如下形式：

${\displaystyle \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}}}$ 

其中${\displaystyle \sigma }$ 是法向约束，而${\displaystyle \tau }$ 是切向约束。

迹${\displaystyle \sigma _{xx}+\sigma _{yy}+\sigma _{zz}}$ 在流体处于平衡态时为0。这等价于流体粒子上的法向力的积分为0。

我们再加上连续性方程：

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

对于处于平衡的液体，${\displaystyle \mathbb {P} }$ 的迹是3p。

其中

p是压强

最后，我们得到：

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

其中${\displaystyle \mathbb {T} }$ 是${\displaystyle \mathbb {P} }$ 的非对角线部分。

闭合问题

这些方程是不完整的。要对它们进行完备化，必须对${\displaystyle \mathbb {P} }$ 的形式作一些假设。例如在理想流体的情况${\displaystyle \tau }$ 分量为0。用于完备方程组的方程是状态方程。

再如，压强可以主要是密度和温度的函数。

要求解的变量是速度的各个分量，流体密度，静压力，和温度。流场假定为可微并连续，使得这些平衡得以用偏微分方程表达。这些方程可以转化为涡度和流函数这些次变量的威尔金森方程组。解依赖于流体的性质（例如粘滞度、比热、和热导率），并且依赖于所研究的区域的边界条件。

${\displaystyle \mathbb {P} }$ 的分量是流体的一个无穷小元上面的约束。它们代表垂直和剪切约束。${\displaystyle \mathbb {P} }$ 是对称的，除非存在非零的自旋密度。

所谓非牛顿流体就是其中该张量没有特殊性质使得方程的特殊解出现的流体

这些是问题的特定的常见简化，有时解是已知的。

牛顿流体

在牛顿流体中，如下假设成立：

${\displaystyle \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)}$ 

其中

${\displaystyle \mu }$ 是液体的粘滞度。

${\displaystyle \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)}$ 

${\displaystyle \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)}$ 

其中为简化书写，对脚标使用了爱因斯坦求和约定。

不采用简化书写的完整形式非常繁琐，分别为：

动量守恒：

${\displaystyle \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]}$ 

${\displaystyle \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]}$ 

${\displaystyle \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]}$ 

质量守恒：

${\displaystyle {\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}$ 

因为密度是一个未知数，我们需要另一个方程。

能量守恒：

${\displaystyle \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 }$ 

其中：

${\displaystyle \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}}$ 

假设一个理想气体：

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

上面是一共6个方程6个未知数的系统。（u，v，w，T，e以及　${\displaystyle \rho }$ ）。

宾汉（Bingham）流体

在宾汉流体中，我们有稍微不同的假设：

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

那些流体在开始流动之前能够承受一定的剪切。牙膏是一个例子。

幂律流体

这是一种理想化的流体，其剪切应力，${\displaystyle \tau }$ ，由下式给出

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

该形式对于模拟各种一般流体有用。

不可压缩流体

其纳维－斯托克斯方程（Navier-Stoke equation）分为动量守恒公式

${\displaystyle \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 )}} ^{\text{惯 性}}=\underbrace {-\nabla p} _{\begin{smallmatrix}{\text{压 强}}\\{\text{梯 度}}\end{smallmatrix}}+\underbrace {\mu \nabla ^{2}\mathbf {v} } _{\text{黏 滞 力}}+\underbrace {\mathbf {f} } _{\begin{smallmatrix}{\text{其 他 力}}\end{smallmatrix}}}$ 

和质量守恒公式

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

其中，对不可压缩牛顿流体来说，只有对流项（convective terms）为非线性形式。对流加速度（convective acceleration）来自于流体流动随空间之变化所产生的速度改变，例如：当流体通过一个渐缩喷嘴（convergent nozzle）时，流体产生加速之情况。由于此项的存在，对于暂态运动中的流体来说，其流场速度变化不再单是时间的函数，亦与空间有关。

另外一个重要的观察重点，在于黏滞力（viscosity）在流场中的以流体速度作拉普拉斯运算来表现。这暗示了在牛顿流体中，黏滞力为动量扩散（diffusion of momentum），与热扩散方程非常类似。

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

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

${\displaystyle \delta _{ij}}$ 是克罗内克记号。

若${\displaystyle \mu }$ 在整个流体上均匀，动量方程简化为

${\displaystyle \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)}$ 

（若${\displaystyle \mu =0}$ 这个方程称为欧拉方程；那里的重点是可压缩流和冲击波）。

如果现在再有${\displaystyle \rho }$ 为常数，我们得到如下系统：

${\displaystyle \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}}$ 

${\displaystyle \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}}$ 

${\displaystyle \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}}$ 

连续性方程（假设不可压缩性）：

${\displaystyle {\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

^[12]

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秒的模拟结果

• 雷诺数
• 马赫数
• 雷诺平均纳维－斯托克斯方程
• 纳维-斯托克斯存在性与光滑性