SIMULATION OF 2D HEAT CONDUCTION IN STEADY AND TRANSIENT CONDITIONS USING MATLAB

OBJECTIVE:

The objective of the project is to solve the 2D heat conduction equation in MATLAB using different iterative solving techniques available. And also to compare the results on the basis of number of iterations to converge, time taken for covergence for each technique. This will be helpful in choosing the optimal solver for this equation. The equation is a PDE and will be solved for both steady state and transient conditions.

Types of Iterative techniques:

1. JACOBI METHOD

2. GAUSS-SEIDAL METHOD

3. SUCCESSIVE OVER RELAXATION METHOD

All the three techniques will be used in both STEADY and TRANSIENT conditions. And for Transient conditions alone, the methods will be used both Implicitly and Explicitly.

2D Heat conduction equation:

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

T = Temperature in K

t = time in sec

x,y = Co-ordinates

`alpha` = Thermal diffusivity in m2/s

Domain Parameters:

The domain taken into consideration is a square,

Length of the domain = 1 m

Breadth of the domain = 1 m

No. of grid points along x axis = 20

No. of grid points along y axis = 20

Boundary conditions:

1. Left:400K

2. Top:600K

3. Right:800K

4. Bottom:900K

SOLVING THE 2D STEADY STATE HEAT CONDUCTION EQUATION:

The steady state can be specified as

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

Applying numerical scheme to each of the terms in the equation to solve it numerically,

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

Now,

Substituting those values in the original steady state equation,

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

Rearranging them,

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

Consider,

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

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

Now the above equation has to be solved in the previously mentioned three iterative techniques.

MATLAB main function:

clear all
close all
clc

%Dimensions of the domain
L = 1;
B = 1;

%No. of grid points
nx = 20;
ny = 20;

%Initiating T matrix
T = ones(nx,ny);
x = linspace(0,L,nx);
y = linspace(0,B,ny);

%BC
T(:,1) = 400;
T(:,ny) = 800;
T(1,:) = 600;
T(nx,:) = 900;

%Refining BC\'s at the corners
T(1,1) = (T(1,2) + T(2,1))/2;
T(1,end) = (T(1,end-1) + T(2,end))/2;
T(end,1) = (T(end-1,1) + T(2,end))/2;
T(end,end) = (T(end-1,end) + T(end,end-1))/2;

T_old = T;

dx = x(2) - x(1);
dy = y(2) - y(1);

error_tol = 1e-4;

%Specifying the iterative method

