2. 流体力学控制方程与工业简化
在工业流场CFD建模中，控制方程是整个求解器的物理内核。无论使用Fluent、CFX、OpenFOAM还是STAR-CCM+，软件在后台求解的都是同一组偏微分方程：连续性方程、动量方程和能量方程。理解这三类方程的物理含义、耦合关系以及工程简化条件，是判断仿真结果是否可信的前提。本节从方程本身出发，逐步过渡到湍流封闭问题和工业适用条件，帮助读者建立“方程—离散—求解—验证”的完整认知链条。
一、三大守恒方程的物理含义
连续性方程描述质量守恒。对不可压缩流体，其微分形式为 ∂u/∂x + ∂v/∂y + ∂w/∂z = 0，即速度场的散度为零。物理意义是：流体微元在运动过程中体积不变，流入某控制体的质量等于流出的质量。在工业管道、泵阀内流、换热器壳程等场景中，连续性方程直接约束了流量分配和速度场形态。若仿真中出现局部速度异常增大而流量不守恒，通常意味着连续性方程离散格式或边界条件设置存在问题。
动量方程描述牛顿第二定律在流体中的表达。以x方向为例，其形式为 ρ(∂u/∂t + u∂u/∂x + v∂u/∂y + w∂u/∂z) = -∂p/∂x + μ(∂²u/∂x² + ∂²u/∂y² + ∂²u/∂z²) + Fx。方程左侧是惯性项，右侧依次为压力梯度、黏性扩散和体积力。动量方程决定了流场的压力分布、涡旋结构和壁面剪切应力，是CFD求解中最核心也最耗时的部分。在工业应用中，压力-速度耦合算法（如SIMPLE、PISO、Coupled）本质上就是在联立求解连续性与动量方程。
能量方程描述热力学第一定律。其通用形式为 ρcp(∂T/∂t + u∂T/∂x + v∂T/∂y + w∂T/∂z) = ∂/∂x(k∂T/∂x) + ∂/∂y(k∂T/∂y) + ∂/∂z(k∂T/∂z) + S。左侧为热惯性项，右侧为导热项和源项。当流动存在显著温差、高速可压缩效应或化学反应时，能量方程不可忽略。工业换热器、燃烧室、涡轮叶片冷却等场景必须求解能量方程；而常温水管流动、低速风道等场景可仅求解连续性与动量方程，以降低计算成本。
二、湍流封闭问题的本质
直接数值模拟（DNS）可以求解所有尺度的湍流涡，但网格量随雷诺数Re的9/4次方增长。工业问题Re通常在10⁵～10⁷量级，DNS计算成本不可接受。因此工程CFD普遍采用雷诺平均（RANS）方法，将瞬时速度分解为时均量和脉动量：u = ū + u'。代入动量方程后，会出现额外的雷诺应力项 -ρu'ᵢu'ⱼ，导致方程组未知量多于方程数，这就是湍流封闭问题。
解决封闭问题的途径是引入湍流模型，用已知量近似表达雷诺应力。常用模型分为涡黏模型和雷诺应力模型两大类。涡黏模型基于Boussinesq假设，将雷诺应力表示为湍流黏度μt与平均应变率的乘积。k-ε模型适用于高雷诺数充分发展湍流，计算稳健、收敛快，是工业外流和管道内流的默认选择；k-ω SST模型在近壁面采用k-ω形式、远场采用k-ε形式，兼顾壁面边界层精度和自由剪切层稳定性，适合换热器、涡轮机械等有逆压梯度的场景；Spalart-Allmaras模型为单方程模型，计算量小，常用于航空外流。LES则直接解析大尺度涡，仅模化小尺度涡，精度更高但网格和计算资源需求显著增加，适合分离流、燃烧等复杂非稳态问题。
三、工程简化条件与判断依据
工业CFD建模并非总是求解完整方程。合理简化可以大幅降低计算成本，但必须建立在明确的物理判断之上。以下给出常见简化条件及其判断依据。
| 简化假设 | 适用条件 | 判断依据 | 典型工业场景 |
|---|---|---|---|
| 不可压缩 | 马赫数 Ma < 0.3 | 密度变化小于5% | 水管、风道、泵阀内流 |
| 稳态求解 | 流动特征时间远小于工况变化时间 | 边界条件不随时间变化 | 连续换热器、稳态管道输送 |
| 忽略能量方程 | 温差小、低速、无化学反应 | 布辛涅斯克近似成立 | 常温通风、水力学计算 |
| 二维/轴对称 | 几何和流动沿某一方向均匀 | 周向梯度可忽略 | 长直管道、轴对称喷嘴 |
| 层流假设 | Re < 2300（圆管） | 黏性力主导 | 微通道、高黏度聚合物流动 |
四、从方程到求解：一个工业管道示例
以DN200工业水管为例，水温20°C，流速2 m/s，管长10 m。首先计算雷诺数：Re = ρuD/μ = 998×2×0.2/0.001 ≈ 3.99×10⁵，远大于2300，属于湍流。马赫数Ma = u/c ≈ 2/1480 ≈ 0.0014，远小于0.3，可作不可压缩处理。温差若小于5°C，能量方程可忽略。因此该问题可简化为不可压缩稳态湍流问题，选用k-ε模型配合标准壁面函数。
网格划分时需关注y+值。标准壁面函数要求第一层网格节点落在对数律层，通常30 < y+ < 300。若y+过小，壁面函数失效；若过大，则无法捕捉边界层速度梯度。可通过经验公式估算第一层网格高度：Δy = y+·μ/(ρ·uτ)，其中uτ为摩擦速度，可由达西摩擦因子估算。对于本例，若目标y+≈50，估算Δy约为0.3 mm。求解设置中，压力-速度耦合选用SIMPLE算法，动量方程采用二阶迎风格式，湍动能和耗散率采用一阶迎风以保证收敛，收敛后再切换二阶。残差监控目标设为1×10⁻⁴，同时监控出口质量流量和壁面剪切应力是否稳定。
五、方程选择与验证的闭环
控制方程的选择不是一次性决策，而应贯穿建模始终。建议在仿真前填写一份“方程与简化决策表”，记录每个简化假设的依据和预期影响；仿真后通过网格无关性验证、时间步无关性验证（瞬态问题）和实验/理论对比来确认简化是否合理。若发现偏差，应回溯到方程层面重新审视假设，而非仅调整求解器参数。只有方程、模型、网格和边界条件四者自洽，工业CFD结果才具备工程参考价值。