【ITBEAR】
為什么無論多高的雨滴落下來都有一個(gè)不大的末速度？普通物理課會告訴學(xué)生，因?yàn)橛甑卧谙侣鋾r(shí)會受到與速度正相關(guān)的空氣阻力，導(dǎo)致雨滴的下落末速度趨于某個(gè)定值。但這個(gè)阻力與速度的關(guān)系是怎樣的呢？一百七十年前，愛爾蘭物理學(xué)家斯托克斯對類似的問題產(chǎn)生了興趣，并提出了在小雷諾數(shù)條件下流體阻力與速度成正比的斯托克斯定律。斯托克斯當(dāng)時(shí)又是如何推導(dǎo)的呢？
回顧雨滴下落速度的求解
假設(shè)有一個(gè)密度為ρ_d的球形雨滴，它在密度為ρ_a的空氣中從靜止開始下落，這個(gè)過程中雨滴會受到三個(gè)力：向下的重力G、向上的空氣浮力F_{浮}和向上的空氣阻力F?，其中空氣阻力滿足斯托克斯定律，將它們依次寫出來是
注意看最后一個(gè)式子所描述的空氣阻力，它正是斯托克斯定律所給出的表達(dá)式，其中μ是空氣的粘滯系數(shù)，R是雨滴半徑，v是雨滴相對空氣的速度，這三個(gè)量的一次方相乘再乘以6π就得到了球形物體在不可壓縮定常流體中所受到的阻力。本節(jié)課的任務(wù)，正是從基本的流體力學(xué)最基本的NS方程出發(fā)，推導(dǎo)出斯托克斯定律。不過在進(jìn)入正題前，先看看斯托克斯定律在雨滴下落情景中發(fā)揮的作用。
由剛剛分析出的三個(gè)力可以寫下雨滴的運(yùn)動方程
為了化簡方程，可以將系數(shù)用兩個(gè)字母代替
這樣就能把雨滴運(yùn)動方程化簡為
可以猜測這個(gè)微分方程的解具有如下形式
代回原方程，得到
再考慮到初始條件，即雨滴在零時(shí)刻靜止，可以得到
所以
和自由落體的速度線性增長不同的是，在考慮空氣阻力的情況下，雨滴的下落速度在t→∞時(shí)無限逼近一個(gè)有限值，它就是雨滴下落的最終速度
為了更加直觀地感受空氣阻力的效果，張朝陽取了一些具體數(shù)據(jù)代入
算出雨滴下落的最終速度為10m/s，和大家日常生活的經(jīng)驗(yàn)十分吻合，說明斯托克斯定律對空氣阻力的描述相當(dāng)?shù)睾谩?/p>
以上計(jì)算只能是針對0.3mm的小雨滴，或者說，它更接近云層中的水霧。現(xiàn)實(shí)中，云層中的水霧會匯聚成毫米量級的大水珠再落下來，水珠下落時(shí)受到空氣阻力還會發(fā)生形變，所以斯托克斯定律并不能完美地描述雨滴實(shí)際受到的阻力。事實(shí)上，此時(shí)占主導(dǎo)地位的阻力是一個(gè)與速度平方成正比的壓降阻力。不過在本堂課，我們更關(guān)心斯托克斯定律的導(dǎo)出，至于斯托克斯定律在雨滴下落問題的適用性，可以認(rèn)為，在雨滴很小的情況下，它總是近似成立的。
關(guān)于斯托克斯力并不是大部分情況下雨滴所受空氣阻力的論述，是課后由孫昌璞院士友情指正的，特此致謝。
尋根問底：從流體力學(xué)最基礎(chǔ)的NS方程出發(fā)
上一節(jié)求解了雨滴在受到與速度成正比的空氣阻力下的運(yùn)動過程，那為什么空氣阻力可以具有與速度的線性關(guān)系？斯托克斯定律為什么能在小雨滴在空氣中下落的情形下成立？它又是否是對所有的流體導(dǎo)致的阻力都適用的呢？
帶著這樣的疑問，張朝陽開啟了本場“馬拉松”課程的正式內(nèi)容。球形雨滴在空氣中下落，這個(gè)過程具有以下落方向?yàn)檩S的旋轉(zhuǎn)對稱性，而雨滴本身又是球?qū)ΨQ的，所以自然地可以想到要把這些拉普拉斯算子在柱坐標(biāo)或者球坐標(biāo)下寫開。張朝陽選擇先以下落方向?yàn)闃O軸，建立了如下所示的球坐標(biāo)系。
建立球坐標(biāo)系
建立好坐標(biāo)系后，就可以寫出空氣的各點(diǎn)性質(zhì)關(guān)于坐標(biāo)的函數(shù)，比如各個(gè)位置上的壓強(qiáng)p(r,θ,?)和速度v(r,θ,?)。物理學(xué)家通常用場的概念描述這種隨空間分布的性質(zhì)，比如壓強(qiáng)就是一個(gè)標(biāo)量場，速度就是一個(gè)矢量場。一般的場還會帶有時(shí)間坐標(biāo)t，但這里研究的是達(dá)到穩(wěn)定狀態(tài)后的現(xiàn)象，所以可以忽略坐標(biāo)t。
NS方程本質(zhì)上就是牛頓第二定律在流體上的應(yīng)用，方程的左邊類似于ma，其中加速度包含速度場對時(shí)間的偏導(dǎo)項(xiàng)，和速度場對空間偏導(dǎo)后再對時(shí)間求偏導(dǎo)，方程的右邊類似于F，其中包含了壓強(qiáng)梯度所導(dǎo)致的正向壓力差，和流體運(yùn)動的粘性所導(dǎo)致的粘滯力。
現(xiàn)在假設(shè)流體已經(jīng)達(dá)到了穩(wěn)態(tài)，所以NS方程左邊第一項(xiàng)，也就是關(guān)于時(shí)間求偏導(dǎo)的項(xiàng)為0。而方程左邊的第二項(xiàng)是一個(gè)速度場關(guān)于空間分布變化率的項(xiàng)，可以進(jìn)一步假設(shè)流體微元在隨著流線運(yùn)動的過程中速度的空間變化率是緩慢的，也就是近似認(rèn)為NS方程左邊第二項(xiàng)為0。經(jīng)過穩(wěn)態(tài)和空間緩變的這樣兩個(gè)假設(shè)，NS方程被簡化為了一個(gè)線性的微分方程
類比電動力學(xué)，巧妙引入渦度
方程的左邊是一階導(dǎo)，右邊是二階導(dǎo)，有沒有可能將右邊降階為一階導(dǎo)呢？回憶矢量微積分中有這樣一條公式
而流體的質(zhì)量守恒給出
這里引入第三個(gè)假設(shè)：在速度比較小時(shí)空氣是不可壓縮的流體。這個(gè)假設(shè)在流速小于0.3倍音速的情形下通常是成立的。在不可壓縮假設(shè)下，空氣密度是一個(gè)常量，所以質(zhì)量守恒導(dǎo)出速度場是無散的
將(1)-(3)式聯(lián)立，得到
這個(gè)等式的右邊看起來還是二階導(dǎo)，但與(1)式不同的是，這里的nabla算子▽是依次以叉乘的形式作用在后面的矢量上的，而(1)式是兩個(gè)nabla算子以點(diǎn)乘成拉普拉斯算子的形式作用到速度矢量上，前者的兩次求導(dǎo)操作是容易拆分的，后者要拆分的話比較困難，需要先作用一次導(dǎo)出二階張量再求散度來縮并回一階矢量。受到(4)式的啟發(fā)，可以定義一個(gè)中間變量
它是速度場的旋度，稱為渦度(vorticity)。這樣就能把(1)式右邊也降階為一階導(dǎo)的形式
這里巧妙地引入了渦度場來代替速度場，從而讓NS方程兩邊的求導(dǎo)階數(shù)相等。而巧合的是，電動力學(xué)中也曾有類似的操作。張朝陽將電動力學(xué)的量類比成三層樓，第一層是電勢Φ和磁矢勢A，第二層是電場E和磁場B，第三層是電荷密度ρ和電流密度j。每上升一層就對應(yīng)求一次導(dǎo)，例如
比較(5)和(7)，不難看出對于同樣是無散的速度場和磁場，都可以利用它們的旋度，也就是渦度和電流密度，來作為中間變量進(jìn)而簡化方程。
類比電動力學(xué)的“三層樓”
經(jīng)過三個(gè)假設(shè)后，NS方程簡化為了(6)式，它是一個(gè)關(guān)于壓強(qiáng)場和渦度場的方程，有沒有可能將這兩個(gè)場分開呢？首先，注意到任何一個(gè)矢量場的旋度無散，即
所以對(6)式兩邊求散度，得到只關(guān)于壓強(qiáng)場的泊松方程
單獨(dú)一個(gè)(8)式只給出了壓強(qiáng)場的信息，顯然它損失了渦度場的信息。張朝陽打比方道，這就像是有兩個(gè)數(shù)3和4，將它們乘在一起得到了12，但單獨(dú)一個(gè)12是還原不出3和4的，它也可能對應(yīng)2和6。類似的道理，對方程(6)求散度得到(8)式是會損失信息的，還需要再從方程(6)導(dǎo)出渦度場的信息。同樣地，注意到任意一個(gè)標(biāo)量場的梯度無旋，即
所以對(6)式兩邊求旋度，得到只關(guān)于渦度場的方程
這里再次利用了矢量微積分的公式(2)。另外，渦度是速度場的旋度，所以它是無散的。因此渦度場也滿足泊松方程
靈活切換坐標(biāo)系，轉(zhuǎn)化渦度場的矢量泊松方程
上一節(jié)將NS方程簡化成了兩個(gè)泊松方程，一個(gè)是關(guān)于壓強(qiáng)標(biāo)量的式(8)，另一個(gè)是關(guān)于渦度矢量的式(9)。標(biāo)量的泊松方程在之前的課上有相當(dāng)多的處理經(jīng)驗(yàn)，但矢量的泊松方程是之前從沒遇到過的，為此需要先研究一下矢量被拉普拉斯算符作用后得到的是什么。
根據(jù)流體運(yùn)動的軸向?qū)ΨQ性，速度矢量應(yīng)該只有徑向分量v_r和極角方向的分量v_θ，而沒有方位角方向上的分量v_?，并且各方向的速度分量對?的求導(dǎo)都為0。基于對速度場的軸對稱性的認(rèn)知，可以在球坐標(biāo)下寫出渦度
可見在流體速度有軸對稱性的情況下，渦度是一個(gè)只有方位角方向分量的矢量。由此(9)式變成了
這是一個(gè)矢量場的泊松方程。首先來計(jì)算第一次nabla算符作用后的結(jié)果，它將被作用的矢量沿不同方向求導(dǎo)，但對求導(dǎo)方向的基矢和被作用后的矢量的基矢這兩個(gè)基矢而言做了張量積，張量積既不是點(diǎn)乘也不是叉乘，而是把兩個(gè)基矢直接放在一起作為二階張量的基底，以三維空間來看，它包含了3×3=9個(gè)系數(shù)和基底。用?代表矢量的張量積，可以寫成
(12)式的兩項(xiàng)分別將nabla算符作用在了系數(shù)和基矢上，它的第一項(xiàng)再被nabla算符點(diǎn)乘后是
上式第二項(xiàng)中被大括號標(biāo)出的部分為0，因?yàn)榍蜃鴺?biāo)的散度公式為
而基矢vec{e}_?就相當(dāng)于g_r=g_θ=0，g_?=1的一個(gè)矢量，代入散度公式可知它等于0。
(12)式的第二項(xiàng)涉及到直接對一個(gè)矢量求“梯度”得到二階張量，展開來寫是
第二個(gè)等式新定義了一個(gè)矢量vec{e}_ρ，它是從vec{e}_?對?求導(dǎo)出來的。方位角基矢與r和θ無關(guān)，只會隨著?變化，且變化的矢量是指向極軸的一個(gè)向內(nèi)的矢量，類似于柱坐標(biāo)下的徑向矢量。
柱坐標(biāo)的切面
為了方便對這個(gè)二階張量求散度，接下來把它切換到直角坐標(biāo)系下，將vec{e}_?和vec{e}_ρ用vec{i}和vec{j}展開
直角坐標(biāo)系下的坐標(biāo)和球坐標(biāo)之間有這樣的關(guān)系
因此，
至此，(12)式的第二項(xiàng)化為了直角坐標(biāo)系下的形式，無論從張量積的形式看還是從矩陣的形式看，它都是一個(gè)二階張量，可以用字母F上加兩個(gè)箭頭來表示這是一個(gè)二階張量。對這個(gè)二階張量再求一次散度
最后一步把第一項(xiàng)又寫回到球坐標(biāo)的形式，第二項(xiàng)則引入了矢量
將(13)和(14)式的結(jié)果代入(11)式的矢量泊松方程，可以將基底vec{e}_?提出到拉普拉斯算符外面，從而得到只和系數(shù)ω_?有關(guān)的泊松方程
巧借氫原子球諧函數(shù)，求解球坐標(biāo)下的壓強(qiáng)場和渦度場
將壓強(qiáng)場的泊松方程(8)式和渦度系數(shù)的方程(15)寫在球坐標(biāo)下，并且注意到它們在?方向上的對稱性，可以得到
壓強(qiáng)場的方程比較簡單，先求解它。假設(shè)壓強(qiáng)場方程的解能分離成與徑向有關(guān)的部分以及與極角方向有關(guān)的部分
代入(16.a)，得到
令ψ=rf，上式可化簡為
這個(gè)方程有兩個(gè)解系，一個(gè)是ψ~r^(l+1)，一個(gè)是ψ~r^(-l)。
對于極角方向部分，取(17)式的第二部分和第三部分
令x=cos(θ),h(x)=g(θ)
這個(gè)方程的解正是勒讓德多項(xiàng)式
結(jié)合徑向和角向的結(jié)論，可以寫下壓強(qiáng)場的通解
剩下的問題是確定展開系數(shù)A和B。為此，張朝陽先從量綱的角度分析了它與哪些物理量有關(guān)，對于流體產(chǎn)生的阻力，不難想象它和在其中運(yùn)動的球狀物體的半徑R、運(yùn)動速度v?、以及流體粘滯系數(shù)μ有關(guān)
一個(gè)樸素的想法是，先將這些量的一次方相乘起來，看看量綱是什么
和壓強(qiáng)的量綱正好差L2，因此可以大膽地假設(shè)，(19)式中只包含r^(-2)的項(xiàng)。當(dāng)然這里會有一個(gè)小問題，R/r是一個(gè)無量綱量，完全可以在表達(dá)式中再乘上這一無量綱量來維持量綱的平衡。針對這個(gè)不確定性，數(shù)學(xué)家們已經(jīng)證明了，在給定足夠的邊界條件下，微分方程的解具有唯一性，因此只要求出了其中一個(gè)滿足邊界條件的解，就可以認(rèn)為得到了完整的解。
經(jīng)過這樣的量綱分析，不難發(fā)現(xiàn)只有l(wèi)=1的情況才滿足假設(shè)
再來看(16.b)的渦度場方程，徑向部分的分析和壓強(qiáng)場是一樣的，而角向部分相比(18)式會多出一項(xiàng)
多出來的第三項(xiàng)，正好是球諧函數(shù)Y在m=1時(shí)方位角分量e^(im?)被拉普拉斯算子作用后出現(xiàn)的項(xiàng)
由此可以確定，(21)式的解是連帶勒讓德多項(xiàng)式
所以渦度場的通解為
渦度是速度的散度，它具有量綱
如果和壓強(qiáng)場一樣，也樸素地假設(shè)它含有R的一次項(xiàng)和v?的一次項(xiàng)，那么可以推斷(22)式只含有r^(-2)的項(xiàng)，也就是l=1。由于
因此
降一層臺階：用斯托克斯流函數(shù)代入邊界條件
上一節(jié)利用之前的物理課解氫原子薛定諤方程的經(jīng)驗(yàn)，成功得出了壓強(qiáng)場(20)和渦度場(23)，但它們各自都有一個(gè)待定系數(shù)，需要用邊界條件來定出。然而渦度的邊界條件是難以給定的，因?yàn)楦鶕?jù)(10)式，它和速度的兩個(gè)分量vr、vθ的導(dǎo)數(shù)都有關(guān)系?；貞浺霚u度時(shí)，出發(fā)點(diǎn)是為了將(1)式右邊的二階導(dǎo)變成一階導(dǎo)，相當(dāng)于視角升了一個(gè)臺階，從而降低微分的階數(shù)。那么有沒有可能再轉(zhuǎn)化一次視角，把vr和vθ統(tǒng)一成另一個(gè)中間量，以便于引入邊界條件呢？
答案是可行的。在引入渦度時(shí)，曾將速度場類比成磁場，因?yàn)樗鼈兌季哂袩o散的性質(zhì)，可以通過(2)式將拉普拉斯算子轉(zhuǎn)化為兩個(gè)散度的依次作用。無散的矢量場還有一個(gè)性質(zhì)，就是它們總可以寫成另一個(gè)矢量場的旋度，正是這一點(diǎn)允許從磁場定義一個(gè)磁矢勢。同理，可以定義一個(gè)矢量場vec{Ψ}，它的旋度等于速度場，從而將視角降一個(gè)臺階
張朝陽類比電動力學(xué)引入斯托克斯流函數(shù)
在球坐標(biāo)下寫開各分量
和一般的在球坐標(biāo)下求散度的情形不同的是，最后一步重新定義了Ψ?，把rsin(θ)吸收了進(jìn)去，并簡記為Ψ，稱為斯托克斯流函數(shù)。再次利用軸對稱的性質(zhì)，可以知道v沒有?方向的分量，且Ψ是一個(gè)與?無關(guān)的函數(shù)，所以
也就是
到這一步，就把兩個(gè)速度分量統(tǒng)一成了一個(gè)斯托克斯流函數(shù)。將它們代入(10)式
再代入(23)式，就可以從渦度的場方程轉(zhuǎn)換到斯托克斯流函數(shù)的方程
這個(gè)偏微分方程相比(16)式要復(fù)雜得多，它很難直接地分離出徑向部分和角向部分，但可以先假設(shè)它是可分量變量的，從而猜測它應(yīng)該滿足什么樣的形式，一個(gè)最簡單的猜解是
代回檢驗(yàn)，不難發(fā)現(xiàn)這個(gè)形式能非常方便地配平(25)式中有關(guān)角向的部分的冪次
由此得到徑向部分的方程
這是一個(gè)變系數(shù)非齊次線性微分方程，可以確定它的通解是這些冪函數(shù)的和
代回原式可以定出a=D?/2。
至此，終于可以考慮速度的邊界條件了。在本問題中，有兩個(gè)邊界條件需要考慮，一個(gè)是無窮遠(yuǎn)處的流體應(yīng)當(dāng)沒有受到球狀雨滴的影響，所以保持靜止；另一個(gè)是雨滴表面的流體相對雨滴靜止，稱為no-slip condtion，這在小雷諾數(shù)情形下是普遍成立的。把(26)代入(24)，并要求無窮遠(yuǎn)處速度趨于0，可以定出b=0。再要求球面處有vec{v}(R)=vec{v}?，對應(yīng)r分量和θ分量分別是
由此可以定出
所有的待定系數(shù)都已經(jīng)得到了確定，最終結(jié)果是
在用邊界條件定下速度場和渦度場后，就可以聯(lián)立(6)和(20)來確定壓強(qiáng)場的待定系數(shù)B?了。
比較上下兩式，可以定出
所以壓強(qiáng)場為
當(dāng)然，空氣應(yīng)該還有個(gè)大氣壓強(qiáng)p?作為背景壓強(qiáng)，但加上一個(gè)常數(shù)并不影響壓強(qiáng)場的梯度，因此這里可以忽略大氣壓。
計(jì)算流體作用在球狀物體上的斯托克斯阻力
最終作用在球體上的力，就是這個(gè)二階的應(yīng)力張量與球面的法向量vec{e}?點(diǎn)乘得到的一階矢量在整個(gè)球面上的積分。正是由于二階張量有非對角項(xiàng)，點(diǎn)乘出的力并不一定垂直于法向，而是會有一定的切向分量。
(應(yīng)力張量與面法向量的點(diǎn)乘)
對于(28)式的第一項(xiàng)，將張量基底寫開并與球面法向點(diǎn)乘是
根據(jù)軸對稱性，壓強(qiáng)在球面上的積分只剩下沿極軸z方向的分量，其他方向的分量會互相抵消，所以可以只對壓強(qiáng)的z向分量做積分
同理，(28)式的第二項(xiàng)與球面法向點(diǎn)乘是
上式最后一步代入r=R時(shí)，r方向上的速度梯度剛好是一個(gè)極值點(diǎn)，所以不貢獻(xiàn)粘滯力，只有θ方向貢獻(xiàn)出來一個(gè)切向的壓強(qiáng)。將這個(gè)壓強(qiáng)的z分量在球面上積分
(28)式的第三項(xiàng)與球面法向點(diǎn)乘是
這說明這一項(xiàng)在球面不會貢獻(xiàn)壓強(qiáng)。其中第二個(gè)等號運(yùn)用了球坐標(biāo)單位基矢的求導(dǎo)關(guān)系
將以上三項(xiàng)加起來，就得到了斯托克斯定律
回顧整堂課的推導(dǎo)過程，雖然要解決的是流體力學(xué)問題，但其中貫穿了各個(gè)領(lǐng)域的知識，比如從電動力學(xué)的矢量微積分引入了渦度、速度、斯托克斯流函數(shù)的“三層樓”，借用量子力學(xué)中的球諧函數(shù)解出了渦度角向分量的泊松方程。當(dāng)然，從歷史進(jìn)程來看，數(shù)學(xué)家們先是在研究流體力學(xué)時(shí)發(fā)展了這些工具，再到后來發(fā)展電動力學(xué)和量子力學(xué)時(shí)，人們就可以直接使用這些已經(jīng)發(fā)展好的工具。而《張朝陽的物理課》研習(xí)路徑則剛好反了過來，是先學(xué)習(xí)了電動力學(xué)和量子力學(xué)，再來研究這個(gè)大部分普通物理課上一筆帶過的斯托克斯定律。張朝陽以此提示大家，學(xué)科的發(fā)展本就沒有領(lǐng)域的限制，只有將一個(gè)領(lǐng)域的內(nèi)容學(xué)扎實(shí)了，才能在面對新問題時(shí)靈活運(yùn)用先前的知識。
本堂課也是在向斯托克斯致敬。斯托克斯在劍橋執(zhí)教期間對流體力學(xué)領(lǐng)域做出很多奠基性的貢獻(xiàn)，并早在1851年就得出了今天所講的斯托克斯定律。他的研究對后世產(chǎn)生了深遠(yuǎn)的影響，比如矢量分析的工具幫助了麥克斯韋在1865年建立起電磁方程組，密立根油滴實(shí)驗(yàn)用到了斯托克斯定律來確定油滴的質(zhì)量，風(fēng)洞試驗(yàn)、汽車風(fēng)阻模擬、機(jī)翼設(shè)計(jì)等工程問題也處處有NS方程的身影。NS方程的解的存在性與光滑性至今仍是數(shù)學(xué)界的研究熱點(diǎn)，并在21世紀(jì)初被列為千禧年七大數(shù)學(xué)難題之一。盡管現(xiàn)在能借助計(jì)算機(jī)數(shù)值求解NS方程，但尋找解析解依舊是幫助人們理解物理圖像的重要手段之一。
