★ FEATURED ARTICLE 第一次把纳维-斯托克斯方程从纸面搬进代码跑出第一条残差单调下降的曲线是在一个二维顶盖驱动方腔的算例上。那台机器不算快128×128 的网格跑了两千步屏幕上涡心慢慢挪到该在的位置我用鼠标量了一下中心线的速度剖面和历史文献里的基准数据能对上才敢相信前面那些乱七八糟的溢出和发散不是白熬的。这就是纳维-斯托克斯方程下文简称 N-S 方程有意思的地方它是流体力学的地基从管道输运到飞行器外流从血液流动到天气预报工程上一天都离不开它可与此同时它在三维情形下的光滑解是否永远存在、是否唯一至今还是一道没被推倒的墙被列进了公认的七个顶级数学难题。同一组方程一边被工程师天天用来做设计决策一边让数学家挠头了几十年。这篇东西写给两类人一类是想弄明白这方程到底在说什么、为什么难的读者另一类是准备动手写一个流体求解器、或者正在被仿真结果折磨的同行。我会把方程的每一项拆开讲清楚它对应的物理和数值后果把数学挑战里最核心的那几个卡点用大白话讲明白然后给出一套可以照着复现的二维方腔实操流程包括网格、时间步、压力泊松方程的处理以及我踩过的那些坑。不追求公式的严密推导追求的是你读完能自己动手跑起来并且知道结果为什么对、为什么错。1. 这对方程到底在算什么从一小团流体说起1.1 一次记账动量守恒在流动中的写法把流体想象成无数个微小控制体每个控制体都在做一件事被周围的流体推、被压力挤、被粘性拖。牛顿第二定律说的是质量乘加速度等于合外力套到一个流体微团上左边是惯性项右边是各种面力和体力写出来就是动量方程。不可压缩、常物性的形式大概长这样ρ(∂u/∂t (u·∇)u) -∇p μ∇²u f四项各自有明确身份。∂u/∂t是当地加速度也就是站在固定点上看流速怎么随时间变(u·∇)u是对流加速度描述流体微团自己带着自己跑时感受到的速度变化。等号右边-∇p是压力梯度产生的推动力方向永远从高压指向低压μ∇²u是粘性扩散它让速度场变得平滑把速度差抹平f是体积力重力、电磁力、旋转参考系里的离心力和科氏力都往这里塞。我第一次真正理解对流项是把它和站在河边看水对照起来。你站在桥上看河水速度可能一直是 2 m/s 不变当地加速度为零但一片叶子从上游漂到下游它经过的每一处速度不同叶子自己感受到的是加速或减速这就是对流项。所以同一时刻同一位置的两个观测者一个是欧拉视角固定点一个是拉格朗日视角随流体走两者差的就是这一项。数值上这一项最讨人厌因为它是非线性的是后面所有麻烦的源头。1.2 连续性方程不是附加条款而是耦合约束光有动量方程不够还得保证质量不凭空产生也不凭空消失。这由连续性方程管∂ρ/∂t ∇·(ρu) 0不可压缩流动里密度恒定它退化成一个极其简洁的约束∇·u 0。不少人初学时会把它当成一个顺便满足一下的条件觉得主力是动量方程。恰恰相反这个约束是整数值方法里最难伺候的部分因为它不含时间导数意味着压力不是独立演化的物理量而是一个拉格朗日乘子式的角色压力存在的意义是在每一时刻调整速度场让散度归零。这个认识直接决定了求解器的架构。你不可能像解温度场那样解压力场因为压力没有自己的输运方程。必须把动量方程和连续性方程当成一个耦合系统来处理要么同时求解耦合式要么把压力当作校正量分步求解分离式。后面第 4 节要讲的投影法本质上就是后者它的核心操作是把压力从动量方程里拆出来变成一个泊松方程。1.3 它为什么同时出现在工程手册和数学悬案清单里同一个方程有两种命运原因在于算和证是两件完全不同的事。工程上我们只需要在给定的边界和初值下针对特定几何、特定雷诺数得到一个足够准的数值解。只要离散格式稳定、网格够密、残差收敛结果拿来用没有问题。工程实践里的 N-S 方程是可用的因为具体问题有具体约束解的行为被边界和耗散牢牢锁住。数学上要的是更强的结论对任意满足一定光滑性的初值三维不可压缩 N-S 方程是否一定存在全局时间上的光滑解如果不是是否会在有限时间内出现奇点速度梯度爆破这个问题到现在没有答案。二维情形倒是早就被证明了全局正则所以你能看到二维算例几乎从不炸三维高雷诺数就很容易出问题。这个维度差异本身就很有意思它提示我们对流项在三维里多出来的涡拉伸机制是麻烦的关键。有一件事要提前说清楚数学上的未解决不代表工程上的数值解不可信。这两件事经常被混为一谈导致一些夸张的说法。实际做仿真时你面对的误差主要来自离散误差、湍流模型误差和边界条件不准确而不是方程的解可能不存在。2. 逐项拆解每一项背后都是一次工程取舍2.1 对流项非线性的全部麻烦都从这里长出来(u·∇)u这一项让方程变成拟线性的它带来三个直接后果。第一叠加原理失效你不能把两个解加起来得到新解第二解可能在空间中形成陡峭梯度网格跟不上就会出现非物理振荡第三它和粘性项之间存在竞争竞争的结果由雷诺数决定。数值上处理对流项有三大流派各有各的脾气。中心差分精度高但会在网格雷诺数较大时产生振荡因为它的数值粘性是负的一阶迎风无条件稳定但数值扩散严重你会看到本该尖锐的界面被抹成一片糊二阶迎风、QUICK、TVD 类格式在两者之间找平衡代价是实现复杂、边界处理麻烦。我的经验是低雷诺数、光滑解的问题中心差分配合足够密的网格完全够用一旦网格雷诺数超过 2 左右就该考虑迎风类格式了。判断网格雷诺数的公式很简单Re_cell |u|Δx/ν它衡量的是在一个网格尺度上对流和粘性谁占主导。这个数如果远大于 1说明网格根本没分辨到粘性起作用的尺度中心差分的振荡几乎是必然的。注意看到速度场出现棋盘状的交替振荡先别急着调时间步八成是对流项格式和网格尺度不匹配。2.2 无量纲化之后雷诺数把参数空间压成一维N-S 方程有一个极其漂亮的性质通过选择合适的特征长度 L、特征速度 U、特征时间 L/U、特征压力 ρU²可以把方程改写成不含任何物性参数的形式只在粘性项前面剩一个系数1/Re其中Re ρUL/μ UL/ν这件事的意义比看上去大得多。原本你要考虑密度、粘度、几何尺寸、来流速度四个量的组合现在只需要一个无量纲数。两个物理上完全不同的系统只要雷诺数相同、几何相似它们的无量纲流场就相同这就是相似准则也是风洞实验和缩比模型的理论基础。常见流动的雷诺数量级可以给个参照一滴水在毛细血管里流动Re 大约在 10⁻³ 量级粘性完全主导流动像糖浆一样规整家用自来水管里大约 10³ 到 10⁴处在过渡区汽车外流是 10⁶ 到 10⁷大型船舶和飞机是 10⁸ 以上。这个跨度意味着同一套方程在不同区间里的性格完全不同低雷诺数下它接近一个扩散方程高雷诺数下它接近一个双曲型方程。雷诺数量级典型流动主导机制数值难点10⁻³ 至 10⁰微通道、润滑膜粘性扩散几乎无难点隐式格式随便跑10² 至 10³方腔、圆柱绕流起涡惯性开始显现定常/非定常转换10⁴ 至 10⁵管道湍流、叶轮机械湍流开始充分发展网格量与模型误差10⁶ 以上飞行器、汽车外流惯性主导、边界层极薄壁面解析成本极高2.3 边界条件才是真正的项目说明书方程本身是通用的真正决定你的问题是什么的是边界条件。这句话我在带新人的时候说过很多次。动量方程是二阶的速度在每一方向需要两个边界条件压力方程是一阶的在投影法里是二阶泊松方程只需要一个。常见的边界类型就那么几种无滑移固壁速度为零、滑移壁法向速度为零、切向自由、入口给定速度剖面、出口给定压力或者零梯度、对称面、周期边界。每一种的实现细节都能坑人。比如出口给定压力很多人直接令出口边界的压力节点等于 0看起来简单实际上如果出口离尾迹太近会人为地把流动吸出去涡结构被压扁。稳妥的做法是把出口推远用零法向梯度加上整体质量守恒修正。还有一个特别容易忽略的点入口给速度剖面时如果给的是均匀流入口角落会立刻形成一个奇异点那里速度从 0 突变到 U梯度无穷大。数值上表现为入口附近的压力出现尖峰。解决办法是在入口附近用几层网格过渡或者直接给一个光滑的解析剖面抛物、双曲正切都行。提示压力泊松方程如果是全 Neumann 边界封闭腔体就是这种情况雅可比矩阵奇异直接解会失败。必须给出一个参考点压力或者在每次迭代后把压力场减去平均值。这一点后面还会细说。2.4 可压缩与不可压缩两套截然不同的数值生态∇·u 0这条约束把不可压缩流体的数学结构变得很特殊压力不再是热力学状态量而是约束的乘子声波被完全滤掉。好处是稳定、好算坏处是不能直接模拟声学、激波、密度剧变。可压缩流则完全相反密度是独立变量需要状态方程闭合压力由能量方程和状态方程给出方程是双曲型的存在声波、激波、接触间断三种特征。数值上必须用守恒型格式加通量分裂或者近似黎曼求解器还得处理激波捕捉和正定性保持。工程上选择哪一套取决于马赫数。一般 M 0.3 用不可压缩近似误差在可接受范围M 在 0.3 到 0.8 之间用低马赫预处理的可压缩求解器再往上就必须用可压缩格式了。中间地带最难受两边都不太适用。我见过一个案例用不可压缩求解器算 M0.6 的叶栅内流密度变化超过 20%结果叶背的压力分布整体偏低和实测差了一截最后只能换求解器重做。维度不可压缩可压缩压力角色约束乘子需解泊松方程热力学量由状态方程给出时间推进常需分步或耦合迭代可直接显式推进声波被滤除必须解析时间步受声速限制典型格式投影法、SIMPLE 类Godunov 类、通量分裂主要难点压力-速度耦合激波捕捉、正定性3. 数学上的硬骨头算不出来和证不出来是两回事3.1 三维光滑解从算得挺稳到证明存在之间有道墙关于三维不可压缩 N-S 方程公开的数学现状是这样的给定任意光滑、无散度的初值和足够好的外力是否存在一个在任意长时间内都保持光滑的解如果不存在那么是否存在一个有限时间点让速度梯度或等价的涡量变得无穷大这两个问题的回答至今空白。已经知道的部分也不少。弱解的存在性早在上世纪三十年代就被证明了思路是对速度场做能量意义上的估计构造一个近似序列取极限得到一个满足积分形式方程的解。这个弱解足够好物理上能用但光滑性不够唯一性证明不出来。二维情形则被完整证明全局弱解存在且光滑解的唯一性成立。所以二维算例里几乎看不到真正的数值爆破除非你的格式本身不稳定。为什么那一维之差这么大关键在于涡拉伸。在三维里涡量方程里有一项(ω·∇)u它描述涡管被流动拉伸和扭曲涡管变细涡量增强涡量增强又进一步加剧拉伸形成一个正反馈。二维里这一项要么为零要么退化成守恒形式涡量只能被对流输运和粘性耗散无法自我放大。所以三维湍流的能量级串能一路传下去而二维湍流会反向传能形成大涡。这个机制对做数值的人有直接指导意义三维高雷诺数模拟里涡拉伸带来的小尺度结构是指数级增长的任何一点数值耗散不足都会在几步之内放大成灾难。这就是为什么三维 LES 和 DNS 对格式的耗散特性极其敏感。3.2 能量不等式给到什么又拿不走什么推导正则性最核心的工具是能量估计。对不可压缩 N-S 方程把动量方程两边点乘速度在区域上积分利用∇·u 0让压力项消失分部积分后压力乘散度为零得到(1/2) d/dt ∫|u|² dx -ν ∫|∇u|² dx ∫f·u dx这个式子说得非常清楚动能的变化率等于粘性耗散永远为负起抑制作用加上外力做功。没有外力时动能单调不增。这是所有定性分析的起点。但这条不等式拿不到的恰恰是最要命的东西。它控制的是速度的 L² 范数和速度梯度的 L² 范数而在三维里这两个量还不足以控制涡量本身的增长。要控制涡量需要更强的一类空间范数目前的技术手段还不够。所有关于三维正则性的进展本质上都是在多要一点假设初值属于某个比能量空间更强的空间或者假设速度在某段时间内满足额外的可积性条件。比较有名的一类结果形如如果在某段时间内速度在某个混合时空范数下足够可积且可积指数满足某个不等式那么这段时间内解就是光滑的。翻译成人话只要流动没有在某个小区域里聚集过多能量解就不会出问题。这给了爆破一个大致的图像但既不能排除也不能确认。所以这个问题到现在还挂在墙上。3.3 尺度不变性与那些差一点点的正则性判据N-S 方程有一个漂亮的对称性如果u(x,t)是解那么对任意 λ 0u_λ(x,t) λ u(λx, λ²t)也是解。也就是说把空间压缩 λ 倍、时间压缩 λ² 倍、速度放大 λ 倍方程形式不变。这个尺度变换在分析里有大用处因为它能确定临界空间什么样的范数在这个变换下保持不变。能量∫|u|²dx在三维里不是尺度不变的它会随 λ 变化。所以能量空间不是临界空间这也解释了为什么能量估计不够用它控制的东西在不同尺度上是偏心的。真正的临界空间是一个光滑性介于 L² 和 H¹ 之间的空间在这个空间里的量在尺度变换下不变。这个观察引出了很多研究路线也带来一个让人牙痒的结论现有方法差的那一点点可能不是技巧问题而是需要一个全新的控制量。我个人的看法偏悲观短期内不太可能被解决因为过去几十年里所有努力都卡在同一个地方——三维的涡拉伸项无法被现有估计吸收。对做工程的人来说这段的意义在于两点。第一别指望完美的湍流理论很快出现湍流模型在可见的将来仍然是必需品。第二数值格式的设计要格外小心不要引入非物理的额外耗散或反耗散因为方程本身处于临界状态你对它的扰动可能被放大。4. 把方程变成代码二维顶盖驱动方腔的可复现实操4.1 方案选型原始变量法、交错网格、投影法要动手算第一件事是选方案。二维方腔这个经典算例可用路径有三条涡量-流函数法、原始变量法、格子玻尔兹曼法。涡量-流函数法很优雅二维下它自动满足连续性把两个动量方程压成一个涡量输运方程压力完全消失只剩流函数满足的泊松方程。缺点是无法直接推广到三维边界条件处理也绕。原始变量法直接解速度和压力通用性强代价是要处理压力-速度耦合。格子玻尔兹曼法介于两者之间介观尺度边界处理简单但可压缩性误差和稳定性需要留心。我选原始变量法理由是它的思路最贴近真实方程的物理并且所有踩过的坑在你转到三维或者商业软件时都用得上。具体到空间离散用有限体积 交错网格。交错网格的布局是这样的压力存在单元中心水平速度 u 存在单元的左右面上竖直速度 v 存在上下面上。为什么这么麻烦因为如果所有变量都在同一套节点上同位网格压力梯度会被中心差分算成隔点之差压力场会出现棋盘振荡而且这种振荡在连续性方程里完全看不见求解器不会自动修正它。交错网格从结构上保证压力梯度算的是相邻两个面的压差不会出现这个问题。方案连续性处理压力振荡三维推广上手难度涡量-流函数自动满足无困难中原始变量交错网格显式约束不出现直接中原始变量同位网格需 Rhie-Chow 插值需专门处理直接较难格子玻尔兹曼自动满足无直接较易4.2 网格与时间步把 CFL 和扩散数算清楚时间步不是随便给的它受两个条件约束。对流主导时看 Courant 数C max|u| · Δt / Δx ≤ 0.5扩散主导时看扩散数D ν · Δt / Δx² ≤ 0.25二维显式方腔算例里顶盖速度设为 1网格 128×128Δx 1/128 ≈ 0.0078。取 C 0.4则 Δt ≤ 0.4 × 0.0078 / 1 0.0031。再看扩散约束Re 1000 时 ν 0.001D 0.001 × 0.0031 / 0.0078² ≈ 0.051远小于 0.25说明对流约束更紧。所以这个算例是按 Courant 数来定时间步的。如果 Re 很低比如 Re 100ν 0.01D 会升到 0.51超限了此时必须按扩散约束把 Δt 压到 0.0015 左右。这就带来一个实际困境低雷诺数下粘性主导显式格式被扩散数死死卡住步数多到离谱。解决办法是用隐式格式处理粘性项或者用算子分裂把扩散项单独隐式求解。我在 Re 50 的算例里一律把扩散项做隐式时间步可以放大五到十倍。心得先用估算公式算一遍两个约束取更严的那个然后在运行的前一百步里打印实际 Courant 数。如果发现它在剧烈跳动说明速度场正在快速增长这时候要么减小时间步要么检查初始条件是不是太突兀。4.3 压力泊松方程整个流程里最费时间的环节投影法的三步走每一本教材都会写但实际写起来细节决定成败。第一步忽略压力梯度用上一步的速度场算一个中间速度u*(u* - uⁿ)/Δt -(uⁿ·∇)uⁿ ν∇²uⁿ f第二步解压力泊松方程让校正后的速度满足散度为零∇²pⁿ⁺¹ (ρ/Δt) ∇·u*第三步用压力梯度校正速度uⁿ⁺¹ u* - (Δt/ρ) ∇pⁿ⁺¹关键在于第二步。这个泊松方程在封闭方腔里是全 Neumann 边界因为壁面法向速度为零压力在壁面上的法向梯度由动量方程给出为零。全 Neumann 的泊松方程有一个相容性条件源项的积分必须等于边界通量的积分否则无解。而数值上由于舍入误差这个条件几乎不可能严格满足。两个处理办法都要用上。第一每次解完压力后把整个压力场减去平均值强制均值为零这样参考点是明确的第二在迭代求解时监控残差如果残差停在某个位置下不去说明相容性条件被破坏了要检查边界上的散度是否真的为零。解泊松方程本身也有讲究。128×128 的网格直接法比如基于快速傅里叶变换的求解器最快因为它利用了这个算例的周期性或者正弦结构几步就能给出机器精度解。但通用几何下就不能用了。常用的迭代法里SOR 实现最简单收敛因子可以通过最优松弛因子 ω 2/(1sin(π/N)) 算出来N 是单方向网格数共轭梯度法收敛更快但需要好的预条件多重网格是终极方案O(N) 复杂度但实现量大。我的建议是先用 SOR 跑通把最优松弛因子算出来跑通了再考虑共轭梯度如果网格上到 512 以上再考虑多重网格别一上来就啃硬骨头。4.4 代码骨架与三个必须校核的量下面是我一般用的骨架结构Python NumPy128×128 网格在普通笔记本上几分钟能跑几千步。为了可读性边界处理部分做了简化说明实际写的时候要单独封装函数。import numpy as np nx, ny 128, 128 Lx, Ly 1.0, 1.0 dx, dy Lx/nx, Ly/ny # 交错网格布局 # u[i,j]: i0..nx-1 (单元面), j0..ny (节点) # v[i,j]: i0..nx (节点), j0..ny-1 (单元面) # p[i,j]: i0..nx-1, j0..ny-1 (单元中心) u np.zeros((nx, ny1)) v np.zeros((nx1, ny)) p np.zeros((nx, ny)) Re 1000.0 nu 1.0 / Re dt 0.4 * dx # 由 CFL 反推先按 max|u|~1 估 U_lid 1.0 # 顶盖速度 def divergence(u, v, dx, dy): div np.zeros((nx, ny)) div[:, :] (u[:, 1:] - u[:, :-1])/dx (v[1:, :] - v[:-1, :])/dy return div def solve_pressure_poisson(src, p, dx, dy, iters200, omega1.85): # 源项 src rho/dt * div(u*), 全 Neumann均值强制为零 rhs src - src.mean() for _ in range(iters): p_new p.copy() p_new[1:-1, 1:-1] (1 - omega) * p[1:-1, 1:-1] omega * 0.25 * ( p[2:, 1:-1] p[:-2, 1:-1] p[1:-1, 2:] p[1:-1, :-2] - dx*dx*rhs[1:-1, 1:-1]) # Neumann 边界镜像等价于法向梯度为零 p_new[0, :] p_new[1, :] p_new[-1, :] p_new[-2, :] p_new[:, 0] p_new[:, 1] p_new[:, -1] p_new[:, -2] p_new - p_new.mean() p p_new return p校核环节有三个量必须每次都看不能只看残差。第一个是散度残差。算完之后全场算一遍max|∇·u|它应该和机器精度同一个量级。如果这个量在 10⁻³ 上下说明压力校正没收敛你的速度场并不是一个合法解后面再算多久都是白算。第二个是涡心位置和中心线速度剖面。方腔算例之所以经典就是因为有大家都认的参考数据可以比对。用 Re 1000 的工况把中垂线上的 u 分量沿高度画出来和公开文献里被引用最多的那组基准曲线叠在一起。如果中间段差了超过三五个百分点先别怀疑物理检查网格分辨率和收敛程度。第三个是质量守恒的累积误差。把每一步全场散度的绝对值加起来随时间画如果它缓慢上升而不是维持在噪声水平说明有系统性的边界泄漏。八成是顶盖的角点处理有问题那里两条边界交汇两个方向的速度都是强加的互相冲突。5. 踩坑实录残差曲线突然竖起来的时候5.1 现象到原因的速查表发散这件事几乎每个做流体数值的人都经历过而且第一次遇到的时候毫无头绪。我把这些年遇到的典型症状整理成一张表方便对照排查。要注意的是同一个症状可能对应多个原因要按概率从高到低依次排除。症状最可能原因排查方法前几步正常几十步后突然溢出对流项格式在网格 Re 高时失稳算网格 Re改迎风或减时间步速度场出现棋盘状正负交替同位网格压力振荡改交错网格或加 Rhie-Chow压力场整体漂移不收敛全 Neumann 缺参考点每步减去压力均值残差降不下去停在定值泊松方程迭代不足提高迭代次数或换求解器边界附近速度异常大角点处边界条件冲突检查拐角节点赋值顺序网格加密后反而更差格式精度不足或湍流模型不匹配换二阶格式重标模型参数长时间后涡结构被抹平数值扩散过大改高阶格式减小网格对称算例结果不对称求解顺序引入的方向偏置检查 Gauss-Seidel 扫描方向我印象最深的一次是一个圆柱绕流的算例Re 200理论上应该出现稳定的涡脱落。结果算了五千步尾迹是稳态的一点周期性都没有。查了两天最后发现是出口边界放在了离圆柱只有 8 倍直径的位置人为地把涡压住了。把出口推到 25 倍直径周期性立刻出来。这个坑教训很深边界条件不是形式上的东西它真的会改变物理。5.2 两个反直觉的经验网格加密不一定更好、残差下降不等于收敛第一个反直觉网格加密不一定让结果更准。这话听起来荒谬但有前提——当你用的是一阶迎风格式网格越密数值扩散越小解会越来越接近真实的、包含更多小尺度结构的解而如果这些结构没有被正确解析反而会看到解在高频上剧烈抖动甚至比粗网格更难看。真正要做的是格式和网格同步升级而不是单方面加密。二阶格式配中等网格通常比一阶格式配细网格更可靠。第二个反直觉残差下降不等于收敛。残差是方程离散后的代数残差它衡量的是当前迭代解距离满足离散方程有多远和离散方程的解距离真实物理解有多远是两回事。我遇到过残差降到 10⁻⁶、但压力场在入口附近有几十个百分点误差的情况原因是最开始几百步的瞬态被当成了定常解。定常计算的正确做法是残差降到 10⁻⁴ 之后同时监控几个关键物理量比如总压降、升力系数只有它们也稳定了才算收敛。提示做定常计算时把关键物理量随迭代步数画成曲线比只看残差曲线靠谱得多。物理量曲线的抖动幅度小于某个阈值并且趋势平坦才算真的定常。5.3 后处理阶段最容易骗自己的三件事第一件把涡量图的色标调得很大。涡量场在壁面附近有极强的峰值如果色标覆盖全量程中间区域的涡结构就全白了如果色标卡在中间壁面又全饱和。看起来漂亮的图可能是因为色标选择性地隐藏了信息。我的做法是同时出两张图一张全量程、一张截断到中间范围两张一起看。第二件用流线判断分离点。流线在分离点附近会显得很挤但网格不够密的时候分离泡可能被数值扩散抹掉看起来像没有分离。判断分离应该直接看壁面剪应力的符号变化那是定量的。壁面剪应力为零的位置就是分离点正负号变化才说明有回流。第三件忽略湍流模型的适用性。低雷诺数修正、壁面函数、可实现的 k-ε每个模型都有自己的适用区间。把壁面函数用在 y 小于 5 的网格上相当于假设第一个节点就在对数律区实际上它还在粘性底层误差能到 30% 以上。每次换模型先用一个简单验证算例标定一遍别直接上复杂几何。最后分享一个我自己坚持的习惯任何一个新算例先算一个低雷诺数的稳态版本把边界条件和网格检查一遍再逐步升到目标工况。这一步花的时间远比在目标工况下发散之后满世界找原因要少。流体这行前面的功夫省不得。
阅读完成 · 觉得有帮助？
