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

1. AIM

To write a program to solve the 2D heat conduction equation by both implicit and explicit schemes

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation

1. METHODOLOGY
• INTRODUCTION

Assumption

• Domain is a unit square and number of grid points along x and y are equal
• Boundary Conditions: (Drichlet type)

1. Left wall = 400K
2. Top wall = 600K
3. Right wall = 800K
4. Bottom wall = 900K

STEADY STATE

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

UNSTEADY STATE

2 Dimensional Heat Conduction Equations Unsteady state:

`(∂^2T)/(∂x^2)+(∂^2T)/(∂y^2) =α (∂T)/(∂t)`

Discretisation of PDE spatial terms using Central Difference and time terms using forward difference then the final equation is :

`( T ^(n + 1) )_ (i , j )=( ( T^ n )_( i , j) + ( K 1 ( T_ ((i − 1 ), j) + T_ (i + 1 , j ) )+ K 2 ( T_( i , j − 1) + T_( i , j + 1) ) )^( n + 1)) /(1 + 2 K 1 + 2 K 2)`

where

`K 1 = α Δ t/ ( Δ x ) ^2`

`K 2 = α Δ t/ ( Δ y ) ^2`

n is the time index

(i,j) is the space index

Given:

omega = 1.2;

alpha=1.4;

dx=dy;

tolerance value = `10^-4`;

Inspite of the fact that conduction is a time dependent phenomena, it finally culminates to a steady state phenomena with time.

The condition for convergence of the 2D coduction equation is given by Courant-Lewy-Friedrichs Condition. Since we have distcretized the PDE there will be error terms which progrates with the solution as it moves form one time step to another in case of transient solver. Thus it is mandatory to minimze this effect or oscillate with damped amplitude so that it gradually alleviates as it moves forward with respect to time. CFL is an important criteria which when defined sets a condition that if the time time step is above a limit the convergence is not possible by the solver.

`CFL = α * ( 1/ ( Δ x ) ^2 + 1 /( Δ y )^ 2 ) ⋅ Δ t<= 0.5`

where:
α = Thermal Diffusivity
Δx=Δy = No. of grid points in the x or y direction
Δt = Time step

1. Unsteady state explicit Iteration Scheme

MATLAB CODE

% MATLAB Program  2D Transient State Heat Conduction Equation-UnSteady State Explicitly
clear all;
close all;
clc;
%Initialize variables
Nx =11;
Ny = 11;
x=linspace(0,1,Nx);
y=linspace(0,1,Ny);
dx=1/(Nx-1);
dy=1/(Ny-1);
alpha=1.4;% faster convergence as suggested in help video
dt=0.001;
%%Temperature field
T=298*ones(Nx,Ny);
T(1,:) = 600; % bottom
T(end,:)=900; % top 
T(:,1) = 400; % left 
T(:,end)=800; % right 
% Corner refining
T(1,1)=(600+400)/2;
T(1,end) = (600+800)/2;
T(end,1) = (400+900)/2;
T(end,end) = (800+900)/2;
%copy T
Told=T;
%%CFL Number
CFL_num=alpha*dt/(dx^2)+ alpha*dt/(dy^2);
%% Solving the 2D unsteady state system 
[xx,yy] = meshgrid(x,y);
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
iter_solver=1; % cannot be kept inside the time loop as for every new time loop the counter will start from 1
tic;
%computing:
for nt=1:1600
 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
 Told = T;
 iter_solver = iter_solver+1;
end
time_required=toc;
figure(1)
[C,h]=contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
set(gca,'YDIR','reverse');
title_text= sprintf('Unsteady state heat conduction using Explicit method \n No. of Iterations: %d \n Computation Time: %0.3f seconds (Time %0.2f s)\n CFL Number=%0.2f & Time step dt = %0.4f',iter_solver,time_required, nt*dt,CFL_num,dt);
title(title_text);
xlabel('xx');
ylabel('yy');
saveas(figure(1),'unsteady_explicit_dt=0.001_cfl=0.25_','jpeg');

Figure: Unsteady State - Explicit Scheme

Inference: It is clear from the results that the solution has good convergence and statble since CFL value is less than 0.5.

