纯数学视角下，纳维-斯托克斯方程的全过程推导

一 、坐标系和基本的力学矢量计算

 

 

图1.2 硬杆受力示意图（xoz面视角）

以杆子的中点作为坐标原点建立一个新的坐标系k’。k‘系的坐标轴同k系的平行且同向，如图1.3所示。

图1.3 k系和k‘系示意图（xoz面视角）

显然k’系相对k系存在一个朝向x轴正方向的加速度ac。在k’系中，杆中心点是静止的。在k’系中看来杆受到一个朝向x‘轴负方向的惯性力，但惯性力的等效作用点是杆的中点，所以它不会使杆产生转动的趋势。只有F会让杆旋转，其旋转角加速度α’为：

以上方程的物理意义是：角加速度等于位矢叉乘力除以转动惯量，这里需要注意旋转参考基点的选取。

所以在k’系中看来，杆上各点的加速度为：

式中：x‘和z‘表示杆上的点在k’系中的空间位置坐标。

式中：x‘和z‘表示杆上的点在k’系中的空间位置坐标。

以上方程的物理意义是：角加速度叉乘位矢等于加速度，这里同样需要注意旋转参考基点的选取。

现在将方程（1-3）切换到k系中，那么表示杆上各点的加速度的方程为：

式中：x和z表示杆上的点在k系中的空间位置坐标。

这里需要处理坐标系之间的变换问题，包括空间位置的变换和加速度的叠加。在相对运动速度远低于光速的情况下使用伽利略变换即可，而不需要用到相对论变换。

在k系中，将杆上的点的空间位置坐标值代入方程（1-4）中，即可求得各点的加速度：

在k系原点处：

若θ=π/2，则：

流体微元的流速u可与其所经过的空间位置的坐标（x,y,z）建立对应关系：

为流体微元的运动速率。

2.2 欧拉法

欧拉法则着眼于坐标系中任意一个固定空间点上的物理量随时间变化而变化过程，它不像拉格朗日法那样跟踪某一个固定的流体微元的运动过程，而是从“场”的视角来考察流体，将流体覆盖的区域视作“流场”。如果存在这样一条光滑曲线，其上任意一点的切线都与该点处流体微元的速度矢量重合，那么就称该曲线为“流线”。

在欧拉法框架内任意流体微元的流速u为关于其所处的时空位置的函数：

那么根据流线的定义有：

3.1 随体导数

在一般流场内，流体微元的流速u为关于其时空坐标的函数：

不单是流速u，流体的其它物理量也满足上述形式，即：

式中：“X”表示流体的任意一种物理量。

3.2 旋度和散度

电磁理论的基础方程组——伟大的麦克斯韦方程组就含有旋度算子，应该说它对于矢量场的研究是至关重要的。

图3.1 速度场的环量计算（xoz面视角）

在微小的范围内，流速关于位置的变化率视为线性。如图3.1所示，在空间中取一个同xoz面平行的微矩形ABCD，矩形的长和宽分别为dx和dz，在点A处流体在x方向和z方向上的速度分量分别为ux和uz，那么速度矢量u顺时针绕矩形周线ABCDA的环量为：

如果将旋度推广到三维情形，则：

从式（3-4）可看出散度为标量，它可用于表征矢量场在空间各点处的发散或吸收的强弱程度。由于散度相较旋度更容易理解，所以这里就不做更多的介绍或说明了。

在时刻t在一个三维流场中取一个长方体流体微元作为考察对象，该流体微元的六个表面均同坐标轴垂直或平行，其长、宽、高分别为dx、dy、dz，如图4.1所示。

图4.1 三维流场中的流体微元（考察对象）

现考察该流体微元的前表面或后表面在xoz面上的投影，将该投影记为ABCD，显然在时刻t它是长方形的，如图4.2所示。

图4.2 流体微元的前表面或后表面在xoz面上的投影示意图

需要注意的是，此处考察的是流体微元的前表面或后表面在xoz面上的投影，而非其前表面或后表面本身，因为在三维流场中随着流体微元的运动其前、后表面并不始终跟xoz面平行。

