Solving the steady and unsteady 2-D Heat Conduction Problem
OBJECTIVE:
Solving the steady and unsteady 2-D heat conduction problem. In this we have to solve the conduction equation using a Transient solver and a Steady state solver using Iterative techniques (Jacobi.Gauss Seidal,SOR).
- Assume that the domain is a unit square (L = 1)
- Assume nx = ny [Number of points along the x direction is equal to the number of points along the y direction]
- Boundary conditions for steady and transient case
- Left:400K
- Top:600K
- Right:800K
- Bottom:900K
Background:
`(deltaT)/(deltat)+alphaDelta(T)=0`
`(T_(i,j)^(n+1)-T_(i,j)^n)/(Deltat)+alpha((delta^2T)/(deltax^2)+(delta^2T)/(deltay^2))^n=0`
`(T_(i,j)^(n+1)-T_(i,j)^n)/(Deltat)+alpha((T_(i+1,j)-2T_(i,j)+T_(i-1,j))/(Deltax^2)+(T_(i,j+1)-2T_(i,j)+T_(i,j-1))/(Deltay^2))^n=0`
----equation 1
By solving uncoupled linear equation 1 we get
`T_(i,j)^(n+1)=T_(i,j)^(n)+alphaDeltat((T_(i+1,j)-2T_(i,j)+T_(i-1,j))/(Deltax^2)+(T_(i,j+1)-2T_(i,j)+T_(i,j-1))/(Deltay^2))^n`
Implicit
`T_(i,j)^(n+1)=T_(i,j)^(n)+alphaDeltat((T_(i+1,j)-2T_(i,j)+T_(i-1,j))/(Deltax^2)+(T_(i,j+1)-2T_(i,j)+T_(i,j-1))/(Deltay^2))^(n+1)`
[A][X]=[B]
`k_1=alpha(Deltat)/(Deltax^2)`
`k_2=alpha(Deltat)/(Deltay^2)`
`T_(i,j)^(n+1)=T_(i,j)^(n)+k_1(T_(i+1,j)^(n+1)-2T_(i,j)^(n+1)+T_(i-1,j)^(n+1))+k_2(T_(i,j+1)^(n+1)-2T_(i,j)^(n+1)+T_(i,j-1)^(n+1))`
`T_(i,j)^(n+1)-k_1T_(i+1,j)^(n+1)+2k_1T_(i,j)^(n+1)-k_1T_(i-1,j)^(n+1))-k_2T_(i,j+1)^(n+1)+2k_2T_(i,j)^(n+1)-k_2T_(i,j-1)^(n+1)=T_(i,j)^(n)`
`T_(i,j)^(n+1)(1+2k_1+2K_2)-k_1T_(i+1,j)^(n+1)-k_1T_(i-1,j)^(n+1)-k_2T_(i,j+1)^(n+1)-k_2T_(i,j-1)^(n+1)=T_(i,j)^(n)`
[A][X]=[B]
Iterative solver
`T_(i,j)^(n+1)=T_(i,j)^(n)+k_1T_(i+1,j)^(n+1)+k_1T_(i-1,j)^(n+1)+k_2T_(i,j+1)^(n+1)+k_2T_(i,j-1)^(n+1)`
Initial Equation for "T" -Temparature
`T_(p)^(n+1)=T_(p)^(n)+k_1(T_(i+1,j)^(n+1)+T_(i-1,j)^(n+1))+k_2(T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1))`
Explicitly
`(delta^2T)/(deltax^2)+(delta^2T)/(deltay^2)=0`
`(2T_(i+1,j)-2T_(i,j)+T_(i-1,j))/(Deltax^2)+(2T_(i,j+1)-2T_(i,j)+T_(i,j-1))/(Deltay^2)=0`
Simplifying
`(T_(i+1,j)+T_(i-1,j))/(Deltax^2)+(T_(i,j+1)+T_(i,j-1))/(Deltay^2)=(2T_(i,j))/(Deltax^2)+(2T_(i,j))/(Deltay^2)`
`(T_(i+1,j)+T_(i-1,j))/(Deltax^2)+(T_(i,j+1)+T_(i,j-1))/(Deltay^2)=T_(i,j)(2((Deltax^2+Deltay^2)/(Deltax^2Deltay^2)))`
`(2((Deltax^2+Deltay^2)/(Deltax^2Deltay^2)))=k`
`T_(i,j)=1/k((T_(i+1,j)+T_(i-1,j))/(Deltax^2)+(T_(i,j+1)+T_(i,j-1))/(Deltay^2))`
`T_k=T_(i+1,j)`
`T_T=T_(i,j+1)`
`T_P=T_(i,j)`
`T_L=T_(i-1,j)`
`T_B=T_(i,j-1)`
Generalized Formula for Jacobi
`T_jac=D^-1(L+N)`
`L_jac=D^-1b`
`X_i^(k+1)=1/a_(i,j)(b_i-sum_(j=1)^(i-1)a_iX_j^k-sum_(i=i+1)^na_(i,j)X_j^k)`
D – Diagonal
L – Lower Diagonal
U – Upper Diagonal
steady state
`(delta^2T)/(deltax^2)+(delta^2T)/(deltay^2)=0`
`T_p=1/k((T_L+T_R)/(Deltax^2))+1/k((T_T+T_B)/(Deltay^2))`
unsteady state or transcient state
`(deltaT)/(deltat)=alpha((delta^2T)/(deltax^2)+(delta^2T)/(deltay^2))`
Discretizating
explicit
`(T_p^(n+1)-T_p^n)/(Deltat)=alpha((T_L-2T_P+T_R)/(Deltax^2)+(T_T-2T_P+T_B)/(Deltay^2))^n`
`T_p^(n+1)=T_p^n+k1(T_L-2T_P+T_R)^n+k_2(T_T-2T_P+T_B)^n`
`k_1=alpha(Deltat)/(Deltax^2)`
`k_2=alpha(Deltat)/(Deltay^2)`
implicit
`(T_p^(n+1)-T_p^n)/(Deltat)=alpha((T_L-2T_P+T_R)/(Deltax^2)+(T_T-2T_P+T_B)/(Deltay^2))^(n+1)`
`T_p^(n+1)=1/(1+2k_1+2k_2)(T_p^n+k_1(T_L+T_R)^(n+1)+k_2(T_T+T_B)^(n+1))`
`k_1=alpha(Deltat)/(Deltax^2)`
`k_2=alpha(Deltat)/(Deltay^2)`
term 1 `1/(1+2k_1+2k_2)`
term 2=`k_1*term1`
term 3=`k_2*term1`
Explicit Methods:
Jacobi:
clear all
close all
clc
%Length of the Domain
L = 1;
%number of gridpoints along x-axis
nx = 10;
ny = nx;
%Grid Spacing
x = linspace(0,L,nx)
y = linspace(0,L,ny)
[X,Y] = meshgrid(x,y);
dx = x(2)-x(1);
dy = y(2)-y(1);
%initialising velocity
T = ones(nx,ny)
%Boundary conditions
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
T(1,1) = (600+400)/2;
T(1,10) = (600+800)/2;
T(10,10) = (800+900)/2;
T(10,1) = (900+900)/2;
Told = T;
K = 2*((dx^2+dy^2)/((dx^2)*(dy^2)));
Tolerance = 1e-4;
error = 9e9;
%Jacobi_Solver
tic
jacobi_iteration_solver = 1;
while (error>Tolerance)
for i = 2:nx-1
for j = 2:ny-1
term_1 = (1/K)*((Told(i-1,j)+Told(i+1,j))/(dx^2))
term_2 = (1/K)*((Told(i,j-1)+Told(i,j+1))/(dy^2))
T(i,j) = term_1+term_2
end
end
error = max(max(abs(Told-T)));
Told = T;
jacobi_iteration_solver = jacobi_iteration_solver+1
end
Iteration_time = toc;
figure(1)
contourf(x,y,T);
[a,b] = contourf(x,y,T);
clabel(a,b);
colormap(jet);
colorbar
xlabel('x-axis')
ylabel('y-axis')
set(gca,'YDIR','reverse')
title_text = sprintf('Number of iterations =%d Iteration time =%d',jacobi_iteration_solver,Iteration_time);
title(title_text)
Guass-Seidel:
clear all
close all
clc
%Length of the Domain
L = 1;
%number of gridpoints along x-axis
nx = 10;
ny = nx;
%Grid Spacing
x = linspace(0,L,nx)
y = linspace(0,L,ny)
[X,Y] = meshgrid(x,y);
dx = x(2)-x(1);
dy = y(2)-y(1);
%initialising velocity
T = ones(nx,ny)
%Boundary conditions
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
T(1,1) = (600+400)/2;
T(1,10) = (600+800)/2;
T(10,10) = (800+900)/2;
T(10,1) = (900+900)/2;
Told = T;
K = 2*((dx^2+dy^2)/((dx^2)*(dy^2)));
Tolerance = 1e-4;
error = 9e9;
%Gauss_Seidel_Solver
tic
Gauss_Seidel_iteration_solver = 1;
while (error>Tolerance)
for i = 2:nx-1
for j = 2:ny-1
term_1 = (1/K)*((T(i-1,j)+Told(i+1,j))/(dx^2))
term_2 = (1/K)*((T(i,j-1)+Told(i,j+1))/(dy^2))
T(i,j) = term_1+term_2
end
end
error = max(max(abs(Told-T)));
Told = T;
Gauss_Seidel_iteration_solver = Gauss_Seidel_iteration_solver+1
end
Iteration_time = toc;
figure(1)
contourf(x,y,T);
[a,b] = contourf(x,y,T);
clabel(a,b);
colormap(jet);
colorbar
xlabel('x-axis')
ylabel('y-axis')
set(gca,'YDIR','reverse')
title_text = sprintf('Number of iterations =%d Iteration time =%d',Gauss_Seidel_iteration_solver,Iteration_time);
title(title_text)
Succesive Over Relaxation:
clear all
close all
clc
%Length of the Domain
L = 1;
%number of gridpoints along x-axis
nx = 10;
ny = nx;
%Grid Spacing
x = linspace(0,L,nx)
y = linspace(0,L,ny)
[X,Y] = meshgrid(x,y);
dx = x(2)-x(1);
dy = y(2)-y(1);
%initialising velocity
T = ones(nx,ny)
%Boundary conditions
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
T(1,1) = (600+400)/2;
T(1,10) = (600+800)/2;
T(10,10) = (800+900)/2;
T(10,1) = (900+900)/2;
Told = T;
K = 2*((dx^2+dy^2)/((dx^2)*(dy^2)));
Tolerance = 1e-4;
error = 9e9;
alpha = 1.2;
%Gauss_Seidel_Solver
tic
Succesive_Over_Relaxation_solver = 1;
while (error>Tolerance)
for i = 2:nx-1
for j = 2:ny-1
term_1 = (1/K)*((T(i-1,j)+Told(i+1,j))/(dx^2))
term_2 = (1/K)*((T(i,j-1)+Told(i,j+1))/(dy^2))
T(i,j) = (1-alpha)*Told(i,j)+(alpha*(term_1+term_2))
end
end
error = max(max(abs(Told-T)));
Told = T;
Succesive_Over_Relaxation_solver = Succesive_Over_Relaxation_solver+1
end
Iteration_time = toc;
figure(1)
contourf(x,y,T);
[a,b] = contourf(x,y,T);
clabel(a,b);
colormap(jet);
colorbar
xlabel('x-axis')
ylabel('y-axis')
set(gca,'YDIR','reverse')
title_text = sprintf('Number of iterations =%d Iteration time =%d',Succesive_Over_Relaxation_solver,Iteration_time);
title(title_text)
Implicit Methods:
Jacobi:
clear all
close all
clc
L = 1;
nx = 10;
ny = nx;
x = linspace(0,L,nx)
y = linspace(0,L,ny)
[X,Y] = meshgrid(x,y);
dx = x(2)-x(1);
dy = y(2)-y(1);
alpha = 1.4;
dt = 0.001;
T = ones(nx,ny)
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
T(1,1) = (600+400)/2;
T(1,10) = (600+800)/2;
T(10,10) = (800+900)/2;
T(10,1) = (900+900)/2;
K1 = alpha*(dt/(dx^2));
K2 = alpha*(dt/(dy^2));
jacobi_iteration_solver = 0;
Tolerance = 1e-4
Told = T;
T_prev_dt = Told;
tic
for nt = 1:1000;
error = 9e9;
while (error>Tolerance)
for i = 2:nx-1
for j = 2:ny-1
term_1 = (1+((2*K1)+(2*K2)))^-1;
term_2 = K1*term_1;
term_3 = K2*term_1;
H = (Told(i-1,j)+Told(i+1,j))
V = (Told(i,j-1)+Told(i,j+1))
T(i,j) = (T_prev_dt(i,j)*term_1)+(term_2*H)+(term_3*V);
end
end
error = max(max(abs(Told-T)));
Told = T;
jacobi_iteration_solver = jacobi_iteration_solver+1
end
T_prev_dt = T;
end
Iteration_time = toc;
figure(1)
contourf(x,y,T);
[a,b] = contourf(x,y,T);
clabel(a,b);
colormap(jet);
colorbar
xlabel('x-axis')
ylabel('y-axis')
set(gca,'YDIR','reverse')
title_text = sprintf('Number of iterations =%d Iteration time =%d',jacobi_iteration_solver,Iteration_time);
title(title_text)
Guass – Seidel:
clear all
close all
clc
L = 1;
nx = 10;
ny = nx;
x = linspace(0,L,nx)
y = linspace(0,L,ny)
[X,Y] = meshgrid(x,y);
dx = x(2)-x(1);
dy = y(2)-y(1);
alpha = 1.4;
dt = 0.001;
T = ones(nx,ny)
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
T(1,1) = (600+400)/2;
T(1,10) = (600+800)/2;
T(10,10) = (800+900)/2;
T(10,1) = (900+900)/2;
K1 = alpha*(dt/(dx^2));
K2 = alpha*(dt/(dy^2));
Guass_Seidel_solver = 0;
Tolerance = 1e-4
Told = T;
T_prev_dt = Told;
tic
for nt = 1:1000;
error = 9e9;
while (error>Tolerance)
for i = 2:nx-1
for j = 2:ny-1
term_1 = (1+((2*K1)+(2*K2)))^-1;
term_2 = K1*term_1;
term_3 = K2*term_1;
H = (T(i-1,j)+T(i+1,j))
V = (T(i,j-1)+T(i,j+1))
T(i,j) = (T_prev_dt(i,j)*term_1)+(term_2*H)+(term_3*V);
end
end
error = max(max(abs(Told-T)));
Told = T;
Guass_Seidel_solver = Guass_Seidel_solver+3
end
T_prev_dt = T;
end
Iteration_time = toc;
figure(1)
contourf(x,y,T);
[a,b] = contourf(x,y,T);
clabel(a,b);
colormap(jet);
colorbar
xlabel('x-axis')
ylabel('y-axis')
set(gca,'YDIR','reverse')
title_text = sprintf('Number of iterations =%d Iteration time =%d',Guass_Seidel_solver,Iteration_time);
title(title_text)
Succesive Over Relaxation:
clear all
close all
clc
L = 1;
nx = 10;
ny = nx;
x = linspace(0,L,nx)
y = linspace(0,L,ny)
[X,Y] = meshgrid(x,y);
dx = x(2)-x(1);
dy = y(2)-y(1);
alpha = 1.4;
w = 1.3;
dt = 0.01;
T = ones(nx,ny)
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
T(1,1) = (600+400)/2;
T(1,10) = (600+800)/2;
T(10,10) = (800+900)/2;
T(10,1) = (900+900)/2;
K1 = alpha*(dt/(dx^2));
K2 = alpha*(dt/(dy^2));
SOR_solver = 0;
Tolerance = 1e-4
Told = T;
T_prev_dt = Told;
tic
for nt = 1:1000;
error = 9e9;
while (error>Tolerance)
for i = 2:nx-1
for j = 2:ny-1
term_1 = (1+((2*K1)+(2*K2)))^-1;
term_2 = K1*term_1;
term_3 = K2*term_1;
H = (T(i-1,j)+T(i+1,j))
V = (T(i,j-1)+T(i,j+1))
G_S = (T_prev_dt(i,j)*term_1)+(term_2*H)+(term_3*V);
T(i,j) = ((1-w)*Told(i,j))+(w*(G_S))
end
end
error = max(max(abs(Told-T)));
Told = T;
SOR_solver = SOR_solver+1
end
T_prev_dt = T;
end
Iteration_time = toc;
figure(1)
contourf(x,y,T);
[a,b] = contourf(x,y,T);
clabel(a,b);
colormap(jet);
colorbar
xlabel('x-axis')
ylabel('y-axis')
set(gca,'YDIR','reverse')
title_text = sprintf('Number of iterations =%d Iteration time =%d',SOR_solver,Iteration_time);
title(title_text)
CONCLUSION:
- The steady-state solution is the fastest schemes to compute the solution. This is also because they are absent we the temperature terms in the equation which leads to faster computations.
- Explicit scheme forms to be the fastest among the transient method. It reuires way less number of iterations and computation time to converge compared to Implicit schemes.
- Relaxation factor(omega)must be maintained within a raange of 1.1-1.3. More or less tha the range it may become over relaxed or under-relaxed