Week 5 - Mid term project - Solving the steady and unsteady 2D heat conduction problem

We have to solve for two dimensional heat conduction equation in this challenge.

First of all the two dimensional heat conduction equation is,

`(del^2T)/(delx^2)+(del^2T)/(dely^2)=0`

When we discritize the above equation by using central differencing scheme,we get the following

`(del^2T)/(delx^2)=(T_(i-1,j)-2*T_(i,j)+T_(i+1,j))/(Deltax^2)`

`(del^2T)/(dely^2)=(T_(i,j-1)-2*T_(i,j)+T_(i,j+1))/(Deltay^2)`

Then the heat conduction becomes,

`[(T_(i-1,j)-2*T_(i,j)+T_(i+1,j))/(Deltax^2)]+[(T_(i,j-1)-2*T_(i,j)+T_(i,j+1))/(Deltay^2)]=0`

Given Data:

• we want to assume the domain length as unit square.[L=1m].
• nx=ny.
• Absolute error=1e-4
• Boundary conditions,

                           Top Boundary = 600 K

                           Bottom Boundary = 900 K

                           Left Boundary = 400 K

                           Right Boundary = 800 K

Steady State:

                              The main program to solve two dimensional heat conduction equation is given below,

%Solving Steady state 2-dimensional equation
clear all
close all
clc


%JACOBI method
%inputs for the problem

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;

%Boundary conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;

T=ones(nx,ny);

% calculating the 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;


% Calculating average temperatures at the matrix corners.
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));


[X,Y] = meshgrid(x,y);
Told = T; 

% Jacobi Iteration.
for iterative_solver = 1
tic;
 if iterative_solver == 1
    jacobi_iteration = 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_iteration =  jacobi_iteration + 1;
    end
 end

 toc;
 iteration_time = toc;
%Plotting 
figure(1)
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('y axis')
hold on
title_text_1 = sprintf('Steady State - JACOBI METHOD ');
title_text_2 = sprintf('iteration number=%d', jacobi_iteration);
title({title_text_1,title_text_2});
end




%GAUSS-SEIDAL method
%inputs for the problem

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;

%Boundary conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;

T=ones(nx,ny);

% calculating the 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;


% Calculating average temperatures at the matrix corners.
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));


[X,Y] = meshgrid(x,y);
Told = T; 

% gauss-seidal Iteration.
for iterative_solver = 2
tic;
 if iterative_solver == 2
    gauss_seidal_iteration = 2;
    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;
        gauss_seidal_iteration =  gauss_seidal_iteration + 1;
    end
 end

 toc;
 iteration_time = toc;
%Plotting 
figure(2)
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('y axis')
hold on
title_text_1 = sprintf('Steady State - GAUSS SEIDAL METHOD ');
title_text_2 = sprintf('iteration number=%d', gauss_seidal_iteration);
title({title_text_1,title_text_2});
end




%SUCESSIVE OVER RELAXATION method
%inputs for the problem

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;
alpha = 1.3; 
omega = 1.2;

%Boundary conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;

T=ones(nx,ny);

% calculating the 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;


% Calculating average temperatures at the matrix corners.
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));


[X,Y] = meshgrid(x,y);
Told = T; 

% sor Iteration.
for iterative_solver = 3
tic;
 if iterative_solver == 3
    sor_iteration = 1;
    while(error > tol) 
        for i = 2:nx-1 
            for j = 2:ny-1
                T(i,j) = ((1-omega)*Told(i,j) + (omega*0.25*(T(i-1,j) + Told(i+1,j) + T(i,j-1) + Told(i,j+1))));
            end
        end
        error = max(max(abs(Told - T)));
        Told = T;
        sor_iteration =  sor_iteration + 1;
    end
 end

 toc;
 iteration_time = toc;
%Plotting 
figure(3)
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('y axis')
hold on
title_text_1 = sprintf('Steady State - SUCESSIVE OVER RELAXATION METHOD ');
title_text_2 = sprintf('iteration number=%d', sor_iteration);
title({title_text_1,title_text_2});
end

EXPLANATION:

• First of all we need to provide the initial values and given datas in the problem.
• Here consider the error=9e9.
• After that need to provide the Boundary conditions which are given in the problem.
• Now using the above boundary conditions we need to calculate the corner temperatures by averaging the boundary values.
• Then provide the K1 and K2 values which we calculated.
• Now we are initializing the command for the meshgrid.
• Then with reference to the coding steps, write the correct sequence of codes and update the old values for each timestep and it is mandatory.
• after when we came to plot command, we want to use coutourf command to acquire the plot as mentioned in the challenge.
• colorbar command is used to display a vertical bar to represent the temperature.
• clabel command this is used to represent the values of temparature inside the contour plot.
• Finally we have to code a line to give title for the plot.
• Follow these steps for various methods like, 
• Gauss seidal method(iterative solver==2)
• Sor method(iterative solver==3) and omega=1.2.                                        

