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

AIM :

To write the code which solves the steady and unsteady 2D heat conduction equation by using the point iterative techniques by implementing the following methods.

• Jacobi
• Gauss-Seidel
• Successive over-relaxation

THEORY :

Considering general 3D heat condution equation which has both time derivative and the spatial derivative terms.

`(delT)/(delt)+alpha*(((del^2T)/(delx^2))+((del^2T)/(dely^2))+((del^2T)/(delz^2)))=0`

In the above equation 1st term represents time derivative and 2nd term is spatial derivative, where alpha is thermal diffusivity.

Now by considering the assumption that the above equation is solved for a thin plate, dimensions in the Z direction are neglected, thereby the above equation is converted into the following form which is as follows:

`(delT)/(delt)+alpha*(((del^2T)/(delx^2))+((del^2T)/(dely^2)))=0`

The above equation can be calculated by considering two approaches:

1.Steady State equation :

Under these conditions, the rate of change of temperature remains 0. Therefore in this condition time, derivative term is equated to zero.

Steady-state 2D heat conduction equation:

`((del^2T)/(delx^2))+((del^2T)/(dely^2))=0`

2.Transient State equation :

Under these conditions, both terms are considered therefore the equation is as follows:

`(delT)/(delt)+alpha*(((del^2T)/(delx^2))+((del^2T)/(dely^2)))=0`

Both the above 2 equations can be solved in 2 ways by 3 iterative methods like Jacobi,Gauss Seidel and SOR.

1.Explicit method

2.Implicit method

We have used transient state equation using both ways.

Assumptions taken for solving the problem:

1.Assume that the domain is a unit square.

2.Assume nx = ny=10(for both steady and unsteady) [Number of points along the x-direction is equal to the number of points along the y direction]

3.Boundary conditions for steady and transient case

Left:400K , Top:600K , Right:800K , Bottom:900K , Tolerance=1e-4 , Starting error=9e9 , Thermal diffusivity(Alpha)=1.4m^2/sec , Domain length=1m

CFL(for stability)= `CFL_x+CFL_y`, =`(alpha*(dt/dx^2)+alpha*(dt/dy^2))`

Solving Steady state equation

Considering Two dimensional Steady State Equation

`((del^2T)/(delx^2))+((del^2T)/(dely^2))=0`

Discretizing the equation by applying central differencing scheme ,we get

`T_(i,j)=1/k*{((T_(old)(i,j))+(T_(old)(i+1,j)))/dx^2+((T_(old)(i,j-1))+(T_(old)(i,j+1)))/dy^2}`

`k=2*((dx^2+dy^2)/(dx^2*dy^2))`

Code for solving the above equation using 3 iterative methods :

main program

%Solving the Steady state equation  by employing
%different solvers like Jacobi,Gauss Seidel,Sucessive Over Relaxation
clear all
close all
clc
%defining the inputs
nx=10;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,nx);
dx=x(2)-x(1);
dy=y(2)-y(1);
[xx,yy] = meshgrid(x,y);

%defining the boundary conditions
T=ones(nx,ny);
T(1,:)=900;
T(end,:)=600;
T(:,1)=400;
T(:,end)=800;
Told=T;
%error_tol,start error
error_tol=1e-4;
start_error=9e9;
k=2*((dx^2+dy^2)/(dx^2*dy^2));
%defining the solvers
iterative_solver = input('Enter the desired number w.r.t solver:');
if iterative_solver==1
    tic;
    T_jacobi= jacobi_iter(start_error,error_tol,xx,yy,nx,ny,dx,dy,T,k,Told);
    time_jacobi_in_sec=toc
elseif iterative_solver==2
    tic;
    T_gauss= gauss_iter(start_error,error_tol,x,y,nx,ny,dx,dy,k,T,Told);
    time_gauss_in_sec=toc
elseif iterative_solver==3
    tic;
    T_sor= sor_iter(start_error,error_tol,x,y,nx,ny,dx,dy,k,T,Told);
    time_sor_in_sec=toc
end

Jacobi function

