How can I plot the Temperature distribution in this code.
Afficher commentaires plus anciens
% 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
Réponse acceptée
% 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 commentaires
Mahendra Yadav
le 14 Mai 2021
Stephan
le 14 Mai 2021
Did you notice that you can accept and/or vote for useful answers?
Mahendra Yadav
le 14 Mai 2021
Plus de réponses (0)
Catégories
En savoir plus sur Creating, Deleting, and Querying Graphics Objects dans Centre d'aide et File Exchange
le 14 Mai 2021
le 14 Mai 2021
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
