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

Aim : Solving the steady and unsteady 2D heat conduction problem

Objective : 

Solve the 2D heat conduction equation by using the point iterative techniques that were taught in the class. The Boundary conditions for the problem are as follows;

Top Boundary = 600 K

Bottom Boundary = 900 K

Left Boundary = 400 K

Right Boundary = 800 K

Theory:

Three dimension heat conduction equation is given as :

`(deltaT)/(deltat) = alpha((delta^2T)/(deltax^2)+(delta^2T)/(deltay^2)+(delta^2T)/(deltaz^2))`

For steady state 2D problem the above equation becomes: `(deltaT)/(deltat) = 0`

For unsteady 2D state the equation becomes : `(deltaT)/(deltat) = alpha((delta^2T)/(deltax^2)+(delta^2T)/(deltay^2))`

1. 2D steady state problem:

We apply central difference method for the above steady state equation. We get:

`T(i) =  (Deltax^2Deltay^2)/(2(Deltax^2+Deltay^2))((T(i+1,j)+T(i-1,j))/(Deltax^2)+(T(i,j+1)+T(i,j-1))/(Deltay^2))`

Here we assign : k = `(2(Deltax^2+Deltay^2))/(Deltax^2Deltay^2)`

Methods used:

1.jacobian method: In jacobian method we use values from previous iterations therefore the right side is already knowsn.Here we use values from nth time step while calculating values of n+1 time.

`(T(i,j))^(n+1) =  (1/k).((T(i+1,j)+T(i-1,j))/(Deltax^2)+(T(i,j+1)+T(i,j-1))/(Deltay^2))^(n)`

2.Gauss sidel method : We update values after each iterations and we use the update values.

`(T(i,j))^(n+1) =  (1/k).((T(i+1,j)+T(i-1,j))/(Deltax^2)+(T(i,j+1)+T(i,j-1))/(Deltay^2))^(n+1)`

3.Successive over relaxation: It is used either over jacobian or Gauss sidel method it enhances the performance of other iterative solver. It uses the scaling factor omega. It either makes solution to converge faster or slower based on the value of Omega. If we choose the value of omega above treshold values may shoot over the required results and would never converege. Therefore we need to be carefull while choosing the value of omega. 

`T(i,j) = Told(i,j) +omega(T(i,j)-Told(i,j))`

Where T(i,j) are the updated values from Gauss sidel method.

Steps in coding:

1. First we input the values of length of x and y coordinates.
2. We enter the number of grid points which are same in both x and u direction
3. Create temperature array and set temperature values.
4. We set tollerance value and error value.
5. We set SOR factor value.
6. We create for loop for each interval of time and inside we create space coordinate loops.
7. after each iteration we compare the values of tollerance and error if the error is more than tollerance the loop continues.
8. We plot the contour plot of x,y and T

Main code:

clear all
close all
clc

%inputs
L = 1;
nx = 20;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
[X, Y]=meshgrid(x,y);
T = 300*(ones(nx,ny));
T(1,:)=600;
T(ny,:)=900;
T(:,1)=400;
T(:,ny)=800;
Told = T;
T_new = T;
dx = L/(nx-1);
dy = L/(ny-1);
k = (2*(dx^2+dy^2))/((dx^2)*(dy^2));
q = k^-1;
q1 = q/dx^2;
q2 = q/dy^2;
e = 8e888;
tol = 1e-4;
e = 1;
omega = 1.73;
iteration_method = 1;

%jacobian method
if iteration_method == 1
    tic
iteration =1;
while e > tol
for i = 2:(nx-1)
    for j = 2:(ny-1)
        
TR = Told(i+1,j);
TL = Told(i-1,j);
Tt = Told(i,j+1);
Tb = Told(i,j-1);

T(i,j) = q1*(TR+TL)+q2*(Tt+Tb);
    end
end
error = max(abs(T-Told));
e = max(abs(error));
Told = T;
iteration = iteration+1;
end
timetaken = toc;
[m, n]=contourf(X,Y,T);
clabel(m,n);
colormap(jet);
colorbar
ax = gca;
ax.YDir = 'reverse';
title_text = sprintf('Number of iteration = %d',iteration);
title({['Jacobian method'];title_text;['time taken = ',num2str(timetaken)]})
xlabel('x-axis')
ylabel('y-axis')
zlabel('Temperature')


