纯数学视角下，纳维-斯托克斯方程的全过程推导
前 言对于有流体参与的反应器而言，其内部的水力学状态是影响反应效率的重要因素之一。实验是评估水力学条件优良与否的途径之一，但这往往需要承担很多的时间和金钱成本。侧重于理论计算的流体模拟之所以能在化工、航天、车辆、船舶、能源、水利、环保和军事等领域都有很大的发挥空间，是因为其有不受客观条件限制、能模拟恶劣条件、便于调整初始和边界条件、能得出高参考价值的结果、无需承担很多时间和资金成本等优势。随着计算机技术的不断发展，个人认为流体模拟的运用空间将会被进一步开拓。 最初知道有流体模拟这一项技术还是在我读研的时候，当时我师父带着我的另一位同学做相关研究，使用的模拟软件是“Fluent”。现在开始写这一系列的文章，以记录自己业余学习流体模拟的历程，并提出对一些内容的看法和见解。在正文开始前先对这一系列文章做一些统领性的约定或说明：（1）若无特别说明，则矢量采用斜体加粗格式，标量采用斜体非加粗格式；（2）若无特别说明，则所讨论的流体均为连续流体，对于诸如混入气泡、颗粒等的非连续流体不在讨论范围内； （3）若无特别说明，则所讨论的流体均为各向同性的牛顿流体，即流体内部任意一点处的黏性正应力或黏性剪应力均为关于该点处的线性形变速率或角形变速率的线性函数，且该线性函数的形式不随坐标系选取的改变而改变； （4）所取的流体微元对于宏观无穷小但对于微观无穷大，即流体微元相对宏观尺度充分小，各宏观强度量在其内部可认为处处相同，但对于分子尺度充分大，其内部包含的微观粒子数量足够多而可以满足统计学规律； （5）若无特别说明，则所建立的三维坐标系满足右手螺旋定则。 一 、坐标系和基本的力学矢量计算如图1.1所示，建立一个三维笛卡尔直角坐标系k，并以i、j、k分别作为x方向、y方向和z方向上的基矢。 图1.1 三维笛卡尔直角坐标系（满足右手螺旋定则）设有一根线密度为ρl、长度为l的硬杆。杆的一头在坐标原点，另一头在位置（lcosθ,0, lsinθ）处并有力F作用其上，力的方向矢量为（1,0,0），如图1.2所示。 图1.2 硬杆受力示意图（xoz面视角）因为是均质硬干，所以其中点（质心）的速度vc满足动量定理： 所以杆中点处的加速度ac为： 以杆子的中点作为坐标原点建立一个新的坐标系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系原点处： 在k系点（lcosθ,0,lsinθ）处： 若θ=π/2，则： 二、 拉格朗日法和欧拉法2.1 拉格朗日法 拉格朗日法着眼于某一固定流体微元在空间中的运动过程，该流体微元在空间中的运动轨迹称为“迹线”。 流体微元的流速u可与其所经过的空间位置的坐标（x,y,z）建立对应关系： 式中：ux、uy和uz分别表示速度u在x方向、y方向和z方向上的分量大小。 如果在迹线上的任意一点（x,y,z）处取一线段微矢dl： 那么在流体微元流经该微矢量时须满足： 其中： 为线段微矢的长度； 为流体微元的运动速率。2.2 欧拉法欧拉法则着眼于坐标系中任意一个固定空间点上的物理量随时间变化而变化过程，它不像拉格朗日法那样跟踪某一个固定的流体微元的运动过程，而是从“场”的视角来考察流体，将流体覆盖的区域视作“流场”。如果存在这样一条光滑曲线，其上任意一点的切线都与该点处流体微元的速度矢量重合，那么就称该曲线为“流线”。在欧拉法框架内任意流体微元的流速u为关于其所处的时空位置的函数： 若于时刻t在流线上的任意一点（x,y,z）处取一段线段微矢dl 那么根据流线的定义有： 其中： 为线段微矢的长度； 为时刻t、点（x,y,z）处的流体微元的运动速率。 因为定常流场中空间任意一点处的物理量均不随时间t的改变而改变，包括速度u，所以这时原本体现在流线方程（2-2）中的时间变量t可略去，方程（2-2）化为如下形式： 方程（2-3）与（2-1）具有完全一样的微分形式。如果在同一个定常流场中给出同一个初始点，则显然方程（2-1）和（2-3）会积分出一样的结果，所以可知在定常流场中迹线和流线是完全重合的。 三 、流场内的随体导数、旋度和散度3.1 随体导数在一般流场内，流体微元的流速u为关于其时空坐标的函数： 过了短时间dt后，该流体微元所处的时空坐标为： 根据全微分思想可知在时刻t+dt时，该流体微元的流速为： 进而可求得该流体微元的速度变化率为： 式中： 为哈密顿算子。 不单是流速u，流体的其它物理量也满足上述形式，即： 式中：“X”表示流体的任意一种物理量。3.2 旋度和散度“旋度”用以度量矢量场中某一点上的旋转的性质。旋度是矢量，其方向表示矢量场在这一点附近旋转度最大的环量的旋转轴，它和该最大环量的旋转方向间满足右手螺旋定则，其大小则是该最大环量与旋转路径围成的面元的面积之比。 电磁理论的基础方程组——伟大的麦克斯韦方程组就含有旋度算子，应该说它对于矢量场的研究是至关重要的。 图3.1 速度场的环量计算（xoz面视角）在微小的范围内，流速关于位置的变化率视为线性。如图3.1所示，在空间中取一个同xoz面平行的微矩形ABCD，矩形的长和宽分别为dx和dz，在点A处流体在x方向和z方向上的速度分量分别为ux和uz，那么速度矢量u顺时针绕矩形周线ABCDA的环量为： 式中：微分变量前用符号“D”是为了同矩形的长“dx”和宽“dz”中的“d”相区分。 根据右手螺旋定则可知方向为j。再结合旋度的定义可得： 如果将旋度推广到三维情形，则： 这边主要是对旋度做一个基本的介绍来为后续的推导做铺垫，关于其更详细的论述就放在下一篇文章中展开了。 在矢量计算中，出镜率同旋度相仿甚至更高的另一个概念是散度，其表达式是哈密顿算子点乘矢量。对于流体的速度场u，其散度为： 从式（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-5）的第三个等号右边也表示“∠DAB变化速率的一半”。方程（4-6）的第三个等号右边也表示“流体微元所在位置处的速度场的旋度在y方向上的分量（方程（3-2））的一半”，即(rotu)y是∠DAB的角平分线AM绕点A的旋转速率的度量，其几何意义是角（β+θ）的变化速率的两倍。 不妨将方程（4-1）和（4-2）改写为： 结合方程（4-3）~（4-6），方程（4-7）和（4-8）又可改写为： 方程（4-9）和（4-10）说明了可将点A和点C之间的相对运动看作是流体微元的前表面或后表面在xoz面上的投影的线性形变、旋转和角形变的合成。将方程（4-7）和（4-8）合并为矢量形式的话： 或 式中： 表示点C相对于点A的位置矢量； 是一个表示投影转动速率的张量； 是一个表示投影角形变速率的张量。 以上结论可以推广到三维流场： 或 式中： 表示两个点之间的相对位置矢量； 是一个斜角反对称张量，表示流体微元的旋转速率； 是一个斜角对称张量，表示流体微元的角形变速率。 方程（4-13）就是著名的亥姆霍兹速度分解定理。如同之前论述的那样，亥姆霍兹速度分解定理将流体微元内两个点之间的相对运动看作是流体微元的线性形变、旋转和角形变的合成。五、 广义牛顿内摩擦定律5.1 黏性剪应力对于图4.1中的流体微元，其六个面都会受到应力作用，应力示意图如图5.1所示（图中仅画出了上、右、后三个互相垂直的面上的力）。对于流体微元的任意一个面，方向平行于面的两个力为剪应力，垂直于面的一个力为正应力。同一个面上的三个力的方向之间满足右手螺旋定则。力均等效作用于面的中心点上。应力符号“τ”的上标表示应力所在的面，下标表示力的方向。在剩余三个没有标出应力的面上，剪应力的方向同对侧面上的相反，正应力的则相同。 图5.1 流体微元所受应力的示意图（仅展示上、右、后三个互相垂直的面上的力） 先来思考最简单的一维流（流体只在一个坐标方向上有分速度）情形：基于图5.1，若流体只在x方向上有分速度，那么描述剪应力的方程在一般的本科物理化学教材中能找到：如方程（5-1）所示，这就是牛顿内摩擦定律。 式中：τx(xy)表示流体微元的上表面受到的在x方向上剪应力，具体可参看图5.1；η表示流体的动力粘度，对于各向同性的流体η不随考察方向的变化而变化；ux为流体微元在x方向上的速度分量。 但是进一步分析可以发现，剪应力不光存在于流体微元的上、下两个面上，左、右两个面上也有，而左、右两个面上的剪应力的存在是容易被忽略的。 位于流体微元上表面的剪应力τx(xy)和位于其下表面的剪应力τx(xy)-dτx(xy)方向相反。如图5.2所示，以同y轴平行的、过流体微元前后表面的中心点的轴作为转动轴，那么τx(xy)和τx(xy)-dτx(xy)产生的转动力矩为： 式中：dx、dy、dz分别表示流体微元的长、宽、高。 图5.2 转动轴示意图流体微元关于该转动轴的转动惯量为： 式中：ρ为流体密度。My(xy)产生的角加速度： 因为dx和dz都是微小量，故而角加速度αy(xy)将充分大，这显然不符合物理实际。要使这流体微元不产生充分大的转动角加速度，必须有一个同My(xy)大小几乎一样但方向相反的转动力矩作用其上，所以在该流体微元的左右两个面上也必须有剪应力。 接下来的分析将反应出流体微元上、下两个面上的x方向上的剪应力同左、右两个面上的z方向上的剪应力间存在特殊的关系。 左右两个面上的z方向上的剪应力对流体微元产生的转动力矩为： 式中：τz(yz)表示流体微元的右表面受到的在z方向上剪应力，具体可参看图5.1。对两个转动力矩产生的转动趋势进行加和： 要使流体微元不会产生充分大的角加速度，分子也必须是无穷小量，即： 式中：o表示一个无穷小量。将无穷小量均移到等号右边，得： 因为上式的等号右边都是无穷小量，故而可认为： 式（5-2）说明了在流体微元的两个互相垂直的表面上，指向共有边的两个剪应力的大小是相同的，这就是所谓的“剪应力互等定理”，它可以推广到三维流场，证明过程同上。 关于一维流的场合就分析到这里。 由于实际的流场往往是二维或三维的，方程（5-1）难以适用，所以必须去寻找度量剪应力的更本质的参数以获取描述剪应力的更通用的方程，即必须将牛顿内摩擦定律进行推广。现依旧着眼于图4.1中的流体微元的前表面或后表面在xoz面上的投影ABCD的形变过程，如图4.3所示，其中射线AM为∠DAB的角平分线： 图4.3 流体微元的前表面或后表面在xoz面上的投影的形变过程示意图若有一个观察者甲在边AD的中点、面向边BC、头顶朝向y轴负方向站立，那么他可以看到边BC被扯向了右侧，边AB往自己这边倒，其绕点A的旋转速率为： 如果有另一个观察者乙站在边AB的中点、面向边CD、头顶朝向y轴正方向站立，那么他可以看到边CD被扯向了右侧，边AD往自己这边倒，其绕点A的旋转速率同样为： 甲和乙观察到的投影ABCD的形变的速率——角形变速率——是完全一样的，流体微元的上表面和右表面上的剪应力τx(xy)和τz(yz)大小也是一样的，故而认角形变速率是度量剪应力大小的更本质的参数。在之前的一维流场景下： 所以方程（5-1）可改写为： 在本三维流场景下描述剪应力的方程维持同上的形式： 因为： 所以将方程（5-4）代入（5-3）可得： 从方程（5-5）可以看出，不光ux在z方向上的变化会对τx(xy)产生影响，uz在x方向上的变化也会对其产生影响，并且两者拥有同等的地位。方程（5-1）没有涵盖uz在x方向上的变化这一影响因素，所以它对剪应力的描述是不完备的。同理可推得流体微元的上、右和后表面上的剩余的剪应力方程，具体如下： 5.2 黏性正应力这是一个非常容易被忽略的力，我最初在自行推导纳维-斯托克斯方程的时候就忽略这一个力，致使最终导得的方程跟正确的方程有一定的差别。据说斯托克斯在初期做流体力学研究的过程中也犯了同样的“错误“。斯托克斯在推导应力方程时先做了三条假设：（1）流体是连续的。类似于胡克定律，流体微元受到的任意一个应力均是关于相关形变速率的线性函数；（2）流体是各向同性的，它的性质与方向无关，不随坐标系选取的不同而不同，在同一时刻无论坐标系如何选取，假设（1）中提及的线性函数均保持同一个形式；（3）流体静止时，流体中只剩下静压力，流体的黏性不会产生任何额外的应力； 其中：假设（1）和（2）已经在前言中给出，如果流体不连续且非各向同性，那么系统的复杂性会急剧增加，本文的大多数方程也会随之失效；而假设（3）更像是一种“自然而然的认定”。 流体的形变速率张量如（4-15）。基于假设（1），列出如下应力张量关于形变速率张量的线性方程： 式中： 为应力张量 为单位张量；a和b为不随坐标系选取的改变而改变的常量； 剪应力方程（5-5）~（5-7）同上述三条假设是自洽的，将方程（5-8）同剪应力方程（5-5）~（5-7）对比可知： 假设（2）要求流体各向同性，那么η值为不随坐标系选取的改变而改变的常量。这一点在推导剪应力方程时也已经提及。 注：个人认为直接给出方程（5-8）不太严谨，因为没有决定性的理由能让方程推导者将线性形变速率和角形变速率前的系数设置为一样，即不能排除方程会呈现出如下形式的可能性： 假设（1）要求黏性正应力是关于相关形变速率的线性函数； 由于剪应力方程（5-5）~（5-7）中不含有线性形变速率项，那么认为黏性正应力方程中不含有角形变速率项； 假设（3）要求静止流体的黏性正应力为0，所以黏性正应力方程中不存在单独的常数项。综合以上三点，以黏性正应力τx(yz)为例，可给出如下线性方程： 式中：c1和c2为比例常数。根据物理定律的对称性，既然考察的是x方向上的黏性正应力，那么y方向和z方向的地位应该是平等的，即∂uy/∂y和∂uz/∂z前的比例常数没有理由不同，所以这里为两者设了同一个值c2。 y方向和z方向的黏性正应力τx(yz)和τx(yz)的方程需要跟τx(yz)的保持对称，故而： 上述三个方程又可写成如下更便于同方程（5-8）对比的形式： 将它们同方程（5-8）对比不难看出： 即为速度的散度。以随体视角看速度散度的话，其物理意义为流体微元的体积变化速率，体积的变化速率显然不会以坐标系选取的改变而改变。以纯数学角度看的话，速度散度为流体微元形变速率张量的“迹”，为不随坐标系选取的改变而改变的不变量。故以上对b的赋值方式是满足假设（1）和（2）提出的要求的。关于张量的坐标变换问题放在下一篇文章中再展开讨论。到了这一步，光凭上述三条假设已无法确定c2的值。斯托克斯当年又追加了一个假设，即： 对于这一假设个人有如下几点的解读： （1）因为考量的是流体微元受到的正应力，所以压强p也须是随体测得的，即：在测压的瞬间测压点跟随流体微元同步运动（测压的瞬间测压点同流体微元间保持相对静止）。依照这一方式测得的压强p大小应各向一样； （2）流场中任意一点处的压强p等于该点处三个正交方向的正应力（包含压强和黏性正应力）的平均值，或者说流场中任意一点处的三个正交方向的黏性正应力的加和始终为零。（3）流场中任意一点处的三个正交方向的正应力（包含压强和黏性正应力）的平均值定义为该点处的压强p；（4）根据假设（2）提出的要求，解读（2）和解读（3）中的“三个正交方向”——即坐标系的选取——是任意的。关于黏性正应力同坐标系选取方面的详细讨论放在后续文章中展开。依托这额外的假设可得 解出： 5.3 本构方程综上，描述流体微元应力的本构方程为： 因为本构方程由牛顿内摩擦定律推广而得，所以它被称为广义牛顿内摩擦定律，而契合本构方程的流体就被称作“牛顿流体”。如果没有特别说明则所讨论的流体均为牛顿流体。 纵观推导本构方程的过程不难看出在一些问题的处理上用了外推、假设等手段，这需要推导者有很好物理直觉和理论敏感性。虽然原则上讲这些处理手段不很严谨，也难以通过实验来直接验证它们的正确性。但大量的实验结果和纳维-斯托克斯方程的整体性预测相符合，这也就间接证明了本构方程的正确性。六 、纳维-斯托克斯方程6.1 质量守恒如图4.1，将流体微元所处的长、宽、高分别为dx、dy、dz的这一微小的长方体空间作为查考对象列出质量守恒方程： 该方程等号左边表示包围该长方体空间的封闭面上的净质量通量，等号右边表示该长方体空间内质量随时间的变化率。 也可以将方程（6-1）写成更为简洁的散度的形式： 对于不可压缩流体，密度关于整个时空为常量，方程（6-2）化为： 速度的散度的物理意义是流体微元的体积变化率，对于不可压缩流体该项自然为0。6.2 动量守恒 随体考察流体微元时其质量为不变量，单位时间内其各向动量的增量如下： 式中： 表示流体微元的体积。 流体微元各向受到的作用力如下: 这三个方程的物理意义均为“合力=压力+黏性正应力+剪应力+剪应力+重力”。 先考察x方向。将表达黏性应力的本构方程（5-11）代入方程（6-7）中得： 根据冲量定理有： 所以联立方程（6-4）和（6-10）可得： 同理可得y方向和z方向上的动量守恒方程，如下： 将方程（6-11）~（6-13）进行矢量加和： 整理得： 方程（6-14）即为描述流体动量守恒的方程。由于它的导出基于本构方程，所以其仅适用于牛顿流体。 也可以直接将方程（6-7）~（6-9）写成矩阵向量的形式： 将 代入方程（6-15）可得： 又因为 联立方程（6-16）和（6-17）并整理即可得到方程（6-14）。 6.3 纳维-斯托克斯方程 综合质量守恒和动量守恒即可得纳维-斯托克斯方程： 如果流体不可压缩，则： 如果流体不可压缩且是定常流，则： 由方程组（6-20）可导出不可压缩流体的伯努利方程。 七结论纵观纳维-斯托克斯方程的推导过程可知它的适用对象为连续的牛顿流体。由于在实际的工程中大多数的流体均可以视作或近似视作连续的牛顿流体，所以该方程的适用范围是非常广泛的。虽然现在纳维-斯托克斯方程已经作为基础内容写进了专业的流体力学教材，但是在它的诞生过程中汇集了物理和数学界的各路大神，这其中有为大家所熟知的科学巨匠牛顿，有学过高数一定听说过的拉格朗日和欧拉，还有就是出现在方程名中的纳维（Claude-Louis Navier）和斯托克斯（George GabrielStokes）。科学理论往往能以简洁而优雅的方程予以描述，但发现它们的道路却总是崎岖而坎坷的。虽然纳维-斯托克斯方程极为基础而重要，但从纯数学视角看它是一个非线性的偏微分方程，受限于其求解，它位列美国克雷数学研究所 (Clay Mathematics Institute) 在2000年提出的7个千禧年大奖难题之一。由于本人目前的数学功力有限，这方面暂无法展开论述。文章转载自公众 号“北斗数学与物理”，作者：周枫来源：多相流在线