2. Unsteady state Implicit Iteration Scheme

MATLAB CODE

% MATLAB Program  2D Transient State Heat Conduction Equation-UnSteady state Implicitly
clear all;
close all;
clc;
%Initialize variables
Nx = 11;
Ny = 11;
x=linspace(0,1,Nx);
y=linspace(0,1,Ny);
dx=1/(Nx-1);
dy=1/(Ny-1);
alpha=1.4;% faster convergence as suggested in help video
dt=0.001;
tol=1e-4;
%%Temperature field
T=298*ones(Nx,Ny);
T(1,:) = 600; % bottom
T(end,:)=900; % top 
T(:,1) = 400; % left 
T(:,end)=800; % right 
% Corner refining
T(1,1)=(600+400)/2;
T(1,end) = (600+800)/2;
T(end,1) = (400+900)/2;
T(end,end) = (800+900)/2;
%copy T
Told=T;
T_prev_dt = Told;
%%CFL Number and omega
CFL_num=alpha*dt/(dx^2)+ alpha*dt/(dy^2);
omega=1.2;
%% 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);
%% Input
iter_solver= input('For solver type: \n 1 for Gauss Jaccobi \n 2 for Gauss Siedel \n 3 for SOR n'); 
jaccobi_iter=1;% cannot be kept inside the time loop as for every new time loop the counter will start from 1
gs_iter=1;
sor_iter=1;
tic;

% computing- time loop:
for nt = 1:1600
    % Jaccobi iteration
    if iter_solver == 1
    term1 = 1/(1+2*k1+2*k2);
    term2 = k1*term1;
    term3 = k2*term1;
    error = 1e9;% error outside this loop then for next time loop error will be less than tolerance and program wont procede
    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 2 max
       Told=T;
       jaccobi_iter = jaccobi_iter +1;
    end
   end
   

%% Gauss siedel iteration scheme
if iter_solver == 2
    term1 = 1/(1+2*k1+2*k2);
    term2 = k1*term1;
    term3 = k2*term1;
    error = 1e9;
    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 2 max
       Told=T;
    gs_iter= gs_iter +1;
    end
end

%% SOR iteration scheme
if iter_solver == 3
    term1 = 1/(1+2*k1+2*k2);
    term2 = k1*term1;
    term3 = k2*term1;
    error = 1e9;
    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 2 max
       Told=T;
    sor_iter= sor_iter +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);
set(gca,'YDIR','reverse');
if iter_solver == 1
title_text= sprintf('unsteady state heat conduction using Gauss Jacobi method \n No. of Iterations: %d \n Computation Time: %0.3f seconds\n CFL_num = %0.3f & dt = %0.4f',jaccobi_iter,time_required,CFL_num, dt);
title(title_text);
saveas(figure(1), 'unsteady_Implict_Jaccobi_dt=0.001_cfl=0.28.jpeg','jpeg');
else if iter_solver == 2
title_text= sprintf('unsteady state heat conduction using Gauss Siedel method \n No. of Iterations: %d \n Computation Time: %0.3f seconds\n CFL_num = %0.3f & dt = %0.4f',gs_iter,time_required,CFL_num,dt);
title(title_text);
saveas(figure(1), 'unsteady_Implict_Gauss_Siede_dt=0.001_cfl=0.28.jpeg','jpeg');
else if iter_solver == 3
title_text= sprintf('unsteady state heat conduction using SOR method \n No. of Iterations: %d \n Computation Time: %0.3f seconds \n CFL_num = %0.3f & dt = %0.4f',sor_iter,time_required,CFL_num, dt);
title(title_text);
saveas(figure(1), 'unsteady_Implict_SOR_dt=0.001_cfl=0.28.jpeg','jpeg');
    end
    end
end

xlabel('xx');
ylabel('yy');

Figure: Unsteady State - Explicit Scheme Jacobi

Inference: It is clear from the results that the solution has good convergence and statble since CFL value is less than 0.5. But Implicit schemes are unconditionally stable even when CFL>0.5 which is discussed in part 2 of challenge.

Figure: Unsteady State - Explicit Scheme Gauss Siedel

