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

Aim : The aim of this project is to solve steady and unsteady 2D heat conduction equation.

Theory : The heat conduction equation is a partial differential equation that describes the distribution of heat (or the temperature field) in a given body over time.

Assumption & Boundary conditions :

Following assumption are made to solve the steady and unsteady equation which as follows:-

1. Assuming that the domain is a unit square.

2. Assuming that nx = ny [Number of points along the x direction is equal to the number of points along the y direction]

3. Boundary conditions for steady and unsteady case

• Left:400K
• Top:600K
• Right:800K
• Bottom:900K

4. Tolerance = 1e-4

Solution :

A. Solving Steady state 2D heat conduction equation :

Govering equation for 2D heat conduction equation for steady state is

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

Discretizing the above equation using Central differencing scheme(CDS), we get

`((T_L-2*T_P+T_R)/(Deltax^2) + (T_T-2*T_P+T_B)/(Deltay^2))=0`

Futher simplifying the above equation, we get

`T_P =(1/K)*((T_L+T_R)/(Deltax^2)+(T_T+T_B)/(Deltay^2))`

`K =2*((Deltax^2+Deltay^2)/(Deltax^2*Deltay^2))`

1. Jacobi Method

`T_P^(n+1) =(1/K)*((T_L+T_R)/(Deltax^2)+(T_T+T_B)/(Deltay^2))`

2. Guass seidel Method

`T_P^(n+1) =(1/K)*((T_L+T_R)/(Deltax^2)+(T_T+T_B)/(Deltay^2))^(n+1)`

3. Successive over relaxation Method

`T_P^(n+1) =T_P^n*(1-omega)+omega*((1/K)*((T_L+T_R)/(Deltax^2)+(T_T+T_B)/(Deltay^2))^(n+1))`

Main Code :

clear all
close all
clc

%Solving the steady 2D heat conduction equation

%Input parameter
nx=10; %no of grids in x-direction
ny=10; %no of grids in y_direction

%Calculation
x=linspace(0,1,nx);
y=linspace(0,1,ny);

dx=x(2)-x(1); %Time step in x-direction
dy=y(2)-y(1); %Time step in y-direction

%Creating temperatures across given region
T=300*ones(nx,ny);
T(:,1)=400;
T(1,:)=600;
T(:,end)=800;
T(end,:)=900;

%Temperature at edges
T(1,1)=(600+400)/2;
T(1,end)=(600+800)/2;
T(end,1)=(400+900)/2;
T(end,end)=(800+900)/2;

% Constant
k=2*((dx^2+dy^2)/(dx^2*dy^2)); %For steady state
alpha=1.4;
dt=1e-3;
k1=alpha*(dt/(dx^2)); %CFLx for transient state
k2=alpha*(dt/(dy^2)); %CFLy for transient state
cfl=k1+k2;

for method=linspace(1,7,7)
if method == 1
    error_1=steady_jacobi(x,y,nx,ny,dx,dy,T,k);
    fprintf('\n Error using Jacobi method(steady state)=%d',error_1)
    
elseif method == 2
    error_2=steady_gs(x,y,nx,ny,dx,dy,T,k);
    fprintf('\n\n Error using Gauss seidel method(steady state)=%d',error_2)
    
elseif method == 3
    error_3=steady_sor(x,y,nx,ny,dx,dy,T,k);
    fprintf('\n\n Error using SOR method(steady state)=%d',error_3)
elseif method == 4
    error_4=transient_explicit(x,y,nx,ny,T,k1,k2);
    fprintf('\n\n Error using Explicit transient state =%d',error_4) 

elseif method == 5
    error_5=transient_implicit_jacobi(x,y,nx,ny,T,k1,k2);
    fprintf('\n\n Error using Jacobi method(implicit transient state)=%d',error_5)
    
elseif method == 6
    error_6=transient_implicit_gs(x,y,nx,ny,T,k1,k2);
    fprintf('\n\n Error using Gauss seidel method(implicit transient state)=%d',error_6)
    
elseif method == 7
    error_7=transient_implicit_sor(x,y,nx,ny,T,k1,k2);    
    fprintf('\n\n Error using SOR method(implicit transient state)=%d',error_7)     
    
end
end

Explanation :

