Solving 2D Heat Conduction Equation for Stedy and Unsteady State using MATLAB
The 3D Heat Conduction is givien by:
`(delT)/(delt)+alpha((del^2T)/(delx^2)+(del^2T)/(dely^2)+(del^2T)/(delz^2))=0`
Similarly the 2D Heat Conduction equation is given by:
`(delT)/(delt)+alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))=0`
In this Project the 2D Heat Conduction Equation was solved for :
- Steady Case
- Unsteady/ Transient case
- The Steady Heat conduction is the special case of heat conduction in which the change of temperature w.r.t. the time is 0. The 2D Heat Conduction equation can be converted to steady form as:`(del^2T)/(delx^2)+(del^2T)/(dely^2)=0`
- The Trainsient heat cinsuction equation is the equation in which there is a temperature gradient which means that the temperature changes w.r.t. time.`(delT)/(delt)+(del^2T)/(delx^2)+(del^2T)/(dely^2)=0`
The Steady/Unsteady heat conduction equation an be solved by different numerical methods.
The Methods, solvers used for solving Steady cases were
- Jacobi Method
- Gauss Siedel Method
- Succsessive over relaxation Method
Similalry, For the transient cases the method used were:
- Implicit Jacobi Method
- Implicit G.S. Method
- Implicit SOR Method
- Explicit Method
The inputs given for the solution are:
1. The domain size was kept as 1x1
2. The domain was discretised with size of 0.02
3. The Boundary condition was set to:
Left : 400K
Top: 600K
Right: 800K
Bottom:900K
Ambient Temp: 300K
4. For the transient cases:
total time : 5 sec
dt : 0.1
The Contour plot for initial condition of temprature distribution is shown below:
The contour shows the defined wall temperature and the initial ambinet temp.
The code for solving the various cases using various methods is given below:
clear all
close all
clc
% Inputs
nx = 10;
ny = 10;
x_domain_length = 1 ;
y_domain_length = 1 ;
omega_sor = 1.7
% Select method : 1. Steady State 2. Unsteady State
state = 2
% Select method : 1. Implicit Jacobi 2. Implicit Gauss-Seidel 3. Implicit
% SOR 4. Explicit Method(only for state 2)
method = 1
% Defining Grid Points
x_domain = linspace(0,x_domain_length,nx);
y_domain = linspace(0,y_domain_length,ny);
dx = x_domain(3)-x_domain(2);
dy = y_domain(3)-y_domain(2);
% Time Domain and Time Steps
total_time = 5;
dt = 0.1;
dt_exp = 0.001;
time_steps = total_time/dt;
time_steps_exp = total_time/dt_exp;
% Constants
ambient_temp = 300;
k1 = (dt/dx^2);
k2 = (dt/dy^2);
term1 = 1/(1+2*k1+2*k2);
% Intialization
for i= 1:length(x_domain)
for j= 1:length(y_domain)
temp(i,j) = ambient_temp;
end
end
% Bottom Wall
temp (1, :) = 900;
% Top Wall
temp (end, :) = 600;
% Left Wall
temp ( :, 1) = 400;
% Right Wall
temp ( :, end) = 800;
% Plot of initial conditions
temp_old = temp;
figure (1)
% Contour Plot for Initial Condition
contourf( x_domain, y_domain, temp)
jac_error = 1e+9;
gs_error= 1e+9;
sor_error= 1e+9;
iteration = 1;
temp_init = temp;
% Steady State Cases
if state ==1
% Constants for Steady State Analysis
c1 = 1/dx^2;
c2 = 1/dy^2;
c3 = c1/(2*(c1+c2));
c4 = c2/(2*(c1+c2));
% 1.1 Steady State analysis (Jacobi Method)
if method == 1
jac_error= 1e+9;
iteration = 1;
temp_t_minus_1 = temp;
while (jac_error > 1e-4)
temp_old = temp;
% Y axis Loop
for j = 2:ny-1
% X axis Loop
for i = 2: nx-1
c1 = 1/dx^2;
c2 = 1/dy^2;
c3 = c1/(2*(c1+c2));
c4 = c2/(2*(c1+c2));
H = temp_old(i-1,j)+temp_old(i+1,j);
V = temp_old(i,j+1)+temp_old(i,j-1);
temp(i,j) = c3*H+c4*V;
%temp(i,j) = ((temp_t_minus_1(i,j))+(k1*temp_old(i-1,j))+(k1*temp_old(i+1,j))+(k2*temp_old(i,j+1))+(k2*temp_old(i,j-1)))/term1;
end
end
iteration = iteration +1;
%Contour Plots
figure (2)
contourf(x_domain, y_domain,temp)
jac_error = max(max(abs(temp-temp_old)));
title(sprintf( \'Steady State analysis (Jacobi Method), Iteration no.= %d\', iteration))
xlabel(\'X Axis\')
ylabel(\'Y Axis\')
colorbar
end
end
% 1.2 Steady State Analysis (Gauss Seidel Method)
if method == 2
gs_error= 1e+9
iteration = 1;
temp_t_minus_1 = temp;
while (gs_error > 1e-4)
temp_old = temp;
% Y axis Loop
for j = 2:ny-1
% X axis Loop
for i = 2: nx-1
c1 = 1/dx^2;
c2 = 1/dy^2;
c3 = c1/(2*(c1+c2));
c4 = c2/(2*(c1+c2));
H = temp(i-1,j)+temp(i+1,j);
V = temp(i,j+1)+temp(i,j-1);
temp(i,j) = c3*H+c4*V;
end
end
iteration = iteration +1;
% Contour Plots
figure (2)
contourf(x_domain, y_domain,temp)
gs_error = max(max(abs(temp-temp_old)));
title(sprintf(\'Steady State Analysis (Gauss Seidel Method), Iteration no.= %d\', iteration));
xlabel(\'X Axis\')
ylabel(\'Y Axis\')
colorbar
end
end
% Steady State Analysis (SOR Method)
if method == 3
sor_error= 1e+9
iteration = 1;
temp_t_minus_1 = temp;
while (sor_error > 1e-4)
temp_old = temp;
% Y axis Loop
for j = 2:ny-1
% X axis Loop
for i = 2: nx-1
c1 = 1/dx^2;
c2 = 1/dy^2;
c3 = c1/(2*(c1+c2));
c4 = c2/(2*(c1+c2));
H = temp(i-1,j)+temp(i+1,j);
V = temp(i,j+1)+temp(i,j-1);
temp(i,j) = omega_sor*(c3*H+c4*V)+(1-omega_sor)*temp_old(i,j);
end
end
iteration = iteration +1;
%Contour Plots
figure (2)
contourf(x_domain, y_domain,temp)
sor_error = max(max(abs(temp-temp_old)));
title(sprintf(\' Steady State Analysis (SOR Method), Iteration no.= %d\', iteration))
xlabel(\'X Axis\')
ylabel(\'Y Axis\')
colorbar
end
end
end
if state ==2
% Transient analysis (Jacobi Method)
if method == 1
tic
% Time Loop
for t = 1:time_steps
jac_error= 1e+9
iteration = 1;
temp_t_minus_1 = temp;
while (jac_error > 1e-4)
temp_old = temp;
% Y axis Loop
for j = 2:ny-1
% X axis Loop
for i = 2: nx-1
term1 = 1/(1+2*k1+2*k2);
term2 = k1*term1;
term3 = k2*term1;
H = temp_old(i-1,j)+temp_old(i+1,j);
V = temp_old(i,j+1)+temp_old(i,j-1);
temp(i,j) = term1*temp_t_minus_1(i,j)+ term2*H+ term3*V;
end
end
iteration = iteration +1;
% Contour Plots
figure (2)
contourf(x_domain, y_domain,temp)
jac_error = max(max(abs(temp-temp_old)));
title(sprintf(\'Transient analysis (Jacobi Method), time step = %d, iteration no.= %d\',t, iteration))
xlabel(\'X Axis\')
ylabel(\'Y Axis\')
colorbar
end
end
toc
end
% Transient analysis (Gauss Seidel Method)
if method == 2
tic
% Time Loop
for t = 1:time_steps
gs_error= 1e+9
iteration = 1;
temp_t_minus_1 = temp;
while (gs_error > 1e-4)
temp_old = temp;
% Y axis Loop
for j = 2:ny-1
% X axis Loop
for i = 2: nx-1
term1 = 1/(1+2*k1+2*k2);
term2 = k1*term1;
term3 = k2*term1;
H = temp(i-1,j)+temp(i+1,j);
V = temp(i,j+1)+temp(i,j-1);
temp(i,j) = term1*temp_t_minus_1(i,j)+ term2*H+ term3*V;
%temp(i,j) = ((temp_t_minus_1(i,j))+(k1*temp(i-1,j))+(k1*temp(i+1,j))+(k2*temp(i,j+1))+(k2*temp(i,j-1)))/term1;
end
end
iteration = iteration +1;
% Contour Plots
figure (2)
contourf(x_domain, y_domain,temp)
gs_error = max(max(abs(temp-temp_old)));
title(sprintf(\'Transient analysis (Gauss Seidel Method), time step = %d, iteration no.= %d\',t, iteration))
xlabel(\'X Axis\')
ylabel(\'Y Axis\')
colorbar
end
end
toc
end
% Transient analysis (SOR Method)
if method == 3
tic
temp = temp_init
% Time Loop
for t = 1:time_steps
sor_error= 1e+9
iteration = 1;
temp_t_minus_1 = temp;
while (sor_error > 1e-4)
temp_old = temp;
sor_old_error = sor_error;
% Y axis Loop
for j = 2:ny-1
% X axis Loop
for i = 2: nx-1
term1 = 1/(1+2*k1+2*k2);
term2 = k1*term1;
term3 = k2*term1;
term4 = 1+omega_sor*(term1-1);
H = temp(i-1,j)+temp(i+1,j);
V = temp(i,j+1)+temp(i,j-1);
temp(i,j) = temp(i,j) + omega_sor*(term1*temp_t_minus_1(i,j)+ term2*H+ term3*V - temp(i,j));
end
end
iteration = iteration +1;
% Contour Plots
figure (2)
contourf(x_domain, y_domain,temp)
colorbar
sor_error = max(max(abs(temp-temp_old)));
title(sprintf(\'Transient analysis (SOR Method), time step = %d, iteration no.= %d\',t, iteration))
end
end
sor_matrix = temp(i,j);
toc
end
% Transient analysis (Explicit Method)
if method == 4
tic
for t = 1:time_steps_exp
exp_error = 1e+9
iteration = 1
temp_t_minus_1 = temp;
while (exp_error > 1e-4)
temp_old = temp;
% Y axis Loop
for j = 2:ny-1
% X axis Loop
for i = 2: nx-1
term1 = 1/(1+2*k1+2*k2);
term2 = k1*term1;
term3 = k2*term1;
H = temp_t_minus_1(i-1,j)+temp_t_minus_1(i+1,j);
V = temp_t_minus_1(i,j+1)+temp_t_minus_1(i,j-1);
temp(i,j) = term1*temp_t_minus_1(i,j)+term2*H+term2*V;
end
end
iteration = iteration +1;
% Contour Plots
figure (2)
contourf(x_domain, y_domain,temp)
exp_error = max(max(abs(temp-temp_old)))
title(sprintf(\'time step = %d, iteration no.= %d\',t, iteration))
end
end
toc
end
end1. Steady Cases:
The contour plots for temperature distribution was obtained for the steady states as shown in the figures below.
- In Steady State the SOR, Method have the least Iterations at selected alpha value.
- The highest iterations was required by Jacobi method.
- After decreasing the dx and dy value the number of iterations increased for each method.
- So the computational time increases with mesh refinement for all the methods.
- The results obtained with all the method are similar with a little deviation.
- The obtained contour shows the variation Temprature in the mesh for steady condition.
As we can see from the number of iteration shown on the plots, the highest number of iteration is for the jacobi method to get the solution with 1e-4 error and the least is for the SOR method which shows that the SOR method is very much faster than the jacobi and Gauss siedel. This is due to the relaxation factor involved in the solution.
2. Transient Cases:
In transient case the equations was solved explicitly as well as implicitly. For solving the equations implicitly the jacobi method, GS method and SOR methods were used.
Explicit method:
As the Explicit method is conditonally stable the value of dt was obtained based on CFL number. The graph obtained from the explicit method is shown in the figure below.
Implicit Method:
The implicit methods are unconditianally stable so the same dt value was kept as 0.01.
The Contour Plots showing temprature distribution was obtained for the transient states solved with explicit method as well as implicit method as shown in the figures below.
After comparing the Contours obtained from the all four methods and the values obtained from the matrix all the methods give very similar data and results.
Comparing the Explicit and Implicit method the Explicit methods are much faster but needs lower time steps according to the CFL Number.
Compairing the three Implicit methods the SOR methods gives the results much faster when keeping the omega more than 1 and gives the results at slower rate when the omega value is less than one. In this case the value of omega is 1.7 so it converges very much faster than other two methods.
The time Comparison for Jacobi, GS and SOR with nx=10 and ny=10 is shown below:
So Convergence Speed (SOR>GS>Jaccobi).