%function for Jacobi method to solve in Steady state equation
function T= jacobi_iter(start_error,error_tol,xx,yy,nx,ny,dx,dy,T,k,Told)
     
     jacobi_iter=1;
     while(start_error > error_tol)
     for i=2:nx-1
         for j=2:ny-1
          term1 = (Told(i-1,j)+Told(i+1,j))/(k*(dx^2));
          term2 = (Told(i,j-1)+Told(i,j+1))/(k*(dy^2));
          T(i,j) = term1+term2;
         end
     end
     start_error=max(max(abs(Told-T)));
     Told=T;
     jacobi_iter=jacobi_iter+1;
     figure(1)  
     [c,d]=contourf(xx,yy,T);
     colormap(jet)
     colorbar
     title_text=sprintf('Jacobi iteration number=%d at steady state',jacobi_iter);
     title(title_text)
     xlabel('X-axis')
     ylabel('Y-axis')
     clabel(c,d)
     pause(0.005)
     end
     fprintf('No. of Iteration in Jacobi=%d',jacobi_iter);
end

Gauss Seidel function

%function for Gauss_Seidel method to solve in Steady state method
function gauss_seidel= gauss_iter(start_error,error_tol,x,y,nx,ny,dx,dy,k,T,Told)
     GS_iter=1;
   
     while(start_error > error_tol)
     for i=2:nx-1
         for j=2:ny-1
         T(i,j)=(((T(i-1,j)+Told(i+1,j))/dx^2)+((T(i,j-1)+Told(i,j+1))/dy^2))/k;
         end
     end
     start_error=max(max(abs(Told-T)));
     Told=T;
     GS_iter=GS_iter+1;
     gauss_seidel=figure(2)   
     [c,d]=contourf(x,y,T);
     colormap(jet)
     colorbar
     title_text=sprintf('GS iteration number=%d at steady state',GS_iter);
     title(title_text)
     xlabel('X-axis')
     ylabel('Y-axis')
     clabel(c,d)
     pause(0.005)
     end
     fprintf('No. of Iteration in GS=%d',GS_iter);
end

SOR function

%function for Successive Over Relaxation method to solve in Steady state method
function T= sor_iter(start_error,error_tol,x,y,nx,ny,dx,dy,k,T,Told)
     sor_iter=1;
     omega=1.2; 
     while(start_error > error_tol)
     for i=2:nx-1
         for j=2:ny-1
         T_gs=(((T(i-1,j)+Told(i+1,j))/dx^2)+((T(i,j-1)+Told(i,j+1))/dy^2))/k;
         T(i,j)= omega*(T_gs)+(1-omega)*(Told(i,j));
         end
     end
     start_error=max(max(abs(Told-T)));
     Told=T;
     sor_iter=sor_iter+1;
     figure(3)   
     [c,d]=contourf(x,y,T);
     colormap(jet)
     colorbar
     title_text=sprintf('SOR iteration number=%d at steady state',sor_iter);
     title(title_text)
     xlabel('X-axis')
     ylabel('Y-axis')
     clabel(c,d)
     pause(0.005)
     end
     fprintf('No. of Iteration in SOR=%d',sor_iter);
end

OUTPUT :

Enter the desired number w.r.t solver:1
No. of Iteration in Jacobi=217
time_jacobi_in_sec =

   52.2490

No. of Iteration in GS=116
time_gauss_in_sec =

   30.6671

Enter the desired number w.r.t solver:3
No. of Iteration in SOR=78
time_sor_in_sec =

   18.7859

From the above plots we can infer that SOR takes both less time and less number of iterations than the other 2 methods. So we can conclude that SOR is the best among the three.

Solving Transient state equation

Transient state 2D Steady state heat conduction equation is solved in both ways i.e explicitly and implicitly. Following is the 2D heat conduction equation

`(delT)/(delt)+alpha*(((del^2T)/(delx^2))+((del^2T)/(dely^2)))=0`

By explicit method

Explicit Approach is used in numerical analysis for obtaining numerical approximations to the solution of time-dependent ordinary or partial differential equations as in required in computer simulations of a physical process. Explicit methods calculate the state of a system at a later time from the state of the system at the current time.

In this method, each differential equation approach contains only one unknown and can be solved explicitly for this one unknown using this method.

Here in this equation time derivative terms are calculated using the first order differential equations and whereas the term 2 has second-order derivative terms and elliptical in nature so information needs to be taken from all the points equally, So, therefore, a central differencing scheme is applied for the discretising the term2.

`(delT)/(delt)approx(T_(i,j)^(n+1)-T_(i,j)^(n))/(Deltat)`

`(del^2T)/(delx^2)approx(T_(i+1,j)^(n)-2T_(i,j)^(n)+T_(i-1,j)^(n))/(Deltax^2)`