OUTPUT:

1)Jacobi Method:-

number of iteration = 869

2)Gauss Seidal Method:-

number of iteration = 464

3)Sucessive Over Relaxation Method:-

number of iteration = 318.

Results and discussion:

                                        From the above plots we came to know that the number of iterations for, Sucessive Over Relaxation method is less when compared to gauss seidal and jacobi method.

                                        Jacobi > gauss seidal >sucessive over relaxation.

                                        Hence proved.

Unsteady/Transient state:

 As mentioned in the challenge we have to do this for two methods namely,

• Explicit
• Implicit

The two dimensional unsteady heat conduction equation is given below,

`(delT)/(delt)=alpha*((del^2T)/(delx^2)+(del^2T)/(dely^2))`

Wehen we discritize the above equation, we want to use forward scheme for LHS and central difference schemefor RHS,

L.H.S

`(delT)/(delt)=(T_(i,j)^(n+1)-T_(i,j)^n)/(Deltat)`

R.H.S:

`(del^2T)/(delx^2)=(T_(i-1,j)-2*T_(i,j)+T_(i+1,j))/(Deltax^2)`

`(del^2T)/(dely^2)=(T_(i,j-1)-2*T_(i,j)+T_(i,j+1))/(Deltay^2)`

Then, the equation becomes,

`(T_(i,j)^(n+1)-T_(i,j)^n)/(Deltat)=alpha*[((T_(i-1,j)-2*T_(i,j)+T_(i+1,j))/(Deltax^2))+((T_(i,j-1)-2*T_(i,j)+T_(i,j+1))/(Deltay^2))]`

Then to make the above equation simple we are providing some constants,

`K_1=alpha*((Deltat)/(Deltax^2))`

`K_2=alpha*((Deltat)/(Deltay^2))`

Hence,

`T_(i,j)^(n+1)=(T_(i,j)^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)))/(1+2*K_1+2*K_2)`

In unsteady condition the temprature changes with time and reaches a steady state.

EXPLIXCIT METHOD:

                         The main program to solve for the explicit method of two dimensional trasient equation is follows,

%Solving UnSteady state 2-dimensional equation(explicit method)
clear all
close all
clc


%JACOBI method
%inputs for the problem

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;
alpha=2.3e-5;
dt=30;

%Boundary conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;

T=ones(nx,ny);

% calculating the 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;


% Calculating average temperatures at the matrix corners.
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 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);


[X,Y] = meshgrid(x,y);
Told = T; 

% Jacobi Iteration.
for iterative_solver = 1
tic;
 if iterative_solver == 1
    iteration = 1;
    while(error > tol) 
        for i = 2:nx-1 
            for j = 2:ny-1
               term1 = Told(i,j);
               term2 = k1*(Told(i+1,j) - 2*Told(i,j) + Told(i-1,j));
               term3 = k2*(Told(i,j+1) - 2*Told(i,j) + Told(i,j-1));
               T(i,j) = (term1 + term2 + term3);
            end
        end
        error = max(max(abs(Told - T)));
        Told = T;
        iteration = iteration + 1;
    end
 end

 toc;
 iteration_time = toc;
%Plotting 
figure(1)
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('y axis')
hold on
title_text_1 = sprintf('UnSteady State - EXPLICIT METHOD ');
title_text_2 = sprintf('iteration number=%d',iteration);
title({title_text_1,title_text_2});


end

EXPLANATION:

• First of all we need to provide the initial values and given datas in the problem.
• Here consider the error=9e9.
• After that need to provide the Boundary conditions which are given in the problem.
• Now using the above boundary conditions we need to calculate the corner temperatures by averaging the boundary values.
• Then provide the K1 and K2 values which we calculated.
• Now we are initializing the command for the meshgrid.
• Then with reference to the coding steps, write the correct sequence of codes and update the old values for each timestep and it is mandatory.
• after when we came to plot command, we want to use coutourf command to acquire the plot as mentioned in the challenge.
• colorbar command is used to display a vertical bar to represent the temperature.
• clabel command this is used to represent the values of temparature inside the contour plot.
• Finally we have to code a line to give title for the plot.

 OUTPUT:

