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

Aim: To solve 2D heat conduction equation for steady state and transient state using various iterative techiniques like Jacobi, Gauss Seidel & Successive over relaxation methods using MATLAB.

Problem Define

Heat conduction equation is define by pde
`del T/delt = alpha*(del^2T/del x + del^2T/del y + del^2T/del z)`

Which is heat distribution in 3D space with period of time.
Here alpha is known to be thermal diffusivity.


Methods used

Explicit Method: When a direct computation is made for an unknown quantity in terms of known quantity This computation is called Explicit. This method requires stability criteria which when leads to stable solution.As the stability criteria small steps has to be taken to obtain stable solutions. With decrease in time step the number of computation increases hence time required for solution to converge is less and accuracy is high.

Implicit Method: When a computation is made for dependent variable defined by set of equations and matrix is involved to find unknown variable. This computation is called Implicit method. This is stable method so stability criteria is not required. In genereal we guess the solution to find aproximate values and find error .Then those error is used to find new overrite data same process follows till we find the exact value or value of error is under tolerance.

In order to solve dependant variable iterative solvers are used. Those are:-

1. Jacobi iterative solver
   In this we assume the initial values at all nodes to calculate the new values, which will update once we find all nodes values and process will go on untill we arrive at converge solution.But for this method we need to make equations into diagonal dominant matrix.
   
   
2. Guass siedel solver
   In this iteration we assume the values obtianed in jacobi iteration for the next equation immidiately without waiting for the commpletion of first cycle.
   
   
3. Succesive over relaxation solver
   This method is an updated version of previous two .In addition we use the relaxtion factor to ease the process. For tis analysis we use the relaxtion factor to be greater than 1 and less than 2.

Transient state of eqation in 2D

`del T/delt = alpha*(del^2T/del x + del^2T/del y)`

solving 2D transient heat conduction equation explicitly we get

`T_P^(n+1) = (1-2*K_1-2*K_2)*T_P^(n)+K_1(T_L+T_R)^(n)+K_2(T_T+T_B)^(n)`

`K_1 = alpha*(trianglet)/(trianglex^2)`

`K_2 = alpha*(trianglet)/(triangley^2)`

Solving for implicitly we get

`T_P^(n+1) = 1/K*[T_P^(n)+K_1(T_L+T_R)^(n+1)+K_2(T_T+T_B)^(n+1) ]`

`K=1+2*K_1+2*K_2`

`K_1 = alpha*(trianglet)/(trianglex^2)`

`K_2 = alpha*(trianglet)/(triangley^2)`

For steady state the temperarure is not changing hence

Steady state equation

`(del^2T/del x + del^2T/del y) =0`

For solving the steady state 2D using different iterative techniques

`(T_L - 2T +T_R)/(Delta x^2) + (T_B- 2T +T_T)/(Delta y^2) =0`

or simplyfing,

assume ,`k = 2(Delta x^2 +Delta y^2)/(Delta x^2*Delta y^2`

we can say that

`T_p = 1/k((T_L +T_R)/(Delta x^2) + (T_B +T_T)/(Delta y^2))`

Boundary conditions

We select boundary of independent variables of 10 X 10 matrix in which we defing given temperature values for bottom, top , left & right. and we guess initial value of 100K.

Main Code

close all
clear all
clc

%% defining space domain
nx=10;
ny=nx;

% specifying length of domain
L=1;
% defining space domain
x= linspace(0,L,nx);
y= linspace(0,L,ny);

%grid size
dx= x(2) - x(1);
dy= y(2) - y(1);

% defining temperature using guess temperature 200K
T= 100*ones(nx,ny);

% Initia boundary conditions at left and right
T(1:end,1) = 400;
T(1:end,end) = 800;
% Initia boundary conditions at bottom and top
T(end,1:end) = 900;
T(1,1:end) = 600;
%average values of T at corner boundaries
T(1,1) = 500;
T(end,1) = 650;
T(1,end) = 700;
T(end,end) = 850;

% Tolerence
tol= 1e-4;

% copy of temperature
told=T;
t_initial=T;

%Coefficient obtained from the derivation for steady state equation
k = 2*((dx^2)+(dy^2))/((dx^2)*(dy^2));

%diffusivity co-efficient in the  2D heat conduction equation
alpha = 1.4;

% Calculation of time step
dt = (0.5)/((alpha/(dx^2))+(alpha/(dy^2)));