sprintf(\"1.Jacobi\\n2.Gauss seidal\\n3.SOR\")

method = input(\'Method = \');

if method == 1
      
%JACOBI METHOD    
tic;
[jacobi, jacobi_iter]  = jacobi_steady(T,T_old,error_tol,nx,ny,dx,dy);
time_jacobi = toc;

%Plotting
figure(1)
[X,Y] = contourf(x,fliplr(y),jacobi,\'ShowText\',\'on\');
xlabel(\'X axis\');
ylabel(\'Y axis\');
title_text1 = sprintf(\'2D Steady-state heat conduction-Jacobi\\nNo.of iterations to converge = %d\\nSimulation time = %f sec\',jacobi_iter,time_jacobi);
title(title_text1);


elseif method == 2

%GAUSS SEIDAL METHOD
tic;
[gs,gs_iter]= gs_steady(T,T_old,error_tol,nx,ny,dx,dy);
time_gs = toc;

%Plotting
figure(2)
[A,C] = contourf(x,fliplr(y),gs,\'ShowText\',\'on\');
xlabel(\'X axis\');
ylabel(\'Y axis\');
title_text2 = sprintf(\'2D Steady-state heat conduction-Gauss seidal\\nNo.of iterations to converge = %d\\nSimulation time = %f sec\',gs_iter,time_gs);
title(title_text2);
    
elseif method == 3
    
%SOR METHOD
tic;
[sor,sor_iter] = sor_steady(T,T_old,error_tol,nx,ny,dx,dy);
time_sor = toc;

figure(3)
[O,Q] = contourf(x,fliplr(y),sor,\'ShowText\',\'on\');
xlabel(\'X axis\');
ylabel(\'Y axis\');
title_text3 = sprintf(\'2D Steady-state heat conduction-SOR\\nNo.of iterations to converge = %d\\nSimulation time = %f sec\',sor_iter,time_sor);
title(title_text3);

end

MATLAB function for steady state Jacobi solver:

function [T_out_jac,jacobi_iter]  = jacobi_steady(T,T_old,error_tol,nx,ny,dx,dy)

jacobi_iter = 1;

error = 5e5;
    %Convergence loop
 while error > error_tol
   
     %Nodal loop
    for i = 2 : (nx-1)
         
        for j = 2 : (ny-1) 
            
            Z = (2*((dx^2)+(dy^2)))/((dx^2)*(dy^2));
            H = (T_old(i-1,j) + T_old(i+1,j))/(dx^2);
            V = (T_old(i,j-1) + T_old(i,j+1))/(dy^2);
            T(i,j) = (H + V)/Z;
            
        end
    end
    error = max(abs(T_old(i,j) - T(i,j)));
    T_old = T;
    jacobi_iter = jacobi_iter + 1;
    
 end

T_out_jac = T;

MATLAB function for steady state Gauss seidal solver:

function [T_out_gs, gs_iter] = gs_steady(T,T_old,error_tol,nx,ny,dx,dy)

gs_iter = 1;

error = 5e5;
    %Convergence loop
 while error > error_tol
   
     %Nodal loop
    for i = 2 : (nx-1)
         
        for j = 2 : (ny-1) 
            
            Z = (2*((dx^2)+(dy^2)))/((dx^2)*(dy^2));
            H = (T(i-1,j) + T_old(i+1,j))/(dx^2);
            V = (T(i,j-1) + T_old(i,j+1))/(dy^2);
            T(i,j) = (H + V)/Z;
            
        end
    end
    error = max(abs(T_old(i,j) - T(i,j)));
    T_old = T;
    gs_iter = gs_iter + 1;
    
 end

T_out_gs = T;

MATLAB function for steady state SOR solver:

function [T_out_sor, sor_iter] = sor_steady(T,T_old,error_tol,nx,ny,dx,dy)

w = 0.999; %relaxation parameter
sor_iter = 1;

error = 5e5;
    %Convergence loop
 while error > error_tol
   
     %Nodal loop
    for i = 2 : (nx-1)
         
        for j = 2 : (ny-1) 
            
            Z = (2*((dx^2)+(dy^2)))/((dx^2)*(dy^2));
            H = (T(i-1,j) + T_old(i+1,j))/(dx^2);
            V = (T(i,j-1) + T_old(i,j+1))/(dy^2);
            T(i,j) = ((H + V)/Z)*w;
            
        end
    end
    error = max(abs(T_old(i,j) - T(i,j)));
    T_old = T;
    sor_iter = sor_iter + 1;
    
 end

T_out_sor = T;

The plots depicting the temperature distribution using the three methods in steady state are as follows,

Inference from steady state heat conduction plots:

The temperature matrix at the end of the computation is same for all the three techniques which is necessary. So as their plots. The key results to be noted is that the time taken for computation and the number of iterations each technique takes to attain convergence.

The number of iterations for the solution to converge is in the following order,

Jacobi > Gauss seidal > SOR

So, the successive over relaxation method\'s solution has been converged in least number of iterations possible among the three. But it mainly depends on the relaxation factor specified. The relaxation factor depends on the spectral radius of the jacobi of the iterative matrix. In this case it is taken as 0.999.

Now the time taken for computation is in the order of,

Jacobi > SOR > Gauss seidal

So for the 2D heat conduction equation, in the steady state condition, the optimal iterative technique is Successive Over Relaxation method.

SOLVING THE 2D TRANSIENT HEAT CONDUCTION EQUATION:

IMPLICIT METHOD:

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

Applying numerical scheme to each of the terms in the equation to solve it numerically,

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

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

Substituting these terms in the original equation,

`((T_(i,j))^(n+1) - (T_(i,j))^(n))/(Deltat) + alpha((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))^(n+1) = 0`

`(T_(i,j))^(n+1) = (T_(i,j))^(n) + alphaDeltat((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))^(n+1)`

Now lets take,

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

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

The terms k1 and k2 are called the Courant numbers. It is also responsible for the stability of the solution obtained.There can be three courant numbers in an PDE for three dimensions.

So by substituting the k1 and k2 and rearranging the form, we get,

`(T_(i,j))^(n+1) = ((T_(i,j))^(n) + k_1(T_(i-1,j) + T_(i+1,j))^(n+1)+ k_2(T_(i,j-1) + T_(i,j+1))^(n+1))/(1+2k_1 + 2k_2`

Solving above equation iteratively in MATLAB will give the results.

MATLAB main function for 2D transient -Implicit method:

clear all
close all
clc

%Dimensions of the domain
L = 1;
B = 1;

%No. of grid points
nx = 20;
ny = 20;

%Initiating T matrix
T = ones(nx,ny);
x = linspace(0,L,nx);
y = linspace(0,B,ny);

%BC
T(:,1) = 400;
T(:,ny) = 800;
T(1,:) = 600;
T(nx,:) = 900;

%Refining BC\'s at the corners
T(1,1) = (T(1,2) + T(2,1))/2;
T(1,end) = (T(1,end-1) + T(2,end))/2;
T(end,1) = (T(end-1,1) + T(2,end))/2;
T(end,end) = (T(end-1,end) + T(end,end-1))/2;

T_old = T;
T_prev_dt = T;

dx = x(2) - x(1);
dy = y(2) - y(1);

%For CFL condition, Courant number < 1
CFL_no = 0.1;
alpha = 1.11e-4; %Thermal diffusivity for copper(m2/s)
dt = (CFL_no*((dx)^2))/alpha;
nt = 1000;

%Courant numbers
a = (alpha*dt)/(dx)^2;
b = (alpha*dt)/(dy)^2;

error_tol = 1e-4;

term1 = 1 + (2*a) + (2*b);
term2 = a/term1;
term3 = b/term1;

%Specifying the iterative method

sprintf(\"1.Jacobi\\n2.Gauss seidal\\n3.SOR\")

method = input(\'Method = \');

if method == 1
      
%JACOBI METHOD    
tic;
[jacobi, jacobi_iter]  = jacobi_transient(T,T_old,T_prev_dt,term1,term2,term3,error_tol,nt,nx,ny);
time_jacobi = toc;

%Plotting
figure(1)
[X,Y] = contourf(x,fliplr(y),jacobi,\'ShowText\',\'on\');
xlabel(\'X axis\');
ylabel(\'Y axis\');
title_text1 = sprintf(\'2D Transient heat conduction-Jacobi(Implicit)\\nNo.of iterations to converge = %d\\nSimulation time = %f sec\',jacobi_iter,time_jacobi);
title(title_text1);


elseif method == 2

%GAUSS SEIDAL METHOD
tic;
[gs,gs_iter]= gs_transient(T,T_old,T_prev_dt,term1,term2,term3,error_tol,nt,nx,ny);
time_gs = toc;

%Plotting
figure(2)
[A,C] = contourf(x,fliplr(y),gs,\'ShowText\',\'on\');
xlabel(\'X axis\');
ylabel(\'Y axis\');
title_text1 = sprintf(\'2D Transient heat conduction-Gauss seidal(Implicit)\\nNo.of iterations to converge = %d\\nSimulation time = %f sec\',gs_iter,time_gs);
title(title_text1);
    
elseif method == 3
    
%SOR METHOD
tic;
[sor,sor_iter] = sor_transient(T,T_old,T_prev_dt,term1,term2,term3,error_tol,nt,nx,ny);
time_sor = toc;

figure(3)
[O,Q] = contourf(x,fliplr(y),sor,\'ShowText\',\'on\');
xlabel(\'X axis\');
ylabel(\'Y axis\');
title_text1 = sprintf(\'2D Transient heat conduction-SOR(Implicit)\\nNo.of iterations to converge = %d\\nSimulation time = %f sec\',sor_iter,time_sor);
title(title_text1);

end

The material considered here is Copper. So the thermal diffusivity of copper is 1.11e-4 m2/s. Also the Courant number for both dimensions are taken as 0.1 to satisfy the CFL condition.

MATLAB function for transient Jacobi solver-Implicit method:

function [T_out_jac,jacobi_iter]  = jacobi_transient(T,T_old,T_prev_dt,term1,term2,term3,error_tol,nt,nx,ny)

jacobi_iter = 1;

%Time loop
for k = 1:nt
error = 5e5;
    %Convergence loop
 while error > error_tol
   
     %Nodal loop
    for i = 2 : (nx-1)
         
        for j = 2 : (ny-1) 
            
            H = T_old(i-1,j) + T_old(i+1,j);
            V = T_old(i,j-1) + T_old(i,j+1);
            T(i,j) = (T_prev_dt(i,j)/term1) + (H*term2) + (V*term3);
            
        end
    end
    error = max(max(abs(T_old - T)));
    T_old = T;
    jacobi_iter = jacobi_iter + 1;
    
 end
T_prev_dt = T;
end

T_out_jac = T;

MATLAB function for transient Gauss seidal solver-Implicit method:

function [T_out_gs, gs_iter] = gs_transient(T,T_old,T_prev_dt,term1,term2,term3,error_tol,nt,nx,ny)

gs_iter = 1;

%Time loop
for k = 1:nt
error = 5e5;
    %Convergence loop
 while error > error_tol
   
     %Nodal loop
    for i = 2 : (nx-1)
         
        for j = 2 : (ny-1) 
            
            H = T(i-1,j) + T_old(i+1,j);
            V = T(i,j-1) + T_old(i,j+1);
            T(i,j) = (T_prev_dt(i,j)/term1) + (H*term2) + (V*term3);
            
        end
    end
    error = max(max(abs(T_old - T)));
    T_old = T;
    gs_iter = gs_iter + 1;
    
 end
T_prev_dt = T;
end

T_out_gs = T;

MATLAB function for transient SOR solver-Implicit method:

function [T_out_sor,sor_iter] = sor_transient(T,T_old,T_prev_dt,term1,term2,term3,error_tol,nt,nx,ny)

w = 0.999; %relaxation parameter
sor_iter = 1;

%Time loop
for k = 1:nt
error = 5e5;
    %Convergence loop
 while error > error_tol
   
     %Nodal loop
    for i = 2 : (nx-1)
         
        for j = 2 : (ny-1) 
            
            H = T(i-1,j) + T_old(i+1,j);
            V = T(i,j-1) + T_old(i,j+1);
            T(i,j) = T_prev_dt(i,j)*(1-w) + w*((T_prev_dt(i,j)/term1) + (H*term2) + (V*term3));
            
            
        end
    end
    error = max(max(abs(T_old - T)));
    T_old = T;
    sor_iter = sor_iter + 1;
    
 end
T_prev_dt = T;
end

T_out_sor = T;

The plots depicting the temperature distribution using the three methods in transient state are as follows,

Inference from transient state heat conduction plots-Implicit method:

The temperature distribution is exactly the same all the three types, but the number of iterations for convergence and the itime taken for it makes them distinguished. The number of iterations for the solution to be converged is in the order of,

Jacobi > Gauss seidal > SOR

The order is same as that of the steady state method,but the ratio of difference here is very low when compared to steady state. In steady state, the number of iterations in GS method is half the time of Jacobi method. But here, it\'s more than 3/4th. The time taken for convergence is in the order of

Jacobi > SOR> Gauss seidal

So, to solve the 2D transient heat conduction implicitly, the optimal method is Gauss-seidal, since Gauss-seidal method takes only 2 iterations more than SOR to converge and the time is also lessser.

EXPLICIT METHOD:

`(T_(i,j))^(n+1) = ((T_(i,j))^(n) + k_1(T_(i-1,j) + T_(i+1,j))^(n)+ k_2(T_(i,j-1) + T_(i,j+1))^(n))/(1+2k_1 + 2k_2`

The difference between the implicit and explicit method is that in explicit method, all the required temperature values for calculating the temperature at later state of the system are taken from the current state of the system. It can be seen in the above equation that all the terms in RHS are at \'n\' timestep and the LHS being \'n+1\'.

MATLAB code for 2D transient -Explicit method:

clear all
close all
clc

%Dimensions of the domain
L = 1;
B = 1;

%No. of grid points
nx = 20;
ny = 20;

%Initiating T matrix
T = ones(nx,ny);
x = linspace(0,L,nx);
y = linspace(0,B,ny);

%BC
T(:,1) = 400;
T(:,ny) = 800;
T(1,:) = 600;
T(nx,:) = 900;

%Refining BC\'s at the corners
T(1,1) = (T(1,2) + T(2,1))/2;
T(1,end) = (T(1,end-1) + T(2,end))/2;
T(end,1) = (T(end-1,1) + T(2,end))/2;
T(end,end) = (T(end-1,end) + T(end,end-1))/2;

T_old = T;
T_prev_dt = T;

dx = x(2) - x(1);
dy = y(2) - y(1);

%For CFL condition, Courant number < 1
CFL_no = 0.2;
alpha = 1.11e-4; %Thermal diffusivity for copper(m2/s)
dt = (CFL_no*((dx)^2))/alpha;
nt = 1000;

%Courant numbers
a = (alpha*dt)/(dx)^2;
b = (alpha*dt)/(dy)^2;

tic;
%Time loop
for k = 1:nt
   
     %Nodal loop
    for i = 2 : (nx-1)
         
        for j = 2 : (ny-1) 
            
            H = T_old(i-1,j) - (2.*(T_old(i,j))) + T_old(i+1,j);
            V = T_old(i,j-1) - (2.*(T_old(i,j))) + T_old(i,j+1);
            T(i,j) = T_old(i,j) + (a*H) + (b*V);
            
        end
    end
    T_old = T;
    
end
time_jacobi = toc;

%Plotting
figure(1)
[X,Y] = contourf(x,fliplr(y),T,\'ShowText\',\'on\');
xlabel(\'X axis\');
ylabel(\'Y axis\');
title_text1 = sprintf(\'2D Transient heat conduction-Explicit\\nSimulation time = %f sec\\nCourant number = %f\',time_jacobi,CFL_no);
title(title_text1);

The plots depicting the temperature distribution using the explicit method for different courant number are as follows,

Inference from transient state heat conduction plot-Explicit method:

It can be clearly seen that the time taken for the computation in Explicit method is nearly 1/3rd of that of the implicit techniques. Also the stability of solution in explicit method depends mainly on the courant number. It should satisfy the CFL condition, which is C<1 for the solution to be stable. But that does not mean all the solutions with C<1 should be stable. The stability phenomenon depends on other parameters too. The timestep size should be as smaller as possible.

In this case, the solution with C=0.2 is stable. But raising it to C=0.3 makes it unstable. So the optimal Courant number for this case is 0.2.

CONCLUSION:

If the 2D heat conduction equation needs to be in steady state condition, then it is optimal to use SOR method. When it comes to transient condition, a number of things should be taken into consideration. The implicit method is computationally costly and also difficult to program. But the advantage is that it gives accurate results and less prone to unstable solutions. On the other hand, explicit method is cheap at computation and also easy to program. But the price is that it is prone to unstable solutions easily. It depends many parameters but mainly the courant number and timestep value should be chosen optimally.