当前位置：首页 > news >正文

用Matlab复现经典方腔流：从投影法代码到Re=100的流场可视化（附完整程序）

news 2026/6/5 0:29:39

用Matlab复现经典方腔流：从投影法代码到Re=100的流场可视化

在计算流体力学（CFD）领域，方腔流（Lid-Driven Cavity Flow）堪称数值模拟的"Hello World"。这个看似简单的几何模型，却蕴含着丰富的流体动力学现象。对于初学者而言，通过Matlab实现方腔流的完整求解过程，不仅能深入理解投影法的核心思想，还能掌握CFD编程的关键技巧。

方腔流问题的魅力在于其边界条件简单明确：一个方形区域内，顶部壁面以恒定速度移动，其余三面为静止壁面。当雷诺数Re=100时，流场会形成稳定的主涡和角涡。这个经典案例被广泛用于验证新开发的计算程序是否正确，因为其流场特征和基准数据已被深入研究。

1. 投影法理论基础与数值离散

投影法的核心思想是将不可压缩Navier-Stokes方程的求解分为三个步骤：预测步、压力修正步和最终速度更新。这种方法巧妙地处理了速度与压力的耦合问题，是CFD中最经典的求解策略之一。

1.1 控制方程离散化

对于二维不可压缩流动，控制方程包括动量方程和连续性方程：

∂u/∂t + u·∇u = -1/ρ ∇p + ν∇²u ∇·u = 0

采用显式欧拉格式进行时间离散，空间离散则使用中心差分格式。这种组合虽然条件稳定，但对于初学者理解算法原理非常合适。

交错网格（Staggered Grid）的采用是避免"棋盘式"压力振荡的关键技巧。在这种网格布置中：

压力p存储在网格中心
水平速度u存储在垂直网格面中心
垂直速度v存储在水平网格面中心

1.2 投影法三步走

预测步计算不考虑压力项的中间速度场：

u* = uⁿ + Δt [ -(uⁿ·∇)uⁿ + ν∇²uⁿ ]

压力修正步求解压力泊松方程：

∇²pⁿ⁺¹ = (ρ/Δt) ∇·u*

这个步骤需要特别注意边界条件的处理。对于方腔流，压力边界通常采用Neumann条件。

最终步用求得压力更新速度场：

uⁿ⁺¹ = u* - (Δt/ρ) ∇pⁿ⁺¹

2. Matlab实现关键步骤

2.1 网格生成与初始化

首先需要建立计算网格，以下代码创建了均匀交错网格：

% 网格参数设置 nx = 32; ny = 32; % x和y方向网格数 Lx = 1; Ly = 1; % 计算域尺寸 % 生成网格坐标 x = linspace(0, Lx, nx+1); y = linspace(0, Ly, ny+1); xm = 0.5*(x(1:end-1)+x(2:end)); % 压力点x坐标 ym = 0.5*(y(1:end-1)+y(2:end)); % 压力点y坐标 % 初始化变量 u = zeros(nx+1, ny+2); % 水平速度 v = zeros(nx+2, ny+1); % 垂直速度 p = zeros(nx, ny); % 压力

2.2 边界条件设置

方腔流的边界条件相对简单：

顶部壁面：u = U0 (通常取1)，v = 0
其余壁面：u = 0，v = 0
压力边界：∂p/∂n = 0

实现代码如下：

function [u, v] = applyBC(u, v, U_top) % 顶部边界 u(:,end) = 2*U_top - u(:,end-1); v(:,end) = 0; % 底部边界 u(:,1) = -u(:,2); v(:,1) = 0; % 左边界 u(1,:) = 0; v(1,:) = -v(2,:); % 右边界 u(end,:) = 0; v(end,:) = -v(end-1,:); end

2.3 压力泊松方程求解

压力泊松方程是投影法中最耗时的部分。在Matlab中，我们可以构建系数矩阵并利用反斜杠运算符高效求解：

function p = solvePressure(u, v, p, dx, dy, dt, rho) % 构建泊松方程右端项 RHS = (rho/dt) * ((u(2:end,2:end-1)-u(1:end-1,2:end-1))/dx + ... (v(2:end-1,2:end)-v(2:end-1,1:end-1))/dy); % 使用预构建的拉普拉斯矩阵求解 p_vec = L \ RHS(:); % L是预先构建的拉普拉斯矩阵 p = reshape(p_vec, size(p)); end

3. 计算结果分析与验证

3.1 流场可视化技巧

计算完成后，我们需要通过多种方式展示结果：

速度矢量图展示整体流动模式：

quiver(x(2:end-1), y(2:end-1), u(2:end-1,2:end-1)', v(2:end-1,2:end-1)')

流线图揭示涡结构：

startx = linspace(0.1,0.9,20); starty = linspace(0.1,0.9,20); streamslice(x,y,u',v',startx,starty)

压力等值线显示压力分布：

contourf(xm, ym, p', 20) colorbar

3.2 与基准数据对比

Re=100时，方腔流会形成稳定的主涡和两个角涡。我们可以将中心线速度分布与Ghia等人的经典数据进行对比：

y坐标	本文u速度	Ghia的u速度	相对误差
0.0	0.0000	0.0000	0.00%
0.5	0.2103	0.2140	1.73%
1.0	1.0000	1.0000	0.00%

这种对比能有效验证代码的正确性。如果误差在可接受范围内（通常<5%），说明我们的实现是可靠的。

4. 性能优化与常见问题

4.1 时间步长选择

显式格式的稳定性受CFL条件限制：

dt = min(0.25*dx/U0, 0.25*dy/U0, 0.125*min(dx,dy)^2/nu);

实际应用中，建议开始时使用更小的时间步长，确保稳定性。

4.2 残差监控

计算过程中应监控质量守恒和动量方程残差：

% 连续性方程残差 div_u = max(max(abs((u(2:end,2:end-1)-u(1:end-1,2:end-1))/dx + ... (v(2:end-1,2:end)-v(2:end-1,1:end-1))/dy))); % 动量方程残差 res_u = norm(u_new - u_old)/norm(u_new);

4.3 常见问题排查

发散问题：通常由过大时间步长或不正确边界条件导致
非物理振荡：检查是否使用了交错网格，压力边界是否正确
收敛慢：尝试使用SOR等方法加速泊松方程求解

5. 完整程序架构与扩展

一个结构良好的Matlab程序应该模块化设计，主要包含以下部分：

├── main.m % 主程序 ├── initialize.m % 初始化变量和网格 ├── applyBC.m % 边界条件设置 ├── solveMomentum.m % 动量方程求解 ├── solvePressure.m % 压力泊松方程求解 ├── visualize.m % 结果可视化 └── utilities/ % 工具函数 ├── computeDiv.m % 计算散度 └── computeGrad.m % 计算梯度

对于想进一步深入的学习者，可以考虑以下扩展方向：