納維-斯托克斯方程

納維爾－斯托克斯方程，以法國工程師兼物理學家克勞德-路易·納維、愛爾蘭物理學和數學家喬治·斯托克斯兩人命名，是一組偏微分方程，描述液體和空氣等流體的運動。納維爾－斯托克斯方程表達了牛頓流體運動時，動量和質量守恆。有時，還連同狀態方程列出，說明流體壓強、溫度、密度之間的關係。方程斷言，流體粒子動量的改變率，來自作用在液體內部的壓力變化、耗散粘滯力、以及重力。其中粘滯力類似於摩擦力，產生於分子的相互作用，越黏的流體，該作用就越強。這樣，納維-斯托克斯方程描述作用於液體任意給定區域的力的動態平衡。學術研究和經濟生活中，許多重要物理過程都可用納維爾－斯托克斯方程描述，因此該些方程有很重要的研究價值。它們可以用於模擬天氣、洋流、管道中的水流、星系中恆星的運動、翼型周圍的氣流，也可以用於設計飛行器和車輛、研究血液迴圈、設計電站、分析汙染效應等等。納－斯方桯組與馬克士威方程組聯立，用於研究磁流體力學。納維－斯托克斯方程依賴微分方程來描述流體的運動。不同於代數方程，這些方程不尋求建立所研究的變數與內部壓力的導數成正比。這表示對於給定的物理問題，至少要用微積分才可以求得其納維-斯托克斯方程的解。實用上，也只有最簡單的情況才能用這種方法獲得已知解。這些情況通常涉及穩定態。對於更復雜的情形，例如厄爾尼諾這樣的全球性氣象系統或機翼的升力，現時僅能借助計算機求出納維-斯托克斯方程的數值解。這個科學領域稱為計算流體力學。雖然紊流是日常經驗中就可以遇到的，但這類非線性問題在理論上極難求解，仍未能證明三維空間中是否總存在光滑解，甚至有界解。此問題稱為納維-斯托克斯存在性與光滑性。克雷數學學院於2000年5月21日列入七大未解難題，懸賞一百萬美元，獎勵找到數學證明或反例的任何人。性質.非線性.納維－斯托克斯方程的一般形式是非線性的偏微分方程，所以在大多數實際情況下仍是如此。在特定情況，如一維流和斯托克斯流下，方程可以簡化為線性方程組。非線性項的出現，令大部分問題變得很難，甚至無法求解。另一方面，方程組能描述湍流，而非線性正是湍流出現的重要因素。方程式中的非線性項是對流加速，流經一個細小並逐漸收窄的。此種流，不論能否確切解出，通常都能透徹研究、理解。湍流.湍流是時變的混沌行為，這種行為常見於許多流體流動中。人們普遍認為，湍流的成因，是整個流體的慣性：時變加速度與對流加速度疊加，以產生亂流告終。因此慣性影響很小的流體，往往是層流。雖然不完全確定，但一般相信納維－斯托克斯方程能夠合理地描述湍流。納維-斯托克斯方程關於湍流的數值解是非常難得到的，而且由於湍流之中，有多個顯著不同的尺度，若要得到穩定解，所需要的解析度要極度精細，於是計算或的時間長得不可行。若試圖用解層流的方法來解決湍流問題，通常會得到時間不穩定的解，而不能適當收斂。為了解決這個問題，計算流體力學中，實際模擬湍流的程式，多採用雷諾平均納維－斯托克斯方程。這種方法比RANS方法，佔用更多計算時間和記憶體空間，但效果較好，因為LES明確解出較大的湍流尺度。適用性.連同其他方程也符合實際觀察結果。納維－斯托克斯方程假定所研究的流體連續，且不具相對論流速。在非常小的尺度或極端條件下，由離散分子組成的真實流體，與由納維－斯托克斯方程描繪的連續流體，將產生不同的結果。例如，大梯度流的流體內層，有毛細現象。對於大克努森數的問題，用統計力學的波茲曼方程式可能更適合 ，甚至要用分子動力學或其他混合方法。另一個限制是方程的複雜性。要刻劃一般常見的流體類，有行之有效的公式，但對於較罕見的類別，應用納維－斯托克斯方程時，往往會得到非常複雜的描述，甚至是未解難題。出於這個原因，這些方程通常用於描述牛頓流體。研究這種液體是“簡單”的，因為粘度模型最終被線性化；真正描述其他型別流體的普遍模型，截至2012年還不存在。基本假設.在解釋納維－斯托克斯方程的細節之前，首先對流體的性質作幾個假設。第一個假設是流體是連續的。這強調其內部沒有空隙，例如，汽水中有氣泡，便不屬此例；同樣，包含霧狀粒子的氣體亦不屬此例。另一個必要的假設是，所有涉及到的場，即壓強、速度、密度、溫度等，全部都可微。該方程從質量守恆、動量守恆、能量守恆等基本原理匯出。推導過程中，經常考慮一個有限的任意體積，稱為控制體積，並對其應用這些原理。該有限體積記為formula_1，而其表面記為formula_2。該控制體積可以在空間中固定，也可能隨著流體運動。這會導致一些特殊的結果，見下文。隨質導數.運動流體的屬性變化，譬如大氣風速的變化，可以有兩種不同的方法來測量：風速儀可以固定在氣象站上，也可以隨氣象氣球飄動。第一種情況下，風速儀測量的速度，是所有運動的粒子經過一個固定點的速度；而第二種情況下，儀器測量的數值，是其隨著流體運動時，速度的變化。同樣的論證對於密度、溫度等物理量的測量也是成立的。因此，當作微分時必須區分兩種情況。第一種情況稱為空間導數或者尤拉導數。第二種情況稱為實質導數或拉格朗日導數。隨質導數定義為運算元：formula_3其中formula_4是流體的速度。方程右邊的第一項是普通的尤拉導數，而第二項表示由於流體的運動帶來的變化。這個效應稱為。L在一個控制體積formula_1上守恆，可用以下積分式表示：formula_6設Ω共動，則它隨著時間而改變，所以不能將時間導數和積分簡單的交換，但有：formula_7因為這個表示式對於所有formula_1成立，它可以簡化為：formula_9，formula_11給出了正確的共動時間導數。守恆定律.在上式中，依次取formula_12為下列各守恆量：並且，就能推匯出NS方程。連續性方程.質量的守恆寫作：formula_13其中formula_14是流體的密度。在不可壓縮流體的情況，formula_14是常數，不取決於時間或位置，於是方程簡化為：formula_16動量守恆.動量守恆寫作：formula_17注意formula_18是一個張量，formula_19代表張量積。可以進一步簡化，利用連續性方程，上式變成：formula_20可以認出這就是通常牛頓第二定律formula_21。方程組.一般形式.方程組的形式.納維-斯托克斯方程的一般形式是：formula_22關於動量守恆。張量formula_23代表施加在一個流體粒子上的表面力。除非流體是由象旋渦這樣的旋轉自由度組成，formula_23是一個對稱張量。一般來講，我們有如下形式：formula_25其中formula_26是法向約束，而formula_27是切向約束。跡formula_28在流體處於平衡態時為0。這等價於流體粒子上的法向力的積分為0。我們再加上連續性方程：formula_29對於"處於平衡"的液體，formula_23的跡是3"p"。其中"p"是壓強最後，我們得到：formula_31其中formula_32是formula_23的非對角線部分。閉合問題.這些方程是不完整的。要對它們進行完備化，必須對formula_23的形式作一些假設。例如在理想流體的情況formula_27分量為0。用於完備方程組的方程是狀態方程。再如，壓強可以主要是密度和溫度的函式。要求解的變數是速度的各個分量，流體密度，靜壓力，和溫度。流場假定為可微並連續，使得這些平衡得以用偏微分方程表達。這些方程可以轉化為渦度和流函式這些次變數的威爾金森方程組。解依賴於流體的性質，並且依賴於所研究的區域的邊界條件。formula_23的分量是流體的一個無窮小元上面的約束。它們代表垂直和剪下約束。formula_23是對稱的，除非存在非零的自旋密度。所謂非牛頓流體就是其中該張量沒有特殊性質使得方程的特殊解出現的流體特殊形式.這些是問題的特定的常見簡化，有時解是已知的。牛頓流體.在牛頓流體中，如下假設成立：formula_38其中formula_39是液體的粘滯度。formula_40formula_41其中為簡化書寫，對腳標使用了愛因斯坦求和約定。不採用簡化書寫的完整形式非常繁瑣，分別為：動量守恆：formula_42formula_43formula_44質量守恆：formula_45因為"密度"是一個未知數，我們需要另一個方程。能量守恆：formula_46其中：formula_47假設一個理想氣體：formula_48上面是一共6個方程6個未知數的系統。。賓漢流體.在賓漢流體中，我們有稍微不同的假設：formula_50那些流體在開始流動之前能夠承受一定的剪下。牙膏是一個例子。冪律流體.這是一種理想化的流體，其剪下應力，formula_27，由下式給出formula_52該形式對於模擬各種一般流體有用。不可壓縮流體.其納維－斯托克斯方程分為動量守恆公式-{zh-hant:formula_53;zh-hans:formula_54;和質量守恆公式formula_55。其中，對不可壓縮牛頓流體來說，只有對流項時，流體產生加速之情況。由於此項的存在，對於暫態運動中的流體來說，其流場速度變化不再單是時間的函式，亦與空間有關。另外一個重要的觀察重點，在於黏滯力，與熱擴散方程式非常類似。formula_56;formula_57是散度，formula_58是克羅內克記號。若formula_39在整個流體上均勻，動量方程簡化為formula_60。如果現在再有formula_14為常數，我們得到如下系統：formula_63formula_64formula_65連續性方程：formula_66N-S方程的簡化版本。採用《不可壓縮流》，Ronald Panton所著第二版注意納維－斯托克斯方程僅可近似描述液體流，而且在非常小的尺度或極端條件下，由離散的分子和其他物質的混合體組成的真實流體，會產生和納維－斯托克斯方程所描述的連續並且齊性的液體不同的結果。依賴於問題的納森數，統計力學可能是一個更合適的方法。但是，納維－斯托克斯方程對於很大範圍的實際問題是有效的，只要記住他們的缺陷是天生的就可以了。程式模擬.參考MIT18086_NAVIERSTOKESfunction 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 numberdt = 1e-2; % time steptf = 4e-0; % final timelx = 1; % width of boxly = 1; % height of boxnx = 90; % number of x-gridpointsny = 90; % number of y-gridpointsnsteps = 10; % number of steps with graphic outputnt = 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 conditionsU = zeros(nx-1,ny); V = zeros(nx,ny-1);% boundary conditionsuN = 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 termsgamma = 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 viscosityrhs = 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 correctionrhs = 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;% visualizationif floor(25*k/nt)&gt;floor(25*(k-1)/nt), fprintf('.'), endif k==1|floor(nsteps*k/nt)&gt;floor(nsteps*(k-1)/nt)% stream functionrhs = 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 oncontour(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))drawnowendendfprintf('\n')function B = avg(A,k)if nargin&lt;2, k = 1; endif size(A,1)==1, A = A'; endif k&lt;2, B = (A(2:end,:)+A(1:end-1,:))/2; else, B = avg(A,k-1); endif size(A,2)==1, B = B'; endfunction 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秒的模擬結果