MATLAB code to solve for the 2D heat conduction equation in different schemes
Aim: To solve for the 2D heat conduction equation in Steady-state and Transient state in explicit and implicit methods using the iterative techniques.
Objectives:
- To write a code in MATLAB to solve for the 2D heat conduction equation in Steady-state for the given boundary conditions using the point iterative techniques.
- To write a code in MATLAB to solve for the 2D heat conduction equation in Transient state in the explicit method.
- To write a code in MATLAB to solve for the 2D heat conduction equation in Transient state in the implicit method using the point iterative techniques.
2D Steady-state heat conduction equation:
The equation is as follows:
`(del^2T)/(delx^2)+(del^2T)/(dely^2)=0`
The boundary conditions for all four sides of our domain are defined initially. Hence, we need to solve for the temperature for all the other points in the domain whose value has an acceptable error value.
The boundary conditions for our problem is as follows:
Top boundary = 600K,
Bottom boundary = 900K,
Left boundary = 400K,
Right boundary = 800K.
Since the boundary conditions are defined on all sides, we can use the central difference scheme that takes values from all four sides. Hence, the simplifies form the equation to solve using the central differencing scheme for the space derivatives is as follows:
`T_(i,j)=(dx^2dy^2)/(2(dx^2+dy^2))[(T_(i-1,j)+(T_(i+1,j)))/dx^2 +(T_(i,j-1)+T_(i,j+1))/dy^2 ]`
We predefine an error term and a guess value for the temperature and then use iterative solvers to arrive at a value for the temperature at each node that has an acceptable value for the error.
The predetermined guess value is taken as 300K in our analysis.
The above equation is solved using three different iterative techniques and their description is as follows:
1. Jacobi method: We use the following equation to solve for the temperature in an iterative loop until the error value, i.e., the absolute difference is less than 1e-4.
`T_(i,j)=(dx^2dy^2)/(2(dx^2+dy^2))[(T_(old(i-1,j))+(T_(old(i+1,j))))/dx^2 +(T_(old(i,j-1))+T_(old(i,j+1)))/dy^2 ]`
Where the subscript ‘old’ refers to the previous temperature at that particular point.
2. Gauss-Seidel method: The following is the equation used in an iterative manner until the error is acceptable.
`T_(i,j)=(dx^2dy^2)/(2(dx^2+dy^2))[(T_((i-1,j))+(T_(old(i+1,j))))/dx^2 +(T_((i,j-1))+T_(old(i,j+1)))/dy^2 ]`
3. Successive Over Relaxation method: The following is the equation used in an iterative manner until the error is acceptable.
`T_(i,j)=T_(old(i,j))*(1-ω)+(ω*T_(GS))`
Where TGS refers to the temperature obtained from the Gauss-Seidel equation and `omega` is the over-relaxation factor.
The following is the MATLAB code to solve for the 2D steady-state heat conduction in the implicit method using iterative solvers.
%MATLAB code to solve the 2D steady state heat conduction equation using iterative solvers.
clear all
close all
clc
%defining the boundary
x = 1:10;
dx = abs(x(1)-x(2));
nx = length(x);
y = 1:10;
dy = abs(y(1)-y(2));
ny = length(y);
%defining the boundary conditions
t_top = 600;
t_bottom = 900;
t_left = 400;
t_right = 800;
%defining the initial conditions and the 1st guess value of temperature.
t = 300*ones(nx,ny);
t(1,:) = t_top;
t(ny,:) = t_bottom;
t(:,1) = t_left;
t(:,nx) = t_right;
t(1,1) = (t_left + t_bottom)/2;
t(1,ny) = (t_top + t_left)/2;
t(nx,ny) = (t_top + t_right)/2;
t(nx,1) = (t_right + t_bottom)/2;
%initialising the t_initial and t_old for future computation.
t_initial = t;
t_old = t;
%defining the k term for future simplification.
k = (dx^2*dy^2)/(2*(dx^2+dy^2));
%defining the tolerance and error values
tolerance = 1e-4;
error = 1;
%starting the iteration count.
iteration = 0;
%starting the while loop that uses Jacobi method and it runs until the error value is acceptable.
while error>tolerance
%starting the spatial loops
for i = 2:(length(x)-1)
for j = 2:(length(y)-1)
%calculating the temperature of a point using Jacobi method
term1 = (t_old(i-1,j)+t_old(i+1,j))/(dx^2);
term2 = (t_old(i,j-1)+t_old(i,j+1))/(dy^2);
t(i,j) = k*(term1 + term2);
end
end
%calculating the error value
error = max(max(abs(t_old - t)));
%updating the t_old
t_old = t;
%updating the iteration count.
iteration = iteration + 1;
end
%ploting for the solution obtained using the Jacobi method.
figure(1)
[c,h] = contourf(x,y,t);
clabel(c,h);
colorbar
colormap(jet)
set(gca,'YDIR','reverse');
xlabel('X Axis');
ylabel('Y Axis');
title({['Plot for Jacobi method'],['Number of iterations = ',num2str(iteration)]});
%reinitialising the temperature, error and iterations count and repeat the above steps for Gauss-Seidel method.
t = t_initial;
t_old = t;
error = 1;
iteration = 0;
while error>tolerance
for i = 2:(length(x)-1)
for j = 2:(length(y)-1)
term1 = (t(i-1,j)+t_old(i+1,j))/(dx^2);
term2 = (t(i,j-1)+t_old(i,j+1))/(dy^2);
t(i,j) = k*(term1 + term2);
end
end
error = max(max(abs(t_old - t)));
t_old = t;
iteration = iteration + 1;
end
figure(2)
[c,h] = contourf(x,y,t);
clabel(c,h);
colorbar
colormap(jet)
set(gca,'YDIR','reverse');
xlabel('X Axis');
ylabel('Y Axis');
title({['Plot for Gauss Seidel method'],['Number of iterations = ',num2str(iteration)]});
%reinitialize the temperature, error and the iteration count and solve for the temperature using SOR method.
t = t_initial;
t_old = t;
error = 1;
iteration = 0;
while error>tolerance
for i = 2:(length(x)-1)
for j = 2:(length(y)-1)
term1 = (t(i-1,j)+t_old(i+1,j))/(dx^2);
term2 = (t(i,j-1)+t_old(i,j+1))/(dy^2);
t(i,j) = k*(term1 + term2);
end
end
error = max(max(abs(t_old - t)));
t_old = t;
iteration = iteration + 1;
end
figure(2)
[c,h] = contourf(x,y,t);
clabel(c,h);
colorbar
colormap(jet)
set(gca,'YDIR','reverse');
xlabel('X Axis');
ylabel('Y Axis');
title({['Plot for Gauss Seidel method'],['Number of iterations = ',num2str(iteration)]});
t = t_initial;
t_old = t;
error = 1;
iteration = 0;
omega = 1.2;
while error>tolerance
for i = 2:(length(x)-1)
for j = 2:(length(y)-1)
term1 = (t(i-1,j)+t_old(i+1,j))/(dx^2);
term2 = (t(i,j-1)+t_old(i,j+1))/(dy^2);
t(i,j) = k*(term1 + term2);
t(i,j) = t_old(i,j)*(1-omega)+(omega*t(i,j));
end
end
error = max(max(abs(t_old - t)));
t_old = t;
iteration = iteration + 1;
end
figure(3)
[c,h] = contourf(x,y,t);
clabel(c,h);
colorbar
colormap(jet)
ax = gca;
ax.YDir = 'reverse';
xlabel('X Axis');
ylabel('Y Axis');
title({['Plot for Successive Over Relaxation method'],['Number of iterations = ',num2str(iteration)]});The following are the results obtained from the 3 different iterative solvers.
2D Transient state heat conduction:
The 2D heat conduction equation in the transient state is as follows:
`(∂T)/(∂t) = α((del^2T)/(delx^2))+((del^2T)/(dely^2))`
Where, `alpha` is the thermal diffusivity.
Explicit method: We use a time marching approach to solve for the temperature at each point and the equation to solve for the temperature at a particular point is as follows:
`(T^(n+1))_(i,j)=(T^n)_(i,j)+αdt[((T_(i-1,j)-2T_(i,j)+(T_(i+1,j)))/dx^2)+((T_(i,j-1)-2T_(i,j)+(T_(i,j+1)))/dy^2)]`
Where the super-script n+1 refers to the next time step and the super-script n refers to the current time step.
The MATLAB code to solve for the transient state heat conduction in the explicit method is as follows:
%MATLAB code to solve for the transient state conduction in explicit method.
clear all
close all
clc
%defining the boundary.
x = 1:10;
dx = abs(x(1)-x(2));
nx = length(x);
y = 1:10;
dy = abs(y(1)-y(2));
ny = length(y);
%defining the time steps
T = 1000;
dt = 0.1;
%defining the boundary conditions.
t_top = 600;
t_bottom = 900;
t_left = 400;
t_right = 800;
%defining the initial conditions.
t = 300*ones(nx,ny);
t(1,:) = t_top;
t(ny,:) = t_bottom;
t(:,1) = t_left;
t(:,nx) = t_right;
t(1,1) = (t_left + t_bottom)/2;
t(1,ny) = (t_top + t_left)/2;
t(nx,ny) = (t_top + t_right)/2;
t(nx,1) = (t_right + t_bottom)/2;
%initialising the t_initial and t_old for future computaions
t_initial = t;
t_old = t;
%defining the thermal diffusivity.
alpha = 1.2;
%initialising the k1 and k2 for future simplifications
k1 = alpha*(dt/(dx^2));
k2 = alpha*(dt/(dy^2));
%initialising the iteration count.
iteration = 0;
%starting the time loop. The loop is predetermined to run for 500 time steps.
for k = 1:500
%starting the spatial loops.
for j = 2:(ny-1)
for i = 2:(nx-1)
%calculating the temperature at a perticular point.
term1 = (t_old(i-1,j) - 2*t_old(i,j) + t_old(i+1,j));
term2 = (t_old(i,j-1) - 2*t_old(i,j) + t_old(i,j+1));
t(i,j) = t_old(i,j) + (term1*k1) + (term2*k2);
end
end
%updating the t_old and the iteration count,
t_old = t;
iteration = iteration + 1;
%plotting the solution at each iteration.
[c,h] = contourf(x,y,t);
clabel(c,h);
colorbar
colormap(jet)
set(gca,'YDIR','reverse');
xlabel('X Axis');
ylabel('Y Axis');
title({['2D Heat Conduction in Transient State.'],['Number of iterations = ',num2str(iteration)]});
pause(0.01);
endThe result for the above code is as follows:
The solution ran for 500 iterations and running it any further did not result in an observable change.
Implicit method: In the implicit method, we take the following equation to solve for the temperature at a particular node particular time step.
`(T^(n+1))_(i,j)=(T^n)_(i,j)+αdt[((T_(i-1,j)-2T_(i,j)+(T_(i+1,j)))/dx^2)+((T_(i,j-1)-2T_(i,j)+(T_(i,j+1)))/dy^2)]^(n+1)`
Here, the RHS also has the n+1 terms.
These are solved using the different iterative methods.
The MATLAB code to solve for the transient state 2D heat conduction equation using different iterative methods is as follows:
%MATLAB code to solve for transient state heat conduction in implicit methods.
clear all
close all
clc
%defining the boundary.
x = 1:10;
dx = abs(x(1)-x(2));
nx = length(x);
y = 1:10;
dy = abs(y(1)-y(2));
ny = length(y);
%defining the time steps.
T = 1000;
dt = 10;
%defining the initial conditions.
t_top = 600;
t_bottom = 900;
t_left = 400;
t_right = 800;
t = 300*ones(nx,ny);
t(1,:) = t_top;
t(ny,:) = t_bottom;
t(:,1) = t_left;
t(:,nx) = t_right;
t(1,1) = (t_left + t_bottom)/2;
t(1,ny) = (t_top + t_left)/2;
t(nx,ny) = (t_top + t_right)/2;
t(nx,1) = (t_right + t_bottom)/2;
%defining the t_initail and t_old for future computations.
t_initial = t;
t_old = t;
%defining the thermal diffusivity.
alpha = 1.2;
%initialising the k1 and k2 for future simplifications
k1 = alpha*(dt/(dx^2));
k2 = alpha*(dt/(dy^2));
%defining the tolerance and error terms for future computations.
tolerance = 1e-4;
error = 1;
%initializing the iteration count.
iteration = 0;
%starting the time loop.
for nt = 1:1500
%starting the loop that runs until the error value is within tolerance.
while error > tolerance
%starting the spatial loops.
for j = 2:(ny-1)
for i = 2:(nx-1)
%calculating the temperature value at a point using Jacobi method.
term1 = 1/(1+(2*k1)+(2*k2));
term2 = k1*term1;
term3 = k2*term1;
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_initial(i,j)*term1)+(h*term2)+(v*term3);
end
end
%calculating the error
error = max(max(abs(t_old - t)));
%updating the t_old and iteration count.
t_old = t;
iteration = iteration + 1;
end
%updating the t_initial.
t_initial = t;
end
%plotting the solution for Jacobi method.
figure(1)
[c,h] = contourf(x,y,t);
clabel(c,h);
colorbar
colormap(jet)
set(gca,'YDIR','reverse');
xlabel('X Axis');
ylabel('Y Axis');
title({['2D Heat Conduction in Transient State - Implicit method using Jacobi method.'],['Number of iterations = ',num2str(iteration)]});
%reinitialising the error value and the iteration count and repeat the above steps for Gauss Seidel method.
error = 1;
iteration = 0;
for nt = 1:1500
while error > tolerance
for j = 2:(ny-1)
for i = 2:(nx-1)
term1 = 1/(1+(2*k1)+(2*k2));
term2 = k1*term1;
term3 = k2*term1;
h = (t((i-1),j)+t_old((i+1),j));
v = (t(i,(j-1))+t_old(i,(j+1)));
t(i,j) =(t_initial(i,j)*term1)+(h*term2)+(v*term3);
end
end
error = max(max(abs(t_old - t)));
t_old = t;
iteration = iteration + 1;
end
t_initial = t;
end
figure(2)
[c,h] = contourf(x,y,t);
clabel(c,h);
colorbar
colormap(jet)
set(gca,'YDIR','reverse');
xlabel('X Axis');
ylabel('Y Axis');
title({['2D Heat Conduction in Transient State - Implicit method using Gauss Seidel method.'],['Number of iterations = ',num2str(iteration)]});
%reinitialising the error value and the iteration count and repeat the above steps for SOR method.
error = 1;
iteration = 0;
%defining the thermal diffusivity.
omega = 1.2;
for nt = 1:1500
while error > tolerance
for j = 2:(ny-1)
for i = 2:(nx-1)
term1 = 1/(1+(2*k1)+(2*k2));
term2 = k1*term1;
term3 = k2*term1;
h = (t((i-1),j)+t_old((i+1),j));
v = (t(i,(j-1))+t_old(i,(j+1)));
t(i,j) =(t_initial(i,j)*term1)+(h*term2)+(v*term3);
t(i,j) = t_old(i,j)*(1-omega)+omega*t(i,j);
end
end
error = max(max(abs(t_old - t)));
t_old = t;
iteration = iteration + 1;
end
t_initial = t;
end
figure(3)
[c,h] = contourf(x,y,t);
clabel(c,h);
colorbar
colormap(jet)
set(gca,'YDIR','reverse');
xlabel('X Axis');
ylabel('Y Axis');
title({['2D Heat Conduction in Transient State - Implicit method using SOR method.'],['Number of iterations = ',num2str(iteration)]});The results for the above code is as follows:
Results and Conclusions:Both in steady state and the transient state analysis, we see that the solution converges in with lesser number of iterations for the SOR method, then the Gauss Seidel method and the Jacobi method converges after the most number of iterations.