2D Heat Conduction Simulation

Steady-State Heat Conduction

Steady-state conduction is the form of conduction that happens when the temperature difference(s) driving the conduction are constant, so that (after an equilibration time), the spatial distribution of temperatures (temperature field) in the conducting object does not change any further. Thus, all partial derivatives of temperature concerning space may either be zero or have nonzero values, but all derivatives of temperature at any point concerning time are uniformly zero. In steady-state conduction, the amount of heat entering any region of an object is equal to the amount of heat coming out (if this were not so, the temperature would be rising or falling, as thermal energy was tapped or trapped in a region).

For example, a bar may be cold at one end and hot at the other, but after a state of steady-state conduction is reached, the spatial gradient of temperatures along the bar does not change any further, as time proceeds. Instead, the temperature remains constant at any given cross-section of the rod normal to the direction of heat transfer, and this temperature varies linearly in space in the case where there is no heat generation in the rod.

2D Steady-State Conduction Equation

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

`(T_L-2T_P+T_R)/(Deltax^2)+(T_T-2T_P+T_B)/(Deltay^2)=0`

`(2T_P)/(Deltax^2)+(2T_P)/(Deltay^2)=(T_L+T_R)/(Deltax^2)+(T_T+T_B)/(Deltay^2)`

`T_P(2(Deltax^2+Deltay^2))/(Deltax^2Deltay^2)=(T_L+T_R)/(Deltax^2)+(T_T+T_B)/(Deltay^2)`

Let  `k=(2(Deltax^2+Deltay^2))/(Deltax^2Deltay^2)`  

`T_P=1/k((T_L+T_R)/(Deltax^2)+(T_T+T_B)/(Deltay^2))`

The temperature at the centre node is dependent upon the four adjoining nodes in the domain.

Declaring the Initial and Boundary Condition

A domain of length one meter is taken in x and y direction and it is discretised taking the number of 10 nodes in x and y direction respectively.

For assigning the boundary condition to the discretised domain, a 10X10  matrix is generated and the values are made unity and then the respective boundary values are multiplied with each respective boundary and at the junction of two boundary average value is assigned.

An error variable is made for making the entry into the loop and the given tolerance is also assigned to a variable.

A variable is used to store the current values of temperature which then acts as old values for the new iteration.

An assumed value of Temperature is used as the initial condition which is solved iteratively.

%Code to solve 2nd order Unsteady State 2D Heat Conduction Equation
% d^2T/dx^2+d^2T/dy^2=0

% Inputs
n_x=10;
n_y=10;

%Domain
x=linspace(0,1,n_x);
dx=x(2)-x(1);
y=linspace(0,1,n_y);
dy=y(2)-y(1);
T=ones(n_x,n_y);

%Initial and Boundary Condition
T=T*350;          %Assumed Value
T(1,:)=900;       %T_Lower=900K
T(n_x,:)=600;     %T_Top=600K
T(:,1)=400;       %T_Left=400K
T(:,n_y)=800;     %T_Right=800K

%Temperature at junction of two boundaries
T(1,1)=(900+400)/2;       
T(n_y,1)=(400+600)/2;
T(1,n_x)=(600+800)/2;
T(n_y,n_x)=(900+800)/2;

T_old=T;
error_jac=error_gs=error_sor=9e9;
tol=1e-4;

Solving using Iterative Technique

We can solve the above steady-state heat conduction equation using 3 different iterative methods, namely Jacobi Method, Gauss-Seidel Method and Successive over-relaxation method. We assign a value to each solver and to choose a solver we create a variable and assign a value corresponding to the solver.

Jacobi Method

`T_(i,j)=1/k[(T_(old_(i-1,j))+T_(old_(i+1,j)))/dx^2+(T_(old_(i,j-1))+T_(old_(i,j+1)))/dy^2]`

Gauss Siedel Method

`T_(i,j)=1/k[(T_((i-1,j))+T_(old_(i+1,j)))/dx^2+(T_((i,j-1))+T_(old_(i,j+1)))/dy^2]`

