纳维-斯托克斯方程
纳维尔-斯托克斯方程(),以法國工程師兼物理學家克劳德-路易·纳维、愛爾蘭物理學和數學家乔治·斯托克斯兩人命名,是一组偏微分方程,描述液体和空气等流体的運動。
纳维尔-斯托克斯方程表達了牛頓流體運動時,動量和質量守恆。有時,還連同状态方程列出,說明流體壓強、溫度、密度之間的關係。[1]方程斷言,流体粒子动量的改变率(力),來自作用在液体内部的压力变化、耗散粘滞力、以及重力。其中粘滞力类似于摩擦力,产生于分子的相互作用,越黏的流體,該作用就越強。这样,纳维-斯托克斯方程描述作用于液体任意给定区域的力的动态平衡。
学术研究和经济生活中,許多重要物理過程都可用纳维尔-斯托克斯方程描述,因此該些方程有很重要的研究价值。它们可以用于模拟天气、洋流、管道中的水流、星系中恒星的运动、翼型周围的气流,也可以用于设计飞行器和车辆、研究血液循环、设计电站、分析污染效应等等。納-斯方桯組與馬克士威方程組聯立,用於研究磁流體力學。
纳维-斯托克斯方程依赖微分方程来描述流体的运动。不同于代数方程,这些方程不寻求建立所研究的变量(譬如速度和壓力)的关系,而寻求建立这些量的变化率或通量之间的关系。用数学术语来讲,这些变化率对应于变量的导数。其中,在零粘滞度的最简单情况下,纳维-斯托克斯方程化為歐拉方程,表明加速度(速度的导数,或者说变化率)與内部压力的导数成正比。
这表示对于给定的物理问题,至少要用微积分才可以求得其纳维-斯托克斯方程的解。实用上,也只有最简单的情况才能用这种方法获得已知解。这些情况通常涉及稳定态(流场不随时间变化)的非紊流,其中流体的粘滞系数很大或者其速度很小(低雷诺数)。
对于更复杂的情形,例如厄尔尼诺这样的全球性气象系统或机翼的升力,現時僅能借助计算机求出纳维-斯托克斯方程的數值解。这个科学领域称为计算流体力学。
虽然紊流是日常经验中就可以遇到的,但这类非线性问题在理論上极难求解,仍未能證明三維空間中是否總存在光滑解,甚至有界解。此問題稱為納維-斯托克斯存在性與光滑性。克雷数学学院于2000年5月21日列入七大未解難題,懸賞一百萬美元,奖励找到數學證明或反例的任何人。[2][3]
性质
非线性
纳维-斯托克斯方程的一般形式是非线性的偏微分方程,所以在大多数实际情况下仍是如此。[4][5]在特定情况,如一维流和斯托克斯流(又稱蠕动流)下,方程可以简化为线性方程组。非线性項的出現,令大部分问题變得很难,甚至無法求解。另一方面,方程組能描述湍流,而非線性正是湍流出現的重要因素。
方程式中的非线性項是对流加速(与点速度变化相关联的加速度),因此,任何对流,无论湍流与否,都会涉及非线性。有對流但仍為層流(而非湍流)的例子如下:考慮黏性液體(如油),流經一個細小並逐漸收窄的噴嘴。此種流,不論能否確切解出,通常都能透徹研究、理解。[6]
湍流
湍流是时变的混沌行为,这种行为常见于许多流体流动中。人们普遍认为,湍流的成因,是整個流體的慣性:时變加速度與对流加速度疊加,以產生亂流告終。因此惯性影响很小的流体,往往是层流(雷诺数量化了流所受惯性的大小)。虽然不完全确定,但一般相信纳维-斯托克斯方程能够合理地描述湍流。[7]
纳维-斯托克斯方程关于湍流的数值解是非常难得到的,而且由于湍流之中,有多個显著不同的混合长度尺度,若要得到穩定解,所需要的分辨率要極度精细,於是计算或直接數值模擬的时间長得不可行。若试圖用解层流的方法来解决湍流问题,通常会得到时间不稳定的解,而不能适当收敛。为了解决这个问题,計算流體力學中,實際模擬湍流的程序,多採用雷諾平均納維-斯托克斯方程(RANS)等时间平均方程,再辅以各湍流模型,如Spalart-Allmaras、k–ω、k–ε、SST,以添加另外的方程。另一种數值解法是大涡模拟(LES)。这种方法比RANS方法,佔用更多計算时间和記憶體空間,但效果較好,因为LES明确解出較大的湍流尺度。
适用性
連同其他方程(如质量守恒定律)和良好的边界条件一併考慮時,纳维-斯托克斯方程似乎是流体运动的精确模型;甚至湍流(平均而言)也符合实际观察结果。
纳维-斯托克斯方程假定所研究的流体連續(无限可分,而不是由粒子组成),且不具相對論流速。在非常小的尺度或极端条件下,由离散分子组成的真实流体,与由纳维-斯托克斯方程描繪的连续流体,将产生不同的结果。例如,大梯度流的流體內層,有毛细现象。[8]對於大克努森數的问题,用统计力学的波茲曼方程式可能更適合[9] ,甚至要用分子动力学或其他混合方法[10]。
另一个限制是方程的复杂性。要刻劃一般常見的流体類,有行之有效的公式,但對於較罕見的類別,應用纳维-斯托克斯方程時,往往会得到非常复杂的描述,甚至是未解難題。出于这个原因,这些方程通常用于描述牛顿流体。研究这种液体是“简单”的,因为粘度模型最终被线性化;真正描述其他类型流体(如血液)的普遍模型,截至2012年还不存在。[11]
基本假设
在解释纳维-斯托克斯方程的细节之前,首先对流体的性质作几个假设。第一个假设是流体是连续的。这强调其内部沒有空隙,例如,汽水中有氣泡,便不屬此例;同樣,包含雾状粒子的氣體亦不屬此例。另一个必要的假设是,所有涉及到的场,即压强、速度、密度、温度等,全部都可微。
该方程从质量守恆、动量守恆、能量守恆等基本原理导出。推導過程中,經常考虑一个有限的任意体积,称为控制体积,並對其應用这些原理。该有限体积记为,而其表面记为。该控制体积可以在空间中固定,也可能随着流体运动。这会导致一些特殊的结果,見下文。
随質导数
运动流体的属性变化,譬如大气风速的变化,可以有两种不同的方法来测量:风速仪可以固定在气象站上,也可以隨气象气球飄動。第一种情况下,风速仪测量的速度,是所有运动的粒子经过一个固定点的速度;而第二种情况下,仪器测量的數值,是其随着流体运动时,速度的变化。同样的论证对于密度、温度等物理量的测量也是成立的。因此,当作微分时必须区分两种情况。第一种情况称为空间导数或者欧拉导数。第二种情况称为實質導數或拉格朗日导数。
随質导数定义为算子:
其中是流体的速度。方程右边的第一项是普通的欧拉导数(也就是在静止参照系中的导数),而第二项表示由于流体的运动带来的变化。这个效应称为移流。
L在一个控制体积上守恆,可用以下积分式表示:
設Ω共动(隨流體移動),則它随着时间而改变,所以不能将时间导数和积分简单的交换,但有:
因为这个表达式对于所有成立,它可以简化为:
(第一個等號用到隨質導數的定義,以及對用到導數的積法則)对于不是密度的量(因而它不必在空间中积分),给出了正确的共动时间导数。
连续性方程
质量的守恒写作:
其中是流体的密度。
在不可压缩流体的情况,是常數,不取決於时间或位置,於是方程简化为:
方程组
特殊形式
这些是问题的特定的常见简化,有时解是已知的。
牛顿流体
在牛顿流体中,如下假设成立:
其中
- 是液体的粘滞度。
其中为简化书写,对脚标使用了爱因斯坦求和约定。
不采用简化书写的完整形式非常繁琐,分别为:
动量守恒:
质量守恒:
因为密度是一个未知数,我们需要另一个方程。
能量守恒:
其中:
假设一个理想气体:
上面是一共6个方程6个未知数的系统。(u,v,w,T,e以及 )。
不可壓缩流體
其纳维-斯托克斯方程(Navier-Stoke equation)分为動量守恆公式
-{zh-hant: ; zh-hans: ; }-
和質量守恆公式
- 。
其中,對不可壓縮牛頓流體來說,只有對流項(convective terms)為非線性形式。對流加速度(convective acceleration)來自於流體流動隨空間之變化所產生的速度改變,例如:當流體通過一個漸縮噴嘴(convergent nozzle)時,流體產生加速之情況。由於此項的存在,對於暫態運動中的流體來說,其流場速度變化不再單是時間的函數,亦與空間有關。
另外一個重要的觀察重點,在於黏滯力(viscosity)在流場中的以流體速度作拉普拉斯運算來表現。這暗示了在牛頓流體中,黏滯力為動量擴散(diffusion of momentum),與熱擴散方程式非常類似。
若在整个流体上均匀,动量方程简化为
如果现在再有为常数,我们得到如下系统:
连续性方程(假设不可压缩性):
- 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秒的模擬結果
参考文献
- McLean, Doug. . . 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...
- , claymath.org, Clay Mathematics Institute, March 27, 2017 [2017-04-02], (原始内容存档于2015-12-22) (英语)
- Fefferman, Charles L. (PDF). claymath.org. Clay Mathematics Institute. [2017-04-02]. (原始内容 (PDF)存档于2015-04-15) (英语).
- Fluid Mechanics(Schaum's Series), M. Potter, D.C. Wiggert, Schaum's Outlines, McGraw-Hill (USA), 2008, ISBN 978-0-07-148781-8
- Vectors, Tensors, and the basic Equations of Fluid Mechanics, R. Aris, Dover Publications, 1989, ISBN (10) 0-486-66110-5
- Parker, C. B. 2nd. 1994. ISBN 0-07-051400-3 (英语).
- 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
- Gorban, A.N.; Karlin, I. V. . Contemporary Physics (Review article). 2016, 58 (1): 70–90. Bibcode:2017ConPh..58...70G. S2CID 55317543. arXiv:1702.00831 . doi:10.1080/00107514.2016.1256123.
- Cercignani, C. . Friedlander, S.; Serre, D. (编). 1. Amsterdam: North-Holland. 2002: 1–70. ISBN 978-0444503305 (英语).
- Nie, X.B.; Chen, S.Y.; Robbins, M.O. . Journal of Fluid Mechanics (Research article). 2004, 500: 55–64 [2021-10-12]. Bibcode:2004JFM...500...55N. doi:10.1017/S0022112003007225. (原始内容存档于2018-08-09) (英语).
- Öttinger, H.C. . Berlin, Heidelberg: Springer Science & Business Media. 2012. ISBN 9783540583530. doi:10.1007/978-3-642-58290-5 (英语).
- Chen, Weijia; Chen, JC. . International Journal for Numerical Methods in Fluids. 2011-06-06, 68 (10). ISSN 0271-2091. doi:10.1002/fld.2602.
外部链接
- 克雷数学研究院纳维-斯托克斯方程大奖
- 纳维-斯托克斯方程的一个推导
- 纳维-斯托克斯方程的推导 (页面存档备份,存于)
- NASA关于纳维-斯托克斯方程的网页
- 纳维-斯托克斯方程(一些精确解) (页面存档备份,存于),位于EqWorld:数学方程的世界