Solving Steady and transient 2D heat conduction equation in matlab

OBJECTIVE:

1.To solve the steady and transient(implicit and explicit) 2D heat conduction equation using jacobi,gauss- seidal and successive over relaxation iterative techiques in matlab.

2.To calculate and compare the simulation time and number of iterations taken to converge for each transient solver.

3. To determine the effect of time step(dt), courant number for each transient solver.

EQUATIONS:

1.Steady 2D heat conduction:

Since, the governing equation is elliptic in nature, central differencing scheme is taken to derive the difference equation.

substituting Central difference scheme in equation,

assuming number of grid points in x and y direction are same,therefore dx=dy,

2.Transient 2D heat conduction:

where, `alpha` is the thermal diffudivity(m2/s)

Thermal diffusivity[m2/s]: Thermal diffusivity is the thermal conductivity divided by density and specific heat capacity at const pressure.In a sense, Thermal diffusivity is a measure of thermal inertia in a material.

we have taken forward time marchining for time derivative and central difference scheme for spatial derivative,

For explicit transient 2D heat conduction:

For implicit transient 2D heat conduction:

by solving the above equation we get,

where,

We have solved Transient 2D heat conduction implicitly by using three types of iterative solvers,namely

1. Jacobi iteration method: Jacobi method is an iterative algorithm for determining the solutions of a strictly diagonally dominant system of linear equations. Each diagonal element is solved for, and an approximate value is plugged in. The process is then iterated until it converges.In jacobi method, values obtained at an iteration (values at i-1 or j-1 grid points) are not updated in further calculations of the iteration.

2. Gauss seidal iteration method: Gauss seidal iteration method is similar to jacobi iteration method, unlike jacobi iterative method, the values obtained at an iteration (values at i-1 or j-1 grid points) are updated in further calculations of the iteration.

3. Successive over relaxation method(sor): Successive over relaxation method is a variant of gauss-seidal method,where relaxation factor is used for usually faster convergence(if relaxation factor is greater than 1).

CFL Condition:

The CFL condition for heat conduction equation is

here, C is known as the courant number and its maximum value is 0.5 in order to get a stable solution.The CFL condition is valid for transient explicit condition.

CASE 1: Solution of steady heat conduction equation using jacobi,gauss-seidal and SOR gauss seidal iterative solvers.

Assumptions taken:

1. Number of grid points in x and y direction, nx=ny=10 also since domain length is same, dx=dy.

2. Boundary conditions: Temperature at left : 400K

Temperature at top : 600K

Temperature at right: 800K

Temperature at bottom: 900K

And Average termeprature is taken at edges of the domain.

3. Absolute error criterion,tolerance= 1e-4

4. SOR relaxation factor, omega= 1.25

PROGRAM1:

%steady heat conduction
%d^2T/dx^2 + d^2 T/dy^2=0
clear all
close all
clc
global nx ny omega tolerance %declaring global varaibles 
nx=10;   %no of grid points in x direction
ny=nx;   %no of grid points in y direction

x=linspace(0,1,nx);   %grid generation
y= linspace(0,1,ny);
omega=1.25;           %relaxation factor

%applying boundary conditions and assuming temperature at remaining domain
%to be 1
t= ones(nx,ny);
t(1,1)=0.5*(900+400);          %average temperature at edges
t(1,nx)=0.5*(900+800);
t(ny,1)=0.5*(400+600);
t(nx,ny)=0.5*(600+800);
t(2:ny-1,1)=400;
t(1,2:nx-1)=900;
t(2:ny-1,nx)=800;
t(ny,2:ny-1)=600;

told=t;                     
tolerance=1e-4;             %absolute error criterion


for iterative_sol=1:3
%jacobi
if iterative_sol==1
   [t_jac,sim_time_jac,jacobi_iter]=steady_jac(t,told) ;
