BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

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

BBhavin Patel16 Jun 2021

Aim-To get the solution of the 2d heat conduction equation.

2d heat conduction equation of Cartesian co-ordinates.

from the figure we can see that we want to compute the value of Tp and our differential equation is elliptic type so we have to take information from four points around the point Tp

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

now for steady state

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

now by numerical descritization we can get

Tp=[(Tl+Tr)/dx^2]/K + [(Tt+Tb)/dy^2]/K

now for jacobi method we use old values for the Tl and Tb which we get from the previous equation.

for the gauss seidal method we use latest equation for Tl and Tb

code for the steady state

Step 1-clearing the previous variables and previous plots

Step2-defining the problem using the comments

Step3-defining the range of x and y and creating the mesh for that values

Step4-creating the array temperature

Step5- defining the boundary conditions

Step6-updating the boundary conditions in the temperature array

Step7-defining the iterative solver variable for executing different iterative method

Step8-section for the Jacobi method

Step9-tolerance and error values are defined for the while loop

Step10-two for loops are defined for updating the i and j values.

Step11-calculating error after the iterations and updating the values of T

Step12- ploting the graph using counterf

Step13- y axis is inverted with using the command reverse

clear all
close all
clc

%% solving the 2d heat conduction equation steady state

%equation d^2T/dx^2 + d^2T/dy^2=0
%inputs
nx=10;
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
[xx,yy]=meshgrid(x,y);

T=300*ones(nx,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
K=2*((dx^2 +dy^2)/(dx^2)*dy^2);

%applying the boundary condition
T(1,1)=500;
T(2:9,1)=400;
T(1,2:9)=600;
T(2:9,10)=800;
T(10,2:9)=900;
T(1,10)=700;
T(10,1)=650;
T(10,10)=850;

Told=T;

iterative_solver=3;

%% jacobi method

tic
if iterative_solver==1
    
    tolerance=0.01;
    error=1e9;
    jacobi_iter=1;
    while(error>tolerance)
        
        
        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(abs(Told(i,j)-T(i,j)));
        Told=T;
        jacobi_iter=jacobi_iter+1;
    end
     elapsed_time= toc;
    
   % surf(x,y,T)
       [t,h]=contourf(xx,yy,T)
       clabel(t,h);
       colorbar;
       colormap(jet)
       set(gca,'ydir','reverse')
       text1=sprintf('numbeer of jacobi iterations %d',jacobi_iter)
       text2=sprintf('elapsed time=%f',elapsed_time)
       text3=sprintf('steady state')
       title({text3;text1;text2})
        
       xlabel('x')
       ylabel('y')
       zlabel('T')
       pause(0.03)
   
end


%% gauss_seidal
tic
if iterative_solver==2
    %gauss_seidal
    tolerance=0.01;
    error=1e9;
    gs_iter=1;
    while(error>tolerance)
        
        
        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(abs(Told(i,j)-T(i,j)));
        gs_iter=gs_iter+1;
        Told=T;
    end
    elapsed_time=toc; 
    %surf(x,y,T)
    [t,h]=contourf(xx,yy,T)
    clabel(t,h);
    colorbar;
    colormap(jet)
    set(gca,'ydir','reverse')
    xlabel('x')
    ylabel('y')
    zlabel('T')
    text1=sprintf('numbeer of gauss iterations %d',gs_iter)
       text2=sprintf('elapsed time=%f',elapsed_time)
       text3=sprintf('steady state')
       title({text3;text1;text2})
       
    pause(0.03)
    
end

%% sor method
tic
if iterative_solver==3
    %sor method
    omega=1.1;
    tolerance=0.01;
    error=1e9;
    SOR_iter=1;
    alpha=2/(1+sin(pi*dx));
    while(error>tolerance)
        
        
        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*dx^2;
                T(i,j)=(1-omega)*Told(i,j)+omega*(term1+term2);
            end
        end
        
        error=max(abs(Told(i,j)-T(i,j)));
        Told=T;
        SOR_iter=SOR_iter+1;
    end
    elapsed_time=toc;
   
    %surf(x,y,T)
    [t,h]=contourf(xx,yy,T)
    clabel(t,h);
    set(gca,'ydir','reverse')
    xlabel('x')
    ylabel('y')
    zlabel('T')
    pause(0.03)
    colorbar;
    colormap(jet)
    text1=sprintf('numbeer of sor iterations %d',SOR_iter)
       text2=sprintf('elapsed time=%f',elapsed_time)
       text3=sprintf('steady state')
       title({text3;text1;text2})
       
end

output for steady state

steady state
method number of iteration
for convergence		 elapsed time
jacobi 100 0.011843
gauss seidal 54 0.010997
sor 45 0.010998

from the table we can see that the convergence rate of sor method is more than the jacobi method,and the computational time for the sor method is least.

unsteady state heat conduction in 2d

the governing equation is as below

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

now by the numerical descritization

for implicitly

term1=1/(1+2*K1+2*K2)
term2=K1*term1
term3=K2*term1
H=TL+ToldR
V=TT+TB
T=(T_initial*term1)+(H*term2)+(V*term2)

now for the gauss seidal we will use latest values for TL and TB computed at iteration

code for the solving the heat conduction implicitly

clear all
close all
clc

%% solving the 2d heat conduction equation unsteady state implicit method 
%equation d^2T/dx^2 + d^2T/dy^2=(1/alpha)*dT/dt
%inputs
nx=10;
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
[xx,yy]=meshgrid(x,y);

T=300*ones(nx,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);


%applying the boundary condition
T(1,1)=500;
T(2:9,1)=400;
T(1,2:9)=600;
T(2:9,10)=800;
T(10,2:9)=900;
T(1,10)=700;
T(10,1)=650;
T(10,10)=850;

Told=T;
T_initial=T;
%time inputs
dt=1e-2;

thermal_diffutivity=1.5;

K1=thermal_diffutivity*(dt/dx^2);
K2=thermal_diffutivity*(dt/dy^2);


iterative_solver=3;



 %% jacobi method
 tic
if iterative_solver==1
   
    tolerance=0.01;
    
    jacobi_iter=1;
    for n=1:1200
        error=1e9;
        while(error>tolerance)
            
            
            
            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)=(T_initial(i,j)*term1)+(H*term2)+(V*term2);
                    
                    
                end
            end
            error=max(max(abs(Told-T)));
            Told=T;
            jacobi_iter=jacobi_iter+1;
            
        end
        T_initial=T;
    end
    

