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

Aim: To solve two dimensional equation for steady state and unsteady state (transient) by using explicit and implicit method with given boundary.

Using implicit method

solving using different iterative solvers i.e. .. Jacobi, Gauss seidal and successive over relaxation for both steady state and transient state heat conduction equation                                                 

 The accuracy of the solution are obtained both explicit and implicit method should be accordance with the tolerance value obtained.

Boundary Conditions:                                                                                                                                                 

  Top boundary= 600k                                                                                                                                                   

  Bottom boundary= 900 k                                                                                                                                           

  Left boundary= 400k                                                                                                                                               

  Right boundary= 800k                                                                                                                                                 

  Tolerance Value =>error=1e-4.                                                                                                                                                             

Introduction:

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 satisfied leads to stable solution. Because of the stability criteria small steps has to be taken to obtain stable solutions. With decrease in time step the number of computation increases and the computational time required for the solution to converge is less and accuracy is high 

Implicit Method:

When a computational are made by different variables are defined by set of equation and matrix is involved to find the unknown variable. This computation is called implicit method. This method is unconditionally stable as a result, no stability criteria is required to obtain a stable result and the computational time required for the solution to converge is high and accuracy of the solution is less relative to explicit method because this method uses simultaneous equation and uses guess value as the initial values to find the approximated value, If the guess value are closer to the desired value the computation time is less. As there is no stable criteria requirement larger time steps can be used.

 

In order to solve for independent variables iterative solvers are used.

Iterative solvers are mentioned as

1. Jacobi iterative solvers
2. Gauss seidal iterative solvers
3. Successive over relaxation solvers

2D Heat conduction equation: 

The time derivative term  is hyperbolic in nature which can be discretized using forward differencing scheme and the special derivative are elliptical in nature which can be discretized using central differencing scheme

Discretization of time derivative using forward differencing scheme since the time derivative in time direction which can be solved by marching method.

Discretization of special derivatives using central differencing scheme since the special derivatives are elliptical in nature, which uses information from neighboring points to obtain the solution. Below are the discretization for x and y special .derivatives

From the PDE Heat conduction Equation is:

Where n is current time step and n+1 is next step                                                                                             

 I,j+1 are neighboring points in X direction ,j,j+1 are the neighboring points in Y direction. is a grid spacing in x direction is a grid spacing in Y direction.

Iterative solver:

1. Jacobi Iterative solver:

In this method the temperature at n+1th iteration is found using old values of the temperature defined at nth iterations. Using the approximation obtained iterations procedure it continued until the desired accuracy of the solution obtained.

2. Gauss seidal iterative solver:

In this type iterative solver, temperature at n+1 iterations uses values of n+1 iterations that have already been computed .By this way old approximated values are overwritten as new approximated values for the iteration performed. Unlike Jacobi, Gauss seidal require one storage vector which can be over written with each computation this helps the solution to coverage faster than Jacobi method

3. Successive over relaxation iterative solver:

In this type of iterative solver, relaxation factor is used which helps in increasing the convergence rate. It goes by the formula.

whereis the relaxation factor, it is range is form

where   Is between 0-1 it is called under relaxation state

 This condition is used when the temperature gradient is high and due to this condition the time step does not over shoot and helps to converge at a faster rate.

where    is Where omega is Between 1-2 it is called over relaxation rate.

This condition is used when the temp gradient is low and due to this condition larger time step are used to obtain the steady state i.e., for the values to converge at faster rate.

Is the values of temperature for previous iterations.

T (GS) or T (Jacobi) is the temperature of value for n+1th iteration using Gauss seidal or Jacobi iterative solver.

Steady state 2D Heat Conduction Equation:

Steady state is where the fluid properties which does not change with respect to time.

From the 2D heat conduction equation since the system is in steady state time derivative is zero and diffusivity is constant.

Hence we have the equation for steady state as:

Discretization for steady state:

Solving we get;

Solving the 2D steady state conduction equation implicitly.                                                                                         

    Formula for Jacobi iterative solver:

Formula for Gauss seidal iterative solver:

Formula for successive over relaxation iterative solver:

Transient state:

Unsteady state is where the fluid properties change with respective time                                                                     

  The formula for 2D transient heat conduction equation is:

