爲什麼雨滴落下不會砸死人？《張朝陽的物理課》推導斯托克斯定律

爲什麼無論多高的雨滴落下來都有一個不大的末速度？普通物理課會告訴學生，因爲雨滴在下落時會受到與速度正相關的空氣阻力，導致雨滴的下落末速度趨於某個定值。但這個阻力與速度的關係是怎樣的呢？一百七十年前，愛爾蘭物理學家斯托克斯對類似的問題產生了興趣，並提出了在小雷諾數條件下流體阻力與速度成正比的斯托克斯定律。斯托克斯當時又是如何推導的呢？

10月6日14時，《張朝陽的物理課》第二百二十五期開播，搜狐創始人、董事局主席兼首席執行官、物理學博士張朝陽坐鎮搜狐視頻直播間，從流體力學最基本的納維爾-斯托克斯方程（Navier-Stokes equations，簡稱NS方程）出發，持續四個多小時“馬拉松”式地一步步推導出斯托克斯定律。

回顧雨滴下落速度的求解

假設有一個密度爲ρ_d的球形雨滴，它在密度爲ρ_a的空氣中從靜止開始下落，這個過程中雨滴會受到三個力：向下的重力G、向上的空氣浮力F_{浮}和向上的空氣阻力Fₛ，其中空氣阻力滿足斯托克斯定律，將它們依次寫出來是

注意看最後一個式子所描述的空氣阻力，它正是斯托克斯定律所給出的表達式，其中μ是空氣的粘滯係數，R是雨滴半徑，v是雨滴相對空氣的速度，這三個量的一次方相乘再乘以6π就得到了球形物體在不可壓縮定常流體中所受到的阻力。本節課的任務，正是從基本的流體力學最基本的NS方程出發，推導出斯托克斯定律。不過在進入正題前，先看看斯托克斯定律在雨滴下落情景中發揮的作用。

由剛剛分析出的三個力可以寫下雨滴的運動方程

爲了化簡方程，可以將係數用兩個字母代替

這樣就能把雨滴運動方程化簡爲

可以猜測這個微分方程的解具有如下形式

代回原方程，得到

再考慮到初始條件，即雨滴在零時刻靜止，可以得到

所以

和自由落體的速度線性增長不同的是，在考慮空氣阻力的情況下，雨滴的下落速度在t→∞時無限逼近一個有限值，它就是雨滴下落的最終速度

爲了更加直觀地感受空氣阻力的效果，張朝陽取了一些具體數據代入

算出雨滴下落的最終速度爲10m/s，和大家日常生活的經驗十分吻合，說明斯托克斯定律對空氣阻力的描述相當地好。

以上計算只能是針對0.3mm的小雨滴，或者說，它更接近雲層中的水霧。現實中，雲層中的水霧會匯聚成毫米量級的大水珠再落下來，水珠下落時受到空氣阻力還會發生形變，所以斯托克斯定律並不能完美地描述雨滴實際受到的阻力。事實上，此時占主導地位的阻力是一個與速度平方成正比的壓降阻力。不過在本堂課，我們更關心斯托克斯定律的導出，至於斯托克斯定律在雨滴下落問題的適用性，可以認爲，在雨滴很小的情況下，它總是近似成立的。

關於斯托克斯力並不是大部分情況下雨滴所受空氣阻力的論述，是課後由孫昌璞院士友情指正的，特此致謝。

尋根問底：從流體力學最基礎的NS方程出發

上一節求解了雨滴在受到與速度成正比的空氣阻力下的運動過程，那爲什麼空氣阻力可以具有與速度的線性關係？斯托克斯定律爲什麼能在小雨滴在空氣中下落的情形下成立？它又是否是對所有的流體導致的阻力都適用的呢？

帶着這樣的疑問，張朝陽開啓了本場“馬拉松”課程的正式內容。球形雨滴在空氣中下落，這個過程具有以下落方向爲軸的旋轉對稱性，而雨滴本身又是球對稱的，所以自然地可以想到要把這些拉普拉斯算子在柱座標或者球座標下寫開。張朝陽選擇先以下落方向爲極軸，建立了如下所示的球座標系。