`(del^2T)/(dely^2)approx(T_(i,j+1)^(n)-2T_(i,j)^(n)+T_(i,j-1)^(n))/(Deltay^2)`

The discretized equation for 2D heat conduction equation for the explicit form is as follows:

`T^(n+1)(i,j)=T^n+k1*(T_(old)^n(i+1,j))-2*(T_(old)^n(i,j))+(T_(old)^n(i-1,j))+k2*{T_(old)^n(i,j-1))-2*(T_(old)^n(i,j))+(T_(old)^n(i,j+1)}`

`k1=2*((alpha*dt)/(dx^2)) ,k2=2*((alpha*dt)/(dy^2))`

code for 2D Transient Heat conduction equation solving explicitly:

main program

%solving 2D Transient heat conduction equation explicitly
clear all
close all
clc
%defining the inputs
nx=10;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,nx);
dx=x(2)-x(1);
dy=y(2)-y(1);
dt=1e-3; 
alpha=1.4; %thermal diffusivity
[xx,yy] = meshgrid(x,y);

%defining the boundary conditions
T=ones(nx,ny);
T(1,:)=900;
T(end,:)=600;
T(:,1)=400;
T(:,end)=800;
Told=T;
%rounding off the corner value by taking the average value
T(1,1) = (400+900)/2;
T(1,ny) = (900+800)/2;
T(nx,1) = (400+600)/2;
T(nx,ny) = (600+800)/2;
%error_tol,start error
error_tol=1e-4;
start_error=9e9;
k1= alpha*(dt/dx^2);
k2= alpha*(dt/dy^2);
CFL_1= alpha*(dt/(dx^2));
CFL_2= alpha*(dt/(dy^2));
CFL=CFL_1+CFL_2;
%defining the solvers
iterative_solver = 1;
if iterative_solver==1
    tic;
    T_jacobi= jacobi_iter_transient_expli(xx,yy,nx,ny,T,k1,k2,Told);
end

Jacobi function :

%function for Jacobi method to solve in transient state equation(explicit)
function T= jacobi_iter_transient_expli(xx,yy,nx,ny,T,k1,k2,Told)
     jacobi_iter=1;
    for k=1:1500
     for i=2:nx-1
         for j=2:ny-1
          term1 = Told(i,j);
          term2 = k1*(Told(i+1,j)-2*Told(i,j)+Told(i-1,j));
          term3= k2*(Told(i,j+1)-2*Told(i,j)+Told(i,j-1));
          T(i,j) = term1+term2+term3 ;
         end
     end
     start_error=max(max(abs(Told-T)));
     Told=T;
     jacobi_iter=jacobi_iter+1;
    end
     time_jacobi_in_sec=toc
     figure(1)  
     [c,d]=contourf(xx,yy,T);
     colormap(jet)
     colorbar
     title_text=sprintf('Jacobi iteration number=%d at transient state explicitly and time taken in sec is=%f',jacobi_iter,time_jacobi_in_sec);
     title(title_text)
     xlabel('X-axis')
     ylabel('y-axis')
     clabel(c,d)
     pause(0.005)
     fprintf('No. of Iteration in Jacobi=%d ',jacobi_iter);
    
end

OUTPUT :

time_jacobi_in_sec =

    0.0171

No. of Iteration in Jacobi=1501

By Implicit method:

Implicit Approach is used in numerical analysis for obtaining numerical approximations to the solution of time-dependent ordinary or partial differential equations as in required in computer simulations of a physical process.

Implicit methods find the solution by solving an equation involving both the current state of the system and the later one. Mathematically if Y(t) is the current system state and Y(T+dt) is the state at the later time dt, then, for an explicit method Y(T+dT)=F(Y(t))

Here in this equation time derivative terms are calculated using the first order differential equations and whereas the term 2 has second-order derivative terms and elliptical in nature so information needs to be taken from all the points equally, So, therefore, a central differencing scheme is applied for the discretising the term2.

`(delT)/(delt)approx(T_(i,j)^(n+1)-T_(i,j)^(n))/(Deltat)`

`(del^2T)/(delx^2)approx(T_(i+1,j)^(n+1)-2T_(i,j)^(n+1)+T_(i-1,j)^(n+1))/(Deltax^2)`

`(del^2T)/(dely^2)approx(T_(i,j+1)^(n+1)-2T_(i,j)^(n+1)+T_(i,j-1)^(n+1))/(Deltay^2)`

