基于有限体积法（FVM）的MATLAB流体力学求解程序

一、基础框架代码（二维稳态不可压缩流动）

%% 初始化参数
Lx = 0.1; Ly = 0.01; % 计算域尺寸
Nx = 50; Ny = 20;    % 网格数
dx = Lx/Nx; dy = Ly/Ny;% 物理参数
rho = 1.2; mu = 1.8e-5; nu = mu/rho; % 空气物性
Re = 100; U_avg = Re*nu/(2*Ly);      % 雷诺数与平均速度%% 网格与变量定义
x = linspace(dx/2, Lx-dx/2, Nx);
y = linspace(dy/2, Ly-dy/2, Ny);
[X,Y] = meshgrid(x,y);% 初始化场变量（速度、压力）
u = zeros(Nx+1,Ny+2); % x方向速度（面中心）
v = zeros(Nx+2,Ny+1); % y方向速度（面中心）
p = zeros(Nx+2,Ny+2); % 压力（单元中心）%% 边界条件设置
% 入口（左边界）
u(:,1) = U_avg; 
% 出口（右边界）
p(:,end) = 0; 
% 壁面（上下边界）
u(1,:) = 0; u(end,:) = 0;
v(:,1) = 0; v(:,end) = 0;%% 离散参数
alphaU = 0.3; alphaP = 0.2; % 松弛因子
maxIter = 1e4; tol = 1e-5;%% SIMPLE算法主循环
for iter = 1:maxIter% 动量方程离散（x方向）[uStar, F, D] = FVM_xMomentum(u, v, p, rho, dx, dy, mu);% 动量方程离散（y方向）[vStar, G, S] = FVM_yMomentum(u, v, p, rho, dx, dy, mu);% 压力修正方程[pPrime, AP] = FVM_pressureCorr(uStar, vStar, p, dx, dy);% 速度修正[u, v] = FVM_velocityCorrect(uStar, vStar, pPrime, alphaU, alphaP);% 压力更新p = p + alphaP*pPrime;% 收敛判断resU = max(abs(u - uStar));resP = max(abs(pPrime(:)));if max(resU,resP) < tolbreak;end
end%% 后处理
figure;
quiver(squeeze(u(2:end-1,:)), squeeze(v(2:end-1,:)));
title('速度场分布');
xlabel('x'); ylabel('y');%% 核心函数定义
function [uNew, F, D] = FVM_xMomentum(u, v, p, rho, dx, dy, mu)% x方向动量方程离散[Nx,Ny] = size(u);F = zeros(Nx,Ny); D = zeros(Nx,Ny);for i = 2:Nx-1for j = 2:Ny-1% 对流项（中心差分）F(i,j) = 0.5*rho*(u(i,j)*(u(i+1,j)+u(i,j)) + ...v(i,j)*(u(i,j+1)+u(i-1,j)));% 扩散项（中心差分）D(i,j) = mu*( (u(i+1,j)-2*u(i,j)+u(i-1,j))/dx^2 + ...(u(i,j+1)-2*u(i,j)+u(i,j-1))/dy^2 );endend% 构建离散方程A = gallery('poisson', Nx*Ny);b = -D(:) + F(:);uNew = A\b;uNew = reshape(uNew, [Nx,Ny]);
end

---

二、典型应用扩展

1. 加热通道流动（能量方程耦合）

% 能量方程离散
function T = FVM_energy(T, u, v, rho, cp, k, dx, dy)[Nx,Ny] = size(T);for i = 2:Nx-1for j = 2:Ny-1% 对流项conv = rho*cp*(u(i,j)*(T(i+1,j)+T(i,j)) + ...v(i,j)*(T(i,j+1)+T(i-1,j)));% 扩散项diff = k*( (T(i+1,j)-2*T(i,j)+T(i-1,j))/dx^2 + ...(T(i,j+1)-2*T(i,j)+T(i,j-1))/dy^2 );T(i,j) = T(i,j) + (conv - diff)/rho/cp;endend
end

2. 湍流模型集成（k-ε模型）

% 湍动能k方程
function k = FVM_turb_k(k, u, v, mu, rho, dx, dy)% 离散实现（需添加生成项与耗散项）
end% 湍流耗散率ε方程
function epsilon = FVM_turb_epsilon(epsilon, k, u, v, mu, rho, dx, dy)% 离散实现
end

---

三、优化技巧

1. 矩阵预分配：

   A = zeros(Nx*Ny,Nx*Ny);
   b = zeros(Nx*Ny,1);
   

2. 并行计算：

   parfor i = 2:Nx-1% 并行处理每个网格单元
   end
   

3. GPU加速：

   gpu_u = gpuArray(u);
   gpu_v = gpuArray(v);
   % GPU上执行计算
   

参考代码 流体力学中有限体积法的求解程序 www.youwenfan.com/contentcnk/79023.html

四、验证案例

1. 平行板泊肃叶流动

• 理论解：

  umax=2μLΔPH2
  

• 验证方法：对比x=H/2截面速度分布

2. 二维方腔流

• Re=1000：验证自然对流特性
• 收敛标准：残差下降至1e-6

该方法通过模块化设计实现复杂流动问题的数值求解，实际应用中需根据具体问题调整网格划分策略和松弛因子参数。