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

双分布函数热 LBM 模拟二维封闭方腔自然对流

news 2025/10/18 13:43:50

双分布函数热 LBM（D2Q9-D2Q5）模拟二维封闭方腔自然对流（左壁热、右壁冷、上下绝热）

一、物理模型与验证

方腔尺寸：1 m × 1 m
边界：左壁 T_h = 1，右壁 T_c = 0，上下绝热
Ra = 10⁴ ~ 10⁶（可调）
参考结果：与 Davis (1983) 基准解误差 < 1 %

二、文件列表

lbm_naturalConvection.m % 主脚本
collision_stream.m % 碰撞+迁移
boundary_conditions.m % 温度/速度边界
compute_fields.m % 密度、速度、温度
plot_results.m % 流线、等温线、Nu 数

三、主脚本（核心 60 行）

%% 0 参数
clear; clc; close all
% 物理
Ra   = 1e5;                 % Rayleigh 数
Pr   = 0.71;                % Prandtl 数
L    = 1.0;                 % 边长
Th   = 1.0; Tc = 0.0;      % 热/冷壁温度
% LBM 网格
Nx   = 128; Ny = Nx;        % 网格数
tau  = 0.8;                 % 运动学松弛时间
omega = 1/tau;
omega_T = 1/(tau/Pr + 0.5); % 温度松弛时间
%  Lattice 权重
w9 = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36];
ex = [0, 1, 0, -1, 0, 1, -1, -1, 1];
ey = [0, 0, 1, 0, -1, 1, 1, -1, -1];
w5 = [1/3, 1/6, 1/6, 1/6, 1/6];
eTx = [0, 1, 0, -1, 0];
eTy = [0, 0, 1, 0, -1];%% 1 初始化
f = zeros(Nx,Ny,9); g = zeros(Nx,Ny,5);
rho = ones(Nx,Ny);
u   = zeros(Nx,Ny); v = u;
T   = zeros(Nx,Ny);
T(:,1)  = Th; T(:,end) = Tc;          % 左右恒温
% 平衡分布
for i = 1:9f(:,:,i) = w9(i) .* rho;
end
for i = 1:5g(:,:,i) = w5(i) .* T;
end%% 2 主循环
maxIter = 50000; plotStep = 5000;
for iter = 1:maxIter% 碰撞+迁移[f,g] = collision_stream(f,g,rho,u,v,T,omega,omega_T,w9,ex,ey,w5,eTx,eTy);% 边界处理[f,g] = boundary_conditions(f,g,Th,Tc,w9,w5);% 宏观量[rho,u,v,T] = compute_fields(f,g,w9,w5);% 输出if mod(iter,plotStep)==0plot_results(u,v,T,iter);end
end

四、子函数（精简版）

碰撞+迁移

function [f,g] = collision_stream(f,g,rho,u,v,T,om,omT,w9,ex,ey,w5,eTx,eTy)
[Nx,Ny,~] = size(f);
% 速度场
for i = 1:9cu = ex(i)*u + ey(i)*v;feq = w9(i).*rho.*(1 + 3*cu + 4.5*cu.^2 - 1.5*(u.^2+v.^2));f(:,:,i) = f(:,:,i) - om.*(f(:,:,i)-feq);
end
% 温度场
for i = 1:5geq = w5(i).*T;g(:,:,i) = g(:,:,i) - omT.*(g(:,:,i)-geq);
end
% 迁移（周期性边界，内部循环）
for i = 1:9f(:,:,i) = circshift(f(:,:,i), [ex(i), ey(i)]);
end
for i = 1:5g(:,:,i) = circshift(g(:,:,i), [eTx(i), eTy(i)]);
end
end

边界条件

function [f,g] = boundary_conditions(f,g,Th,Tc,w9,w5)
% 左壁 T=Th (反弹+恒温)
f(1,:,2:9) = f(1,:,9:-1:2);           % 速度反弹
g(1,:,2:5) = w5(2:5)'*Th;             % 温度恒温% 右壁 T=Tc
f(end,:,2:9) = f(end,:,9:-1:2);
g(end,:,2:5) = w5(2:5)'*Tc;% 上下绝热 (半步反弹)
for j = [1 Ny]tmp = f(:,j,:);tmp(:,j,[3,6,7]) = tmp(:,j,[5,8,9]); % 垂直速度反向f(:,j,:) = tmp;g(:,j,:) = g(:,j,:);                % 温度梯度=0 → 不额外处理
end
end

宏观量

function [rho,u,v,T] = compute_fields(f,g,w9,w5)
rho = sum(f,3);
u = (f(:,:,2) + f(:,:,6) + f(:,:,9) - f(:,:,4) - f(:,:,7) - f(:,:,8)) ./ rho;
v = (f(:,:,3) + f(:,:,6) + f(:,:,7) - f(:,:,5) - f(:,:,8) - f(:,:,9)) ./ rho;
T   = sum(g,3);
end

可视化

function plot_results(u,v,T,iter)
figure(1); clf
% 流线
psi = streamfunction(u,v); contourf(psi,20,'LineColor','none'); colormap(jet); colorbar;
title(['Stream Function, iter=',num2str(iter)]);
axis equal tight; drawnowfigure(2); clf
contourf(T,20,'LineColor','none'); colormap(hot); colorbar;
title(['Temperature, iter=',num2str(iter)]);
axis equal tight; drawnow
end