% Co-efficient  multipied for x diffuse term
k1 = dt*alpha/dx^2;
% Co-efficient  multipied for y diffuse term
k2 = dt*alpha/dy^2;

%relaxation factor for SOR for steady state heat conduction
a = 1.5;
%relaxation factor for SOR for transiet heat conduction
a1 = 1.15;


%Function to calculate the temperature of Steady state using jacobi iterative solver
[error_j,itr_jacobi] = jacobi(nx,ny,T,told,k,dx,dy,tol);

%Function to calculate the temperature of Steady state using Gauss Seidal iterative solver    
[error_gs,itr_gs] = gs(nx,ny,T,told,k,dx,dy,tol);

%Function to calculate the temperature Steady state using Successive over relaxation iterative solver    
[error_sor_gs,itr_sor] = sor_gs(nx,ny,T,told,k,dx,dy,tol,a); 

%Function to calculate the temperature of Transient heat conduction explicitly    
[error_dt,iter_dt]=transient_heat_conduction(nx,ny,k1,k2,told,T,tol,dt);   

%Function to calculate the temperature of Unsteady state using jacobi iterative solver
[error_jac,iter_jac]=transient_heat_conduction_jacobi(nx,ny,k1,k2,told,T,tol,t_initial,dt);

%Function to calculate the temperature of Unsteady state using Gauss Seidal iterative solver
[error_gs_imp,iter_gs_imp]=transient_heat_conduction_gs(nx,ny,k1,k2,told,T,tol,t_initial,dt);

%Function to calculate the temperature of Unsteady state using Successive over relaxation iterative solver
[error_sor_gs_imp,iter_sor_gs_imp]=transient_heat_conduction_sor_gs(nx,ny,k1,k2,told,T,tol,t_initial,a1,dt);

Explanation of Main Code

1. Specifying the no of grid points in the spacial domain x & y In this case 10 grid points in x& y direction.
2. Specifying of length of domain in x & y in this case equals to 1.
3. Defining spacial domain for the length and grid points specified.
4. Calculating the grid spacing in x & y space domain .
5. Temperaure guess value as 100K fo no. of grid points defined.
6. Applying boundary condition as mentioned in problems.
7. The accuracy of the solution is measured with computed approximated solution and compared with the tolerance value so that the converged solution is in the desired accuracy.
8. A copy of initial temperature matrix is made in told & t_initial which is used while calculating the temperature at n+1 time steps explicitly and implicitly.
9. Calculating `K,K_1 & K_2`. here `alpha`is taken 1.4.
10. Function is called to calculate temperature value at n+1 time step for both steady and unsteady using differnt iterative solvers and explicitly method is also used for calculating temperature of unsteady heat conduction equation.

Explanation of FUNCTIONS in Main Code

• In functions error dummy is assigned to error variable so the while loop condition would satisfy .
• Inside while there is 2 for loop so we calculate data for all x and y directions and absolute value of error is calculated.
• At end of for loop we assign told to t so for next iteration the value is updated.
• Iterative value is assigned to its value +1 so the no. of iterarions used can be found.
• "T" stores the temperature data using formulas we mention for transient and steady state respectively.
• After all while loop and for loop end we plot contour for desired data plot.
• In Transient state the time marching step is added nt= T/dt . It is done by adding for loop before the while loop and after end of while loop assign told =T so the final value of 1st iteration is used for next time step.

FUNCTION CODES

Steady state jacobi method

function [error_j,itr_jacobi] = jacobi(nx,ny,T,told,k,dx,dy,tol)

% defining iteration
itr_jacobi=0;
% dummy error value 
error_j= 900;

tic
%convergence loop
while (error_j > tol)
    %iteration loop
    for i= 2:nx-1
        
        for j= 2:ny-1
            % formula to calculate the temperature
                TL = told(i-1,j);
                TR = told(i+1,j);
                TB = told(i,j-1);
                TT = told(i,j+1);
                
                H = (TL +TR)/(dx^2);
                V = (TB +TT)/(dy^2);
                T(i,j)= (H + V)* 1/k;
        end
    end
    error_j= max(max(abs(told - T)));
    told= T;
    % iteration updation 
    itr_jacobi= itr_jacobi+1;
end
time_jacobi=toc;
figure(1)
[C,h]=contourf(T);
%Labelling the plot here C reperesents the temperature values and h rotates the temperature values as it moves.
 clabel(C,h);
 axis ij;