在时刻t，点A在x方向和z方向的速度分量分别为ux和uz，那么点C在x方向和z方向的速度分量ux‘和uz’分别为：

虽然点A和点C并非流体微元上的点，但根据投影法则它们在x方向和z方向上的运动速度同流体微元上对应的真实的点的是一样，所以ux、uz、 ux‘和uz’其实表示的就是流体微元上对应的真实的点的运动速度。

因为流体微元在运动过程中受到黏性剪应力的作用而会发生形变，所以投影ABCD也会跟着发生形变。在时刻t，投影ABCD是一个长方形，在经过了短时间dt后，它会变成菱形，如图4.3所示。

图4.3 流体微元的前表面或后表面在xoz面上的投影的形变过程示意图

将边AD的长度记为LAD、边AB的长度记为LAB的，在时刻t，边AD和边AB的线性形变速率（单位长度的长度变化速率）分别为：

两式中：dlAD/dt和dlAB/dt分别表示边AD和边AB的线性形变速率。

射线AM为∠DAB的角平分线，在时刻t，角θ的变化速率和射线AM绕点A的旋转速率分别为：

不妨将方程（4-1）和（4-2）改写为：

方程（4-9）和（4-10）说明了可将点A和点C之间的相对运动看作是流体微元的前表面或后表面在xoz面上的投影的线性形变、旋转和角形变的合成。

将方程（4-7）和（4-8）合并为矢量形式的话：

是一个表示投影转动速率的张量；

或

是一个斜角反对称张量，表示流体微元的旋转速率；

方程（4-13）就是著名的亥姆霍兹速度分解定理。如同之前论述的那样，亥姆霍兹速度分解定理将流体微元内两个点之间的相对运动看作是流体微元的线性形变、旋转和角形变的合成。

5.1 黏性剪应力

对于图4.1中的流体微元，其六个面都会受到应力作用，应力示意图如图5.1所示（图中仅画出了上、右、后三个互相垂直的面上的力）。对于流体微元的任意一个面，方向平行于面的两个力为剪应力，垂直于面的一个力为正应力。同一个面上的三个力的方向之间满足右手螺旋定则。力均等效作用于面的中心点上。应力符号“τ”的上标表示应力所在的面，下标表示力的方向。在剩余三个没有标出应力的面上，剪应力的方向同对侧面上的相反，正应力的则相同。

先来思考最简单的一维流（流体只在一个坐标方向上有分速度）情形：

基于图5.1，若流体只在x方向上有分速度，那么描述剪应力的方程在一般的本科物理化学教材中能找到：如方程（5-1）所示，这就是牛顿内摩擦定律。

位于流体微元上表面的剪应力τx(xy)和位于其下表面的剪应力τx(xy)-dτx(xy)方向相反。如图5.2所示，以同y轴平行的、过流体微元前后表面的中心点的轴作为转动轴，那么τx(xy)和τx(xy)-dτx(xy)产生的转动力矩为：

图5.2 转动轴示意图

流体微元关于该转动轴的转动惯量为：

式中：ρ为流体密度。

My(xy)产生的角加速度：

式中：τz(yz)表示流体微元的右表面受到的在z方向上剪应力，具体可参看图5.1。

对两个转动力矩产生的转动趋势进行加和：

要使流体微元不会产生充分大的角加速度，分子也必须是无穷小量，即：

式中：o表示一个无穷小量。

将无穷小量均移到等号右边，得：

由于实际的流场往往是二维或三维的，方程（5-1）难以适用，所以必须去寻找度量剪应力的更本质的参数以获取描述剪应力的更通用的方程，即必须将牛顿内摩擦定律进行推广。

现依旧着眼于图4.1中的流体微元的前表面或后表面在xoz面上的投影ABCD的形变过程，如图4.3所示，其中射线AM为∠DAB的角平分线：

图4.3 流体微元的前表面或后表面在xoz面上的投影的形变过程示意图

