2D Steady and Un-Steady Heat Conduction
2D Conduction
Assuming that after the domain of uniform discretization, the step length of x-axis is `Δx = x_(i+1) – x_i` and that of y-axis is `Δy = y_(j+1) - y_j` obviously, `x_i = iΔx,y_j = jΔy,and i,j =0,1,2,3 ....`
`n(n = 0,1,2,3 ...)` is used to uniformly discretize the time domain and the step length between two time moments `Δt = t_(n+1 )- t_n,t_n= nΔt` where `(i,j)⋅ ⋅ ⋅ (x_i,y_j ),n ⋅ ⋅ ⋅ t_n; t_(i,j)^n⋅ ⋅ ⋅n`
and the temperature at node `(i,j)` in the time moment `n` is `T(x_i,y_j,t_n).`
The explicit difference array of the two-dimensional unsteady-state heat conduction without internal heat source is expressed as follows.
Applying the first heat conduction equation to node `(i,j)` at the time moment of `n`, the equation can be rewritten as
`((∂^2 t)/(∂x^2 )+(∂^2 t)/(∂y^2 ))_(i,j)^n =1/α*((∂T)/(∂t))_(i,j)^n; -----(1`
The partial differential in the two sides of can be approximated by difference quotient. The temperature item in the right of the equal sign can be approximated by first order time forward difference quotient.
`((∂T)/(∂t))_(i,j)^n=((T_(i,j)^(n+1)-T_(i,j)^n ))/(∆t)` -----(2)
The second-order partial differential in the left of the equal sign can be approximated by the central difference quotient.
`((∂^2 T)/(∂x^2 ))_(i,j)^n=((T_(i-1,j)^n-2T_(i,j)^n+T_(i+1,j)^n ))/((∆x)^2 )`---(3)
`((∂^2 T)/(∂y^2 ))_(i,j)^n=((T_(i-1,j)^n-2T_(i,j)^n+T_(i+1,j)^n ))/((∆y)^2 )`---(4)
Substituting equations (2),(3),(4) into (1) we get
Explicit Form
`T_(i,j)^(n+1)=T_(i,j)^n+α*∆t(((T_(i-1,j)^n-2T_(i,j)^n+T_(i+1,j)^n ))/((∆x)^2 )+((T_(i-1,j)^n-2T_(i,j)^n+T_(i+1,j)^n ))/((∆y)^2 ))`
Similarly, applying the first heat conduction equation to node `(i,j)` at the time moment of `n+1` and approximating the second-order partial differential in the left of the equal sign by the central difference quotient we get,
Implicit Form
`T_(i,j)^(n+1)=T_(i,j)^n+α*∆t(((T_(i-1,j)^(n+1)-2T_(i,j)^n+T_(i+1,j)^(n+1) ))/((∆x)^2 )+((T_(i-1,j)^(n+1)-2T_(i,j)^n+T_(i+1,j)^(n+1) ))/((∆y)^2 ))`
For two dimensional steady state, with no heat generation, the Laplace equation can be applied and approximating the second-order partial differential in the left of the equal sign by the central difference quotient we get,
Steady State
`T_(i,j)^n=(((∆x)^2*(∆y)^2)/(2*(∆x)^2+(∆y)^2))(((T_(i-1,j)^n+T_(i+1,j)^n ))/((∆x)^2 )+((T_(i,j-1)^n+T_(i,j+1)^n ))/((∆y)^2 ))`
Main Code
clear all
close all
clc
nx=20;
ny=nx;
x=linspace(0,5,nx);
y=linspace(0,5,ny);
dx=x(2)-x(1);
dy=dx;
dt=0.005;
nt=80;
alpha=1;
T=300*ones(20,20);
T(:,1)=400;
T(:,end)=800;
T(1,:)=900;
T(end,:)=600;
Told=T;
Tprev=T;
T_initial=T;
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
k3=(1+2*k1+2*k2);
k4=(1-2*k1-2*k2);
w=0.9;
tol=1e-4;
% 2D Transient Heat Conduction (Explicit Form)
for k=1:nt
for j=2:ny-1
for i=2:nx-1
T(i,j)= k4*Told(i,j)+k1*(Told(i-1,j)+Told(i+1,j))+k2*(Told(i,j-1)+Told(i,j+1));
end
end
Told=T;
end
Told=T_initial;
figure(1)
contourf(x,y,T)
[C,h]=contourf(x,y,T);
clabel(C,h)
colormap;
colorbar;
xlabel('X axis')
ylabel('Y axis')
title ('2D Transient Heat Conduction (Explicit Form) at T=0.4s')
T=T_initial;
Tprev=T_initial;
%2D Transient Heat Conduction(Implicit Forms)
[T, jacobian_iter,gauss_iter,SOR_iter] = Transient(T,Told,Tprev,T_initial,k1,k2,k3,nt,nx,ny,tol,x,y,w)
Told=T_initial;
Tprev=T_initial;
T=T_initial;
%2D Steady State Heat Conduction
[T,jacobi_steady_iter,gauss_steady_iter,SOR_steady_iter] =steady(T,Told,T_initial,Tprev,x,y,nx,ny,nt,dx,dy,tol,w)
Un-Steady State
function [T, jacobian_iter,gauss_iter,SOR_iter] = Transient(T,Told,Tprev,T_initial,k1,k2,k3,nt,nx,ny,tol,x,y,w)
jacobian_iter=1;
for k=1:nt
err=1;
while (err>tol)
for j=2:ny-1
for i=2:nx-1
T(i,j)= (Tprev(i,j)+k1*(Told(i-1,j)+Told(i+1,j))+k2*(Told(i,j-1)+Told(i,j+1)))/k3
end
end
err=max(max(abs(Told-T)));
Told=T;
jacobian_iter= jacobian_iter+1;
end
Tprev=T;
end
figure(2)
contourf(x,y,T);
[C,h]=contourf(x,y,T);
clabel(C,h);
colormap;
colorbar;
xlabel('X axis')
ylabel('Y axis')
title ('2D Transient Heat Conduction (Jacobi Implicit) at T=0.4s')
T=T_initial;
Told= T_initial;
Tprev=T_initial;
gauss_iter=1;
for k=1:nt
err=1;
while (err>tol)
for j=2:ny-1
for i=2:nx-1
T(i,j)= (Tprev(i,j)+k1*(T(i-1,j)+T(i+1,j))+k2*(T(i,j-1)+T(i,j+1)))/k3
end
end
err=max(max(abs(T-Told)));
Told=T;
gauss_iter= gauss_iter+1;
end
Tprev=T;
end
figure(3)
contourf(x,y,T)
[C,h]=contourf(x,y,T);
clabel(C,h)
colormap;
colorbar;
xlabel('X axis')
ylabel('Y axis')
title ('2D Transient Heat Conduction (Gauss Siedel Implicit) at T=0.4s')
T=T_initial;
Told= T_initial;
Tprev=T_initial;
Tgauss=T_initial;
SOR_iter=1;
for k=1:nt
err=1;
while (err>tol)
for j=2:ny-1
for i=2:nx-1
Tgauss(i,j)= (Tprev(i,j)+k1*(T(i-1,j)+T(i+1,j))+k2*(T(i,j-1)+T(i,j+1)))/k3
T(i,j) = (1-w)*Tprev(i,j)+ w*Tgauss(i,j)
end
end
err=max(max(abs(T-Told)));
Told=T;
SOR_iter= SOR_iter+1;
end
Tprev=T;
end
figure(4)
contourf(x,y,T)
[C,h]=contourf(x,y,T);
clabel(C,h)
colormap;
colorbar;
xlabel('X axis')
ylabel('Y axis')
title ('2D Transient Heat Conduction (SOR Implicit) at T=0.4s')
T=T_initial;
Told= T_initial;
Tprev=T_initial;
end
Steady State
function [T,jacobi_steady_iter,gauss_steady_iter,SOR_steady_iter] =steady(T,Told,T_initial,Tprev,x,y,nx,ny,nt,dx,dy,tol,w)
k1=0.5*((dy^2)/(dx^2+dy^2));
k2=0.5*((dx^2)/(dx^2+dy^2));
jacobi_steady_iter=1;
err=1;
while(err>tol)
for j=2:ny-1
for i=2:nx-1
T(i,j)= k1*(Told(i-1,j)+Told(i+1,j))+k2*(Told(i,j-1)+Told(i,j+1))
end
end
err= max(max(abs(T-Told)));
Told=T;
jacobi_steady_iter=jacobi_steady_iter+1;
end
figure(5)
contourf(x,y,T)
[C,h]=contourf(x,y,T);
clabel(C,h);
colormap;
colorbar;
xlabel('X axis')
ylabel('Y axis')
title ('2D Steady State Heat Conduction (Jacobi) at T=0.4s')
T=T_initial;
Told= T_initial;
Tprev=T_initial;
gauss_steady_iter=1;
err=1;
while(err>tol)
for j=2:ny-1
for i=2:nx-1
T(i,j)= k1*(T(i-1,j)+T(i+1,j))+k2*(T(i,j-1)+T(i,j+1))
end
end
err= max(max(abs(T-Told)));
Told=T;
gauss_steady_iter=gauss_steady_iter+1;
end
figure(6)
contourf(x,y,T)
[C,h]=contourf(x,y,T);
clabel(C,h);
colormap;
colorbar;
xlabel('X axis')
ylabel('Y axis')
title ('2D Steady State Heat Conduction (Gauss Siedal) at T=0.4s')
T=T_initial;
Told= T_initial;
Tprev=T_initial;
Tgauss=T_initial;
SOR_steady_iter=1;
err=1;
while(err>tol)
for j=2:ny-1
for i=2:nx-1
Tgauss(i,j)= k1*(T(i-1,j)+T(i+1,j))+k2*(T(i,j-1)+T(i,j+1));
T(i,j)= (1-w)*Told(i,j)+w*Tgauss(i,j);
end
end
err= max(max(abs(T-Told)));
Told=T;
SOR_steady_iter=SOR_steady_iter+1;
end
figure(7)
contourf(x,y,T)
[C,h]=contourf(x,y,T);
clabel(C,h);
colormap;
colorbar;
xlabel('X axis')
ylabel('Y axis')
title ('2D Steady State Heat Conduction (SOR) at T=0.4s')
T=T_initial;
Told= T_initial;
Tprev=T_initial;
endPLOTS
Explicit Form
For explicit form,the solution is obtained at t=0.4seconds by time marching using time step value of 0.005 seconds. CFL number (0.1444) is well within the range of `CFL <= 0.25` for nx=ny.Hence , the obtained solution is stable.
Implicit Forms
Number of Iterations =1219 iterations.
Number of Iterations =938 iterations.
Number of Iterations =914 iterations. (`omega`=0.9 i.e it is under relaxed)
Number of Iterations =1166 iterations. (`omega`=1.5 i.e it is over relaxed)
For Implicit forms there is no effect of CFL number on the stability of the solution.
Steady State
Number of Iterations =826 iterations.
Number of Iterations =441 iterations.
Number of Iterations =530 iterations (`omega`=0.9 i.e it is under relaxed)
Number of Iterations =155 iterations (`omega`=1.5 i.e it is over relaxed)
In steady state there is no need of time marching as the temperature at all points remains unchanged.
SOR (over relaxed) is the fastest iterative method for calculating the steady state. Value of `omega` should be optimum so that the solution converges faster.