Discretization for unsteady state heat conduction equation:

Solving 2D transient heat conduction equation explicitly we get:

Stability criteria:

K1+K20.5 stability criteria for 2D heat conduction equation when solving explicitly.

Solving for 2D Transient Heat conduction equation implicitly we get:

 

Here n+1 is the next time step and n is the preceding time step

Formula for Jacobi iterative solver:

Formula for Gauss seidal method:

Formula for successive over relaxation method:

Program for steady state heat conduction:

% program to solve 2D steady state heat conduction using iterative solvers
% Steady state heat conduction

clear all
close all
clc

% initializing  variables
% specifyingthe length ofdomain , L=1
% Defining space domain
nx=40; % number of grid points along in x-axis
ny=40; % number of grid points along in y-axis
x=linspace(0,1,nx);
y=linspace(0,1,ny);
% grid size
dx=1/(nx-1);
dy=1/(ny-1);

% dummy error value
error=1e9;
% tolerance
tol=1e-4;
% coefficient obtained from the derivative for steady state equation
k=2*((dx^2+dy^2)/(dx^2*dy^2));

% temperature fieldnames
T=ones(nx,ny);
T(1,:)=600; % bottom
T(end,:)=900; % top
T(:,1)=400; % left
T(:,end)=800; % right


% relaxation factor for successsive over relaxation (SOR) method, for steady stateheat condition
omega=1.1;
% making acopy of the initial temperature matrix_
Told=T;

% solving 2D steady state system
[xx,yy]=meshgrid(x,y);
disp('For solver type: n=1 for jacobi; n =2 for Gauss seidal; n=3 for SOR n')
iterative_solvers = input('Enter the solver number=')
tic


% jacobi iterative method:

if iterative_solvers==1
  jacobi_iterative=1; % defining iterative variable to store and display the number of iteratations computed for the values of convergence
  
  % convergence loop
  while(error>tol)
       % Iteratation loop.
  for i=2:nx-1;
  for j=2:ny-1;
    % calculating tempat each grid points
    term1=(Told(i-1,j)+Told(i+1,j))/(dx^2);
    term2=(Told(i,j-1)+Told(i,j+1))/(dy^2);
    T(i,j)=(term1+term2)/k;
    end
    end
  error=max(max(abs(Told-T))); % 2D problem so two max
  
  % updating the temperature value for the next iteration
  Told=T;
  
  % incremnting iterative values by 1 to display the number of iteratations taken for the solution
  jacobi_iterative=jacobi_iterative+1;
  end
end

% gauss sidal iterative method:

if iterative_solvers==2
  gs_iterative=1; % defining iterative variable to store and display the number of iteratations computed for the values of convergence
  
  % convergence loop
  while(error>tol)
       % Iteratation loop.
  for i=2:nx-1;
  for j=2:ny-1
    % calculating tempat each grid points
    term1=(T(i-1,j)+Told(i+1,j))/(k*(dx^2));
    term2=(T(i,j-1)+Told(i,j+1))/(k*(dy^2));
    T(i,j)=term1+term2;
    end
    end
  error=max(max(abs(Told-T))); % 2D problem so two max
  
  % updating the temperature value for the next iteration
  Told=T;
  
  % incremnting iterative values by 1 to display the number of iteratations taken for the solution
  gs_iterative=gs_iterative+1;
  end
end


% SOR iterative method:

if iterative_solvers==3
  sor_iterative=1; % defining iterative variable to store and display the number of iteratations computed for the values of convergence
  
  % convergence loop
  while(error>tol)
       % Iteratation loop.
  for i=2:nx-1;
  for j=2:ny-1;
    % calculating tempat each grid points
    term1=(T(i-1,j)+Told(i+1,j))/(k*(dx^2));
    term2=(T(i,j-1)+Told(i,j+1))/(k*(dy^2));
    T(i,j)=(1-omega)*(Told(i,j))+omega*(term1+term2);
    end
    end
  error=max(max(abs(Told-T))); % 2D problem so two max
  
  % updating the temperature value for the next iteration
  Told=T;
  
  % incremnting iterative values by 1 to display the number of iteratations taken for the solution
  sor_iterative=sor_iterative+1;
  end
