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

solution for steady and transient state 2D heat conduction equation using matlab programming

SSaurabh Kumar Sharma14 Sept 2019

Steady state solution:-

We make the assumption as there is no heat convection, no internal heat generation, and no change of itemerature gradient with resect to time. therefore the generalised heat equation in 2D becomes as:-

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

Solving the steady state equation does not tinvolve time marching therefore the temperature gradient does not change with respect to the time.

We discritize the spatial term using central differencing:-

`(∂^2 T)/(∂x^2) =(T_(i-1,j)-2*T_(i,j)+T_(i+1,j))/(∆x^2)`

`(∂^2 T)/(∂y^2) =(T_(i,j-1)-2*T_(i,j)+T_(i,j+1))/(∆y^2)`

Therefore we get

`[(T_(i-1,j)-2*T_(i,j)+T_(i+1,j))/(∆x^2) ]+[(T_(i,j-1)-2*T_(i,j)+T_(i,j+1))/(∆y^2) ]=0`

Following is the main script for performing the steady 2D heat conduction:-

%solving the steady and unsteady 2D heat conduction equation
%dT/dt + alpha*(d^T/dx^2 + d^2T/dy^2)=0
%x=0 T=400K
%x=1 T=800K
%y=0 T=900K
%y=1 T=600K

clear all 
close all 
clc

%inputs for the iterative solvers
nx=20; ny=20;
T=298*ones(nx,ny);
T_L=400;
T_R=800;
T_T=600;
T_B=900;
L=1;
x=linspace(0,L,nx);
y=linspace(L,0,ny);
[X Y]=meshgrid(x,y);
dx=x(2)-x(1);
dy=y(2)-y(1);

%assigning the boundary condition 
T(2:ny-1,1)=T_L;    %left wall temp.
T(2:ny-1,nx)=T_R;   %right wall temp.
T(1,2:ny-1)=T_T;    %top wall temp.
T(nx,2:ny-1)=T_B;   %bottom wall temp.

%temperature at the corner can be taken by average temperature value
T(1,1)=(T_L+T_T)/2;
T(nx,ny)=(T_R+T_B)/2;
T(1,ny)=(T_T+T_R)/2;
T(nx,1)=(T_L+T_B)/2;

%creating a copy of T
T_old=T;

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

