# 纳维-斯托克斯方程

## 基本假设

### 随質导数

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

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 {\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 {\frac {\partial \rho }{\partial t}}+\nabla \cdot (\rho \mathbf {v} )=0}$

${\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 \rho {\frac {D\mathbf {v} }{Dt}}=\sum \rho \mathbf {f} ,}$

## 方程组

### 一般形式

#### 方程组的形式

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

${\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 _{xx}+\sigma _{yy}+\sigma _{zz}}$ 在流体处于平衡态时为0。这等价于流体粒子上的法向力的积分为0。

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

p是压强

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

#### 闭合问题

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

### 宾汉（Bingham）流体

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

### 幂律流体

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

### 不可壓缩流體

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

${\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 \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所著第二版

## 程式模擬

[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;


## 参考文献

1. ^ McLean, Doug. Continuum Fluid Mechanics and the Navier-Stokes Equations. Understanding Aerodynamics: Arguing from the Real Physics. John Wiley & Sons. 2012: 13–78 （英语）. The main relationships comprising the NS equations are the basic conservation laws for mass, momentum, and energy. To have a complete equation set we also need an equation of state relating temperature, pressure, and density...
2. ^ Millennium Prize Problems—Navier–Stokes Equation, claymath.org, Clay Mathematics Institute, March 27, 2017 [2017-04-02], （原始内容存档于2015-12-22） （英语）
3. ^ Fefferman, Charles L. Existence and smoothness of the Navier–Stokes equation (PDF). claymath.org. Clay Mathematics Institute. [2017-04-02]. （原始内容 (PDF)存档于2015-04-15） （英语）.
4. ^ Fluid Mechanics（Schaum's Series）, M. Potter, D.C. Wiggert, Schaum's Outlines, McGraw-Hill (USA), 2008, ISBN 978-0-07-148781-8
5. ^ Vectors, Tensors, and the basic Equations of Fluid Mechanics, R. Aris, Dover Publications, 1989, ISBN (10) 0-486-66110-5
6. ^ Parker, C. B. McGraw Hill Encyclopaedia of Physics 2nd. 1994. ISBN 0-07-051400-3 （英语）.
7. ^ Encyclopaedia of Physics (2nd Edition), R.G. Lerner, G.L. Trigg, VHC publishers, 1991, ISBN (Verlagsgesellschaft) 3-527-26954-1, ISBN（VHC Inc.）0-89573-752-3
8. ^ Gorban, A.N.; Karlin, I. V. Beyond Navier–Stokes equations: capillarity of ideal gas. Contemporary Physics (Review article). 2016, 58 (1): 70–90. Bibcode:2017ConPh..58...70G. S2CID 55317543. . doi:10.1080/00107514.2016.1256123.
9. ^ Cercignani, C. The Boltzmann equation and fluid dynamics. Friedlander, S.; Serre, D. (编). Handbook of mathematical fluid dynamics 1. Amsterdam: North-Holland. 2002: 1–70. ISBN 978-0444503305 （英语）.
10. ^ Nie, X.B.; Chen, S.Y.; Robbins, M.O. A continuum and molecular dynamics hybrid method for micro-and nano-fluid flow. Journal of Fluid Mechanics (Research article). 2004, 500: 55–64 [2021-10-12]. Bibcode:2004JFM...500...55N. doi:10.1017/S0022112003007225. （原始内容存档于2018-08-09） （英语）.
11. ^ Öttinger, H.C. Stochastic processes in polymeric fluids. Berlin, Heidelberg: Springer Science & Business Media. 2012. ISBN 9783540583530. doi:10.1007/978-3-642-58290-5 （英语）.
12. ^ Chen, Weijia; Chen, JC. Combined compact difference method for solving the incompressible Navier-Stokes equations. International Journal for Numerical Methods in Fluids. 2011-06-06, 68 (10). ISSN 0271-2091. doi:10.1002/fld.2602.
• 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.