%Gauss sidel method
elseif iteration_method == 2
    tic
iteration =1;
while e > tol
for i = 2:(nx-1)
    for j = 2:(ny-1)
        
TR = T_new(i+1,j);
TL = T_new(i-1,j);
Tt = T_new(i,j+1);
Tb = T_new(i,j-1);

T(i,j) = q1*(TR+TL)+q2*(Tt+Tb);
T_new = T;
    end
end
error = max(abs(T-Told));
e = max(abs(error));
Told = T;
iteration = iteration+1;
end
timetaken = toc;
[m, n]=contourf(X,Y,T);
clabel(m,n);
colormap(jet);
colorbar
ax = gca;
ax.YDir = 'reverse';
title_text = sprintf('Number of iteration = %d',iteration);
title({['Gauss sidel method'];title_text;['time taken = ',num2str(timetaken)]})
xlabel('x-axis')
ylabel('y-axis')
zlabel('Temperature')

%successive over relaxation method applyed to Gauss sidel


elseif iteration_method == 3
tic
iteration =1;
while e > 1e-4
for i = 2:(nx-1)
    for j = 2:(ny-1)
        
TR = T_new(i+1,j);
TL = T_new(i-1,j);
Tt = T_new(i,j+1);
Tb = T_new(i,j-1);

T(i,j) = q1*(TR+TL)+q2*(Tt+Tb);
T(i,j) = Told(i,j)+omega*(T(i,j)-Told(i,j));
T_new = T;
    end
end
error = max(abs(T-Told));
e = max(abs(error));
Told = T;
iteration = iteration+1;
end
timetaken = toc
[m, n]=contourf(X,Y,T);
clabel(m,n);
colormap(jet);
ax = gca;
ax.YDir = 'reverse';
colorbar
title_text = sprintf('Number of iteration = %d',iteration);
title({['SOR method'];title_text;['time taken = ',num2str(timetaken)]})
xlabel('x-axis')
ylabel('y-axis')
zlabel('Temperature')

end

Results:

Conclusion : Based on number of iterations

Jacobian > Gauss sidel > SOR

Therefore we can see that jacobian takes lot of iterations to arrive at the solution where as SOR takes least and Gauss sidel lies in between them.

2. 2D unsteady state:

For 2D unsteady state transient 3D heat conduction equation becomes:

`(deltaT)/(deltat) = alpha((delta^2T)/(deltax^2)+(delta^2T)/(deltay^2))`

The above equation can be discretized in two ways:

1. Explicit

2.Implicit

1.Explicit: In this method the temperature values of previous time steps are inputed on the right side i.e:

`T(i,j)^(n+1) = T(i,j)^(n)+k1(T(i+1,j)+T(i-1,j))^n+k2(T(i,j+1)+T(i,j-1))^n`

where:

`k1 = (alphaDeltat)/(Deltax^2` `k2 = (alphaDeltat)/(Deltay^2`

`alpha =`Thermal diffusivity

 Here we know the all the values on the right hand side.

Main code:

clear all
close all
clc

nx = 20;
ny = 20;
dt = 1e-4;
L = 1;
alpha = 1.4;
dx= L/(nx-1); 
dy= L/(ny-1);
x = linspace(0,L,nx);
y = linspace(0,L,ny);
[X Y] = meshgrid(x,y);
T = 300*(ones(nx,ny));
T(1,:)=600;
T(ny,:)=900;
T(:,1)=400;
T(:,ny)=800;
Told = T;
tol = 1e-4;
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
Term1 = (1-2*k1-2*k2);
Term2 = k1;
Term3 = k2;
iteration = 0;
error = 8e88;   

   
        while error > tol
            for i = 2:(nx-1)
                for j = 2:(ny-1)
                    T(i,j)=Told(i,j)*Term1+Term2*(Told(i+1,j)+Told(i-1,j))+Term3*(Told(i,j+1)+Told(i,j-1));
                end
            end
           
            error = max(max(abs(T-Told)));
             Told = T;
             iteration = iteration+1;
                      
        end

 [m, n]=contourf(X,Y,T);
clabel(m,n);
colormap(jet);
colorbar
ax = gca;
ax.YDir = 'reverse';
title_text = sprintf('Number of iteration = %d',iteration);
title({['explicit method for Transient heat conduction'];[title_text]})
xlabel('x-axis')
ylabel('y-axis')
zlabel('Temperature')

