Solving the steady and unsteady 2D heat conduction problem
OBJECTIVE: To solve steady and unsteady 2D heat conduction equations using different iterative techniques for implicit and explicit methods in MATLAB.
PROBLEM STATEMENT:
1. Assume that the domain is a unit square.
2. Assume nx = ny [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
1. Left:400K
2. Top:600K
3. Right:800K
4. Bottom:900K
Now, using the point iterative techniques solve the 2D heat conduction equation to plot the temperature contour of steady-state conditions using these three iterative methods, namely Jacobi method, Gauss seidel method, Successive over-relaxation (SOR) method. Also, plot the temperature contour for transient state condition describing the effect of implicit and explicit methods of the conduction equations. Use the above three iterative methods while illustrating the implicit method of transient state conditions.
THEORY:
The Heat equation is a partial differential equation that describes how the distribution of heat evolves in a solid medium, as it spontaneously flows from places where it is higher towards places where it is lower. It is a special case of the diffusion equation.
`{partial T}/{partial t}= alpha({partial ^{2}T}/{partial x^{2}}+{partial ^{2}T}/{partial y^{2}}+{partial ^{2}T}/{partial z^{2}})`
where `T`is the temperature and `(x,y,z,t)` is space and time coordinate respectively.
By the second law of thermodynamics, heat will flow from hotter bodies to adjacent colder bodies, in proportion to the difference of temperature and the thermal conductivity of the material between them. When heat flows into (or out of) a material, its temperature increases (respectively, decreases), in proportion to the amount of heat divided by the amount (mass) of material, with a proportionality factor called the specific heat capacity of the material.
Therefore, the equation says that the rate `dot(T)` at which the material at a point will heat up (or cool down) is proportional to how much hotter (or cooler) the surrounding material is. The coefficient `alpha` in the equation takes into account the thermal conductivity, the specific heat, and the density of the material.
GOVERNING EQUATIONS:
The governing equations consisting of 2D heat conduction equations in cartesian coordinates which are required to solve the problem statements are shown below:
For, Steady-state condition:
`{partial ^{2}T}/{partial x^{2}}+{partial ^{2}T}/{partial y^{2}} = 0`
For, Transient state condition:
`{partial T}/{partial t}= alpha({partial ^{2}T}/{partial x^{2}}+{partial ^{2}T}/{partial y^{2}})`
Here, `alpha` in the above equations represents the thermal diffusivity of the medium.
DISCRETIZATION:
Partial differential equations (PDEs) constitute by far the biggest source of sparse matrix problems. The typical way to solve such equations is to discretize them, i.e., to approximate them by equations that involve a finite number of unknowns. The matrix problems that arise from these discretizations are generally large and sparse; i.e., they have very few nonzero entries. There are several different ways to discretize a PDE, but we will use Finite Difference approximations to solve our problem statement.
We will use the central differencing method to discretize our second-order partial differential equation of `T` w.r.t. `x` and `y` because it is elliptical or we can say it takes information from all direction in the mesh and will use Forward differencing method to discretize our time derivative term of `T` as it is hyperbolic or we can say it takes information from left to right.
Let `T_{i,j}` be the temperature of the point of intersection of `i^{th}` row and `j^(th)` column in the generated square mesh. The discretized form our PDEs are shown below,
`{partial ^{2}T}/{partial x^{2}} = (T_{i-1,j} - 2*T_{i,j}+ T_{i+1,j})/(Deltax^2)`
`{partial ^{2}T}/{partial y^{2}} = (T_{i,j+1} - 2*T_{i,j}+ T_{i+1,j-1})/(Deltay^2)`
`{partial T}/{partial t} = (T_{i,j}^{n+1} - T_{i,j}^{n})/(Deltat)`
STEADY-STATE ANALYSIS:
Since there is no change of temperature gradient with respect to time, so our governing equation becomes:
`{partial ^{2}T}/{partial x^{2}}+{partial ^{2}T}/{partial y^{2}} = 0`
The combined discretization of this equation is goes as follow,
`(T_{i-1,j} - 2*T_{i,j}+ T_{i+1,j})/(Deltax^2) + (T_{i,j+1} - 2*T_{i,j}+ T_{i,j-1})/(Deltay^2) = 0`
As we know that, `nx=ny=> Deltax=Deltay` , so the above equation reduces to :
`T_[i,j] = 1/4(T_[i-1,j]+T_[i+1,j]+T_[i,j-1]+T_[i,j+1])`
Thus above equation gives the temperature value at each node as a function of the temperature of its neighboring nodes. These equations can now be converged to the final solution using iterative schemes like Jacobi, gauss seidel, successive over-relaxation(SOR) using the MATLAB code written below:
MATLAB code for Steady state:
% Script for solving 2d Steady state heat conduction equation
close all
clear all
clc
%% Inputs
L=1; %length of domain is a unit square
%no of grid points along X & Y
nx=30;
ny=nx;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy= y(2)-y(1);
error = 9e9;
tol = 1e-4;
alpha=1.8;%For SOR
%% Defining boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
%% calculating average temperature at the matrix corners
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;
%% other factors to be considered
[X,Y] = meshgrid(x,y);
Told = T;
%% Calculation of 2d Steady state heat conduction eqn
iterative_solver=1:3;
%% Jacobi Iteration
for iterative_solver=1
tic;
if iterative_solver==1
jacobi_iter=1;
while(error > tol)
for i=2:nx-1
for j=2:ny-1
T(i,j) = 0.25*((Told(i-1,j)+Told(i+1,j)+Told(i,j-1)+Told(i,j+1)));
end
end
error = max(max(abs(Told-T)));
Told = T;
jacobi_iter=jacobi_iter+1;
end
end
F=toc;
%plot
figure(1)
[M,N]=contourf(X,Y,T);
clabel(M,N)
colormap(jet)
colorbar
title_text=sprintf('Jacobi Steady iteration number =%d ,Total Time F=%fs',jacobi_iter,F);
title(title_text)
xlabel('X-axis')
ylabel('Y-axis')
set(gca,'YDIR', 'reverse')
fprintf('No of Iteration in jacobi=%d\n',jacobi_iter);
end
%% GS Iter
for iterative_solver=2
%update values
L=1;
nx=30;
ny=nx;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy= y(2)-y(1);
error = 9e9;
tol = 1e-4;
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;
[X,Y] = meshgrid(x,y);
Told = T;
tic;
if iterative_solver==2
GS_iter=1;
while(error > tol)
for i=2:nx-1
for j=2:ny-1
T(i,j) = 0.25*(T(i-1,j)+Told(i+1,j)+T(i,j-1)+Told(i,j+1));
end
end
error = max(max(abs(Told-T)));
Told = T;
GS_iter=GS_iter+1;
end
end
F=toc;
% plot
figure(2)
[M,N]=contourf(X,Y,T);
clabel(M,N)
colormap(jet)
colorbar
title_text=sprintf('GS steady iteration number =%d ,Total Time F=%fs',GS_iter,F);
title(title_text)
xlabel('X-axis')
ylabel('Y-axis')
set(gca,'YDIR', 'reverse')
fprintf('No of Iteration in GS=%d\n',GS_iter);
end
%% SOR Iter
for iterative_solver=3
%update values
L=1;
nx=30;
ny=nx;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy= y(2)-y(1);
error = 9e9;
tol = 1e-4;
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;
[X,Y] = meshgrid(x,y);
Told = T;
tic;
if iterative_solver==3
SOR_iter=1;
while(error > tol)
for i=2:nx-1
for j=2:ny-1
T(i,j) = ((1-alpha)*Told(i,j)+(alpha*0.25*(T(i,j-1)+Told(i,j+1)+T(i-1,j)+Told(i+1,j))));
end
end
error = max(max(abs(Told-T)));
Told = T;
SOR_iter=SOR_iter+1;
end
end
F=toc;
% plot
figure(3)
[M,N]=contourf(X,Y,T);
clabel(M,N)
colormap(jet)
colorbar
title_text=sprintf('SOR steady iteration number =%d ,Total Time F=%fs',SOR_iter,F);
title(title_text)
xlabel('X-axis')
ylabel('Y-axis')
set(gca,'YDIR', 'reverse')
fprintf('No of Iteration in SOR=%d\n',SOR_iter);
endRESULTS:
TRANSIENT STATE ANALYSIS:
The generalised 2D temperature Equation is,
`α ((∂ ^2 T)/( ∂ x ^2) +α( ∂^ 2 T)/( ∂ y ^2) )= (∂ T)/( ∂ t)`
We discretize the eqn by Central Differencing & time gradient Eqn with Forward Differencing Method
- The Eqn with respect to X becomes,
`(∂ ^2 T)/( ∂ x^ 2) = (T( i − 1 , j )− 2 T (i , j) + T( i + 1 , j))/ (Δ x^ 2)`
- The Eqn with respect to Y becomes,
`(∂^2T)/( ∂y^2) = (T (i , j-1) − 2 T( i , j) + T (i , j+1 ))/(Δ y ^2)`
- The Eqn with respect to time gradient becomes
`(∂ T)/( ∂ t) = (T^( n + 1)( i , j )− T n ( i , j ))/(Δ t)`
We can define 2 constants namely k1 and k2
`k1=α(Δt)/(Δx^2)`
`k2=α(Δt)/(Δy^2)`
Since we assumed nx=ny then we can take `Δ x ^2=Δy ^2`
then the Eqn becomes as
`T ^(n + 1 )( i , j )= (T^ n ( i , j) + k1 ( T (i − 1 , j) + T (i + 1 , j )) ^(n+1)+ k2 ( T (i , j − 1) + T (i , j + 1 ))^(n + 1))/( 1 + 2 K 4 + 2 K 5)`
MATLAB code for Transient state(Explicit):
% Script for solving 2d Transient state heat conduction equation(Explicit)
close all
clear all
clc
%% Inputs
L=1; %length of domain is a unit square
%no of grid points along X & Y
nx=30;
ny=nx;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy= y(2)-y(1);
error = 9e9;
tol = 1e-4;
dt = 0.01;
alpha=1.88e-5;%thermal diffusivity for steel
%% Defining boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
%% calculating average temperature at the matrix corners
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;
%% other factors to be considered
%we know the values for constants k
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[X,Y] = meshgrid(x,y);
Told = T;
%% Calculation of 2d Transient state heat conduction eqn in Explicit
iterative_solver=1;
tic;
if iterative_solver==1
solver_counter=1;
while(error > tol)
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
error = max(max(abs(Told-T)));
Told = T;
solver_counter=solver_counter+1;
end
end
F=toc;
%% plot
figure(1)
[M,N]=contourf(X,Y,T);
clabel(M,N)
colormap(jet)
colorbar
title_text=sprintf('Transient iteration number thru explicit=%d ,Total Time F=%fs',solver_counter,F);
title(title_text)
xlabel('X-axis')
ylabel('Y-axis')
set(gca,'YDIR', 'reverse')
fprintf('No of Iteration in Transient(explicit)=%d\n',solver_counter);RESULTS:
MATLAB code for Transient state(Implicit):
% Script for solving 2d Transient state heat conduction equation(Implicit)
close all
clear all
clc
%% Inputs
L=1; %length of domain is a unit square
%no of grid points along X & Y
nx=30;
ny=nx;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy= y(2)-y(1);
error = 9e9;
tol = 1e-4;
dt = 2000;
alpha=1.88e-5;%thermal diffusivity for steel
alpha1=1.8;%For SOR
%% Defining boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
%% calculating average temperature at the matrix corners
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;
%% other factors to be considered
%we know the values for constants k
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[X,Y] = meshgrid(x,y);
Told = T;
T_Prev_dt = Told;
%% Calculation of 2d Transient state heat conduction eqn(Implicit)
iterative_solver=1:3;
%% Jacobi Iteration
term1=(1+2*k1+2*k2)^(-1);
term2=k1*term1;
term3=k2*term1;
for iterative_solver=1
tic;
if iterative_solver==1
jacobi_iter=1;
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)));
Told = T;
jacobi_iter=jacobi_iter+1;
end
end
F=toc;
%plot
figure(1)
[M,N]=contourf(X,Y,T);
clabel(M,N)
colormap(jet)
colorbar
title_text=sprintf('Jacobi UnSteady(Implicit) iteration number =%d ,Total Time F=%fs',jacobi_iter,F);
title(title_text)
xlabel('X-axis')
ylabel('Y-axis')
set(gca,'YDIR', 'reverse')
fprintf('No of Iteration in jacobi=%d\n',jacobi_iter);
end
%% GS Iter
%update values
L=1;
nx=30;
ny=nx;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy= y(2)-y(1);
error = 9e9;
tol = 1e-4;
dt = 2000;
alpha=1.88e-5;%thermal diffusivity for steel
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[X,Y] = meshgrid(x,y);
Told = T;
T_Prev_dt = Told;
term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
for iterative_solver=2
tic;
if iterative_solver==2
GS_iter=1;
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)));
Told = T;
GS_iter=GS_iter+1;
end
end
F=toc;
% plot
figure(2)
[M,N]=contourf(X,Y,T);
clabel(M,N)
colormap(jet)
colorbar
title_text=sprintf('GS Transient(Implicit) iteration number =%d ,Total Time F=%fs',GS_iter,F);
title(title_text)
xlabel('X-axis')
ylabel('Y-axis')
set(gca,'YDIR', 'reverse')
fprintf('No of Iteration in GS=%d\n',GS_iter);
end
%% SOR Iter
%update values
L=1;
nx=30;
ny=nx;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy= y(2)-y(1);
error = 9e9;
tol = 1e-4;
dt = 2000;
alpha=1.88e-5;%thermal diffusivity for steel
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[X,Y] = meshgrid(x,y);
Told = T;
T_Prev_dt = Told;
term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
for iterative_solver=3
tic;
if iterative_solver==3
SOR_iter=1;
while(error > tol)
for i=2:nx-1
for j=2:ny-1
H=(alpha1*T(i-1,j)+(1-alpha1)*Told(i-1,j)+Told(i+1,j));
V=(alpha1*T(i,j-1)+(1-alpha1)*Told(i,j-1)+Told(i,j+1));
T(i,j)=T_Prev_dt(i,j)*term1+(term2*H+term3*V);
end
end
error = max(max(abs(Told-T)));
Told = T;
SOR_iter=SOR_iter+1;
end
end
F=toc;
% plot
figure(3)
[M,N]=contourf(X,Y,T);
clabel(M,N)
colormap(jet)
colorbar
title_text=sprintf('SOR Transient(Implicit) iteration number =%d ,Total Time F=%fs',SOR_iter,F);
title(title_text)
xlabel('X-axis')
ylabel('Y-axis')
set(gca,'YDIR', 'reverse')
fprintf('No of Iteration in SOR=%d\n',SOR_iter);
endRESULTS:
CONCLUSIONS:
- From the above plots it is clear that Jacobi Iterations > Gauss Seidel Iterations > Successive Over Relaxation(SOR).
- Successive over-relaxation Iterative method proves to be solving the equation converging towards the solution at the small-time as well as with fewer numbers of iterative cycles with minimum error values.
- So the SUCCESSIVE OVER RELAXATION Method proves to be the Cost & time effective method when compared to Jacobi & Gauss-Seidel iterative Methods.