建立球座標系

建立好座標系後，就可以寫出空氣的各點性質關於座標的函數，比如各個位置上的壓強p(r,θ,ϕ)和速度v(r,θ,ϕ)。物理學家通常用場的概念描述這種隨空間分佈的性質，比如壓強就是一個標量場，速度就是一個矢量場。一般的場還會帶有時間座標t，但這裡研究的是達到穩定狀態後的現象，所以可以忽略座標t。

流體力學中的壓強場和速度場，就好像質點力學中的力和速度。上一節課將雨滴看成一個質點，用牛頓第二定律寫下了它的運動方程。本節要研究的是，空氣作爲一個黏滯流體，它會給在其中做相對運動的雨滴一個什麼樣的阻力（有時也叫流體曳力，drag）。爲了求解空氣阻力，需要先求出空氣自己的運動狀態，而空氣作爲一個流體，它的運動方程可以由NS方程給出

NS方程本質上就是牛頓第二定律在流體上的應用，方程的左邊類似於ma，其中加速度包含速度場對時間的偏導項，和速度場對空間偏導後再對時間求偏導，方程的右邊類似於F，其中包含了壓強梯度所導致的正向壓力差，和流體運動的粘性所導致的粘滯力。

現在假設流體已經達到了穩態，所以NS方程左邊第一項，也就是關於時間求偏導的項爲0。而方程左邊的第二項是一個速度場關於空間分佈變化率的項，可以進一步假設流體微元在隨着流線運動的過程中速度的空間變化率是緩慢的，也就是近似認爲NS方程左邊第二項爲0。經過穩態和空間緩變的這樣兩個假設，NS方程被簡化爲了一個線性的微分方程

類比電動力學，巧妙引入渦度

方程的左邊是一階導，右邊是二階導，有沒有可能將右邊降階爲一階導呢？回憶矢量微積分中有這樣一條公式

而流體的質量守恆給出

這裡引入第三個假設：在速度比較小時空氣是不可壓縮的流體。這個假設在流速小於0.3倍音速的情形下通常是成立的。在不可壓縮假設下，空氣密度是一個常量，所以質量守恆導出速度場是無散的

將(1)-(3)式聯立，得到

這個等式的右邊看起來還是二階導，但與(1)式不同的是，這裡的nabla算子▽是依次以叉乘的形式作用在後面的矢量上的，而(1)式是兩個nabla算子以點乘成拉普拉斯算子的形式作用到速度矢量上，前者的兩次求導操作是容易拆分的，後者要拆分的話比較困難，需要先作用一次導出二階張量再求散度來縮並回一階矢量。受到(4)式的啓發，可以定義一箇中間變量

它是速度場的旋度，稱爲渦度(vorticity)。這樣就能把(1)式右邊也降階爲一階導的形式

這裡巧妙地引入了渦度場來代替速度場，從而讓NS方程兩邊的求導階數相等。而巧合的是，電動力學中也曾有類似的操作。張朝陽將電動力學的量類比成三層樓，第一層是電勢Φ和磁矢勢A，第二層是電場E和磁場B，第三層是電荷密度ρ和電流密度j。每上升一層就對應求一次導，例如

比較(5)和(7)，不難看出對於同樣是無散的速度場和磁場，都可以利用它們的旋度，也就是渦度和電流密度，來作爲中間變量進而簡化方程。

類比電動力學的“三層樓”

經過三個假設後，NS方程簡化爲了(6)式，它是一個關於壓強場和渦度場的方程，有沒有可能將這兩個場分開呢？首先，注意到任何一個矢量場的旋度無散，即

所以對(6)式兩邊求散度，得到只關於壓強場的泊松方程

單獨一個(8)式只給出了壓強場的信息，顯然它損失了渦度場的信息。張朝陽打比方道，這就像是有兩個數3和4，將它們乘在一起得到了12，但單獨一個12是還原不出3和4的，它也可能對應2和6。類似的道理，對方程(6)求散度得到(8)式是會損失信息的，還需要再從方程(6)導出渦度場的信息。同樣地，注意到任意一個標量場的梯度無旋，即