title('number of iterations for Jacobi method for steady state heat conduction',itr_jacobi);
end
          

Steady state Guass siedel method

function [error_gs,itr_gs] = gs(nx,ny,T,told,k,dx,dy,tol)

% defining iteration
itr_gs=0;
% dummy error value 
error_gs= 9e9;

%convergence loop
while (error_gs > tol)
    %iteration loop
    for i= 2:(nx-1)
        for j= 2:(ny-1)
            % formula to calculate the temperature
                TL = T(i-1,j);
                TR = told(i+1,j);
                TB = T(i,j-1);
                TT = told(i,j+1);
                
                H= (TL +TR)/dx^2;
                V = (TB +TT)/dy^2;
                T(i,j)= (1/k)*(H + V);
        end
    end
    error_gs= max(max(abs(told - T)));
    told=T;
    % iteration updation 
    itr_gs= itr_gs+1;
end
 % Plot representing the converged results.
figure(2)
[C,h]=contourf(T);
%Labelling the plot here C reperesents the temperature values and h rotates the temperature values as it moves.
clabel(C,h);
axis ij;
title('Number of iterations for Gauss seidal method steady state heat conduction',itr_gs)
end

Steady state Sor method

function [error_sor_gs,itr_sor] = sor_gs(nx,ny,T,told,k,dx,dy,tol,a)
% defining iteration
itr_sor=0;
% dummy error value 
error_sor_gs= 9e9;

%convergence loop
while (error_sor_gs > tol)
    %iteration loop
    for i= 2:(nx-1)
        for j= 2:(ny-1)
            % formula to calculate the temperature
                TL = T(i-1,j);
                TR = told(i+1,j);
                TB = T(i,j-1);
                TT = told(i,j+1);
                
                H= (TL +TR)/dx^2;
                V = (TB +TT)/dy^2;
                T(i,j)= (told(i,j))*(1-a) + a*((H + V)/k);
               
        end
    end
    error_sor_gs= max(max(abs(told - T)));
    told=T;
    % iteration updation 
    itr_sor= itr_sor+1;
end

 % Plot representing the converged results.
  figure(3)
  [C,h]=contourf(T);
  %Labelling the plot here C reperesents the temperature values and h rotates the temperature values as it moves.
  clabel(C,h);
  axis ij;
  title('number of iteration for Successive over relaxation method for steady state heat conduction', itr_sor)
end

Unsteady state Explicit method

Difference form above is time marching loop before while for n+1 time step.

function [error_dt,iter_dt]=transient_heat_conduction(nx,ny,k1,k2,told,T,tol,dt)

    %dummy error
    error_dt = 9e9;
    %defining iteration variable to store and display the number of
   %iterations computed for the values to converge.
    iter_dt =0;
    %convergence loop
    %time marching
    nt= T/dt;
    tic
    for m= 1:nt
    while (error_dt>tol)
            
            %Iteration loop
     for i = 2:nx-1
         for j = 2:ny-1
             
              H = told(i-1,j)-2*told(i,j)+told(i+1,j);%<-- This line represent temperature values in x direction
              V = told(i,j-1)-2*told(i,j)+told(i,j+1);%<-- This line represent temperature values in y direction
                %Formula to calculate the temperature at each grid
                    %point
                    T(i,j)= told(i,j)+k1*H+k2*V;      
                    
                end
            end

           %error for each iteration
        
            error_dt = max(max(abs(told-T)));
            
            %updating the temperature values for the next iteration. 
            told = T;
            
           %incrementing iter variable by 1 to display the no. of
           %iterations taken for the solution to converge.
           iter_dt = iter_dt+1;
           


    end 
    told=T;
    end
    time1=toc;
           %Plot representing the converged results.
           figure(4)       
           [C,h]=contourf(T);
           %Labelling the plot here C reperesents the temperature values and h rotates the temperature values as it moves.
           clabel(C,h);
           axis ij
           %title('Number of iterations using transient heat conduction explicitly',iter_dt);
           title_text1 = sprintf('Transient state Explicit using Jacobi = %d', iter_dt);
        title_text2 = sprintf('Simulation time  +%f s', time1);
        title_text3 = sprintf('Final error = %f' , error_dt);
        title({title_text1; title_text2;title_text3});
end