elapsed_time=toc;
figure(1)
%surf(xx,yy,T)
[t,h]=contourf(xx,yy,T);
clabel(t,h);
colorbar;
colormap(jet)
set(gca,'YDIR','reverse')
text1=sprintf('number of jacobi iteration= %d',jacobi_iter);
text2=sprintf('elapsed time= %f',elapsed_time);
text3=sprintf('unsteady state ');
title({text1,text2,text3})
xlabel('x')
ylabel('y')
zlabel('T')
end
%% gauss_seidal method

if iterative_solver==2
    tic
        tolerance=0.01;
    gs_iter=1;
    for n=1:1200
        error=1e9;
        while(error>tolerance)
            
            
            
            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)=(T_initial(i,j)*term1)+(H*term2)+(V*term2);
                    
                    
                end
            end
            error=max(max(abs(Told-T)));
            Told=T;
            gs_iter=gs_iter+1;
            
        end
        T_initial=T;
    end
    

elapsed_time=toc;
%surf(xx,yy,T)
figure(2)
[t,h]=contourf(xx,yy,T)
clabel(t,h);
colorbar;
colormap(jet)
set(gca,'YDIR','reverse')
text1=sprintf('number of gauss seidal iteration= %d',gs_iter)
text2=sprintf('elapsed time= %f',elapsed_time)
text3=sprintf('unsteady state ')
title({text3,text1,text2})


xlabel('x')
ylabel('y')
zlabel('T')
pause(0.03)
end
%% sor method
tic
if iterative_solver==3
        tolerance=0.01;
        omega=1.1;
     
    
    sor_iter=1;
    for n=1:1200
        error=1e9;
        while(error>tolerance)
            
            
            
            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*((T_initial(i,j)*term1)+(H*term2)+(V*term2));
                    
                    
                end
            end
            error=max(max(abs(Told-T)));
            Told=T;
            sor_iter=sor_iter+1;
            
        end
        T_initial=T;
    end
    
