First of all I am just a starter in Matlab, so the program is not optimized. For now I just want to get it running. This is a 3D heat problem with boundary conditions. Fixed temperature from the north, insulated from the bottom, and we have convection on the 4 surrounding sides with the same h. when I run it it does not give me the right convergences. I have the analytical and the solidworks approximations, so I know the results I get in Matlab are wrong. Could someone advise me please. Your help is appreciated.
==============================
%===Solution of the problem of 2D in the bookmarks:Case2a=====
clc;
close all;
clear all;
DX=0.01; % step size
DY=0.01;
DZ=0.01;
Lx= 0.6; %length along x-axis in cmn
Ly=0.4; %length along y-axis
Lz=0.3; %length along z-axis
Nx=Lx/DX+1; %Number of nodes along x-axis
Nx=floor(Nx)+1; % I shouldn't add 1. I should get it automatically
Ny=(Ly/DY+1); %Number of nodes along y-axis
Ny=floor(Ny);
Nz=(Lz/DZ+1);
Nz=floor(Nz)+1;
K=2;
h=500;
T_inf=20;
X=0:DX:Lx;
Y=0:DY:Ly;
Z=0:DZ:Lz;
T=zeros(Nx,Ny,Nz);
% I added imaginary nodes outside then it gave me numbers
tol=1;
s=0; %s=2000 could be set as the maximum number of allowed iteration for
example
err=1;
T_old=1;
while err >=tol %&& s<2001
s=s+1;
%------------boundary conditions------------------------------------------
%====North temperature (Boundary)
T(:,Ny,:)=90;
%===South Boundary "insulation" ==============
for k=1:Nz
for i=2:Nx-1
T(i,1,k)=0.25*[2*T(i,2,k)+T(i-1,1,k)+T(i+1,1,k)];
%T(i-1,1,k)=4* T(i,1,k)-2*T(i,2,k)-T(i+1,1,k);
end
end
%================Back Boundary "Convection" ================
for j=2:Ny-1 %changed
for i=2:Nx-1
T(i,j,1)=0.5*K/(h*DX+2*K)*[2*T(i,j,2)+T(i-1,j,1)+T(i+1,j,1)+2*h*DX*T_inf/K];
end
end
%================Front Boundary "Convection" ================
for j=2:Ny-1
for i=2:Nx-1
T(i,j,Nz)=0.5*K/(h*DX+2*K)*[2*T(i,j,Nz-1)+T(i-1,j,Nz-1)+T(i+1,j,Nz-1)+2*h*DX*T_inf/K];
end
end
%================West Boundary "Convection" ================
for j=2:Ny-1
for k=2:Nz-1
T(1,j,k)=0.5*K/(h*DX+2*K)*[2*T(2,j,k)+T(1,j,k-1)+T(1,j,k+1)+2*h*DX*T_inf/K];
end
end
%================East Boundary "Convection" ================
for j=2:Ny-1
for k=2:Nz-1
T(Nx,j,k)=0.5*K/(h*DX+2*K)*[2*T(Nx-1,j,k)+T(Nx,j,k-1)+T(Nx,j,k+1)+2*h*DX*T_inf/K];
end
end
% The interior nodes
for k = 2:Nz-1 % I think it needs to go loop by loop(pairs of loops)
for j = 2:Ny-1
for i=2:Nx-1
T(i,j,k)=0.25*(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1));
end
end
end
err=max(max(abs(T-T_old)));
T_old=T;
end
T=rot90(T)