所以對(6)式兩邊求旋度，得到只關於渦度場的方程

這裡再次利用了矢量微積分的公式(2)。另外，渦度是速度場的旋度，所以它是無散的。因此渦度場也滿足泊松方程

靈活切換座標系，轉化渦度場的矢量泊松方程

上一節將NS方程簡化成了兩個泊松方程，一個是關於壓強標量的式(8)，另一個是關於渦度矢量的式(9)。標量的泊松方程在之前的課上有相當多的處理經驗，但矢量的泊松方程是之前從沒遇到過的，爲此需要先研究一下矢量被拉普拉斯算符作用後得到的是什麼。

根據流體運動的軸向對稱性，速度矢量應該只有徑向分量v_r和極角方向的分量v_θ，而沒有方位角方向上的分量v_ϕ，並且各方向的速度分量對ϕ的求導都爲0。基於對速度場的軸對稱性的認知，可以在球座標下寫出渦度

可見在流體速度有軸對稱性的情況下，渦度是一個只有方位角方向分量的矢量。由此(9)式變成了

這是一個矢量場的泊松方程。首先來計算第一次nabla算符作用後的結果，它將被作用的矢量沿不同方向求導，但對求導方向的基矢和被作用後的矢量的基矢這兩個基矢而言做了張量積，張量積既不是點乘也不是叉乘，而是把兩個基矢直接放在一起作爲二階張量的基底，以三維空間來看，它包含了3×3=9個係數和基底。用⊗代表矢量的張量積，可以寫成

(12)式的兩項分別將nabla算符作用在了係數和基矢上，它的第一項再被nabla算符點乘後是

上式第二項中被大括號標出的部分爲0，因爲球座標的散度公式爲

而基矢\vec{e}_ϕ就相當於g_r=g_θ=0，g_ϕ=1的一個矢量，代入散度公式可知它等於0。

(12)式的第二項涉及到直接對一個矢量求“梯度”得到二階張量，展開來寫是

第二個等式新定義了一個矢量\vec{e}_ρ，它是從\vec{e}_ϕ對ϕ求導出來的。方位角基矢與r和θ無關，只會隨着ϕ變化，且變化的矢量是指向極軸的一個向內的矢量，類似於柱座標下的徑向矢量。

柱座標的切面

爲了方便對這個二階張量求散度，接下來把它切換到直角座標系下，將\vec{e}_ϕ和\vec{e}_ρ用\vec{i}和\vec{j}展開

直角座標系下的座標和球座標之間有這樣的關係

因此，

或者更顯式地，把它寫成一個2×2的矩陣形式（因爲ω沒有z方向分量，可以把問題簡化在xy平面上）

至此，(12)式的第二項化爲了直角座標系下的形式，無論從張量積的形式看還是從矩陣的形式看，它都是一個二階張量，可以用字母F上加兩個箭頭來表示這是一個二階張量。對這個二階張量再求一次散度

最後一步把第一項又寫回到球座標的形式，第二項則引入了矢量

它是從極軸出發指向空間中某個位置的位矢。注意到f是一個和ϕ無關的函數，雖然它的分子上有ω_ϕ，但從(10)式可以看出ω_ϕ是和ϕ無關的，這一切都來源於流速場具有繞軸旋轉的對稱性，因此▽f在xy平面上的分量只能是平行於ρ的，所以(14)式中被大括號包住的部分爲0。

將(13)和(14)式的結果代入(11)式的矢量泊松方程，可以將基底\vec{e}_ϕ提出到拉普拉斯算符外面，從而得到只和係數ω_ϕ有關的泊松方程

巧借氫原子球諧函數，求解球座標下的壓強場和渦度場

將壓強場的泊松方程(8)式和渦度係數的方程(15)寫在球座標下，並且注意到它們在ϕ方向上的對稱性，可以得到

壓強場的方程比較簡單，先求解它。假設壓強場方程的解能分離成與徑向有關的部分以及與極角方向有關的部分

代入(16.a)，得到