end
time_required=toc;
% plot representing the converged results
figure(1)
contourf(x,y,T);
colorbar;
colormap(jet);
%clabel(c,h);
pause(0.03) % labelling the contour plot,here c represents the temperature value and h rotates the temperature value as it moves

if iterative_solvers==1
  title_text=sprintf('steady state heat conduction using jacobi method  n  NO of iteratations: %d n computation time: %0.3f seconds',jacobi_iterative,time_required);
  title(title_text);
  saveas(figure(1),'steady_explicit_jacobi.jpeg');
  
 else if iterative_solvers==2
  title_text=sprintf('steady state heat conduction using gauss siedel method  n  NO of iteratations: %d n computation time: %0.3f seconds',gs_iterative,time_required);
  title(title_text);
  saveas(figure(1),'steady_explicit_gauss_siedel.jpeg');
  
  else if iterative_solvers==3
  title_text=sprintf('steady state heat conduction using SOR method  n  NO of iteratations: %d n computation time: %0.3f seconds',sor_iterative,time_required);
  title(title_text);
  saveas(figure(1),'steady_explicit_SOR.jpeg');
  end
  end
end
xlabel('X Length Domain');
ylabel('Y Length Domain');

 

When the equation was solved using the gauss seidal formula .It took 1725 iterations to get converged results. This is much lower compared to number of iterations taken to solve the same equation using the Jacobi method. (When the equation was solved using the Jacobi formula. It took 3227 iterations to get the converged results). Due to this the time required for computation reduces significantly compared to the Jacobi method. This is because the terms substituted in the discretized equation are the latest terms computed, rather substituting the terms obtained in the previous iteration.

 

From the above result, we can see that SOR has the fastest convergence rate, followed by Gauss seidal and then Jacobi. The SOR method gave converged solutions in 1436 iterations when compared to Gauss seidal with 1725  iterations and Jacobi with  3227 iterations. Consequently, the simulation time of the SOR method was times faster than the Gauss seidal and times faster than the Jacobi.

This is because the SOR method is an optimized from of the Gauss seidal which in turn, is an advanced form of Jacobi, It scales –up the convergence rate of the Gauss seidal method by a factor called the relaxation factor (omega), The relaxation factor can generally take values between 0and 2.But for the method Omega=1, the SOR method is the same as the Gauss seidal method. At omega-0, the solution does not update at all and at omega=2, the solution blows up.

The Convergence rate of SOR also depends on the choice of the relaxation factor, so, we must find an optimum value for the relaxation factor that gives the fastest convergence rate.

 

 Program for Transient state heat conduction(unsteady state Explicit method):

 

 

clear all
close all
clc


% initializing variables
% specfying length of the domain, L=1
% defining space domain
nx=40; % no of grid points in x axis
ny=40; % no of grid points in y axis
x=linspace(0,1,nx);
y=linspace(0,1,ny);
% grid size
dx=1/(nx-1);
dy=1/(ny-1);
alpha=1.4; % diffusivity coefficient in the 2D heat conduction equation
dt=0.0001; % time step
T=ones(nx,ny);
T(1,:)=600; % bottom
T(end,:)=900; % top
T(:,1)=400; % left
T(:,end)=800; % right

% copy T
Told=T;
% errors
error_tol=1e-4;
start_error=9e9;

%% CFL number
CFL_num=alpha*(dt/dx^2)+alpha*(dt/dy^2);

% solving 2D  unsteady state system
[xx,yy]=meshgrid(x,y);
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
iterative_solvers=1 % cannot  be kept inside the time loop as for every new time loopthe counter will start from1
tic;