1. Initally the number of nodes, boundary conditions (temperatures at sides and edges) are assigned.
2. x, dx, y, dy and dt values are assigned.
3. The constants k (for steady state), alpha, k1,k2 (CFL numbers for transient analysis) are assigned.
4. Using for loop and elseif commands, each method (steady and transient state) is executed one by one.
5. Under function steady_jacobi, the error, tol, Told (old temperature) and initial iteration value are asigned.
6. Convergence loop (while loop) is executed in which T(i,j) is calculated under for loop, and after for loop the error value and iteration number are updated and error value is calculated.
7. While loop comes to end when error falls behind tolerance value and the T(i,j) is plotted using contourf and colormap commands.
8. No of iterations is displayed on title. When the function ends, it returns error value to main program.
9. The above procedure is repeated for steady states Gauss seidel and SOR methods in which only Temperature T(i,j) value is altered for each methods.
10. In explicit method, iterative solver (convergence loop) is not used. Time loop (for nt=1:1400) is used, then temperature T(i,j) is found after which error and old temperature are updated. When time loop ends, the T(i,j) is plotted.
11. In implicit scheme, the Told, Tpre (previous temp) and iteration number are assigned. Then the time loop (for nt=1:1400) is used. Then error and tolerance value are assigned.
12. Convergence loop is used and T(i,j) value is calculated, then error value is found out and Told and iteration number are updated. When while loop comes to end, the Tpre (previous temp) is updated (after each time loop).
13. After 1400 time step, time loop comes to end and the temperature T(i,j) is plotted. Same procedure is repeated for Gauss seidel and SOR methods (implicit scheme) with their respective T(i,j) values are calculated and plotted.

Function code :

Steady State

1. Jacobi Method

% Steady state (Jacobi)
function [error]=steady_jacobi(x,y,nx,ny,dx,dy,T,k)

error=9e9; 
tol=1e-4;
Told=T;
jacobi_iter=1;
tic
while(error>tol)
    for i=2:nx-1
        for j=2:ny-1
        term1=(Told(i-1,j)+Told(i+1,j))/(k*(dx^2));
        term2=(Told(i,j-1)+Told(i,j+1))/(k*(dy^2));
        T(i,j)=term1+term2;
        end
    end
    
    error=max(max(abs(Told-T)));
    Told=T;
    jacobi_iter=jacobi_iter+1;
end 
 time=toc;
    figure(1)
    meshgrid(x,y);
    [a, b]=contourf(x,y,T);
    clabel(a,b)
    colorbar 
    colormap(jet)
    set(gca,'ydir','reverse')
    xlabel('X Domain Length')
    ylabel('Y Domain Length')
    title_text1=sprintf('No of iterations (steady state jacobi)=%d',jacobi_iter);
    title_text2=sprintf('Simulation time =%f s',time);
    title(title_text1,title_text2)  
    
end

2. Guass seidel Method

% Steady state (Gauss seidel)
function [error]=steady_gs(x,y,nx,ny,dx,dy,T,k)

error=9e9; 
tol=1e-4;
Told=T;
gs_iter=1;
tic
while(error>tol)
    for i=2:nx-1
        for j=2:ny-1
        term1=(T(i-1,j)+Told(i+1,j))/(k*dx^2);
        term2=(T(i,j-1)+Told(i,j+1))/(k*dy^2);
        T(i,j)=term1+term2;
        end
    end
    
    error=max(max(abs(Told-T)));
    Told=T;
    gs_iter=gs_iter+1;
end
time=toc;
    figure(2)
    meshgrid(x,y);
    [a, b]=contourf(x,y,T);
    clabel(a,b)
    colorbar 
    colormap(jet)
    set(gca,'ydir','reverse')
    xlabel('X Domain Length')
    ylabel('Y Domain Length')
    title_text1=sprintf('No of iterations (steady state gauss seidel)=%d',gs_iter);
    title_text2=sprintf('Simulation time =%f s',time);
    title(title_text1,title_text2)  
end

3. Successive over relaxation Method

% Steady state (SOR)
function [error]=steady_sor(x,y,nx,ny,dx,dy,T,k)
error=9e9; 
tol=1e-4;
omega=1.2; %Relaxation factor for SOR
Told=T;
sor_iter=1;
tic
while(error>tol)
    for i=2:nx-1
        for j=2:ny-1
        term1=(T(i-1,j)+Told(i+1,j))/(k*dx^2);
        term2=(T(i,j-1)+Told(i,j+1))/(k*dy^2);
        T(i,j)=((Told(i,j))*(1-omega))+(omega*(term1+term2));
        end
    end

    error=max(max(abs(Told-T)));
    Told=T;
    sor_iter=sor_iter+1;
end
 time=toc;
    figure(3)
    meshgrid(x,y);
    [a, b]=contourf(x,y,T);
    clabel(a,b)
    colorbar 
    colormap(jet)
    set(gca,'ydir','reverse')
    xlabel('X Domain Length')
    ylabel('Y Domain Length')
    title_text1=sprintf('No of iterations (steady state SOR)=%d',sor_iter);
    title_text2=sprintf('Simulation time =%f s',time);
    title(title_text1,title_text2)
end

Result :

Observation :