若有一个观察者甲在边AD的中点、面向边BC、头顶朝向y轴负方向站立，那么他可以看到边BC被扯向了右侧，边AB往自己这边倒，其绕点A的旋转速率为：

所以方程（5-1）可改写为：

从方程（5-5）可以看出，不光ux在z方向上的变化会对τx(xy)产生影响，uz在x方向上的变化也会对其产生影响，并且两者拥有同等的地位。方程（5-1）没有涵盖uz在x方向上的变化这一影响因素，所以它对剪应力的描述是不完备的。

同理可推得流体微元的上、右和后表面上的剩余的剪应力方程，具体如下：

5.2 黏性正应力

这是一个非常容易被忽略的力，我最初在自行推导纳维-斯托克斯方程的时候就忽略这一个力，致使最终导得的方程跟正确的方程有一定的差别。据说斯托克斯在初期做流体力学研究的过程中也犯了同样的“错误“。

斯托克斯在推导应力方程时先做了三条假设：

（1）流体是连续的。类似于胡克定律，流体微元受到的任意一个应力均是关于相关形变速率的线性函数；

（2）流体是各向同性的，它的性质与方向无关，不随坐标系选取的不同而不同，在同一时刻无论坐标系如何选取，假设（1）中提及的线性函数均保持同一个形式；

流体的形变速率张量如（4-15）。基于假设（1），列出如下应力张量关于形变速率张量的线性方程：

剪应力方程（5-5）~（5-7）同上述三条假设是自洽的，将方程（5-8）同剪应力方程（5-5）~（5-7）对比可知：

假设（3）要求静止流体的黏性正应力为0，所以黏性正应力方程中不存在单独的常数项。

综合以上三点，以黏性正应力τx(yz)为例，可给出如下线性方程：

上述三个方程又可写成如下更便于同方程（5-8）对比的形式：

即为速度的散度。以随体视角看速度散度的话，其物理意义为流体微元的体积变化速率，体积的变化速率显然不会以坐标系选取的改变而改变。以纯数学角度看的话，速度散度为流体微元形变速率张量的“迹”，为不随坐标系选取的改变而改变的不变量。故以上对b的赋值方式是满足假设（1）和（2）提出的要求的。关于张量的坐标变换问题放在下一篇文章中再展开讨论。

到了这一步，光凭上述三条假设已无法确定c2的值。斯托克斯当年又追加了一个假设，即：

（2）流场中任意一点处的压强p等于该点处三个正交方向的正应力（包含压强和黏性正应力）的平均值，或者说流场中任意一点处的三个正交方向的黏性正应力的加和始终为零。

（3）流场中任意一点处的三个正交方向的正应力（包含压强和黏性正应力）的平均值定义为该点处的压强p；

（4）根据假设（2）提出的要求，解读（2）和解读（3）中的“三个正交方向”——即坐标系的选取——是任意的。

关于黏性正应力同坐标系选取方面的详细讨论放在后续文章中展开。

依托这额外的假设可得

5.3 本构方程

综上，描述流体微元应力的本构方程为：

纵观推导本构方程的过程不难看出在一些问题的处理上用了外推、假设等手段，这需要推导者有很好物理直觉和理论敏感性。虽然原则上讲这些处理手段不很严谨，也难以通过实验来直接验证它们的正确性。但大量的实验结果和纳维-斯托克斯方程的整体性预测相符合，这也就间接证明了本构方程的正确性。

6.1 质量守恒

如图4.1，将流体微元所处的长、宽、高分别为dx、dy、dz的这一微小的长方体空间作为查考对象列出质量守恒方程：

速度的散度的物理意义是流体微元的体积变化率，对于不可压缩流体该项自然为0。

随体考察流体微元时其质量为不变量，单位时间内其各向动量的增量如下：

根据冲量定理有：

整理得：

代入方程（6-15）可得：

如果流体不可压缩，则：

文章转载自公众 号“北斗数学与物理”，作者：周枫