上式第一個部分只與f(r)有關，第二個部分只與g(θ)有關，所以它們應該都等於一個常數，注意到這個方程和《張朝陽的物理課（第一卷）》中氫原子薛定諤方程的形式一模一樣，不妨就借鑑那裡的思路，令這個常數等於角量子數l(l+1)。對於徑向部分，取(17)式的第一部分和第三部分做比奈（Binet）變換

令ψ=rf，上式可化簡爲

這個方程有兩個解系，一個是ψ~r^(l+1)，一個是ψ~r^(-l)。

對於極角方向部分，取(17)式的第二部分和第三部分

令x=cos(θ),h(x)=g(θ)

這個方程的解正是勒讓德多項式

結合徑向和角向的結論，可以寫下壓強場的通解

剩下的問題是確定展開係數A和B。爲此，張朝陽先從量綱的角度分析了它與哪些物理量有關，對於流體產生的阻力，不難想象它和在其中運動的球狀物體的半徑R、運動速度v₀、以及流體粘滯係數μ有關

一個樸素的想法是，先將這些量的一次方相乘起來，看看量綱是什麼

和壓強的量綱正好差L²，因此可以大膽地假設，(19)式中只包含r^(-2)的項。當然這裡會有一個小問題，R/r是一個無量綱量，完全可以在表達式中再乘上這一無量綱量來維持量綱的平衡。針對這個不確定性，數學家們已經證明了，在給定足夠的邊界條件下，微分方程的解具有唯一性，因此只要求出了其中一個滿足邊界條件的解，就可以認爲得到了完整的解。

經過這樣的量綱分析，不難發現只有l=1的情況才滿足假設

再來看(16.b)的渦度場方程，徑向部分的分析和壓強場是一樣的，而角向部分相比(18)式會多出一項

多出來的第三項，正好是球諧函數Y在m=1時方位角分量e^(imϕ)被拉普拉斯算子作用後出現的項

由此可以確定，(21)式的解是連帶勒讓德多項式

所以渦度場的通解爲

渦度是速度的散度，它具有量綱

如果和壓強場一樣，也樸素地假設它含有R的一次項和v₀的一次項，那麼可以推斷(22)式只含有r^(-2)的項，也就是l=1。由於

因此

降一層臺階：用斯托克斯流函數代入邊界條件

上一節利用之前的物理課解氫原子薛定諤方程的經驗，成功得出了壓強場(20)和渦度場(23)，但它們各自都有一個待定係數，需要用邊界條件來定出。然而渦度的邊界條件是難以給定的，因爲根據(10)式，它和速度的兩個分量vr、vθ的導數都有關係。回憶引入渦度時，出發點是爲了將(1)式右邊的二階導變成一階導，相當於視角升了一個臺階，從而降低微分的階數。那麼有沒有可能再轉化一次視角，把vr和vθ統一成另一箇中間量，以便於引入邊界條件呢？

答案是可行的。在引入渦度時，曾將速度場類比成磁場，因爲它們都具有無散的性質，可以通過(2)式將拉普拉斯算子轉化爲兩個散度的依次作用。無散的矢量場還有一個性質，就是它們總可以寫成另一個矢量場的旋度，正是這一點允許從磁場定義一個磁矢勢。同理，可以定義一個矢量場\vec{Ψ}，它的旋度等於速度場，從而將視角降一個臺階

張朝陽類比電動力學引入斯托克斯流函數

在球座標下寫開各分量

和一般的在球座標下求散度的情形不同的是，最後一步重新定義了Ψᵩ，把rsin(θ)吸收了進去，並簡記爲Ψ，稱爲斯托克斯流函數。再次利用軸對稱的性質，可以知道v沒有ϕ方向的分量，且Ψ是一個與ϕ無關的函數，所以

也就是

到這一步，就把兩個速度分量統一成了一個斯托克斯流函數。將它們代入(10)式

再代入(23)式，就可以從渦度的場方程轉換到斯托克斯流函數的方程

這個偏微分方程相比(16)式要複雜得多，它很難直接地分離出徑向部分和角向部分，但可以先假設它是可分量變量的，從而猜測它應該滿足什麼樣的形式，一個最簡單的猜解是

