% LBM 3D - D3Q15 Heat Diffusion in a Plate
% Problem - 5.10
clc
clear variables
close all
alpha = 0.25;
Cs = 1/sqrt(3);
omega = 1/(alpha/Cs^2 + 0.5);
% Domain Discretization
L = 40; M = 40; N = 40;
xE = 40; yE = 40; zE = 40;
dx = 1; dy = 1; dz = 1;
x = 0:dx:xE;
y = 0:dy:yE;
z = 0:dz:zE;
tfinal = 200;
dt = 1.0;
twall = 1.0;
% Weights
w0 = 2/9;
w = [1/9,1/9,1/9,1/9,1/9,1/9,1/72,1/72,1/72,1/72,1/72,1/72,1/72,1/72];
% Distribution Functions
f0 = zeros(L+1,M+1,N+1);
f0eq = zeros(L+1,M+1,N+1);
f = zeros(L+1,M+1,N+1,14);
feq = zeros(L+1,M+1,N+1,14);
T = zeros(L+1,M+1,N+1);
% Initially When system is at rest
for i = 1:L+1
for j = 1:M+1
for k = 1:N+1
f0(i,j,k) = w0*T(i,j,k);
for p = 1:14
f(i,j,k,p) = w(p)*T(i,j,k);
end
end
end
end
%-----------------------------------------------------------------------------------------
% LBM Simulation
for t = 1:tfinal
% Collision
for i = 1:L+1
for j = 1:M+1
for k = 1:N+1
f0eq(i,j,k) = w0*T(i,j,k);
f0(i,j,k) = (1-omega)*f0(i,j,k) + omega*f0eq(i,j,k);
for p = 1:14
feq(i,j,k,p) = w(p)*T(i,j,k);
f(i,j,k,p) = (1-omega)*f(i,j,k,p) + omega*feq(i,j,k,p);
end
end
end
end
% Streaming
f(:,:,:,1) = circshift(squeeze(f(:,:,:,1)),[+1,0,0]);
f(:,:,:,2) = circshift(squeeze(f(:,:,:,2)),[-1,0,0]);
f(:,:,:,3) = circshift(squeeze(f(:,:,:,3)),[0,+1,0]);
f(:,:,:,4) = circshift(squeeze(f(:,:,:,4)),[0,-1,0]);
f(:,:,:,5) = circshift(squeeze(f(:,:,:,5)),[0,0,+1]);
f(:,:,:,6) = circshift(squeeze(f(:,:,:,6)),[0,0,-1]);
f(:,:,:,7) = circshift(squeeze(f(:,:,:,7)),[+1,+1,+1]);
f(:,:,:,8) = circshift(squeeze(f(:,:,:,8)),[-1,-1,-1]);
f(:,:,:,9) = circshift(squeeze(f(:,:,:,9)),[+1,+1,-1]);
f(:,:,:,10) = circshift(squeeze(f(:,:,:,10)),[-1,-1,+1]);
f(:,:,:,11) = circshift(squeeze(f(:,:,:,11)),[-1,+1,-1]);
f(:,:,:,12) = circshift(squeeze(f(:,:,:,12)),[+1,-1,+1]);
f(:,:,:,13) = circshift(squeeze(f(:,:,:,13)),[-1,+1,+1]);
f(:,:,:,14) = circshift(squeeze(f(:,:,:,14)),[+1,-1,-1]);
% Boundary Conditions
% Left Wall
for j = 1:M+1
for k = 1:N+1
f(1,j,k,1) = w(1)*twall + w(2)*twall - f(1,j,k,2);
f(1,j,k,5) = w(5)*twall + w(6)*twall - f(1,j,k,6);
f(1,j,k,7) = w(7)*twall + w(8)*twall - f(1,j,k,8);
f(1,j,k,9) = w(9)*twall + w(10)*twall - f(1,j,k,10);
f(1,j,k,12) = w(12)*twall + w(11)*twall - f(1,j,k,11);
f(1,j,k,14) = w(14)*twall + w(13)*twall - f(1,j,k,13);
end
end
% Bottom Wall (Bounce Back) {Adiabatic Boundary Conditions)
for i = 1:L+1
for k = 1:N+1
f(i,1,k,3) = f(i,2,k,3);
f(i,1,k,5) = f(i,2,k,5);
f(i,1,k,7) = f(i,2,k,7);
f(i,1,k,9) = f(i,2,k,9);
f(i,1,k,11) = f(i,2,k,11);
f(i,1,k,13) = f(i,2,k,13);
end
end
% Right Wall (Constant Temperature T = 0.0)
for j= 1:M+1
for k = 1:N+1
f(L+1,j,k,2) = -f(L+1,j,k,1);
f(L+1,j,k,6) = -f(L+1,j,k,5);
f(L+1,j,k,8) = -f(L+1,j,k,7);
f(L+1,j,k,10) = -f(L+1,j,k,9);
f(L+1,j,k,11) = -f(L+1,j,k,12);
f(L+1,j,k,13) = -f(L+1,j,k,14);
end
end
% Top Wall (Constant Temperature T = 0.0)
for i = L+1
for k = 1:N+1
f(i,M+1,k,4) = -f(i,M+1,k,3);
f(i,M+1,k,6) = -f(i,M+1,k,5);
f(i,M+1,k,8) = -f(i,M+1,k,7);
f(i,M+1,k,10) = -f(i,M+1,k,9);
f(i,M+1,k,12) = -f(i,M+1,k,11);
f(i,M+1,k,14) = -f(i,M+1,k,13);
end
end
% Front Wall (Constant Temperature T = 0.0)
for i = 1:L+1
for j = 1:M+1
f(i,j,1,5) = -f(i,j,1,6);
f(i,j,1,7) = -f(i,j,1,8);
f(i,j,1,10) = -f(i,j,1,9);
f(i,j,1,12) = -f(i,j,1,11);
f(i,j,1,13) = -f(i,j,1,14);
% f(i,j,1,5) = -f(i,j,1,6);
end
end
% Backward Wall(Constant Temperature T = 0.0)
for i = 1:L+1
for j = 1:M+1
f(i,j,N+1,6) = -f(i,j,N+1,5);
f(i,j,N+1,8) = -f(i,j,N+1,7);
f(i,j,N+1,9) = -f(i,j,N+1,10);
f(i,j,N+1,11) = -f(i,j,N+1,12);
f(i,j,N+1,14) = -f(i,j,N+1,13);
end
end
% Final Temperature
for i = 1:L+1
for j = 1:M+1
for k = 1:N+1
sum = 0.0;
for p = 1:14
sum = sum + f(i,j,k,p);
end
T(i,j,k) = f0(i,j,k) + sum;
end
end
end
end
figure
slice(x,y,z,T,0:8:40,[],[],'nearest')
figure
slice(x,y,z,T,[],5,[],'nearest')
figure
slice(x,y,z,T,[10 30],[],20,'nearest')