Successive over-relaxation Method

`T_(i,j)=(1-alpha)(T_(old_(i,j)))+alpha T_((i,j) _(Jacobi or Gauss))`

To find the solution two for loops are formulated inside a while loop. The while loop checks the error against the given tolerance and stops the loops if the error is less than tolerance. The two for loops is also called the special loop as it depends upon the grid points in the x and y-direction. The loops can be used to find all the temperature along the y-direction by fixing the x value or vice-versa. At the end of the spacial loop, we again find the new error and updates the temperature values to a variable that acts as the old values for the next iteration.

A contour plot is used within the while loop to plot the values of Temperature along x and y-direction.

Similarly, the code is repeated for the Gauss Siedel and Successive over-relaxation method.

Program for Steady-State Solution 

clear all
close all 
clc

%Code to solve 2nd order Unsteady State 2D Heat Conduction Equation
% d^2T/dx^2+d^2T/dy^2=0

% Inputs
n_x=10;
n_y=10;

%Domain
x=linspace(0,1,n_x);
dx=x(2)-x(1);
y=linspace(0,1,n_y);
dy=y(2)-y(1);
T=ones(n_x,n_y);

%Initial and Boundary Condition
T=T*350;          %Assumed Value
T(1,:)=900;       %T_Lower=900K
T(n_x,:)=600;     %T_Top=600K
T(:,1)=400;       %T_Left=400K
T(:,n_y)=800;     %T_Right=800K

%Temperature at junction of two boundaries
T(1,1)=(900+400)/2;       
T(n_y,1)=(400+600)/2;
T(1,n_x)=(600+800)/2;
T(n_y,n_x)=(900+800)/2;

T_old=T;
error_jac=error_gs=error_sor=9e9;
tol=1e-4;

%Calculating the constant k
k=2*(dx^2+dy^2)/(dx^2*dy^2);

%Choosing the solver
iterative_solver=3;

%Jacobi Method
if iterative_solver==1
  jacobi_ct=1;
  while(error_jac>tol)
    for i=2:(length(x)-1)
        for j=2:(length(y)-1)
          T(i,j)=(1/k)*(((T_old(i-1,j)+T_old(i+1,j))/dx^2)+((T_old(i,j-1)+T_old(i,j+1))/dy^2));
        end
    end 
    error_jac=max(max(abs(T_old.-T)));
    T_old=T;
    jacobi_ct=jacobi_ct+1;
    figure(1)
    contourf(x,y,T)
    colorbar
    xlabel('X-Axis')
    ylabel('Y-Axis')
    zlabel('Temperature')
    title_text=sprintf('Steady State Solution by Jacobi Method at Iteration Number=%d',jacobi_ct);
    title(title_text)
    pause(0.009)
  end
end

%Gauss Siedel Method
if iterative_solver==2
  gs_ct=1;
  while(error_gs>tol)
    for i=2:(length(x)-1)
        for j=2:(length(y)-1)
          T(i,j)=(1/k)*(((T(i-1,j)+T_old(i+1,j))/dx^2)+((T(i,j-1)+T_old(i,j+1))/dy^2));
        end
    end 
    error_gs=max(max(abs(T_old.-T)));
    T_old=T;
    gs_ct=gs_ct+1;
    figure(2)
    contourf(x,y,T)
    colorbar
    xlabel('X-Axis')
    ylabel('Y-Axis')
    zlabel('Temperature')
    title_text=sprintf('Steady State Solution by Gauss Siedel Method at Iteration Number=%d',gs_ct);
    title(title_text)
    pause(0.009)
  end
end

%SOR Method
if iterative_solver==3
  sor_ct=1;
  omega=1.5;  %constant
  while(error_sor>tol)
    for i=2:(length(x)-1)
        for j=2:(length(y)-1)
          T(i,j)=(1/k)*(((T(i-1,j)+T_old(i+1,j))/dx^2)+((T(i,j-1)+T_old(i,j+1))/dy^2));
          T(i,j)=(1-omega)*T_old(i,j)+omega*(T(i,j));
        end
    end 
    error_sor=max(max(abs(T_old.-T)));
    T_old=T;
    sor_ct=sor_ct+1;
    figure(3)
    contourf(x,y,T)
    colorbar
    xlabel('X-Axis')
    ylabel('Y-Axis')
    zlabel('Temperature')
    title_text=sprintf('Steady State Solution by SOR Method at Iteration Number=%d',sor_ct);
    title(title_text)
    pause(0.009)
  end
