%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A LATTICE BOLTZMANN FORMULATION FOR THE ANALYSIS OF
% RADIATIVE HEAT TRANSFER PROBLEMS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lattice Boltzmann sample, written in matlab
% Copyright (C) 2010-12-10 14:53:49 Zhao Lei .
% Address:
% E-mail: haorizi2008.cool@163.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clf
% RADIATION CONSTANS
ly = 51;
lx = 501;
delta_x= 1./(lx-1);
delta_y= 1./(ly-1);
delta_t= 1.0e-004;
sigma = 5.67e-008; % constant s is the geometric distance
Thot = 1;
Tcold = 0;
Ib = sigma*Thot.^4/pi ; % blackbody intensity
G = 4*pi.*Ib ; % volumetric absorption
IEq = G/(4*pi); % equilibrinm PDFs
beta = 1.0; % extinction coefficient
n =8;
alpha = pi/4; %
tPlot = 100; % iterations between successive graphical outputs
tStatistics = 10; % iterations between successive file accesses
Kn = delta_x*beta <= 0.05; % Kn number is smaller than a threshold value that ensures stability and accuracy
% D2Q8 LATTICE CONSTANTS
cx = [ 1, 0, -1, 0, 1, -1, -1, 1 ];
cy = [ 0, 1, 0, -1, 1, 1, -1, -1 ];
t = ones(1,8)*alpha/(2*pi); % weight
az = [0:2*pi/n:2*pi];
[y,x] = meshgrid(1:ly,1:lx);
ux = delta_x/delta_t; % velocity
uy = delta_y/delta_t;
omega = sqrt((cx*ux)*(cx*ux)'+(cx*uy)*(cx*uy)')*beta; % relaxation frequency
% INITIAL CONDITION FOR RADIATION:
IIn = ones(8,lx,ly); %8 is the direction number,2 is r
IOut = ones(8,lx,ly);
ts = 0;
% MAIN LOOP
while (ts<80000 && Kn <= 0.05 && abs((abs(IOut)-abs(IIn))./abs(IIn)))<1e-6 || ts<100
% BOUNDARY CONDITIONS FOR RADIATION:
for i=1:8
IIn (i , : ,1) = sigma*Thot^4/pi; % sigma is Stefan-Boltzmann constant
IOut(i , 1 ,:) = 0;
IOut(i ,lx ,:) = 0;
IOut(i ,: ,ly) = 0;
end
% COLLISION STEP
for i=1:8 % IEq = zeros(1,lx,ly); %让平衡分布函数IEq初始为同维的零矩阵
IEq (i,:,:) = IEq(i,:,:)+ reshape(t*reshape(IIn,8,lx*ly),1,lx,ly); % equilibrinm PDFs % IOut(i,:,:) = zeros(1,lx,ly);
IOut(i,:,:) = IIn(i,:,:) - omega .* (IIn(i,:,:)-IEq(i,:,:));
end
% STREAMING STEP RADIATION
for i=1:8
IIn(i,:,:) = circshift(IOut(i,:,:), [0,cx(i),cy(i)]);
tx = pi*sin(alpha/2)*cos(az(i)); %weight x direction
ty = pi*sin(alpha/2)*sin(az(i)); %weight y direction
end
ts = ts+1;
end
% radiative heat fluxes
Qx = sum(tx.*IOut);
Qy = sum(ty.*IOut);
% VISUALIZATION
if (mod(ts,tPlot)==0)
Q = reshape(sqrt(Qx.^2+Qy.^2),lx,ly);
imagesc(Q');colorbar
title('heat flux');
axis equal off; drawnow
end