iterative_solver=input(\'enter the desired value\')
error=9e9;
tol=1e-4;

%jacobi method
if iterative_solver==1
    jacobi_iter=1
    while(error>tol)
       for i=2:(nx-1)
           for j=2:(ny-1)
             T(i,j)=k1*(T_old(i-1,j)+T_old(i+1,j))+k2*(T_old(i,j+1)+T_old(i,j-1));
           end  
       end
       error=max(max(abs(T_old-T)));
       T_old=T;
       jacobi_iter=jacobi_iter+1;
    end  
       
       figure(1)
       contourf(x,y,T);
       [a,b]=contourf(x,y,T);
       clabel(a,b);
       xlabel(\'X\')
       ylabel(\'Y\')
       colorbar
       colormap(jet)
       title_text=sprintf(\'solution of steady state 2-D conduction equation(jacobi iterative solver) \\n iteration number=%d\' ,jacobi_iter);
       title(title_text)
end


    %GS method
if iterative_solver==2
    gs_iter=1
    while(error>tol)
       for i=2:(nx-1)
           for j=2:(ny-1)
             T(i,j)=k2*(T_old(i+1,j)+T(i-1,j))+k1*(T_old(i,j+1)+T(i,j-1));
           end  
       end
       error=max(max(abs(T_old-T)));
       T_old=T;
       gs_iter=gs_iter+1;
    end
       figure(2)
       contourf(x,y,T);
       [a,b]=contourf(x,y,T);
       clabel(a,b);
       xlabel(\'X\')
       ylabel(\'Y\')
       colorbar
       colormap(jet)
       title_text=sprintf(\'solution of steady state 2-D conduction equation(gauss seidel iterative solver) \\n iteration number=%d\',gs_iter);
       title(title_text)
end


   %SOR method
if iterative_solver==3
    SOR_iter=1
    alpha=149e-6;
    while(error>tol)
       for i=2:(nx-1)
           for j=2:(ny-1)
             T(i,j)=(1-alpha)*T_old(i,j)+alpha*(k2*(T_old(i+1,j)+T(i-1,j))+k1*(T_old(i,j+1)+T(i,j-1)));
           end  
       end
       error=max(max(abs(T_old-T)));
       T_old=T;
       SOR_iter=SOR_iter+1;
    end
       figure(3)
       contourf(x,y,T);
       [a,b]=contourf(x,y,T);
       clabel(a,b);
       xlabel(\'X\')
       ylabel(\'Y\')
       colorbar
       colormap(jet)
       title_text=sprintf(\'solution of 2-D steady state conduction equation(SOR iterative solver) \\n iteration number=%d\',SOR_iter);
       title(title_text)
    
end

OUTPUT:-

1. Jacobi method:-

number of iteration =827

\"\"

2.Gauss seidel method:-

number of iteration=443

\"\"

3. SOR method:-

number of iteration=1235745

\"\"

Discussion(steady state):-

From the above results we can deduce that among the three iterative methods we get the faster results and with less error in the case of successive over relaxation method than the gauss seidel method which is faster than tha jacobi method.

Unsteady/Transient heat conduction:-

The generalised transient heat conduction equation in 2D can be written as:-

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

The above equation is generalised by assuming that there is no heat convection, no internal heat generaation or dissipation and alpha is constant.

We discritize spatial terms using central differencing:-

`(∂^2 T)/(∂x^2) =(T_(i-1,j)-2*T_(i,j)+T_(i+1,j))/(∆x^2)`

`(∂^2 T)/(∂y^2) =(T_(i,j-1)-2*T_(i,j)+T_(i,j+1))/(∆y^2)`

We dicritize temporal terms by forward differencing:-

`(∂T)/(∂t)=(T_(i,j)^(n+1)-T_(i,j)^n)/(∆t)`

hence

`(T_(i,j)^(n+1)-T_(i,j)^n)/(∆t)=alpha*[(T_(i-1,j)-2*T_(i,j)+T_(i+1,j))/(∆x^2)+(T_(i,j-1)-2*T_(i,j)+T_(i,j+1))/(∆y^2) ]`

We have defined some constants nammely k1 and k2

`k1=alpha*((/_\\t)/(/_\\x^2))`

`k2=alpha*((/_\\t)/(/_\\y^2))`

Therefore

`T_(i,j)^(n+1)=(T_(i,j)^n+k1*(T_(i-1,j)^(n+1)+T_(i+1,j)^(n+1) )+k2*(T_(i,j-1)^(n+1)+T_(i,j+1)^(n+1)))/(1+2*k1+2*k2)`

In transient condition temperature changes with time and reaches a steady state.The above equation is iterated by both methods explicit and implicit.

Explicit:-

In explicit method the temperature at next level is calculated by means of surrounding teperature values and increase the time steps.

Following is the code to solve the 2d transient heat conduction equation explictly:-

%solving the transient 2D heat conduction equation explicitly
%dT/dt + alpha*(d^2T/dx^2 + d^2T/dy^2)=0
%x=0 T=400K
%x=1 T=800K
%y=0 T=900K
%y=1 T=600K

clear all 
close all 
clc

%inputs for the iterative solvers
nx=20; ny=20;L=1;
T=298*ones(nx,ny);
x=linspace(0,L,nx);
y=linspace(L,0,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
dt=1e-2;
nt=500;
T_L=400;
T_R=800;
T_T=600;
T_B=900;

%assigning the boundary condition 
T(2:ny-1,1)=T_L;    %left wall temp.
T(2:ny-1,nx)=T_R;   %right wall temp.
T(1,2:ny-1)=T_T;    %top wall temp.
T(nx,2:ny-1)=T_B;   %bottom wall temp.

%temperature at the corner can be taken by average temperature value
T(1,1)=(T_L+T_T)/2;
T(nx,ny)=(T_R+T_B)/2;
T(1,ny)=(T_T+T_R)/2;
T(nx,1)=(T_L+T_B)/2;

%creating a copy of T
T_old=T;
T_prev_dt=T_old;


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

CFL=alpha*(dt/dx^2);


iterative_solver=input(\'enter the desired number=\')
error=9e9;
tol=1e-4;

%jacobi method
if iterative_solver==1
    jacobi_iter=1
   for k=1:nt 
    while(error>tol)
       for i=2:(nx-1)
           for j=2:(ny-1)
               term1=T_old(i,j);
               term2=k1*(T_old(i-1,j)+T_old(i+1,j)-2*(term1));
               term3=k2*(T_old(i,j+1)+T_old(i,j-1)-2*(term1));
               T(i,j)=term1+term2+term3;
           end  
       end
       error=max(max(abs(T_old-T)));
       T_old=T;
       jacobi_iter=jacobi_iter+1;
     end  
       T_prev_dt=T;
   end
       figure(1)
       contourf(x,y,T);
       [a b]=contourf(x,y,T);
       clabel(a,b);
       xlabel(\'X\')
       ylabel(\'Y\')
       colorbar
       colormap(jet)
       title_text=sprintf(\'solution of Transient 2-D conduction equation explicitly(jacobi iterative solver) \\n iteration number=%d\',jacobi_iter);
       title(title_text)

end


    %GS method
if iterative_solver==2
    gs_iter=1
   for k=1:nt 
    while(error>tol)
       for i=2:(nx-1)
           for j=2:(ny-1)
               term4=T_old(i,j);
               term5=k1*(T_old(i+1,j)+T(i-1,j)-2*(term4));
               term6=k2*(T_old(i,j+1)+T(i,j-1)-2*(term4));
             T(i,j)=term4+term5+term6;
           end  
       end
       error=max(max(abs(T_old-T)));
       T_old=T;
       gs_iter=gs_iter+1;
    end
    T_prev_dt=T;
   end
       figure(2)
       contourf(x,y,T);
       [a,b]=contourf(x,y,T);
       clabel(a,b);
       xlabel(\'X\')
       ylabel(\'Y\')
       colorbar
       colormap(jet)
       title_text=sprintf(\'solution of Transient 2-D conduction equation expicitly(Gauss seidel iterative solver) \\n iteration number=%d\',gs_iter);
       title(title_text)
end




   %SOR method
if iterative_solver==3
    SOR_iter=1;
    alpha=1.2;
   for k=1:nt 
    while(error>tol)
       for i=2:(nx-1)
           for j=2:(ny-1)
               term_1=T_old(i,j);
               term_2=k1*(T_old(i+1,j)+T(i-1,j)-2*(term_1));
               term_3=k2*(T_old(i,j+1)+T(i,j-1)-2*(term_1));
               T(i,j)=(1-alpha)*term_1+alpha*(term_1+term_2+term_3);
           end  
       end
       error=max(max(abs(T_old-T)));
       T_old=T;
       SOR_iter=SOR_iter+1;
    end
    T_prev_dt=T;
   end
       figure(3)
       contourf(x,y,T);
       [a,b]=contourf(x,y,T);
       clabel(a,b);
       xlabel(\'X\')
       ylabel(\'Y\')
       colorbar
       colormap(jet)
       title_text=sprintf(\'solution of Transient 2-D conduction equation explicitly(SOR iterative solver) \\n iteration number=%d\',SOR_iter);
       title(title_text)
end

OUTPUT:-

1.jacobi method:-

number of iteration:-10748

\"\"

2. Gauss seidel method:-

number of iteration=10479

\"\"

3. SOR method:-

number of iteration=8874

\"\"

Implicit method:-

Following is the code to solve the 2D transient heat conduction equation implicitly:-

%solving the transient 2D heat conduction equation implicitly
%dT/dt + alpha*(d^2T/dx^2 + d^2T/dy^2)=0
%x=0 T=400K
%x=1 T=800K
%y=0 T=900K
%y=1 T=600K

clear all 
close all 
clc

%inputs for the iterative solvers
nx=20; ny=20;L=1;
T=298*ones(nx,ny);
x=linspace(0,L,nx);
y=linspace(L,0,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
T_L=400;
T_R=800;
T_T=600;
T_B=900;

%assigning the boundary condition 
T(2:ny-1,1)=T_L;    %left wall temp.
T(2:ny-1,nx)=T_R;   %right wall temp.
T(1,2:ny-1)=T_T;    %top wall temp.
T(nx,2:ny-1)=T_B;   %bottom wall temp.

%temperature at the corner can be taken by average temperature value
T(1,1)=(T_L+T_T)/2;
T(nx,ny)=(T_R+T_B)/2;
T(1,ny)=(T_T+T_R)/2;
T(nx,1)=(T_L+T_B)/2;

%creating a copy of T
T_old=T;
T_prev_dt=T_old;

dt=300;
nt=50000;
alpha=0.04;
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);


iterative_solver=input(\'enter the desired number=\')
error=9e9;
tol=1e-4;

%jacobi method
if iterative_solver==1
    jacobi_iter=1
   for k=1:nt
    while(error>tol)
       for i=2:(nx-1)
           for j=2:(ny-1)
               term1=T_old(i,j);
               term2=k1*(T_old(i-1,j)+T_old(i+1,j));
               term3=k2*(T_old(i,j+1)+T_old(i,j-1));
               T(i,j)=(T_prev_dt(i,j)+term2+term3)/(1+2*k1+2*k2);
           end  
       end
       error=max(max(abs(T_old-T)));
       T_old=T;
       jacobi_iter=jacobi_iter+1;
     end  
       T_prev_dt=T;
   end
       figure(1)
       contourf(x,y,T);
       [a,b]=contourf(x,y,T);
       clabel(a,b);
       xlabel(\'X\')
       ylabel(\'Y\')
       colorbar
       colormap(jet)
       title_text=sprintf(\'solution of Transient 2-D conduction \\n equation implicitly(jacobi iterative solver) \\n jacobi iteration number=%d\',jacobi_iter);
       title(title_text)
end


    %GS method
if iterative_solver==2
    gs_iter=1
   for k=1:nt 
    while(error>tol)
       for i=2:(nx-1)
           for j=2:(ny-1)
               term4=T_old(i,j);
               term5=k1*(T_old(i+1,j)+T(i-1,j));
               term6=k2*(T_old(i,j+1)+T(i,j-1));
             T(i,j)=(T_prev_dt(i,j)+term5+term6)/(1+2*k1+2*k2);
           end  
       end
       error=max(max(abs(T_old-T)));
       T_old=T;
       gs_iter=gs_iter+1;
    end
    T_prev_dt=T;
   end
       figure(2)
       contourf(x,y,T);
       [a,b]=contourf(x,y,T);
       clabel(a,b);
       xlabel(\'X\')
       ylabel(\'Y\')
       colorbar
       colormap(jet)
       title_text=sprintf(\'solution of Transient 2-D conduction \\n equation implicitly(Gauss seidel iterative solver) \\n gauss iteration number=%d\',gs_iter);
       title(title_text)
end




   %SOR method
if iterative_solver==3
    SOR_iter=1
   for i=1:nt 
       error = 10;
    while(error>tol)
       for i=2:(nx-1)
           for j=2:(ny-1)
               term_1=T_old(i,j);
               term_2=k1*(T_old(i+1,j)+T(i-1,j));
               term_3=k2*(T_old(i,j+1)+T(i,j-1));
               T(i,j)=(1-alpha)*T_prev_dt(i,j)+alpha*((T_prev_dt(i,j)+term_2+term_3)/(1+2*k1+2*k2));
           end  
       end
       error=max(max(abs(T_old-T)));
       T_old=T;
       SOR_iter=SOR_iter+1;
    end
    T_prev_dt=T;
   end
       figure(3)
       contourf(x,y,T);
       [a,b]=contourf(x,y,T);
       clabel(a,b);
       xlabel(\'X\')
       ylabel(\'Y\')
       colorbar
       colormap(jet)
       title_text=sprintf(\'solution of Transient 2-D conduction \\n equation implicitly(SOR iterative solver) \\n SOR iteration number=%d\',SOR_iter);
       title(title_text)
end

OUTPUT:-

1. Jacobi method:-

number of iteration = 823

\"\"

2. Gauss seidel:-

number of iteration:-441

\"\"

3. SOR method:-

number of iterations=72530

\"\"

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.