有限差分法实战：手把手教你用Python求解不可压NS方程（附完整代码）

原创

于 2026-03-11 02:21:53 发布 · 554 阅读

标签

有限差分法实战：用Python从零构建不可压NS方程求解器

计算流体力学（CFD）的魅力，在于它能将现实世界中那些看似混沌、难以捉摸的流动现象，转化为计算机屏幕上清晰演化的数值解。对于许多工程师和科研新手来说，从理论公式到一行行可运行的代码，这中间的路似乎总隔着一层迷雾。市面上不乏优秀的理论推导，也不缺成熟的商业软件，但恰恰是那个“自己动手，从网格开始一步步搭建求解器”的过程，最能让人深刻理解流动的每一个细节。这篇文章，就是为你准备的这样一份“手工作坊”指南。我们不依赖任何现成的CFD库，只用最基础的Python和NumPy，从最经典的二维方腔顶盖驱动流问题出发，完整地走一遍有限差分法求解不可压Navier-Stokes方程的全流程。你会看到网格如何生成，边界条件如何“说话”，压力与速度如何迭代耦合，以及那些教科书上很少提及、却能让你的程序崩溃或给出荒谬结果的“编程陷阱”。准备好了吗？让我们开始这段从方程到动画的旅程。

1. 物理问题与数学模型的建立

在动手写代码之前，我们必须清晰地定义要解决什么问题，以及描述这个问题的数学语言是什么。选择一个经典且边界条件相对简单的算例，有助于我们聚焦于数值方法本身，而不是被复杂的几何或边界条件所困扰。

1.1 方腔顶盖驱动流：一个CFD的“Hello World”

方腔顶盖驱动流，堪称计算流体力学领域的标准测试案例。想象一个正方形的空腔，里面充满流体。空腔的顶部壁面以一个恒定的速度水平运动，就像一块滑盖。其余三个壁面（左、右、下）则保持静止。顶部壁面的运动将带动腔内的流体产生复杂的旋涡运动。

注意：选择这个案例并非因为它简单，而是因为它边界条件明确（所有壁面速度已知），且存在大量公开的基准解，便于我们验证自己代码的正确性。

这个物理问题可以用不可压缩Navier-Stokes方程组来描述。对于二维流动，方程组如下：

连续性方程（质量守恒）： ∂u/∂x + ∂v/∂y = 0

动量方程（x方向）： ∂u/∂t + u ∂u/∂x + v ∂u/∂y = - (1/ρ) ∂p/∂x + ν (∂²u/∂x² + ∂²u/∂y²)

动量方程（y方向）： ∂v/∂t + u ∂v/∂x + v ∂v/∂y = - (1/ρ) ∂p/∂y + ν (∂²v/∂x² + ∂²v/∂y²)

其中：

u, v 分别是流体在 x 和 y 方向的速度分量。
p 是压力。
ρ 是流体密度（常数）。
ν 是流体的运动粘度。
t 是时间。

我们的目标就是求解出计算域内每一个点上的 u, v, p 随时间的变化。

1.2 投影法：解耦压力与速度的钥匙

直接求解上述耦合的方程组非常困难。投影法（Projection Method）是一种高效的时间推进策略，其核心思想是将一个时间步内的更新分解为两个子步：

中间速度步：忽略压力梯度项，只考虑对流和粘性项，计算出一个不满足连续性方程的“中间速度场”。
投影修正步：将上一步得到的中间速度场，投影到满足连续性方程（无散度）的子空间上，这个投影过程自然地引出了关于压力的泊松方程。求解这个压力泊松方程，得到压力场，并用压力梯度去修正速度场，使其恢复无散度。

这种方法将压力和速度解耦，使得我们在每个时间步只需要依次求解几个相对简单的方程。最经典的投影法是Chorin的投影法，我们将在代码中实现它的一种形式。

2. 计算网格与变量布置策略

数值计算是在离散的网格点上进行的，如何摆放我们的未知量（u, v, p），是有限差分法第一个关键的设计决策，它直接影响着离散方程的推导和边界条件的实施。

2.1 交错网格：为什么它优于同位网格

初学者很自然地会想：把 u, v, p 都定义在同一个网格点（单元格中心）不就好了吗？这种布局称为同位网格。但在实践中，同位网格会导致一个严重问题：压力场的空间振荡（棋盘式压力场）。这是因为动量方程和连续性方程对压力的离散形式，在同位网格上可能无法耦合，允许出现高低压交替的虚假解。

交错网格（Staggered Grid）是解决这一问题的经典方案。在交错网格中：

压力 p 定义在网格单元的中心。
x方向速度 u 定义在网格单元垂直边的中心（即，u 点的位置在 p 点的左右两侧）。
y方向速度 v 定义在网格单元水平边的中心（即，v 点的位置在 p 点的上下两侧）。

这种布置有两大优势：

自然避免棋盘压力场：压力梯度项（如 ∂p/∂x）的离散化天然地使用相邻压力点的差值，而 u 点正好位于这两个压力点的中间，离散精度高且稳定。
便于实施连续性方程：在单元中心计算速度散度 ∂u/∂x + ∂v/∂y 时，所需要的 u 和 v 正好位于该单元的四条边上，离散形式简洁。

2.2 Python中的网格实现与虚拟网格层

为了统一处理内部点和边界点，我们通常在物理计算域的外围添加一层“虚拟网格”（Ghost Cells）。这些虚拟网格点上的值不由控制方程决定，而是根据边界条件由相邻的内部点值“镜像”或“外推”得到。这能极大简化边界条件处理的代码逻辑。

假设我们的方腔物理尺寸是 Lx = Ly = 1.0，我们将其离散为 Nx 和 Ny 个内部单元。由于采用交错网格并添加虚拟层，数组的尺寸会略有不同。

让我们用代码来初始化这些场：

import numpy as np

# 参数定义
Lx = 1.0       # 方腔长度
Ly = 1.0       # 方腔高度
Nx = 41        # x方向内部单元数 (建议为奇数，便于定位中心线)
Ny = 41        # y方向内部单元数
dx = Lx / (Nx - 1)  # 网格间距 (基于内部节点数)
dy = Ly / (Ny - 1)

# 创建交错网格变量
# 压力 p 定义在单元中心，包括一层虚拟网格
# 因此，p 的数组形状为 (Nx+2

最低0.47元/天解锁文章