Solving the steady and unsteady 2D heat conduction problem
Objective:
Solve the 2D heat conduction equation by using the point iterative techniques for both Steady State and Transient State Analysis.
Following method should be implemented for solving Implicit Equations
- Jacobi
- Gauss-Seidel
- Successive over-relaxation
Input:
Boundary conditions:
Top Boundary = 600 K
Bottom Boundary = 900 K
Left Boundary = 400 K
Right Boundary = 800 K
Error Criteria = 1e-4
STEADY STATE:
To define a steady state for a system, you’d observe that whatever affects the system’s behaviour are constant or unchanging in time Also, in terms of continuous-time, this means that for those properties of the system, the partial derivative in regards to time remains zero.
For our analysis we have to solve 2D heat conduction problem. Steady state for the problem can be defined as-
`(δ^2T)/(δx^2) + (δ^2T)/(δy^2) = 0`
TRANSIENT STATE:
A transient state is essentially the time between the beginning of the event and the steady state.
Transient state can be defined by stating that if a quantity is at rest and a change in time takes place, thus changing the current state, a transient has occurred.
`(δT)/(δt) = α*((δ^2T)/(δx^2)+(δ^2T)/(δy^2))`
JACOBI ITERATION METHOD:
In numerical linear algebra, the Jacobi method is an iterative algorithm for determining the solutions of a strictly diagonally dominant system of linear equations. Each diagonal element is solved for, and an approximate value is plugged in. The process is then iterated until it converges.
The formula for Jacobi Method for Steady State can be written as:
`T(i,j) = (1/k)*(T_(old)(i-1,j)+T_(old)(i+1,j))/(Δx^2) + (1/k)*(T_(old)(i,j-1)+T_(old)(i,j+1))/(Δy^2)`
where `k = 2*(Δx^2+Δy^2)/(Δx^2*Δy^2)`; `T_(old)` is last iteration's value.
The formula for Jacobi Method for Transient Implicit State can be written as:
`T(i,j)^(n+1) = term_1*T_(prev)(i,j)^n + term_2*(T_(old)(i-1,j)^(n+1)+T_(old)(i+1,j)^(n+1)) + term_3*(T_(old)(i,j-1)^(n+1)+T_(old)(i,j+1)^(n+1))`
where `k1 = (α.ΔT)/(Δx^2)` ; `k2 = (α.ΔT)/(Δy^2)` ; `term_1 = 1/(1+2k1+2k2)` ; `term_2 = k1*term_1` ; `term_3 = k2*term_1`;`T_(old)`is last iteration's value.; `T_(prev)` is last Time Step's value; the value n and n+1, gives the time step
GAUSS-SEIDEL ITERATION METHOD:
Gauss–Seidel method is an iterative method to solve a set of linear equations and very much similar to Jacobi’s method.
Here the better solution is x = (x1, x2, … , xn), if x1(k+1) is a better approximation to the value of x1 than x1(k) is, then it would better that we have found the new value x1(k+1) to use it (rather than the old value that isx1(k)) in finding x2(k+1), … , xn(k+1). So x1(k+1) is found as in Jacobi’s Method, but in finding x2(k+1), instead of using the old value of x1(k) and old values of x3(k),…, xn(k), we then use the new value x1(k+1) and the old values x3(k), … , xn(k), and similarly for finding x3(k+1), … , xn(k+1)
The formula for Gauss-Seidel Method for Steady State can be written as:
`T(i,j) = (1/k)*(T(i-1,j)+T(i+1,j))/(Δx^2) + (1/k)*(T(i,j-1)+T(i,j+1))/(Δy^2)`
where `k = 2*(Δx^2+Δy^2)/(Δx^2*Δy^2)`; `T_(old)` is last iteration's value.
The formula for Transient Implicit State can be written as:
`T(i,j)^(n+1) = term_1*T_(prev)(i,j)^n + term_2*(T(i-1,j)^(n+1)+T(i+1,j)^(n+1)) + term_3*(T(i,j-1)^(n+1)+T(i,j+1)^(n+1))`
where `k1 = (α.ΔT)/(Δx^2)` ; `k2 = (α.ΔT)/(Δy^2)` ; `term_1 = 1/(1+2k1+2k2)` ; `term_2 = k1*term_1` ; `term_3 = k2*term_1`;`T_(old)` is last iteration's value.; `T_(prev)` is last Time Step's value; the value n and n+1, gives the time step
SUCCESSIVE OVER RELAXATION ITERATION METHOD:
The method of successive over-relaxation (SOR) is a variant of the Gauss–Seidel method for solving a linear system of equations, resulting in faster convergence.
For any iterative method, in finding x(k+1) from x(k), we move a certain amount in a particular direction from x(k) to x(k+1). This direction is the vector x(k+1) − x(k), since x(k+1) = x(k) + (x(k+1) − x(k)). If we assume that the direction from x(k) to x(k+1) is taking us closer, but not all the way, to the true solution x , then it would make sense to move in the same direction x(k+1) − x(k), but farther along that direction.
SOR Method uses the value from Gauss_Seidel to converge the results faster with a parameter ω">ω
`T_(SOR)(i,j) = (1-ω)*T_(old)(i,j)+ω*T_(GS)(i,j)`
For ω>1 we have Over-Relaxation and ω<1 we have Under-Relaxation.
CODE FOR STEADY STATE:
%Solving the Steady State 2D heat conduction problem
clear all
close all
clc
%Defining GRID and INPUTS_GRID - 10*10 matrix
nx = 10;
ny = 10;
t = ones(nx, ny);
t(1:nx, 1) = 400;
t(1, 1:ny) = 600;
t(1:nx, ny) = 800;
t(nx, 1:ny) = 900;
t(1,1) = (t(1,2)+t(2,1))/2; %Assigning Corner values
t(nx,1) = (t(nx-1,1)+t(nx,2))/2; %Assigning Corner values
t(1,ny) = (t(1,ny-1)+t(2,ny))/2; %Assigning Corner values
t(nx,ny) = (t(nx-1,ny)+t(nx,ny-1))/2; %Assigning Corner values
t(2:nx-1,2:ny-1) = 300; %Assuming grid values
error = 9e9;
tolerance = 1e-4;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = 1/(nx-1);
dy = 1/(ny-1);
[xx yy] = meshgrid(x,y);
k = 2*(dx^2 + dy^2)/(dx^2 * dy^2);
told = t;
tinitial = t;
jacobi_iteration = 0;
gs_iteration = 0;
sor_iteration = 0;
w = 1.5;
%Selecting method to solve Steady State Problem
method = input('Select the method to solve the Steady State problem n 1. Jacobi n 2. Gauss-seidel n 3. Successive over-relaxation n')
%Code for JACOBI METHOD
if method==1
while error>tolerance
%Space Loop
for i = 2:nx-1
for j = 2:ny-1
H = (told(i-1,j)+told(i+1,j))/(dx^2);
V = (told(i,j-1)+told(i,j+1))/(dy^2);
t(i,j) = (H+V)/k;
end
end
diff = abs(told-t);
error = max(max(diff));
told = t;
jacobi_iteration = jacobi_iteration +1;
end
fprintf('Number of Iterations for Jacobi Method = %d', jacobi_iteration)
%Plot for JACOBI METHOD
contourf(x,y,t,'showtext','on')
colormap(jet)
colorbar
set(gca,'YDIR','reverse')
xlabel('X axis')
ylabel('Y axis')
title(sprintf('Steady Analysis using Jacobi Method'))
end
%Code for Gauss-seidel
if method==2
while error>tolerance
%Space Loop
for i = 2:nx-1
for j = 2:ny-1
H = (t(i-1,j)+told(i+1,j))/(dx^2);
V = (t(i,j-1)+told(i,j+1))/(dy^2);
t(i,j) = (H+V)/k;
end
end
diff = abs(told-t);
error = max(max(diff));
told = t;
gs_iteration = gs_iteration +1;
end
fprintf('Number of Iterations for Gauss-Seidel Method = %d', gs_iteration)
%Plot for Gauss-seidel
contourf(x,y,t,'showtext','on')
colormap(jet)
colorbar
set(gca,'YDIR','reverse')
xlabel('X axis')
ylabel('Y axis')
title(sprintf('Steady Analysis using Gauss-seidel Method'))
end
%Code for Successive over-relaxation
if method==3
while error>tolerance
%Space Loop
for i = 2:nx-1
for j = 2:ny-1
H = (t(i-1,j)+told(i+1,j))/(dx^2);
V = (t(i,j-1)+told(i,j+1))/(dy^2);
tgs(i,j) = (H+V)/k;
t(i,j) = (1-w)*told(i,j) + w*tgs(i,j);
end
end
diff = abs(told-t);
error = max(max(diff));
told = t;
sor_iteration = sor_iteration +1;
end
fprintf('Number of Iterations for Successive over-relaxation method = %d', sor_iteration)
%Plot for Successive over-relaxation
contourf(x,y,t,'showtext','on')
colormap(jet)
colorbar
set(gca,'YDIR','reverse')
xlabel('X axis')
ylabel('Y axis')
title(sprintf('Steady Analysis using Successive Over Relaxation Method'))
end
OUTPUT:
1. Jacobi Method
2. Gauss-Seidel Method
3. SOR Method
OUTPUT PLOTS:
1. Jacobi Method
2. Gauss-Seidel Method
3. SOR Method
CODE FOR TRANSIENT STATE:
%Solving the Transient State 2D heat conduction problem
clear all
close all
clc
%Defining GRID and INPUTS_GRID - 10*10 matrix
nx = 10;
ny = 10;
t = ones(nx, ny);
t(1:nx, 1) = 400;
t(1, 1:ny) = 600;
t(1:nx, ny) = 800;
t(nx, 1:ny) = 900;
t(1,1) = (t(1,2)+t(2,1))/2; %Assigning Corner values
t(nx,1) = (t(nx-1,1)+t(nx,2))/2; %Assigning Corner values
t(1,ny) = (t(1,ny-1)+t(2,ny))/2; %Assigning Corner values
t(nx,ny) = (t(nx-1,ny)+t(nx,ny-1))/2; %Assigning Corner values
t(2:nx-1,2:ny-1) = 300; %Assuming grid values
told = t;
tprev = t;
error = 9e9;
tolerance = 1e-4;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = 1/(nx-1);
dy = 1/(ny-1);
[xx yy] = meshgrid(x,y);
dt = 1e-4;
iteratoin_explicit = 0;
alpha = 1.4;
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
%Selecting method to solve Transient State Problem
method = input('Select the method to solve the Steady State problem n 1. Explicit n 2. Implicit Jacobi n 3. Implicit Gauss-seidel n 4. Implicit Successive over-relaxation n')
% For Explicit Method
if method == 1
for nt = 1:1400
for i = 2:nx-1
for j = 2:ny-1
t(i,j) = told(i,j) + k1*(told(i-1,j)-2*told(i,j)+told(i+1,j)) + k2*(told(i,j-1)-2*told(i,j)+told(i,j+1));
end
end
told = t;
end
contourf(x,y,t,'showtext','on')
colormap(jet)
colorbar
set(gca,'YDIR','reverse')
xlabel('X axis')
ylabel('Y axis')
title(sprintf('Transient Analysis using Explicit Method'))
end
% For Implicit Jacobi Method
if method == 2
iteration_implicit_jacobi = 1;
term1 = 1/(1+2*k1+2*k2);
term2 = k1*term1;
term3 = k2*term1;
for nt = 1:1400
error = 9e9;
while(error>tolerance)
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*tprev(i,j) + term2*H + term3*V;
end
end
diff = abs(told-t);
error = max(max(diff));
told = t;
iteration_implicit_jacobi = iteration_implicit_jacobi + 1;
end
tprev = t;
end
fprintf('Number of Iterations for Implicit Jacobi Method = %d', iteration_implicit_jacobi)
contourf(x,y,t,'showtext','on')
colormap(jet)
colorbar
set(gca,'YDIR','reverse')
xlabel('X axis')
ylabel('Y axis')
title(sprintf('Transient Analysis using Implicit Jacobi Method'))
end
% For Implicit Gauss-Seidel Method
if method == 3
iteration_implicit_GS = 1;
term1 = 1/(1+2*k1+2*k2);
term2 = k1*term1;
term3 = k2*term1;
for nt = 1:1400
error = 9e9;
while(error>tolerance)
for i = 2:nx-1
for j = 2:ny-1
H = t(i-1,j)+t(i+1,j);
V = t(i,j-1)+t(i,j+1);
t(i,j) = term1*tprev(i,j) + term2*H + term3*V;
end
end
iteration_implicit_GS = iteration_implicit_GS + 1;
diff = abs(told-t);
error = max(max(diff));
told = t;
end
tprev = t;
end
fprintf('Number of Iterations for Implicit Gauss-Seidel Method = %d', iteration_implicit_GS)
contourf(x,y,t,'showtext','on')
colormap(jet)
colorbar
set(gca,'YDIR','reverse')
xlabel('X axis')
ylabel('Y axis')
title(sprintf('Transient Analysis using Implicit Gauss-Seidel Method'))
end
% For Implicit SOR Method
if method == 4
iteration_implicit_SOR = 1;
term1 = 1/(1+2*k1+2*k2);
term2 = k1*term1;
term3 = k2*term1;
w = 1.015;
for nt = 1:1400
error = 9e9;
while(error>tolerance)
for i = 2:nx-1
for j = 2:ny-1
H = t(i-1,j)+t(i+1,j);
V = t(i,j-1)+t(i,j+1);
t(i,j) = term1*tprev(i,j) + term2*H + term3*V;
t(i,j) = (1-w)*told(i,j) + w*t(i,j);
end
end
iteration_implicit_SOR = iteration_implicit_SOR + 1;
diff = abs(told-t);
error = max(max(diff));
told = t;
end
tprev = t;
end
fprintf('Number of Iterations for Implicit SOR = %d', iteration_implicit_SOR)
contourf(x,y,t,'showtext','on')
colormap(jet)
colorbar
set(gca,'YDIR','reverse')
xlabel('X axis')
ylabel('Y axis')
title(sprintf('Transient Analysis using Implicit SOR Method'))
endOUTPUT:
1. Implicit Jacobi Method
2. Implicit Gauss-Seidel Method
3. Implicit SOR Method
OUTPUT PLOTS:
1. Explicit Method
2. Implicit Jacobi Method
3. Gauss-Seidel Method
4. SOR Method
RESULT EXPLANATION:
It can be seen that for the Steady State there is significant drop in the number of iterations to converge the results from assumed values for following,
Jacobi (207) > Gauss-Seidel (111) > SOR(28).
For SOR the ω selected was 1.5. For ω value we can check for optimum value by trial and error method, which was employed here for the accuracy of 0.1 and found that 1.5 was optimum owing to its quicker iteration and convergence. Higher values for ω can lead to higher iterations and convergence time.
Similarly, for Transient Implicit State it can be seen that it follows similar trend,
Jacobi (5455) > Gauss-Seidel (4926) > SOR(4391).
For SOR the ω values selected is 1.015. For Transient Implicit Analysis it seems that the SOR is highly sensitive with the ω value as the iterations values starts increasing exponentially with the increase in ω value.