2D Steady state Unsteady state heat conduction equation solve using implicit and explicit methods

Steady state & unsteady state 2D - heat conduction equation solve using implicit & explicit method

Diffusion – It’s the process of moving from higher concentration to lower concentration.

Heat diffusion – If the temperature increases the diffusion process will starts, Then the heat diffuses from high temperature region to low temperature region.

2D -Steady state heat conduction equation

del^2T/delx^2 +del^2T /dely^2 = 0

Discretize the above equation using taylor table method

2D – unsteady heat conduction equation.

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

The above equation is discretized by using taylors table method.

Solving the Steady state & unsteady state by implicit method

I.Gauss jacobi method

In this method, assume the guess value from previous level iteration to compute the all the unknowns.

II.Gauss seidel method

In this method we will use the compute the first unknown as same as jacobi method, but for the upcoming unknown we have used the updated value.

III.Sucessive over relaxation

It can used on the jacobi method (or) Gauss seidel method.It is used to enhance the iterative solver,which is used to enhance the performance of jacobi/gauss seidel method.

Matlab code for Unsteady state & steady state:

%Solving 2D unsteady state heat conduction equation using implicit method
%========================================================================

%du/dt = alpha*(d^T/dx^2 + d^T/dy^2)

clear all
close all
clc


%length 
L= 1;

%Defining square domain
nx=10;
ny=nx;

%Grids and nodes
%================
%Grids & nodes in X-axis 
x=linspace(0,L,nx);
dx=x(2)-x(1);

%Grids & nodes in Y-axis
y=linspace(0,L,nx);
dy=y(2)-y(1);

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

%nx X ny Matrix with boundary condition
T=ones(nx,ny);
T(:,1)=400;
T(1,:)=600;
T(:,end)=800;
T(end,:)=900;

dt = 0.08;%1e-4;

%Previous temperature
Told = T;