end

 Results

The temperature distribution at a steady state is for the three iterative methods are shown below. The temperature distribution for the steady-state is independent of the time thus the profile of temperature remain constant across the domain and this is found using the iterative solver.  

Steady-State Solution using the Jacobi method took 205 iterations with an error less than `1e^-4` .

Steady-State Solution using Gauss Siedel method took 110 iterations with an error less than `1e^-4` .             

Steady-State Solution using successive over-relaxation method took 28 iterations with an error less than `1e^-4` . 

Transient Heat Conduction 

In a time period if the temperatures changes in time at any place within an object, the mode of thermal energy flow is termed transient conduction. Another term is "non-steady-state" conduction, referring to the time-dependence of temperature fields in an object. Non-steady-state situations appear after an imposed change in temperature at a boundary of an object. They may also occur with temperature changes inside an object, as a result of a new source or sink of heat suddenly introduced within an object, causing temperatures near the source or sink to change in time.

2D Transient-State Conduction Equation

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

Explicit Method   

`T_(i,j)^(n+1)=T_(i,j)^n+k_1(T_(i,j+1)^n-2T_(i,j)^n+T_(i,j-1)^n)+k_2(T_(i+1,j)^n-2T_(i,j)^n+T_(i-1,j)^n)`  where

`k_1=alpha(Deltat)/(Deltax^2` 

`k_2=alpha(Deltat)/(Deltay^2`    

Declaring the Initial and Boundary Condition

A domain of length one meter is taken in x and y direction and it is discretised taking the number of 10 nodes in x and y direction respectively.

For assigning the boundary condition to the discretised domain, a 10X10  matrix is generated and the values are made unity and then the respective boundary values are multiplied with each respective boundary and at the junction of two boundary average value is assigned.

An error variable is made for making the entry into the loop and the given tolerance is also assigned to a variable.

A variable is used to store the current values of temperature which then acts as old values for the new iteration.

An assumed value of Temperature is used as the initial condition which is solved iteratively.

clear all
close all 
clc

%Code to solve 2nd order Steady State 2D Heat Conduction Equation
% dT/dt+d^2T/dx^2+d^2T/dy^2=0

% Inputs
n_x=10;
n_y=10;

%Domain
x=linspace(0,1,n_x);
dx=x(2)-x(1);
y=linspace(0,1,n_y);
dy=y(2)-y(1);
T=ones(n_x,n_y);

%Initial and Boundary Condition
T=T*350;
T(1,:)=900;       %T_Lower=900K
T(n_x,:)=600;     %T_Top=600K
T(:,1)=400;       %T_Left=400K
T(:,n_y)=800;     %T_Right=800K

T(1,1)=(900+400)/2;
T(n_y,1)=(400+600)/2;
T(1,n_x)=(600+800)/2;
T(n_y,n_x)=(900+800)/2;

%Assumed Inputs
dt=0.002;
alpha=1.4;
omega=1.5;

T_old=T;
T_pre=T;
error_ex=9e9;
error_gs=9e9;
error_jac=9e9;
error_sor=9e9;
tol=1e-4;
k1=(alpha*dt)/dx^2;
k2=(alpha*dt)/dy^2;

Implicit Method

We can solve the above steady-state heat conduction equation using 3 different iterative methods, namely Jacobi Method, Gauss-Seidel Method and Successive over-relaxation method. We assign a value to each solver and to choose a solver we create a variable and assign a value corresponding to the solver.

