Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
The formula for steady state heat conduction
`(del^2T)/(delx^2)+(del^2T)/(dely^2)=0`
`(TL-2TP-TR)/(delx^2)+(T T-2TP-TB)/(dely^2)`
`(del^2T)/(delx^2)+(del^2T)/(dely^2)`=`(TL-2TP-TR)/(delx^2)+(T T-2TP-TB)/(dely^2)`
Transient 2D heat conduction :
`(delT)/(delt)=alpha[(del^2T)/(delx^2)+(del^2T)/(dely^2)]`
Time interval
n is for explicit
n+1 is for implicit
for explicit:
`Tp^(n+1) = Tp^n+k1[TL-2Tp+TR]^n+k2[T T-2Tp+TB]^n`
`k1=alpha(delt)/(delx^2)`
`k2=alpha(delt)/(dely^2)`
for implicit method:
`Tp^(n+1)=1/(1+2k1+2k2)[Tp^n+k1[TL+TR]^(n+1)+k2[T T+TB]^(n+1)]`
close all
clear all
clc
% inputs
L = 1;
nx = 10;
ny = 10;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = 1/(nx-1);
dy = 1/(ny-1);
tol = 1e-4;
error = 9e9;
omega = 1.2;
T = ones(nx,ny);
iterative_solver = 1;
% boundary conditions
T(:,1) = 400;
T(:,end) = 800;
T(1,:) = 900;
T(end,:) = 600;
Told = T;
k = 2/(dx^2)+2/(dy^2);
tic
if iterative_solver==1
jacobi_iter = 1;
while (error>tol)
for i = 2:nx-1
for j = 2:ny-1
V= (Told(i-1,j)+Told(i+1,j))/(k*(dx^2));
H = (Told(i,j-1)+Told(i,j+1))/(k*(dy^2));
T(i,j) = V+H;
end
end
error = max(max(abs(Told-T)));
Told = T;
jacobi_iter = jacobi_iter+1;
figure(1)
[X, Y] = contourf(x,y,T);
colorbar
colormap(jet)
clabel(X,Y);
text_iter = sprintf('jacobi_iter = %d',jacobi_iter)
title(text_iter)
ylabel('Y')
xlabel('X')
pause(0.03)
end
end
if iterative_solver == 2
sor_iter = 1;
while (error>tol)
for i = 2:nx-1
for j = 2:ny-1
V= (T(i-1,j)+T(i+1,j))/(k*(dx^2));
H = (T(i,j-1)+T(i,j+1))/(k*(dy^2));
T(i,j) = V+H;
end
end
error = max(max(abs(Told-T)));
Told = T;
sor_iter = sor_iter+1;
figure(1)
[X, Y] = contourf(x,y,T);
colorbar
colormap(jet)
clabel(X,Y);
text_iter = sprintf('gs_iter = %d',sor_iter)
title(text_iter)
ylabel('Y')
xlabel('X')
pause(0.03)
end
end
if iterative_solver == 3
sor_iter = 1;
while (error>tol)
for i = 2:nx-1
for j = 2:ny-1
V= (T(i-1,j)+Told(i+1,j))/(k*(dx^2));
H = (T(i,j-1)+Told(i,j+1))/(k*(dy^2));
T(i,j) = (1-omega)*Told(i,j)+omega*(H+V);
end
end
error = max(max(abs(Told-T)));
Told = T;
sor_iter = sor_iter+1;
figure(1)
[X, Y] = contourf(x,y,T);
colorbar
colormap(jet)
clabel(X,Y);
text_iter = sprintf('sor_iter = %d',sor_iter)
title(text_iter)
ylabel('Y')
xlabel('X')
pause(0.03)
end
end
toc
time_taken = toc
close all
clear all
clc
L = 1;
nx = 10;
ny = 10;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = L/(nx-1);
dy = L/(ny-1);
dt = 0.001;
t = 0.1;
nt = t/dt;
error = 9e9;
tol = 1e-4;
T = 298*ones(nx,ny);
% boundary conditions
T(:,1) = 400;
T(:,end) = 800;
T(1,:) = 900;
T(end,:) = 600;
Told = T;
alpa = 1.4;
k1 = alpa*(dt/dx^2);
k2 = alpa*(dt/dy^2);
CFL = k1+k2
Told = T;
explicit_iter = 1;
tic
for nt = 1:1400
while (error>tol)
for i = 2:nx-1
for j = 2:ny-1
H = Told(i-1,j)-2*Told(i,j)+Told(i+1,j);
V = Told(i,j-1)-2*Told(i,j)+Told(i,j+1);
T(i,j) = Told(i,j)+k1*(H)+k2*(V);
end
end
error = max(max(abs(Told-T)));
Told = T;
explicit_iter = explicit_iter+1
figure(1)
[A,B] = contourf(x,y,T);
colorbar
colormap("jet")
text_iter = sprintf('explicit_iter= %d',explicit_iter)
title(text_iter)
clabel(A,B);
pause(0.03)
end
end
time_taken = toc
close all
clear all
clc
L = 1;
nx = 11;
ny = 11;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = L/(nx-1);
dy = L/(ny-1);
dt = 0.001;
t = 0.1;
nt = t/dt;
error = 9e9;
tol = 1e-4;
T = 298*ones(nx,ny);
% boundary conditions
T(:,1) = 400;
T(:,end) = 800;
T(1,:) = 900;
T(end,:) = 600;
Told = T;
alpa = 1.4;
k1 = alpa*dt/(dx^2);
k2 = alpa*dt/(dy^2);
CFL = k1+k2
Told = T;
T_old_dt = T;
iteratve_solver = 2;
omega = 1.1;
t1 = (1+2*k1+2*k2)^(-1);
t2 = k1*t1;
t3 = k2*t1;
tic
for nt = 1:1400
if iteratve_solver == 1
jacobi_iter = 1;
while (error>tol)
for i = 2:nx-1
for j = 2:ny-1
H = Told(i-1,j)+Told(i+1,j);
V = Told(i,j-1)+Told(i,j+1);
T(i,j) = (T_old_dt(i,j)*t1)+(H*t2)+(V*t3);
end
end
error = max(max(abs(Told-T)));
Told = T;
T_old_dt = T;
jacobi_iter = jacobi_iter+1
figure(1)
[A,B] = contourf(x,y,T);
colorbar
colormap(jet)
text_iter = sprintf('implicit_Jacobi_iter= %d',jacobi_iter)
ylabel('Y')
title(text_iter)
clabel(A,B);
pause(0.003)
end
end
if iteratve_solver == 2
Gas_iter = 1
t1 = (1+2*k1+2*k2)^(-1);
t2 = k1*t1;
t3 = k2*t1;
while (error>tol)
for i = 2:nx-1
for j = 2:ny-1
H = T(i-1,j)+Told(i+1,j);
V = T(i,j-1)+Told(i,j+1);
T(i,j) = (T_old_dt(i,j)*t1)+(H*t2)+(V*t3);
end
end
error = max(max(abs(Told-T)));
Told = T;
T_old_dt = T;
Gas_iter = Gas_iter+1;
figure(1)
[A,B] = contourf(x,y,T);
colorbar
colormap(jet)
text_iter = sprintf('implicit_Gas_iter= %d',Gas_iter)
ylabel('Y')
title(text_iter)
clabel(A,B);
pause(0.03)
end
end
if iteratve_solver == 3
sor_iter = 1
t1 = (1+2*k1+2*k2)^(-1);
t2 = k1*t1;
t3 = k2*t1;
while (error>tol)
for i = 2:nx-1
for j = 2:ny-1
H = T(i-1,j)+Told(i+1,j);
V = T(i,j-1)+Told(i,j+1);
T(i,j) = (1-omega)*Told(i,j)+omega*((T_old_dt(i,j)*t1)+(H*t2)+(V*t3));
end
end
error = max(max(abs(Told-T)));
Told = T;
T_old_dt = T;
sor_iter = sor_iter+1;
figure(1)
[A,B] = contourf(x,y,T);
colorbar
colormap(jet)
text_iter = sprintf('implicit_sor_iter= %d',sor_iter)
ylabel('Y')
title(text_iter)
clabel(A,B);
pause(0.03)
end
end
end
time_taken = toc
Results;
In steady-state conditions, there is no time to be included so results would converge with fewer iterations.
Unsteady state conditions, take much time to compute when compared to steady-state
SOR gives answers with fewer iterations and takes no longer to compute