Solving the Steady State and Transient state of heat conduction equation using matlab
The Heat conduction equation for the steady state ( dt = 0)
`d^2/dx^2 (T) + d^2/dy^2(T) = 0`
Temperature is the function of (x,y).
Solving the steady state heat conduction in the matlab using the taylor function
`(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) = 0`
using the given condition
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
Temperature analysis by steady state
Steady State temperature analysis using jacobi iterative
clear all
close all
clc
% input parameter
nx = 10;
ny = nx;
x = linspace(0,1,nx);
y = linspace(1,0,nx);
dx = x(2) - x(1);
dy = y(2) - y(1);
[X,Y] = meshgrid(x,y);
Temp = 400*ones(nx,ny);
Temp(1,:) = 600;
Temp(10,:) = 900;
Temp((2:9),1) = 400;
Temp((2:9),10) = 800;
Temp(1,1) = (600+400)/2;
Temp(10,1)= (400+900)/2;
Temp(1,10) = (600+800)/2;
Temp(10,10) = (800+900)/2;
Tempold = Temp;
% for steady state analysis
% d^2T/dx^2 + d^2T/dy^2 = 0
% (TL-2Tp+Tr)/dx^2 + (Tt -2Tp +Tb)/dy^2 = 0
error = 9e9;
tolerance = 1e-4;
jacobi_itera = 1;
while(error>tolerance)
for i = 2:nx-1
for j = 2:ny-1
term1 = 0.25*(Tempold(i-1,j)+ Tempold(i+1,j));
term2 = 0.25*(Tempold(i,j-1)+ Tempold(i,j+1));
Temp(i,j) = (term1)+ (term2);
end
end
error = max(max(abs(Tempold - Temp)));
Tempold = Temp;
jacobi_itera = (jacobi_itera + 1);
end
figure(1)
[C,M] = contourf(X,Y,Temp);
clabel(C,M);
colormap(jet);
colorbar;
xlabel(\'xaxis\')
ylabel(\'yaxis\')
title_text1 = sprintf(\'iteration_number1 = %d\',jacobi_itera);
title_text1a = sprintf(\'Steady state = jacobi iteration\');
title({title_text1a,title_text1});
pause(0.003);
2. Steady state temperature analysis using gauss-siedel iteration
clear all
close all
clc
% input parameter
nx = 10;
ny = nx;
x = linspace(0,1,nx);
y = linspace(1,0,nx);
dx = x(2) - x(1);
dy = y(2) - y(1);
[X,Y] = meshgrid(x,y);
Temp = 400*ones(nx,ny);
Temp(1,:) = 600;
Temp(10,:) = 900;
Temp((2:9),1) = 400;
Temp((2:9),10) = 800;
Temp(1,1) = (600+400)/2;
Temp(10,1)= (400+900)/2;
Temp(1,10) = (600+800)/2;
Temp(10,10) = (800+900)/2;
Tempold = Temp;
% for steady state analysis
% d^2T/dx^2 + d^2T/dy^2 = 0
% (TL-2Tp+Tr)/dx^2 + (Tt -2Tp +Tb)/dy^2 = 0
error = 9e9;
tolerance = 1e-4;
gauss_sidel = 1;
while(error>tolerance)
for k = 2:nx-1
for m = 2:ny-1
term1 = 0.25*(Temp(k-1,m)+ Temp(k+1,m));
term2 = 0.25*(Temp(k,m-1) + Temp(k,m+1));
Temp(k,m) = term1 + term2;
end
end
error = max(max(abs(Temp - Tempold)));
Tempold = Temp;
gauss_sidel = gauss_sidel +1;
end
figure(2)
[V,M] = contourf(X,Y,Temp);
clabel(V,M);
colormap(jet);
colorbar;
xlabel(\'xaxis\')
ylabel(\'yaxis\')
title_text2 = sprintf(\'iteration_number = %d\',gauss_sidel);
title_text2a = sprintf(\'Steady state = gauss iteration\');
title({title_text2a,title_text2});
pause(0.003);
3. Steady State temperature analysis using SOR iterative
clear all
close all
clc
% input parameter
nx = 10;
ny = nx;
x = linspace(0,1,nx);
y = linspace(1,0,nx);
dx = x(2) - x(1);
dy = y(2) - y(1);
[X,Y] = meshgrid(x,y);
Temp = 400*ones(nx,ny);
Temp(1,:) = 600;
Temp(10,:) = 900;
Temp((2:9),1) = 400;
Temp((2:9),10) = 800;
Temp(1,1) = (600+400)/2;
Temp(10,1)= (400+900)/2;
Temp(1,10) = (600+800)/2;
Temp(10,10) = (800+900)/2;
Tempold = Temp;
error = 9e9;
tolerance = 1e-4;
% SOR ( successive over relaxation)
alpha = 1.2;
SOR_iterative = 1 ;
while(error>tolerance)
for u = 2:nx-1
for v = 2:ny-1
SOR_term1 = (1-alpha)*Tempold(u,v);
SOR_term2 = alpha*0.25*(Temp(u-1,v) + Tempold(u+1,v));
SOR_term3 = alpha*0.25*(Temp(u,v-1) + Tempold(u,v+1));
Temp(u,v) = SOR_term1 + SOR_term2 + SOR_term3;
end
end
error = max(max(abs(Tempold - Temp)));
Tempold = Temp;
SOR_iterative = SOR_iterative + 1;
end
figure(3)
[M,N] = contourf(X,Y,Temp);
clabel(M,N);
colormap(jet);
colorbar;
xlabel(\'Xaxis\');
ylabel(\'Yaxis\');
title_text3 = sprintf(\'iteration_number3 = %d\',SOR_iterative\');
title_text3a = sprintf(\'Steady state\');
title({title_text3a,title_text3});
pause(0.03);
The figure for the for the following condition
1. Steady State Temperture analysis using jacobi iteration technique
The number of iteration in jacobi iteration is 203.
2. Steady state temperature analysis using Gauss-Siedel technique
The number of iteration is reduced in the gauss siedel and it reached from 203 in jacobi iteration to 110 by gauss-Siedel technique.
3. Steady state temperature analysis using SOR Technique
In SOR(Successive Over relaxation) technique the iteration number is reduced by a large number compared to the Jacobi and Gauss - Siedel technique.
The iteration number in SOR is 74.
The Transient Heat Conduction equation
`d^2/dx^2 (T) + d^2/dy^2(T) = alpha. d(T)/dt`
alpha = heat diffusion coefficient
Temperature is function of temperature and space.
expand the equation using the taylor table
`(T_(i,j)^ (n+1) - T_(i,j)^n )/(Delta t) = alpha. ((T_(i-1,j) - T_(i,j) + T_(i+1,j))/(Deltax^2) + (T_(i,j-1) - T_(i,j) + T_(i,j+1))/(Deltay^2))^(n or (n+1))`
n = for explicit scheme
n+1 = for implicit scheme
using the same the given condition
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
Transient heat analysis using explicit scheme
Transient State of heat conduction analysis using explicit scheme
clear all
close all
clc
% input parameter
nx = 10;
ny = nx;
x = linspace(0,1,nx);
y = linspace(1,0,nx);
dx = x(2) - x(1);
dy = y(2) - y(1);
t = 0.5;
dt = 0.001;
nt = t/dt;
[X,Y] = meshgrid(x,y);
Temp = 100*ones(nx,ny);
Temp(1,:) = 600;
Temp(10,:) = 900;
Temp((2:9),1) = 400;
Temp((2:9),10) = 800;
Temp(1,1) = (600+400)/2;
Temp(10,1)= (400+900)/2;
Temp(1,10) = (600+800)/2;
Temp(10,10) = (800+900)/2;
Tempold = Temp;
temp_coeff = 1.4;
k1 = (temp_coeff * dt)/dx^2;
k2 = (temp_coeff*dt)/dy^2;
error = 6e6;
tolerance = 1e-4;
iterative_solver = 1
while(error>tolerance);
for k = 1:nt
for i = 2:nx-1
for j = 2:ny-1
term1 = Tempold(i,j);
term2 = k1*(Tempold(i-1,j) - 2*Tempold(i,j)+Tempold(i+1,j));
term3 = k2*(Tempold(i,j-1) - 2*Tempold(i,j)+ Tempold(i,j+1));
Temp(i,j) = term1 + term2 + term3;
end
end
error = max(max(abs(Temp - Tempold)));
Tempold = Temp;
iterative_solver = iterative_solver +1 ;
end
end
figure(1)
[C,h] = contourf(X,Y,Temp);
clabel(C,h);
colormap(jet);
colorbar;
xlabel(\'Xaxis\');
ylabel(\'Yaxis\');
title_text = sprintf(\'explicit iter = %d\',iterative_solver);
title(title_text);
pause(0.003)
The output of the above program is shown by the figure
The EXPLICIT method: you use the known y0 (beginning of time interval) as the value of y in this period. Very simple method but not very stable.
Solving the Transient heat analysis using the implicit method
1.Transient State heat analysis by implicit scheme using jacobi iteration
clear all
close all
clc
% input parameter
nx = 10;
ny = nx;
x = linspace(0,1,nx);
y = linspace(1,0,nx);
dx = x(2) - x(1);
dy = y(2) - y(1);
dt = 1;
t = 5;
nt = t/dt
[X,Y] = meshgrid(x,y);
Temp = 100*ones(nx,ny);
Temp(1,:) = 600;
Temp(10,:) = 900;
Temp((2:9),1) = 400;
Temp((2:9),10) = 800;
Temp(1,1) = (600+400)/2;
Temp(10,1)= (400+900)/2;
Temp(1,10) = (600+800)/2;
Temp(10,10) = (800+900)/2;
Tempold = Temp;
Tprev = Temp;
error = 6e6;
tolerance = 1e-4;
temp_coeff = 2.4;
k1 = (temp_coeff*dt)/dx^2;
k2 = (temp_coeff*dt)/dy^2;
k = (1+2*k1 +2*k2)^-1 ;
jacobi_iterative = 1;
while(error>tolerance)
for g = 1:nt
for i = 2 :nx-1
for j = 2:ny-1
term1 = k*Tprev(i,j);
term2 = k*(k1*(Tempold(i-1,j)+ Tempold(i+1,j)));
term3 = k*(k2*(Tempold(i,j-1)+ Tempold(i,j+1)));
Temp(i,j) = term1+term2 +term3;
end
end
error = max(max(abs(Tempold - Temp)));
Tempold = Temp;
jacobi_iterative = jacobi_iterative +1 ;
end
Tprev = Temp;
end
figure(1)
[T,V] = contourf(X,Y,Temp);
clabel(T,V);
colormap(jet);
colorbar;
xlabel(\'xaxis\');
ylabel(\'yaxis\');
title_text1 = sprintf(\'iteration number = %d\',jacobi_iterative);
title_text2 = sprintf(\'Transient state = Jacobi iteration\');
title({title_text1,title_text2});
pause(0.003);
2. Transient State heat analysis by implicit scheme using gauss iteration
clear all
close all
clc
% input parameter
nx = 10;
ny = nx;
x = linspace(0,1,nx);
y = linspace(1,0,nx);
dx = x(2) - x(1);
dy = y(2) - y(1);
dt = 1;
t = 5;
nt = t/dt
[X,Y] = meshgrid(x,y);
Temp = 100*ones(nx,ny);
Temp(1,:) = 600;
Temp(10,:) = 900;
Temp((2:9),1) = 400;
Temp((2:9),10) = 800;
Temp(1,1) = (600+400)/2;
Temp(10,1)= (400+900)/2;
Temp(1,10) = (600+800)/2;
Temp(10,10) = (800+900)/2;
Tempold = Temp;
Tprev = Temp;
error = 6e6;
tolerance = 1e-4;
temp_coeff = 2.4;
k1 = (temp_coeff*dt)/dx^2;
k2 = (temp_coeff*dt)/dy^2;
k = (1+2*k1 +2*k2)^-1 ;
gauss_iterative = 1;
while(error>tolerance)
for g = 1:nt
for i = 2 :nx-1
for j = 2:ny-1
term1 = k*Tprev(i,j);
term2 = k*(k1*(Temp(i-1,j)+ Temp(i+1,j)));
term3 = k*(k2*(Temp(i,j-1)+ Temp(i,j+1)));
Temp(i,j) = term1+term2 +term3;
end
end
error = max(max(abs(Tempold - Temp)));
Tempold = Temp;
gauss_iterative = gauss_iterative +1 ;
end
Tprev = Temp;
end
figure(1)
[T,V] = contourf(X,Y,Temp);
clabel(T,V);
colormap(jet);
colorbar;
xlabel(\'xaxis\');
ylabel(\'yaxis\');
title_text1 = sprintf(\'iteration number = %d\',gauss_iterative);
title_text2 = sprintf(\'Transient state = Jacobi iteration\');
title({title_text1,title_text2});
pause(0.003);
3. Transient heat analysis by Implicit scheme using SOR iteration.
clear all
close all
clc
% input parameter
nx = 10;
ny = nx;
x = linspace(0,1,nx);
y = linspace(1,0,nx);
dx = x(2) - x(1);
dy = y(2) - y(1);
dt = 1;
t = 5;
nt = t/dt
[X,Y] = meshgrid(x,y);
Temp = 100*ones(nx,ny);
Temp(1,:) = 600;
Temp(10,:) = 900;
Temp((2:9),1) = 400;
Temp((2:9),10) = 800;
Temp(1,1) = (600+400)/2;
Temp(10,1)= (400+900)/2;
Temp(1,10) = (600+800)/2;
Temp(10,10) = (800+900)/2;
Tempold = Temp;
Tprev = Temp;
error = 6e6;
tolerance = 1e-4;
alpha= 1.1;
temp_coeff = 2.4;
k1 = (temp_coeff*dt)/dx^2;
k2 = (temp_coeff*dt)/dy^2;
k = (1+2*k1 +2*k2)^-1 ;
SOR_iterative = 1;
while(error>tolerance)
for g = 1:nt
for i = 2 :nx-1
for j = 2:ny-1
term1 = k*Tprev(i,j);
term2 = k*(k1*(Temp(i-1,j)+ Tempold(i+1,j)));
term3 = k*(k2*(Temp(i,j-1)+ Tempold(i,j+1)));
term4 = Tempold(i,j);
Temp(i,j) = (1-alpha)*term4 + alpha*(term1+term2 +term3);
end
end
error = max(max(abs(Tempold - Temp)));
Tempold = Temp;
SOR_iterative = SOR_iterative +1 ;
end
Tprev = Temp;
end
figure(1)
[T,V] = contourf(X,Y,Temp);
clabel(T,V);
colormap(jet);
colorbar;
xlabel(\'xaxis\');
ylabel(\'yaxis\');
title_text1 = sprintf(\'iteration number = %d\',SOR_iterative);
title_text2 = sprintf(\'Transient state = SOR iteration\');
title({title_text1,title_text2});
pause(0.003);
The result of the implicit scheme is shown by the figure
An implicit solution contains information obtained from solving simultaneous equations for the full grid for each time step. This is computationally more demanding but allows for larger time steps and better stability.
(i) Jacobi Iteration
(ii) Gauss - Siedel
(iii) SOR iteration