How can I plot the Temperature distribution in this code.
4 views (last 30 days)
Show older comments
Mahendra Yadav
on 14 May 2021
Commented: Mahendra Yadav
on 14 May 2021
% 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
0 Comments
Accepted Answer
Stephan
on 14 May 2021
Edited: Stephan
on 14 May 2021
% 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')
3 Comments
More Answers (0)
See Also
Categories
Find more on Schedule Model Components in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!