elapsed_time=toc;
%surf(xx,yy,T)
[t,h]=contourf(xx,yy,T);
clabel(t,h);
colorbar;
colormap(jet)
set(gca,'YDIR','reverse')
text1=sprintf('number of sor iteration= %d',sor_iter);
text2=sprintf('elapsed time= %f',elapsed_time);
text3=sprintf('unsteady state ');
title({text1,text2,text3})


xlabel('x')
ylabel('y')
zlabel('T')
pause(0.03)
end

output of the code

unsteadystate(implicitly)
method number of iteration
for convergence		 elapsed time
jacobi 1741 0.0811036
gauss seidal 1549 0.073998
sor 1499 0.069868

here we are seeing the same result as steady state that the convergence rate of the sor method is high so the computational time is also low for the sor.

unsteady state 2d heat conduction using the explicit method

main code

clear all
close all
clc

%% unsteady state heat conduction equation by explicit method
%equation d^2T/dx^2 + d^2T/dy^2=(1/alpha)*dT/dt
%inputs
nx=10;
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
[xx,yy]=meshgrid(x,y);

T=300*ones(nx,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);


%applying the boundary condition
T(1,1)=500;
T(2:9,1)=400;
T(1,2:9)=600;
T(2:9,10)=800;
T(10,2:9)=900;
T(1,10)=700;
T(10,1)=650;
T(10,10)=850;

Told=T;
T_initial=T;
%time inputs
dt=1e-2;

thermal_diffutivity=1.5;

K1=thermal_diffutivity*(dt/dx^2);
K2=thermal_diffutivity*(dt/dy^2);


iterative_solver=3;
%% jacobi method
tic
if iterative_solver==1
    
    
    
    
    for n=1:120000
        
        
        
        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)=(T_initial(i,j)*term1)+(H*term2)+(V*term2);
                
                
            end
        end
        Told=T;
        
        
        T_initial=T;
    end
    
    
    elapsed_time=toc;
    figure(1)
    %surf(xx,yy,T)
    [t,h]=contourf(xx,yy,T);
    clabel(t,h);
    colorbar;
    colormap(jet)
    set(gca,'YDIR','reverse')
    
    text2=sprintf('elapsed time= %f',elapsed_time);
    text3=sprintf('unsteady state explicitly jacobi');
    title({text2,text3})
    xlabel('x')
    ylabel('y')
    zlabel('T')
    
end

%% gauss seidal method
if iterative_solver==2
    tic
    for n=1:1200
        
        
        
        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)=(T_initial(i,j)*term1)+(H*term2)+(V*term2);
                
                
            end
        end
        Told=T;
        
        
        T_initial=T;
    end
    
    
    elapsed_time=toc;
    %surf(xx,yy,T)
    figure(2)
    [t,h]=contourf(xx,yy,T)
    clabel(t,h);
    colorbar;
    colormap(jet)
    set(gca,'YDIR','reverse')
    text2=sprintf('elapsed time= %f',elapsed_time)
    text3=sprintf('unsteady state explicitly gauss seidal ')
    title({text3,text2})
    
    
    xlabel('x')
    ylabel('y')
    zlabel('T')
    pause(0.03)
end
%% sor method
tic
if iterative_solver==3
    omega=1.1;
    
    
    for n=1:1200
        
        
        
        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*((T_initial(i,j)*term1)+(H*term2)+(V*term2));
                
                
            end
        end
        Told=T;
        
        
        T_initial=T;
    end
    
    elapsed_time=toc;
    %surf(xx,yy,T)
    [t,h]=contourf(xx,yy,T);
    clabel(t,h);
    colorbar;
    colormap(jet)
    set(gca,'YDIR','reverse')
    
    text2=sprintf('elapsed time= %f',elapsed_time);
    text3=sprintf('unsteady state explicitly by sor ');
    title({text2,text3})
    
    
    xlabel('x')
    ylabel('y')
    zlabel('T')
    pause(0.03)
end

output of the code

Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.