RESULTS AND DISCUSSIONS:

                                                  number of iteration = 872.

IMPLICIT METHOD:

                            The main program to solve for the implicit method of two dimensional trasient equation is follows,

%Solving Steady state 2-dimensional equation
clear all
close all
clc


%JACOBI method
%inputs for the problem

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;
dt=30;
alpha=1.3;


%Boundary conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;

T=ones(nx,ny);

% calculating the 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;


% Calculating average temperatures at the matrix corners.
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 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);


[X,Y] = meshgrid(x,y);
Told = T; 
T_prev_dt = Told;
iterative_solver = 1;

% Jacobi method
term1 = (1+ 2*k1 + 2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
tic;
 if iterative_solver == 1
    jacobi_iteration = 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_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
            end    
        end
        error = max(max(abs(Told - T)));
        Told = T;
        jacobi_iteration =   jacobi_iteration + 1;
    end
    T_prev_dt = T;
 end 
toc;
time_required = toc;
%Plotting 
figure(1)
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('y axis')
title_text_1 = sprintf(' IMPLICT - JACOBI METHOD');
title_text_2 = sprintf('iteration number=%d', jacobi_iteration);
title({title_text_1,title_text_2});




%GAUSS-SEDIAL METHOD
%inputs for the problem

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;
dt=30;
alpha=1.3;


%Boundary conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;

T=ones(nx,ny);


% calculating the 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;


% Calculating average temperatures at the matrix corners.
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 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);


[X,Y] = meshgrid(x,y);
Told = T;
T_prev_dt = Told;
iterative_solver = 2;

%Gauss seidal method
term1 = (1+ 2*k1 + 2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
tic,
  if iterative_solver == 2
    gauss_seidal_iter = 1;
    while(error > tol) 
        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;
       gauss_seidal_iter = gauss_seidal_iter + 1;
    end
    T_prev_dt = Told;
 end 
toc;
time_required = toc;
%Plotting 
figure(2)
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('y axis')
title_text_1 = sprintf('IMPLICIT - GAUSS SEDIAL METHOD');
title_text_2 = sprintf('iteration number=%d', gauss_seidal_iter);
title({title_text_1,title_text_2});  








%SUCESSIVE OVER RELAXATION METHOD
%inputs for the problem

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;
dt=30;
alpha=1.3;
omega =1.2;


%Boundary conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;

T=ones(nx,ny);


% calculating the 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;


% Calculating average temperatures at the matrix corners.
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 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);


[X,Y] = meshgrid(x,y);
Told = T;
T_prev_dt = Told;
iterative_solver = 3;

%SOR method
term1 = (1+ 2*k1 + 2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
tic,
  if iterative_solver == 3
    SOR_iteration = 1;
    while(error > tol) 
        for i = 2:nx-1 
            for j = 2:ny-1
               H = (omega*T(i-1,j) + (1-omega)*Told(i-1,j) + Told(i+1,j));
               V = (omega*T(i,j-1) + (1-omega)* 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;
       SOR_iteration = SOR_iteration + 1;
    end
    T_prev_dt = Told;
 end 
toc;
time_required = toc;
%Plotting 
figure(3)
[a,b] = contourf(X,Y,T);
colorbar;
clabel(a,b);
xlabel('x axis')
ylabel('y axis')
title_text_1 = sprintf('IMPLICIT - SUCESSIVE OVER RELAXATION METHOD');
title_text_2 = sprintf('iteration number=%d', SOR_iteration);
title({title_text_1,title_text_2}); 

EXPLANATION AND RESULTS:

• First of all we need to provide the initial values and given datas in the problem.
• Here consider the error=9e9.
• After that need to provide the Boundary conditions which are given in the problem.
• Now using the above boundary conditions we need to calculate the corner temperatures by averaging the boundary values.
• Then provide the K1 and K2 values which we calculated.
• Now we are initializing the command for the meshgrid.
• Then with reference to the coding steps, write the correct sequence of codes and update the old values for each timestep and it is mandatory.
• after when we came to plot command, we want to use coutourf command to acquire the plot as mentioned in the challenge.
• colorbar command is used to display a vertical bar to represent the temperature.
• clabel command this is used to represent the values of temparature inside the contour plot.
• Finally we have to code a line to give title for the plot.
• In the SOR method  we take alpha=1.3 which is the relaxation factor.

OUTPUT:

1)Jacobi Method:-

number of iteration = 868

2)Gauss Seidal Method:-

number of iteration = 462

3)Sucessive Over Relaxation Method:-

number of iteration = 376.