`k1=2*((alpha*dt)/(dx^2)) ,k2=2*((alpha*dt)/(dy^2))`

The above equation is solved iteratively using different iterative methods such as Jacobi, Gauss-Seidel and SOR method.

In Jacobi, the unknowns on the right-hand side are assumed to have an initial value, i.e. T at 'n+1' at all nodes is assumed. Using the equation, then the LHS is calculated. On the next value of i,j similar method is applied to calculate LHS. Thus in the first iteration, value of T for time 'n+1' is calculated at all nodes. The difference between T at 'n+1' which is obtained and T at 'n+1' which was assumed previously is taken. If it is greater than a specified tolerance value, iterations are carried out to calculate T at 'n+1'. But in the second iteration, T at 'n+1' is assumed as the value of T which was obtained in the previous iteration. Once the value of T is found out iteratively whose error is less than the tolerance value, the value of T is calculated at the next time step 'n+2' using the same iterative method.

The Gauss Seidel method is similar to the Jacobi method, but the only difference is- For the first iteration when the value of T at a particular node is calculated, this updated value at that node is used to calculate the value at the next node. But, the error term is calculated similar to Jacobi method, i.e. T at 'n+1' calculated minus the T at 'n+1' previously assumed. Once the value of T is found out iteratively whose error is less than the tolerance value, the value of T is calculated at the next time step 'n+2' using the same iterative method.

SOR method can be applied either by using Jacobi or Gauss seidal but the only difference is consideration of factor called omega(scaling factor) which in turn makes the process faster by which the number of iterations required to solve will be decreased

`omega<1 ,` under-relaxation

`omega>1 ,` over-relaxation

code for 2D Transient Heat conduction equation solving implicitly:

main program

%solving 2D Transient heat conduction equation implicitly
clear all
close all
clc
%defining the inputs
nx=10;
ny=nx;
L=1;
x=linspace(0,1,nx);
y=linspace(0,1,nx);
dx=x(2)-x(1);
dy=y(2)-y(1);
dt=1e-3;   
alpha=1.4;%thermal diffusivity
[xx,yy] = meshgrid(x,y);

%defining the boundary conditions
T=ones(nx,ny);
T(1,:)=900;
T(end,:)=600;
T(:,1)=400;
T(:,end)=800;
%rounding off the corner value by taking the average value
T(1,1) = (400+900)/2;
T(1,ny) = (900+800)/2;
T(nx,1) = (400+600)/2;
T(nx,ny) = (600+800)/2;
Told=T;
Tprev=Told;
%error_tol,start error
error_tol=1e-4;
start_error=9e9;
k1= alpha*(dt/(dx^2));
k2= alpha*(dt/(dy^2));
omega=1.2;
CFL_1= alpha*(dt/(dx^2));
CFL_2= alpha*(dt/(dy^2));
CFL=CFL_1+CFL_2;
%defining the solvers
iterative_solver = input('input the appropriate number to select required the method:');
if iterative_solver==1
    tic;
    T=jacobi_iter_trans_impli(error_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev);
elseif iterative_solver==2
    tic;
    T_gs= gs_iter_trans_impli(error_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev);
elseif iterative_solver==3
    tic;
    T_sor= sor_iter_trans_impli(error_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev,omega);
end
 

Jacobi function :

%function for Jacobi method to solve in Transient state equation(implicit)     
function   T=jacobi_iter_trans_impli(error_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev)   
    jacobi_iter=1; 
for k= 1:1500 %time loop
    start_error=9e9;
     while(start_error > error_tol)%convergence loop
     for i=2:nx-1  %nodal loop
         for j=2:ny-1
          term1 = (1+(2*k1)+(2*k2))^(-1);
          term2= k1*term1;
          term3= k2*term1; 
          H=(Told(i-1,j)+Told(i+1,j));
          v=(Told(i,j-1)+Told(i,j+1));
          T(i,j)= (Tprev(i,j)*term1)+term2*H+term3*v;
         end
     end
     start_error=max(max(abs(Told-T)));
     Told=T;
     jacobi_iter=jacobi_iter+1; 
     end
    Tprev=T;
