1. 背景介绍
1.1 方腔驱动流问题
方腔驱动流（Lid-driven Cavity Flow）是计算流体力学（CFD）中的经典基准问题。该问题描述了一个二维正方形腔体，其顶部壁面以恒定速度运动，而其他三个壁面保持静止。这是一个理想化的流动问题，具有以下特点：
-
几何简单：正方形计算域（1×1单位正方形），边界条件明确 -
物理丰富：包含边界层、涡旋、分离、再附着等复杂流动现象 -
验证标准：有大量文献数据可供验证，是CFD算法验证的黄金标准 -
数值挑战：高雷诺数下存在数值稳定性问题，需要精细的数值处理
1.1.1 几何配置
U = 1 (top lid velocity)
↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓
┌─────────────────────────────┐
│ │
│ │
│ ● │ ← y = Ly = 1
│ │
│ │
│ │
└─────────────────────────────┘
0 x = Lx = 1
1.1.2 边界条件
# 顶壁面 (y = Ly): u = 1, v = 0 (驱动边界)
# 底壁面 (y = 0): u = 0, v = 0 (无滑移壁面)
# 左壁面 (x = 0): u = 0, v = 0 (无滑移壁面)
# 右壁面 (x = Lx): u = 0, v = 0 (无滑移壁面)
1.2 控制方程
不可压缩流体的运动由Navier-Stokes方程组描述：
连续性方程（质量守恒）：
∂u/∂x + ∂v/∂y = 0
动量方程：
∂u/∂t + u(∂u/∂x) + v(∂u/∂y) = -∂p/∂x + (1/Re)(∂²u/∂x² + ∂²u/∂y²)
∂v/∂t + u(∂v/∂x) + v(∂v/∂y) = -∂p/∂y + (1/Re)(∂²v/∂x² + ∂²v/∂y²)
其中：
-
u, v：x和y方向的速度分量 [LT⁻¹] -
p：压力 [ML⁻¹T⁻²] -
Re：雷诺数，Re = ρUL/μ [无量纲] -
ρ：密度 [ML⁻³]，U：特征速度 [LT⁻¹]，L：特征长度 [L]，μ：动力粘度 [ML⁻¹T⁻¹]
1.2.1 无量纲化处理
代码中采用以下无量纲化方案：
-
长度尺度：L = 1（方腔边长） -
速度尺度：U = 1（顶盖速度） -
时间尺度：T = L/U = 1 -
压力尺度：P = ρU² = 1
1.3 流动特征与雷诺数效应
1.3.1 雷诺数定义
Re = ρUL/μ = UL/ν
其中 ν = μ/ρ 是运动粘度。
1.3.2 不同雷诺数下的流动特征
**低雷诺数 (Re < 100)**：
-
流动为层流，流场结构简单 -
存在一个主要的涡旋结构 -
数值求解相对容易
**中等雷诺数 (100 < Re < 1000)**：
-
流场变得复杂，出现多个涡旋结构 -
在底部角落出现小的分离涡 -
本代码模拟的 Re = 1000 属于此范围
**高雷诺数 (Re > 1000)**：
-
流动可能变为非定常或湍流 -
边界层变薄，对网格分辨率要求更高 -
数值稳定性成为主要挑战
2. 应用领域
2.1 工程应用
-
微流体器件：芯片实验室、微混合器设计 -
血液检测芯片中的流场控制 -
药物混合微设备的优化设计 -
搅拌设备：反应器、混合器的流场分析 -
化工反应器内的混合效率评估 -
食品加工中的搅拌器设计 -
润滑系统：轴承、密封件中的流动分析 -
轴承间隙内的润滑油流动 -
机械密封面的流体力学分析 -
涂层工艺：旋涂、刮涂过程中的流动控制 -
半导体制造中的光刻胶旋涂 -
汽车喷漆工艺的流场优化
2.2 科学研究
-
湍流研究：层流到湍流的转捩机理 -
临界雷诺数的确定 -
不稳定性发展机制研究 -
数值方法验证：CFD算法的准确性验证 -
新开发数值格式的benchmark test -
网格收敛性分析标准 -
稳定性分析：流动失稳机理研究 -
线性稳定性分析 -
非线性动力学行为研究 -
传热传质：耦合传输现象研究 -
自然对流与强制对流的耦合 -
传质过程中的流动影响
2.3 教学价值
-
CFD入门：理想的教学案例 -
编程实践：数值算法实现训练 -
物理理解：流体力学现象的直观展示 -
验证方法：科学计算结果验证的重要性
3. 数据说明
3.1 参考数据来源
代码中使用了两组重要的基准数据，用于验证数值结果的准确性：
3.1.1 Ghia等（1982）数据
# 沿中心线的位置坐标
x_Ghia = [0.0000, 0.0625, 0.0703, 0.0781, 0.0938, 0.1563, 0.2266,
0.2344, 0.5000, 0.8047, 0.8594, 0.9063, 0.9453, 0.9531,
0.9609, 0.9688, 1.0000]
y_Ghia = [0.0000, 0.0547, 0.0625, 0.0703, 0.1016, 0.1719, 0.2813,
0.4531, 0.5000, 0.6172, 0.7344, 0.8516, 0.9531, 0.9609,
0.9688, 0.9766, 1.0000]
# 对应的速度值
u_Ghia = [0.00000, -0.18109, -0.20196, -0.22220, -0.29730, -0.38289,
-0.27805, -0.10648, -0.06080, 0.05702, 0.18719, 0.33304,
0.46604, 0.51117, 0.57492, 0.65928, 1.00000]
v_Ghia = [0.00000, 0.27485, 0.29012, 0.30353, 0.32627, 0.37095,
0.33075, 0.32235, 0.02526, -0.31966, -0.42665, -0.51500,
-0.39188, -0.33714, -0.27669, -0.21388, 0.00000]
特点：
-
经典的高精度数值解 -
广泛用于CFD验证 -
基于涡度-流函数方法 -
数据点选择在关键流动区域
3.1.2 Erturk等（2005）数据
# 更高分辨率的坐标点
x_Erturk = [0.0000, 0.015, 0.030, 0.045, 0.060, 0.075, 0.090, 0.105,
0.120, 0.135, 0.150, 0.500, 0.850, 0.865, 0.880, 0.895,
0.910, 0.925, 0.940, 0.955, 0.970, 0.985, 1.000]
y_Erturk = [0.0000, 0.020, 0.040, 0.060, 0.080, 0.100, 0.120, 0.140,
0.160, 0.180, 0.200, 0.500, 0.900, 0.910, 0.920, 0.930,
0.940, 0.950, 0.960, 0.970, 0.980, 0.990, 1.000]
# 对应的速度分布
u_Erturk = [0.0000, -0.0757, -0.1392, -0.1951, -0.2472, -0.2960,
-0.3381, -0.3690, -0.3854, -0.3869, -0.3756, -0.0620,
0.3838, 0.3913, 0.3993, 0.4101, 0.4276, 0.4582,
0.5102, 0.5917, 0.7065, 0.8486, 1.0000]
v_Erturk = [0.0000, 0.1019, 0.1792, 0.2349, 0.2746, 0.3041,
0.3273, 0.3460, 0.3605, 0.3705, 0.3756, 0.0258,
-0.4028, -0.4407, -0.4803, -0.5132, -0.5263, -0.5052,
-0.4417, -0.3400, -0.2173, -0.0973, 0.0000]
特点：
-
更高分辨率的计算结果 -
提供了更详细的流场信息 -
使用现代数值方法计算 -
在边界层区域有更密集的数据点
3.2 数据可视化
代码在开始时就展示了这两组基准数据的对比：
plt.figure(figsize=(8, 4))
# u速度沿垂直中心线的分布
plt.subplot(1, 2, 1)
plt.scatter(u_Ghia, y_Ghia, marker="1", alpha=.7, label="Ghia et al. 1982")
plt.scatter(u_Erturk, y_Erturk, marker="2", alpha=.7, label="Erturk et al. 2005")
plt.legend(loc="lower right")
plt.xlabel("$u$")
plt.ylabel("$y$")
plt.grid(alpha=.3)
# v速度沿水平中心线的分布
plt.subplot(1, 2, 2)
plt.scatter(x_Ghia, v_Ghia, marker="1", alpha=.7, label="Ghia et al. 1982")
plt.scatter(x_Erturk, v_Erturk, marker="2", alpha=.7, label="Erturk et al. 2005")
plt.legend(loc="upper right")
plt.xlabel("$x$")
plt.ylabel("$v$")
plt.grid(alpha=.3)
3.3 数据验证意义
-
算法验证：确保数值方法的正确性 -
精度量化：评估数值解的准确程度 -
网格收敛：判断网格是否足够精细 -
物理一致性：验证流动现象的物理合理性
4. 算法和数学公式原理
4.1 投影方法（Projection Method）
本代码采用Chorin的投影方法（也称为分数步方法），将速度和压力的耦合求解分为两个步骤。这种方法的核心思想是：
-
预测步：暂时忽略压力梯度，预测中间速度场 -
投影步：通过求解压力泊松方程，将中间速度场投影到无散度空间
4.1.1 算法流程
步骤1：预测步（不考虑压力梯度）
û = u^n + Δt[-∇·(uu) + (1/Re)∇²u]
v̂ = v^n + Δt[-∇·(uv) + (1/Re)∇²v]
在代码中的实现：
# 时间积分（显式欧拉）
u_hat[2:-2, 2:-2] = u_old[2:-2, 2:-2] + dt * (- conv_u + diff_u)
v_hat[2:-2, 2:-2] = v_old[2:-2, 2:-2] + dt * (- conv_v + diff_v)
步骤2：投影步（压力修正）
∇²p = (∇·û)/Δt
u^(n+1) = û - Δt∇p
v^(n+1) = v̂ - Δt∇p
在代码中的实现：
# 计算中间速度的散度
div_u_hat = u_hat_x + v_hat_y
b[2:-2, 2:-2] = div_u_hat / dt # 压力泊松方程右端项
# 求解压力泊松方程（迭代求解）
p[2:-2, 2:-2] = 1. / (2. * (dx**2 + dy**2)) * (
- b[2:-2, 2:-2] * dx**2 * dy**2
+ (p_old[2:-2, 3:-1] + p_old[2:-2, 1:-3]) * dy**2
+ (p_old[3:-1, 2:-2] + p_old[1:-3, 2:-2]) * dx**2
)
# 速度修正
u[2:-2, 2:-2] = u_hat[2:-2, 2:-2] - dt * p_x
v[2:-2, 2:-2] = v_hat[2:-2, 2:-2] - dt * p_y
4.1.2 投影方法的数学基础
投影方法基于Helmholtz分解定理：任意向量场都可以分解为无旋部分和无散度部分：
F = ∇φ + ∇×ψ
对于不可压缩流动，速度场必须满足 ∇·u = 0，投影方法正是利用这一约束条件。
4.2 有限差分格式
4.2.1 对流项处理
代码使用4阶中心差分 + 人工粘性来处理对流项，这是处理高雷诺数流动的有效方法。
4阶中心差分公式：
∂u/∂x ≈ (-u[i+2] + 8u[i+1] - 8u[i-1] + u[i-2])/(12Δx)
人工粘性项：
β|u|Δx³(∂⁴u/∂x⁴) ≈ β|u|Δx³(u[i+2] - 4u[i+1] + 6u[i] - 4u[i-1] + u[i-2])/Δx⁴
在代码中的具体实现：
# u方向对流项：u∂u/∂x
u_u_x = u_old[2:-2, 2:-2] * (- u_old[2:-2, 4:] + 8. * u_old[2:-2, 3:-1]
- 8. * u_old[2:-2, 1:-3] + u_old[2:-2, :-4]) / (12. * dx) \
+ beta * np.abs(u_old[2:-2, 2:-2]) * dx**3 * (u_old[2:-2, 4:] - 4. * u_old[2:-2, 3:-1]
+ 6. * u_old[2:-2, 2:-2] - 4. * u_old[2:-2, 1:-3] + u_old[2:-2, :-4]) / dx**4
# v方向对流项：v∂u/∂y
v_u_y = v_old[2:-2, 2:-2] * (- u_old[4:, 2:-2] + 8. * u_old[3:-1, 2:-2]
- 8. * u_old[1:-3, 2:-2] + u_old[:-4, 2:-2]) / (12. * dy) \
+ beta * np.abs(v_old[2:-2, 2:-2]) * dy**3 * (u_old[4:, 2:-2] - 4. * u_old[3:-1, 2:-2]
+ 6. * u_old[2:-2, 2:-2] - 4. * u_old[1:-3, 2:-2] + u_old[:-4, 2:-2]) / dy**4
# 总对流项
conv_u = u_u_x + v_u_y
人工粘性的物理意义：
-
抑制高频数值振荡 -
模拟亚网格尺度的耗散效应 -
保持数值稳定性
4.2.2 扩散项处理
扩散项使用标准的2阶中心差分：
∂²u/∂x² ≈ (u[i+1] - 2u[i] + u[i-1])/Δx²
代码实现：
# 二阶中心差分计算拉普拉斯算子
u_xx = (u_old[2:-2, 3:-1] - 2. * u_old[2:-2, 2:-2] + u_old[2:-2, 1:-3]) / dx**2
u_yy = (u_old[3:-1, 2:-2] - 2. * u_old[2:-2, 2:-2] + u_old[1:-3, 2:-2]) / dy**2
lap_u = u_xx + u_yy # ∇²u
diff_u = k * lap_u # 扩散项 = (1/Re) * ∇²u
4.2.3 压力泊松方程求解
使用5点差分格式进行迭代求解：
∇²p ≈ (p[i+1,j] + p[i-1,j] + p[i,j+1] + p[i,j-1] - 4p[i,j])/(Δx²)
代码中的Gauss-Seidel迭代格式：
for it in range(0, it_max+1):
p_old = p.copy()
# 5点差分格式
p[2:-2, 2:-2] = 1. / (2. * (dx**2 + dy**2)) \
* (- b[2:-2, 2:-2] * dx**2 * dy**2
+ (p_old[2:-2, 3:-1] + p_old[2:-2, 1:-3]) * dy**2
+ (p_old[3:-1, 2:-2] + p_old[1:-3, 2:-2]) * dx**2)
# 压力边界条件（Neumann边界条件）
p[1, :] = p[2, :] # ∂p/∂y = 0 在 y = 0
p[-2, :] = p[-3, :] # ∂p/∂y = 0 在 y = Ly
p[:, 1] = p[:, 2] # ∂p/∂x = 0 在 x = 0
p[:, -2] = p[:, -3] # ∂p/∂x = 0 在 x = Lx
p[2, 2] = 0. # 固定参考点消除压力常数
# 检查收敛性
p_res = np.sqrt(np.sum((p - p_old)**2)) / np.sqrt(np.sum(p_old**2))
if p_res < p_tol:
break
4.3 数值稳定性条件
4.3.1 CFL条件（对流稳定性）
对于显式时间积分，必须满足Courant-Friedrichs-Lewy条件：
Δt₁ = h/U
其中 h = min(dx, dy) 是最小网格尺寸。
4.3.2 扩散稳定性条件
对于扩散项的显式处理，稳定性条件为：
Δt₂ = h²/(2ν) = h²Re/2
4.3.3 实际时间步长选择
代码中采用保守的时间步长：
U = 1. # 特征速度
h = min(dx, dy) # 最小网格尺寸
k = 1. / Re # 运动粘度
dt1 = h / U # 对流时间步长
dt2 = 1. / 2. * h**2 / k # 扩散时间步长
dt = min(dt1, dt2) * 0.2 # 安全系数0.2
安全系数的作用：
-
确保数值稳定性 -
考虑非线性项的影响 -
为人工粘性项预留安全余量
4.4 边界条件实施
4.4.1 Dirichlet边界条件
代码采用虚拟网格点（ghost points）方法实施边界条件：
# 物理边界点
u[0, :], v[0, :] = 0., 0. # 底壁面：无滑移条件
u[-1, :], v[-1, :] = 1., 0. # 顶壁面：u=1（驱动速度），v=0
u[:, 0], v[:, 0] = 0., 0. # 左壁面：无滑移条件
u[:, -1], v[:, -1] = 0., 0. # 右壁面：无滑移条件
# 虚拟边界点（确保高阶格式的实施）
u[1, :], v[1, :] = 0., 0.
u[-2, :], v[-2, :] = 1., 0.
u[:, 1], v[:, 1] = 0., 0.
u[:, -2], v[:, -2] = 0., 0.
4.4.2 压力边界条件
压力采用自然边界条件（Neumann边界条件）：
∂p/∂n = 0 （在所有壁面上）
这确保了压力场的唯一性（除了一个常数）。
4.5 网格设计
4.5.1 网格生成
Lx, Ly = 1., 1. # 计算域尺寸
dx, dy = 5e-3, 5e-3 # 网格间距
x = np.arange(0. - dx, Lx + 2. * dx, dx) # 包含虚拟点的x坐标
y = np.arange(0. - dy, Ly + 2. * dy, dy) # 包含虚拟点的y坐标
4.5.2 网格质量评估
代码中包含了Kolmogorov尺度的估算：
eta = Lx / Re**(3. / 4.) # Kolmogorov长度尺度估计
print(f"dx: {dx:.3e}, eta: {eta:.3e}")
物理意义：
-
η 是湍流中最小涡旋的特征尺度 -
dx < η 确保能够分辨最小的流动结构 -
对于 Re = 1000，η ≈ 5.6×10⁻⁴，而 dx = 5×10⁻³
4.6 收敛判断准则
4.6.1 速度收敛判断
u_res = np.sqrt(np.sum((u - u_old)**2)) / np.sqrt(np.sum(u_old**2))
这是相对L²范数误差，物理意义为速度场变化的相对幅度。
4.6.2 压力收敛判断
p_res = np.sqrt(np.sum((p - p_old)**2)) / np.sqrt(np.sum(p_old**2))
4.6.3 收敛标准设置
u_tol = 1e-6 # 速度场收敛容差
p_tol = 1e-5 # 压力场收敛容差
这些值的选择平衡了计算精度和效率。
5. 代码逐行解释和维度分析
5.1 导入库和设置参数
import os
import numpy as np
import matplotlib.pyplot as plt
# 设置图形参数
plt.rcParams["mathtext.fontset"] = "cm" # 使用Computer Modern字体
plt.rcParams["legend.framealpha"] = 1. # 图例背景透明度
plt.rcParams["savefig.dpi"] = 300 # 保存图形的分辨率
说明：
-
os：用于文件系统操作和目录管理 -
numpy：核心数值计算库 -
matplotlib：绘图和可视化库 -
plt.rcParams：设置全局绘图参数，确保输出图形的专业质量
5.2 物理和数值参数设置
Re = 1000. # 雷诺数 [无量纲]
Lx, Ly = 1., 1. # 计算域尺寸 [L]
dx, dy = 5e-3, 5e-3 # 网格间距 [L]
维度分析：
-
Re：无量纲数，Re = ρUL/μ -
Lx, Ly：长度量纲 [L] -
dx, dy：长度量纲 [L]
网格分辨率分析：
-
总网格点数：Nx × Ny = 204 × 204 = 41,616 -
网格长宽比：dx/Lx = dy/Ly = 0.005（0.5%） -
这是一个相对精细的网格，能够捕捉Re=1000时的主要流动特征
5.3 网格生成
x = np.arange(0. - dx, Lx + 2. * dx, dx) # 包含边界的x坐标 [L]
y = np.arange(0. - dy, Ly + 2. * dy, dy) # 包含边界的y坐标 [L]
Nx, Ny = len(x), len(y) # 网格点数 [无量纲]
X, Y = np.meshgrid(x, y) # 二维网格 [L]
网格设计分析：
-
采用均匀结构化网格，便于高阶差分格式实施 -
包含两层虚拟边界点（ghost points）： -
内层虚拟点：用于边界条件实施 -
外层虚拟点：用于4阶差分格式 -
实际计算域：x[2:-2], y[2:-2]
网格坐标范围：
x ∈ [-0.005, 0, 0.005, ..., 0.995, 1.000, 1.005]
y ∈ [-0.005, 0, 0.005, ..., 0.995, 1.000, 1.005]
5.4 Kolmogorov尺度估算
eta = Lx / Re**(3. / 4.)
print(f"dx: {dx:.3e}, eta: {eta:.3e}")
物理意义：
-
eta：Kolmogorov长度尺度的粗略估计 -
对于Re=1000：η ≈ 1.0 / 1000^(3/4) ≈ 5.6×10⁻⁴ -
网格分辨率：dx/η ≈ 5×10⁻³ / 5.6×10⁻⁴ ≈ 9 -
这表明网格可能无法完全分辨最小的湍流尺度
5.5 时间步长确定
U = 1. # 特征速度 [LT⁻¹]
h = min(dx, dy) # 最小网格尺寸 [L]
k = 1. / Re # 运动粘度 [L²T⁻¹]
dt1 = h / U # 对流时间步长 [T]
dt2 = 1. / 2. * h**2 / k # 扩散时间步长 [T]
dt = min(dt1, dt2) * 0.2 # 安全系数0.2 [T]
print(f"dt1: {dt1:.3e}, dt2: {dt2:.3e}, dt: {dt:.3e}")
维度分析验证：
-
dt1 = [L]/[LT⁻¹] = [T]✓ -
dt2 = [L²]/[L²T⁻¹] = [T]✓ -
k = 1/Re = μ/(ρUL) = [ML⁻¹T⁻¹]/[ML⁻³·LT⁻¹·L] = [L²T⁻¹]✓
数值结果示例：
dt1: 5.000e-03 (对流限制)
dt2: 1.250e-05 (扩散限制) ← 控制因子
dt: 2.500e-06 (实际时间步长)
可见，对于Re=1000，扩散限制比对流限制更严格。
5.6 人工粘性和收敛参数
beta = 1. / 4. # 人工粘性系数 [无量纲]
u_tol = 1e-6 # 速度收敛容差 [无量纲]
p_tol = 1e-5 # 压力收敛容差 [无量纲]
it_max = int(1e3) # 压力迭代最大次数
人工粘性参数选择：
-
beta = 0.25：经验最优值 -
过小：无法抑制数值振荡 -
过大：过度耗散，影响解的精度
5.7 结果输出目录设置
dir_res = f"Re{Re:.0f}"
dir_fig_velocity_norm = os.path.join(dir_res, "velocity_norm")
dir_fig_velocity_u = os.path.join(dir_res, "velocity_u")
dir_fig_velocity_v = os.path.join(dir_res, "velocity_v")
# ... 其他目录设置
# 创建目录
os.makedirs(dir_res, exist_ok=True)
os.makedirs(dir_fig_velocity_norm, exist_ok=True)
# ... 其他目录创建
输出组织结构：
Re1000/
├── velocity_norm/ # 速度模长分布
├── velocity_u/ # u速度分量
├── velocity_v/ # v速度分量
├── divergence_hat/ # 中间速度散度
├── pressure/ # 压力分布
├── divergence/ # 最终速度散度
├── vorticity/ # 涡度分布
├── u/ # u速度中心线分布
└── v/ # v速度中心线分布
5.8 变量初始化
u = np.zeros(shape=(Ny, Nx)) + 1e-3 # u速度 [LT⁻¹]
v = np.zeros_like(u) + 1e-3 # v速度 [LT⁻¹]
p = np.zeros_like(u) + 1e-3 # 压力 [L²T⁻²]
b = np.zeros_like(u) + 1e-3 # 压力泊松方程右端项 [T⁻²]
初始化策略：
-
添加小的非零初值（1e-3）避免数值奇异 -
所有变量采用相同的数组形状：(Ny, Nx) -
注意：numpy数组索引为[行, 列] = [y, x]
5.9 主时间循环结构
n = 0 # 时间步计数器
t = 0. # 当前时间
u_res = 9999. # 初始残差
while u_res > u_tol:
u_old = u.copy() # 保存前一时间步的值
v_old = v.copy()
# 预测步、投影步、边界条件...
n += 1
t += dt
u_res = np.sqrt(np.sum((u - u_old)**2)) / np.sqrt(np.sum(u_old**2))
if u_res < u_tol:
print(f" >>> main converged")
break
if t > 200.:
print(f" >>> maximum simulation time reached")
break
5.10 对流项计算详解
5.10.1 u-动量方程的对流项
# u∂u/∂x项
u_u_x = u_old[2:-2, 2:-2] * (- u_old[2:-2, 4:] + 8. * u_old[2:-2, 3:-1]
- 8. * u_old[2:-2, 1:-3] + u_old[2:-2, :-4]) / (12. * dx) \
+ beta * np.abs(u_old[2:-2, 2:-2]) * dx**3 * (u_old[2:-2, 4:] - 4. * u_old[2:-2, 3:-1]
+ 6. * u_old[2:-2, 2:-2] - 4. * u_old[2:-2, 1:-3] + u_old[2:-2, :-4]) / dx**4
分解分析：
-
4阶中心差分部分：
u_old[2:-2, 2:-2] * (- u_old[2:-2, 4:] + 8. * u_old[2:-2, 3:-1]
- 8. * u_old[2:-2, 1:-3] + u_old[2:-2, :-4]) / (12. * dx) -
系数：[-1, 8, -8, 1] / (12dx) -
模板：u[i-2], u[i-1], u[i+1], u[i+2] -
人工粘性部分：
beta * np.abs(u_old[2:-2, 2:-2]) * dx**3 * (u_old[2:-2, 4:] - 4. * u_old[2:-2, 3:-1]
+ 6. * u_old[2:-2, 2:-2] - 4. * u_old[2:-2, 1:-3] + u_old[2:-2, :-4]) / dx**4 -
系数：[1, -4, 6, -4, 1] / dx⁴ -
这是4阶差分算子，用于计算∂⁴u/∂x⁴
5.10.2 人工粘性的量纲分析
人工粘性项 = β|u|Δx³(∂⁴u/∂x⁴)
量纲：[无量纲] × [LT⁻¹] × [L³] × [LT⁻¹L⁻⁴] = [LT⁻²]
与加速度项 ∂u/∂t 的量纲 [LT⁻²] 一致 ✓
5.10.3 v∂u/∂y项
v_u_y = v_old[2:-2, 2:-2] * (- u_old[4:, 2:-2] + 8. * u_old[3:-1, 2:-2]
- 8. * u_old[1:-3, 2:-2] + u_old[:-4, 2:-2]) / (12. * dy)
索引解释：
-
u_old[4:, 2:-2]：u[i+2, j] -
u_old[3:-1, 2:-2]：u[i+1, j] -
u_old[1:-3, 2:-2]：u[i-1, j] -
u_old[:-4, 2:-2]：u[i-2, j]
5.11 扩散项计算
# 拉普拉斯算子的x方向分量
u_xx = (u_old[2:-2, 3:-1] - 2. * u_old[2:-2, 2:-2] + u_old[2:-2, 1:-3]) / dx**2
# 拉普拉斯算子的y方向分量
u_yy = (u_old[3:-1, 2:-2] - 2. * u_old[2:-2, 2:-2] + u_old[1:-3, 2:-2]) / dy**2
# 总拉普拉斯算子
lap_u = u_xx + u_yy # ∇²u [T⁻¹]
# 扩散项
diff_u = k * lap_u # (1/Re) × ∇²u [LT⁻²]
维度验证：
-
u_xx：[LT⁻¹]/[L²] = [L⁻¹T⁻¹] -
lap_u：[L⁻¹T⁻¹] -
diff_u：[L²T⁻¹]·[L⁻¹T⁻¹] = [LT⁻²]✓
5.12 预测步实现
u_hat[2:-2, 2:-2] = u_old[2:-2, 2:-2] + dt * (- conv_u + diff_u)
v_hat[2:-2, 2:-2] = v_old[2:-2, 2:-2] + dt * (- conv_v + diff_v)
物理意义：
-
这是显式欧拉时间积分 -
暂时忽略压力梯度项 -
只考虑对流和扩散效应
维度验证：
-
dt * conv_u：[T] × [LT⁻²] = [LT⁻¹]✓ -
dt * diff_u：[T] × [LT⁻²] = [LT⁻¹]✓
5.13 中间速度散度计算
# 计算中间速度的散度
u_hat_x = (u_hat[2:-2, 3:-1] - u_hat[2:-2, 1:-3]) / (2. * dx)
v_hat_y = (v_hat[3:-1, 2:-2] - v_hat[1:-3, 2:-2]) / (2. * dy)
div_u_hat = u_hat_x + v_hat_y # ∇·û [T⁻¹]
b[2:-2, 2:-2] = div_u_hat / dt # 泊松方程右端项 [T⁻²]
物理意义：
-
div_u_hat：中间速度场的散度 -
理想情况下应该为零（不可压缩条件） -
实际计算中由于数值误差不为零 -
b：压力泊松方程的源项
5.14 压力泊松方程迭代求解
for it in range(0, it_max+1):
p_old = p.copy()
# 5点差分格式
p[2:-2, 2:-2] = 1. / (2. * (dx**2 + dy**2)) \
* (- b[2:-2, 2:-2] * dx**2 * dy**2
+ (p_old[2:-2, 3:-1] + p_old[2:-2, 1:-3]) * dy**2
+ (p_old[3:-1, 2:-2] + p_old[1:-3, 2:-2]) * dx**2)
# 压力边界条件（Neumann边界条件）
p[1, :] = p[2, :] # ∂p/∂y = 0 在 y = 0
p[-2, :] = p[-3, :] # ∂p/∂y = 0 在 y = Ly
p[:, 1] = p[:, 2] # ∂p/∂x = 0 在 x = 0
p[:, -2] = p[:, -3] # ∂p/∂x = 0 在 x = Lx
p[2, 2] = 0. # 固定参考点
# 收敛性检查
p_res = np.sqrt(np.sum((p - p_old)**2)) / np.sqrt(np.sum(p_old**2))
if it % int(it_max / 5) == 0:
print(f" >>> PPE -> it: {it}, p_res: {p_res:.6e}")
if p_res < p_tol:
print(f" >>> PPE converged")
break
算法细节：
-
Gauss-Seidel迭代：使用最新的压力值进行迭代 -
边界条件：自然边界条件（零法向梯度） -
参考点：p[2,2] = 0 固定压力常数 -
收敛监控：每20%的迭代输出一次残差
5.15 投影步（速度修正）
# 计算压力梯度
p_x = (p[2:-2, 3:-1] - p[2:-2, 1:-3]) / (2. * dx) # ∂p/∂x [LT⁻²]
p_y = (p[3:-1, 2:-2] - p[1:-3, 2:-2]) / (2. * dy) # ∂p/∂y [LT⁻²]
# 速度修正
u[2:-2, 2:-2] = u_hat[2:-2, 2:-2] - dt * p_x # [LT⁻¹]
v[2:-2, 2:-2] = v_hat[2:-2, 2:-2] - dt * p_y # [LT⁻¹]
物理意义：
-
从中间速度场中减去压力梯度的贡献 -
确保最终速度场满足不可压缩条件 -
这是投影方法的核心步骤
维度验证：