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

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

 

 


Projects by Saurabh sharma

Challenge on Combustion
Saurabh sharma · 2020-02-16 18:45:07

  Challenge on Combustion Objective:- Perform a combustion simulation on the combuster model and plot the variation of the mass fraction of the different species in the simulation using line probes at different locations of the combuster. You need to plot Read more

 

Cyclone Separator Challenge

 

Objective:-

We have to perform an analysis on the cyclone seperator model for  different number of inectioned particles.

1. The four imperical models use Read more

  Gate Valve Parametric Study A gate valve, also known as a sluice valve, is a valve that opens by lifting a barrier(gate) out of the path of the fluid.Gate valves are used to shut off the flow of liquids rather than for flow regulation. Gate valve require Read more

Conjugate Heat Transfer Analysis on a graphics card. Description of a conjugate heat transfer analysis:- Conjugate heat transfer represents the combination of heat transfer in liquids and solids simulteneously. The mode of heat transfer in solids is conduction and con Read more

    Rayleigh Taylor Instability Challenge SOLUTION-1:- Models used for Rayleigh-Taylor instability:- Early-time growth Potential flow Buoyancy-drag balance Energy budget   1. Early-time growth:- It was noted in the taylor (1950) that t Read more

Exhaust Port Challenge
Saurabh sharma · 2019-12-19 08:14:22

              Exhaust Port Challenge SOLUTION1:- Why and Where a CHT analysis is used? Conjugate Heat Transfer analysis:- Conjugate heat transfer analysis is generally used when there is a teperature variation during heating and co Read more

Ahmed Body Challenge
Saurabh sharma · 2019-12-17 16:08:31

              Ahmed Body Challenge  SOLUTION1: Ahmed body:- The Ahmed body was originally described By S.R. Ahmed in 1984. Ahmed body is a simplified car body. The shape provides a  model to study geometric effects Read more

    Simulation of Flow through a pipe in                                OpenFoam Objective of the project:- 1. Calculation of flow variables. 2. Generation of the blockMeshDict Read more

Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method Introduction:- nozzles are most commonly used in every fields in order to study the liquids and air properties, because those two elements plays an important role in the enginnering fields. H Read more

             FVM:: Interpolation and Gradient Schemes literature review FVM-Flux limiters and interpolation schemes What is FVM? FVM is a finite volume method for representing and evaluating PDE\'s in the form of algebraic equations Read more


Loading...

The End