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

                                                                                                                                 Date:  30-01-2021

 

STEADY STATE AND TRANSIENT STATE ANALYSIS OF 2D-HEAT CONDUCTION EQUATION

 

Main Objective: -

1. Solve steady state two-dimensional heat conduction equation by Using the three iterative techniques in implicit scheme.
2. Solve the Transient state two-dimensional heat conduction equation by using both implicit and explicit schemes.

 

Given input values: -

        Boundary conditions are

        Top Boundary = 600 K

       Bottom Boundary = 900 K

      Left Boundary = 400 K

      Right Boundary = 800 K

     Absolute error criteria is 1e-4

 

Description of domain: -

• Consider, domain as a square with equal number of points in both ‘x’ and ‘y’ directions.
• Here we are considering 10 points in both directions.
• Apply boundary conditions mentioned in problem as an input.
• Take average for corner points.
• For applying this iterative technique, we require initial guess of T-matrix.

 

                                                          

                                                                    Fig(1): - Domain discretisation

I) Steady state -Implicit Scheme-Iterative solvers: -

2D-Heat conduction equation in spatial coordinates or simply Laplacian heat conduction equation in 2D as follows

                     ꝺ²T/ ꝺx² + ꝺ²T/ ꝺy² = ∇²T =0

 

Mathematical formulation: -

 

• Now we have to discretise PDE to algebraic equations for numerical steady. This process is called Discretisation.
• We are discretising by the help of one of the numerical techniques called Finite difference method (FDM).
• Here we are using 2^nd order central difference scheme for 2^nd order derivatives in both ‘x’ and ‘y’ directions.

 

(T_R -2 T_P+ T_L) / Δx² + (T_T -2 T_P+ T_B) / Δy² = 0 ---------------1

Where,

       T_P = T (i, j) = our main aim to find the value of this point for different i,j values  

        T_R = T(i+1,j)

        T_L = T(i-1,j)

        T_T = T(i,j+1)

        T_B = T(i,j-1)

        Δx – distance between two successive points in x-direction which is constant .

        Δy – distance between two successive points in y-direction which is constant

 

We can understand equation1 by help of below figure

 

                                   

                                   

                                                  Fig(2): - Domain discretisation

• From equation 1 we can write that

 T_P = (1/k) *[(T_L + T_R/ Δx²) + (T_T + T_B/ Δy²) ]  ------2

Where ,

                     k=2*(Δx²+ Δy²)/(Δx²* Δy²)

1. Jacobi-iterative solver: -

• In this solver by the help of the old values we can march forward till we attain convergence.
• If the error value is less than tolerance (or) absolute error criteria, then we can say solution is converged.
• This technique is very helpful to find desired output in implicit scheme.
• Follow same mathematical formulation described above and apply this solving technique to equation2.
• For each iteration give old values as input and we can end up with output as new ‘T’ values. Again, for next iteration take new ‘T’ values as a input and find next iteration values. this process is continuing till solution meet converge criteria.

 

                        (T_P ) _new= (1/k) *[ (T_L + T_R/ Δx²)_old + (T_T + T_B/ Δy²)_old ]

• ‘Computational time’ or in other word ‘number of iterations’ is measuring factor for iterative solvers.
• We can see application of this method in MATLAB code and then we end up with steady state solution.

 

OUTPUT:-

 

 Fig (3): -Jacobi iterative solver steady state at 208 iterations

 

2. Gauss seidel-iterative solver: -

 

• This iterative solving technique is upgrade version to Jacobi technique.
• Instead of using old values in equation we use updated values as we can see in below equation.
• There is updated values in T_R and T_T .That’s why in that points (T_R and T_T) for each iteration ,we are still using old values. 

(T_P ) _new= (1/k) *[ (T_L + T_R/ Δx²) + (T_T + T_B/ Δy²) ]

(Or)

T(i,j)= (1/k) *[ (T(i-1,j) + T(i+1,j)/ Δx²) + (T(i,j+1)+ T(i,j-1)/ Δy²)]

     While applying Gauss Seidel Method