% computing:
for nt=1:1600;
 %while(start_error>error_tol);
    for i=2:nx-1;
    for j=2:ny-1;
    term1=(1-2*k1-2*k2);
    term2=k1;
    term3=k2;
    h=Told(i-1,j)+Told(i+1,j);
    v=Told(i,j-1)+Told(i,j+1);;
    T(i,j)=Told(i,j)*term1+term2*h+term3*v;
    end
    end
    Told=T;
    %iterative_solvers=iterative_solvers+1;
  %end
  end
  time_required=toc;
  figure(1)
  contourf(xx,yy,T);
  colorbar;
  colormap(jet);
  title_text=sprintf('unsteady atate heat conduction using Explicit method n No.of Iterattations: %d ncomputation time : %0.3f seconds n CFL Number=%0.2f n time step dt=%0.4f',iterative_solvers,time_required,CFL_num,dt);
  title(title_text);
  xlabel('X Length Domain');
  ylabel('Y Length Domain');
  saveas(figure(1),'unsteady_explicit_dt=0.001_cf1=0.25.jpeg');
  pause(0.3)

 

 

  Program for Transient state heat conduction(unsteady state implicit method):

% program for 2D transient state heat conduction equation -unsteady state Implicity
clear all;
close all;
clc;

% initialize variables
nx=11; % no of grid points in x -axis
ny=11; % no of grid points in y-axis
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=1/(nx-1);
dy=1/(ny-1);
alpha=1.4; % diffusivity coefficient in the 2D heat conduction equation
% clacualting the time step
% dt=(0.5)/((alpha/dx^2)+(alpha/dy^2));
dt=0.001;
% tolerance
tol=1e-4;

% % temperature 
T=ones(nx,ny);
T(1,:)=600; % bottom
T(end,:)=900; % top
T(:,1)=400; % left
T(:,end)=800; % right

% copy T
Told=T;
T_prev_dt=Told;

% CFL number and omega
CFL_num=alpha*dt/(dx^2)+alpha*dt/(dy^2);
omega=1.1; % relaxation factor for successive over relaxation (SOR) method for transient heat conduction

% solving the 2D unsteady state system Implicity using iterative_solvers
[xx,yy]=meshgrid(x,y);
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
disp('For solver type: n=1 for jacobi; n =2 for Gauss seidal; n=3 for SOR n')
 iterative_solvers = input('Enter the solver number=')
tic;
% input 
%iterative_solvers=input('For solver type :n 1 for jacobi n 2 for gauss siedel n 3for SOR n')
jacobi_iterative=1; % cannot be kept inside the time loop as for every new time lop the counter will start from 1
gs_iterative=2;
sor_iterative=3;


% computing the time loop
for nt=1:1600;
  % jacobi iterative solver
  if iterative_solvers==1;
    term1=1/(1+2*k1+2*k2);
    term2=term1*k1;
    term3=term1*k2;
    error=1e9; % error outside this loop then next time loop error will be less than tolerance and program wont proceed
    while(error>tol);
    for i=2:nx-1;
    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_prev_dt(i,j)+term2*H+term3*V;
  end
end
error=max(max(abs(Told-T))); %  2D problem  so two max
Told=T;
jacobi_iterative=jacobi_iterative+1;
end
end

% Gauss siedel iterative solver
  if iterative_solvers==2;
    term1=1/(1+2*k1+2*k2);
    term2=term1*k1;
    term3=term1*k2;
    error=1e9; % error outside this loop then next time loop error will be less than tolerance and program wont proceed
    while(error>tol);
    for i=2:nx-1;
    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_prev_dt(i,j)+term2*H+term3*V;
  end
end
error=max(max(abs(Told-T))); %  2D problem  so two max
Told=T;
gs_iterative=gs_iterative+1;
end
end

% SOR iterative solver
  if iterative_solvers==3;
    term1=1/(1+2*k1+2*k2);
    term2=term1*k1;
    term3=term1*k2;
    error=1e9; % error outside this loop then next time loop error will be less than tolerance and program wont proceed
    while(error>tol);
    for i=2:nx-1;
    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)=(1-omega)*Told(i,j)+omega*(term1*T_prev_dt(i,j)+term2*H+term3*V);
  end
end
error=max(max(abs(Told-T))); %  2D problem  so two max
Told=T;
sor_iterative=sor_iterative+1;
end
end
T_prev_dt=T;
end
time_required=toc;
figure(1);
[c,h]=contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(c,h);
if iterative_solvers==1
  title_text=sprintf('unsteady state heat conduction using jacobi method n No.of iteratation :%d n computation Time: %0.3f seconds n CFL number=%0.3f& dt=%0.4f',jacobi_iterative,time_required,CFL_num,dt);
  title(title_text)
  saveas(figure(1),'unsteady_Implicit_jacobi_dt=0.001_cf1=0.28.jpeg');
  