figure(1)
contourf(x,y,t_jac,\'showtext\',\'on\')
hold on
colorbar
xlabel(\'X-Direction\')
ylabel(\'Y-Direction\')
title({\'Steady heat conduction-Jacobi\',[\'iteration no=\',num2str(jacobi_iter)],[\'simulation time=\',num2str(sim_time_jac)]})

  

end


%gauss-Seidal method
if iterative_sol==2
 [t_gs,sim_time_gs,gs_iter]=steady_gs(t,told) ;   


figure(2)
contourf(x,y,t_gs,\'showtext\',\'on\')
hold on
colorbar
xlabel(\'X-Direction\')
ylabel(\'Y-Direction\')
title({\'Steady heat conduction-gauss seidal\',[\'iteration no=\',num2str(gs_iter)],[\'simulation time=\',num2str(sim_time_gs)]})
  

end


%SOR method
if iterative_sol==3
 [t_sor,sim_time_sor,sor_iter]=steady_sor(t,told) ;  

figure(3)
contourf(x,y,t_sor,\'showtext\',\'on\')
hold on
colorbar
xlabel(\'X-Direction\')
ylabel(\'Y-Direction\')
title({\'Steady heat conduction-SOR\',[\'iteration no=\',num2str(sor_iter)],[\'simulation time=\',num2str(sim_time_sor)]})
  

end
end

1. if iterative_sol=1.function steady_jac is called with input arguments(t, told) and the output of the function is temperatue profile, simulation time and number of iteration of jacobi iterative method.

2. if iterative_sol=2 and 3, function steady_gs and steady_sor is called repectively and the output of the function is temperatue profile, simulation time and number of iteration of gauss seidal and SOR iterative method repectively.

Function steady_jac:

function [t,sim_time,jacobi_iter]=steady_jac(t,told) 
global nx ny tolerance
tic
error1=5;
jacobi_iter=1;
  while (error1>tolerance)          %iteration loop
    for i=2:nx-1                     %spacial loop
      for j=2:ny-1
       t(i,j)= 0.25*(told(i-1,j)+t(i+1,j)+told(i,j-1)+told(i,j+1));
      end
    end
   error1=max(max(abs(told-t)));    %error updation
   told=t;                          %temperature profile update for next iteration
   jacobi_iter=jacobi_iter+1;       %iteration count
  end
  
  sim_time=toc;                     %simulation time count
end

Function steady_gs:

function [t,sim_time,gs_iter]=steady_gs(t,told)
global nx ny tolerance
tic
error1=5;


gs_iter=1;
  while (error1>tolerance)           %iteration loop
    for i=2:nx-1                    %spacial loop
      for j=2:ny-1
       t(i,j)= 0.25*(t(i-1,j)+told(i+1,j)+t(i,j-1)+told(i,j+1));
      end
    end
   error1=max(max(abs(told-t)));  %error updation
   told=t;                       %temperature profile update for next iteration
   gs_iter=gs_iter+1;          %iteration count 
   
  end
  sim_time=toc;            %simulation time count
end

Function steady_sor:

function [t,sim_time,sor_iter]=steady_sor(t,told)
global nx ny tolerance omega
tic
error1=5;


sor_iter=1;
  while (error1>tolerance)          %iteration loop
    for i=2:nx-1                    %spacial loop
      for j=2:ny-1
        t(i,j)= told(i,j)*(1-omega)+omega*0.25*(t(i-1,j)+told(i+1,j)+t(i,j-1)+told(i,j+1));
      end
    end
   error1=max(max(abs(told-t)));    %error updation
   told=t;                          %temperature profile update for next iteration
   sor_iter=sor_iter+1;              %iteration count
   
  end
  sim_time=toc;                     %simulation time count

end

PROGRAM 1 OUPUT:

1. Temperature profile contour using jacobi iteravtive method

2. Temperature profile contour using gauss seidal iteravtive method

3. Temperature profile contour using SOR gauss seidal iteravtive method

Plotting conclusion:

1.The temperature profile for steady state is same for all three iterative method but the number of iterations taken for solution to converge is different for different iterative method.

2. The number of iterations for gauss seidal method is lesser (almost half) than jacobi iteration method.The number of iterations in gauss seidal is lesser because gauss seidal updates the temperature values calculated in an iteration (temperature values at i-1 and j-1 grid points) whereas jacobi iterative method keeps the old temperature value(at i-1 and j-1 grid points) for calculation of temperature at i grid point thereby taking more number of iterations to coverge.

3. The number of iterations for SOR gauss seidal is lesser than gauss seidal iterative method. The number of iterations for SOR gauss seidal is lesser because relaxation factor(omega) gives weightage to value calculated at previous iteration and current iteration value.If relaxation factor considered is greater than 1, SOR takes lesser iterations than gauss seidal to converge.If relaxation factor is equal to 1,SOR takes same number of iterations as gauss seidal.And if relaxation is lesser than 1, SOR takes more iterations to converge than gauss.Since lesser iterations to converge is desirable,usually relaxation factor is usually considered greater than 1.

4. Simulation time is least for gauss seidal as compared to jacobi and SOR in case of steady heat conduction.

CASE 2: Solution of transient 2D heat conduction equation by explicit approach

Assumptions taken:

1. Number of grid points in x and y direction, nx=ny=10 also since domain length is same, dx=dy.

2. Boundary conditions: Temperature at left : 400K

Temperature at top : 600K

Temperature at right: 800K

Temperature at bottom: 900K

And Average termeprature is taken at edges of the domain.

3. Absolute error criterion,tolerance= 1e-4

4. Time step values,dt=[1e-4 1e-3 1e-2 1e-1]

5. Thermal diffusivity=0.25

6. total time, T=50

PROGRAM 2:

clear all
close all
clc
%transient 2D heat conduction
nx=10;    %number of grid points
ny=nx;
T=50;     %total time

x=linspace(0,1,nx);     %grid generation
y= linspace(0,1,ny);
alpha=0.25;             %thermal diffusivity
dx= x(2)-x(1);          
dy=dx;                   %dx=dy
dt=[1e-4 1e-3 1e-2 1e-1];     %time step values

%applying boundary condition and initial condition
t= ones(nx,ny);         
t(1,1)=0.5*(900+400);                 %average temperature at edges
t(1,nx)=0.5*(900+800);
t(ny,1)=0.5*(400+600);
t(nx,ny)=0.5*(600+800);
t(2:ny-1,1)=400;
t(1,2:nx-1)=900;
t(2:ny-1,nx)=800;
t(ny,2:ny-1)=600;

told=t;
tolerance=1e-4;             %absolute error criterion

for l=1:4
    C(l)=(alpha*dt(l))*(1/(dx)^2+1/(dy)^2); %courant number
nt(l)=T/dt(l);                              %number of time steps
k1(l)= (alpha*dt(l))/(dx^2);                
k2(l)= (alpha*dt(l))/(dy^2);
tic
for time=1:nt(l)         %time loop
    for i=2:nx-1         %spacial loop
         for j=2:ny-1
             H=k1(l)*(told(i+1,j)-2*told(i,j)+told(i-1,j));
             V=k2(l)*(told(i,j+1)-2*told(i,j)+told(i,j-1));
            t(i,j)=told(i,j)+H+V;
         end
      end
       told=t;          %temperature update for each time loop
 end
t_exp=t;                %final temperature update
sim_time(l)=toc;        %simulation time calculation

figure(l)
   contourf(x,y,t_exp,\'showtext\',\'on\')
   hold on
   colorbar
   xlabel(\'X-Direction\')
   ylabel(\'Y-Direction\')
   title({\'transient heat conduction-explicit method\',[\'simulation time=\',num2str(sim_time(l))],[\'dt=\',num2str(dt(l))],[\'courant no=\',num2str(C(l))]})

end
   

PROGRAM 2 OUPUT:

1. Temperature profile for dt= 0.0001 and courant no= 0.00405

2. Temperature profile for dt=0.001 and courant no=0.0405

3. Temperature profile for dt=0.01 and courant no=0.405

4. Temperature profile for dt=0.1 and courant no=4.05

PROGRAM 2 Conclusion:

1. The solution obtained from transient explicit is similar to the steady state solution.

2. In transient explicit method,as the time step size increases the simulation time decreases considerably which is an desirable effect.

3.Increase in time step size leads to increase in courant number.As the courant number increases more than 0.5, the solution blows up or there is no solution.Therefore in order to get a stable solution, time step(dt) cannot be increased beyond a certain value.

CASE 3: Solution of transient heat conduction equation Implicitly using jacobi,gauss-seidal and SOR gauss seidal iterative solvers.

Assumptions taken:

1. Number of grid points in x and y direction, nx=ny=10 also since domain length is same, dx=dy.

2. Boundary conditions: Temperature at left : 400K

Temperature at top : 600K

Temperature at right: 800K

Temperature at bottom: 900K

And Average termeprature is taken at edges of the domain.

3. Absolute error criterion,tolerance= 1e-4

4. Time step values,dt=[ 1e-1 2e-1 1e+1 2e+1];

5. Thermal diffusivity=0.25

6. total time, T=500

PROGRAM 3:

clear all
close all
clc
global nx ny tolerance omega  %declaring global varaibles
nx=10;                       % number of grid points
ny=nx;
T=500;                       %total time

x=linspace(0,1,nx);            %grid generation
y= linspace(0,1,ny);
alpha=0.25;                    %thermal diffusivity
dx= x(2)-x(1);                   %grid size
dy=dx;      
dt=[ 1e-1 2e-1 1e+1 2e+1];       %time step size
omega=1.25;                       %relaxation factor

%applying boundary conditions and assuming temperature at remaining domain
%to be 1
t= ones(nx,ny);                        
t(1,1)=0.5*(900+400);           %average temperature at edges
t(1,nx)=0.5*(900+800);
t(ny,1)=0.5*(400+600);
t(nx,ny)=0.5*(600+800);
t(2:ny-1,1)=400;
t(1,2:nx-1)=900;
t(2:ny-1,nx)=800;
t(ny,2:ny-1)=600;

told=t;
t_pre_dt=t;
tolerance=1e-4;                 % absolute error criterion

for k=1:length(dt)              % for every value of dt

C(k)=(alpha*dt(k))*(1/(dx)^2+1/(dy)^2);    %courant number
nt(k)=T/dt(k);                             %number of grid points
k1(k)= (alpha*dt(k))/(dx^2);
k2(k)= (alpha*dt(k))/(dy^2);
term1(k)=1/(1+2*k1(k)+2*k2(k));
term2(k)=(k1(k))/(1+2*k1(k)+2*k2(k));
term3(k)=(k2(k))/(1+2*k1(k)+2*k2(k));

for iterative_sol=1:3
if iterative_sol==1      %for jacobi iterative method
    [T_jac,sim_time_jac,jacobi_iter]=transient_imp_jac(t,told,t_pre_dt,nt(k),term1(k),term2(k),term3(k));
   
    figure(1) 
    subplot(2,2,k)      
   contourf(x,y,T_jac,\'showtext\',\'on\')
   hold on
   colorbar
   xlabel(\'X-Direction\')
   ylabel(\'Y-Direction\')
   title({\'transient heat conduction-jacobi\',[\'iteration no=\',num2str(jacobi_iter)],[\'dt=\',num2str(dt(k))],[\'simulation time\',num2str(sim_time_jac)],[\'courant no=\',num2str(C(k))]})

  
   
end



if iterative_sol==2     %for gauss seidal iterative method
    
  [T_gs,sim_time_gs,gs_iter]=transient_imp_gs(t,told,t_pre_dt,nt(k),term1(k),term2(k),term3(k));
  
  figure(2)
  subplot(2,2,k)
   contourf(x,y,T_gs,\'showtext\',\'on\')
   hold on
   colorbar
   xlabel(\'X-Direction\')
   ylabel(\'Y-Direction\')
   title({\'transient heat conduction-gauss seidal\',[\'iteration no=\',num2str(gs_iter)],[\'dt=\',num2str(dt(k))],[\'simulation time\',num2str(sim_time_gs)],[\'courant no=\',num2str(C(k))]})

  
   
   
end

if iterative_sol==3 %for SOR gauss seidal iteravtive method
    
   [T_sor,sim_time_sor,sor_iter]=transient_imp_sor(t,told,t_pre_dt,nt(k),term1(k),term2(k),term3(k));
 
   figure(3)
   subplot(2,2,k)
   contourf(x,y,T_sor,\'showtext\',\'on\')
   hold on
   colorbar
   xlabel(\'X-Direction\')
   ylabel(\'Y-Direction\')
   title({\'transient heat conduction-SOR gauss seidal\',[\'iteration no=\',num2str(sor_iter)],[\'dt=\',num2str(dt(k))],[\'simulation time\',num2str(sim_time_sor)],[\'courant no=\',num2str(C(k))]})

  

   
end


end

end

1. dt is input as array,and the solution is calculated for every value of dt in an for loop.

2. if iterative_sol=1.function transient_imp_jac is called with input arguments(t, told,t_pre_dt) and the output of the function is temperatue profile, simulation time and number of iteration of jacobi iterative method.

3. if iterative_sol=2 and 3, function transient_imp_gs and transient_imp_sor is called repectively and the output of the function is temperatue profile, simulation time and number of iteration of gauss seidal and SOR iterative method repectively.

Function transient_imp_jac:

function [t,sim_time,jacobi_iter]= transient_imp_jac(t,told,t_pre_dt,nt,term1,term2,term3)
global nx ny  tolerance               %global varaibles
%jacobi
tic
 jacobi_iter=0;                        %iteration count
 error1=5;                            %initial error
    for time=1:nt                      %time loop
      while (error1>tolerance)         %iteration loop
       for i=2:nx-1                    %spatial loop
         for j=2:ny-1
          H=(told(i-1,j)+told(i+1,j));   
          V=(told(i,j-1)+told(i,j+1));
          t(i,j)= term1*t_pre_dt(i,j)+term2*H+term3*V;
         end
       end
       error1=max(max(abs(told-t)));
       told=t;                           % temperature upddate for every iteration
       jacobi_iter=jacobi_iter+1;      %iteration count
     
   
      end
      t_pre_dt=t;         %temperatue update for every time step
    end

sim_time=toc;              %simulation time
end

1. told is used to update the temperature profile for every iteration and t_pre_dt is used to update the temperature profile for every time loop.

Function transient_imp_gs:

function [t,sim_time,gs_iter]= transient_imp_gs(t,told,t_pre_dt,nt,term1,term2,term3)
global nx ny  tolerance
%gauss seidal
tic
 gs_iter=0;         %iteration count
 error1=5;           %initial error
    for time=1:nt          %time loop
      while (error1>tolerance)   %iteration loop
       for i=2:nx-1               %spatial loop
         for j=2:ny-1
             H=(t(i-1,j)+told(i+1,j));
             V=(t(i,j-1)+told(i,j+1));
             
          t(i,j)= term1*t_pre_dt(i,j)+term2*H+term3*V;
         end
       end
       error1=max(max(abs(told-t)));
       told=t;                        % temperature upddate for every iteration
       gs_iter=gs_iter+1;             %iteration count
     
   
      end
      t_pre_dt=t;                     %temperatue update for every time step
    end

sim_time=toc;                         %simulation time
end

Function transient_imp_sor:

function [t,sim_time,sor_iter]= transient_imp_sor(t,told,t_pre_dt,nt,term1,term2,term3)
global nx ny  tolerance omega               %global varaibles
%SOR 
tic
 sor_iter=0;                             %iteration count
 error1=5;                                %initial error
    for time=1:nt                         %time loop
      while (error1>tolerance)            %iteration loop
       for i=2:nx-1                       %spatial loop
         for j=2:ny-1
              H=(t(i-1,j)+told(i+1,j));
             V=(t(i,j-1)+told(i,j+1));
          tgs=(omega)*(term1*t_pre_dt(i,j)+term2*H+term3*V);   
          t(i,j)= (told(i,j))*(1-omega)+tgs;
         end
       end
       error1=max(max(abs(told-t)));
       told=t;                        % temperature upddate for every iteration
       sor_iter=sor_iter+1;           %iteration count
     
   
      end
      t_pre_dt=t;                      %temperatue update for every time step
    end

sim_time=toc;                         %simulation time
end

Program 3 Output:

1.Temperature profiles for transient heat conduction using jacobi iterative method for dt=0.1,0.2,10 and 20

2.Temperature profiles for transient heat conduction using gauss seidal iterative method for dt=0.1,0.2,10 and 20

3.Temperature profiles for transient heat conduction using SOR gauss seidal iterative method for dt=0.1,0.2,10 and 20

Program 3 conclusion:

1.For a particular value of time step(dt), all the three iterative methods(i.e. jacobi,gauss seidal and SOR gauss seidal) yeild similar solution.but the number of iterations and simulation time is different for the three iteravtive methods.

2. For a particular value of time step, SOR gauss seidal takes least number of iterations to converge and jacobi takes highest.And gauss seidal iteravtive method has the least simulation time while jacobi has the highest simulation time.

3. For a particular iteravtive method, As dt increases the simulation time decreases and also leads to increase in courant number.

Project conclusion:

1. In case of steady heat conduction and transient heat conduction, similar solutions are observed.Therefore we can conclude that transient method yeilds a solution similar to steady solution after a particular point in time marching. In other words, Time marching is also a method of attaining a steady solution in case of heat conduction.

2. In both steady and transient methods, SOR gauss seidal took least number of iterations to converge and jacobi takes the highest number of iterations to converge.The number of iterations in gauss seidal is lesser because gauss seidal updates the temperature values calculated in an iteration (temperature values at i-1 and j-1 grid points) whereas jacobi iterative method keeps the old temperature value(at i-1 and j-1 grid points) for calculation of temperature at i grid point thereby taking more number of iterations to coverge.

3. Comparing the implicit and explicit transient heat conduction, explcit method does not yeild a stable solution when courant number of the solution increases more than 0.5.Hence Explicit method has limitation value of time step and increasing time step beyond the limit leads to unstable solution.Whereas implicit solution does have a limitation in time step size(dt)(i.e implicit method is not limited to CFL condition), Hence higher time step values can be choosen in order to get quicker results.