代回檢驗，不難發現這個形式能非常方便地配平(25)式中有關角向的部分的冪次

由此得到徑向部分的方程

這是一個變係數非齊次線性微分方程，可以確定它的通解是這些冪函數的和

代回原式可以定出a=D₁/2。

至此，終於可以考慮速度的邊界條件了。在本問題中，有兩個邊界條件需要考慮，一個是無窮遠處的流體應當沒有受到球狀雨滴的影響，所以保持靜止；另一個是雨滴表面的流體相對雨滴靜止，稱爲no-slip condtion，這在小雷諾數情形下是普遍成立的。把(26)代入(24)，並要求無窮遠處速度趨於0，可以定出b=0。再要求球面處有\vec{v}(R)=\vec{v}₀，對應r分量和θ分量分別是

由此可以定出

所有的待定係數都已經得到了確定，最終結果是

在用邊界條件定下速度場和渦度場後，就可以聯立(6)和(20)來確定壓強場的待定係數B₁了。

比較上下兩式，可以定出

所以壓強場爲

當然，空氣應該還有個大氣壓強p₀作爲背景壓強，但加上一個常數並不影響壓強場的梯度，因此這裡可以忽略大氣壓。

計算流體作用在球狀物體上的斯托克斯阻力

在本文最開始介紹完整的NS方程時，就說道等號右側代表流體受力的兩項分別是壓強梯度導致的正向壓力差和流體粘性導致的粘滯力。流體對雨滴的作用力也一樣有這兩個來源，所以對於流體的應力張量，它的表達式是

最終作用在球體上的力，就是這個二階的應力張量與球面的法向量\vec{e}ᵣ點乘得到的一階矢量在整個球面上的積分。正是由於二階張量有非對角項，點乘出的力並不一定垂直於法向，而是會有一定的切向分量。

(應力張量與面法向量的點乘)

對於(28)式的第一項，將張量基底寫開並與球面法向點乘是

根據軸對稱性，壓強在球面上的積分只剩下沿極軸z方向的分量，其他方向的分量會互相抵消，所以可以只對壓強的z向分量做積分

同理，(28)式的第二項與球面法向點乘是

上式最後一步代入r=R時，r方向上的速度梯度剛好是一個極值點，所以不貢獻粘滯力，只有θ方向貢獻出來一個切向的壓強。將這個壓強的z分量在球面上積分

(28)式的第三項與球面法向點乘是

這說明這一項在球面不會貢獻壓強。其中第二個等號運用了球座標單位基矢的求導關係

將以上三項加起來，就得到了斯托克斯定律

回顧整堂課的推導過程，雖然要解決的是流體力學問題，但其中貫穿了各個領域的知識，比如從電動力學的矢量微積分引入了渦度、速度、斯托克斯流函數的“三層樓”，借用量子力學中的球諧函數解出了渦度角向分量的泊松方程。當然，從歷史進程來看，數學家們先是在研究流體力學時發展了這些工具，再到後來發展電動力學和量子力學時，人們就可以直接使用這些已經發展好的工具。而《張朝陽的物理課》研習路徑則剛好反了過來，是先學習了電動力學和量子力學，再來研究這個大部分普通物理課上一筆帶過的斯托克斯定律。張朝陽以此提示大家，學科的發展本就沒有領域的限制，只有將一個領域的內容學紮實了，才能在面對新問題時靈活運用先前的知識。

本堂課也是在向斯托克斯致敬。斯托克斯在劍橋執教期間對流體力學領域做出很多奠基性的貢獻，並早在1851年就得出了今天所講的斯托克斯定律。他的研究對後世產生了深遠的影響，比如矢量分析的工具幫助了麥克斯韋在1865年建立起電磁方程組，密立根油滴實驗用到了斯托克斯定律來確定油滴的質量，風洞試驗、汽車風阻模擬、機翼設計等工程問題也處處有NS方程的身影。NS方程的解的存在性與光滑性至今仍是數學界的研究熱點，並在21世紀初被列爲千禧年七大數學難題之一。儘管現在能借助計算機數值求解NS方程，但尋找解析解依舊是幫助人們理解物理圖像的重要手段之一。