else if iterative_solvers==2
  title_text=sprintf('unsteady state heat conduction using gauss seidel method n No.of iteratation :%d n computation Time: %0.3f seconds n CFL number=%0.3f& dt=%0.4f',gs_iterative,time_required,CFL_num,dt);
  title(title_text)
  saveas(figure(1),'unsteady_Implicit_gauss_seidel_dt=0.001_cf1=0.28.jpeg');
  
else if iterative_solvers==3
  title_text=sprintf('unsteady state heat conduction using SOR method n No.of iteratation :%d n computation Time: %0.3f seconds n CFL number=%0.3f& dt=%0.4f',sor_iterative,time_required,CFL_num,dt);
  title(title_text)
  saveas(figure(1),'unsteady_Implicit_SOR_dt=0.001_cf1=0.28.jpeg');
  end
  end
end

xlabel('X Length in Domain');
ylabel('Y Length in domain');

 

 

From the above results, we can see that the explicit method gives the solution fastest. This is because it does not require guess values and iterations for convergence as the explicit solution is solely based on the previous state values which had already been completed.

The implicit method requires guess value and iterations for convergence as it is also based on the future state value, which are not known at the time of computation.it checks convergence for each time step. This cause a delay in simulation.

From the above points .It can be found that the explicit scheme is the fastest to solve the transient state 2D heat equation because we are not bothered about iterations and convergence, the scheme uses already defined previous state value to compute present state value. So there is no need for guess values and convergence criteria. However the explicitly scheme is strictly governed by the CFL stability condition.  For a scheme to be stable.

Where alpha is called the thermal diffusivity and it is taken as1.4 in this case for a faster solution.

This is a huge disadvantage for the explicit scheme because we must use extremely small values of dt as compared to dx for the solution to be stable, which means longer computational time .The above explicit solution as the CFL number is  0.28 which is well below 0.5 . And the solution value for dt=1e-3(CFL=0.3) as the courant’s number approaches the limit.

The simplicity scheme are slower and more complicated than the explicit scheme, since were using guess values and iterating for convergence in each time step, However , they are unconditionally stable for any choice of dt with respect to dx. They offer better stability and larger time step size.

Below table gives number of iterations required by iterative solvers to converge.

 

Steady state Heat Conduction Equation:

Iterative solver	No of iterations
Jacobi	3227
Gauss seidal	1725
SOR	1436

 Unsteady state Heat Conduction Equation:

 

Iterative solver	No of iterations
jacobi	4384
Gauss seidal	3650
SOR	3216

Conclusion:

Steady state conduction is in the form of conduction that happens when the temperature difference driving the conduction is constant with respect to time. On the other hand, in transient –state condition, the temperature at the boundaries change with time, and eventually reach a steady state if the temperature is not changed further.

If the error term obtained is below the tolerance value, the solution is said to be converged. In the above problem which is solved in steady and unsteady – state condition, we obtain an error term below the tolerance (1e-4) for each case hence the solution is converged

In the steady state analysis and transient state analysis using the implicit method, the solution is obtained using three methods namely Jacobi, gauss seidal, successive over relaxation method. Out of this three methods, the third method, SOR convergence to the solution with a low number of iterations compared to other two methods. Provided value omega is approximate> the Jacobi method in both cases converges to the solution with the highest number of iterations compared to the other two methods.

From this we can conclude that the SOR method is the most effective method, since it converges the solution for less number of iterations and hence less time is taken for computing to Jacobi and Gauss seidal method. This largely depends on omega value , if an approximate value chosen, then the solution converges faster, else the solution may take longer to converge

Implicit method can be used to problems when converged results are needed in less computation time as they are unconditionally stable. With increase in time step the convergence rate is increased which leads to steady state with less no of iterations

Explicit method results time accurate solutions as it is governed by stability criteria, when solution not satisfied it leads to unstable solution and requires more iterations to attend steady state.

But in the code provided above time step size remains constant which is calculated with respect to time stability criteria. Hence for the time step size the explicit scheme leads to faster convergent rate and as the time step small the number of iterations increases for implicit scheme which results slower convergence rate.