%Assumed thermal diffusivity
alpha = input(\'Enter the value of alpha=\');

%k1 & k2
k1 =alpha*(dt/dx^2);
k2 = alpha*(dt/dy^2);

%Assume Error value  and tolerance value% 
error = 9e9;
tol = 1e-4;
n=5;

%Jacobi method
%=============
iteration_no = input(\'Enter the iteration no = \');

%%Time loop
for nt = 1:n
    if iteration_no == 1
        jacobi_iter = 1;
      while (tol < error)  
            for i=2:nx-1
                for j=2:ny-1
                    Tl=Told(i-1,j);  %Temperature on left side
                    Tr=Told(i+1,j);  %Temperature on Right side
                    Tt=Told(i,j+1);  %Temperatuer on Top side 
                    Tb=Told(i,j-1);  %Temperature on Bottom side
                    term1 = 1+2*k1+2*k2;
                    
                    T(i,j) = (Told(i,j)+(k1*(Tr+Tl))+(k2*(Tt+Tb)))/term1;
                end
            end
                    error_n = max(abs(Told - T));
                    error = max(error_n);
                    Told=T;
                    
                    colormap(jet)
                    [C,h] =contourf(x,y,T);
                    clabel(C,h,\'FontSize\',8);
                    jacobi_iter = jacobi_iter+1;
                    title_text = sprintf(\'Jacobi Iteration Number = %d\',jacobi_iter);
                    title(title_text)
                    colorbar
                    pause(0.03)
                                     
      end
    end
                  %updating time step
                  dt_old = dt; 
end
 

%%Gauss Seidel method
%====================
 
for nt = 1:n
    if iteration_no == 2
        GS_iter = 1;
      while (tol < error)  
            for i=2:nx-1
                for j=2:ny-1
                    Tl=T(i-1,j);
                    Tr=Told(i+1,j);
                    Tt=T(i,j+1);
                    Tb=Told(i,j-1);
                    term1 = 1+2*k1+2*k2;
                    
                    T(i,j) = (Told(i,j)+(k1*(Tr+Tl))+(k2*(Tt+Tb)))/term1;
                end
            end
                    error_n = max(abs(Told - T));
                    error = max(error_n);
                    Told=T;
                    
                    colormap(jet)
                    [C,h] = contourf(X,Y,T);
                    clabel(C,h,\'FontSize\',8);
                    GS_iter = GS_iter+1;
                    title_text = sprintf(\'GS Iteration Number = %d\',GS_iter);
                    title(title_text)
                    colorbar
                    pause(0.03)    
      end
    end
                  %updating time step
                  dt_old = dt; 
end
 

%Sucessive Over Relaxation method
%================================

for nt = 1:n
    if iteration_no == 3
        SOR_iter = 1;
        while (tol < error)  
            for i=2:nx-1
                for j=2:ny-1
                    Tl=T(i-1,j);
                    Tr=Told(i+1,j);
                    Tt=T(i,j+1);
                    Tb=Told(i,j-1);
                    term1 = 1+2*k1+2*k2;
                    
                    T(i,j) = (Told(i,j)*(1-alpha))+alpha*((Told(i,j)+(k1*(Tr+Tl))+(k2*(Tt+Tb)))/term1);
                end
            end
                    error_n = max(abs(Told - T));
                    error = max(error_n);
                    Told=T;
                    
                    colormap(jet)
                    [C,h] =contourf(x,y,T);
                    clabel(C,h,\'FontSize\',8)
                    SOR_iter = SOR_iter+1;
                    title_text = sprintf(\'SOR Iteration Number = %d\',SOR_iter);
                    title(title_text)
                    colorbar
                    pause(0.03)
                                     
      end
    end
                  %updating time step
                  dt_old = dt; 
end

Results of Jacobi method:

Result of Gauss- seidel method

Result of SOR method(Over relaxed condition)

Result of SOR method(Under relaxed method)

Matlab code for unsteady state heat conduction using explicit method:

%Unsteady state 2D heat conduction equation using explicit method
%================================================================


clear all
close all
clc


%length 
L= 1;

%Defining square domain
nx=10;
ny=nx;

n=200;
%Grids and nodes
%---------------
%Grids & nodes in X-axis 
x=linspace(0,L,nx);
dx=x(2)-x(1);

%Grids & nodes in Y-axis
y=linspace(0,L,nx);
dy=y(2)-y(1);


%nx X ny Matrix with boundary condition
T=ones(nx,ny);
T(:,1)=400;
T(1,:)=600;
T(:,end)=800;
T(end,:)=900;


dt = 0.01;

%Previous temperature
Told = T;

%Assumed thermal diffusivity
alpha =0.25;


k1 =alpha*dt/(dx^2);
k2 = alpha*dt/(dy^2);

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


%Assume Error value  and tolerance value% 

k = 1;
%%Time loop
for nt = 1:n  
            for i=2:nx-1
                for j=2:ny-1
                    term1 = Told(i,j);
                    term2 = k1*(T(i-1,j)-2*Told(i,j)+Told(i+1,j));
                    term3 = k2*(T(i,j-1)-2*Told(i,j)+Told(i,j+1));
                    
                    T(i,j) = term1+term2+term3;
                end
            end
                   
                     Told=T;
                    
                     colormap(jet);
                     [C,h]=contourf(x,y,T);
                     clabel(C,h);
                     colorbar
                     title(sprintf(\'No.of iteration =%d\',k))
                     k=k+1;
                     pause(0.03)
                     
                     
                 
      end

Result for explicit method:

Matlab code for steady state heat conduction using implicit method:

%Solving steady state equation using implicit methods
%====================================================

clear all;close all;clc


%Domain length
L=1;

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

% grid spacing
x= linspace(0,L,nx);
dx = x(2)-x(1);

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


%Boundary condition
T=ones(nx,ny);

T(:,1)=400;
T(1,:)=600;
T(:,end)=800;
T(end,:)=900;

Told=T;

%Absolute error
error = 9e9;
tol = 1e-4;

iteration_method =input(\'Enter the iteration method =\'); 

%%Jacobo oteration
%------------------
if iteration_method ==1;
jacobi_iter = 1;
while (tol<error)
    for i=2:nx-1
        for j=2:ny-1
             T(i,j) = 0.25*(Told(i-1,j)+Told(i+1,j)+Told(i,j+1)+Told(i,j-1));
        end
    end
    
    error_n= max(abs(Told - T));
    error = max(error_n);
    Told = T;
    colormap(jet);
    [C,h]=contourf(x,y,T);
    clabel(C,h,\'Fontsize\',8);
    colorbar;
    jacobi_iter =jacobi_iter +1;
    title_text = sprintf(\'Jacobi iteration number = %d\',jacobi_iter);
    title(title_text);
    pause(0.03);
    
end
end
%plot(x,T)

%Gauss seidel method
%-------------------
if iteration_method ==2;
GS_iter = 1;
while (tol<error)
    for i=2:nx-1
        for j=2:ny-1
             T(i,j) = 0.25*(Told(i-1,j)+T(i+1,j)+Told(i,j+1)+T(i,j-1));
        end
    end
    
    error_n= max(abs(Told - T));
    error = max(error_n);
    Told = T;
    colormap(jet);
    [C,h]=contourf(x,y,T);
    clabel(C,h,\'Fontsize\',8);
    colorbar;
    GS_iter =GS_iter +1;
    title_text = sprintf(\'GS iteration number = %d\',GS_iter);
    title(title_text);
    pause(0.03);
    
end
end

%SOR method 
%-----------
if iteration_method ==3;
SOR_iter = 1;
alpha =input(\'Enter the value of alpha=\');
while (tol<error)
    for i=2:nx-1
        for j=2:ny-1
             T(i,j) = Told(i,j)*(1-alpha)+ (alpha*0.25*(Told(i-1,j)+T(i+1,j)+Told(i,j+1)+T(i,j-1)));
        end
    end
    
    error_n= max(abs(Told - T));
    error = max(error_n);
    Told = T;
    colormap(jet);
    [C,h]=contourf(x,y,T);
    clabel(C,h,\'Fontsize\',8);
    colorbar;
    SOR_iter =SOR_iter +1;
    title_text = sprintf(\'SOR iteration number = %d\',SOR_iter);
    title(title_text);
    pause(0.03);
    
end
end

Results for Jacobi method:

Result for Gauss seidel method:

Result for SOR method:

Observation:

1.Implicit method

From both code and plot,it is clearly observed that the number of iteration is getting reduced with different iterative methods used in both unsteady state and steady state.

And if timestep value decreases the number of iteration will increases and vice versa.The iteration process will go on until get the converged result.

2.Explicit method

For explicit method we are using the old level value,it doesnt uses any updated value.

The result looks as same as implicit method but it gives a limited no of iteration what we gave as a time level.

Then with the increase of dt time step,solution gets moving gradually from right to left.

If the dt value below 0.01 the solution is converging with given time level.