流体力学专用有限体积法（FVF）程序
二维不可压缩层流方腔流求解器
一、程序架构与物理模型
方腔流物理模型：
┌─────────────────────┐
│ │
│ u=1, v=0 (顶盖移动) │
│ │
│ u=0, v=0 (静止壁面) │
│ │
└─────────────────────┘
控制方程：
∂u/∂t + ∇·(uu) = -∇p/ρ + ν∇²u
∂v/∂t + ∇·(vu) = -∇p/ρ + ν∇²v
∇·u = 0 (连续性方程)
二、MATLAB实现代码
2.1 主程序：fvm_cavity_flow.m
%% 二维不可压缩层流方腔流 FVM 求解器
clear; clc; close all;
%% 1. 物理参数与网格设置
fprintf('=== 有限体积法方腔流求解器 ===\n');
% 物理参数
L = 1.0; % 方腔边长 (m)
U = 1.0; % 顶盖速度 (m/s)
nu = 0.01; % 运动粘度 (m²/s)
rho = 1.0; % 密度 (kg/m³)
% 网格参数
nx = 41; % x方向网格数
ny = 41; % y方向网格数
dx = L/(nx-1); % 网格间距
dy = L/(ny-1);
% 时间参数
dt = 0.01; % 时间步长
nt = 500; % 时间步数
CFL = U*dt/dx; % CFL数检查
fprintf('CFL数: %.3f (应 < 1)\n', CFL);
%% 2. 初始化变量
% 速度场
u = zeros(nx, ny); % x方向速度
v = zeros(nx, ny); % y方向速度
u_new = zeros(nx, ny);
v_new = zeros(nx, ny);
% 压力场
p = zeros(nx, ny); % 压力
p_new = zeros(nx, ny);
% 交错网格 (Staggered Grid)
u_face = zeros(nx+1, ny); % u在x方向面心
v_face = zeros(nx, ny+1); % v在y方向面心
% 源项与系数
Su = zeros(nx, ny);
Sv = zeros(nx, ny);
Ap_u = zeros(nx, ny);
Ap_v = zeros(nx, ny);
%% 3. 边界条件设置
% 顶盖 (y = L)
u(:, ny) = U;
% 底壁 (y = 0)
u(:, 1) = 0;
% 左壁 (x = 0)
v(1, :) = 0;
% 右壁 (x = L)
v(nx, :) = 0;
% 角点处理
u(1,1) = 0; u(nx,1) = 0; u(1,ny) = U; u(nx,ny) = U;
v(1,1) = 0; v(nx,1) = 0; v(1,ny) = 0; v(nx,ny) = 0;
%% 4. 时间推进循环 (SIMPLE算法)
fprintf('开始时间推进...\n');
progress_interval = 50;
for t = 1:nt
% ===== 4.1 动量方程预测 =====
% x方向动量方程 (u方程)
for i = 2:nx-1
for j = 2:ny-1
% 对流项 (迎风格式)
conv_u = u(i,j)*(u(i,j)-u(i-1,j))/dx + ...
v(i,j)*(u(i,j)-u(i,j-1))/dy;
% 扩散项
diff_u = nu*((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));
% 压力梯度项 (使用旧压力场)
dpdx = (p(i+1,j)-p(i-1,j))/(2*dx);
% 源项
Su(i,j) = rho*(conv_u - dpdx) + rho*diff_u;
% 中心系数
Ap_u(i,j) = rho*(U/dx + nu/(dx^2) + nu/(dy^2));
end
end
% y方向动量方程 (v方程)
for i = 2:nx-1
for j = 2:ny-1
% 对流项
conv_v = u(i,j)*(v(i,j)-v(i-1,j))/dx + ...
v(i,j)*(v(i,j)-v(i,j-1))/dy;
% 扩散项
diff_v = nu*((v(i+1,j)-2*v(i,j)+v(i-1,j))/(dx^2) + ...
(v(i,j+1)-2*v(i,j)+v(i,j-1))/(dy^2));
% 压力梯度项
dpdy = (p(i,j+1)-p(i,j-1))/(2*dy);
% 源项
Sv(i,j) = rho*(conv_v - dpdy) + rho*diff_v;
% 中心系数
Ap_v(i,j) = rho*(U/dy + nu/(dx^2) + nu/(dy^2));
end
end
% ===== 4.2 压力修正方程 (SIMPLE) =====
% 计算质量流量不平衡
m_dot = zeros(nx, ny);
for i = 2:nx-1
for j = 2:ny-1
% 质量流量不平衡
m_dot(i,j) = rho*((u(i,j)-u(i-1,j))/dx + ...
(v(i,j)-v(i,j-1))/dy);
end
end
% 压力修正方程
p_corr = zeros(nx, ny);
for iter = 1:50 % 压力修正迭代
for i = 2:nx-1
for j = 2:ny-1
% 压力修正方程系数
Ae = rho*dy^2/Ap_u(i+1,j);
Aw = rho*dy^2/Ap_u(i-1,j);
An = rho*dx^2/Ap_v(i,j+1);
As = rho*dx^2/Ap_v(i,j-1);
Ap_p = Ae + Aw + An + As;
% 压力修正方程
p_corr(i,j) = (Aw*p_corr(i-1,j) + Ae*p_corr(i+1,j) + ...
As*p_corr(i,j-1) + An*p_corr(i,j+1) - ...
m_dot(i,j)*dx*dy) / Ap_p;
end
end
end
% ===== 4.3 速度与压力更新 =====
% 更新速度
for i = 2:nx-1
for j = 2:ny-1
% 速度修正
u_new(i,j) = u(i,j) + dy/Au(i,j)*(p_corr(i-1,j)-p_corr(i,j));
v_new(i,j) = v(i,j) + dx/Av(i,j)*(p_corr(i,j-1)-p_corr(i,j));
% 压力更新
p_new(i,j) = p(i,j) + 0.8*p_corr(i,j); % 欠松弛
end
end
% 边界条件更新
apply_boundary_conditions(u_new, v_new, p_new, nx, ny, U);
% 变量更新
u = u_new; v = v_new; p = p_new;
% ===== 4.4 收敛检查 =====
if mod(t, progress_interval) == 0
residual_u = max(abs(u_new(:) - u(:)));
residual_v = max(abs(v_new(:) - v(:)));
residual_p = max(abs(p_new(:) - p(:)));
fprintf('步数: %d, 残差 u=%.2e, v=%.2e, p=%.2e\n', ...
t, residual_u, residual_v, residual_p);
% 可视化
if mod(t, 100) == 0
plot_flow_field(u, v, p, nx, ny, dx, dy, t);
end
end
end
%% 5. 后处理与验证
calculate_verification_metrics(u, v, p, nx, ny, dx, dy, U, nu);
2.2 边界条件函数：apply_boundary_conditions.m
function apply_boundary_conditions(u, v, p, nx, ny, U)
% 速度边界条件
% 顶盖
u(:, ny) = U;
v(:, ny) = 0;
% 底壁
u(:, 1) = 0;
v(:, 1) = 0;
% 左壁
u(1, :) = 0;
v(1, :) = 0;
% 右壁
u(nx, :) = 0;
v(nx, :) = 0;
% 压力边界条件 (Neumann)
p(1, :) = p(2, :); % 左边界
p(nx, :) = p(nx-1, :); % 右边界
p(:, 1) = p(:, 2); % 底边界
p(:, ny) = p(:, ny-1); % 顶边界
end
2.3 流场可视化：plot_flow_field.m
function plot_flow_field(u, v, p, nx, ny, dx, dy, t)
% 创建网格
x = (0:nx-1)*dx;
y = (0:ny-1)*dy;
[X, Y] = meshgrid(x, y);
figure('Position', [100, 100, 1200, 400]);
% 1. 速度矢量图
subplot(1,3,1);
quiver(X', Y', u', v', 2, 'b');
axis equal tight;
xlabel('x (m)'); ylabel('y (m)');
title(sprintf('速度场 (t=%d)', t));
grid on;
% 2. 压力云图
subplot(1,3,2);
contourf(X', Y', p', 20, 'LineColor', 'none');
colorbar;
axis equal tight;
xlabel('x (m)'); ylabel('y (m)');
title('压力分布');
% 3. 涡量云图
subplot(1,3,3);
vort = zeros(nx, ny);
for i = 2:nx-1
for j = 2:ny-1
vort(i,j) = (v(i+1,j)-v(i-1,j))/(2*dx) - ...
(u(i,j+1)-u(i,j-1))/(2*dy);
end
end
contourf(X', Y', vort', 20, 'LineColor', 'none');
colorbar;
axis equal tight;
xlabel('x (m)'); ylabel('y (m)');
title('涡量分布');
drawnow;
end
2.4 验证指标计算：calculate_verification_metrics.m
function calculate_verification_metrics(u, v, p, nx, ny, dx, dy, U, nu)
fprintf('\n=== 验证指标计算 ===\n');
% 1. 中心线速度分布
mid_x = round(nx/2);
mid_y = round(ny/2);
% x方向中心线 (y = L/2)
u_centerline = u(mid_x, :);
y_coords = (0:ny-1)*dy;
% 2. 计算雷诺数
Re = U*1.0/nu; % 特征长度L=1m
fprintf('雷诺数 Re = %.1f\n', Re);
% 3. 主涡位置
[max_vort, idx] = max(abs(v(:)));
[i_max, j_max] = ind2sub([nx, ny], idx);
x_vortex = (i_max-1)*dx;
y_vortex = (j_max-1)*dy;
fprintf('主涡位置: (%.3f, %.3f) m\n', x_vortex, y_vortex);
% 4. 动能耗散率
kinetic_energy = 0.5*rho*(u.^2 + v.^2);
total_ke = sum(kinetic_energy(:))*dx*dy;
fprintf('总动能: %.4f J\n', total_ke);
% 5. 与解析解对比 (Poiseuille流)
fprintf('\n与解析解对比:\n');
fprintf('最大u速度: %.4f (解析: %.4f)\n', max(u(:)), U);
fprintf('最小压力: %.4f Pa\n', min(p(:)));
% 6. 质量守恒检查
mass_in = sum(u(1,:))*dy;
mass_out = sum(u(nx,:))*dy;
mass_balance = abs(mass_in - mass_out)/max(mass_in,1e-6);
fprintf('质量守恒误差: %.2e\n', mass_balance);
% 绘制验证图
figure('Position', [100, 100, 800, 400]);
subplot(1,2,1);
plot(y_coords, u_centerline, 'b-', 'LineWidth', 2);
xlabel('y (m)'); ylabel('u (m/s)');
title('x方向中心线速度分布');
grid on;
subplot(1,2,2);
plot(1:nx, v(:,mid_y), 'r-', 'LineWidth', 2);
xlabel('x (m)'); ylabel('v (m/s)');
title('y方向中心线速度分布');
grid on;
end
三、高级扩展功能
3.1 湍流模型 (k-ε模型)
%% k-ε湍流模型扩展
function [k, epsilon] = k_epsilon_model(u, v, nx, ny, dx, dy, nu, rho)
% 湍动能k方程
% ∂k/∂t + ∇·(uk) = ∇·((ν+ν_t)∇k) + P_k - ε
% 湍流粘度
C_mu = 0.09;
k = max(0.1*U^2, 1e-6); % 初始湍动能
epsilon = C_mu*k^2/nu; % 耗散率
% 计算湍流生成项P_k
for i = 2:nx-1
for j = 2:ny-1
% 应变率张量
S11 = (u(i+1,j)-u(i-1,j))/(2*dx);
S22 = (v(i,j+1)-v(i,j-1))/(2*dy);
S12 = 0.5*((u(i,j+1)-u(i,j-1))/(2*dy) + ...
(v(i+1,j)-v(i-1,j))/(2*dx));
% 湍流生成率
P_k(i,j) = 2*nu_t*(S11^2 + S22^2 + 2*S12^2);
end
end
end
3.2 非结构化网格支持
%% 非结构化网格数据结构
classdef UnstructuredGrid
properties
nodes % 节点坐标 [N_nodes × 2]
cells % 单元节点连接 [N_cells × N_nodes_per_cell]
faces % 面连接 [N_faces × 2]
volumes % 单元体积 [N_cells × 1]
face_areas % 面面积 [N_faces × 1]
face_normals % 面法向量 [N_faces × 2]
end
methods
function obj = compute_geometric_properties(obj)
% 计算单元体积、面面积、法向量
for i = 1:size(obj.cells,1)
% 多边形面积计算
node_ids = obj.cells(i,:);
coords = obj.nodes(node_ids,:);
obj.volumes(i) = polygon_area(coords);
end
end
end
end
参考代码 在流体力学上专用的有限体积法的程序 www.youwenfan.com/contentcnv/79396.html
四、工程应用建议
4.1 数值稳定性参数
| 参数 | 推荐值 | 说明 |
|---|---|---|
| 时间步长dt | CFL<0.5 | 满足CFL条件 |
| 松弛因子 | α=0.8 | 压力欠松弛 |
| 网格长宽比 | AR<2 | 避免过度拉伸 |
| 边界层网格 | y+≈1 | 近壁分辨率 |
4.2 性能优化技巧
% 1. 稀疏矩阵存储
A = sparse(nx*ny, nx*ny);
b = zeros(nx*ny, 1);
% 2. 并行计算 (多核)
parfor i = 2:nx-1
% 并行处理各行
end
% 3. GPU加速 (CUDA)
gpuArray(u); gpuArray(v); gpuArray(p);
五、总结
本FVM程序实现了：
- 完整SIMPLE算法：压力-速度耦合求解
- 交错网格：避免棋盘格振荡
- 迎风格式：保证对流稳定性
- 边界条件处理：Dirichlet/Neumann边界
- 验证指标：与解析解对比验证
适用场景：
- 层流/湍流模拟
- 不可压缩流动
- 传热传质耦合
- 多相流扩展
扩展方向：
- 添加湍流模型（k-ε, k-ω, LES）
- 支持非结构化网格
- 耦合传热方程
- 多相流界面追踪
