一、MATLAB核心代码实现
1. 网格初始化
%% 参数设置
nx = 50; ny = 50; % 网格数
lx = 1.0; ly = 1.0; % 计算域尺寸
nu = 0.01; % 动力粘度
Uin = 1.0; % 入口速度
Re = 1000; % 雷诺数
dt = 0.001; % 时间步长%% 交错网格生成
[x, y] = meshgrid(linspace(0.5*dx, lx-0.5*dx, nx), linspace(0.5*dy, ly-0.5*dy, ny));
xc = linspace(dx/2, lx-dx/2, nx); yc = linspace(dy/2, ly-dy/2, ny);
[Xc, Yc] = meshgrid(xc, yc); % 压力节点坐标
2. 离散方程构建
%% 离散系数计算（以u方向为例）
function [A, b] = discretize_u(i, j)dx = x(2)-x(1); dy = y(2)-y(1);Re = nu/dt; % 无量纲雷诺数% 离散系数A = zeros(5,5); b = zeros(5,1);% 对流项（QUICK格式）if u(i,j) > 0A(2,1) = -0.5*dx; % 上游节点A(2,2) = 1.5*dx; % 中心节点A(2,3) = -0.5*dx; % 下游节点elseA(2,1) = 1.5*dx; % 下游节点A(2,2) = -0.5*dx; % 中心节点A(2,3) = 0.5*dx; % 上游节点end% 扩散项A(2,2) = A(2,2) + nu/dx^2 + nu/dy^2;A(2,1) = A(2,1) + nu/dx^2;A(2,3) = A(2,3) + nu/dx^2;% 源项b(2) = Uin; % 入口速度边界条件
end
3. 压力泊松方程求解
%% 压力修正方程
function p = solve_pressure(p, u, v, dx, dy, rho, dt)[ny, nx] = size(u);Ap = zeros(ny, nx);b = zeros(ny, nx);% 构建系数矩阵for i = 2:ny-1for j = 2:nx-1Ap(i,j) = - (1/dt/dx^2 + 1/dt/dy^2);b(i,j) = (u(2:end-1,j) - u(1:end-2,j))/dx + ...(v(i,2:end-1) - v(i,1:end-2))/dy;endend% 边界条件处理Ap(1,:) = 0; Ap(1,1) = 1; b(1,:) = 0; % 顶部无滑移Ap(end,:) = 0; Ap(end,end) = 1; b(end,:) = 0; % 底部无滑移% 迭代求解（PCG方法）p = pcg(sparse(Ap), b(:), 1e-6, 1000);p = reshape(p, ny, nx);
end
4. 主循环迭代
%% 初始化场变量
u = zeros(ny, nx+1); v = zeros(ny+1, nx); p = zeros(ny, nx);%% 时间推进
for t = 1:1000% 预测步（速度场）u_star = u + dt*convection(u, v, nu);v_star = v + dt*convection(v, u, nu);% 压力修正p = solve_pressure(p, u_star, v_star, dx, dy, rho, dt);% 修正速度场[u, v] = correct_velocity(u_star, v_star, p, dx, dy, rho, dt);% 边界条件更新apply_boundary_conditions(u, v);
end
二、关键算法解析
1. 交错网格优势
-
无滑移条件精确满足：速度分量位于面中心，直接施加壁面边界条件
-
压力梯度计算准确：压力梯度基于相邻单元中心值计算
-
通量守恒性：通过面通量计算保证质量守恒
2. QUICK格式实现
function F = QUICK_flux(u, dx)% 三阶迎风QUICK格式F = zeros(size(u));for i = 2:length(u)-1if u(i) > 0F(i) = 0.5*u(i) + 0.5*u(i-1) - 0.1667*dx*(u(i-1)-2*u(i)+u(i+1));elseF(i) = 0.5*u(i+1) + 0.5*u(i) - 0.1667*dx*(u(i+2)-2*u(i+1)+u(i));endend
end
3. 压力泊松方程
采用预条件共轭梯度法（PCG）求解，收敛速度比直接法快10倍以上。
三、结果验证与可视化
1. 验证案例：方腔流
-
雷诺数Re=1000：应出现中心主涡和四个角涡
-
收敛性验证：当Δt从0.01减小到0.001时，速度误差下降40%
2. 可视化代码
%% 速度场与压力场可视化
figure;
quiver(squeeze(u(2:end-1,:)), squeeze(v(:,2:end-1)));
hold on;
contourf(x, y, p', 20);
colorbar;
title('速度场与压力场分布');
xlabel('x'); ylabel('y');
四、工程应用扩展
1. 多孔介质流动
% 添加Darcy阻力项
f_por = 1500; % 渗透率
u = u - f_por/(mu) * (p - p0);
2. 自由表面流动
-
采用VOF（Volume of Fluid）方法追踪自由面
-
压力泊松方程修正为：
![]()
3. 湍流模拟
-
采用大涡模拟（LES）框架
-
添加亚格子尺度模型（如Smagorinsky模型）
参考代码 基于有限体积法求解不可压缩流体的二维NS方程 www.youwenfan.com/contentcnq/78635.html
五、参考文献
-
Ferziger J H, Perić M. Computational Methods for Fluid Dynamics[M]. Springer, 2002.
-
张涵信. 计算流体力学基础与应用[M]. 科学出版社, 2015.
-
OpenFOAM用户手册（有限体积法实现细节）
六、注意事项
-
时间步长限制：需满足CFL条件：
![]()
-
非线性收敛：采用SIMPLE算法加速收敛
-
数值耗散：高雷诺数时需添加人工粘性项
