Solving the steady and unsteady 2D heat conduction problem

1. Introduction:

1.1 Steady State:

A system or process is said to be steady when the state variables don't change along the time.

`(∂p)/(∂t) = 0`

where p is any property of the system.

In the steady-state heat conduction problem, the temperature (T) doesn't change with respect to time.

`(∂T)/(∂t) = 0`

Thus, the equation for 2D study-state heat conduction becomes

`(∂^2T)/(∂x^2) +(∂^2T)/(∂y^2)= 0` ..........(1)

Note: Here the spatial derivative is elliptical, thus the information for computation is taken from all the sides. Thus Central Differencing Method is used.

Discretizing the equation 1,

`T_p = (1/k)(((T_L +T_R)/dx^2)+((T_T +T_B)/dy^2))` .........(2)

where, `k = 2((dx^2+dy^2)/(dx^2dy^2))`

and subscript represent the left(L),right(R), top(T) and bottom(B) the points location with respect to present point (p).

1.2 Unsteady State:

A system or process is said to be in unsteady-state when the state variables change with respect to time.

Thus,

`(∂p)/(∂t) != 0`

where p is any property of the system.

In the steady-state heat conduction problem, the temperature (T) doesn't change with respect to time.

`(∂T)/(∂t) != 0`

Thus, the equation for 2D study-state heat conduction becomes

`(∂T)/(∂t)+alpha((∂^2T)/(∂x^2) +(∂^2T)/(∂y^2))= 0` ............(3)

Note: Here the spatial derivative is elliptical form, the information for computation is taken from all the sides. Thus we will be using Central Differencing Method.

Explicit Form

Discretizing the equation 3 and forming an explicit equation,

`T_p^(n+1) =T_p^(n)+[(k_1)((T_L -2T_p+T_R)/dx^2)+(k_2)((T_T - 2T_p +T_B)/dy^2)]^n` .................(4)

where `k_1 = (alpha*dt)/dx^2` and `k_1 = (alpha*dt)/dy^2`

and n = represents time(t)

Convergence condition by Courant–Friedrichs–Lewy

The convergence condition by Courant–Friedrichs–Lewy for this condition is,

if `CFL > 0.5`then the solution becomes unstable.

then`CFL le 0.5`

implies `CFL_x +CFL_y le 0.5`

`k_1 +k_2 le 0.5`

Implicit Form

Discretizing the equation 3 and forming an implicit equation,

