MATLAB Scripting of steady and unsteady 2D heat conduction equation
Objective
Write a MATLAB Script to solver 2D Heat Diffusion Equation for Steady-state & Transient State using Jacobi, Gauss-seidel & Successive over-relaxation iterative method for Steady-state & Implicit and Explicit schemes for Transient state.
Input
2D, Plate with negligible thickness
Length of Plate = 1 meter, Width of Plate = 1 meter.
Convegence Criteria is 1e-4.
Boundary Condition:-
Top Boundary = 600 K, Bottom Boundary = 900 K, Left Boundary = 400 K, Right Boundary = 800 K
Uniform grid spacing along x & y direction with 10 nodes.
Thermal Diffusivity - 1.4 m2/s
Over-relaxation factor - 1.1
Timestep - 1000 (for transient simulation)
Discretization of 2D Diffusion Equation (Steady Form)
`(d^2T)/dx^2 + (d^2T)/dy^2 = 0`
Discretization scheme for diffusion term: Central Difference Scheme
`(T_(i+1,j) - 2T_(i,j) + T_(i-1,j))/(deltax)^2 + (T_(i,j+1) - 2T_(i,j) + T_(i,j-1))/(deltay)^2`
Iterarive method such as Jacobi, Gauss-seidel & Successive over-relaxation can be used to solve the above equation to obtain the steady state temperature field.
Jacobi Method - Calculates the temperature of nodes at the current iteration using temperature values of the previous iteration.
`T_(i,j) = (1/k)((Told_(i+1,j) + Told_(i-1,j))/(deltax)^2 + (Told_(i,j+1) + Told_(i,j-1))/(deltay)^2)`
T - Temperature of the current iteration
Told - Temperature of the previous iteration
`K = 2((deltax^2 + deltay^2)/(deltax^2deltay^2))`
Gauss-Seidel Method - Calculates the temperature of nodes at current iteration using temperature values of the latest updated iteration during the iterative procedures.
`T_(i,j) = (1/k)((Told_(i+1,j) + T_(i-1,j))/(deltax)^2 + (Told_(i,j+1) + T_(i,j-1))/(deltay)^2)`
T - Temperature of the current iteration
Told - Temperature of the previous iteration
`K = 2((deltax^2 + deltay^2)/(deltax^2deltay^2))`
Successive over-relaxation Method - Calculates the temperature of nodes at current iteration using temperature values of the previous iteration + Rate of change of temperature with respect to iteration multiplied with an over-relaxation factor.
`T_(i,j) = Told_(i,j) + alpha((dT)/(di))`
where `((dT)/(di)) = (T_(GS) - T_(old))/((n+1) - n)`
`alpha` - over-relaxation factor
n+1 - current iteration
n - old iteration
`T_(i,j)` - Temperature at current iteration
`Told_(i,j)` - Temperature at old iteration
`T_(GS)` - Temperature calculated at current iteration using Gauss-Seidel Method
Finally substituting `(dT)/(di)` in main formula and doing some re-arrangement, The final form we will get as follow
`T_(i,j) = Told_(i,j)(1 -alpha ) + alphaT_(GS)`
Discretization of 2D Diffusion Equation (Transient Form)
`(dT)/(dt) = alpha((d^2T)/dx^2 + (d^2T)/dy^2)`
`alpha` - Thermal Diffusivity
Explicit Method - Temperature of nodes at current time step is calculated based on previous time step nodes temperature.
`((T_(i,j) )^(n+1)-(T_(i,j) )^n)/(deltat)=alpha(((T_(i+1,j))^n - 2(T_(i,j) )^n+ (T_(i-1,j))^n)/(deltax)^2 + ((T_(i,j+1))^n - 2(T_(i,j))^n + (T_(i,j-1))^n)/(deltay)^2)`
Right hand side, all values are known so it can be directly solved for next timestep using time marching approach.
Using Jacobi Method for space marching, the final discretized equation is as follow:-
`(T_(i,j) )^(n+1)=(T_(i,j) )^n + (alphadeltat)/(deltax)^2((T_(i+1,j))^n - 2(T_(i,j) )^n+ (T_(i-1,j))^n) + (alphadeltat)/(deltay)^2((T_(i,j+1))^n - 2(T_(i,j))^n + (T_(i,j-1))^n)`
`(T_(i,j) )^n` - Temperature at previous time step
`(T_(i,j) )^(n+1)` - Temperature at current time step
Stability criteria:-
`(alphadeltat)/(deltax)^2+(alphadeltat)/(deltay)^2 <= 0.5` OR `(CFL_x + CFL_y <= 0.5)`
Implicit Method - Temperature of nodes at the current time step is calculated based on the current time step nodes temperature itself. This leads to a system of linear algebirc equations.
`((T_(i,j) )^(n+1)-(T_(i,j) )^n)/(deltat)=alpha(((T_(i+1,j))^(n+1) - 2(T_(i,j) )^(n+1)+ (T_(i-1,j))^(n+1))/(deltax)^2 + ((T_(i,j+1))^(n+1) - 2(T_(i,j))^(n+1) + (T_(i,j-1))^(n+1))/(deltay)^2)`
Right hand side all values are unknown, Hence Iterative technique such as Jacobi, Gauss-seidel & Successive over-relaxation has to be used followed by time marching.
Using different iterative methods(Jacobi, Gauss-Seidel, SOR), Solution will end up with a different number of iterations to reach specified solution time(0.5 seconds in this case) due to the iterative process nature.
Jacobi Method
`(T_(i,j) )^(n+1)=(T_(i,j) )^n + (alphadeltat)/(deltax)^2((T_(i+1,j))^(n+1) - 2(T_(i,j) )^(n+1)+ (T_(i-1,j))^(n+1)) + (alphadeltat)/(deltay)^2((T_(i,j+1))^(n+1) - 2(T_(i,j))^(n+1) + (T_(i,j-1))^(n+1))`
`(alphadeltat)/(deltax)^2 - CFL_x`
`(alphadeltat)/(deltay)^2 - CFL_y`
`(T_(i,j) )^(n+1)=(T_(i,j) )^n + CFL_x((T_(i+1,j))^(n+1) - 2(T_(i,j) )^(n+1)+ (T_(i-1,j))^(n+1)) + CFL_y((T_(i,j+1))^(n+1) - 2(T_(i,j))^(n+1) + (T_(i,j-1))^(n+1))`
`(T_(i,j) )^(n+1) (1 + 2CFL_x + 2CFL_y) = (T_(i,j) )^n + CFL_x((T_(i+1,j))^(n+1)+ (T_(i-1,j))^(n+1)) + CFL_y((T_(i,j+1))^(n+1) + (T_(i,j-1))^(n+1))`
`(T_(i,j) )^(n+1) = 1/((1 + 2CFL_x + 2CFL_y))((T_(i,j) )^n + CFL_x((T_(i+1,j))^(n+1)+ (T_(i-1,j))^(n+1)) + CFL_y((T_(i,j+1))^(n+1) + (T_(i,j-1))^(n+1)))`
`1/((1 + 2CFL_x + 2CFL_y)) - term_1`
`CFL_x*term_1 - term_2`
`CFL_y*term_1 - term_3`
We can write final equation for jacobi method as:-
`(T_(i,j) )^(n+1) = term_1(Tprev_(i,j) )^n + term_2((Told_(i+1,j))^(n+1)+ (Told_(i-1,j))^(n+1)) + term_3((Told_(i,j+1))^(n+1) + (Told_(i,j-1))^(n+1)))`
`(T_(i,j) )^(n+1)` - Temperature at current time step
`(Tprev_(i,j) )^n` - Temperature at previous time step(This will remain constant with in an iterative loop and will be updated for next time loop)
`(Told_(i,j))^(n+1)` - Temperature at previous iteration(This will change with in an iterative loop)
For the Gauss-Seidel method above equation can be used but with the latest updated temperature values as shown below
`(T_(i,j) )^(n+1) = term_1(Tprev_(i,j) )^n + term_2((Told_(i+1,j))^(n+1)+ (T_(i-1,j))^(n+1)) + term_3((Told_(i,j+1))^(n+1) + (T_(i,j-1))^(n+1)))`
For the SOR method, We will use equation as follow:-
`T_(i,j) = Told_(i,j)(1 -alpha ) + alphaT_(GS)`
`T_(GS) -`Temperature calculated at current iteration using Gauss-Seidel Method
`alpha` - over-relaxation factor
MATLAB Script
I have written separate code for different type of simulations i.e. Steady simulation with Jacobi Method, Steady simulation with Gauss-Seidal Method, Steady simulation with SOR Method, Transient simulation using Explicit Method, Transient simulation using Implicit Method.
I have written a main code where I have defined all the inputs and I have scripted user defined input to select the type of simulation steady or unsteady, type of solver Jacobi, Gauss-Seidal, SOR, Explicit or Implicit.
Main Code
clc
clear all
close all
L = 1; %length of plate
W = 1; %width of plate
nx = 10; %mesh points along x direction
ny = 10; %mesh points along y direction
nt = 500; %number of time step for transient simulation
dt = 0.001; %time step size
%generating the grid
x = linspace(0,L,nx);
y = linspace(0,W,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
error = 5e9; %initial error
convergence_criteria = 1e-4; %convergence criteria
T = 300*ones(nx,ny); %initializing the temeperature field
%defining boundary condition
T_top = 600;
T_bottom = 900;
T_right = 800;
T_left = 400;
T(:,1) = T_left;
T(:,end) = T_right;
T(1,:) = T_top;
T(end,:) = T_bottom;
%corner nodes temperature
T(1,1) = (T_top + T_left)/2;
T(end,1) = (T_bottom + T_left)/2;
T(1,end) = (T_top + T_right)/2;
T(end,end) = (T_bottom + T_right)/2;
%Selectthe solver
type = sprintf("Steady Simulation - 1 nTransient Simulation - 2n");
disp(type)
simulation_type = input("Select the simulation type: 1 or 2: ");
if simulation_type == 1 %steady
solver = sprintf("nJacobi Solver - 1 nGauss-Seidel Solver - 2nSOR Solver - 3n");
disp(solver)
solver_type = input("Select the solver: 1 or 2 or 3: ");
if solver_type == 1 %steady jacobi
jacobi_steady_solver(T,nx,ny,dx,dy,error,convergence_criteria,x,y)
elseif solver_type == 2 %steady gauss-seidel
gauss_seidel_steady_solver(T,nx,ny,dx,dy,error,convergence_criteria,x,y)
elseif solver_type == 3 %steady SOR
SOR_steady_solver(T,nx,ny,dx,dy,error,convergence_criteria,x,y)
end
elseif simulation_type == 2 %transient
solver = sprintf("nExplicit Solver - 1 nImplicit Solver - 2n");
disp(solver)
solver_type = input("Select the solver: 1 or 2: ");
if solver_type == 1 %explicit method
jacobi_unsteady_explicit_solver(T,nx,ny,dx,dy,dt,nt,x,y,convergence_criteria)
elseif solver_type == 2 %implicit method
sub_solver = sprintf("nJacobi Solver - 1 nGauss-Seidel Solver - 2nSOR Solver - 3n");
disp(sub_solver)
sub_solver_type = input("Select the sub solver: 1 or 2 or 3: ");
if sub_solver_type == 1 %implicit jacobi method
jacobi_unsteady_implicit_solver(T,nx,ny,dx,dy,dt,nt,x,y,convergence_criteria,error)
elseif sub_solver_type == 2 %implicit gauss-seidel method
gauss_seidel_unsteady_implicit_solver(T,nx,ny,dx,dy,dt,nt,x,y,convergence_criteria,error)
elseif sub_solver_type == 3 %implicit SOR method
SOR_unsteady_implicit_solver(T,nx,ny,dx,dy,dt,nt,x,y,convergence_criteria,error)
end
end
endCode for Steady-State Jacobi Method
function out = jacobi_steady_solver(T,nx,ny,dx,dy,error,convergence_criteria,x,y)
T_old = T; %temperature of previous iteration
jacobi_iteration = 1; %initial iteration
k = (2*(dx^2+dy^2))/((dx^2)*(dy^2));
%solving
while (max(error)>convergence_criteria)
for i = 2:nx-1 %march along x direction
for j = 2:ny-1 %march along y direction
T(i,j) = (1/k)*((T_old(i+1,j)+T_old(i-1,j))/(dx^2) + (T_old(i,j+1)+T_old(i,j-1))/(dy^2));
end
end
error = (max(abs(T_old - T))); %calculating the error
T_old = T; %updating temperature field
jacobi_iteration = jacobi_iteration + 1; %updating the iteration no.
end
%Ploting temperature contour
figure(1);
contourf(x,y,T);
clabel(contourf(x,y,T));
colorbar;
colormap(jet);
xlabel('X-Axis');
ylabel('Y-Axis');
title_text = sprintf('2D Steady Heat Diffusion, Jacobi Solver, Iterations = %d',jacobi_iteration);
title(title_text);
iter = sprintf("Number of Iterations to converge: %d",jacobi_iteration);
disp(iter)
set(gca,'ydir','reverse')
endCode for Steady-State Gauss-Seidel Method
function out = gauss_seidel_steady_solver(T,nx,ny,dx,dy,error,convergence_criteria,x,y)
T_old = T; %temperature of previous iteration
gauss_seidel_iteration = 1; %initial iteration
k = (2*(dx^2+dy^2))/((dx^2)*(dy^2));
%solving
while (max(error)>convergence_criteria)
for i = 2:nx-1 %march along x direction
for j = 2:ny-1 %march along y direction
T(i,j) = (1/k)*((T_old(i+1,j)+T(i-1,j))/(dx^2) + (T_old(i,j+1)+T(i,j-1))/(dy^2));
end
end
error = (max(abs(T_old - T))); %calculating the error
T_old = T; %updating temperature field
gauss_seidel_iteration = gauss_seidel_iteration + 1; %updating the iteration no.
end
%Ploting temperature contour
figure(1);
contourf(x,y,T);
clabel(contourf(x,y,T));
colorbar;
colormap(jet);
xlabel('X-Axis');
ylabel('Y-Axis');
title_text = sprintf('2D Steady Heat Diffusion, Gauss-Seidel Solver, Iterations = %d',gauss_seidel_iteration);
title(title_text);
iter = sprintf("Number of Iterations to converge: %d",gauss_seidel_iteration);
disp(iter)
set(gca,'ydir','reverse')
endCode for Steady-State SOR Method
function out = SOR_steady_solver(T,nx,ny,dx,dy,error,convergence_criteria,x,y)
T_old = T; %temperature of previous iteration
SOR_iteration = 1; %initial iteration
k = (2*(dx^2+dy^2))/((dx^2)*(dy^2));
alpha = 1.6; %Over-relaxation factor
%solving
while (max(error)>convergence_criteria)
for i = 2:nx-1 %march along x direction
for j = 2:ny-1 %march along y direction
T(i,j) = T(i,j)*(1-alpha) + (alpha*((1/k)*((T_old(i+1,j)+T(i-1,j))/(dx^2) + (T_old(i,j+1)+T(i,j-1))/(dy^2))));
end
end
error = (max(abs(T_old - T))); %calculating the error
T_old = T; %updating temperature field
SOR_iteration = SOR_iteration + 1; %updating iteration no.
end
%Ploting temperature contour
figure(1);
contourf(x,y,T);
clabel(contourf(x,y,T));
colorbar;
colormap(jet);
xlabel('X-Axis');
ylabel('Y-Axis');
title_text = sprintf('2D Steady Heat Diffusion, SOR Solver, Iterations = %d',SOR_iteration);
title(title_text);
iter = sprintf("Number of Iterations to converge: %d",SOR_iteration);
disp(iter)
set(gca,'ydir','reverse')
endCode for Transient Explicit Method
function out = jacobi_unsteady_explicit_solver(T,nx,ny,dx,dy,dt,nt,x,y,convergence_criteria)
T_old = T; %temperature of previous time step
alpha = 1.4; %thermal diffusivity
CFL_x = ((2*alpha*dt)/dx^2); %CFL along x
CFL_y = ((2*alpha*dt)/dy^2); %CFL along y
solution_time = dt; %tracking the solution time
for k = 1:nt %time marching
for i = 2:nx-1 %space marching along x
for j = 2:ny-1 %space marching along y
T(i,j) = T_old(i,j) + CFL_x*(T_old(i+1,j)-2*T_old(i,j)+T_old(i-1,j)) + CFL_y*(T_old(i,j+1)-2*T_old(i,j)+T_old(i,j-1));
end
end
error = max(max(abs(T_old - T))); %calculating the error
T_old = T; %updating the temperature
%ploting temperature contour
figure(1);
contourf(x,y,T);
clabel(contourf(x,y,T));
colorbar;
colormap(jet);
xlabel('X-Axis');
ylabel('Y-Axis');
solution_time = dt + solution_time;
title_text = sprintf('2D Explicit Unsteady Heat Diffusion, Jacobi, Solution Time: %.3f s',solution_time);
title(title_text);
set(gca,'ydir','reverse')
if (error < convergence_criteria);
disp("Simulation has reached steady state!")
break
end
endCode for Transient Implicit Jacobi Method
function out = jacobi_unsteady_implicit_solver(T,nx,ny,dx,dy,dt,nt,x,y,convergence_criteria,error)
T_old = T; %previous iteration temperature
T_prev = T; %previous time step temperature
alpha = 1.4; %thermal diffusivity
CFL_x = ((2*alpha*dt)/dx^2); %CFL along x
CFL_y = ((2*alpha*dt)/dy^2); %CFL along y
term_1 = (1 + 2*CFL_x + 2*CFL_y)^(-1);
term_2 = term_1*CFL_x;
term_3 = term_1*CFL_y;
solution_time = dt; %tracking the solution time
jacobi_iteration = 1;
for k = 1:nt %time marching
while (error > convergence_criteria) %Jacobi Iterative Method
for i = 2:nx-1 %space marching along x
for j = 2:ny-1 %space marching along y
H = (T_old(i+1,j) + T_old(i-1,j));
V = (T_old(i,j+1) + T_old(i,j-1));
T(i,j) = ((T_prev(i,j)*term_1) + (term_2*H) + (term_3*V));
end
end
error = max(max(abs(T_old - T))); %calculating the error
T_old = T; %updating the temperature for next iteration
jacobi_iteration = 1 + jacobi_iteration;
end
error = 1; %re-assigning the error
T_prev = T; %updating the temperature for the next time step
%ploting temperature contour
figure(1);
contourf(x,y,T);
clabel(contourf(x,y,T));
colorbar;
colormap(jet);
xlabel('X-Axis');
ylabel('Y-Axis');
solution_time = dt + solution_time;
title_text = sprintf('2D Implicit Unsteady Heat Diffusion, Jacobi, Solution Time: %.3f s',solution_time);
title(title_text);
set(gca,'ydir','reverse')
end
endCode for Transient Implicit Gauss-Seidel Method
function out = gauss_seidel_unsteady_implicit_solver(T,nx,ny,dx,dy,dt,nt,x,y,convergence_criteria,error)
T_old = T; %previous iteration temperature
T_prev = T; %previous time step temperature
alpha = 1.4; %thermal diffusivity
CFL_x = ((2*alpha*dt)/dx^2); %CFL along x
CFL_y = ((2*alpha*dt)/dy^2); %CFL along y
term_1 = (1 + 2*CFL_x + 2*CFL_y)^(-1);
term_2 = term_1*CFL_x;
term_3 = term_1*CFL_y;
solution_time = dt; %tracking the solution time
gauss_seidel_iteration = 1;
for k = 1:nt %time marching
while (error > convergence_criteria) %Gauss-Seidel Iterative Method
for i = 2:nx-1 %space marching along x
for j = 2:ny-1 %space marching along y
H = (T_old(i+1,j) + T(i-1,j));
V = (T_old(i,j+1) + T(i,j-1));
T(i,j) = ((T_prev(i,j)*term_1) + (term_2*H) + (term_3*V));
end
end
error = max(max(abs(T_old - T))); %calculating the error
T_old = T; %updating the temperature for next iteration
gauss_seidel_iteration = 1 + gauss_seidel_iteration;
end
error = 1; %re-assigning the error
T_prev = T; %updating the temperature for the next time step
%ploting temperature contour
figure(1);
contourf(x,y,T);
clabel(contourf(x,y,T));
colorbar;
colormap(jet);
xlabel('X-Axis');
ylabel('Y-Axis');
solution_time = dt + solution_time;
title_text = sprintf('2D Implicit Unsteady Heat Diffusion, Gauss-Seidel, Solution Time: %.3f s',solution_time);
title(title_text);
set(gca,'ydir','reverse')
end
endCode for Transient Implicit SOR Method
function out = SOR_unsteady_implicit_solver(T,nx,ny,dx,dy,dt,nt,x,y,convergence_criteria,error)
T_old = T; %previous iteration temperature
T_prev = T; %previous time step temperature value
alpha = 1.4; %thermal diffusivity
omega = 1.1; %over-relaxation factor
CFL_x = ((2*alpha*dt)/dx^2); %CFL along x
CFL_y = ((2*alpha*dt)/dy^2); %CFL along y
term_1 = (1 + 2*CFL_x + 2*CFL_y)^(-1);
term_2 = term_1*CFL_x;
term_3 = term_1*CFL_y;
solution_time = dt; %tracking the solution time
SOR_iteration = 1;
for k = 1:nt %time marching
while (error > convergence_criteria) %SOR Iterative Method
for i = 2:nx-1 %space marching along x
for j = 2:ny-1 %space marching along y
H = (T_old(i+1,j) + T(i-1,j));
V = (T_old(i,j+1) + T(i,j-1));
T_gauss_seidel = ((T_prev(i,j)*term_1) + (term_2*H) + (term_3*V));
T(i,j) = T_old(i,j)*(1-omega) + (T_gauss_seidel*omega);
end
end
error = max(max(abs(T_old - T))); %calculating the error
T_old = T; %updating the temperature for next iteration
SOR_iteration = 1 + SOR_iteration;
end
error = 1; %re-assigning the error
T_prev = T; %updating the temperature for the next time step
%ploting temperature contour
figure(1);
contourf(x,y,T);
clabel(contourf(x,y,T));
colorbar;
colormap(jet);
xlabel('X-Axis');
ylabel('Y-Axis');
solution_time = dt + solution_time;
title_text = sprintf('2D Implicit Unsteady Heat Diffusion, SOR, Solution Time: %.3f s',solution_time);
title(title_text);
set(gca,'ydir','reverse')
end
endResults for Steady Cases
| Solver Type | Number of iterations to reach convergence |
| Jacobi | 208 |
| Gauss-Seidel | 112 |
| SOR (over relexation factor = 1.1) | 38 |
SOR takes the least iterations to reach a steady-state solution whereas Jacobi takes the maximum number of iterations.
Results for Transient Cases
Explicit Scheme with Jacobi Solver
The temperature field reaches a steady-state in 0.228 seconds.
Implicit Scheme with Jacobi, Gauss-Seidel & SOR Solver
For the implicit method, the solution time is 0.5 seconds.
The temperature field has reached its steady state.
Jacobi method takes a maximum of 2318 iterations to reach 0.5 seconds of solution time.
Gauss-Seidel method takes 1790 iterations to reach 0.5 seconds of solution time.
SOR method takes a minimum of 1548 iterations to reach 0.5 seconds of solution time.