%% Steady state solver
close all
clear all
clc
%Temperature value at centre with various solver for Rectangular BC
%Input conditions
L = 10;
nx = 11;
ny = nx;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
T = 300*ones(nx,ny);
Term = ((dx^2*dy^2) / (2*(dx^2+dy^2)));
%Assignig periphery anotations
T_L = 400;
T_R = 600;
T_T = 800;
T_B = 900;
% Periphery BC
T(:,1) = 400;
T(:,nx) = 600;
T(1, :) = 800;
T(ny, :) = 900;
% Corner BC
T(1,1) = (T_L + T_T) / 2;
T(nx , 1) = (T_L + T_B) / 2;
T(nx, ny) = (T_R + T_B) / 2;
T(1, ny) = (T_T + T_R) / 2;
%Temperature inlet conditions
T_old = T;
Tol = 1e-4;
Error = 9e9;
%%%% Jacobi method
Thermal_Itr = 3;
Iteration_no = 1;
if (Thermal_Itr == 1)
while (Error > Tol)
for i = 2 : nx-1
for j = 2 : ny-1
T(i,j) = Term*(T_old(i-1,j) + T_old(i+1,j) + T_old(i,j-1) + T_old(i,j+1));
end
end
Error = (max(max(abs(T - T_old))));
T_old = T;
Iteration_no = Iteration_no + 1;
figure(1)
plot(x,T)
title_axis = sprintf('Iteration no = %d' , Iteration_no)
title(title_axis)
xlabel('dx space distribution')
ylabel('Temperature profile')
grid on
pause(0.1)
end
figure (2)
contourf(x, y, T)
clabel(contourf(x, y, T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('X axis')
ylabel('Y axis')
title(title_axis)
grid on
pause(0.03)
%%%% Gauss Siedel method
else
if (Thermal_Itr == 2)
while (Error > Tol)
for i = 2 : nx-1
for j = 2 : ny-1
T(i,j) = Term*(((T(i-1,j) + T_old(i+1,j) )/ dx^2) + ((T(i,j-1) + T_old(i,j+1)) / dy^2));
end
end
Error = (max(max(abs(T - T_old))));
T_old = T;
Iteration_no = Iteration_no + 1;
figure(3)
plot(x,T)
title_axis = sprintf('Iteration no = %d' , Iteration_no)
title(title_axis)
xlabel('dx space distribution')
ylabel('Temperature profile')
grid on
pause(0.1)
end
figure (4)
contourf(x, y, T)
clabel(contourf(x, y, T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('X axis')
ylabel('Y axis')
title(title_axis)
grid on
pause(0.03)
end
%%%% Successive over relaxation method
%else
if (Thermal_Itr == 3)
omega = 1.8;
while (Error > Tol)
for i = 2 : nx-1
for j = 2 : ny-1
T(i,j) = Term * (((T(i-1,j) + T_old(i+1,j)) / dx^2) + ((T(i,j-1) + T_old(i,j+1)) / dy^2));
T(i,j) = T_old(i,j)*(1-omega) + T(i,j)*omega;
end
end
Error = (max(max(abs(T - T_old))));
T_old = T;
Iteration_no = Iteration_no + 1;
figure(5)
plot(x,T)
title_axis = sprintf('Iteration no = %d' , Iteration_no)
title(title_axis)
xlabel('dx space distribution')
ylabel('Temperature profile')
grid on
pause(0.1)
end
figure (6)
contourf(x, y, T)
clabel(contourf(x, y, T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('X axis')
ylabel('Y axis')
title(title_axis)
grid on
pause(0.03)
end
end
%%% Transient Implicit solver
close all
clear all
clc
%Temperature value at centre with various solver for Rectangular BC
%Input conditions
L = 10;
nx = 11;
ny = nx;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
T = 300*ones(nx,ny);
% Time step definition
t = 100;
dt = 10;
nt = t/dt;
%Assignig periphery anotations
T_L = 400;
T_R = 600;
T_T = 800;
T_B = 900;
% Periphery BC
T(:,1) = 400;
T(:,nx) = 600;
T(1, :) = 800;
T(ny, :) = 900;
% Corner BC
T(1,1) = (T_L + T_T) / 2;
T(nx , 1) = (T_L + T_B) / 2;
T(nx, ny) = (T_R + T_B) / 2;
T(1, ny) = (T_T + T_R) / 2;
%%%% Jacobi method
Iteration_method = 3;
Itr = 1;
Itr_2 = 1;
T_old = T;
T_prev = T;
alfa = 1.75e-4;
K1 = alfa * dt / dx^2;
K2 = alfa * dt / dy^2;
term1 = 1 / (1 + 2*K1 + 2*K2);
term2 = K1 * term1;
term3 = K2 * term1;
Error = 9e9;
Tol = 1e-4;
if (Iteration_method == 1)
tic
for k = 1 : nt
while(Error > Tol)
for j = 2 : ny-1
for i = 2 : nx-1
H = T_old(i-1,j) + T_old(i+1,j);
V = T_old(i,j-1) + T_old(i,j+1);
T(i,j) = T_prev(i,j)*term1 + H*term2 + V*term3;
end
end
Error = max(max(abs(T - T_old)));
T_old = T;
figure(1)
plot(x,T)
xlabel('X axis')
ylabel('Temperature profile')
title_axis = sprintf('Iteration no = %d',Itr );
title(title_axis)
Itr = Itr + 1;
end
T_prev = T;
figure (2)
contourf(x, y, T)
clabel(contourf(x, y, T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('X axis')
ylabel('Y axis')
title_axis_2 = sprintf('Iteration no = %d',Itr_2 );
title(title_axis_2)
grid on
pause(0.03)
Itr_2 = Itr_2 + 1;
end
toc
end
%%%% Gauss siedel method
if (Iteration_method == 2)
tic
for k = 1 : nt
while(Error > Tol)
for j = 2 : ny-1
for i = 2 : nx-1
H = T(i-1,j) + T_old(i+1,j);
V = T(i,j-1) + T_old(i,j+1);
T(i,j) = T_prev(i,j)*term1 + H*term2 + V*term3;
end
end
Error = max(max(abs(T - T_old)));
T_old = T;
figure(3)
plot(x,T)
xlabel('X axis')
ylabel('Temperature profile')
title_axis = sprintf('Iteration no = %d',Itr );
title(title_axis)
Itr = Itr + 1;
end
T_prev = T;
figure (4)
contourf(x, y, T)
clabel(contourf(x, y, T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('X axis')
ylabel('Y axis')
title_axis_2 = sprintf('Iteration no = %d',Itr_2 );
title(title_axis_2)
grid on
pause(0.03)
Itr_2 = Itr_2 + 1;
end
toc
end
%%%% SOR method
omega = 1.1;
if (Iteration_method == 3)
tic
for k = 1 : nt
while(Error > Tol)
for j = 2 : ny-1
for i = 2 : nx-1
H = T(i-1,j) + T_old(i+1,j);
V = T(i,j-1) + T_old(i,j+1);
T(i,j) = T_prev(i,j)*term1 + H*term2 + V*term3;
T(i,j) = T_old(i,j)*(1-omega) + T(i,j)*omega;
end
end
Error = max(max(abs(T - T_old)));
T_old = T;
figure(5)
plot(x,T)
xlabel('X axis')
ylabel('Temperature profile')
title_axis = sprintf('Iteration no = %d',Itr );
title(title_axis)
Itr = Itr + 1;
end
T_prev = T;
figure (6)
contourf(x, y, T)
clabel(contourf(x, y, T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('X axis')
ylabel('Y axis')
title_axis_2 = sprintf('Iteration no = %d',Itr_2 );
title(title_axis_2)
grid on
pause(0.03)
Itr_2 = Itr_2 + 1;
end
toc
end
%%% Transient Explicit solver
close all
clear all
clc
%Temperature value at centre with various solver for Rectangular BC
%Input conditions
L = 10;
nx = 11;
ny = nx;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
T = 300*ones(nx,ny);
% Time step definition
t = 100;
dt = 10;
nt = t/dt;
%Assignig periphery anotations
T_L = 400;
T_R = 600;
T_T = 800;
T_B = 900;
% Periphery BC
T(:,1) = 400;
T(:,nx) = 600;
T(1, :) = 800;
T(ny, :) = 900;
% Corner BC
T(1,1) = (T_L + T_T) / 2;
T(nx , 1) = (T_L + T_B) / 2;
T(nx, ny) = (T_R + T_B) / 2;
T(1, ny) = (T_T + T_R) / 2;
%%%% SOR method
Iteration_method = 1;
Itr = 1;
Itr_2 = 1;
T_old = T;
T_prev = T;
alfa = 1.75e-4;
K1 = alfa * dt / dx^2;
K2 = alfa * dt / dy^2;
Error = 9e9;
Tol = 1e-4;
if (Iteration_method == 1)
tic
for k = 1 : nt
while(Error > Tol)
for j = 2 : ny-1
for i = 2 : nx-1
t1 = T_old(i,j);
t2 = K1 * (T_old(i+1,j) - 2*T_old(i,j) + T_old(i-1,j));
t3 = K2 * (T_old(i,j+1) - 2*T_old(i,j) + T_old(i,j-1));
T(i,j) = t1 + t2 + t3;
end
end
Error = max(max(abs(T - T_old)));
T_old = T;
figure(1)
plot(x,T)
xlabel('X axis')
ylabel('Temperature profile')
title_axis = sprintf('Iteration no = %d',Itr );
title(title_axis)
Itr = Itr + 1;
end
T_prev = T;
figure (2)
contourf(x, y, T)
clabel(contourf(x, y, T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('X axis')
ylabel('Y axis')
title_axis_2 = sprintf('Iteration no = %d',Itr_2 );
title(title_axis_2)
grid on
pause(0.03)
Itr_2 = Itr_2 + 1;
end
toc
end