• Tough we are getting same result for different iterative method so, the different between them is number of iteration and computation time taken by iterative solver to get convergence result.
• By observing,above plots its quite clear that jacobi_iteration takes more number of iteration and computation time as compared to other two iterative solver this is because jacobi_iteration takes the old values.
• While other two iterative solver(Gauss seidel & Successive over relaxiation) takes the updated values.Hence the new values is passed in new iteration as the solution proceeds.Thus in this type of iterative solver we get convergence result in less number of iteration as well as in less computation time.

2. Solving 2D heat conduction equation for unsteady state :

Govering equation for 2D heat conduction equation for unsteady state is

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

We will solve above equation by using two scheme,

• Explict scheme
• Implicit scheme

A. Explict scheme :

Discretising the above equations using Forward differencing scheme(FDS) for time derivative and Central differencing scheme(CDS) for space derivative, we get,

`(T_P^(n+1) -T_P^n)/(Deltat)=alpha*((T_L-2*T_P+T_R)/(Deltax^2) + (T_T-2*T_P+T_B)/(Deltay^2))^n`

`T_P^(n+1) =T_P^n+K_1*(T_L-2T_P+T_R)^n+K_2*(T_T-2T_P+T_B)^n`

`T_P^(n+1) =T_P^n*(1-2K_1+2K_2)+K_1*(T_L+T_R)^n+K_2*(T_T+T_B)^n`

where `K_1=alpha*(delT)/(delx^2)` and `K_2=alpha*(delT)/(dely^2)`

where `CFL_x = K_1` and `CFL_y=K_2`

For explicit diffusion problems, `CFL_x +CFL_y <= 0.5` for the solution to be stable.

Code :

function [error]=transient_explicit(x,y,nx,ny,T,k1,k2)

Told=T;
jacobi_iter=1;
tic

for nt=1:1400
    for i=2:nx-1
        for j=2:ny-1
        term1=(Told(i-1,j)+Told(i+1,j));
        term2=(Told(i,j-1)+Told(i,j+1));
        T(i,j)=(Told(i,j)*(1-(2*k1)-(2*k2)))+(k1*term1)+(k2*term2);
        end
    end
    error=max(max(abs(Told-T)));
    Told=T;
    jacobi_iter=jacobi_iter+1;
end  
    time=toc;
    figure(4)
    [a, b]=contourf(x,y,T);
    clabel(a,b)
    colorbar 
    colormap(jet)
    set(gca,'ydir','reverse')
    xlabel('X Domain Length')
    ylabel('Y Domain Length')
    title_text1=sprintf('No of iterations (Explicit state)=%d',jacobi_iter);
    title_text2=sprintf('Simulation time =%f s',time);
    title(title_text1,title_text2)
   
end

Result :

B. Implicit scheme :

Govering equation for 2D heat conduction equation for unsteady state is,

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

Discretising the above equations using Forward differencing scheme(FDS) for time derivative and Central differencing scheme(CDS) for space derivative, we get,

`(T_P^(n+1) -T_P^n)/(Deltat)=alpha*((T_L-2*T_P+T_R)/(Deltax^2) + (T_T-2*T_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))`

1. Jacobi Method

`T_P^(n+1) =1/((1-2K_1+2K_2))*(T_P^n+K_1*(T_L+T_R)^(n)+K_2*(T_T+T_B)^(n))`

2. Guass seidel Method

`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))`

3. Successive over relaxation Method

`T_P^(n+1) =T_P^n*(1-omega)+omega*(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)))`

Code :

Implicit scheme :

1. Jacobi Method

% Transient state implcit scheme (Jacobi)
function [error]=transient_implicit_jacobi(x,y,nx,ny,T,k1,k2)

Told=T;
Tpre=T;
jacobi_iter=1;
tic

for nt=1:1400
    error=9e9;
    tol=1e-4;
 while(error>tol)
    for i=2:nx-1
        for j=2:ny-1
        term1=1/(1+(2*k1)+(2*k2));
        term2=k1*term1;
        term3=k2*term1;
        H=(Told(i-1,j)+Told(i+1,j));
        V=(Told(i,j-1)+Told(i,j+1));
        T(i,j)=(Tpre(i,j)*term1)+(H*term2)+(V*term3);
        end
    end
    
    error=max(max(abs(Told-T)));
    Told=T;
    jacobi_iter=jacobi_iter+1;
 end
Tpre=T;
end
time=toc;
    figure(5)
    meshgrid(x,y);
    [a, b]=contourf(x,y,T);
    clabel(a,b)
    colorbar 
    colormap(jet)
    set(gca,'ydir','reverse')
    xlabel('X Domain Length')
    ylabel('Y Domain Length')
    title_text1=sprintf('No of iterations (transient implicit jacobi)=%d',jacobi_iter);
    title_text2=sprintf('Simulation time =%f s',time);
    title(title_text1,title_text2)
end

2. Guass seidel Method