`T_p^(n+1) = (1/(1+2k_1 +2k_2))(T_p^(n) + k_1(T_L+T_R)^(n+1) +k_2(T_T+T_B)^(n+1))` ..................(5`

1.5 Iterative method:

It is a mathematical procedure where the initial values are guessed to generate a sequence of improving approximate solutions, thus the `n^(th)` approximation is derived from the previous approximations.

We have 3 iterative methods:

1. Jacobi method
2. Gauss-Seidel method
3. The Successive over-relaxation method

1.5.1 Jacobi Method

Here, considering a system of linear equation,

`a_11*x_1 + a_12*x_2 + a_13*x_3 =b_1`

`a_21*x_1+ a_22*x_2 + a_23*x_3 =b_2`

`a_31*x_1 + a_32*x_2 + a_33*x_3 =b_3`

In matrix form

`[[a_11, a_12, a_13],[a_21, a_22, a_23],[a_31, a_32, a_33]][[x_1],[x_2],[x_3]] = [[b_1],[b_2],[b_3]]`

i.e [A][X] = [B]

Let `i` be the number of iteration

For the first iteration the values of `[[x_1],[x_2],[x_3]]` is assumed.

For the next iteration(`i+1`),

`x_1^(i+1) = 1/(a_11)(b_1 - (a_12*x_2^i) - (a_13*x_3^i))`

`x_2^(i+1) = 1/(a_22)(b_2 - (a_21*x_1^i) - (a_23*x_3^i))`

`x_3^(i+1) = 1/(a_33)(b_3 - (a_31*x_1^i) - (a_32*x_2^i))`

So, the next iterations are computed using only the previous iteration value of x.

1.5.2 Gauss-Seidel Method

Considering the previous system of the linear equation and assuming the initial values of x.

Next iteration (`i+1`) is computed as follows

`x_1^(i+1) = 1/(a_11)(b_1 - (a_12*x_2^i) - (a_13*x_3^i))`

here the `x_1` value of this iteration is used for the computation of the next `x_2 and x_3` values

`x_2^(i+1) = 1/(a_22)(b_2 - (a_21*x_1^(i+1)) - (a_23*x_3^i))`

here the `x_2` value of this iteration is used for the computation of the next `x_3` values

`x_3^(i+1) = 1/(a_33)(b_3 - (a_31*x_1^(i+1)) - (a_32*x_2^(i+1)))`

This reduces the number of iterations significantly compared to Jacobi method

1.5.3 The Successive over-relaxation method

It is a variant of the Gauss-Seidel method for solving linear equations resulting in faster convergence. Considering the same linear equation system

[A][X] = [B].

The Gauss-Seidel method is used for iteration but after 1 iteration, equation(6) is used and the values of this equation are used for the next iteration. This goes on till it converges.

The equation is `x_(1)^(i+1) =omega*x_1^i +(1-omega)(x_1^((i+1) of GS))` ...............(6)

The same is computed for `x_2 and x_3`.

If `omega > 1` it is called Over-Relaxation.

If `omega < 1` it is called under-Relaxation.

`omega` value should not be very high as it can overshoot the value and increases the number of iterations.

1.6 Problem Statement:

Simulating the 2D heat conduction for both steady and unsteady state using the boundary condition

Top Boundary = 600 K

Bottom Boundary = 900 K

Left Boundary = 400 K

Right Boundary = 800 K

over the domain `1 m^2`.

1.7 Objective

1. Simulate the steady-state 2D heat conduction for the given domain using iterative solvers. And also analyze the number of iterations for different solvers.
2. Simulate the unsteady-state 2D heat conduction for the given domain in both explicit and implicit form.
3. Analyze how CFL number affects the stability and finding the optimum dt and dx at `alpha = 1.4` (Thermal diffusivity)
4. Using iterative solvers for implicit equation and also analyze the number of iterations for different solvers.
5. For the SOR method, omega value is taken by analyzing both steady and unsteady state for giving a minimum number of iteration and the solution being accurate.

2. Preparation of the domain using the Boundary conditions, Initial conditions and few assumptions

Boundary Conditions:

Top Boundary = 600 K

Bottom Boundary = 900 K

Left Boundary = 400 K

Right Boundary = 800 K

The value of corner points is equal to the average of the 2 BC which affects that point.

Initial Conditions:

Temperature = 500k (the closer the value to the solution lesser the iterations)

Domain:

over the domain `1 m^2`.

Assumptions:

1. The mesh at x and y direction are equal. i.e `n_x=n_y = 10`
2. Thermal diffusivity (`alpha`) = 1.4
3. Tolerance = 1e-5 (for convergence)

2.1 CODE 1

clear all
close all
clc

% creating the domain
Lx = 1;
Ly = 1;
nx=ny=10;
x=linspace(0,Lx,nx);
y=linspace(0,Lx,ny);
dx=Lx/(nx-1);
dy=Ly/(ny-1);

%initial and boundary conditions for the domain is prepared
T = ones(ny,nx);

%Boundary condition
T(1,:) = 600;
T(ny,:) = 900;
T(:,1) = 400;
T(:,nx) = 800;

%boundary conditions of the corners
T(1,1) = (T(1,2)+T(2,1))/2;
T(1,nx) = (T(1,(nx-1))+T(2,nx))/2;
T(ny,1) = (T((ny-1),1)+T(ny,2))/2;
T(ny,nx) = (T((ny-1),nx)+T(ny,(nx-1)))/2;
%Initial conditions inside the domain
for i=2:(nx-1)
  for j = 2:(ny-1)
    T(i,j) = 300;
  endfor
endfor

3. Simulate the steady-state 2D heat conduction for the given domain:

Using equation 2,

`T_p = (1/k)(((T_L +T_R)/dx^2)+((T_T +T_B)/dy^2))`

Different iterative solvers are used to solve the equation.

3.1. Jacobi method:

Code 2: Creating a Jacobi method function

% Creating a function to execute jacobi method 
function [T,jacobi_ite] = jacobi_iter(T,Told,nx,ny,dx,dy,error,tol)
  
  % To know the number of iterations
  jacobi_ite=0;
  
  %Convergence loop
  while error >tol
    %Sapce loop
    for i=2:(nx-1)
      for j = 2:(ny-1)
       k = 2*(dx^2+dy^2)/((dx^2)*(dy^2));
       T(i,j) = (((Told(i-1,j)+Told(i+1,j))/dx^2)+((Told(i,j-1)+Told(i,j+1))/dy^2))/k;
      endfor
    endfor

  jacobi_ite=jacobi_ite+1;

  % finding the max error
  error1 = max(abs(Told-T));
  error = max(error1);
  
  % Copying the T values to the old
  Told=T;
  endwhile 
endfunction

3.2 Gauss-Seidel Method

Code 3: Creating a Gauss-Seidel Method function

% Creating a function to execute GS method 
function [T,GS_ite] = GS_iter(T,Told,nx,ny,dx,dy,error,tol)
  
  % To know the number of iterations
  GS_ite=0;
  
  %Convergence loop
  while error >tol
    
    %Sapce loop
  for i=2:(nx-1)
    for j = 2:(ny-1)
      k = 2*(dx^2+dy^2)/((dx^2)*(dy^2));
      T(i,j) = (((T(i-1,j)+T(i+1,j))/dx^2)+((T(i,j-1)+T(i,j+1))/dy^2))/k;
    endfor
  endfor
  GS_ite=GS_ite+1;
  
  % finding the max error
  error1 = max(abs(Told-T));
  error = max(error1);
  
  % Copying the T values to the old
  Told=T;
  endwhile 
endfunction

3.3 The Successive over-relaxation method

Code 4: Creating a Successive over-relaxation method function

% Creating a function to execute SOR method 
function [T,SOR_ite] = SOR_iter(T,Told,nx,ny,dx,dy,error,tol,omega)
  
  % To know the number of iterations
  SOR_ite=0;
  
  %Convergence loop
  while error >tol
    
    %Sapce loop
    for i=2:(nx-1)
      for j = 2:(ny-1)
        % GS method
       k = 2*(dx^2+dy^2)/((dx^2)*(dy^2));
       T(i,j) = (((T(i-1,j)+T(i+1,j))/dx^2)+((T(i,j-1)+T(i,j+1))/dy^2))/k;
       % using equation 6
       T(i,j) = (1-omega)*Told(i,j)+omega*(T(i,j));
      endfor
    endfor
   SOR_ite=SOR_ite+1;
   
   % finding the max error
   error1 = max(abs(Told-T));
   error = max(error1);
   
   % Copying the T values to the old
   Told=T;
   endwhile 
endfunction

Code 5: Finding the optimum `omega`

clear all
close all 
clc

% creating the domain
Lx = 1;
Ly = 1;
nx=ny=10;
x=linspace(0,Lx,nx);
y=linspace(0,Lx,ny);
dx=Lx/(nx-1);
dy=Ly/(ny-1);

%initial and boundary conditions for the domain is prepared
T = ones(ny,nx);

%Boundary condition
T(1,:) = 600;
T(10,:) = 900;
T(:,1) = 400;
T(:,10) = 800;

%boundary conditions of the corners
T(1,1) = (T(1,2)+T(2,1))/2;
T(1,10) = (T(1,9)+T(2,10))/2;
T(10,1) = (T(9,1)+T(10,2))/2;
T(10,10) = (T(9,10)+T(10,9))/2;

%Initial conditions inside the domain
for i=2:9
  for j = 2:9
    T(i,j) = 300;
  endfor
endfor

%making a copy of the temperature profile
T_initial=Told=T;

% defining the Error and tolerance 
error = 9e9;
tol = 1e-5; 

% Range of omega
n = 10
omega=linspace(1,1.9,n);

% Calculating SOR iterations
for h = 1:10
 error = 9e9;
[SOR_ite(h)] = SOR_iter(T,Told,nx,ny,dx,dy,error,tol,omega(h));
end

% Ploting omega vs SOR iterantions 
plot (omega,SOR_ite,'x')
xlabel("omega")
ylabel("Number of iterations")

Note:

Make sure you pass the right variables and receive the arguments correctly for the function (when you call the function)

3.4 Main Code

3.4.1 CODE 6

clear all
close all
clc

% creating the domain
Lx = 1;
Ly = 1;
nx=ny=10;
x=linspace(0,Lx,nx);
y=linspace(0,Lx,ny);
dx=Lx/(nx-1);
dy=Ly/(ny-1);

%initial and boundary conditions for the domain is prepared
T = ones(ny,nx);

%Boundary condition
T(1,:) = 600;
T(10,:) = 900;
T(:,1) = 400;
T(:,10) = 800;

%boundary conditions of the corners
T(1,1) = (T(1,2)+T(2,1))/2;
T(1,10) = (T(1,9)+T(2,10))/2;
T(10,1) = (T(9,1)+T(10,2))/2;
T(10,10) = (T(9,10)+T(10,9))/2;

%Initial conditions inside the domain
for i=2:9
  for j = 2:9
    T(i,j) = 300;
  endfor
endfor

%making a copy of the temperature profile
T_initial=Told=T;

% defining the Error tolerance and omega 
error = 9e9;
tol = 1e-5;
omega=1.5;

%Calling the respective user-definerd functions by inputing all the parameters and expecting the iteration number and the updated temperature profile
[T_jacobi,jacobi_ite] = jacobi_iter(T,Told,nx,ny,dx,dy,error,tol);
[T_GS,GS_ite] = GS_iter(T,Told,nx,ny,dx,dy,error,tol);
[T_SOR,SOR_ite] = SOR_iter(T,Told,nx,ny,dx,dy,error,tol,omega);

%Ploting the temperature contour
figure (1)
[c,h]=contourf(x,y,T_jacobi,10)
clabel(c,h)
colorbar
colormap(jet)
title("Jacobi Method")
set(gca,'YDIR','reverse')

figure (2)
[c,h]=contourf(x,y,T_GS,10)
clabel(c,h)
colorbar
colormap(jet)
title("Gauss-Seidel Method")
set(gca,'YDIR','reverse')

figure (3)
[c,h]=contourf(x,y,T_SOR,10)
clabel(c,h)
colorbar
colormap(jet)
title("Successive over-relaxation method")
set(gca,'YDIR','reverse')


% creating the bar graph to showcase the number of iterrations of different solvers
figure(4)
y = [jacobi_ite 0 0
0 GS_ite 0
0 0 SOR_ite];
h = bar (y);
set (h, "ydata");
set (h(1), "facecolor", "r");
title("Variations of iterations for diiferent solvers")
legend("Jacobi method","Gauss-Seidel Method","Successive over-relaxation method")

4. Simulate the unsteady-state 2D heat conduction for the given domain

4.1 Explicit Form

From equation (4) ,

`T_p^(n+1) =T_p^(n)+[(k_1)((T_L -2T_p+T_R)/dx^2)+(k_2)((T_T - 2T_p +T_B)/dy^2)]^n`

where `k_1 = (alpha*dt)/dx^2` and `k_1 = (alpha*dt)/dy^2`

and n = represents time(t)

4.1.1 Calculating dt

We know that ,`alpha = 1.4 and dx=dy=10`

From CFL condition,

`CFL_x +CFL_y le 0.5`

i.e `k_1 +k_2 le 0.5`

As `n_x = n_y` and we have a unit sq domain so `L_x = L_y`

then `k_1=k_2=k=alpha*dt/dx^2`

imples, `2k = 0.5` consider the max of CFL possible.

`k=0.25`

`alpha*dt/dx^2=0.25`

`1.4*dt/(0.111)^2=0.25`

`dt ~= 0.002`

4.1.2 Main Code 7:

clear all
close all
clc

% creating the domain
Lx = 1;
Ly = 1;
nx=ny=10;
x=linspace(0,Lx,nx);
y=linspace(0,Lx,ny);
dx=Lx/(nx-1);
dy=Ly/(ny-1);

%initial and boundary conditions for the domain is prepared
T = ones(ny,nx);

%Boundary condition
T(1,:) = 600;
T(10,:) = 900;
T(:,1) = 400;
T(:,10) = 800;

%boundary conditions of the corners
T(1,1) = (T(1,2)+T(2,1))/2;
T(1,10) = (T(1,9)+T(2,10))/2;
T(10,1) = (T(9,1)+T(10,2))/2;
T(10,10) = (T(9,10)+T(10,9))/2;

%Initial conditions inside the domain
for i=2:9
  for j = 2:9
    T(i,j) = 300;
  endfor
endfor

%making a copy of the temperature profile
T_initial=Told=T;

% defining the Error and tolerance 
alpha=1.4;
dt=0.002;
k = (alpha*dt)/dx^2

% Time loop
for nt = 1:1400
  
  %Space tool
  for i=2:(nx-1)
    for j=2:(ny-1)
      % explicit equation
      T(i,j) = Told(i,j) + (k*(Told(i-1,j) - 2*Told(i,j) +Told(i+1,j)))+ k*(Told(i,j-1)-2*Told(i,j)+T(i,j+1));
    endfor
  endfor
  
  % equating old to new values
  Told = T;
endfor

figure (1)
[c,h]=contourf(x,y,T,10)
clabel(c,h)
title("Explicit Form")
colorbar
colormap(jet)
set(gca,'YDIR','reverse')

4.2 Implicit Form

Using the equation (5)

`T_p^(n+1) = (1/(1+2k_1 +2k_2))(T_p^(n) + k_1(T_L+T_R)^(n+1) +k_2(T_T+T_B)^(n+1))`

where `k_1 = (alpha*dt)/dx^2` and `k_1 = (alpha*dt)/dy^2`

and n = represents time(t)

We know that ,`alpha = 1.4 and dx=dy=10` and through calculation dt = 0.002

We will be using iterative solvers and compare the number of iteration to achieve the result.

4.2.1 Jacobi method:

Code 8: Creating a Jacobi method function

%Creating jacobi method
function [T, Jacobi_iter] = jacobi_unsteady(T,Told,Tprev,k,nx,ny,dx,dy,error,tol)
  
  %Time loop
  for nt = 1:1400
    Jacobi_iter = 0;
    
    % convergence loop
    while error>tol 
      
      %space loop
     for i = 2:(nx-1)
        for j = 2:(ny-1)
          
          % jacobi method
         term_1 = 1/(1+(4*k));
          term_2 = k*term_1;
          term_3 = k*term_1;
          T(i,j)= term_1*(Tprev(i,j))+term_2*(Told(i-1,j)+Told(i+1,j))+term_3*(Told(i,j-1)+Told(i,j+1));
        end
      end
      
      % calculating the error
      error1 = max(abs(Told-T));
      error = max(error1);
      
      % updating the old values
      Told =T;
      Jacobi_iter = Jacobi_iter+1;
    end
    
    % updating the previous time value
     Tprev =T;
     
     %reseting the error
     error = 9e9;
  end
end

4.2.2 Gauss-Seidel Method

Code 9: Creating a Gauss-Seidel method function

%Creating GS method
function [T, GS_iter] = GS_unsteady(T,Told,Tprev,k,nx,ny,dx,dy,error,tol)
  
  %Time loop
  for nt = 1:1400
    GS_iter(nt) = 0;
    
    % convergence loop
    while error>tol 
      
      %space loop
     for i = 2:(nx-1)
        for j = 2:(ny-1)
          
          % GS method
         term_1 = 1/(1+(4*k));
          term_2 = k*term_1;
          term_3 = k*term_1;
          T(i,j)= term_1*(Tprev(i,j))+term_2*(T(i-1,j)+T(i+1,j))+term_3*(T(i,j-1)+T(i,j+1));
        end
      end
      
      % calculating the error
      error1 = max(abs(Told-T));
      error = max(error1);
      
      % updating the old values
      Told =T;
      GS_iter(nt) = GS_iter(nt)+1;
    end
    
    % updating the previous time value
     Tprev =T;
     
     %reseting the error
     error = 9e9;
  end
end

4.2.3 The Successive over-relaxation method

Code 10: Creating a Successive over-relaxation method function

%Creating GS method
function [T,SOR_iter] = SOR_unsteady(T,Told,Tprev,k,nx,ny,dx,dy,error,tol,omega)
  
  %Time loop
  for nt = 1:1400
    SOR_iter(nt) = 0;
    
    % convergence loop
    while error>tol 
      
      %space loop
     for i = 2:(nx-1)
        for j = 2:(ny-1)
          % GS method
          term_1 = 1/(1+(4*k));
          term_2 = k*term_1;
          term_3 = k*term_1;
          T(i,j)= term_1*(Tprev(i,j))+term_2*(T(i-1,j)+T(i+1,j))+term_3*(T(i,j-1)+T(i,j+1));
          
          % SOR method
          T(i,j) = (1-omega)*Told(i,j)+omega*(T(i,j));
        end
      end
      
      % calculating the error
      error1 = max(abs(Told-T));
      error = max(error1);
      
      % updating the old values
      Told =T;
      SOR_iter(nt) = SOR_iter(nt)+1;
    end
    
     % updating the previous time value
     Tprev =T;
     
     %reseting the error
     error = 9e9;
  end
  %used for finding omega value can bew ignored for simulating but pass the value if finding omega
  ans = max(SOR_iter);
end

Code 11: Finding optimum `omega` value

clear all
close all
clc
% creating the domain
Lx = 1;
Ly = 1;
nx=ny=10;
x=linspace(0,Lx,nx);
y=linspace(0,Lx,ny);
dx=Lx/(nx-1);
dy=Ly/(ny-1);
%initial and boundary conditions for the domain is prepared
T = ones(ny,nx);
%Boundary condition
T(1,:) = 600;
T(10,:) = 900;
T(:,1) = 400;
T(:,10) = 800;
%boundary conditions of the corners
T(1,1) = (T(1,2)+T(2,1))/2;
T(1,10) = (T(1,9)+T(2,10))/2;
T(10,1) = (T(9,1)+T(10,2))/2;
T(10,10) = (T(9,10)+T(10,9))/2;
%Initial conditions inside the domain
for i=2:9
  for j = 2:9
    T(i,j) = 500;
  endfor
endfor
%making a copy of the temperature profile
T_initial=Told=Tprev=T;
alpha=1.4;
dt=0.002;
k = (alpha*dt)/dx^2
error = 9e9;
tol = 1e-5;
omega = linspace(1,1.9,10);

for h=1:10
[SOR_ite(h)] = SOR_unsteady(T,Told,Tprev,k,nx,ny,dx,dy,error,tol,omega(h));
end
figure(1)
plot(omega,SOR_ite,"x")
xlabel("omega")
ylabel("Number of SOR iterations")

4.3 Main Code

Code 12:

clear all
close all
clc

% creating the domain
Lx = 1;
Ly = 1;
nx=ny=10;
x=linspace(0,Lx,nx);
y=linspace(0,Lx,ny);
dx=Lx/(nx-1);
dy=Ly/(ny-1);

%initial and boundary conditions for the domain is prepared
T = ones(ny,nx);

%Boundary condition
T(1,:) = 600;
T(10,:) = 900;
T(:,1) = 400;
T(:,10) = 800;

%boundary conditions of the corners
T(1,1) = (T(1,2)+T(2,1))/2;
T(1,10) = (T(1,9)+T(2,10))/2;
T(10,1) = (T(9,1)+T(10,2))/2;
T(10,10) = (T(9,10)+T(10,9))/2;

%Initial conditions inside the domain
for i=2:9
  for j = 2:9
    T(i,j) = 500;
  endfor
endfor

%making a copy of the temperature profile
T_initial=Told=Tprev=T;
alpha=1.4;
dt=0.002;
k = (alpha*dt)/dx^2
error = 9e9;
tol = 1e-5;
omega = 1.1;

%Calling the function for different solvers and the number of iteration is taken for one particular time when max iteration is visibale 
[T_jacobi,jacobi_ite] = jacobi_unsteady(T,Told,Tprev,k,nx,ny,dx,dy,error,tol);
[T_GS,GS_iter] = GS_unsteady(T,Told,Tprev,k,nx,ny,dx,dy,error,tol);
[T_SOR,SOR_ite] = SOR_unsteady(T,Told,Tprev,k,nx,ny,dx,dy,error,tol,omega);

%ploting temperature profile
figure (1)
[c,h]=contourf(x,y,T_jacobi,10,"linestyle","none")
clabel(c,h)
colorbar
colormap(jet)
set(gca,'YDIR','reverse')
pause(0.03)

figure (2)
[c,h]=contourf(x,y,T_GS,10,"linestyle","none")
clabel(c,h)
colorbar
colormap(jet)
set(gca,'YDIR','reverse')
pause(0.03)

figure (2)
[c,h]=contourf(x,y,T_SOR,10,"linestyle","none")
clabel(c,h)
colorbar
colormap(jet)
set(gca,'YDIR','reverse')
pause(0.03)

%ploting the number of iterations for different solvers
figure(4)
y = [max(jacobi_ite) 0 0
0 max(GS_iter) 0
0 0 max(SOR_ite)];
h = bar (y);
set (h, "ydata");
set (h(1), "facecolor", "r");
title("Variations of iterations for diiferent solvers")
legend("Jacobi method","Gauss-Seidel Method","Successive over-relaxation method")

5. RESULTS

CODE 1:

Fig 1: Initial temperature values

Fig 2: Countour of the initial profile

CODE 5:

Fig 3: Omega vs Number of iterations for steady-state

The above graph is drawn with omega varying from 1-1.9. The least number of iteration is found at 1.5. Thus `omega = 1.5` is used for simulation in further programs.

CODE 2,3,4,6:

Fig 4: Simulation of 2D steady-state heat conduction with different iterative solvers.

Fig 5: Number of iterations for different iterative methods -Steady-state.

The fig 5, indicates that the number of iteration is least in the Successive over-relaxation method compared to Jacobi or Gauss-Seidel method.

i.e Number of iteration for

Successive over-relaxation method

Code 7:

Fig 6: Simulation of 2D unsteady-state heat conduction using the explicit form.

Code 11:

Fig 7:Omega vs Number of iterations for unsteady-state

The fig(7), shows the least number of iterations for SOR at`omega`=1.1. Thus `omega` = 1.1 is taken in the further programs.

Code 8,9,10,12:

Fig 8: Simulation of 2D Unsteady-state heat conduction with different iterative solvers.

Fig 9: Number of iterations for different iterative methods -Unsteady-state.

The fig 9, indicates that the number of iteration is least in the Successive over-relaxation method compared to Jacobi or Gauss-Seidel method.

i.e Number of iteration for

Successive over-relaxation method

6.CONCLUSION

The steady-state is simulated using iterative methods such as Jacobi, Gauss-Seidel and Successive over-relaxation(SOR) method. For SOR from fig(3), `omega` is found when the iteration is minimum. Thus we get `omega = 1.5`. Fig(4) shows the simulation. From fig(5) we can prove that number of iterations for Successive over-relaxation method

For Unsteady state, using a CFL condition, the time step is calculated.

i.e `dt~=0.002`

Fig (8), shows the simulation for Unsteady-state using iterative methods such as Jacobi, Gauss-Seidel and Successive over-relaxation(SOR) method. For SOR from fig(7), `omega` is found when the iteration is minimum. Thus we get `omega = 1.1`. From fig(9) we can prove that number of iterations for Successive over-relaxation method

Thus for both steady and unsteady state, the number of iterations for Successive over-relaxation method is the least and `omega = 1.5 and 1.1` respectively.