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

Aim:

The aim of the code is to display the temperature distribution acrros a 2D plate which has been subjected to constant temperatures across the boundary. This would show the conduction of heat across the plate

Design of code:

• The temperature distribution for considering steady as well as transient state(impicit & explicit) has been displayed
• A contour depicting the temperature in all three cases has also been provided.
• Since steady state as well as implicit transient state uses iterative methods, a switch case statement has been provided in order to chose between SOR, Jacobi and Gauss Seidel techniques

Formulae

For steady state

`(d^2T)/dx^2 + (d^2T)/dy^2 = 0`

`(Tr - 2Tp + Tl )/h^2 + (Tt - 2Tp + Tb )/k^2 = 0`

`(Tr - 2Tp + Tl )*k^2/h^2 + (Tt - 2Tp + Tb ) = 0``c^2 * (Tr + Tl ) + (Tt + Tb ) -2* Tp*(1 + c^2) = 0`

`(c^2 * (Tr + Tl ) + (Tt + Tb ) )/(2*(1 + c^2)) = Tp`

For transient state its the normal formula as discussed in the lecture

Variables:

Most of the variables are self explanatory or have already been mentioned aongside

• err = maximum row or column error in temperature matrix between 2 iteration
• iter = iteration
• a = themal conductivity of metal
• Told = Temperature matrix for previous iteration for a given time step
• Tp = Temperature matrix for the previous time step
• nt = number of time steps

  close all
  clear
  clc
  % equation
  % d^2T/dx^2 + d^T/dy^2 = 0; steady state
  % d^2T/dx^2 + d^T/dy^2 = (dT/dt)/a; transient state
  
  % initial condition
  nx = 10; ny = 10;
  h = 1/nx; k = 1/ny ;% values of dx and dy respectively
  T = 400 * ones(nx,ny);
  T(1,1:end) = 600; T(1:end,1) = 400;
  T(end,1:end) = 900; T(1:end,end) = 800;
  T(1,1) = avg(400,600); T(1,end) = avg(600,800); % avg is a function to determine average for the corner points
  T(end,1) = avg(400,900); T(end,end) = avg(800,900);
  Told = T; 
  c = h/k;
  err = 99;
  choice = "SOR";
  w = 1.5; %optimal over relaxation parameter
  iter = 0;
  %steady state
  
  while err > 1e-4
      for i = 2:nx-1
          for j = 2:ny-1
              switch choice
                  case{"J"}
                      T(i,j) = (c^2*(Told(i-1,j) + Told(i+1,j)) + Told(i,j-1) + Told(i,j+1))/(2*(1+c^2));
                  case{"GS"}
                      T(i,j) = (c^2*(T(i-1,j) + Told(i+1,j)) + T(i,j-1) + Told(i,j+1))/(2*(1+c^2));
                  case{"SOR"}
                      T(i,j) = (1-w)*Told(i,j) + w*(c^2*(T(i-1,j) + Told(i+1,j)) + T(i,j-1) + Told(i,j+1))/(2*(1+c^2));
              end 
          end
      end
       err = max(max(abs(T-Told)));
       Told = T;
       iter = iter + 1;
  end
  figure(1)
  contourf(T)
  
        
  
  %transient state
  T = ones(nx,ny);%considering the plat to have a temperature of 1K before boundary conditions are applied 
  T(1,1:end) = 600; T(1:end,1) = 400;
  T(end,1:end) = 900; T(1:end,end) = 800;
  T(1,1) = avg(400,600); T(1,end) = avg(600,800); % avg is a function to determine average for the corner points
  T(end,1) = avg(400,900); T(end,end) = avg(800,900);
  Tp = T; a = 1.4; dt = 0.002; time = 5;
  nt = time/dt;
  k1 = a*dt/h^2; k2 = a*dt/k^2;
  choice  = "J"; w = 1.09; iter1 = [];%iter1 is the matrix containing no. of iterations required per time step 
  Te = T;% for explicit scheme later on
  
  %implicit form
  Told = Tp;
  
  for t = 1:nt
      iter = 0;
      err = 99;
      while err > 1e-4
          for i = 2:nx-1
              for j = 2:ny-1
                  switch choice
                      case{"J"}
                          T(i,j) = (Tp(i,j) + k1*( Told(i,j-1) + Told(i,j+1) ) + k2*(Told(i+1,j) + Told(i-1,j) )) /(1+ 2*(k1+k2) );
                      case{"GS"}
                          T(i,j) = (Tp(i,j) + k1*( T(i,j-1) + Told(i,j+1) ) + k2*(Told(i+1,j) + T(i-1,j) )) /(1+ 2*(k1+k2) );
                      case{"SOR"}
                          T(i,j) = (1-w)*Told(i,j) + w*((Tp(i,j) + k1*( T(i,j-1) + Told(i,j+1) ) + k2*(Told(i+1,j) + T(i-1,j) )) /(1+ 2*(k1+k2) ));
                  end
              end
          end
          err = max(max(abs(T - Told))); 
          Told = T;
          iter = iter + 1;
      end
      Tp = T;
      iter1(t) = iter;
  end
  iter = max(iter1);
  figure(2)
  contourf(T)
  
  % explicit form
  Tp = Te;
  for t = 1:35 % no. of iterations uptill which value is converging properly
      for i = 2:nx-1
          for j = 2:ny-1
              T(i,j) = Tp(i,j) + k1*(Tp(i,j-1) + Tp(i,j+1) ) + k2*(Tp(i+1,j) + Tp(i-1,j) ) 
  - 2*Tp(i,j)*(k1 + k2);
          end
      end
      Tp = T;
      figure(3)
      contourf(T)
      pause(0.1) %in order to see the value converge
  end
  
  

Observations:

From the temperature distribution in general

• The temperature along the bottom boundary having Temp = 900K is higher than the rest of the plate
• The temperature along the top boundary having Temp = 400K is the lowest
• Initially in the transient state, the lowest temmperature is initially at the centre, and the overall temperture is also low
• As time passes, lower temperature shifts towards the boundary and the boundary temperature zone starts becoming broader

Regarding the techniques used

• SOR method when used with the right value of w, has least number of iterations( 28 in steady state, 12 in transient)
• w having value higher or lower than 1.5 yields greater number of iterations in steady state
• w having value higher or lower than 1.1 yields greater number of iterations in transient state
• In the explicit transient scheme, higher number of time steps yield non convergent answer
• Around 40 time steps provides a convergent response, anything beyond that makes the results fluctuates
• The value of dt(which affects the CFL number k1,k2) drastically affects the end value and number of iterations
• As dt increases, the number of iterations increase in implicit technique, but still yields a convergent solution
• However in the explicit scheme, a non convergent solution is obtained if dt is high

Inference:

• SOR method has least iterations
• Implicit method is inherently stable because value is converged in each time step
• Explicit scheme is highly volatile