`T_(i,j)^(n+1)=T_(i,j)^(n)+k_1(T_(i,j+1)^(n+1)-2T_(i,j)^(n+1)+T_(i,j-1)^(n+1))+k_2(T_(i+1,j)^(n+1)-2T_(i,j)^(n+1)+T_(i-1,j)^(n+1))`      

`T_(i,j)^(n+1)=[1/(1+2k_1+2k_2)][T_(i,j)^(n)+k_1(T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1))+k_2(T_(i+1,j)^(n+1)+T_(i-1,j)^(n+1))`   

Solving using Iterative Technique

Jacobi Method

`T_(i,j)^(n+1)=[1/(1+2k_1+2k_2)][T_(i,j)^(n)+k_1(T_(old_(i,j+1))^(n+1)+T_(old_(i,j-1))^(n+1))+k_2(T_(old_(i+1,j))^(n+1)+T_(old_(i-1,j))^(n+1))`

Gauss Siedel Method

`T_(i,j)^(n+1)=[1/(1+2k_1+2k_2)][T_(i,j)^(n)+k_1(T_(old_(i,j+1))^(n+1)+T_(i,j-1)^(n+1))+k_2(T_(old_(i+1,j))^(n+1)+T_(i-1,j))(n+1))`

Successive over-relaxation Method

`T_(i,j)^(n+1)=(1-alpha)(T_(old_(i,j))^(n+1))+alpha T_((i,j) _(Jacobi or Gauss))^(n+1)`

To find the solution two for loops are formulated inside a while loop. The while loop checks the error against the given tolerance and stops the loops if the error is less than tolerance. The two for loops is also called the special loop as it depends upon the grid points in the x and y-direction. The loops can be used to find all the temperature along the y-direction by fixing the x value or vice-versa. At the end of the spacial loop, we again find the new error and updates the temperature values to a variable that acts as the old values for the next iteration.

A contour plot is used within the while loop to plot the values of Temperature along x and y-direction.

Similarly, the code is repeated for the Gauss Siedel and Successive over-relaxation method.

iterative_solver=1;

%Explicit Method
if iterative_solver==1
  t=0+dt;
  ct=1;
  while(error_ex>tol)
    for i=2:(length(x)-1)
        for j=2:(length(y)-1)
            T(i,j)=T_pre(i,j)+k1*(T_pre(i,j+1)-2*T_pre(i,j)+T_pre(i,j-1))+k2*(T_pre(i+1,j)-2*T_pre(i,j)+T_pre(i-1,j));
        end
    end 
    ct=ct+1;
    error_ex=max(max(abs(T_pre-T)));
    T_pre=T;
    t=t+dt;
    figure(1)
    contourf(x,y,T)
    colorbar
    xlabel('X-Axis')
    ylabel('Y-Axis')
    zlabel('Temperature')
    title_text=sprintf('Transient state solution using Explicit scheme is obatined at time=%d seconds',t)
    title(title_text)
    pause(0.009)
  end
 end

 %Implicit Method

%Jacobi Method
if iterative_solver==2
  jacobi_ct=1;
  t=0+dt;
  while (error_jac>tol)
    error_jac=9e9;
    while(error_jac>tol)
      for i=2:(length(x)-1)
          for j=2:(length(y)-1)
           T(i,j)=(1/(1+2*k1+2*k2))*(T_pre(i,j)+k1*(T_old(i,j+1)+T_old(i,j-1))+k2*(T_old(i+1,j)+T_old(i-1,j)));
          end
      end 
      error_jac=max(max(abs(T_old-T)));
      T_old=T;   
    end
    jacobi_ct=jacobi_ct+1;
    error_jac=max(max(abs(T_pre-T)));
    T_pre=T;
    t=t+dt;
    figure(1)
    contourf(x,y,T)
    colorbar
    xlabel('X-Axis')
    ylabel('Y-Axis')
    zlabel('Temperature')
    title_text=sprintf('Transient state solution using Implicit scheme with Jacobi method is obatined at time=%d seconds',t);
    title(title_text)
    pause(0.009)
  end
end

