Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Aim:Solving the steady and unsteady 2D heat conduction problem.
Transient 2D heat conduction (implicit)
Steady state
Transient 2D heat conduction (explicit)
Code:
clear all
close all
clc
number_node = 10;
domain_length = 1;
x = linspace(0,domain_length,number_node);
y = linspace(0,domain_length,number_node);
dx = domain_length/(number_node-1);
dy = domain_length/(number_node-1);
T = 300*ones(number_node,number_node);
T(1:end,1) = 400;
T(1:end,end) = 800;
T(end,1:end) = 900;
T(1,1:end) = 600;
T(1,1) = 500;
T(end,1) = 650;
T(1,end) = 700;
T(end,end) = 850;
T_old = T;
K = (2*(dx^2+dy^2))/(dx^2*dy^2);
omegha = 1.2;
error = 10;
tol = 1e-4;
no_iteration = 0;
method = input('Enter(1) for jacobi,(2) for gauss siediel,(3) for sos');
if method == 1
while(error>tol)
for i = 2:number_node-1
for j = 2:number_node-1
term1 = (T_old(i+1,j)+T_old(i-1,j))/dx^2;
term2 = (T_old(i,j+1)+T_old(i,j-1))/dy^2;
T(i,j) = (term1+term2)/K;
end
end
error = max(max(abs(T_old-T)));
T_old = T;
no_iteration = no_iteration+1;
end
elseif method == 2
while(error>tol)
for i = 2:number_node-1
for j = 2:number_node-1
term1 = (T_old(i+1,j)+T(i-1,j))/dx^2;
term2 = (T_old(i,j+1)+T(i,j-1))/dy^2;
T(i,j) = (term1+term2)/K;
end
end
error = max(max(abs(T_old-T)));
T_old = T;
no_iteration = no_iteration+1;
end
else
while(error>tol)
for i = 2:number_node-1
for j = 2:number_node-1
term1 = (T_old(i+1,j)+T(i-1,j))/dx^2;
term2 = (T_old(i,j+1)+T(i,j-1))/dy^2;
T(i,j) = ((1-omegha)*T_old(i,j))+(omegha*((term1+term2)/K));
end
end
error = max(max(abs(T_old-T)));
T_old = T;
no_iteration = no_iteration+1;
end
end
[X,Y] = meshgrid(x,y);
contourf(X,Y,T)
colorbar
set(gca,'ydir','reverse')
Unsteady:
clear all
close all
clc
% Domain and dicretization length
number_node = 10;
domain_length = 1;
%value of x and y at nodes in a discretization plane
x = linspace(0,domain_length,number_node);
y = linspace(0,domain_length,number_node);
dx = (domain_length/number_node-1);
dy = (domain_length/number_node-1);
dt = 1/10;
% Let define a temperature value at all nodes with initial guess 300k
T = 300*ones(number_node,number_node);
% Initia boundary conditions at left and right
T(1:end,1) = 400;
T(1:end,end) = 800;
% Initia boundary conditions at bottom and top
T(end,1:end) = 900;
T(1,1:end) = 600;
%average values of T at corner boundaries
T(1,1) = 500;
T(end,1) = 650;
T(1,end) = 700;
T(end,end) = 850;
%copy of Temperature to compute for previous space terms
T_old = T;
% copy of T for previous time_step
T_prev = T;
%constants
alpha = 1.5;
omegha = 1.15;
K1 = (alpha*dt)/(dx^2);
K2 = (alpha*dt)/(dy^2);
iteration = 0;
approach = input('Enter 1 for implicit or Enter 2 for explicit')
if approach == 1
method = input('Enter 1 for jacobi or 2 for gauss siedel or 3 for sos')
if method == 1
for nt = 1:1500
error = 10;
tol = 1e-4;
while(error>tol)
for i = 2:number_node-1
for j = 2:number_node-1
term1 = 1/(1+(2*K1)+(2*K2));
term2 = K1*(T_old(i+1,j)+T_old(i-1,j));
term3 = K2 *(T_old(i,j+1)+T_old(i,j-1));
T(i,j) = term1*(T_prev(i,j)+term2+term3);
end
end
iteration = iteration+1;
error = max(max((abs(T_old-T))));
T_old = T;
end
T_prev = T;
end
elseif method ==2
for nt = 1:1500
error = 10;
tol = 1e-4;
while(error>tol)
for i = 2:number_node-1
for j = 2:number_node-1
term1 = 1/(1+(2*K1)+(2*K2));
term2 = K1*(T_old(i+1,j)+T(i-1,j));
term3 = K2 *(T_old(i,j+1)+T(i,j-1));
T(i,j) = term1*(T_prev(i,j)+term2+term3);
end
end
iteration = iteration+1;
error = max(max(abs(T_old-T)));
T_old = T;
end
T_prev = T;
end
else
for nt = 1:1500
error = 10;
tol = 1e-4;
while(error>tol)
for i = 2:number_node-1
for j = 2:number_node-1
term1 = 1/(1+(2*K1)+(2*K2));
term2 = K1*(T_old(i+1,j)+T(i-1,j));
term3 = K2 *(T_old(i,j+1)+T(i,j-1));
term4 = term1*(T_prev(i,j)+term2+term3);
T(i,j)=(1-omegha)*T_old(i,j)+(omegha*term4);
end
end
iteration = iteration+1;
error = max(max(abs(T_old-T)));
T_old = T;
end
T_prev = T;
end
end
else
for nt = 1:1500
for i = 2:number_node-1
for j = 2:number_node-1
term1 = K1*(T_prev(i+1,j)-2*T_prev(i,j)+T_prev(i-1,j));
term2 = K2*(T_prev(i,j+1)-2*T_prev(i,j)+T_prev(i,j-1));
T(i,j) = T_prev(i,j)+term1+term2;
end
end
T_prev = T;
iteration = iteration+1;
end
end
[X,Y] = meshgrid(x,y);
figure(1)
contourf(X,Y,T)
colorbar
set(gca,'ydir','reverse')
Result:
contour
unsteady
explicit
Conclusion:
The temperature distribution plot obtained for both the steady and unsteady conditions is the same.Because in Transient analysis also, the solution approaches to steady-state.The top, bottom, left, and right side of the graphs representing their respective boundary conditions.