T(i,j)_new= (1/k) *[ (T(i-1,j)_new+ T(i+1,j)_old/ Δx²) + (T(i,j+1)_old+ T(i,j-1)_new/ Δy²)]

 

• We can see application of this method in MATLAB code and then we end up with steady state solution.

OUTPUT:-

Fig (4): -Gauss seidel iterative solver steady state at 112 iterations

3.Successive Over Relaxation-iterative solver: -

• This technique is advanced to both ‘Jacobi’ and ‘GS’ techniques.
• There is one factor called relaxation factor (omega ‘ω’)
• Based on choice of ‘ω’, solution accuracy and computational time to converge (number of iterations) is depend.

      (T(i,j))_new = (1- ω) T_old  + ω *{GS (or) jacobi}   ------------3

• As I mentioned above , based on ω value solution will change .

                              If ω>1 over relaxation

                              If ω<1 under relaxation

• By the help of jacobi and gauss seidal method we calculate sor That’s why in equation 3, in the place of GS and jacobi we have to any one of the this method. Generally ,we prefer GS method.
• ‘T_old’ in equation 3 will change for each iteration till solution reach steady state.
• We can see application of this method in MATLAB code and then we end up with steady state solution.

 

OUTPUT:-

   Fig (5): - SOR solver steady state at 29 iterations

 

Conclusion form above graphs: -

From above three plots we can say that regarding number of iterations

Jacobi > Gauss seidel > SOR

So, that means by the help of SOR iterative technique in implicit scheme, we can obtain steady state solution with less number of iterations (or) in other words ,SOR iterative technique took less computational time to achieve steady state solution.

 

 

II) TRANSIENT ANALYSIS: -

 

2D-Heat conduction equation in both spatial and temporal coordinates or simply transient heat conduction equation in 2D as follows

                    ꝺT/ ꝺt +α (ꝺ²T/ ꝺx² + ꝺ²T/ ꝺy²) = 0

                                              or

                                ꝺT/ ꝺt +α ∇²T =0

• Where ‘α’ is coefficient of thermal diffusivity.

Mathematical formulation: -

• Spatial part already we discuss above. Regarding temporal part we are using forward difference method with 1^st order approximation.

    ((T_pⁿ⁺¹ -T_pⁿ)/Δt) + α* [(T_R -2 T_P+ T_L) / Δx²] + α* [(T_T -2 T_P+ T_B) / Δy²] = 0 -----4

Where ‘n’ indicates number of time steps

• For each value of time step(n+1), find changes in space coordinates by help of previous timestep values(n).
• We can march forward by using two schemes explicit and implicit schemes.

Transient analysis by using Explicit scheme: -

• Generally, in this scheme we know all the values of RHS and by using them we can find LHS

Form equation 4 we can write

  T_pⁿ⁺¹(LHS) = {T_pⁿ + α *Δt* [(T_R -2 T_P+ T_L)ⁿ / Δx²]+ α *Δt* [(T_T -2 T_P+ T_B)ⁿ / Δy²]}(RHS)---5

Where ‘n’ indicates number of time steps

• Here we march forward by finding ‘T’ values at ‘n+1’ time step by the help of values of ‘T’ at ‘n’ time step.
• We can write equation 5 as

       T_pⁿ⁺¹ = T_pⁿ + k₁* [(T_R -2 T_P+ T_L)ⁿ / Δx²]+ k₂* [(T_T -2 T_P+ T_B)ⁿ / Δy²]

 

     Where

                      k₁ = α *(Δt/ Δx²)

                      k₂ = α *(Δt/ Δy²)

• This k₁ and k₂ are the CFL numbers or courant numbers along x and y directions.
• For this method, stability of solution is very important.
• Depends on CFL numbers only we can attain stable solution or converged solution .
• Here our criteria is K₁ + k2 <= 0.5
• We have to choose ‘Δt’ and ‘α’ such a way that that should meet above criteria.

 

 

OUTPUT: -

Fig (6): - At t=0.44 seconds we achieved steady state in transient analysis of 2d- heat conduction by using explicit method

Conclusion form above graphs: -