FUnction for Unsteady state using jacobi method.

As we are now solveing for transient state its achieved by certain added timesteps.

1st while loop is used to check for steady state solution by finding error between temperature succesive time steps comapring with tolerancing values . @nd loop while finds the converged value of temperature for particular time step

Before while loop "for" loop for time marching is calculated nt time steps. at end we assign told to T so for new time step updated value is used.
Tic & toe is used for calculating iteration solving time.

unction [error_jac,iter_jac]=transient_heat_conduction_jacobi(nx,ny,k1,k2,told,T,tol,t_initial,dt)

    % dummy error
    error_jac=9e9;
    error_f = 9e9;
   %defining iteration variable to store and display the number of
   %iterations computed for the values to converge.
    iter_jac=0;
     %time marching
    nt= T/dt;
    tic
    for m= 1:nt
    %While loop to stop iterations when the temperature values converge
    %whle computing for different time steps.
    while (error_f>tol)
        
        %Convergence loop each time step
        while (error_jac>tol)

                %Iteration loop
                for i = 2:nx-1

                    for j = 2:ny-1
                        
                        %Co-efficient obtained while deriving for transient heat conduction equation
                        term1= (1+(2*k1)+(2*k2))^-1;
                        %k1 and k2 represent the cfl value for x and y
                        %domain repectively.
                        H = told(i-1,j)+told(i+1,j);%<-- This line represent temperature values in x direction
                        V = told(i,j-1)+told(i,j+1);%<-- This line represent temperature values in y direction
                        %Formula to calculate the temperature at each grid
                        %point
                        T(i,j)= (t_initial(i,j)*term1)+(term1*k1*H)+(term1*k2*V);


                    end
                end

                %error for each iteration
              
                error_jac = max(max(abs(told-T)));
                %updating the temperature values for the next iteration. 
                told = T;
                
                %incrementing iter variable by 1 to display the no. of
                %iterations taken for the solution to converge.
                iter_jac = iter_jac+1;

 
        end 
      
      % error defined to obtain desired accuracy for the convergenc rate.
      error_f = max(max(abs(t_initial-T)));
      t_initial = T; 
      %Reseting the error value for the covergence loop to compute for the
      %next time step.
      error_jac = 9e9;

    end
    told=T;
    end
    time2=toc;
                %Plot representing the converged results.
                figure(5)       
               [C,h]=contourf(T);
               %Labelling the plot here C reperesents the temperature values and h rotates the temperature values as it moves.
               clabel(C,h)
               axis ij
               %title('Number of iterations for implicit Jacobi',iter_jac);
               title_text1 = sprintf('Transient state Implicit using Jacobi = %d', iter_jac);
        title_text2 = sprintf('Simulation time  +%f s', time2);
        title_text3 = sprintf('Final error = %f' , error_f);
        title({title_text1; title_text2;title_text3});

end

Transient Guass siedel method

function [error_gs_imp,iter_gs_imp]=transient_heat_conduction_gs(nx,ny,k1,k2,told,T,tol,t_initial,dt)

    %Definition of dummy error
    error_f = 9e9;
    error_gs_imp =9e9;
   %defining iteration variable to store and display the number of
   %iterations computed for the values to converge.
    iter_gs_imp=0;
     %time marching
    nt= T/dt;
    tic
    for m= 1:nt
    %While loop to stop iterations when the temperature values converge
    %whle computing for different time steps.
    while(error_f >tol)
        
        %Convergence loop each time step
        while (error_gs_imp>tol)

                %Iteration loop
                for i = 2:nx-1

                    for j = 2:ny-1
                        
                        %Co-efficient obtained while deriving for transient heat conduction equation
                        term1= (1+(2*k1)+(2*k2))^-1;
                        %k1 and k2 represent the cfl value for x and y
                        %domain repectively.
                        H = T(i-1,j)+told(i+1,j);%<-- This line represent temperature values in x direction
                        V = T(i,j-1)+told(i,j+1);%<-- This line represent temperature values in y direction
                        %Formula to calculate the temperature at each grid
                        %point
                        T(i,j)= (t_initial(i,j)*term1)+(term1*k1*H)+(term1*k2*V);


                    end
                end

                %error for each iteration
               
                error_gs_imp =max( max(abs(told-T)));
                
                %updating the temperature values for the next iteration.  
                told = T;
                
               %incrementing iter variable by 1 to display the no. of
               %iterations taken for the solution to converge.
               iter_gs_imp = iter_gs_imp+1;

 
        end 
      
      % error defined to obtain desired accuracy for the convergenc rate.
      error_f = max(max(abs(t_initial-T)));
      %Updating the previous temperature values to computed values.
      t_initial = T; 
      %Reseting the error value for the covergence loop to compute for the
      %next time step.
      error_gs_imp = 9e9;

    end
    told=T;
    end
      time3=toc;
      %Plot representing the converged results.
                figure(6)       
               [C,h]=contourf(T);
               %Labelling the plot here C reperesents the temperature values and h rotates the temperature values as it moves.
               clabel(C,h)
               axis ij
               %title('Number of iterations for implicit Gauss seidal',iter_gs_imp);
                title_text1 = sprintf('Transient state Implicit using Gauss SIEDEL = %d', iter_gs_imp);
        title_text2 = sprintf('Simulation time  +%f s', time3);
        title_text3 = sprintf('Final error = %f' , error_f);
        title({title_text1; title_text2;title_text3});