Inference: It is clear from the results that the solution has good convergence and statble since CFL value is less than 0.5. But Implicit schemes are unconditionally stable even when CFL>0.5 which is discussed in part 2 of challenge.

Figure: Unsteady State - Explicit Scheme SOR

Inference: It is clear from the results that the solution has good convergence and statble since CFL value is less than 0.5. But Implicit schemes are unconditionally stable even when CFL>0.5 which is discussed in part 2 of challenge. Convergence rate is highest compared to other 2 shcemes.

Convergence rate is of the order SOR>>Gauss Siedel >> Jacobi

3. Steady State 2D Head Conduction

STEADY STATE

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

MATLAB CODE

% MATLAB Program to generate temperature profile of 2D Transient State Heat Conduction Equation 
%Steady State 
clear all;
close all;
clc;
%Initialize variables
Nx = 21;
Ny = 21;
x=linspace(0,1,Nx);
y=linspace(0,1,Ny);
dx=1/(Nx-1);
dy=1/(Ny-1);
error = 1e9;
tol=1e-4;
k=2*((1/(dx^2))+ (1/(dy^2)));
%%Temperature field
T=298*ones(Nx,Ny);
T(1,:) = 600; % bottom
T(end,:)=900; % top 
T(:,1) = 400; % left 
T(:,end)=800; % right 
% Corner refining
T(1,1)=(600+400)/2;
T(1,end) = (600+800)/2;
T(end,1) = (400+900)/2;
T(end,end) = (800+900)/2;
omega = 1.2;
%copy T
Told=T;
%% Solving the 2D steady state system 
[xx,yy] = meshgrid(x,y);
iter_solver= input('For solver type: \n 1 for Gauss Jaccobi \n 2 for Gauss Siedel \n 3 for SOR n'); 
tic;
%Jacobi Iteration Scheme:
if iter_solver == 1
jaccobi_iter=1;
while (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
error = max(max(abs(Told-T))); % 2d problem so 2 max
Told=T;
jaccobi_iter= jaccobi_iter +1;
end
end
% Gauss siedel
if iter_solver == 2
gs_iter=1;
while (error > tol)
for i = 2:Nx-1
for j = 2:Ny-1
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)));
Told=T;
gs_iter= gs_iter +1;
end
end
%SOR
if iter_solver == 3
sor_iter=1;
while (error > tol)
for i = 2:Nx-1
for j = 2:Ny-1
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)));
Told=T;
sor_iter= sor_iter +1;
end
end
time_required=toc;
figure(1)
[C,h]=contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
set(gca,'YDIR','reverse');
if iter_solver == 1
title_text= sprintf('Steady state heat conduction using Gauss Jacobi method \n No. of Iterations: %d \n Computation Time: %0.3f seconds',jaccobi_iter,time_required);
title(title_text);
saveas(figure(1), 'unsteady_Explicit_Jaccobi.jpeg','jpeg');
else if iter_solver == 2
title_text= sprintf('Steady state heat conduction using Gauss Siedel method \n No. of Iterations: %d \n Computation Time: %0.3f seconds',gs_iter,time_required);
title(title_text);
saveas(figure(1), 'unsteady_Explicit_Gauss_siedel.jpeg','jpeg');
else if iter_solver == 3
title_text= sprintf('Steady state heat conduction using SOR method \n No. of Iterations: %d \n Computation Time: %0.3f seconds',sor_iter,time_required);
title(title_text);
saveas(figure(1), 'unsteady_Explicit_SOR.jpeg','jpeg');
    end
    end
end
xlabel('xx');
ylabel('yy');

Figure: Steady State - Gauss Jacobi

Inference: It is clear from the results that the solution has good convergence and statble.

Figure: Steady State - Gauss Siedel

Inference: It is clear from the results that the solution has good convergence and statble.

Figure: Steady State - SOR

Inference: It is clear from the results that the solution has good convergence and statble. for given omeg =1.2 SOR has higher convergence rate when compared to Gauss Jacobi and Siedel iterative solvers

Convergence rate is of the order SOR>>Gauss Siedel >> Jacobi.

3. CONCLUSION

The program to study steady and unsteady state 2D Heat conduction equation was developed and executed successfully.