8.1 有限体积法与SIMPLE算法
8.1 有限体积法与SIMPLE算法
本节摘要：CFD 的主流骨架是有限体积法——把第 3 章的积分形式守恒方程逐控制体记账，通量进出的差等于内部储量变化，守恒性天然保持。真正的拦路虎是压力速度耦合：速度藏在对流项里、压力藏在梯度里，没有独立的压力方程。SIMPLE 算法用"猜压力、解速度、算修正、更新、再迭代"的循环解开这道死结。本节手写一个一维扩喷管求解器，把整条流水线走通。
第 3 章写下的纳维-斯托克斯方程，在多数几何里没有解析解——泊肃叶、库埃特那几个精确解只是幸运的特例。工程上的出路是把连续的方程变成一大堆代数方程让计算机解。有限体积法是商用与开源 CFD 软件共同的选择，原因只有一条：它对每一个小控制体执行的就是守恒律本身，进多少、出多少、存多少，账目平衡到机器精度。
一、从记账思想到离散方程
把第 3 章的控制体记账法（质量账、动量账）套到每个网格单元上：单元面上有通量进出，单元内部有源项。对一维定常问题，通用输运方程离散后形如
a_P · φ_P = a_E · φ_E + a_W · φ_W + b
其中 φ 是待求变量（速度、温度、浓度都行），下标 P 是本单元，E、W 是东、西邻居，系数 a 由对流通量与扩散通量组成。整片网格上的方程连成一个稀疏线性方程组，迭代或直接法求解。中心差分与迎风格式的差别就藏在系数里：迎风把 a_E、a_W 的权重交给上游单元，牺牲一点精度换来对流占优时的稳定性——网格粗、Pe 数大时，中心差分会给出振荡的非物理解，这是新手最常撞的第一堵墙。
import numpy as np # 一维对流扩散方程 迎风离散 系数结构演示 n, L, rho, U, gamma = 40, 1.0, 1.0, 0.5, 0.01 dx = L / n Pe = rho * U * dx / gamma # 网格佩克莱数 aE = max(-rho*U, 0) + gamma/dx # 东侧系数 对流迎风加扩散 aW = max( rho*U, 0) + gamma/dx # 西侧系数 aP = aE + aW # 定常无源 print(f"网格 Pe = {Pe:.1f} (>2 时中心差分易振荡, 迎风无条件稳定)")
网格 Pe 超过 2 时，中心差分的系数会出现负值，解随之振荡；迎风格式永远用正系数，稳定但多了一点"数值扩散"。商用软件默认的"高阶有界格式"就是在两者之间做聪明折中。
二、压力速度耦合这道死结
不可压缩流动里，压力不出现在连续性方程中，只藏在动量方程的梯度项里——你没法直接解出压力，也没法撇开压力解速度：猜一个速度场，连续性方程多半不满足；改速度满足连续性，动量方程又被破坏。两本账互相牵制。
SIMPLE（半隐式压力耦合方程组的英文缩写）的破局思路朴素得像会计对账：
- 猜一个压力场 p*；
- 用 p* 解动量方程，得速度场 u*（此时连续性不满足）；
- 由连续性方程的欠账推出压力修正方程，解出 p'；
- 用 p' 修正速度与压力：u = u* + u'，p = p* + ω·p'（欠松弛因子 ω 小于 1 防振荡）；
- 回到第 2 步，直到残差收敛。
三、手写求解器：一维扩喷管等温流动
把整套思想落到一个能跑的模型上：截面积渐变的管道，密度恒定（等温、低速近似），求沿程速度与压力。质量守恒给出各截面速度，动量方程积分给出压降。用松弛迭代模拟 SIMPLE 的外循环：
import numpy as np n = 60 x = np.linspace(0, 1, n) A = 1.0 - 0.5*x # 面积从 1 收缩到 0.5 单位任取 mdot = 0.6 # 质量流量 rho U A 恒定 rho 取 1 U = mdot / A # 连续性直接定速度 dp, K, omega = np.zeros(n), 0.1, 0.7 # 逐段损失系数与松弛因子 for it in range(300): dp_new = np.zeros(n) for i in range(1, n): dU = U[i] - U[i-1] # 动量账 压差 抵抗 加速所需的动量流变化 加摩擦损失 dp_new[i] = dp_new[i-1] + mdot*dU + 0.5*K*U[i-1]*abs(U[i-1])*A[i] dp = omega*dp_new + (1-omega)*dp # 欠松弛更新 if it % 100 == 0: print(f"迭代 {it:3d} 出口压降 {dp[-1]:.4f}") print(f"收敛后入口到出口压降 = {dp[-1]:.3f}") print(f"出口流速 {U[-1]:.2f} 为入口 {U[0]:.2f} 的 {U[-1]/U[0]:.1f} 倍 与面积比一致")
出口压降由两部分组成：加速流体做的功（动量流变化）与壁面摩擦。把 K 从 0.1 改成 0.2 重跑，你会看到压降近乎线性上抬——这就是管路设计中"沿程损失"与"加速损失"分开记账的意义：改造弯头降摩擦，救不了收缩加速的开销。
现场解读：收敛不是算完的代名词
残差曲线下降到三个数量级只是底线，还要查三项：出口流量是否还在漂（守恒性）、网格加密一档结果变多少（网格无关性）、峰值变量的位置是否物理（如收缩段压降梯度集中）。三关都过，输出才有资格进报告。
⚠️ 常见坑：欠松弛因子开太大——迭代振荡不收敛；开太小——收敛慢到怀疑人生。工程默认压力 0.3、动量 0.7 起步，收敛困难再降。
变式思考
把上面的动量账从一维换成准三维（截面平均加形状系数），就得到管网计算的常规做法；把密度换成温度的函数重跑，就摸到了可压缩一维流动的门槛——第 5 章的拉瓦尔喷管分析在数值上正是这条路线的延伸。想体验真正的 SIMPLE，可把本例改写为交错网格上的两方程耦合（速度方程加压力修正方程），结构与商用求解器内核同源。
本节要点回顾
- 有限体积法逐控制体记账，守恒性天然保持，是 CFD 主流离散框架；
- 迎风格式换稳定、中心差分换精度，网格 Pe 数是选择开关；
- 压力速度耦合是不可压缩求解的核心困难，SIMPLE 用猜测-修正-欠松弛循环解开；
- 收敛不等于可信：守恒性、网格无关性、物理性三关都要过；
- 欠松弛是稳定器：压力 0.3、动量 0.7 是工程惯用的起步组合；
- 手写一维求解器是理解商用软件黑盒的最短路径。
方法骨架立好，下一节处理里面最不确定的那一项：湍流。