为什么雨滴落下不会砸死人？《张朝阳的物理课》推导斯托克斯定律

为什么无论多高的雨滴落下来都有一个不大的末速度？普通物理课会告诉学生，因为雨滴在下落时会受到与速度正相关的空气阻力，导致雨滴的下落末速度趋于某个定值。但这个阻力与速度的关系是怎样的呢？一百七十年前，爱尔兰物理学家斯托克斯对类似的问题产生了兴趣，并提出了在小雷诺数条件下流体阻力与速度成正比的斯托克斯定律。斯托克斯当时又是如何推导的呢？

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方程，但寻找解析解依旧是帮助人们理解物理图像的重要手段之一。

据了解，《张朝阳的物理课》于每周日中午12时在搜狐视频直播，网友可以在搜狐视频APP“关注流”中搜索“张朝阳”，观看直播及往期完整视频回放；关注“张朝阳的物理课”账号，查看课程中的“知识点”短视频；此外，还可以在搜狐新闻APP的“搜狐科技”账号上，阅览每期物理课程的详细文章。