计算流体力学原理经典教材MATLAB配套源码

资 源 简 介

本项目完整复现了《Principles of Computational Fluid Dynamics》经典教材中涉及的核心数值计算方法与流体力学算法。该代码库旨在帮助学习者通过MATLAB编程深入理解计算流体力学（CFD）的数学原理与离散化技术。具体功能包括：1. 基础方程求解：提供一维及二维对流-扩散方程（Advection-Diffusion Equation）的数值解算器，用于演示标量传输过程；2. 离散化方案实现：涵盖多种空间离散格式，包括中心差分（Central Differencing）、一阶迎风（First-order Upwind）、二阶迎风及QUICK格式，并支持通过数值实验对比不同格式的数值扩散与振荡特性；3. 时间积分策略：实现了显式欧拉、隐式欧拉、Crank-Nicolson以及Runge-Kutta等时间推进方案，并包含Von Neumann稳定性分析工具；4. Navier-Stokes方程求解：针对二维不可压缩流动，实现了交错网格上的MAC（Marker-and-Cell）算法和基于同位网格的 SIMPLE 算法或投影法，能够模拟经典的顶盖驱动方腔流（Lid-driven Cavity Flow）等基准算例；5. 线性方程组求解与加速：内置了雅可比（Jacobi）、高斯-赛德尔（Gauss-Seidel）迭代法，以及几何多重网格（Geometric Multigrid）算法来加速泊松方程的求解收敛；6. 双曲型方程处理：包含浅水方程（Shallow Water Equations）或欧拉方程的求解示例，展示激波捕捉能力；7. 后处理可视化：集成专门的绘图模块，用于动态展示流场演化、误差分析及网格依赖性测试结果。

详 情 说 明

Principles of Computational Fluid Dynamics 经典教材配套源码

项目简介

本项目针对经典教材《Principles of Computational Fluid Dynamics》中的核心理论进行了代码复现。该项目旨在通过 MATLAB 编程，帮助学习者深入直观地理解计算流体力学（CFD）中的数值离散化技术、稳定性分析以及不可压缩流动的解算流程。

代码的核心是一个综合性的数值求解器，涵盖了从基础的一维对流-扩散方程特性分析，到复杂的二维不可压缩 Navier-Stokes 方程求解。

功能特性

本项目基于 main.m 文件的实际实现，包含以下核心功能：

• 一维对流-扩散方程求解与格式对比

* 构建了一维有限体积法求解器。 * 实现了 中心差分格式 (CDS)：展示其数值振荡特性。 * 实现了 一阶迎风格式 (UDS)：展示其数值扩散（人工粘性）特性。 * 包含 QUICK 格式 的逻辑演示框架。 * 支持与解析解（Exact Solution）进行实时绘图对比。

• 二维不可压缩 Navier-Stokes 求解器

* 模拟经典的 顶盖驱动方腔流 (Lid-driven Cavity Flow) 算例。 * 基于 MAC (Marker-and-Cell) 算法，采用交错网格（Staggered Grid）布局，防止压力失耦。 * 使用 投影法 (Projection Method) 进行速度-压力解耦求解。

• 高级线性方程组求解

* 在压力泊松方程的求解中，集成了 几何多重网格 (Geometric Multigrid) 算法（代码中通过 wrapper 调用）以加速收敛。 * 同时也保留了 SOR/Jacobi 迭代法的接口逻辑。

• 可视化系统

* 动态监测计算残差（散度误差）。 * 实时绘制流场、速度矢量及标量分布对比图。

系统要求

• MATLAB R2016b 或更高版本
• 无需额外的工具箱（Toolbox），代码基于 MATLAB 基础函数编写。

使用方法

1. 将所有文件下载到同一目录。
2. 在 MATLAB 中打开 main.m 文件。
3. 直接运行 main.m 即可。
4. 程序将依次执行：

* 阶段一：运行 1D 对流扩散演示，弹出窗口展示不同离散格式与解析解的对比。在控制台按任意键继续。 * 阶段二：启动 2D 顶盖驱动方腔流模拟，弹出新窗口动态展示流场演化过程。

详细功能与算法实现分析

以下内容基于 main.m 的源代码逻辑进行详细解析：

1. 全局控制逻辑

程序通过布尔变量控制运行流程：

• 1D 演示模块：由 run_1d_demo 控制，默认开启。
• 2D NS 求解模块：由 run_2d_ns 控制，默认开启。

2. 一维对流-扩散方程 (solve_1d_advection_diffusion)

该函数用于演示高 Peclet 数（Pe = 50）下的数值离散特性。

• 物理模型：求解域长度 L=1，边界条件 $phi(0)=0, phi(1)=1$。
• 离散格式实现：

* 中心差分 (CDS)：构建系数矩阵，处理除边界外的内部节点，展示在对流占优时的非物理振荡。 * 一阶迎风 (UDS)：根据流速方向选择上游节点，展示其虽然稳定但在粗网格下带来的人工扩散误差。 * QUICK 格式：代码中保留了 QUICK 格式的实现框架，展示高阶格式的基本混合逻辑。

• 结果验证：计算精确解（Exact Solution），并将所有数值解绘制在同一坐标系中进行直观对比。

3. 二维 Navier-Stokes 求解器 (solve_2d_boussinesq_cavity)

该函数是 CFD 求解器的核心，针对雷诺数 Re=100 的方腔流进行模拟。

• 网格系统：

* 采用 交错网格 (Staggered Grid)。 * 压力 $p$ 存储在网格中心。 * 水平速度 $u$ 存储在垂直网格面上（东/西面）。 * 垂直速度 $v$ 存储在水平网格面上（南/北面）。

• 求解流程 (基于投影法)：

1. 边界条件处理： * 对 $u$ 和 $v$ 施加 Dirichlet 边界条件。 * 利用 Ghost Cells（虚拟网格）处理无滑移壁面和移动顶盖（$U_{top}=1$）。 2. 动量预测步 (Predictor Step)： * 调用 solve_momentum 函数。 * 使用显式欧拉法进行时间推进。 * 离散策略：扩散项采用中心差分，对流项采用守恒型格式。 * 计算出临时速度场 ($u^*, v^*$)，此时速度场不满足连续性方程。 3. 压力泊松方程 (Pressure Poisson Equation)： * 计算中间速度场的散度：$R = frac{rho}{Delta t} nabla cdot mathbf{u}^*$。 * 调用 multigrid_vcycle_wrapper（多重网格 V 循环包装器）求解压力修正量 $p$。此步骤是算法中最耗时的部分，多重网格技术显著提升了求解效率。 4. 速度修正步 (Correction Step)： * 调用 correct_velocity 函数。 * 利用计算出的压力梯度对中间速度场进行修正，使最终速度场满足无散度条件（$nabla cdot mathbf{u} = 0$）。 5. 误差监测与可视化： * 记录每一步的最大散度误差 (div_max)。 * 使用 mesh 或 quiver 等函数定期绘制流场。

4. 关键算法函数

• solve_momentum：显式求解 u 和 v 动量方程，包含对流项（非线性项）和粘性扩散项的离散。
• compute_divergence：基于有限体积通量概念，计算每个控制体的质量源项残差。
• multigrid_vcycle_wrapper / poisson_sor：压力方程线性系统的求解器。代码结构显示优先使用多重网格方法。
• correct_velocity：标准的投影法修正步骤，将压力梯度的影响更新到速度分量上。