Solving the steady and unsteady 2D heat conduction problem
Steady State Solution:
We make the assumptions as there is no convection, no internal heat generation and no change of temperature gradient with respect to time, therefore the generalised energy equation in 2D becomes as:
`(del^2T)/(delx^2)+(del^2T)/(dely^2) = 0`
Solving the steady state equations does not involve time marching therefore the temperature gradient do not change with respect to time.
We discritize spatial term by central differencing:
`(del^2T)/(delx^2) = (T_(i-1,j)-2T_(i,j)+T_(i+1,j))/(Deltax^2)`
`(del^2T)/(dely^2) = (T_(i,j-1)-2T_(i,j)+T_(i,j+1))/(Deltay^2)`
therefore we get
`(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) = 0`
we take `Deltax = Deltay`
`T_i = 1/4(T_(i-1,j)+T_(i+1,j)+T_(i,j-1)+T_(i,j+1))`
Following is the Script and Functions to perform steady state analysis:
%Steady State Solution
clear all
close all
clc
%Step 1:
%First lets Initial the given condition
L = 1;
nx = 20;
ny = 20;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
error = 9e9;
tol = 1e-4;
% Step 2:
% Define boundary condition
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = ones(nx,ny);
% Temperature at Boundary
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
% Temperature at the corner by average value temperature
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;
k1 = (dx^2)/(2*(dx^2+dy^2));
k2 = (dy^2)/(2*(dx^2+dy^2));
i_s=1; %Selection of method
if i_s==1;
jacobi_iter(T,x,y,nx,ny,tol,error,k1,k2);
elseif i_s==2;
gauss_iter(T,x,y,nx,ny,tol,error,k1,k2);
elseif i_s==3;
sor_iter(T,x,y,nx,ny,tol,error,k1,k2);
end
%Jacobi Iteration
function jacobi_iter(T,x,y,nx,ny,tol,error,k1,k2)
[X, Y] = meshgrid(x,y);
Told = T;
tic
jacobi_iterr = 1;
while(error>tol)
for i = 2:nx-1
for j = 2:ny-1
T(i,j)=k1*(Told(i-1,j)+Told(i+1,j))+k2*(Told(i,j-1)+Told(i,j+1))
end
end
error = max(max(abs(Told-T)));
Told = T;
jacobi_iterr = jacobi_iterr+1;
end
time_taken=toc;
figure(1)
contourf(x,y,T);
[a b] = contourf(x,y,T);
clabel(a,b);
xlabel(\'X\')
ylabel(\'Y\')
colorbar
title_text = sprintf(\'Jacobi No of iteration = %d\' , jacobi_iterr);
title(title_text)
file_text = sprintf(\'Jacobi.png\');
saveas(gca,file_text);
fprintf(\'Number of iterations done are %d \\n\',jacobi_iterr);
fprintf(\'Time taken for solving equation is %f \\n\',time_taken);
end
%Gauss Iteration
function gauss_iter(T,x,y,nx,ny,tol,error,k1,k2)
[X, Y] = meshgrid(x,y);
Told = T;
gs_iter = 1;
while(error>tol)
for i = 2:nx-1
for j = 2:ny-1
T(i,j) = k1*(T(i,j-1)+Told(i,j+1))+k2*(T(i-1,j)+Told(i+1,j));
end
end
error =max(max(abs(Told-T)));
Told = T;
gs_iter = gs_iter+1;
end
time_taken=toc;
figure(1)
contourf(x,y,T);
[a b] = contourf(x,y,T);
clabel(a,b);
colorbar
title_text = sprintf(\'Gauss no of iteration = %d\' , gs_iter);
title(title_text)
file_text = sprintf(\'G_S.png\');
saveas(gca,file_text);
fprintf(\'Number of iterations done are %d \\n\',gs_iter);
fprintf(\'Time taken for solving equation is %f \\n\',time_taken);
end
%SOR Iteration
function sor_iter(T,x,y,nx,ny,tol,error,k1,k2)
% relaxation factor calculation
z=cos(pi/nx)+cos(pi/ny);
omega = 4/(2+sqrt(4-z^2));
Told = T;
[X Y] = meshgrid(x,y);
sor_iter = 1;
while(error>tol)
for i = 2:nx-1
for j = 2:ny-1
T(i,j) = omega*(k1*(T(i,j-1)+T(i,j+1))+k2*(T(i-1,j)+T(i+1,j)))+(1-omega)*Told(i,j);
end
end
error = max(max(abs(Told-T)));
Told = T;
sor_iter = sor_iter+1;
end
time_taken=toc;
figure(1)
contourf(x,y,T);
[a b] = contourf(x,y,T);
clabel(a,b);
colorbar;
title_text = sprintf(\'Sor no of iteration = %d\' , sor_iter);
title(title_text)
file_text = sprintf(\'S-SOR.png\');
saveas(gca,file_text);
fprintf(\'Relaxation factor considered is %d \\n\',omega);
fprintf(\'Number of iterations done are %d \\n\',sor_iter);
fprintf(\'Time taken for solving equation is %f \\n\',time_taken);
end
Output:
1.Jacobi Method:
Number of iterations done are 869
Time taken for solving equation is 361.586736
2.Gauss Seidel Method:
Number of iterations done are 463
Time taken for solving equation is 0.005830
3.Successive Over Relaxation Method:
Relaxation factor considered is 1.729454e+00
Number of iterations done are 61
Time taken for solving equation is 0.002747
Discussion (Steady State):
From the above results we can deduce that amongst the three iterative methods we get the results faster and with less error in the case of Successive Over Relaxation method than the Gauss Seidel method which is faster than the Jacobi method .
Unsteady State Solution:
The Generalised energy equation in 2D becomes as :
`(del^2T)/(delx^2)+(del^2T)/(dely^2) = alpha(delT)/(delt)`
the above equation is obtained by making assumptions that there is no convection , no internal heat generation or dessipation and aplha being constant .
We discritize spatial term by central differencing:
`(del^2T)/(delx^2) = (T_(i-1,j)-2T_(i,j)+T_(i+1,j))/(Deltax^2)`
`(del^2T)/(dely^2) = (T_(i,j-1)-2T_(i,j)+T_(i,j+1))/(Deltay^2)`
we discritize temporal terms by Forward differencing:
`(delT)/(delt) = ((T^(n+1))_(i,j)-(T^(n))_(i,j))/(Deltat)`
hence ,
`(T_(i-1,j)-2T_(i,j)+T_(i+1,j))^(n+1)/(Deltax^2) + (T_(i,j-1)-2T_(i,j)+T_(i,j+1))^(n+1)/(Deltay^2) = ((T^(n+1))_(i,j)-(T^(n))_(i,j))/(Deltat)`
we have defined 2 constants namely k4 and k5
`K4 = alpha(Deltat)/((Deltax)^2)`
`K5 = alpha(Deltat)/((Deltay)^2)`
as `Deltax = Deltay`
hence:`(T^(n+1))_(i,j)=(T^(n))_(i,j)+K4((T^(n+1))_(i-1,j)-2(T^(n+1))_(i,j)+(T^(n+1))_(i+1,j))+K5((T^(n+1))_(i,j-1)-2(T^(n+1))_(i,j)+(T^(n+1))_(i,j+1))`
therefore:
`(T^(n+1))_(i,j)=((T^(n))_(i,j)+(K4(T_(i-1,j)+T_(i+1,j))+K5(T_(i,j-1)+T_(i,j+1)))^(n+1))/(1+2K4+2K5)`
In Transient condition Temperature changes with time and reaches a steady state. The above equation is iterated both explicitly and implicitly.
EXPLICIT METHOD:
In explicit method the temperature at next level is calculated by means of surrounding temperature values and increasing the time steps i.e dt shall result in oscillation of result and the solution diverges from the original result.
Following is the code for
1. 2 D Transient Unsteady equation solved explicitly:
% 2D Transient State Heat Conduction Equation explicitly:
clear all
close all
clc
%Step 1:
%First lets Initial the given condition
L=1;
nx= 20;
ny = 20;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy= y(2)-y(1);
dt = 30;
nt = 1000;
tol = 1e-4;
error = 9e9;
% Step 2:
% Define boundary condition
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T =ones(nx,ny);
% Temperature at Boundary
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
% Temperature at corner (average temperature)
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;
alpha=2.3e-5;
CFL = alpha*dt/(dx^2);
k4 = (alpha*dt)/(dx^2);
k5 = (alpha*dt)/(dy^2);
[X,Y] = meshgrid(x,y);
Told = T;
T_prev_dt = Told;
tic;
counter = 0;
for k = 1:nt
% Jacobi explicit Iteration Scheme
for i = 2:nx-1
for j = 2:ny-1
term1 = Told(i,j);
term2 = k4*(Told(i-1,j)-2*Told(i,j)+Told(i+1,j));
term3 = k5*(Told(i,j-1)-2*Told(i,j)+Told(i,j+1));
T(i,j) = term1+term2+term3;
end
end
Told = T;
counter = counter+1;
end
time_required = toc;
figure(1)
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
title(sprintf(\'2D Unsteady state heat conduction-- Explicit Method n Number of iterations: %d\',counter));
xlabel(\'X\');
ylabel(\'Y\');
hold on
saveas(gca,\'2D Unsteady state heat conduction -Explicit Method .png\');
fprintf(\'Number of iterations done are %d \\n\',counter);
fprintf(\'Time taken for solving equation is %f \\n\',time_required);
fprintf(\'cfl number %f\',CFL)
::Graph 2D Unsteady Transient Heat conduction by Explicit Method::
Number of iterations done are 1000
Time taken for solving equation is 0.017363
cfl number 0.249090
IMPLICIT METHOD:
2. 2 D Transient Unsteady equation solved Implicitly:
% TwoD_Unsteady_Transient_state_Heat_Conduction_Equation_implicitly:
clear all
close all
clc
%Step 1:
%First lets Initial the given condition
L=1;
nx= 20;
ny =20;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy= y(2)-y(1);
dt = 30;
nt = 1000;
tol = 1e-4;
error = 9e9;
% Step 2:
% Define boundary condition
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T =ones(nx,ny);
% Temperature at Boundary
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
% Temperature at corner (average temperature)
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;
alpha=2.3e-5;
k4=alpha*dt/(dx^2);
k5=alpha*dt/(dy^2);
i_s = 1;%Selection of method
if i_s==1;
jacobi_unsteady_implicit(T,x,y,nx,ny,tol,k4,k5,nt);
elseif i_s==2;
gauss_unsteady_implicit(T,x,y,nx,ny,tol,k4,k5,nt);
elseif i_s==3;
sor_unsteady_implicit(T,x,y,nx,ny,tol,k4,k5,nt);
end
%Jacobi Iteration Implicit
function jacobi_unsteady_implicit(T,x,y,nx,ny,tol,k4,k5,nt)
Told = T;
T_prev_dt = Told;
term1=(1+2*k4+2*k5)^(-1);
term2=k4*term1;
term3=k5*term1;
tic
counter = 1;
for k = 1:nt %Time loop
error = 9e9;
while (error>tol) % Convergence loop
% Nodal loop
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_prev_dt(i,j)*term1)+(H*term2)+(V*term3));
end
end
error=max(max(abs(Told-T)));
Told = T;
counter = counter+1;
end
T_prev_dt = T;
end
time_required=toc;
figure(1)
[X, Y] = meshgrid(x,y);
contourf(x,y,T);
[a b] = contourf(x,y,T);
clabel(a,b);
title(sprintf(\'Plot for solving Transient state heat conduction equation implicitly using Jacobi method \\n No. of Iterations: %d\\n\',counter));
xlabel(\'X\');
ylabel(\'Y\');
colorbar
hold on;
saveas(gca,\'2D Transient state heat conduction Jacobi method IMPLICITLY1.png\');
fprintf(\'Number of iterations done are %d \\n\',counter);
fprintf(\'Time taken for solving equation is %f \\n\',time_required);
end
%Guass Seidel
function gauss_unsteady_implicit(T,x,y,nx,ny,tol,k4,k5,nt)
term1=(1+(2*k4)+(2*k5))^-1;
term2=k4*term1;
term3=k5*term1;
Told=T;
T_prev_dt = Told;
gs_iter = 1;
tic
for k =1:nt %Time loop
error = 9e9;
while (error>tol) % Convergence loop
% Nodal loop
for i = 2:nx-1
for j = 2:ny-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)*term1)+(H*term2)+(V*term3));
end
end
error=max(max(abs(Told-T)));
Told = T;
gs_iter = gs_iter+1;
end
T_prev_dt = Told;
end
time_required=toc;
figure(1)
[X, Y] = meshgrid(x,y);
figure(1)
contourf(x,y,T);
[a b] = contourf(x,y,T);
clabel(a,b);
title(sprintf(\'Plot for solving Transient state heat conduction equation implicitly By using gauss method. \\n No. of Iterations: %d\\n\',gs_iter));
xlabel(\'X\');
ylabel(\'Y\');
colorbar
hold on;
saveas(gca,\'2D Transient state heat conduction Gauss sidel method IMPLICITLY.png\');
fprintf(\'Number of iterations done are %d \\n\',gs_iter);
fprintf(\'Time taken for solving equation is %f \\n\',time_required);
end
%Successive Over Relaxation
function sor_unsteady_implicit(T,x,y,nx,ny,tol,k4,k5,nt)
Told = T;
T_prev_dt = Told;
z=cos(pi/nx)+cos(pi/ny);
omega = 4/(2+sqrt(4-z^2));
term1=(1+(2*k4)+(2*k5))^-1;
term2=k4*term1;
term3=k5*term1;
sor_iter = 1;
tic
for k = 1:nt %Time loop
error = 9e9;
while (error>tol) % Convergence loop
% Nodal loop
for i = 2:nx-1
for j = 2:ny-1
H=(T(i-1,j)+T(i+1,j));
V=(T(i,j-1)+T(i,j+1));
T(i,j) =omega*((T_prev_dt(i,j)*term1)+(H*term2)+(V*term3))-(omega-1)*T_prev_dt(i,j);
end
end
error=max(max(abs(Told-T)));
Told = T;
sor_iter = sor_iter+1;
end
T_prev_dt = Told;
end
time_required=toc;
figure(1)
contourf(x,y,T);
[a b] = contourf(x,y,T);
clabel(a,b);
title(sprintf(\'Plot for solving Transient state heat conduction implicitly using SOR method. \\n No. of Iterations: %d\\n\',sor_iter));
xlabel(\'X\');
ylabel(\'Y\');
colorbar
hold on;
saveas(gca,\'2D Transient state heat conduction SOR method IMPLICITLY1.png\');
fprintf(\'Relaxation Factor %f\\n\',omega);
fprintf(\'Number of iterations done are %d \\n\',sor_iter);
fprintf(\'Time taken for solving equation is %f \\n\',time_required);
endGraphs:
1.Jacobi Unsteady Implicit
Number of iterations done are 8045
Time taken for solving equation is 0.079604
2.Gauss Seidel Unsteady Implicit
Number of iterations done are 5862
Time taken for solving equation is 0.058604
3.SOR Unsteady Implicit
Relaxation Factor 1.729454
Number of iterations done are 4222
Time taken for solving equation is 0.047837
Discussion(Unsteady state):
Implicit Methods are found to be unconditionally stable . Successive Over relaxation method takes less time and iterations amongst the methods.