Result:

2.Implicit method:

`T(i,j)^(n+1) = T(i,j)^(n)+k1(T(i+1,j)+T(i-1,j))^(n+1)+k2(T(i,j+1)+T(i,j-1))^(n+1)`

where:

`k1 = (alphaDeltat)/(Deltax^2` `k2 = (alphaDeltat)/(Deltay^2`

`alpha =`Thermal diffusivity

Here the right side values are unknown.

We have 3 iteration methods to solve for explicit type.

1. jacobian method

2.gauss sidel method

3.Successive over relaxation method

Main code:

clear all
close all
clc

nx = 20;
ny = 20;
nt = 1400;
dt = 0.01;
L = 1;
alpha = 1;
dx= L/(nx-1); 
dy= L/(ny-1);
x = linspace(0,L,nx);
y = linspace(0,L,ny);
[X Y] = meshgrid(x,y);
T = 300*(ones(nx,ny));
T(1,:)=600;
T(ny,:)=900;
T(:,1)=400;
T(:,ny)=800;
Told = T;
Tprevious = T;
tol = 1e-4;
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
Term1 = ((1+2*k1+2*k2))^-1;
Term2 = k1*Term1;
Term3 = k2*Term1;
iteration = 0;
solver = 1;

if solver == 1
   
    for k = 1:nt
        error = 8e88
        while error > tol
            for i = 2:(nx-1)
                for j = 2:(ny-1)
                    T(i,j)=Tprevious(i,j)*Term1+Term2*(Told(i+1,j)+Told(i-1,j))+Term3*(Told(i,j+1)+Told(i,j-1));
                end
            end
           
            error = max(max(abs(T-Told)));
             Told = T;
             iteration = iteration+1;
            
        end
        Tprevious = T;
[m, n]=contourf(X,Y,T);
clabel(m,n);
colormap(jet);
colorbar
ax = gca;
ax.YDir = 'reverse';
title_text = sprintf('Number of iteration = %d',iteration);
title({['implicit jacobian method for Transient heat conduction'];[title_text]})
xlabel('x-axis')
ylabel('y-axis')
zlabel('Temperature')
    end
end

if solver == 2
   T_update = T;
    for k = 1:nt
        error = 8e88
        while error > tol
            for i = 2:(nx-1)
                for j = 2:(ny-1)
                    T(i,j)=Tprevious(i,j)*Term1+Term2*(T_update(i+1,j)+T_update(i-1,j))+Term3*(T_update(i,j+1)+T_update(i,j-1));
                end
                T_update = T;
            end
           
            error = max(max(abs(T-Told)));
             Told = T;
             iteration = iteration+1;
            
        end
        Tprevious = T;
[m, n]=contourf(X,Y,T);
clabel(m,n);
colormap(jet);
colorbar
ax = gca;
ax.YDir = 'reverse';
title_text = sprintf('Number of iteration = %d',iteration);
title({['implicit Gauss sidel method for Transient heat conduction'];[title_text]})
xlabel('x-axis')
ylabel('y-axis')
zlabel('Temperature')
    end
end

%SOR method


if solver == 3
   T_update = T;
   omega = 1.2;
    for k = 1:nt
        error = 8e88
        while error > tol
            for i = 2:(nx-1)
                for j = 2:(ny-1)
                    T(i,j)=Tprevious(i,j)*Term1+Term2*(T_update(i+1,j)+T_update(i-1,j))+Term3*(T_update(i,j+1)+T_update(i,j-1));
                    T(i,j)=Told(i,j)+omega*(T(i,j)-Told(i,j));
                end
                T_update = T;
            end
           
            error = max(max(abs(T-Told)));
             Told = T;
             iteration = iteration+1;
            
        end
        Tprevious = T;
[m, n]=contourf(X,Y,T);
clabel(m,n);
colormap(jet);
colorbar
ax = gca;
ax.YDir = 'reverse';
title_text = sprintf('Number of iteration = %d',iteration);
title({['implicit SOR method for Transient heat conduction'];[title_text]})
xlabel('x-axis')
ylabel('y-axis')
zlabel('Temperature')
    end
end

Results:

Conclusion : SOR method converges faster then the gauss sidel and jacobian method.

Number of iterations: Jacobian>Gauss sidel>SOR