Week 5 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Here 2D heat conduction equation is solved for steady and unsteady cases. Various Iterative methods are choosen in order to understand which method converges faster towards the precise result.
Steady state equation:
∂^2T/∂x^2 + ∂^2T/∂y^2 = 0
Unsteady state equation:
α(∂^2T/∂x^2 + ∂^2T/∂y^2) = ∂T/∂t
all the second order terms are converted into algebraic form by central differencing.
For the steady state equation the converted algebraic form for various methods is as folows:
Jacobi method:
T(i,j)=0.25*(Told(i−1,j)+Told(i+1,j)+Told(i,j−1)+Told(i,j+1))
where all the terms are in same time step and all the RHS time variables take values from previous 'j' and update only after all the values are computed along x.
Gauss Seidal method:
T(i,j) = 0.25*(T(i-1,j) + Told(i+1,j) + T(i,j-1) + Told(i,j+1))
Here the previous node values are updated as soon as iteration is done for a (i,j) point for next (i+1,j) point.
Successive over relaxation method:
T(i,j) =(1-omega)*T(i,j)+omega*0.25*(T(i-1,j)+Told(i+1,j)+T(i,j-1)+Told(i,j+1))
omega is found by using this equation:
c = cos(π/nx)+cos(π/ny)
omega = 4/(2+(√(4−c^2)))
For the unsteady state equation there are two more divisions namely:
1. Explicit method
2. Implicit method
An explicit calculation becomes easier due to all know parameters and it doesnt require any kind of itrative solvers.
While the implicit method has many unknowns we require the iterative methods as above mentioned.
Jacobi method:
T_1(i,j) = (T1(i,j)*(1+4*K)^(-1)) + ((T_old1(i-1,j)+T_old1(i+1,j))*K*(1+4*K)^(-1)) + ((T_old1(i,j-1)+T_old1(i,j+1))*K*(1+4*K)^(-1))
Where K is CFL for x and y, since dx and dy are the same.
Gauss Seidal method:
T_1(i,j) = (T1(i,j)*((1+4*K)^(-1))) + ((T_1(i-1,j)+T_old1(i+1,j))*K*((1+4*K)^(-1))) + ((T_1(i,j-1)+T_old1(i,j+1))*K*((1+4*K)^(-1)))
Successive Over relaxation method:
T_1(i,j) = (1-omega)*T1(i,j)+omega*(T_1(i,j)/(1+4*K)+(K/(1+4*K))*(T1(i,j-1)+T_old1(i,j+1)+T1(i-1,j)+T_old1(i+1,j)))
omega is calculated by,
Trial and error method and it is found to be 1.18 which is optimal for this condition.
........................................................................Steady State.............................................................................................
clear all
close all
clc
nx = 11;
ny = 11;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
[X,Y] = meshgrid(x,y);
T = zeros(nx,nx);
T(1,:)=600;
T(nx,:)=900;
T(:,1)=400;
T(:,nx)=800;
T_old = T;
error1 = 10;
error2 = 10;
error3 = 10;
T_real = T;
jacobi_iter = 0;
k1 = (dx^2)/(2*((dx^2)+(dy^2)));
k2 = (dy^2)/(2*((dx^2)+(dy^2)));
T_old1 = T_real;
T1 = T_real;
while error1 > 1e-4
jacobi_iter = jacobi_iter + 1;
for j = 2:(ny-1)
for i = 2:(nx-1)
T1(i,j) = 0.25*(T_old1(i-1,j) + T_old1(i+1,j) + T_old1(i,j-1) + T_old1(i,j+1));
end
end
error1 = max(max(abs(T_old1-T1)));
T_old1 = T1;
figure(1)
[a,b] = contourf(x,y,T_old1);
colorbar;
colormap(jet);
title_2 = sprintf('Jacobi Iteration number = %d',jacobi_iter);
title(title_2);
xlabel('X');
ylabel('Y');
end
T_old2 = T_real;
T2 = T_real;
gs_iter = 0;
while error2 > 1e-4
gs_iter = gs_iter + 1;
for j = 2:(ny-1)
for i = 2:(nx-1)
T2(i,j) = 0.25*(T2(i-1,j) + T_old2(i+1,j) + T2(i,j-1) + T_old2(i,j+1));
end
end
error2 = max(max(abs(T_old2-T2)));
T_old2 = T2;
figure(2)
[a,b] = contourf(x,y,T_old2);
colorbar;
colormap(jet);
title_2 = sprintf('Gauss Seidal Iteration number = %d',gs_iter);
title(title_2);
xlabel('X');
ylabel('Y');
end
T_old3 = T_real;
T3 = T_real;
SOR_iter = 0;
c = cos(pi/nx)+cos(pi/ny);
omega = 4/(2+sqrt(4-(c^2)));
while error3 > 1e-4
SOR_iter = SOR_iter + 1;
for j = 2:(ny-1)
for i = 2:(nx-1)
T3(i,j) =(1-omega)*T3(i,j)+omega*0.25*(T3(i-1,j)+T_old3(i+1,j)+T3(i,j-1)+T_old3(i,j+1));
end
end
error3 = max(max(abs(T_old3-T3)));
T_old3 = T3;
figure(3)
[a,b] = contourf(x,y,T_old3);
colorbar;
colormap(jet);
title_3 = sprintf('SOR Iteration number = %d',SOR_iter);
title(title_3);
xlabel('X');
ylabel('Y');
end
Results for Steady state:
1. Jacobi method took 264 iterations;
2. Gauss Seidal method took 140 iterations
3. SOR method took just 35 iterations
..................................................................Unsteady Explicit solution.................................................................................
clear all
close all
clc
%unsteady stade 2D heat conduction solved using explicit method
alpha = 0.002;
nx = 11;
ny = 11;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dt = 1;
dx = x(2) - x(1);
dy = y(2) - y(1);
[X,Y] = meshgrid(x,y);
K = (alpha*dt)/(dx^2);
T = zeros(nx,nx);
T(1,:)=600;
T(nx,:)=900;
T(:,1)=400;
T(:,nx)=800;
T_old = T;
for t = 1:1400
for j = 2:(nx-1)
for i = 2:(ny-1)
T(i,j) = T_old(i,j)+(K*((T_old(i-1,j)-2*T_old(i,j)+T_old(i+1,j)) + K*(T_old(i,j-1)-2*T_old(i,j)+T_old(i,j+1))));
end
end
T_old = T;
figure(1)
[b,c] = contourf(x,y,T);
colorbar;
title_1 = sprintf('Time step = %d',t);
title(title_1);
xlabel('X');
ylabel('Y');
pause(0.03),
endResults:
Here the explicit solution arrived at the minimal error after 1400 iterations
................................................................................Unsteady Implicit Sloution..................................................................
close all
clear all
clc
alpha = 0.002;
nx = 11;
ny = 11;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
dt = 1;
[X,Y] = meshgrid(x,y);
K = (alpha*dt)/(dx^2);
time = 1400;
T = zeros(nx,ny);
T(1,:)=600;
T(nx,:)=900;
T(:,1)=400;
T(:,nx)=800;
T_old = T;
T_real = T;
jacobi_iter = 0;
T_old1 = T_real;
T_1 = T_real;
T1 = T_real;
for t = 1:time
error1 = 9e9;
while error1 > 1e-4
jacobi_iter = jacobi_iter + 1;
for j = 2:(ny-1)
for i = 2:(nx-1)
T_1(i,j) = (T1(i,j)*(1+4*K)^(-1)) + ((T_old1(i-1,j)+T_old1(i+1,j))*K*(1+4*K)^(-1)) + ((T_old1(i,j-1)+T_old1(i,j+1))*K*(1+4*K)^(-1));
end
end
error1 = max(max(abs(T_1-T_old1)));
T_old1 = T_1;
figure(1)
[a1,b1] = contourf(x,y,T_old1);
colorbar;
colormap(jet);
title_2 = sprintf('Jacobi Iteration number = %d',jacobi_iter);
title(title_2);
xlabel('X');
ylabel('Y');
end
T1 = T_1;
end
T_old2 = T_real;
T_2 = T_real;
gs_iter = 0;
T2 = T_real;
for t = 1:time
error2 = 9e9;
while error2 > 1e-4
gs_iter = gs_iter + 1;
for j = 2:(ny-1)
for i = 2:(nx-1)
T_2(i,j) = (T2(i,j)*((1+4*K)^(-1))) + ((T_2(i-1,j)+T_old2(i+1,j))*K*((1+4*K)^(-1))) + ((T_2(i,j-1)+T_old2(i,j+1))*K*((1+4*K)^(-1)));
end
end
error2 = max(max(abs(T_2-T_old2)));
T_old2 = T_2;
figure(2)
[a2,b2] = contourf(x,y,T_old2);
colorbar;
colormap(jet);
title_2 = sprintf('Gauss Seidal Iteration number = %d',gs_iter);
title(title_2);
xlabel('X');
ylabel('Y');
end
T2 = T_2;
end
T_old3 = T_real;
T3 = T_real;
T_3 = T_real;
SOR_iter = 0;
omega = 1.19;
for t = 1:time
error3 = 9e9;
while error3 > 1e-4
SOR_iter = SOR_iter + 1;
for j = 2:(ny-1)
for i = 2:(nx-1)
T_3(i,j) = (1-omega)*T3(i,j)+omega*(T_3(i,j)/(1+4*K)+(K/(1+4*K))*(T3(i,j-1)+T_old3(i,j+1)+T3(i-1,j)+T_old3(i+1,j)));
end
end
error3 = max(max(abs(T_3-T_old3)));
T_old3 = T_3;
figure(3)
[a3,b3] = contourf(x,y,T_old3);
colorbar;
colormap(jet);
title_3 = sprintf('SOR Iteration number = %d',SOR_iter);
title(title_3);
xlabel('X');
ylabel('Y');
end
T3 = T_3;
endResults:
1. Jacobi method arrived at the solution for t = 1400 after 3893 iterations
2. Gauss Seidal mathod arrived at the solution for t =1400 after 3172 iterations
3. SOR method arrived at the solution after 3059 iterations
Conclusion:
We have finally observed that SOR or Successive over relaxation method arrives at the solution earlier than jacobi or gauss seidal method when the omega value is perfectly selected. Here the omega is taken as 1.18 after several trail and error methods and found to be an optimal value for the solution to converge faster.