end

Transient SOR method

function [error_sor_gs_imp,iter_sor_gs_imp]=transient_heat_conduction_sor_gs(nx,ny,k1,k2,told,T,tol,t_initial,a1,dt)
        
        %Definition of dummy error
        error_f = 9e9;
        error_sor_gs_imp = 9e9;
        %defining iteration variable to store and display the number of
         %iterations computed for the values to converge.
        iter_sor_gs_imp =0;
        %time marching
    nt= T/dt;
    tic
    for m= 1:nt
    %While loop to stop iterations when the temperature values converge
    %whle computing for different time steps.
    while (error_f>tol)
        
        %Convergence loop each time step
        
        while (error_sor_gs_imp>tol)

                %Iteration loop
                for i = 2:nx-1

                    for j = 2:ny-1
                        %Co-efficient obtained while deriving for transient heat conduction equation
                        term1= (1+(2*k1)+(2*k2))^-1;
                        H = T(i-1,j)+told(i+1,j);%<-- This line represent temperature values in x direction
                        V = T(i,j-1)+told(i,j+1);%<-- This line represent temperature values in y direction
                        %Formula to calculate the temperature at each grid
                        %point
                        T(i,j)= (t_initial(i,j)*term1)+(term1*k1*H)+(term1*k2*V);
                        %Successive over relaxation Formula to calculate the
                        %temperature at each grid point.
                        T(i,j)= (told(i,j))*(1-a1)+ a1*T(i,j);
                        


                    end
                end
                
               %error for each iteration
               
               error_sor_gs_imp = max(abs(told - T));
               
               %updating the temperature values for the next iteration.
               told = T;
               
                %incrementing iter variable by 1 to display the no. of
               %iterations taken for the solution to converge.
               iter_sor_gs_imp = iter_sor_gs_imp+1;

 
        end 
        
      %error defined to obtain desired accuracy for the convergenc rate.
      error_f = max(max(abs(t_initial-T)));
      %Updating the previous temperature values to computed values.
      t_initial = T;
      %Reseting the error value for the covergence loop to compute for the
      %next time step.
      error_sor_gs_imp = 9e9;
      
    end
    told=T;
    end
    time4=toc;
               %Plot representing the converged results.
               figure(7)       
               [C,h]=contourf(T);
               %Labelling the plot here C reperesents the temperature values and h rotates the temperature values as it moves.
               clabel(C,h)
               axis ij
               %title('Number of iteration for implicit SOR',iter_sor_gs_imp);
               title_text1 = sprintf('Transient state Implicit using SOR = %d', iter_sor_gs_imp);
        title_text2 = sprintf('Simulation time  +%f s', time4);
        title_text3 = sprintf('Final error = %f' , error_f);
        title({title_text1; title_text2;title_text3});

               

end

OUTPUT GRAPH

CONCLUSION

Implicit method can be used to problems when converged results are needed in less computational time as they are unconditionaly stable, With increase in timestep size the convergence rate is increased which leads to steady state with less. nno of iterations.

Explicit method results time accurate method as it is governed by stability critereia when not satisfied it leads to unstable solution and require more iterations to attain steady state.

But code provided above time step size remains constant so its stable criteria .Hence the explicit leads to faster convergence rate and as the time step is small no. of implicit scheme iterations increaes which shows slower convergence rate.