Clear Filters
Clear Filters

How can I plot the Temperature distribution in this code.

4 views (last 30 days)
% 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

Accepted Answer

Stephan
Stephan on 14 May 2021
Edited: Stephan on 14 May 2021
Play arround with slice:
% 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')

More Answers (0)

Categories

Find more on Schedule Model Components in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!