Solving the steady and unsteady 2D heat conduction problem
Given :-
Boundary conditions :- Left:400K
Top:600K
Right:800K
Bottom:900K
tolerance : 1e-4
nx = ny
Steady State Case :
2D heat conduction equation for steady state:
`(∂^2T)/(∂x^2) + (∂^2T)/(∂y^2) = 0`
`((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)) = 0`
after solving the eqaution we get,
`T(i,j) = (T(i-1,j) + T(i+1,j) + T(i,j-1) + T(i,j+1)) / 4`
Program for steaty state equation using jacobi method :-
clear all
close all
clc
% inputs
L = 1;
nx = 20;
ny = nx;
omega = 1.3;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
% BC
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_B + T_R)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_l + T_B)/2;
[X,Y] = meshgrid(x,y);
told = T;
error = 9e9;
tol = 1e-4;
for iterative_solver = 1
% jacobi method
if iterative_solver ==1
tic;
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
p = toc;
figure(1)
[M,N] = contourf(X,Y,T);
clabel(M,N);
colormap(jet);
colorbar
title_text = sprintf('No. of iterations(jacobi) = %d \t simulation time = %fs',jacobi_iter,p);
xlabel('x axis');
ylabel('y axis');
title(title_text);
set(gca,'ydir','reverse')
end
end
Program for steaty state equation using Gauss Seidal method :-
clear all
close all
clc
% inputs
L = 1;
nx = 20;
ny = nx;
omega = 1.3;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
% BC
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_B + T_R)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_l + T_B)/2;
[X,Y] = meshgrid(x,y);
told = T;
error = 9e9;
tol = 1e-4;
for iterative_solver = 1
% Gauss Seidal method
if iterative_solver ==1
tic;
Gauss_Seidal_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;
Gauss_Seidal_iter = Gauss_Seidal_iter + 1;
end
p = toc;
figure(1)
[M,N] = contourf(X,Y,T);
clabel(M,N);
colormap(jet);
colorbar
title_text = sprintf('No. of iterations = %d \t simulation time = %fs', Gauss_Seidal_iter,p);
xlabel('x axis');
ylabel('y axis');
title(title_text);
set(gca,'ydir','reverse')
end
endProgram for steaty state equation using Successive over relaxation method :-
clear all
close all
clc
% Inputs
L = 1;
nx = 20;
ny = nx;
omega = 1.3;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
% BC
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_B + T_R)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_l + T_B)/2;
[X,Y] = meshgrid(x,y);
told = T;
error = 9e9;
tol = 1e-4;
for iterative_solver = 1
% Successive over relaxation method
if iterative_solver ==1
tic;
successive_over_relaxation = 1;
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
T(i,j) = (1-omega)*told(i,j)+(omega*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;
successive_over_relaxation = successive_over_relaxation + 1;
end
p = toc;
figure(1)
[M,N] = contourf(X,Y,T);
clabel(M,N);
colormap(jet);
colorbar
title_text = sprintf('No. of iterations(SOR) = %d \t simulation time = %fs', successive_over_relaxation,p);
xlabel('x axis');
ylabel('y axis');
title(title_text);
set(gca,'ydir','reverse')
end
endOutput for steady state equation using implicit methods:-
- Output for jacobi method
- Output for Gauss Seidal method
- Output for Successive over relaxation
- In the above heat conduction equation for steady state case output, we can see that the no of iterations required vary for every method used.
- Jacobi method requires high no. of iterations beacause this method uses same old values for entire calculations till it gets stable.
- Gauss Seidal method requires less no. of iterations for same given inputs because it uses updated values of T i.e. after calculating the value of T at certain node then it uses that updated value of T for next node.
- Successive over relaxation requires relatively very less no. of iterations to get stable, as it uses same equation as gauss seidal but it also uses Relaxation Factor ω which makes it to calculate faster. Also the relaxation factor affects the speed of calculation i.e if we use higher ω which will cause it to overshoot and make SOR slower causing it high no. of iterations.
Transient Case :
2D heat conduction equation for transient case :
`(∂T)/(∂t) = α ((∂^2T)/(∂x^2) + (∂^2T)/(∂y^2))`
Explicit method (for transient case) :-
solving the above equation explicitly :-
Using dicretization method & using Time Marching solution we get,
`( 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`
`T^(n+1)(i,j) = T^n(i,j) + k1(T(i+1,j) - 2T(i,j) + T(i-1,j))^n + k2(T(i,j+1) - 2T(i,j) + T(i,j-1)) ^n`
where `k1 = (∝*dt)/dx^2` & `k2 = (∝*dt)/dy^2`
Program for Transient equation using Explicit method :-
clear all
close all
clc
% Inputs
L = 1;
nx = 20;
ny = nx;
omega = 1.3;
alpha = 1.4;
dt = 1/3000;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
% BC
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_B + T_R)/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);
nt = 800;
CFL = alpha*((1/dx^2) + (1/dy^2))*(dt)
[X,Y] = meshgrid(x,y);
told = T;
tol = 1e-4;
tic;
% Time Loop
for g = 1:nt
% Nodal Loop
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
told = T;
end
p = toc;
figure(1)
[M,N] = contourf(X,Y,T);
clabel(M,N);
colormap(jet);
colorbar
title_text = sprintf('No. of iterations = %d \t time taken = %fs',nt,p);
xlabel('x axis');
ylabel('y axis');
title(title_text);
set(gca,'ydir','reverse')
- Output for Transient equation using explicit method:-
- Here as we using explicit method, so we need to make sure that CFL (i.e. alpha*dt/dx^2) should be < or = 0.5 to make the solution stable. Here CFL = k1 + k2 because we are dealing with 2d heat conduction equation, as k1 = alpha*dt/dx^2 & k2 = alpha*dt/dy^2.
- after assuming nt as 800 and dt as 1/3000 we get the output stable as CFL is < 0.5.
- In explicit method, dt should be very less to maintain stability which results in higher computing time to make calculations than usual.
- Also the stability depends upon numerical errors i.e. if it amplifies from one martching step to next then the calculations will become unstable and if the numerical error decreases after every next marching step then the calculations will be stable.
Implicit method (for Transient case) :-
2D heat conduction equation for transient case :
`(∂T)/(∂t) = α ((∂^2T)/(∂x^2) + (∂^2T)/(∂y^2))`
solving this equation using dicretization we get,
`T^(n+1)(i,j) = (T^n(i,j)) / (1+2k1+2k2) + (k1( T(i+1,j) + T(i-1,j))^(n+1)) / (1+2k1+2k2) + (k2( T(i,j+1) + T(i,j-1))^(n+1)) / (1+2k1+2k2)`
where `k1 = (∝*dt)/dx^2` & `k2 = (∝*dt)/dy^2`
Program for Transient equation using Jacobi method :-
clear all
close all
clc
% Inputs
L = 1;
nx = 20;
ny = nx;
omega = 1.3;
alpha = 1.4;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
% BC
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_B + T_R)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_l + T_B)/2;
dt = 20;
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
nt = 1000;
CFL = alpha*((1/dx^2) + (1/dy^2))*(dt);
[X,Y] = meshgrid(x,y);
told = T;
t_prev = told;
tol = 1e-4;
error = 9e9;
for iterative_solver = 1
% Jacobi method
if iterative_solver ==1
tic;
jacobi_iter = 1;
% Time Loop
for g = 1:nt
while (error > tol)
% Nodal Loop
for i = 2:nx-1
for j = 2:ny-1
term1 = (1+(2*k1)+(2*k2))^(-1);
term2 = k1*term1;
term3 = k2*term1;
H = (told(i-1,j) + told(i+1,j));
V = (told(i,j-1) + told(i,j+1));
T(i,j) = ((t_prev(i,j)*term1) + (H*term2) + (V*term3));
end
end
error = max(max(abs(told - T)));
told = T;
jacobi_iter = jacobi_iter + 1;
end
t_prev = T;
end
p = toc;
figure(1)
[M,N] = contourf(X,Y,T);
clabel(M,N);
colormap(jet);
colorbar
title_text = sprintf('No. of iterations(Jacobi) = %d simulation time = %fs',jacobi_iter,p);
xlabel('x axis');
ylabel('y axis');
title(title_text);
set(gca,'ydir','reverse')
end
end
Program for Transient equation using Gauss Seidal method :-
clear all
close all
clc
% Inputs
L = 1;
nx = 20;
ny = nx;
omega = 1.3;
alpha = 1.4;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
% BC
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_B + T_R)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_l + T_B)/2;
dt = 20;
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
nt = 1000;
CFL = alpha*((1/dx^2) + (1/dy^2))*(dt);
[X,Y] = meshgrid(x,y);
told = T;
t_prev = told;
tol = 1e-4;
error = 9e9;
for iterative_solver = 1
% Gauss Seidal method
if iterative_solver ==1
tic;
Gauss_seidal_iter = 1;
% Time Loop
for g = 1:nt
while (error > tol)
% Nodal Loop
for i = 2:nx-1
for j = 2:ny-1
term1 = (1+(2*k1)+(2*k2))^(-1);
term2 = k1*term1;
term3 = k2*term1;
H = (T(i-1,j) + told(i+1,j));
V = (T(i,j-1) + told(i,j+1));
T(i,j) = ((t_prev(i,j)*term1) + (H*term2) + (V*term3));
end
end
error = max(max(abs(told - T)));
told = T;
Gauss_seidal_iter = Gauss_seidal_iter + 1;
end
t_prev = T;
end
p = toc;
figure(1)
[M,N] = contourf(X,Y,T);
clabel(M,N);
colormap(jet);
colorbar
title_text = sprintf('No. of iterations(GS) = %d simulation time = %fs',Gauss_seidal_iter,p);
xlabel('x axis');
ylabel('y axis');
title(title_text);
set(gca,'ydir','reverse')
end
endProgram for Transient equation using Successive Over Relaxation method :-
clear all
close all
clc
% Inputs
L = 1;
nx = 20;
ny = nx;
omega = 1.3;
alpha = 1.4;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
% BC
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_B + T_R)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_l + T_B)/2;
dt = 20;
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
nt = 1000;
CFL = alpha*((1/dx^2) + (1/dy^2))*(dt);
[X,Y] = meshgrid(x,y);
told = T;
t_prev = told;
tol = 1e-4;
error = 9e9;
for iterative_solver = 1
% Successive Over Relaxation method
if iterative_solver ==1
tic;
Successive_over_relaxation_iter = 1;
% Time Loop
for g = 1:nt
while (error > tol)
% Nodal Loop
for i = 2:nx-1
for j = 2:ny-1
term1 = (1+(2*k1)+(2*k2))^(-1);
term2 = k1*term1;
term3 = k2*term1;
H = (T(i-1,j) + told(i+1,j));
V = (T(i,j-1) + told(i,j+1));
T(i,j) = ((1-omega)*told(i,j))+omega*((t_prev(i,j)*term1) + (H*term2) + (V*term3));
end
end
error = max(max(abs(told - T)));
told = T;
Successive_over_relaxation_iter = Successive_over_relaxation_iter + 1;
end
t_prev = T;
end
p = toc;
figure(1)
[M,N] = contourf(X,Y,T);
clabel(M,N);
colormap(jet);
colorbar
title_text = sprintf('No. of iterations(SOR) = %d simulation time = %fs',Successive_over_relaxation_iter,p);
xlabel('x axis');
ylabel('y axis');
title(title_text);
set(gca,'ydir','reverse')
end
end
Output for transient equation using implicit methods:-
- output for Jacobi method
- output for Gauss Seidal method
- output for Successive over relaxation method
- In the above heat conduction equation for transient case output, we can see that the no of iterations required vary for every method used.
- Jacobi method requires high no. of iterations beacause this method uses same old values for entire calculations till it gets stable.
- Gauss Seidal method requires less no. of iterations for same given inputs because it uses updated values of T i.e. after calculating the value of T at certain node then it uses that updated value of T for next node.
- Successive over relaxation requires relatively very less no. of iterations to get stable, as it uses same equation as gauss seidal but it also uses Relaxation Factor ω which makes it to calculate faster. Also the relaxation factor affects the speed of calculation i.e if we use higher ω which will cause it to overshoot and make SOR slower causing it high no. of iterations.
- Here we can also see that the dt values are higher than explicit method which states that, in implicit method the dt value don't matter much. That means stability can be maintained at higher values of dt using fewer time steps and calculated in vey less computational time, as shown in figure below.
- For the output shown below, where dt = 40 and time steps = 900, this output still shows stability with a minor increment in time taken for the calculation than the one with dt = 20 which is shown above using same iterative method (i.e. Jacobi ).
- In implicit method, as the dt increases so does the truncation error, which may show less accuracy.
Conclusion :- based on the outputs,
- steady state equation requires less time for computaion than transient case.
- Also the SOR method gives solution in less iterations and less time for implicit method.
- whereas the computational time for explicit and implicit depends upon dt where explicit needs smaller dt values and implicit can compute at higher dt values.
- Explicit method depends upon CFL conditions for stable solution (i.e. CFL < 0.5).