Week 5 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Objective
To solve the steady and unsteady (implicit, explicit methods) 2-D heat conduction equation by using finite difference discretization, and comment on accuracy, stability and time step.
Introduction
The governing partial differential equation for steady case
`{d^2T}/dx^2 +{d^2T}/dy^2 = 0`
The governing partial differential equation for unsteady case
`{dT}/dt = alpha*({d^2T}/dx^2 + {d^2T}/dy^2)`
The discretized equation can was
for unsteady case
`T^{n+1}(i,j) = T^n(i,j) + (alpha*dt)/(dx)^2[T(i+1,j)-2*T(i,j)+T(i-1,j)]+(alpha*dt)/(dy)^2[T(i,j+1)-2*T(i,j)+T(i,j-1)]`
for steady case
`{T(i+1,j)-2*T(i,j)+T(i-1,j)}/(dx)^2 +{T(i,j+1)-2*T(i,j)+T(i,j-1)}/(dx)^2`
In steady case for getting a required convergence criterion, following iterative methods are used
- Jacobi method.
- Gauss-Seidel method.
- Successive over-relaxation
Basically, in Jacobi method present iteration value computed by using old iteration values, initially we are guessing the values, then by using those values, we iterate further until achieving required convergent values.
Basically, in Gauss-Seidel method present iteration values computed by using old iteration values and also already computed present iteration values, remaining procedure same as Jacobi method.
Basically, SOR method derived from extrapolating the Gauss-Seidel method, it will improve the rate of convergence, i.e. by providing some relaxation factors.
STEADY STATE CODE;
% 2-D heat conduction problem
% Inputs
L = 1; % LENGTH
W = 1; % WIDTH
nx = 10; % NUMBER OF DIVI....
ny = 10;
dx = L/nx;
dy = W/ny;
alpha = 1.4;
w = 1.6 % relaxation
k = (dx^2+dy^2)/(dx^2)*(dy^2)
x = linspace(0,1,nx+1);
y = linspace(0,1,ny+1);
X = repmat(x,ny+1,1);
Y = repmat(y,nx+1,1);
%plot(X,Y,'k');
%grid on
% initial and boundery condition
T = ones(nx+1,ny+1);
T(1,2:ny) = 600; % top side
T(nx+1,2:ny) = 900; % bottom side
T(2:ny,1) = 400; % left side
T(2:ny,nx+1) = 800; % right side
T(1,1) = (T(1,2)+T(2,1))*0.5;
T(ny+1,1)= (T(ny+1,2)+T(ny,1))*0.5;
T(ny+1,nx+1)= (T(ny,nx+1)+T(ny+1,nx))*0.5;
T(1,nx+1) = (T(1,nx)+T(2,nx+1))*0.5;
Told = T ;
T_initial = T;
% discritization
error = 16;
for m = 1:3
Told = T_initial; % initialising
T_prev = T_initial;
error = 16; % let inirial error , without this it will not entered into convegence loop
itter = 1;
while (error > 1e-04) % convergence loop
for i = 2:nx
for j = 2:ny
if m == 1
T(i,j) = 0.25*(Told(i+1,j)+Told(i-1,j)+Told(i,j+1)+Told(i,j-1)); % jacobi
elseif m == 2
T(i,j) = 0.25*(Told(i+1,j)+T(i-1,j)+Told(i,j+1)+T(i,j-1)); % GAUSS-seidal
elseif m == 3
T(i,j) = (1-w)*T_prev(i,j) +w*(0.25*(Told(i+1,j)+T(i-1,j)+Told(i,j+1)+T(i,j-1))); % SOR
end
end
end
error = max(max(abs(T-Told))); % error value
Told = T; % update the Told value
T_prev = T;
itter = itter + 1;
end
if m == 1 % REQUIRED PLOTS
figure(1);
contourf(x,y,T)
set(gca,'YDIR','reverse')
colorbar;
colormap(jet);
xlabel('X-axis');
ylabel('y-axis');
set(gca,'Xtick',(0:dx:1))
set(gca,'Ytick',(0:dy:1))
Title_txt = sprintf('number of itter for jacobi = %d', itter);
title(Title_txt);
grid on
elseif m == 2
figure(2);
contourf(x,y,T)
set(gca,'YDIR','reverse')
set(gca,'Xtick',(0:dx:1))
set(gca,'Ytick',(0:dy:1))
colorbar;
colormap(jet);
xlabel('X-axis');
ylabel('y-axis');
Title_txt = sprintf('number of itter for gause seidal = %d', itter);
title(Title_txt);
grid on
elseif m == 3
figure(3);
contourf(x,y,T)
set(gca,'Xtick',(0:dx:1))
set(gca,'Ytick',(0:dy:1))
set(gca,'YDIR','reverse')
colorbar;
colormap(jet);
xlabel('X-axis');
ylabel('y-axis');
Title_txt = sprintf('number of itter sor = %d', itter);
title(Title_txt);
grid on
end
end
The following result are obtained
the above graph shows that to get the required convergence, the number of iterations run for jacobi method more as compared to seidal and SOR.
The order of rate of convergence
JACOBI < GAUSS-SEIDAL < SOR
(2) FOR UNSTEADY EXPLICIT METHODE CODE
The time domain can be discretized by using Euler scheme, in this problem we have to consider the time step such a way to satisfies the required courant number for conditional stability criteria.
For diffusive problems the courant number should be less than 0.5, then only the solution will be stable.
clear all
close all
clc
% inputs
L = 1 ;
W = 1;
nx = 10;
ny = 10;
dx = L/nx;
dy = W/ny;
w = 1.59; % relaxation
% meshing
x = linspace(0,1,nx+1);
y = linspace(0,1,ny+1);
X = repmat(x,ny+1,1);
Y = repmat(y,nx+1,1);
alpha = 1.4;
dt = 0.001;
run_time = 0.1;
K_1 = (alpha*dt)/(dx^2);
K_2 = (alpha*dt)/(dy^2);
T = ones(nx+1,ny+1);
T(nx+1,2:ny) = 900; % bottom side
T(1,2:ny) = 600;
T(2:ny,1) = 400; % left side
T(2:ny,nx+1) = 800; % right side
T(1,1) = (T(1,2)+T(2,1))*0.5;
T(ny+1,1)= (T(ny+1,2)+T(ny,1))*0.5;
T(ny+1,nx+1)= (T(ny,nx+1)+T(ny+1,nx))*0.5;
T(1,nx+1) = (T(1,nx)+T(2,nx+1))*0.5;
T_initial = T;
T_rev = T;
Told = T;
% discritization
for k = 1:3 % loop for required methode
T_rev = T_initial; % initialising
Told = T_initial;
time = 1;
for t = 1:500 % run for required number of time steps
T(1,2:ny) = 600; % top side
T(nx+1,2:ny) = 900; % bottom side
T(2:ny,1) = 400; % left side
T(2:ny,nx+1) = 800; % right side
T(1,1) = (T(1,2)+T(2,1))*0.5;
T(ny+1,1)= (T(ny+1,2)+T(ny,1))*0.5;
T(ny+1,nx+1)= (T(ny,nx+1)+T(ny+1,nx))*0.5;
T(1,nx+1) = (T(1,nx)+T(2,nx+1))*0.5;
for i = 2:nx % space loop
for j = 2:ny
if k == 1
T(i,j) = T_rev(i,j) + (K_1)*(Told(i-1,j)-2*Told(i,j)+Told(i+1,j))+(K_2)*(Told(i,j-1)-2*Told(i,j)+Told(i,j+1)) ;
elseif k == 2
T(i,j) = T_rev(i,j) + (K_1)*(T(i-1,j)-2*Told(i,j)+Told(i+1,j))+(K_2)*(T(i,j-1)-2*Told(i,j)+Told(i,j+1)) ;
elseif k == 3
T(i,j) = (1-w)*Told(i,j) + w*(T_rev(i,j) + (K_1)*(T(i-1,j)-2*Told(i,j)+Told(i+1,j))+(K_2)*(T(i,j-1)-2*Told(i,j)+Told(i,j+1)));
end
end
end
error = abs(max(max(T-Told))); % error value
T_rev = T; % updating the values
Told = T;
time = time + 1;
end
if k == 1 %for required plots after the time loop
figure(1);
contourf(x,y,T);
set(gca,'YDIR','reverse')
colorbar;
colormap(jet);
xlabel('X-axis');
ylabel('y-axis');
set(gca,'Xtick',(0:dx:1))
set(gca,'Ytick',(0:dy:1))
title_text = sprintf('number of time steps run for jacobi = %d, error = %d',time,error);
title(title_text);
elseif k == 2
figure(2);
contourf(x,y,T);
set(gca,'YDIR','reverse')
colorbar;
colormap(jet);
xlabel('X-axis');
ylabel('y-axis');
set(gca,'Xtick',(0:dx:1))
set(gca,'Ytick',(0:dy:1))
title_text = sprintf('number of time steps run for seidal = %d,error = %d',time,error);
title(title_text);
elseif k == 3
figure(3);
contourf(x,y,T);
set(gca,'YDIR','reverse')
colorbar;
colormap(jet);
xlabel('X-axis');
ylabel('y-axis');
set(gca,'Xtick',(0:dx:1))
set(gca,'Ytick',(0:dy:1))
title_text = sprintf('number of time steps run for SOR = %d, error = %d',time,error);
title(title_text);
grid on
end
end
The following results are obtained
Above problem simulated for time of 0.5 sec, simulated result shows that for Jacobi method has more error value compared to Gauss –Seidel and SOR,
The order of error value
JACOBI > GAUSS-SEIDAL > SOR.
(3) FOR UNSTEADY IMPLICIT
The time domain discretized by using Euler scheme ,
For implicit method their no stability criterion, basically implicit methods are unconditionally stable method’s, basically in this method the right-hand side of discretized equations, are computed at present iteration step, where as in explicit method computed by using previous iteration values. This will make a difference between implicit and explicit methods.
As we are computing right hand side by using present iteration values, but those values not computed at present time step, so we are using the iterative technique to converge the required values in space loop at a particular timestep.
clear all
close all
clc
% input
L = 1;
W = 1;
w = 1.59;
nx = 10;
ny = 10;
dx = L/nx
dy = L/ny
alpha = 1.4;
dt = 0.1
K_1 = alpha*dt/(dx^2);
K_2 = K_1;
x = linspace(0,1,nx+1);
y = linspace(0,1,ny+1);
[X,Y] = meshgrid(x,y);
grid on
term_1 = (1+(2*K_1)+(2*K_2))^-1;
term_2 = K_1*term_1;
term_3 = K_2*term_1;
% boundery conditions
T = ones(nx+1,ny+1); %
T(1,2:ny) = 600; % top side
T(nx+1,2:ny) = 900; % bottom side
T(2:ny,1) = 400; % left side
T(2:ny,nx+1) = 800; % right side
T(1,1) = (T(1,2)+T(2,1))*0.5;
T(ny+1,1)= (T(ny+1,2)+T(ny,1))*0.5;
T(ny+1,nx+1)= (T(ny,nx+1)+T(ny+1,nx))*0.5;
T(1,nx+1) = (T(1,nx)+T(2,nx+1))*0.5;
T_initial = T;
T_old = T;
% descritization
time = 1; % loop for required methodes
for k = 1:3
Told = T_initial ;
T_prev = T_initial;
time =1
for t = 1:100 % time loop
error = 1e09;
iter = 1
while ( error>1e-04 ) % convergence loop
for i = 2:nx
for j = 2:ny
if k == 1
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)*term_1+term_2*H+term_3*V;
elseif k == 2
H = Told(i+1,j)+T(i-1,j);
V = Told(i,j+1)+T(i,j-1);
T(i,j) = T_prev(i,j)*term_1+term_2*H+term_3*V;
elseif k == 3
H = Told(i+1,j)+T(i-1,j);
V = Told(i,j+1)+T(i,j-1);
T(i,j) =(1-w)*Told(i,j) +(w)* (T_prev(i,j)*term_1+term_2*H+term_3*V);
end
end
end
error = max(max(abs(T-Told)));
Told = T; % updating the temperature
iter = iter +1;
end
error_time = max(max(abs(T-T_prev))); % error with respect to present and previous time temp.. value
T_prev = T;
time = time + 1 ;
end
if k == 1 % reqiured plots
figure(1);
contourf(x,y,T);
set(gca,'YDIR','reverse')
set(gca,'Xtick',(0:dx:1))
set(gca,'Ytick',(0:dy:1))
colorbar;
colormap(jet);
xlabel('X-axis');
ylabel('y-axis');
title_text = sprintf(' error_time jacobi = %d',error_time);
title(title_text);
elseif k == 2
figure(2);
contourf(x,y,T);
set(gca,'YDIR','reverse')
set(gca,'Xtick',(0:dx:1))
set(gca,'Ytick',(0:dy:1))
colorbar;
colormap(jet);
xlabel('X-axis');
ylabel('y-axis');
title_text = sprintf('error_time GAUSS-SEIDAL = %d',error_time);
title(title_text);
elseif k == 3
figure(3);
contourf(x,y,T);
set(gca,'YDIR','reverse')
set(gca,'Xtick',(0:dx:1))
set(gca,'Ytick',(0:dy:1))
colorbar;
colormap(jet);
xlabel('X-axis');
ylabel('y-axis');
title_text = sprintf('error_time SOR = %d ',error_time);
title(title_text);
grid on
end
endRESULT
The error value obtained for Jacobi was more than Gauss-Seidel and SOR, and also the solution was run for lesser time steps as compared to explicit method for same convergence values.
CONCLUSIONS
Based on above result we can conclude that
- If we are simulating steady state case, using SOR with proper relaxation factor will give accurate results with lesser number of iterations.
- If we are simulating unsteady case with explicit method, we should compromise between time step size and mesh grid size to obtain satisfied courant number for a conditional stability criterion, and use of SOR with proper relaxation factor will give accurate result with a smaller number of time steps.
- If we are simulating unsteady case with implicit method, we are not bother about stability criterion, because implicit methods are unconditionally stable, and we can use large time step as compared to explicit method.