In explicit scheme we do time marching until we achieve steady state. Courant number (CFL) place an important role in explicit method to check stability criteria of solution for each time step.

 

Transient analysis by using Implicit scheme: -

• Instead of time marching by help of previous values like explicit scheme, we use iterative solvers to achieve steady state in this scheme.
• So, CFL numbers (along x and y direction) are in-effective while finding time marching solution.
• We can choose any value for ‘Δt’ and ‘α’ no stability restrictions.
• ‘Achieve steady state in less number of iterations’. We have to keep this statement in our mind while selecting ‘Δt’ and ‘α’ values.

Mathematical formulation: -

 T_pⁿ⁺¹ = T_pⁿ + k₁* [(T_R -2 T_P+ T_L)ⁿ⁺¹ / Δx²]+ k₂* [(T_T -2 T_P+ T_B)ⁿ⁺¹ / Δy²] -----6

• In this equation both sides we have (n+1) terms which are unknows
• From equation 6 we can write as

           T_pⁿ⁺¹=(term1*T_tⁿ +(term2*H)+(term3*V) -----7

Where,

 

                     T_pⁿ⁺¹ – new temperature values at n+1 time step

                      T_tⁿ - temperature values at previous time step (or) at ‘n’ time step

                      H=T_Lⁿ⁺¹ + T_Rⁿ⁺¹

                      V= T_Tⁿ⁺¹ + T_Bⁿ⁺¹

                      term1=1/(1+(2*k₁) + (2*k₂))

                      term2=k₁*term1

                      term3=k₂*term1

1. Jacobi-iterative solver: -

• Our iterative solver start with equation 7.We star with some gauss value of ‘Tⁿ⁺¹’
• Then new value for ‘Tⁿ⁺¹’ will come then we check convergence criteria which is same as steady state.
• By seeing MATLAB code which is shown last , we can able know what is space and time and convergence loop.
• This space loop will run until it won’t obey convergence criteria.
• Then after it will move to next time step and use previous time step as the initial guess for next time step.
• We can see application of this method in MATLAB code and then we end up with unsteady state solution .

OUTPUT:-

Fig (7): -Jacobi iterative solver achieves steady state in transient simulation at 832 iterations

 

2. Gauss seidel-iterative solver: -

• we can say this technique is updated technique of jacobi.
• same process like we did in steady state analysis.
• Here instead of using old values in space loop we use update values.
• Here we have another loop which is called time loop and after every space loop we have to update previous time step values.
• We can see application of this method in MATLAB code and then we end up with unsteady state solution.

 OUTPUT:-

Fig (8): - GS iterative solver achieves steady state in transient simulation at 473 iterations

3. Successive Over Relaxation-iterative solver: -

• This technique is advanced to both ‘Jacobi’ and ‘GS’ techniques.
• same process like we did in steady state analysis.
• By the help of jacobi and gauss seidal method we calculate sor That’s why in equation 3, in the place of GS and jacobi we have to any one of the this method. Generally ,we prefer GS method.
• We can see application of this method in MATLAB code and then we end up with unsteady state solution.

OUTPUT: -

Fig (9): - SOR iterative solver achieves steady state in transient simulation at 178 iterations

Conclusion form above graphs: -

From above three plots we can say that regarding number of iterations

Jacobi > Gauss seidel > SOR

So, that means by the help of SOR iterative technique in implicit scheme, we can obtain unsteady state (or) transient solution with less number of iterations (or) in other words, SOR iterative technique took less computational time to achieve steady state solution.

At same ‘Δt’ ,‘α’ values ,implicit scheme took more computational time compare to explicit scheme.

 

 

 

MATLAB code for both steady state and transient state analysis: -

 

% this below programmes devided into two categories steady state and transient
%steady-state have one main programme and 3 function programme. files should be in same folder
%transient have one main programme and four function programmes including explicit function programme. files should be in same folder

%%  FIRST PART

%% steady state 2d-heat conduction equation Numerical solution
% Main programme
clear all;
close all;
clc;
%% Assigning Geometric values
nx=10;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
[xx,yy] = meshgrid(x,y); %using meshgrid for future use in the contourf command
T=300*ones(10);% Here intial guess value is 300kelvin.Our convergence time will depend on our intial guess value  
%% Given values of Boundary conditions
T_lb = 400;
T_rb = 800;
T_tb = 600;
T_bb = 900;
error = 9e9;
tol=1e-4;
%% Assigning Boundary conditions 
T(1,:) = T_tb;
T(:,1) = T_lb;
T(:,end) = T_rb;
T(end,:) = T_bb;
T(1,1) = (T_tb+T_lb)/2;
T(1,10) = (T_tb+T_rb)/2;
T(10,1) = (T_bb+T_lb)/2;
T(10,10) = (T_bb+T_rb)/2;
%%
Told=T;
dx= x(2)-x(1);
dy=y(2)-y(1);
k=2*(dx^2+dy^2)/(dx*dy)^2;
%% Jacobi Method   
[jac_iter] = jacobi(dx,dy,nx,ny,k,tol,T,error,xx,yy);   %function calling
%% GS Method
[gs_iter] = gs(dx,dy,nx,ny,k,tol,T,error,xx,yy); %function calling
%% SOR Method
q=1.5;
[sor_iter] = sor(dx,dy,nx,ny,k,tol,T,error,xx,yy,q);%function calling
%% from here onwards we have to make different files for functions programme and main programme and put them in a single folder 

%% Jacobi method function programme

function [jac_iter] = jacobi(dx,dy,nx,ny,k,tol,T,error,xx,yy)
         Told=T;
         iterations=1;
      while(error > tol)%convergence loop
        for i=2:nx-1     %space loop along x-direcion
            for j=2:ny-1  %space loop along y-direction
            s1 =(1/k)*((Told(i-1,j)+Told(i+1,j))/(dx^2));
            s2 =(1/k)*((Told(i,j-1)+Told(i,j+1))/(dy^2));
            T(i,j) = s1+s2;
            end
        end
      error= max(max(abs(Told-T)));
      Told=T; 
      iterations=iterations+1;
      jac_iter = iterations;
      figure(1)
      [c,d]=contourf(xx,yy,T);
      colormap(jet)
      colorbar
      title_text=sprintf('Jacobi iteration number=%d at steady state',iterations);
      title(title_text)
      xlabel('X-axis')
      ylabel('Y-axis')
      clabel(c,d)
      set(gca,'YDIR','reverse') %required to correct the plot as by default it reverses the top and bottom values 
      end
      
end
%% Gauss seidel method function programme
function [gs_iter] = gs(dx,dy,nx,ny,k,tol,T,error,xx,yy)
         Told=T;
         iterations=1;
      while(error>tol)  %convergence loop 
        for i=2:nx-1     %space loop along x-direction
            for j=2:ny-1  %space loop along y-direction
            s1 =(1/k)*((T(i-1,j)+Told(i+1,j))/(dx^2));
            s2 =(1/k)*((T(i,j-1)+Told(i,j+1))/(dy^2));
            T(i,j) = s1+s2;
            end
        end
      error=max(max(abs(Told-T)));
      Told=T; 
      iterations=iterations+1;
      gs_iter =iterations;
      figure(2)
      [c,d]=contourf(xx,yy,T);
      colormap(jet)
      colorbar
      title_text=sprintf('gauss seidel iteration number=%d at steady state',iterations);
      title(title_text)
      xlabel('X-axis')
      ylabel('Y-axis')
      clabel(c,d)
      set(gca,'YDIR','reverse') %required to correct the plot as by default it reverses the top and bottom values 
      end
end
%% Successive Over Relaxation method function programme
function [sor_iter] = sor(dx,dy,nx,ny,k,tol,T,error,xx,yy,q)
         Told=T;
         iterations=1;
      while (error>tol)   %convergence loop
        for i=2:nx-1      %space loop along x-direction
            for j=2:ny-1   %space loop along y-direction
            s1 =(1/k)*((T(i-1,j)+Told(i+1,j))/(dx^2));
            s2 =(1/k)*((T(i,j-1)+Told(i,j+1))/(dy^2));
            z=s1+s2;
            T(i,j) = ((1-q).*Told(i,j))+(q.*z);  % q is the omega which is sor factor
            end
        end
      error= max(max(abs(Told-T)));
      Told=T; 
      iterations=iterations+1;
      sor_iter = iterations;
      figure(3)
      [c,d]=contourf(xx,yy,T);
      colormap(jet)
      colorbar
      title_text=sprintf('Successive Over Relaxation iteration number=%d at steady state',iterations);
      title(title_text)
      xlabel('X-axis')
      ylabel('Y-axis')
      clabel(c,d)
      set(gca,'YDIR','reverse') %required to correct the plot as by default it reverses the top and bottom values 
      end

end


%%  SECOND PART

% Explicit function programme

function [explicit_transient_time]=explicit(nx,ny,T,xx,yy,nt,Told,k1,k2,dt,t)

for k=1:nt % time loop
     for i=2:nx-1%space loop
            for j=2:ny-1%space loop
                z1=Told(i-1,j)-2*Told(i,j)+Told(i+1,j);
                z2=Told(i,j+1)-2*Told(i,j)+Told(i,j-1);
                T(i,j)=Told(i,j)+k1*z1+k2*z2;
            end
     end
      Told=T;
      figure(4)
      [c,d]=contourf(xx,yy,T);
      colormap(jet)
      colorbar
      title_text=sprintf('Transient behaviour of 2D-heat conduction at=%d seconds',k*dt);
      title(title_text)
      xlabel('X-axis')
      ylabel('Y-axis')
      clabel(c,d)
      set(gca,'YDIR','reverse') %required to correct the plot as by default it reverses the top and bottom values    
   
end
explicit_transient_time = t; %final time we need to achive steady state
end

%Implicit jacobi transient state function programme
function [implicit_transient_jacobi_iter]=implicit_jac_transient(nx,ny,tol,T,xx,yy,nt,Told,term1,term2,term3,T_t) 
iteration=1;
 for k=1:nt %time loop
     error=100;
     while(error > tol) % convergence loop
     for i=2:nx-1  % space loop
            for j=2:ny-1   % space loop
                H=Told(i-1,j)+Told(i+1,j);
                V=Told(i,j-1)+Told(i,j+1);
                T(i,j)=(term1*T_t(i,j))+(term2*H)+(term3*V);
            end
     end
     error= max(max(abs(Told-T)));
     Told=T;
     iteration=iteration+1;
     figure(5)
     [c,d]=contourf(xx,yy,T);
      colormap(jet)
      colorbar
      title_text=sprintf('Jacobi iteration number=%d at unsteady state',iteration);
      title(title_text)
      xlabel('X-axis')
      ylabel('Y-axis')
      clabel(c,d)
      set(gca,'YDIR','reverse') %required to correct the plot as by default it reverses the top and bottom values  
     implicit_transient_jacobi_iter = iteration;
     end
      T_t=T;
 end
end

%Implicit Gauss seidel transient state function programme

function [implicit_transient_gs_iter]=implicit_gs_transient(nx,ny,tol,T,xx,yy,nt,Told,term1,term2,term3,T_t) 
iteration=1;
 for k=1:nt % time loop
     error=100;
     while(error > tol)  % convergence loop
     for i=2:nx-1       % time loop
            for j=2:ny-1 % time loop
                H=T(i-1,j)+Told(i+1,j);
                V=T(i,j-1)+Told(i,j+1);
                T(i,j)=(term1*T_t(i,j))+(term2*H)+(term3*V);
            end
     end
     error= max(max(abs(Told-T)));
     Told=T;
     iteration=iteration+1;
     figure(6)
     [c,d]=contourf(xx,yy,T);
      colormap(jet)
      colorbar
      title_text=sprintf('GS iteration number=%d at unsteady state',iteration);
      title(title_text)
      xlabel('X-axis')
      ylabel('Y-axis')
      clabel(c,d)
      set(gca,'YDIR','reverse') %required to correct the plot as by default it reverses the top and bottom values  
     implicit_transient_gs_iter = iteration;
     end
      T_t=T;
 end
end

%Implicit SOR transient state function programme

function [implicit_transient_sor_iter]=implicit_sor_transient(nx,ny,tol,T,xx,yy,nt,Told,term1,term2,term3,T_t,q) 
iteration=1;
 for k=1:nt   % time loop
     error=100;
     while(error > tol)   % convergence loop
     for i=2:nx-1          % space loop
            for j=2:ny-1  % space loop
                H=T(i-1,j)+Told(i+1,j);
                V=T(i,j-1)+Told(i,j+1);
                z=(term1*T_t(i,j))+(term2*H)+(term3*V);
                T(i,j) = ((1-q).*Told(i,j))+(q.*z);  % q is the SOR factor for over relaxation
            end
     end
     error= max(max(abs(Told-T)));
     Told=T;
     iteration=iteration+1;
     figure(7)
     [c,d]=contourf(xx,yy,T);
      colormap(jet)
      colorbar
      title_text=sprintf('SOR iteration number=%d at unsteady state',iteration);
      title(title_text)
      xlabel('X-axis')
      ylabel('Y-axis')
      clabel(c,d)
      set(gca,'YDIR','reverse') %required to correct the plot as by default it reverses the top and bottom values  
     implicit_transient_sor_iter = iteration;
     end
      T_t=T;
 end
end

% Function programme '.m'files and main programme '.m'file should be in same folder

%%Main programme

clear all;
close all;
clc;
%% Assigning Geometric values
nx=10;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
[xx,yy] = meshgrid(x,y); %using meshgrid for future use in the contourf command
T=300*ones(10);% Here intial guess value is 300kelvin.Our convergence time will depend on our intial guess value  
%% Given values of Boundary conditions
T_lb = 400;
T_rb = 800;
T_tb = 600;
T_bb = 900;
error = 9e9;
tol=1e-4;
%% Assigning Boundary conditions 
T(1,:) = T_tb;
T(:,1) = T_lb;
T(:,end) = T_rb;
T(end,:) = T_bb;
T(1,1) = (T_tb+T_lb)/2;
T(1,10) = (T_tb+T_rb)/2;
T(10,1) = (T_bb+T_lb)/2;
T(10,10) = (T_bb+T_rb)/2;
%%
Told=T;
dx= x(2)-x(1);
dy=y(2)-y(1);
%% Transient simulation EXPLICIT Scheme
dt=0.001;
t=0.44;
nt=t/dt;
alpha=1.4;
k1=alpha*(dt/dx^2);%CFL1 number along x-direction.Which is very important for explicit method to maintain stability
k2=alpha*(dt/dy^2);%CFL2 number along y-direction.Which is very important for explicit method to maintain stability
% CFl1+CFL2 <= 0.5 is the stability criteria for this explicit scheme
%%
[explicit_transient_time]=explicit(nx,ny,T,xx,yy,nt,Told,k1,k2,dt,t);%calling explicit function
%% Transient simulation IMPLICIT Scheme 
% Implicit scheme is independent of CFL number.
% dt,nt,t and alpha values are such away that we can achive steady state.
dt=0.1;
t=1;
nt=t/dt;
alpha=1.4;
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);
T_t=T;
term1=1/(1+(2*k1)+(2*k2));
term2=k1*term1;
term3=k2*term1;
 %% Jacobi_solver
 [implicit_transient_jacobi_iter]=implicit_jac_transient(nx,ny,tol,T,xx,yy,nt,Told,term1,term2,term3,T_t);%calling functio
 %% GS_solver
 [implicit_transient_gs_iter]=implicit_gs_transient(nx,ny,tol,T,xx,yy,nt,Told,term1,term2,term3,T_t);% calling function
 %% SOR_solver
 q=1.5; % q is noting but omega,which is main reason for over relaxation
 [implicit_transient_sor_iter]=implicit_sor_transient(nx,ny,tol,T,xx,yy,nt,Told,term1,term2,term3,T_t,q);%calling function