function [error]=transient_implicit_gs(x,y,nx,ny,T,k1,k2)

Told=T;
Tpre=T;
gs_iter=1;
tic

for nt=1:1400
    error=9e9;
    tol=1e-4;
 while(error>tol)
    for i=2:nx-1
        for j=2:ny-1
        term1=1/(1+(2*k1)+(2*k2));
        term2=k1*term1;
        term3=k2*term1;
        H=(T(i-1,j)+Told(i+1,j));
        V=(T(i,j-1)+Told(i,j+1));
        T(i,j)=(Tpre(i,j)*term1)+(H*term2)+(V*term3);
        end
    end
    
    error=max(max(abs(Told-T)));
    Told=T;
    gs_iter=gs_iter+1;
 end
Tpre=T;
end  
time=toc;
    figure(6)
    meshgrid(x,y);
    [a, b]=contourf(x,y,T);
    clabel(a,b)
    colorbar 
    colormap(jet)
    set(gca,'ydir','reverse')
    xlabel('X Domain Length')
    ylabel('Y Domain Length')
    title_text1=sprintf('No of iterations (transient implicit gauss seidel)=%d',gs_iter);
    title_text2=sprintf('Simulation time =%f s',time);
    title(title_text1,title_text2)
end

3. Successive over relaxation Method

% Transient state implcit scheme (SOR)
function [error]=transient_implicit_sor(x,y,nx,ny,T,k1,k2)

Told=T;
Tpre=T;
sor_iter=1;
omega=1.2; %Relaxation factor for SOR
tic

for nt=1:1400
  error=9e9;
  tol=1e-4;
 while(error>tol)
    for i=2:nx-1
        for j=2:ny-1
        term1=1/(1+(2*k1)+(2*k2));
        term2=k1*term1;
        term3=k2*term1;
        H=(T(i-1,j)+Told(i+1,j));
        V=(T(i,j-1)+Told(i,j+1));
        T(i,j)=(1-omega)*(Told(i,j))+(omega*((Tpre(i,j)*term1)+(H*term2)+(V*term3)));
        end
    end
    
    error=max(max(abs(Told-T)));
    Told=T;
    sor_iter=sor_iter+1;
 end
Tpre=T;
end  
time=toc;
    figure(7)
    meshgrid(x,y);
    [a, b]=contourf(x,y,T);
    clabel(a,b)
    colorbar 
    colormap(jet)
    set(gca,'ydir','reverse')
    xlabel('X Domain Length')
    ylabel('Y Domain Length')
    title_text1=sprintf('No of iterations (transient implicit SOR)=%d',sor_iter);
    title_text2=sprintf('Simulation time =%f s',time);
    title(title_text1,title_text2)
end

Result :

Observation :

• Tough we are getting same result for different iterative method so, the different between them is number of iteration and computation time taken by iterative solver to get convergence result.
• By observing,above plots its quite clear that jacobi_iteration takes more number of iteration and computation time as compared to other two iterative solver this is because jacobi_iteration takes the old values.
• While other two iterative solver(Gauss seidel & Successive over relaxiation) takes the updated values.Hence the new values is passed in new iteration as the solution proceeds.Thus in this type of iterative solver we get convergence result in less number of iteration as well as in less computation time.

Command Window :

Explanation :

In each plot, the domain is considered as square of side length 1m. In Y direction the order is reversed (from 1 to 0) because the command "set(gca,'ydir','reverse')" is used. In title the method type and number of iterations are displayed. Each color in the plot denotes specific temperature value (in Kelvin).

S.no	Method	No. of iterations	Simulation Time	Error value
1.	Steady state (Jacobi)	208	0.004988 s	9.486811e-05
2.	Steady state (Gauss seidel)	112	0.004115 s	9.586452e-05
3.	Steady state (Successive over relaxation)	76	0.004214 s	8.264460e-05
4.	Transient state (Explicit scheme)	1401	0.009893 s	0
5.	Transient state (Implicit scheme, Jacobi)	3666	0.022560 s	1.250555e-12
6.	Transient state (Implicit scheme, Gauss seidel)	3117	0.019852 s	0
7.	Transient state (Implicit scheme, Successive over relaxation)	2991	0.020352 s	5.684342e-14

From the table, we can clearly see that the Jacobi method takes lot of iterations to converge and successive over relaxation method takes the least, whereas iterations took by gauss seidel lies between the other two methods.

Order of method based on iterations,

Jacobi>Gauss seidel>Successive over relaxation

Conclusion :

• Thus in matlab, the 2D heat conduction equation is solved by using the point iterative techniques using the methods like Jacobi, Gauss-seidel, Successive over-relaxation for both implicit and explicit schemes.
• The temperature distribution plots are plotted for each methods.
• The number of iterations and error values for each methods are calculated and displayed in the above table.