说实话能把不可压缩Navier-Stokes方程一步步完整推导出来的人比我预期中少很多。这不是因为它难——对理工科背景来说数学门槛也就是微积分、偏导数和一点线性代数——而是因为大部分教材喜欢把最终公式端出来把中间的物理考量藏进“经过推导可得”这句话里。我早年做管道流动计算时套公式、查经验系数都没问题直到遇到一个非标准的突缩段流动速度分布、压降和教科书给的估算对不上我才回过头老老实实把不可压缩Navier-Stokes方程自己推一遍。推完会有一种“心里踏实”的感觉方程里每一项都不是魔法都有明确的守恒律根源而且以后再看到各种CFD求解器里那些离散格式也能大概知道它到底在解什么。需要先说明的是Navier-Stokes方程不是单一条公式而是一族方程。可压缩流、不可压缩流、变密度流、多相流各自有不同形态。我们今天推的是最基础也最常见的形态单相、等温、不可压缩的牛顿流体。之所以叫“推导”意思是不直接接受结果而是从质量守恒和牛顿第二定律这两个更基本的物理事实出发一步步把流场的偏微分方程写出来。这篇文章适合正在学流体力学、准备做CFD、或者只是想把Fma变成连续介质方程的人建议你拿张纸跟着写比光看要有用得多。1. 为什么要亲手推一遍NS方程1.1 从套公式到信方程最开始用NS方程做计算的时候我完全是“工具人”心态模型选不可压、边界设好、丢给求解器跑。但那一次突缩段算出来压降和实验差挺远导师问我“你用的方程在突缩段入口还成立吗”我愣住了。后来我意识到套公式的人只能处理“已经被别人抽象好的标准场景”一旦遇到几何突变、低雷诺数、非牛顿或者温度耦合你就不知道手里的方程该不该加项、哪一项该剪掉。亲手推导一遍最大的收获不是会背方程而是建立了一种判断力每一项都能量出身的判断力。另外推导本身能帮你避开大量使用中的低级错误。比如很多人分不清“不可压缩”和“定常流动”以为流速不随时间变就是不可压缩还有人把动量方程里的压力梯度方向写反解出来流体会从低压流向高压。这些错如果从推导源头捋一遍几乎不会犯因为你知道压力项为什么带负号、粘性项为什么是速度的二阶导数、对流项为什么是那个让人头疼的非线性形式。1.2 控制体法和微元体法怎么选推导不可压缩NS方程通常有两条路。一条是控制体法把注意力放在某个固定或随流运动的有限体积上写出这个体积内质量、动量的收支平衡得到积分形式的方程再通过高斯公式变成微分形式。另一条是微元体法直接取一个无限小的流体微元用泰勒展开写出每个面的流入流出量直接得到偏微分方程。这篇文章主要走微元体法的路线。为什么选微元体法因为它的每一步都建立在“某个体积元上某种物理量的净流入等于变化量”这种直观关系上。泰勒展开虽然看起来数学味重但好处是能让“对流项从哪里冒出来”“粘性项为什么带拉普拉斯算子”“压力项为什么是负梯度”这些疑问当场解决。你如果先在控制体上理解守恒关系再用微元体法写出微分方程两条路互相印证推导就基本不会出错。我在给学生讲的时候也习惯让他们先用微元体法写一遍再用控制体法反推一次两次答案对上才算真的通了。2. 推导前的三个物理前提2.1 连续介质假定流体不是一锅分子汤流体力学里几乎所有偏微分方程的第一步都默认流体可以被看成一种连续填充空间的介质。这个假定不是理所当然的。空气在标准状态下分子平均自由程大约几十纳米如果管道直径远大于这个尺度那么把空气当作连续介质去定义“某一点的密度、速度”是合理的但如果流道缩小到微纳米量级或者气体极度稀薄单个分子的离散碰撞效应就不能忽略了这时候NS方程这套连续介质模型本身会失效。工程上常用努森数Kn λ/L来判据Kn小于0.01左右时连续介质假定足够可靠Kn再大就要考虑滑移边界甚至分子运动论模型。这个前提放在推导之前说是因为很多人把NS方程当成“水、空气的一切流动规律”其实它只覆盖连续介质范围内的牛顿流体。你推的时候可以把“流体微元”想象成一个比分子大很多、又比流场尺度小很多的“小水团”这个小水团里有大量分子所以密度、速度、压力这些宏观量都有明确意义同时它又足够小小到可以用泰勒展开只取一阶项。2.2 不可压缩条件到底压缩了什么“不可压缩”这四个字经常被误解成“密度不变”。严格来说不可压缩流体的定义是流体质点的密度沿运动轨迹不变也就是随体导数Dρ/Dt 0。对单组分均匀流体这通常就退化成密度处处为常数ρ const。但需要留个心眼海洋里盐水密度随深度变化如果忽略温度效应、只保留密度分层这类“变密度不可压缩流”也用不可压缩NS方程配合某种状态关系只是多了浮力项。所以不可压缩不等于密度处处相等而是密度的随体变化为零。另一个容易搞混的是“不可压缩”和“定常流动”。定常流动指∂/∂t 0是流场不随时间变化不可压缩指Dρ/Dt 0是流体微团密度不随时间变化两者完全不是一回事。低速气体流动经常被当作不可压缩近似理论依据是当流速远小于声速时压力变化带来的密度相对变化很小这时候把密度当常数带来的误差可以接受。这一条在推导时非常实用因为一旦接受Dρ/Dt 0质量守恒方程就会从完整的连续性方程大幅简化为一个漂亮的约束∇·u 0也就是速度场散度为零体积不膨胀也不收缩。2.3 牛顿流体本构关系是最后一块拼图在正式推导动量方程之前还需要知道流体内部的应力到底怎么描述。理想流体只有压力没有粘性所以面上受到的力处处垂直于表面现实流体多了粘性面上除了法向还有切向应力。牛顿流体的核心假设是粘性应力与速度梯度应变率成线性关系水的粘性系数在常见工况下基本恒定所以水是最典型的牛顿流体血液在细血管里、聚合物熔体、泥浆这些往往是非牛顿流体本构关系更复杂推出来的方程形态也会不同。很多初学者会问为什么动量方程里要单独讨论应力因为NS方程本质上就是连续介质版本的牛顿第二定律而力的来源除了重力这类体积力就是流体内部互相推挤和拉扯的表面力。表面力的精确描述就直接关系到最后方程里的压力项和粘性项长什么样。所以推导动量方程时最核心的环节是搞清楚应力张量怎么写这件事值得单独一章说透。3. 质量守恒不可压缩约束从哪来3.1 微元体上的流入流出平衡取一个固定在空间中的矩形微元体边长dx、dy、dz。质量守恒说的是单位时间内微元体内质量的变化等于从六个面净流入的质量。以x方向为例左侧面流入的质量流量是ρu dydz右侧面流出的质量流量用泰勒展开近似为ρu dydz ∂(ρu)/∂x · dxdydz于是x方向净流出量就是∂(ρu)/∂x dxdydz。同理y方向和z方向分别贡献∂(ρv)/∂y、∂(ρw)/∂z。把三个方向加起来得到微元体内质量的总净流出率∂(ρu)/∂x ∂(ρv)/∂y ∂(ρw)/∂z这个净流出必须等于微元体内质量的减少率-∂ρ/∂t。移项后就是完整的连续性方程∂ρ/∂t ∂(ρu)/∂x ∂(ρv)/∂y ∂(ρw)/∂z 0用矢量记号写就是∂ρ/∂t ∇·(ρu) 0。只要你每一步都盯着“流入—流出—积累”这三个词连续方程基本不会推错。这里注意我们取的微元体是固定不动的所以∂ρ/∂t是控制体内密度的局部变化率如果你取的是随流体运动的微元体得到的形式会不一样但物理本质相同。3.2 随体导数和不可压缩化简把连续性方程按乘积法则展开可以得到另一个很有启发性的形式∂ρ/∂t u·∇ρ ρ∇·u 0前面两项合起来正好是密度对时间的随体导数Dρ/Dt ∂ρ/∂t u·∇ρ。随体导数的物理意思是你坐在一个流体质点上跟着它漂看到的密度变化率。它由两部分组成一部分是某固定点密度随时间的变化局部变化另一部分是你跟随流体从一个地方走到另一个地方、因为空间密度分布不同而“感受”到的变化对流变化。这就像你站在河边看水温变化和坐着船顺流而下感受水温变化是两个不同的视角。于是连续性方程变成一句话跟随流体运动的微团密度变化率加上密度乘以体积膨胀率等于零。对不可压缩流体Dρ/Dt 0所以立刻得到∇·u 0。这个式子简洁到很多人会忽略它的出身但它其实是质量守恒在不可压缩条件下的完整表达。以后做CFD时看到求解器方程列表里第一行写着“continuity: div(U)0”你就知道它在解什么了。4. 动量方程从牛顿第二定律到应力张量4.1 作用在流体微元上的力接下来是动量方程它对应牛顿第二定律F ma在流体微元上的表达。流体微元的质量是ρdxdydz加速度是速度的随体导数Du/Dt所以左边是ρ(Du/Dt)dxdydz。右边要统计所有外力分成两类体积力和表面力。体积力最常见的是重力写成ρf dxdydzf是单位质量力表面力来自周围流体通过微元表面推挤、拉扯这个微元是推导中最需要小心的部分。以x方向为例。微元体有六个面每个面上都有一组应力它们对x方向动量的贡献不同。法向为x的两个面上正应力σ_xx通过压力或拉伸推动微元法向为y的两个面上切应力σ_yx在x方向产生剪切作用法向为z的两个面上切应力σ_zx也在x方向有贡献。把每个方向上的表面力用泰勒展开求净力(∂σ_xx/∂x ∂σ_yx/∂y ∂σ_zx/∂z) dxdydz注意这里的下标约定第一个下标表示应力的方向第二个下标表示作用面的法线方向。不同教材可能用相反约定但由于后面会用到应力张量的对称性最终结果一致。这里推荐把下标含义写在推导纸的角落否则符号搞反了会越推越乱。4.2 应力张量为什么不是三个数而是九个分量现在关键来了描述流体微元内部应力状态的为什么是一个二阶张量而不是一个简单的向量因为“一个面上的力”本身是个向量三个分量而经过微元体上任意一点可以取无数个不同朝向的面每个面上的应力向量都可能不同。柯西应力定理告诉我们如果知道三个互相垂直面比如x面、y面、z面上的应力向量任意朝向的面上的应力都能通过线性组合算出来。三个面各三个分量合起来就是九个分量记为σ_ij。这种“九个分量描述一个物理量”的对象就是二阶张量。可以类比一下你去挤一块海绵只用一个数字描述“压强”是远远不够的因为海绵的不同方向可能有的被压缩、有的被剪切你需要知道在每个截面上的正应力和切应力各是多少才能完整描述它受力的状态。流体的“内力状态”也一样九个分量并不算多这是连续介质力学的必然要求。应力张量还有一个重要性质它是对称的即σ_ij σ_ji。原因是角动量守恒。如果切应力不配对微元体会产生净扭矩而开始自转这在无外加力矩的连续介质里不允许。这个对称性在后面代入本构关系时很有用它让我们可以放心地把两个下标互换而不改变结果。4.3 速度梯度分解与牛顿本构关系有了应力张量还需要知道“应力如何和流场速度挂钩”。首先看速度梯度∂u_i/∂x_j它同样是一个二阶张量可以分解成对称部分和反对称部分∂u_i/∂x_j E_ij Ω_ij其中E_ij (∂u_i/∂x_j ∂u_j/∂x_i)/2是应变率张量描述流体微团的变形Ω_ij (∂u_i/∂x_j - ∂u_j/∂x_i)/2是旋转张量描述微团的整体转动。物理上粘性应力只在微团变形时产生纯旋转不产生剪切应力就像你搅动一杯蜂蜜如果整杯蜂蜜像刚体一样旋转没有相对变形就不会有额外阻力。这个直觉决定了本构关系只能和E_ij线性相关。牛顿流体的本构关系可以写成σ_ij -pδ_ij λδ_ij(∂u_k/∂x_k) 2μE_ij这里p是热力学压力δ_ij是克罗内克符号μ是动力粘度λ是第二粘度系数。第二粘度系数对应“体积变形引起的粘性阻力”在可压缩流中很重要但对不可压缩流体由于∇·u 0这一项整项消失。于是本构关系漂亮地简化为σ_ij -pδ_ij 2μE_ij这一步简化是很多人推导时容易漏掉的有人把可压缩的完整本构直接代入不可压缩方程结果多出根本不该有的项费半天劲又不得不划掉。先定死∇·u 0再写本构能少踩这种坑。5. 完成组装不可压缩NS方程的完整推导5.1 代入本构关系逐项展开把应力的净力项用本构关系代进去。以x方向为例∂σ_xx/∂x ∂(-p)/∂x ∂(2μE_xx)/∂x∂σ_yx/∂y ∂(-pδ_yx)/∂y ∂(2μE_yx)/∂y由于δ_yx 0压力部分消失只剩∂(2μE_yx)/∂y∂σ_zx/∂z ∂(2μE_zx)/∂z三个方向合成压力部分只贡献-∂p/∂x。粘性部分要仔细算2μ[∂E_xx/∂x ∂E_yx/∂y ∂E_zx/∂z]用E_ij的定义展开这一整坨可以化简成μ(∂²u/∂x² ∂²u/∂y² ∂²u/∂z²)。化简的关键正是利用了不可压缩条件∇·u 0消去了交叉项。这一步如果你自己动手写过就知道指标记法在这里有巨大优势能避免矢量公式展开时丢项。最终x方向动量方程可以写成ρ(∂u/∂t u∂u/∂x v∂u/∂y w∂u/∂z) -∂p/∂x μ(∂²u/∂x² ∂²u/∂y² ∂²u/∂z²) ρf_xy、z方向完全同理。合并成矢量形式就得到了大名鼎鼎的不可压缩NS方程ρ(∂u/∂t (u·∇)u) -∇p μ∇²u ρf如果两边同除以密度写成运动学形式用ν μ/ρ表示运动粘度就是CFD软件里最常见的样子∂u/∂t (u·∇)u -(1/ρ)∇p ν∇²u f到这里整个推导的“主体工程”已经完成。5.2 把方程“读”出来每一项的物理含义推导完不是结束要能把方程读成一句话。左侧第一项∂u/∂t是局部加速度描述某固定点速度随时间的变化第二项(u·∇)u是对流加速度描述流体微团在空间不均匀的速度场中移动时感受到的加速度。两者合起来是Du/Dt也就是随体加速度。右侧第一项-∇p是压力梯度力负号表示流体从高压区被推向低压区方向是压力减小的方向这个符号如果你从前面的σ_ij -pδ_ij出发会非常自然地得到绝不会写反。右侧第二项μ∇²u是粘性扩散力让速度场趋于一致好比在速度不均匀的地方“抹平”∇²是拉普拉斯算子说明粘性力依赖速度的二阶导数也就是速度梯度的梯度。右侧最后一项ρf是体积力重力是最常见的例子。对流项(u·∇)u是非线性的因为两个未知量相乘。这一项正是NS方程难解的核心线性方程可以用叠加原理非线性方程不行做数值模拟时这一项也直接决定了格式的稳定性和迭代策略。可以说理解了这项非线性你就理解了为什么NS方程直到今天仍然是数学和工程上的硬骨头。5.3 边界条件怎么接方程本身给出的是偏微分方程真正求解还需要边界条件。最常见的壁面边界是无滑移条件在固体壁面上流体相对壁面的速度为零u u_wall。这个条件是NS方程区别于无粘欧拉方程的关键没有它粘性边界层问题根本不成立。自由表面边界处理的是液面要求法向应力满足表面张力跳跃条件、切向应力连续或给定进出口边界则根据实际问题可能是给定速度分布Dirichlet或给定压力Neumann。不少初学者会问NS方程推出来之后是不是配上边界条件就能解析求解了现实是除了少数简单几何和层流工况平板、圆管、旋转流动绝大多数流动没有解析解只能数值离散。这也正是CFD存在的意义但所有数值方法本质上都在求解你刚推出来的这套方程边界条件的数学类型Dirichlet还是Neumann直接决定离散方程组的性质。所以这一步值得和推导过程放在一起理解。6. 推导后必做的两件事验证、改错6.1 用可解析的流动验证方程对不对推完方程第一件事不是急着背下来而是找一个自己能解的特殊流动把方程“验算”一遍。我最推荐库埃特流两平行平板下板静止上板以速度U沿x方向运动板间充满粘性流体。假设流动定常、沿x方向单向发展且压力沿x方向没有梯度那么连续性方程自动满足动量方程只剩0 μ d²u/dy²配上边界条件y0处u0yh处uU解得u Uy/h速度剖面是线性分布。这个结果和实验完全吻合整个验证过程大概只要十分钟却能让你确信前面推导的每一步没有原则性错误。另一个更好的验证是泊肃叶流两板静止由x方向恒定压力梯度驱动流动同样假设定常单向流。动量方程变成0 -dp/dx μ d²u/dy²解得速度剖面是抛物线最大速度出现在板间中线上流量与压力梯度成正比。这就是管道流动压降公式的理论根基。当你推完一个高度抽象的方程再亲手从它出发推出这个熟悉的抛物线速度分布那种感觉完全是“闭环了”。6.2 符号习惯和守恒型转换的坑推导中出现最多的错误基本都来自符号习惯。一是压力前的负号有的书把应力张量定义为σ_ij pδ_ij 2μE_ij正号约定p是机械压力这时动量方程右端压力梯度项就变成∇p方向不对要仔细分辨。为了统一我建议始终使用连续介质力学里σ_ij -pδ_ij 2μE_ij的形式压力是压应力所以取负号。二是应力张量下标的顺序虽然对称性让最终方程不受影响但如果混用两种约定中间推导时很容易丢项最好全程固定一种写法。还有一个实际工程中会遇到的坑动量方程有两种写法非守恒型ρ(∂u/∂t u·∇u)和守恒型∂(ρu)/∂t ∇·(ρu⊗u)。在连续方程精确满足的前提下两者数学等价但数值离散时如果离散格式不能严格保证∇·u 0两种形式会有细微差别。低马赫数可压缩求解器里这个差别有时会被放大。所以我建议做CFD的人至少要知道自己用的求解器到底解的是守恒型还是非守恒型这会影响压力震荡和稳定性。7. 无量纲化从方程到Re数再到CFD选型7.1 特征尺度操作推出来的NS方程虽然完整但信息密度太大很难一眼看出不同流动之间的相似性。无量纲化是流体力学从“解方程”走向“理解流动”的关键一步。选特征速度U、特征长度L、特征时间L/U、动压基准ρU²定义无量纲量u* u/Ux* x/Lt* t/(L/U)p* p/(ρU²)。代入动量方程两边同乘L/U²得到∂u*/∂t* (u*·∇*)u* -∇p (1/Re)∇²u (1/Fr²)f*其中Re ρUL/μ是雷诺数Fr是弗劳德数重力与惯性力之比。这是整篇推导中最有“性价比”的一步你不需要具体求解就能从Re的大小判断方程里哪一项主导。7.2 两个极限与CFD选型Re趋于无穷大时粘性项(1/Re)∇²u趋向零方程退化为无粘欧拉方程。理想流体模型适合描述大部分远离壁面的流动但壁面附近必须接边界层否则会丢掉摩擦阻力。Re趋于零时对流项与粘性项相比可忽略方程退化为线性斯托克斯方程这时流动是时间可逆的微流控、生物流体、页岩油渗流里的低雷诺数流都适用这个极限。这个判断对CFD选型非常实用。高雷诺数流动要考虑湍流模型或者大涡模拟因为直接数值模拟的网格量随Re增长极快低雷诺数流动反而可以用更简单的求解器甚至解析办法。我做项目时的习惯是先量纲分析、估计Re量级再决定用哪种数值框架这比拿到模型就套网格高效得多。说到底NS方程推导的价值不只是得到一组公式而是让这些判断都有了依据。最后再分享一条实际经验我每次带新人学这块内容都要求他们把方程里的每一项用自己的语言复述一遍不准照抄教材。有些人复述得磕磕绊绊但磕绊的过程恰恰能暴露他哪里没懂。你自己推一遍、再验证一遍、再尝试验证一个简单流动这套流程走完NS方程就真的变成你自己的工具了。