%Gauss Siedel Method
if iterative_solver==3
  gs_ct=1;
  t=0+dt;
  while (error_gs>tol)
    error_gs=9e9;
    while(error_gs>tol)
      for i=2:(length(x)-1)
          for j=2:(length(y)-1)
           T(i,j)=(1/(1+2*k1+2*k2))*(T_pre(i,j)+k1*(T_old(i,j+1)+T(i,j-1))+k2*(T_old(i+1,j)+T(i-1,j)));
          end
      end 
      error_gs=max(max(abs(T_old-T)));
      T_old=T;   
    end
    gs_ct=gs_ct+1;
    error_gs=max(max(abs(T_pre-T)));
    T_pre=T;
    t=t+dt;
    figure(1)
    contourf(x,y,T)
    colorbar
    xlabel('X-Axis')
    ylabel('Y-Axis')
    zlabel('Temperature')
    title_text=sprintf('Transient state solution Implicit scheme with Gauss Siedel method is obatined at time=%d seconds',t);
    title(title_text)
    pause(0.009)
  end
end

%SOR Method
if iterative_solver==4
  sor_ct=1;
  t=0+dt;
  while (error_sor>tol)
    error_sor=9e9;
    while(error_sor>tol)
      for i=2:(length(x)-1)
          for j=2:(length(y)-1)
           T(i,j)=(1/(1+2*k1+2*k2))*(T_pre(i,j)+k1*(T_old(i,j+1)+T(i,j-1))+k2*(T_old(i+1,j)+T(i-1,j)));
           T(i,j)=(1-omega)*T_old(i,j)+omega*(T(i,j));
          end
      end 
      error_sor=max(max(abs(T_old-T)));
      T_old=T;   
    end
    sor_ct=sor_ct+1;
    error_sor=max(max(abs(T_pre-T)));
    T_pre=T;
    t=t+dt;
    figure(1)
    contourf(x,y,T)
    colorbar
    xlabel('X-Axis')
    ylabel('Y-Axis')
    zlabel('Temperature')
    title_text=sprintf('Transient state solution using Implicit scheme with SOR method is obatined at time=%d seconds',t);
    title(title_text)
    pause(0.009)
  end
end

Program for Transient-State Solution

clear all
close all 
clc

%Code to solve 2nd order Steady State 2D Heat Conduction Equation
% dT/dt+d^2T/dx^2+d^2T/dy^2=0

% Inputs
n_x=10;
n_y=10;

%Domain
x=linspace(0,1,n_x);
dx=x(2)-x(1);
y=linspace(0,1,n_y);
dy=y(2)-y(1);
T=ones(n_x,n_y);

%Initial and Boundary Condition
T=T*350;
T(1,:)=900;       %T_Lower=900K
T(n_x,:)=600;     %T_Top=600K
T(:,1)=400;       %T_Left=400K
T(:,n_y)=800;     %T_Right=800K

T(1,1)=(900+400)/2;
T(n_y,1)=(400+600)/2;
T(1,n_x)=(600+800)/2;
T(n_y,n_x)=(900+800)/2;

%Assumed Inputs
dt=0.002;
alpha=1.4;
omega=1.5;

T_old=T;
T_pre=T;
error_ex=9e9;
error_gs=9e9;
error_jac=9e9;
error_sor=9e9;
tol=1e-4;
k1=(alpha*dt)/dx^2;
k2=(alpha*dt)/dy^2;


iterative_solver=1;

%Explicit Method
if iterative_solver==1
  t=0+dt;
  ct=1;
  while(error_ex>tol)
    for i=2:(length(x)-1)
        for j=2:(length(y)-1)
            T(i,j)=T_pre(i,j)+k1*(T_pre(i,j+1)-2*T_pre(i,j)+T_pre(i,j-1))+k2*(T_pre(i+1,j)-2*T_pre(i,j)+T_pre(i-1,j));
        end
    end 
    ct=ct+1;
    error_ex=max(max(abs(T_pre-T)));
    T_pre=T;
    t=t+dt;
    figure(1)
    contourf(x,y,T)
    colorbar
    xlabel('X-Axis')
    ylabel('Y-Axis')
    zlabel('Temperature')
    title_text=sprintf('Transient state solution using Explicit scheme is obatined at time=%d seconds',t)
    title(title_text)
    pause(0.0000009)
  end
 end

 %Implicit Method

