Simulation of Steady and Unsteady heat conduction problem with a custom build 2D solver using MATLAB
Objectives of the project:-
- Understanding the temperature distribution along the 2D body under both steady and unsteady state conditions.
- Writing a general Matlab code to 2D Heat conduction solver which includes Steady and unsteady states with finite differences method.
- In the generic code, three different iterative methods such as Jacobi, Gauss-Seidel and Successsiveover relaxation methods are used in both explicit and implicit ways to solve, System of linearly uncoupled and coupled equations respectively for both steady and unsteady states.
Literature:-
Heat transfer from one place to another by three methods, namely conduction in solids, convection in liquids and gases and radiation, through anything that will allow the radiation to pass. If there is a temperature difference in a system, heat will always move from higher to lower temperatures.
Heat conduction occurs when two objects are at different temperatures are in contact with each other, until both bodies reached the same temperature. Conduction is the movement of heat through a substance by the collision of molecules. At the place where the two object touch, the faster-moving molecules of the warmer object collide with the slower moving molecules of the cooler object. As they collide, the faster molecules give up some of their energy to the slower molecules. The slower molecules gain more thermal energy and collide with other molecules in the cooler object. This process continues until heat energy from the warmer object spreads throughout the cooler object. Thermal conductivity is the ability to conduct the heat.
General heat conduction equation is
`∇.(k∇T) +Q\'\' +βT(Dp)/(Dt)+μ∅= ρCp(DT)/(Dt)`
Following are the assumptions made to formulate setup:-
- Heat conduction in the Z direction is negligible.
- Viscous dissipation (μ∅) is neglected.
- Compressibility effects are neglected as it\'s an incompressible fluid which means volumetric dilation term (βT(Dp)/(Dt))neglected.
- Source or heat generation(Q\'\') within the body is absent.
- Heat is transferring through molecular vibrations only, i.e. conduction.
- Thermal conductivity is Independent of temperature, i.e K = constant
The final reduced form of the equations is of the form
`(∂T)/(∂t)=α((∂^2T)/(∂x^2) +(∂^2T)/(∂x^2))`
where `α = k/(ρCp)`
Finite difference approach:-
Applying Forward in Time and Central in space (FTCS) scheme which means Temporal terms(parabolic in time) discrete in Forwarding difference scheme whereas Spatial terms(Elliptic in space) discrete in Central difference scheme.
lets consider the numerical stencil at gridpoints(x,y) be (i-1,j),(i,j),(i+1,j),(i,j-1),(i,j+1).
at (i,j)th node at n th time step,
`((delT)/(delt))=(T^(n+1)(i,j) - T^(n)(i,j) )/dt`
`alpha(((del^2T)/(delx^2)) + ((del^2T)/(dely^2))) =alpha((T^(n)(i-1,j) -2T^(n)(i,j) +T^(n)(i+1,j) )/dx)+alpha((T^(n)(i,j-1) -2T^(n)(i,j) +T^(n)(i,j-1) )/dy)`
`(T^(n+1)(i,j) - T^(n)(i,j) )/dt = alpha((T^(n)(i-1,j) -2T^(n)(i,j) +T^(n)(i+1,j) )/dx)+alpha((T^(n)(i,j-1) -2T^(n)(i,j) +T^(n)(i,j-1) )/dy)`
`T^(n+1)(i,j) = T^(n)(i,j) +R1(T^(n)(i-1,j) -2T^(n)(i,j) +T^(n)(i+1,j)) +R2(T^(n)(i,j-1) -2T^(n)(i,j) +T^(n)(i,j-1) )`
where R1 = alpha*dt/dx^2 , R2 = alpha*dt/dy^2 (Fourier numbers in x and y-direction)
Methods of solving
these equations can be represented in Both Explicit and Implicit forms so that uncoupled and coupled Linear equations formed respectively, which can be solved with many methods
Such as 1) Direct analytical solving
2) Iterative - Indirect techniques, Jacobi, Gauss- Siedal, Successive overrelaxation methods.
In this scheme, Explicit forms are Conditionally stable with the limitation on Time step, on the other hand, Implicit method is Unconditionally stable with no effect of time step on a solution.
Steps in FDM Approach:-
- Step 1:- Discretize the spatial domain into finite gridpoints with grid sizes in x and y-direction.
- Step 2:- Choosing the appropriate governing p.d.e to represent the problem.
- Step 3:- Define the initial and boundary conditions of the Situation. Boundary conditions can be
1)Dirichlet Boundary conditions - The value of the dependent variable specified.
2) Neuman boundary conditions - The value of the gradient dependent variable specified.
3) Mixed boundary conditions - The value of the dependent variable is expressed as a function of the gradient.
4) Periodic boundary conditions - The value of the cyclic boundary conditions specified.
- Step 4:- Choose the Appropriate scheme to apply based on time step dependency with factors.
1) Consistency - Truncation error dependency over grid size and step size.
2) Stability - Round off error propagation parameter.
3) Convergence - The error in exact and numerical solution.
Whatever the scheme chose it must be met following requirements
1)Conservativeness.
2) Boundedness.
3) Transportiveness.
- Step 5:- Solve the linearly coupled or uncoupled equations using appropriate Numeriwith converging and tolerance limits in the setup.
- Step 6:- Plot the results of the solution.
Methodology:-
1) In the present study, setup is formulated to the specific geometry of 2D plate with (X, Y) as the dimensions in X and Y- directions.
2) Assumed that the object is set at the boundary conditions at t = 0 where the internal grid points are T = 0 deg.
3) wrote the Generic solver code(Given at end) In MATLAB with User-defined dimensions, Grid sizing, material properties, and Type of study and method of solving.
4) In the piece of code provide the steady-state - explicit method with Jacobi, Gauss, SOR
methods were used to solve the problem.
5)In the unsteady state, Both explicit and implicit methods are used to Solve the problem where Jacobi, Gauss, SOR methods are employed individually.
for the SOR relaxation parameter (0<omega<2) inputted as the optimal value.,
`omega.optimal = 2/((1+sin(pi*dx))`= 2/(1+sin(pi*dx));
Validation of the solver:-
Physical properties:
Material: Copper
Density: 896Kg/m^3
Thermal Conductivity: 400 W/m.k
Specific Heat : 376.812 J/kg.k
Dimensions:
length: 1mt
Breath: 1mt
Boundary Conditions: (Dirichlet Boundary conditions)
Temperature on top surface = 600 k
Temperature on bottom surface = 900 k
Temperature on left surface = 400 k
Temperature on right surface = 800 k
Tolerance limit: 1e-4
Results and plots:-
(A) Steady state analysis:-
Time-dependent terms are neglected in steady state analysis.
`((∂^2T)/(∂x^2) +(∂^2T)/(∂x^2)) = 0`
Case 1:- Number of grid points in the X-direction = 10
Number of grid points in the Y-direction = 10
case:- 2
Number of grid points in the X-direction = 20
Number of grid points in the Y-direction = 20
Comparison between different number of grid points:-
(a) For Nx= Ny =10
(b) For Nx= Ny =20
(c) For Nx= Ny =30
Observations:-
1) With the change in grid size number of iterations increases in multiple folds.
2) Computational time heavily depends on the grid size.
3) In comparing three Linear solving techniques SOR suits and fastest solver.
4)The steady state temperature distribution is grid size independent.
(B) unsteady state analysis:-
Number of grid points in the x-direction: 10
Number of grid points in y-direction: 10
from the steady-state analysis, it is clear that an increase in the number of grid points increases the Computational time rapidly.
Here comparison needs to make between the time step chosen and Convergence rate along with computational time.
Explicit Formulation:-
As Explicit formulation is conditionally stable, for
CFL <= .25, it\'s important to understand the effect of the change in time step on the solution in various methods.
Maximum stable time step for this configuration is 2.1108 sec
When time step = 1 sec.
And results captured at various Times such as T = 1,5,10,25,50sec.
For T = 50sec results are given below.
(A) Comparing the results of T = 1, 5,10,25,50 at time step = 1
(B) Comparing the results of T = 1, 5,10,25,50 at time step = 2
| Steady-state Explicit Method | ||||||
| Method | Nx=Ny=10 | Nx=Ny=20 | Nx=Ny=30 | |||
| No. of iterations | Computation Time | No.of iterations | Computation Time | No. of iterations | Computation Time | |
| Jacobi method | 217 | 1.17 | 869 | 9.72 | 1884 | 50.673 |
| Gauss-seidal method | 116 | 0.96 | 462 | 6.092 | 1003 | 27.758 |
| SOR method | 32 | 0.84 | 60 | 1.94 | 90 | 4.806 |
Observations: -
1.As the time step increased the rate of convergence increased as well\'
2.In both the cases for the same number of iterations SOR method converging to steady state compared to others.
Implicit Formulation:-
As implicit formulation is unconditionally stable, to understand the effect of the change in time step is given below.
Maximum stable time step for this configuration is 2.1108 sec
When time step = 1 sec.
And results captured at various Times , T = 1,5,10,25,50sec.
For T = 50sec results are given below.
Comparing the results at T= 1,5,10,25,50sec for time step 1 and 2 to see the time step
effect on the implicit method.
(A) Comparing the results at T= 1,5,10,25,50sec for time step -1sec:
(B) Comparing the results at T= 1,5,10,25,50sec for time step -2 sec:
observation:-
As the implicit method is unconditionally stable the effect of time steps on the solution is negligible, which can be observed from the above plots.
Whereas the Number of iterations and Computation time have its impact which can be seen in the below table.
| Implicit Method ( dt = 1sec) | ||||||
| Time (T) | Jacobi | Gauss-seidal | SOR | |||
| No .of iterations | Computation Time | No .of iterations | Computation Time | No .of iterations | Computation Time | |
| 1 | 12 | 0.786 | 9 | 0.79 | 9 | 0.828 |
| 5 | 60 | 0.878 | 45 | 0.862 | 45 | 0.87 |
| 10 | 118 | 0.986 | 90 | 0.895 | 89 | 1.018 |
| 25 | 283 | 1.312 | 213 | 1.172 | 209 | 1.43 |
| 50 | 558 | 1.934 | 413 | 1.602 | 386 | 1.994 |
| Implicit Method ( dt = 2sec) | ||||||
| Time (T) | Jacobi | Gauss-seidal | SOR | |||
| No .of iterations | Computation Time | No .of iterations | Computation Time | No .of iterations | Computation Time | |
| 1 | 1 | 0.779 | 1 | 0.762 | 1 | 0.769 |
| 5 | 36 | 0.806 | 25 | 0.795 | 25 | 0.833 |
| 10 | 88 | 0.916 | 61 | 0.846 | 61 | 0.908 |
| 25 | 207 | 1.176 | 145 | 1.035 | 139 | 1.167 |
| 50 | 420 | 1.582 | 291 | 1.316 | 269 | 1.584 |
as time step increases computation time and Number of iterations are decreasing.
A) Matlab Code for the solver:-
clear all;
close all;
clc;
tic;
disp(sprintf(\" \"))
display(\" 2D heat conduction solver for steady and unsteady state \")
display(\" - Using Finite Difference FTCS Scheme \")
disp(sprintf(\" \"))
display(sprintf(\"************Domain defining and discretization**************\"))
lx = input(\'Length along x direction: \'); %length of X domain
ly = input(\'Length along y direction: \'); %length of Y domain
nx = input(\'Number of grids in X direction: \'); %Number of grids in X Domain
ny = input(\'Number of grids in Y direction: \');; %Number of grids in Y Domain
x = linspace(0,lx,nx); %Discretizarion of the X damain with equal step size
y = linspace(0,lx,ny); %Discretizarion of the Y damain with equal step size
dx = lx/nx; %Discretizedd step size in X direction
dy = ly/ny; %Discretizedd step size in Y direction
disp(sprintf(\"Discretizedd step size in X direction: %d \",dx));
disp(sprintf(\"Discretizedd step size in Y direction: %d \",dy));
%%%%%%%%%%%%%%%%%%%%%%%% Boundary Conditions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
display(sprintf(\"****************Boundary conditions*************************\"))
Ttop = input(\'Temperature on the Top surface(K): \');
Tbottom = input(\'Temperature on the Bottom surface(K): \');
Tleft = input(\'Temperature on the Left surface(K): \');
Tright = input(\'Temperature on the Right surface(K): \');
Temp = [Ttop ,Tbottom, Tleft ,Tright];
%%%%%%%%%%%%%%%%%%% Material properties %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
display(\'*********************Material properties*****************************\')
rho = input(\'Density of the Material : \');
k = input (\'Thermal conductivity of the material: \');
Cp = input(\'Specific heat of the material : \');
alpha = k/(rho*Cp) ; % Thermal diffusivity,(m^2/s);
disp(sprintf(\'Thermal diffusivity of the materail: %.5f \',alpha));
deltaT_stable_max =((dx^2*dy^2)/(dx^2+dy^2))*(1/(2*alpha));
%%%%%%%%%%%%%%%%%%%%% Convergence parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$
err = 9e9; %%% Error Limit
tol =input(\'Enter tolerance limit:\'); %%% Tolerance Limit
display(\'************Defining Time dependency and study situation************ \')
display(\'For Steady state analysis: 1, For transient state analysis: 2 \')
situation =input(\'Time dependency(Steady or unsteady):\');
if situation == 2;
%%%%%%%%%%%%%%%%%%%%%%%% Time setting for Transient analysis %%%%%%%%%%%%%%%%%%%
display(\'****************Time limit and Time step parameters setting**********\')
T =input(\'Time to perform analysis (T):\'); %Stoping Time for Time loop,(s);
display(\'Solution with Explict method or implicit method.\')
display(\'For Explict method enter :0,For Implicit method enter :1\')
Approch =input(\'Explicit or Implicit way : \')
if Approch == 0;
display(\'Explicit Method is conditionally stable, for stability choose \')
disp(sprintf(\" time step less than maximum time step allowable for stability i.e:%.6fs\",deltaT_stable_max));
dt = input(\'Time step for the iteration(dt):\') %Time step for the iteration;
elseif Approch == 1;
disp(sprintf(\" Implicit Method is unconditionally stable \"))
dt = input(\'Time step for the iteration(dt):\') %Time step for the iteration;
else
disp(sprintf(\'Enter correct value by running program again.\'))
return
end
end;
display(\" \")
disp(sprintf(\'Solution can be solved solved with three methods.\' ))
display(\" \")
disp(sprintf(\'Choose any:Jacobi, Gauss-seidal, Successive overrelaxation.\' ))
display(\" \")
disp(sprintf(\'Jacobi method: 1 \\n Gauss-seidal method : 2 \\n SOR Method: 3\'))
method = input(\'Solving method setting:\') % 1- Jacobi,2-Gauss seidal ,3-SOR method
if method > 3;
disp(sprintf(\'No method is allocated to this value\'))
disp(sprintf(\'Enter correct value by running program again\'))
return
end
if (situation ==1)&& method <= 3; %%%% Steady state study
[Tcurrent,iterative_step,str_array,execution_time,err] = Steady_state(method,...
lx,ly,x,y,nx,ny,dx,dy,Temp,err,tol);
elseif (situation == 2)&& method <= 3; %%%% unsteady state study
[Tcurrent,iterative_step,execution_time,str_array,deltaT_stable_max]= Transient_state(Approch,method,...
lx,ly,nx,ny,x,y,dx,dy,Temp,tol,alpha,T,dt);
else
disp(sprintf(\'Enter correct value by running program again\'))
return
end
display(\" \")
display(\" \")
display(\'************************ Results ************************************\')
display(\" \")
display(\" \")
disp(sprintf(\'Length along x direction :%d\',lx));
disp(sprintf(\'Length along x direction :%d\',ly));
disp(sprintf(\'Density of material :%f\',rho));
disp(sprintf(\'Thermal conductivity of material :%d\',k));
disp(sprintf(\'Specific heat for the material :%d\',Cp));
if situation == 1;
disp(sprintf(\'Type of the problem : Steady state heat conduction\'));
else
disp(sprintf(\'Type of the problem : Unsteady state heat conduction\'));
end
disp(sprintf(\"Method of implementation : %s\",str_array));
disp(sprintf(\"Number of Iteration : %d\",iterative_step));
Execution_time = toc;
disp(sprintf(\"Computation Time : %d \",Execution_time));Code for Function Call of Steady State:-
function [Tcurrent,iterative_step,str_array,execution_time,err] = Steady_state(method,lx,ly,x,y,...
nx,ny,dx,dy,Temp,err,tol);
%%%%%%%%%%%%%%%%% Boundary conditions setup %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[xx,yy] = meshgrid(x,y); % Mesh Generation
Tcurrent= ones(nx,ny); % Initialize interanal temperature
%%% internal nodes
Tcurrent(end,:) = Temp(1); % Top face temperaturee setting
Tcurrent(1,:) = Temp(2); % Bottom face temperature setting
Tcurrent(:,1) = Temp(3); % left face temperature setting
Tcurrent(:,end) = Temp(4); % Right face temperature setting
%%% corner nodes
Tcurrent(1,1) = (Temp(1)+Temp(3))/2;
Tcurrent(1,end) = (Temp(1)+Temp(4))/2;
Tcurrent(end,1) = (Temp(3)+Temp(2))/2;
Tcurrent(end,end) = (Temp(4)+Temp(2))/2;
%%% Creation of Temperature at nth time step
Told = Tcurrent;
%%% Optimal value for Relaxation parameter
Wi = 2/(1+sin(pi*dx)); %Value of omega (0<omega<2);
%%%%%%%%%%%%%%%% solving for Temperature with various Methods %%%%%%%%%%%%%%%%%%
iterative_step = 1;
tic;
while (err > tol); %%Convergene loop initiation
for i = 2:nx-1;
for j = 2:ny-1;
if method == 1; %%% Jacobi method explicit method
Tcurrent(i,j) = (1/4)*(Told(i+1,j)+Told(i-1,j)...
+Told(i,j+1)+Told(i,j-1));
elseif method == 2; %%% Gauss- seidal explicit method
Tcurrent(i,j) = (1/4)*(Told(i+1,j)+Tcurrent(i-1,j)...
+Told(i,j+1)+Tcurrent(i,j-1));
elseif method ==3; %%% Successive over relaxation method
Tgauss(i,j)=((1/4)*(Told(i+1,j)+Tcurrent(i-1,j)...
+Told(i,j+1)+Tcurrent(i,j-1)));
Tcurrent(i,j) = Told(i,j)*(1-Wi)+Wi*Tgauss(i,j);
end
end
end
Error_matrix= max(abs(Told-Tcurrent));
err = max(Error_matrix);
iterative_step = iterative_step + 1;
Told = Tcurrent;
endwhile
%%%%%%%%%%%%%%%%%%%% plotting results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
[C,h] = contourf(xx,yy,Tcurrent);
colormap(jet); % Sets colors used for contour plot
clabel(C,h);
execution_time = toc;
H =colorbar; % Adds a scale to the plot
ylabel(H,\'Temperature in kelvin\',\'Linewidth\',8);
axis equal tight % Makes the axes have equal length
grid on;
xlabel(\'x(m)\');
ylabel(\'y(m)\');
set(gca,\'XTick\',0:.1:1) % Sets the x-axis tick mark locations
set(gca,\'YTick\',0:.1:1) % Sets the y-axis tick mark locations
if method == 1;
str_array= [\"Jacobi Explicit Method\"];
elseif method == 2 ;
str_array = [\"Gauss-seidal Explicit Method\"];
elseif method ==3 ;
str_array = [\"SOR Explicit Method\"];
end
formatSpec = \"2D Steady state Heat Conduction solver: %s \\n Iteration Counter: %d \\n Computation time: %0.3fs\";
title(sprintf(sprintf(formatSpec,str_array,iterative_step,...
execution_time)));Code For function Call of Unsteady state:-
function [Tcurrent,iterative_step,execution_time,str_array,deltaT_stable_max] = Transient_state(Approch,method,...
lx,ly,nx,ny,x,y,dx,dy,Temp,tol,alpha,T,dt);
%%%%%%%%%%%%%%%%%%%%% Boundary conditions setup %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Tcurrent= ones(nx,ny); % Initialize interanal temperature
[xx,yy] = meshgrid(x,y); % Mesh generation
%%% internal nodes
Tcurrent(end,:) = Temp(1); % Top face temperaturee setting
Tcurrent(1,:) = Temp(2); % Bottom face temperature setting
Tcurrent(:,1) = Temp(3); % left face temperature setting
Tcurrent(:,end) = Temp(4); % Right face temperature setting
%%% corner nodes--- Edge refining
Tcurrent(1,1) = (Temp(1)+Temp(3))/2;
Tcurrent(1,end) = (Temp(1)+Temp(4))/2;
Tcurrent(end,1) = (Temp(3)+Temp(2))/2;
Tcurrent(end,end) = (Temp(4)+Temp(2))/2;
%%% Creation of Temperature at nth time step
Told = Tcurrent;
T_Constant = Told;
%%% Optimal value for Relaxation parameter
Wi = 2/(1+sin(pi*dx));
%%%%%%%%%%%%%%%%%% Stability criterion %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
deltaT_stable_max =((dx^2*dy^2)/(dx^2+dy^2))*(1/(2*alpha));
%Maximum stable time step,(s);
r1 = (alpha*dt)/(dx^2); %Fourier number for x direction.
r2= (alpha*dt)/(dy^2); %Fourier number for Y direction.
r = (1+(2*r1)+(2*r2))^-1;
R1 = r1*r;
R2 = r2*r;
R = r;
%%%%%%%%%%%%%%%%%%%%%%% solving for Temperature %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
iterative_step = 1;
for i =1:(T/dt); %%Time loop initiation
tic;
if Approch == 0;
%%%%%%%%%%%%%%%%%%%%%%%% Explicit Formulation %%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 2:nx-1;
for j = 2:ny-1;
if method == 1; %Jacobi Explicit method
Tcurrent(i,j) = Told(i,j) + r1*(Told(i+1,j)-(2*Told(i,j))+...
Told(i-1,j))+r2*(Told(i,j+1)-(2*Told(i,j))+Told(i,j-1));
elseif method == 2; %Gauss seigal Explicit method
Tcurrent(i,j) = Told(i,j) + r1*(Told(i+1,j)-(2*Told(i,j))+...
Tcurrent(i-1,j))+r2*(Told(i,j+1)-(2*Told(i,j))+Tcurrent(i,j-1));
elseif method == 3; % SOR Explicit method
Tgauss(i,j) = (Told(i,j) + r1*(Told(i+1,j)-(2*Told(i,j))+...
Tcurrent(i-1,j))+r2*(Told(i,j+1)-(2*Told(i,j))+Tcurrent(i,j-1)));
Tcurrent(i,j) =Told(i,j)*(1-Wi) + Wi*Tgauss(i,j);
end
end
end
Told = Tcurrent;
exectution_time= toc;
iterative_step = iterative_step + 1
elseif Approch == 1;
%%%%%%%%%%%%%%%%%%%% Implicit Formulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
err=6e9;
while (err > tol);
for i = 2:nx-1;
for j = 2:ny-1;
if method == 1; %Jacobi Implicit method
Tcurrent(i,j) = R*T_Constant(i,j) + R1*(Told(i+1,j)+...
Told(i-1,j))+ R2*(Told(i,j+1)+ Told(i,j-1));
elseif method == 2; % Gauss-seidal Implicit method
Tcurrent(i,j) = R*T_Constant(i,j) + R1*(Told(i+1,j)+...
Tcurrent(i-1,j))+ R2*(Told(i,j+1)+ Tcurrent(i,j-1));
elseif method == 3; % SOR Implicit method
Tgauss(i,j) = R*T_Constant(i,j) + R1*(Told(i+1,j)+...
Tcurrent(i-1,j))+ R2*(Told(i,j+1)+ Tcurrent(i,j-1));
Tcurrent(i,j) =Told(i,j)*(1-Wi) + Wi*Tgauss(i,j);
end
end
end
err = max(max(abs(Told-Tcurrent)));
Told = Tcurrent;
exectution_time= toc;
iterative_step = iterative_step + 1;
endwhile
T_Constant = Tcurrent;
end
endfor
%%%%%%%%%%%%%%%%%% Plotting results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
[C,h] = contourf(xx,yy,Tcurrent);
colormap(jet); % Sets colors used for contour plot
clabel(C,h);
grid on;
execution_time = toc;
H =colorbar; % Adds a scale to the plot
ylabel(H,\'Temperature in kelvin\',\'Linewidth\',8);
axis equal tight % Makes the axes have equal length
xlabel(\'x(m)\');
ylabel(\'y(m)\');
set(gca,\'XTick\',0:.1:1) % Sets the x-axis tick mark locations
set(gca,\'YTick\',0:.1:1) % Sets the y-axis tick mark locations
if Approch == 0;
if method == 1;
str_array= [\"Jacobi Explicit Method\"];
elseif method == 2 ;
str_array = [\"Gauss-seidal Explicit Method\"];
elseif method ==3 ;
str_array = [\"SOR Explicit Method\"]
end
elseif Approch == 1;
if method == 1;
str_array= [\"Jacobi Implicit Method\"];
elseif method == 2 ;
str_array = [\"Gauss-seidal Implicit Method\"];
elseif method ==3 ;
str_array = [\"SOR Implicit Method\"]
end
end
formatSpec = \"2D UnSteady state Heat Conduction solver: %s \\n Iteration Counter: %d.Time step:%.4fs \\n Computation time: %0.3fs. Maximum stable Time step: %0.4fs. Time :%.3fs\";
title(sprintf(sprintf(formatSpec,str_array,iterative_step,...
dt,execution_time,deltaT_stable_max,T)));