end
     time_jacobi=toc;
     contourf(xx,yy,T,'ShowText','on')
     colormap(jet)
     colorbar
     xlabel('X-axis')
     ylabel('Y-axis')
     pause(0.005)
     colorbar
     title_text=sprintf('Solution for 2D dimensional heat conduction equation solved implicitly using jacobi method with no of iteration=%d and time to solve=%f sec',jacobi_iter,time_jacobi);
     title(title_text)
     fprintf('No of Iteration in Jacobi =%d and time in sec=%f\n',jacobi_iter,time_jacobi);
end

GS function :

%function for GS method to solve in Transient state equation(implicit)
function T= gs_iter_trans_impli(error_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev)
 gs_iter=1;
for k=1:1500
         start_error=9e9;
    while(start_error > error_tol)  
     for i=2:nx-1
         for j=2:ny-1
          term1 = (1+(2*k1)+(2*k2))^-1;
          term2= k1*term1;
          term3= k2*term1;
          H=(T(i-1,j)+Told(i+1,j));
          v=(T(i,j-1)+Told(i,j+1));
          T(i,j)= Tprev(i,j)*term1+term2*H+term3*v;
         end
     end
     start_error=max(max(abs(Told-T)));
     Told=T;
     gs_iter=gs_iter+1; 
    end
    Tprev=T;
end
    time_gs=toc;
    contourf(xx,yy,T,'ShowText','on')
    colormap(jet)
    colorbar
    xlabel('X-axis')
    ylabel('Y-axis')
    pause(0.005)
    colorbar
    title_text=(sprintf('Solution for 2D dimensional heat conduction equation solved implicitly using GS method with no of iteration=%d and time to solve=%f sec',gs_iter,time_gs));
    title(title_text)
    fprintf('No of Iteration in GS=%d and time in sec=%f ',gs_iter,time_gs);
end

SOR function :

%function for SOR method to solve in Transient state equation(implicit)
function T=sor_iter_trans_impli(error_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev,omega)
     sor_iter=1;
 for k=1:1500
     start_error=9e9;
    while(start_error > error_tol)
     for i=2:nx-1
         for j=2:ny-1
          term1 = (1+(2*k1)+(2*k2))^-1;
          term2= k1*term1;
          term3= k2*term1;
          H=(T(i-1,j)+Told(i+1,j));
          v=(T(i,j-1)+Told(i,j+1));
          T(i,j)= omega*(Tprev(i,j)*term1+term2*H+term3*v)+(1-omega)*Tprev(i,j);
         end
     end
     start_error=max(max(abs(Told-T)));
     Told=T;
     sor_iter=sor_iter+1; 
    end
     Tprev=T;
 end
     time_sor=toc;
     contourf(xx,yy,T,'ShowText','on')
     colormap(jet)
     colorbar
     xlabel('X-axis')
     ylabel('Y-axis')
     pause(0.005)
     colorbar
     title_text=(sprintf('Solution for 2D dimensional heat conduction equation solved implicitly using SOR method with no of iteration=%d and time to solve=%f sec',sor_iter,time_sor));
     title(title_text)
     fprintf('No of Iteration in SOR=%d and time in sec=%f',sor_iter,time_sor);
    
end

OUTPUT :

input the appropriate number to select required the method:1
No of Iteration in Jacobi =3983 and time in sec=0.047428

input the appropriate number to select required the method:2
No of Iteration in GS=3371 and time in sec=0.037069 

input the appropriate number to select required the method:3
No of Iteration in SOR=3157 and time in sec=0.038242

We can see from the plots of jacobi method in both explicit and implicit that ,explicit method has lesser number of iterations(1501<3983), also evident from time taken in sec that explicit takes lesser time due to lesser iterations(0.047>0.017 seconds) present.This extra time is due to the fact that we are checking for an additional loop i.e (starting error > error tolerance), also 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 computed.

Also we can infer that in general SOR method > GS method > Jacobi method in implicit form as one can see from their no. of iterations, which is least in SOR-hence faster convergence rate .

The plots shown above are the transition of explicit plot to completely blown out state (due to instability) which can be inferred from CFL(Courant–Friedrichs–Lewy or Courant number) values , it has to be less than 0.5 for stability criteria to be satisfied(ie `CFL_x+CFL_y`) (we noticed it from 0.5216).

For the implicit solution, it was noted that the solution remained stable for all values of dt, irrespective of node step size. This reassures that the implicit solution is unconditionally stable. (we even checked for CFL value of 200)

Even though an implicit solution is unconditionally stable, it has several disadvantages like more complexity, larger error and more computational time.