%Jacobi Method
if iterative_solver==2
  jacobi_ct=1;
  t=0+dt;
  while (error_jac>tol)
    error_jac=9e9;
    while(error_jac>tol)
      for i=2:(length(x)-1)
          for j=2:(length(y)-1)
           T(i,j)=(1/(1+2*k1+2*k2))*(T_pre(i,j)+k1*(T_old(i,j+1)+T_old(i,j-1))+k2*(T_old(i+1,j)+T_old(i-1,j)));
          end
      end 
      error_jac=max(max(abs(T_old-T)));
      T_old=T;   
    end
    jacobi_ct=jacobi_ct+1;
    error_jac=max(max(abs(T_pre-T)));
    T_pre=T;
    t=t+dt;
    figure(1)
    contourf(x,y,T)
    colorbar
    xlabel('X-Axis')
    ylabel('Y-Axis')
    zlabel('Temperature')
    title_text=sprintf('Transient state solution using Implicit scheme with Jacobi method is obatined at time=%d seconds',t);
    title(title_text)
    pause(0.009)
  end
end

%Gauss Siedel Method
if iterative_solver==3
  gs_ct=1;
  t=0+dt;
  while (error_gs>tol)
    error_gs=9e9;
    while(error_gs>tol)
      for i=2:(length(x)-1)
          for j=2:(length(y)-1)
           T(i,j)=(1/(1+2*k1+2*k2))*(T_pre(i,j)+k1*(T_old(i,j+1)+T(i,j-1))+k2*(T_old(i+1,j)+T(i-1,j)));
          end
      end 
      error_gs=max(max(abs(T_old-T)));
      T_old=T;   
    end
    gs_ct=gs_ct+1;
    error_gs=max(max(abs(T_pre-T)));
    T_pre=T;
    t=t+dt;
    figure(1)
    contourf(x,y,T)
    colorbar
    xlabel('X-Axis')
    ylabel('Y-Axis')
    zlabel('Temperature')
    title_text=sprintf('Transient state solution Implicit scheme with Gauss Siedel method is obatined at time=%d seconds',t);
    title(title_text)
    pause(0.009)
  end
end

%SOR Method
if iterative_solver==4
  sor_ct=1;
  t=0+dt;
  while (error_sor>tol)
    error_sor=9e9;
    while(error_sor>tol)
      for i=2:(length(x)-1)
          for j=2:(length(y)-1)
           T(i,j)=(1/(1+2*k1+2*k2))*(T_pre(i,j)+k1*(T_old(i,j+1)+T(i,j-1))+k2*(T_old(i+1,j)+T(i-1,j)));
           T(i,j)=(1-omega)*T_old(i,j)+omega*(T(i,j));
          end
      end 
      error_sor=max(max(abs(T_old-T)));
      T_old=T;   
    end
    sor_ct=sor_ct+1;
    error_sor=max(max(abs(T_pre-T)));
    T_pre=T;
    t=t+dt;
    figure(1)
    contourf(x,y,T)
    colorbar
    xlabel('X-Axis')
    ylabel('Y-Axis')
    zlabel('Temperature')
    title_text=sprintf('Transient state solution using Implicit scheme with SOR method is obatined at time=%d seconds',t);
    title(title_text)
    pause(0.009)
  end
end

 Results

The temperature distribution for the transient state varies with respect to time in the domain and with the time the transient state changes to a steady-state.

Transient-State Solution using Explicit scheme takes 0.45 seconds to reach steady-state.

Transient-State Solution using Implicit scheme with the Jacobi method takes 0.464 seconds to reach steady-state.

Transient-State Solution using Implicit scheme with the Gauss Siedel method takes 0.466 seconds to reach steady-state.

Transient-State Solution using Implicit scheme with the SOR method takes 0.474 seconds to reach steady-state.