Solving the steady and unsteady 2D heat conduction problem
Solving steady and unsteady state 2D heat conduction equations
Objective:
Solving steady and unsteady state 2D heat conduction equations by using iterative solvers:
1) Explicit method (for unsteady only)
2) Implicit method (for steady and unsteady)
- Gauss Jacobi
- Gauss Seidel
- Successive Over Relaxation
SOLVING STEADY STATE 2D HEAT CONDUCTION EQUATION:
Boundary conditions for steady and transient case
1. Left:400K
2. Top:600K
3. Right:800K
4. Bottom:900K
we assume nx = ny
For steady state we know that time derivative is zero
then we get the 2D heat equation
here we use central differencing
We take \'omega\' as successive over relaxation parameter, which is defined as
Input values for the problem:
Thermal diffusivity(alpha): 1e-4
Grid value (nx=ny) : 61
Tolerance: 1e-4
MATLAB CODE FOR STEADY STATE (IMPLICIT METHOD):
clear all
close all
clc
% Inputs
nx = 61;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = 1/(nx-1);
dy = 1/(ny-1);
error = 1e9;
tol = 1e-4;
% Initializing the temparature matrix
T = 298*ones(nx,ny);
T(1,:) = 900;
T(end,:) = 600;
T(:,1) = 400;
T(:,end) = 800;
% edge refining
T(1,1) = (900+400)/2;
T(1,end) = (900+800)/2;
T(end,1) = (400+600)/2;
T(end,end) = (800+600)/2;
% Calcution of over relaxation parameter omega
omega = 2/(1+sin(pi*dx));
% solving 2D steady state implicitly
[xx,yy] = meshgrid(x,y);
Told = T;
solver = input(\'jacobi1 siedel2 SOR3 :\')
tic;
% Jacobi iteration
if solver == 1
count_jac = 1;
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = (Told(i-1,j)+Told(i+1,j))/(4);
term2 = (Told(i,j-1)+Told(i,j+1))/(4);
T(i,j) = term1+term2;
end
end
error = max(max(abs(Told-T)));
Told = T;
count_jac = count_jac + 1;
end
end
% Gauss seidel iteration
if solver == 2
count_seidel = 1;
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = (T(i-1,j)+Told(i+1,j))/(4);
term2 = (T(i,j-1)+Told(i,j+1))/(4);
T(i,j) = term1+term2;
end
end
error = max(max(abs(Told-T)));
Told = T;
count_seidel = count_seidel + 1;
end
end
% SOR iteration
if solver == 3
count_sor = 1;
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = (T(i-1,j)+Told(i+1,j))/(4);
term2 = (T(i,j-1)+Told(i,j+1))/(4);
T(i,j) = omega*(term1+term2)-(omega-1)*Told(i,j);
end
end
error = max(max(abs(Told-T)));
Told = T;
count_sor = count_sor + 1;
end
end
time_required = toc;
figure(1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
if solver == 1
title(sprintf(\'Temparature profile of steady state heat conduction [Gauss Jacobi] n No. of Iterations: %d n Computation Time: %0.3f s\',count_jac,time_required));
else if solver == 2
title(sprintf(\'Temparature profile of steady state heat conduction [Gauss Siedel] n No. of Iterations: %d n Computation Time: %0.3f s\',count_seidel,time_required));
else if solver == 3
title(sprintf(\'Temparature profile of steady state heat conduction [SOR] n No. of Iterations: %d n Computation Time: %0.3f s n Relaxation Parameter: %f\',count_sor,time_required,real(omega)));
end
end
end
xlabel(\'xx\');
ylabel(\'yy\');
PLOT1(GAUSS JACOBI):
PLOT2(GAUSS SEIDEL):
PLOT3(SOR):
SOLVING UNSTEADY(TRANSIENT) STATE 2D HEAT CONDUCTION EQUATION:
The general equation for 2d unsteady state heat conduction equation is
we discretize it we get,
By using the above equation in Matlab code we solve the problem
Grid size:51
time step:0.01
MATLAB CODE FOR UNSTEADY STATE (EXPLICIT METHOD):
clear all
close all
clc
% Inputs
nx = 51;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = 1/(nx-1);
dy = 1/(ny-1);
alpha = 1e-4;
dt = 0.01;
error = 6e9;
tol = 1e-4;
% Initializing the temparature grid
T = 298*ones(nx,ny);
T(1,:) = 900;
T(end,:) = 600;
T(:,1) = 400;
T(:,end) = 800;
% edge refining
T(1,1) = (900+400)/2;
T(1,end) = (900+800)/2;
T(end,1) = (400+600)/2;
T(end,end) = (800+600)/2;
% CFL number calculation and adaptive Time_step control
CFL_Number = (alpha*dt)/(dx^2);
% calcution of teparature distribution explicitly using iterative solvers
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[xx,yy] = meshgrid(x,y);
Told = T;
tic;
counter = 1;
for k = 1:20000
% gauss jacobi explicit iteration
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;
counter = counter+1;
end
time_counter = toc;
figure(1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf(\'2D transient state heat conduction explicit method & No.of iterations:%d (Time: %0.2fs) & Computation time: %0.4f & CFL Number: %0.4f\',counter-1,k*dt,time_counter,CFL_Number));
xlabel(\'xx\');
ylabel(\'yy\');PLOT(GAUSS JACOBI EXPLICIT):
MATLAB CODE FOR UNSTEADY STATE (IMPLICIT METHOD):
Omega is successive over relaxation parameter which is defined as
grid size = 51
time step = 0.01
clear all
close all
clc
% Inputs
nx = 51;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = 1/(nx-1);
dy = 1/(ny-1);
alpha = 1e-4;
t_start = 0;
t_end = 1/alpha;
dt = 0.01;
nt = [t_start:dt*50:t_end];
error = 6e9;
tol = 1e-4;
% Initializing the temparature grid
T = 298*ones(nx,ny);
T(1,:) = 900;
T(end,:) = 600;
T(:,1) = 400;
T(:,end) = 800;
% edge refining
T(1,1) = (900+400)/2;
T(1,end) = (900+800)/2;
T(end,1) = (400+600)/2;
T(end,end) = (800+600)/2;
% calculation of optimal value of over_relaxation parameter(omega)
omega = 2/(1+sin(pi*dx));
% calcution of teparature distribution explicitly using iterative solvers
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[xx,yy] = meshgrid(x,y);
Told = T;
T_prev_dt = Told;
solver = input(\'enter the number for the dedicated solver type: n 1 for gauss jacobi n 2 for gauss siedel n 3 for SOR n\');
tic;
count_jac = 1;
count_seidel = 1;
count_sor = 1;
for k = 1:20000
if solver == 1
term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
error = 1e9;
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;
count_jac = count_jac+1;
end
end
if solver == 2
term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
error = 1e9;
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;
count_seidel = count_seidel+1;
end
end
if solver == 3
term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
error = 1e9;
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
H = (omega*T(i-1,j)+(1-omega)*Told(i-1,j)+Told(i+1,j));
V = (omega*T(i,j-1)+(1-omega)*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;
count_sor = count_sor+1;
end
end
T_prev_dt = T;
end
execution_time = toc;
figure(1);
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
if solver == 1
title(sprintf(\'Temparature profile of unsteady state heat conduction [Gauss Jacobi] Time : %d s & No. of Iterations: %d & Computation Time: %f s\',dt*k,count_jac,execution_time));
else if solver == 2
title(sprintf(\'Temparature profile of unsteady state heat conduction [Gauss Siedel] Time : %d s & No. of Iterations: %d & Computation Time: %f s\',dt*k,count_seidel,execution_time));
else if solver == 3
title(sprintf(\'Temparature profile of unsteady state heat conduction [SOR] Time : %d s & No. of Iterations: %d & Computation Time: %f s & Relaxation Parameter: %f\',dt*k,count_sor,execution_time,real(omega)));
end
end
end
xlabel(\'xx\');
ylabel(\'yy\');
PLOT1(GAUSS JACOBI):
PLOT2(GAUSS SEIDEL):
PLOT3(SOR):
from above plots we can compare the several parameters such as iterations, computation time.