Week 5 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Week 5 Challenge
Week 5 - Mid term project
Solving the steady and unsteady 2D heat conduction problem
- Objective
Steady state analysis & Transient State Analysis
Solve the 2D heat conduction equation by using the point iterative techniques that were taught in the class.
- 2D heat conduction Equation
The heat equation is a partial differential equation that describes how the distribution of some quantity (such as heat) evolves over time in a solid medium, as it spontaneously flows from places where it is higher towards places where it is lower. It is a special case of the diffusion equation.
`(∂T)/(∂t) =α ( (∂^2 T)/(∂x^2 ) + (∂^2 T)/(∂y^2 ))`
- Steady state analysis
Steady-state conduction is the form of conduction that happens when the temperature difference driving the conduction are constant, so that, the spatial distribution of temperatures in the conducting object does not change any further. Thus, all partial derivatives of temperature concerning space may either be zero or have nonzero values, but all derivatives of temperature at any point concerning time are uniformly zero. In steady-state conduction, the amount of heat entering any region of an object is equal to the amount of heat coming out.
For 2D, steady state and without heat generation, the above equation reduces to:
`(∂^2 T)/(∂x^2 ) + (∂^2 T)/(∂y^2 ) = 0`
This is Laplace’s equation.
- Transient analysis
The equation for transient/unsteady heat conduction in 2 dimension is as below
`(∂T)/(∂t) =α ( (∂^2 T)/(∂x^2 ) + (∂^2 T)/(∂y^2 ))`
Equation needs 2 boundary conditions in each direction. We want to find steady state solutions to that satisfy the Dirichlet boundary conditions
T(x, 0) = f1,
T(x, b) = f2,
0 < x < a
T(0, y) = g1,
T(a, y) = g2,
0 < y < b
The problem of finding a solution is known as the Dirichlet problem.
- Difference equation and discretization
- The difference equation
The difference equation for steady state and transient is based on the first derivative for time and second derivative for space variables
For steady state equation, the difference equation is
`(T_(i-1,j)^n - 2 T_(i,j )^n+ T_(i+1,j )^n)/(Δx^2 ) + (T_(i-1,j )^n- 2 T_(i,j )^n +T_(i+1,j)^n)/(Δy^2 ) = 0`
For transient/unsteady equation, the difference equation is
`(T_(i,j)^(n+1) - T_(i,j )^n )/(Δt) = (T_(i-1,j)^n - 2 T_(i,j )^n+ T_(i+1,j )^n)/(Δx^2 ) + (T_(i-1,j )^n- 2 T_(i,j )^n + T_(i+1,j)^n)/(Δy^2 )`
- Space discretization for steady state equation
To discretize this equation in space and time we will Central Difference scheme for space derivative. The partial derivative can be approximated using the derivative definition, removing the limit operator:
`((∂^2 T)/(∂x^2 ) ) = (T_(i-1,j)^n - 2 T_(i,j )^n+ T_(i+1,j )^n)/(Δx^2 )`
`((∂^2 T)/(∂y^2 ) ) = (T_(i-1,j )^n- 2 T_(i,j )^n + T_(i+1,j)^n)/(Δy^2 )`
`(∂^2 T)/(∂x^2 ) + (∂^2 T)/(∂y^2 ) = 0 = (T_(i-1,j)^n - 2 T_(i,j )^n+ T_(i+1,j )^n)/(Δx^2 ) + (T_(i-1,j )^n- 2 T_(i,j )^n + T_(i+1,j)^n)/(Δy^2 )`
Assuming `k = (2*(Δx^2+Δy^2))/(Δx^2*Δy^2 )`
`T_(i,j )^n= (k*(T_(i-1,j)^n + T_(i+1,j )^n))/(Δx^2 ) +(k*( T_(i-1,j )^n + T_(i+1,j)^n))/(Δy^2 )`
This means that for a given iteration for T, numerical values for 4 points around that point is needed.
- Space discretization for transient equation
To discretize this equation in space and we will Central Difference scheme for the second derivative for space variables.
`((∂^2 T)/(∂x^2 ) ) =≈ (T_(i-1,j)^n - 2 T_(i,j )^n+ T_(i+1,j )^n)/(Δx^2 )`
`((∂^2 T)/(∂y^2 ) ) =≈ (T_(i-1,j )^n- 2 T_(i,j )^n + T_(i+1,j)^n)/(Δy^2 )`
`(∂^2 T)/(∂x^2 ) + (∂^2 T)/(∂y^2 ) = (T_(i-1,j)^n - 2 T_(i,j )^n+ T_(i+1,j )^n)/(Δx^2 ) + (T_(i-1,j )^n- 2 T_(i,j )^n + T_(i+1,j)^n)/(Δy^2 )`
- Time discretization for transient equation
To discretize this equation in space and time we will forward difference scheme for time derivative.
`(∂T)/(∂t ) = (T_(i,j)^(n+1) - T_(i,j )^n )/(Δt)`
- Boundary conditions
Top Boundary = 600 K
Bottom Boundary = 900 K
Left Boundary = 400 K
Right Boundary = 800 K
The temperature at the corners are taken as average of the intersecting values.
- Explicit technique and Implicit methods
In an explicit approach, each difference equation contains only one unknown and there fore can be solved explicitly.
In an implicit method, the unknowns must be obtained by means of simultaneous solution of the difference equation applied at all grid points at any particular time level.
For steady state equation, the difference equation is
`(T_(i-1,j)^n - 2 T_(i,j )^n+ T_(i+1,j )^n)/(Δx^2 ) + (T_(i-1,j )^n- 2 T_(i,j )^n + T_(i+1,j)^n)/(Δy^2 ) = 0`
`T_(i,j )^(n+1)= (Δx^2*Δy^2)/(2*(Δx^2+Δy^2 ) )*( ((T_(i-1,j)^n + T_(i+1,j )^n))/(Δx^2 ) +(( T_(i-1,j )^n + T_(i+1,j)^n ))/(Δy^2 ) )`
Assuming k =`(Δx^2*Δy^2)/(2*(Δx^2+Δy^2 ) )`
- Explicit scheme for solving steady state equation
`T_(i,j )^(n)= k*((T_(i-1,j)^n + T_(i+1,j )^n))/(Δx^2 ) + k* (( T_(i-1,j )^n + T_(i+1,j)^n ))/(Δy^2 )`
This means that for a given iteration for T, numerical values for 4 points around that point from the previous iteration is needed all of which are knowns.
For transient/unsteady equation, the difference equation is
`(T_(i,j)^(n+1) - T_(i,j )^n )/(Δt)= (T_(i-1,j)^n - 2 T_(i,j )^n+ T_(i+1,j )^n)/(Δx^2 ) + (T_(i-1,j )^n- 2 T_(i,j )^n + T_(i+1,j)^n)/(Δy^2 )`
`k1 = (α*(Δt))/(Δx^2 )`
`k2 = (α*(Δt))/(Δy^2 )`
`k3 = 1/(1+2*k1+2*k2)`
- Explicit scheme for solving transient equation
`T_(i,j )^(n+1)= T_(i,j )^n+ k1*(T_(i-1,j)^n - 2*T_(i,j )^n+ T_(i+1,j )^n) + k2*( T_(i-1,j )^n - 2*T_(i,j )^n+ T_(i+1,j)^n)`
This means that for a given iteration for T, numerical values for 4 points around that point from the previous iteration is needed and also the previous value for T.
- Implicit scheme for solving transient equation
`T_(i,j )^(n+1)= k3*T_(i,j )^n+k3*k1*(T_(i-1,j)^(n+1) + T_(i+1,j )^(n+1)) +k3*k2* ( T_(i-1,j )^(n+1) + T_(i+1,j)^(n+1))`
This means that for a given iteration for T, numerical values for 4 points around that point from the present iteration is needed and also the previous value for T.
- Comparison between explicit and implicit methods
|
Explicit scheme |
Implicit scheme |
|
Relatively easy to set up and program |
More complicated to set up and program. Massive matrix operations are need to run every iteration so the time taken per time step may be more. |
|
Δt and Δx must be set according to the stability requirements. If Δt is very small for a problem, the computational time can be very long. |
Stability can be maintained over much larger values of Δt. However, if Δt is large, the truncation error can also be large. |
|
Good for transient analysis |
Good for steady state analysis |
- Methods for solving implicit numerical equations
- Jacobi method
The unknown value of a function f(x,y) at a given grid point for the n+1 iteration depends only on the known values of f(x,y) from the neighbouring grid points from the previous iteration. In Jacobi method, the value of function at all the grid points are evaluated for an iteration for (n+1) iteration before moving on the iteration (n+2).
- Gauss-seidel method
The unknown value of a function f(x,y) at a given grid point for the n+1 iteration depends only on the known values of f(x,y) from the neighbouring grid points from the previous iteration and also from the present iteration (if value is available in the present iteration) In Guass-Siedel method, the values of the present iteration is also used to find the value at the neighbouring point.
- Successive over-relaxation method
Let us assume that convergence is obtained when the error is within tolerance at
`f_(i,j )^(n+1)= f_(i,j )^n`
A relaxation factor is used to enhance the convergence of a solution fro Jacobi Method of Guass-Siedel method.
`f_(i,j )^(n+1)= f_(i,j )^n +ω * (f_(i,j )^(n+1)- f_(i,j )^n )`
Where `f_(i,j )^(n+1)` is an intermediate value obtained from the difference equation.
`f_(i,j )^(n+1)= ω * (f_(i,j )^(n+1) )-(1-ω )* f_(i,j )^n`
If w > 1 ,the technique is called successive over relaxation and id w<1 , the technique is called under relaxation. Generally 1 < ω <2
`ω= 2/((1+sin(pi*dx)) )`
- Methodology / Coding for 2D Steady state Heat equation
We will implement the explicit methods for difference equations to solve second order 2 dimensional steady state heat conduction equation
The space domain is considered as1 m x 1m and is discretised with points each using linspace command.
nx = ny = n = 21
dx= dy =L/(n-1) = 1/(21-1)=1/20=0.05
The temperature field is initialized assuming the initial temperature is 273K across a 21 X 21 matrix that represents the grid. Boundary conditions are applied.
Top Boundary = 600 K
Bottom Boundary = 900 K
Left Boundary = 400 K
Right Boundary = 800 K
The temperature at the corners are taken as average of the intersecting values.
The input profile is plotted using ‘meshgrid’ and ‘contourf’ commands with labels and colour bars.
Figure 1: input temperature profile
Figure 1 shows a value of 273k inside the field whereas the value at the boundaries are as per the boundary conditions
Temperature field values are initiated assuming that the old values of T is same as the input values. It is understood that the old values are to be updated after each iteration, otherwise the old values even after each iteration will remain the same inside the code.
Initial error is assumed as 9e9. And a tolerance of 1e-4 is considered acceptable.
- Jacobi method
Jacobi method for 2D problem is implemented using 2 ‘for’ loop for x and y direction. The value of temperature filled is computed using the old value of the temperature field.
`k = (Δx^2*Δy^2)/(2*(Δx^2+Δy^2 )`
`term1= k* ((T_(i-1,j)^n + T_(i+1,j )^n))/(Δx^2 )`
`term2= k*(( T_(i-1,j )^n + T_(i+1,j)^n ))/(Δy^2 )`
`T_(i,j )^(n)= term1 +term2`
For the first iteration, the old value is taken as the initial value. The temperature filed is updated using a ‘while’ loop until the error which is the maximum of the absolute difference between T computed and old value of T. once the error is found, the old value of T is updated with the new values and the iteration is run again until the error becomes less than the tolerance.
Once a converged solution is obtained the solution is plotted using ‘meshgrid’ and ‘contourf’ commands with labels and colour bars. The plot is saved along with the number of iterations it took to achieve convergence.
For Explicit method, a solution is achieved after 913 iterations within the tolerance for error 1e-4 when we use 21 x 21 grid points
Figure 2: temperature profile from Explicit method for steady state
Figure 2 shows thermal diffusion into the field according to the boundary conditions and shows that the temperature field reaching to a steady state after 913 iterations. The temperature at the boundary is still preserved but the temperature inside the field is redistributed compared to the input profile
iterative_solvers_comparison_steady_state.m
% code to solve second order 2 dimensional pde
%steady state heat conduction
% %d^T/dx^2 + d^T/dy^2 = 0
%x=0 T =400
%x=1 T =800
%y=0 T =900
%y=1 T =600
clear all;
close all;
clc;
% space domain inputs
Lx=1; %Ddomain length in x
Ly=1; %Ddomain length in y
nx= 21; % no. of points in x
ny =21; % no. of points in y
x=linspace(0,Lx,nx);
y=linspace(0,Ly,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
k=2*(dx^2+dy^2)/(dx^2*dy^2);
%initiallizing temeprature field and setting boundary conditions
T=273.*ones(nx,ny);
T_t = 600; % Temperature at top
T_b = 900; % Temperature at bottom
T_l = 400; % Temperature at left
T_r = 800; % Temperature at right
T(2:ny-1,1) = T_l ;
T(2:ny-1,nx)= T_r ;
T(1,2:nx-1) = T_b ;
T(ny, 2:nx-1) = T_t ;
T(1,1) = (T_b + T_l)/2;
T(nx,ny) = (T_t + T_r)/2;
T(nx,1) = (T_t + T_l)/2;
T(1,ny) = (T_r + T_b)/2;
[X,Y] = meshgrid(x,y);
Told=T;
##%plotting input profile
figure(1)
[c,h] = contourf(X, Y, T);
xlabel('x domain (T at top = 600k, T at bottom = 900k)')
ylabel(' domain (T at left = 400k, T at right = 800k)')
title('input temperature profile')
colormap(jet)
colorbar
saveas(1,"input temperature profile.png")
iterative_solver = 1;
error =9e9;
tol = 1e-4;
%Explicit method
if iterative_solver ==1
iter = 1;
while(error > tol)
for i=2:nx-1
for j=2:ny-1
term1= (Told(i-1,j) + Told(i+1,j))/(k*dx^2);
term2= (Told(i,j-1) + Told(i,j+1))/(k*dy^2);
T(i,j)= term1 + term2;
endfor
endfor
error=max(max(abs(Told - T)));
Told=T;
iter = iter + 1
end
figure(2)
[c,h] = contourf(X, Y, T);
xlabel('x domain')
ylabel('y domain')
clabel (c, h,"fontsize", 12);
colormap(jet)
colorbar
title_text=sprintf('Explicit Method for 2D steady state heat conduction iteration number = %d', iter);
title(title_text)
saveas(2,"temperature profile with explicit method for steady state.png")
pause(0.03)
endif- Methodolgy /coding for 2D Transient Heat equation explicitly
We will implement the explicit methods for second order 2 dimensional transient heat conduction equation
The time domain is taken with dt=0.01. The space domain is considered as1 m x 1m and is discretised with points each using linspace command.
nx = ny = n = 21
dx= dy =L/(n-1) = 1/(21-1)=1/20=0.05
The temperature field is initialized assuming the initial temperature is 273K across a 21 X 21 matrix that represents the grid. Boundary conditions are applied.
Top Boundary = 600 K
Bottom Boundary = 900 K
Left Boundary = 400 K
Right Boundary = 800 K
The temperature at the corners are taken as average of the intersecting values.
The input profile is plotted using ‘meshgrid’ and ‘contourf’ commands with labels and colour bars.
Figure 3: input temperature profile
Figure 4 shows a value of 273 k inside the field whereas the value at the boundaries are as per the boundary conditions
Temperature field values are initiated assuming that the old values of T is same as the input values. It is understood that the old values are to be updated after each iteration, otherwise the old values even after each iteration will remain the same inside the code.
Initial error is assumed as 9e9. And a tolerance of 1e-4 is considered acceptable.
- Explicit method
Explicit method for 2D problem is implemented using 2 ‘for’ loop for x and y direction. The value of temperature filled is computed using the old value of the temperature field.
`term1= T_(i,j)^n`
`term2= (T_(i-1,j)^n - 2*T_(i,j )^n+ T_(i+1,j )^n)`
`term3=( T_(i-1,j )^n - 2*T_(i,j )^n+ T_(i+1,j)^n)`
`k1 = (α*(dt))/(dx)^2 , k2 = (α*(dt))/(dy)^2`
`T_(i,j )^(n+1)= term1 +k1*term2+ k2*term3`
For the first iteration, the old value is taken as the initial value. The temperature filed is updated using a ‘while’ loop until the error which is the maximum of the absolute difference between T computed and old value of T. once the error is found, the old value of T is updated with the new values and the iteration is run again until the error becomes less than the tolerance.
Once a converged solution is obtained the solution is plotted using ‘meshgrid’ and ‘contourf’ commands with labels and colour bars. The plot is saved along with the number of iterations it took to achieve convergence.
With Jacobi method for 2D unsteady heat conduction equation, an explicit solution is achieved after 4798 iterations within the tolerance for error 1e-4 when we use 21 x 21 grid points
Figure 4: temperature field with explicit method for 2D unsteady heat conduction
Figure 4 shows thermal diffusion into the field according to the boundary conditions and shows that the temperature field reaching to a steady state after 4798 iterations. The temperature at the boundary is still preserved but the temperature inside the field is redistributed compared to the input profile
iterative_solvers_comparison_transient_explicit.m
% code to solve second order 2 dimensional pde
%transient heat conduction
% %alpha*(d^T/dx^2 + d^T/dy^2) = dT/dt
%x=0 T =400
%x=1 T =800
%y=0 T =900
%y=1 T =600
alpha=1e-4
clear all;
close all;
clc;
% space domain inputs
Lx=1; %Ddomain length in x
Ly=1; %Ddomain length in y
nx= 21; % no. of points in x
ny =21; % no. of points in y
x=linspace(0,Lx,nx);
y=linspace(0,Ly,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
dt=1;
alpha = 1e-4;
k1=alpha*(dt)/(dx^2);
k2=alpha*(dt)/(dy^2);
%initiallizing temeprature field and setting boundary conditions
T=273.*ones(nx,ny);
T_t = 600; % Temperature at top
T_b = 900; % Temperature at bottom
T_l = 400; % Temperature at left
T_r = 800; % Temperature at right
T(2:ny-1,1) = T_l ; % Temperature at left
T(2:ny-1,nx)= T_r ; % Temperature at right
T(1,2:nx-1) = T_b ; % Temperature at bottom
T(ny, 2:nx-1) = T_t ; % Temperature at top
T(1,1) = (T_b + T_l)/2;
T(nx,ny) = (T_t + T_r)/2;
T(nx,1) = (T_t + T_l)/2;
T(1,ny) = (T_r + T_b)/2;
[X,Y] = meshgrid(x,y);
Told=T;
##%plotting input profile
figure(1)
[c,h] = contourf(X, Y, T);
xlabel('x domain (T at top = 600k, T at bottom = 900k)')
ylabel(' domain (T at left = 400k, T at right = 800k)')
title('input temperature profile')
colormap(jet)
colorbar
saveas(1,"input temperature profile.png")
error =9e9;
tol = 1e-4;
%explicit_method
explicit_iter =1;
while(error > tol)
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;
endfor
endfor
error=max(max(abs(Told - T)));
Told=T;
explicit_iter = explicit_iter + 1
end
figure(2)
[c,h] = contourf(X, Y, T);
xlabel('x domain')
ylabel('y domain')
clabel (c, h,"fontsize", 12);
colormap(jet)
colorbar
title_text=sprintf('Explicit method for 2D unsteady heat conduction iteration number = %d', explicit_iter);
title(title_text)
saveas(2,"temperature profile with explicit method for 2D unsteady heat conduction.png")
pause(0.03)
- Methodolgy /coding for 2D Transient Heat equation implicitly
We will implement the implicit technique for second order 2 dimensional transient heat conduction equation
The time domain is taken with dt=0.01. The space domain is considered as1 m x 1m and is discretised with points each using linspace command.
nx = ny = n = 21
dx= dy =L/(n-1) = 1/(21-1)=1/20=0.05
The temperature field is initialized assuming the initial temperature is 273K across a 21 X 21 matrix that represents the grid. Boundary conditions are applied.
Top Boundary = 600 K
Bottom Boundary = 900 K
Left Boundary = 400 K
Right Boundary = 800 K
The temperature at the corners are taken as average of the intersecting values.
The input profile is plotted using ‘meshgrid’ and ‘contourf’ commands with labels and colour bars.
Figure 5: input temperature profile
Figure 5 shows a value of 273 k inside the field whereas the value at the boundaries are as per the boundary conditions
Temperature field values are initiated assuming that the old values of T is same as the input values. It is understood that the old values are to be updated after each iteration, otherwise the old values even after each iteration will remain the same inside the code.
Initial error is assumed as 9e9. And a tolerance of 1e-4 is considered acceptable.
- Jacobi method
Jacobi method for 2D problem is implemented using 2 ‘for’ loop for x and y direction. The value of temperature filled is computed using the old value of the temperature field.
`k1 = (α*(dt))/(dx)^2`
`k2 = (α*(dt))/(dy)^2 &`
`k3 = 1/((1+2*k1+2*k2)`
`term1= k3*k1*(T_(i-1,j)^n + T_(i+1,j )^n)`
`term2=k3*k2* ( T_(i-1,j )^n + T_(i+1,j)^n)`
`T_(i,j )^(n+1)= k3* (T_(i,j)^n)prev + term1 +term2`
The important point for implicit scheme is to notice that the final equation contains the value of T at the previous iteration not the old value. For the first iteration, the old value is taken as the initial value. The temperature filed is updated using two ‘while’ loop until the error which is the maximum of the absolute difference between T computed and old value of T and the error for time domain previous value of T and T. once the 1st error is found, the old value of T is updated with the new values and the iteration is run again until the error becomes less than the tolerance. then 2nd error is found, the previous value of T is updated with the new values and the iteration is run again until the error becomes less than the tolerance.
Once a converged solution is obtained the solution is plotted using ‘meshgrid’ and ‘contourf’ commands with labels and colour bars. The plot is saved along with the number of iterations it took to achieve convergence.
With Jacobi method for 2D unsteady heat conduction equation, an explicit solution is achieved after 23105 iterations within the tolerance for error 1e-4 when we use 21 x 21 grid points
Figure 7: temperature profile with implicit Jacobi method for 2D unsteady heat conduction
Figure 7 shows thermal diffusion into the field according to the boundary conditions and shows that the temperature field reaching to a steady state after 23105 iterations. The temperature at the boundary is still preserved but the temperature inside the field is redistributed compared to the input profile
- Gauss-seidel method
Guass-Siedel method for 2D problem is implemented using 2 ‘for’ loop for x and y direction. The values of temperature filled is computed using the values from present iteration of the temperature field. If present values are not available, then old values are taken for those grids.
`k1 = (α*(dt))/(dx)^2`
`k2 = (α*(dt))/(dy)^2 &`
`k3 = 1/((1+2*k1+2*k2)`
`term1= k3*k1*(T_(i-1,j)^(n+1) + T_(i+1,j )^n)`
`term2=k3*k2* ( T_(i-1,j )^(n+1) + T_(i+1,j)^n)`
`T_(i,j )^(n+1)= k3* (T_(i,j)^n)prev + term1 +term2`
The important point for implicit scheme is to notice that the contains the value of T at the previous iteration not the old value. For the first iteration, the old value is taken as the initial value. The temperature filed is updated using two ‘while’ loop until the error which is the maximum of the absolute difference between T computed and old value of T and the error for time domain previous value of T and T. once the 1st error is found, the old value of T is updated with the new values and the iteration is run again until the error becomes less than the tolerance. then 2nd error is found, the previous value of T is updated with the new values and the iteration is run again until the error becomes less than the tolerance.
Once a converged solution is obtained the solution is plotted using ‘meshgrid’ and ‘contourf’ commands with labels and colour bars. The plot is saved along with the number of iterations it took to achieve convergence.
For Guass-seidel method, a solution is achieved after 20696 iterations within the tolerance for error 1e-4 when we use 21 x 21 grid points
Figure 8: temperature profile with implicit Guass-siedel method for 2D unsteady heat conduction
Figure 8 shows thermal diffusion into the field according to the boundary conditions and shows that the temperature field reaching to a steady state after 20696 iterations. The temperature at the boundary is still preserved but the temperature inside the field is redistributed compared to the input profile
- Successive over-relaxation method
Guass-Siedel method for 2D problem is implemented using 2 ‘for’ loop for x and y direction. The value of temperature filed is computed using the old value of the temperature field. Successive over-relaxation method is applied by modifying the final sum as
`k1 = (α*(dt))/(dx)^2 ,`
`k2 = (α*(dt))/(dy)^2`
`k3 = 1/((1+2*k1+2*k2)) &`
`ω= 2/((1+sin(pi*dx))`
`term1= k3*k1*(T_(i-1,j)^(n+1) + T_(i+1,j )^n)`
`term2=k3*k2* ( T_(i-1,j )^(n+1) + T_(i+1,j)^n)`
`T_(i,j )^(n+1)= ω*(k3* (T_(i,j)^n)prev + term1 +term2)+ (ω-1)*(T_(i,j )^n)prev`
The important point for implicit scheme is to notice that the contains the value of T at the previous iteration not the old value. For the first iteration, the old value is taken as the initial value. The temperature filed is updated using two ‘while’ loop until the error which is the maximum of the absolute difference between T computed and old value of T and the error for time domain previous value of T and T. once the 1st error is found, the old value of T is updated with the new values and the iteration is run again until the error becomes less than the tolerance. then 2nd error is found, the previous value of T is updated with the new values and the iteration is run again until the error becomes less than the tolerance.
Once a converged solution is obtained the solution is plotted using ‘meshgrid’ and ‘contourf’ commands with labels and colour bars. The plot is saved along with the number of iterations it took to achieve convergence.
For Guass-seidel method, a solution is achieved after 13040 iterations within the tolerance for error 1e-4 when we use 21 x 21 grid points
Figure 9: temperature profile with implicit Guass-siedel method with SOR for 2D unsteady heat conduction
Figure 9 shows thermal diffusion into the field according to the boundary conditions and shows that the temperature field reaching to a steady state after 13040 iterations. The temperature at the boundary is still preserved but the temperature inside the field is redistributed compared to the input profile
iterative_solvers_comparison_transient_implicit.m
% code to solve second order 2 dimensional pde
%transient heat conduction
% %alpha*(d^T/dx^2 + d^T/dy^2) = dT/dt
%x=0 T =400
%x=1 T =800
%y=0 T =900
%y=1 T =600
alpha=1e-4
clear all;
close all;
clc;
% space domain inputs
Lx=1; %Ddomain length in x
Ly=1; %Ddomain length in y
nx= 21; % no. of points in x
ny =21; % no. of points in y
x=linspace(0,Lx,nx);
y=linspace(0,Ly,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
dt=1;
nt=200;
alpha = 1e-4;
k1=alpha*(dt)/(dx^2);
k2=alpha*(dt)/(dy^2);
k=1/(1+ 2*k1 + 2*k2);
%initiallizing temeprature field and setting boundary conditions
T=273.*ones(nx,ny);
T_t = 600; % Temperature at top
T_b = 900; % Temperature at bottom
T_l = 400; % Temperature at left
T_r = 800; % Temperature at right
T(2:ny-1,1) = T_l ; % Temperature at left
T(2:ny-1,nx)= T_r ; % Temperature at right
T(1,2:nx-1) = T_b ; % Temperature at bottom
T(ny, 2:nx-1) = T_t ; % Temperature at top
T(1,1) = (T_b + T_l)/2;
T(nx,ny) = (T_t + T_r)/2;
T(nx,1) = (T_t + T_l)/2;
T(1,ny) = (T_r + T_b)/2;
[X,Y] = meshgrid(x,y);
Told = T;
T_prev = Told;
T2 = T;
##%plotting input profile
figure(1)
[c,h] = contourf(X, Y, T);
xlabel('x domain (T at top = 600k, T at bottom = 900k)')
ylabel(' domain (T at left = 400k, T at right = 800k)')
title('input temperature profile')
colormap(jet)
colorbar
saveas(1,"input temperature profile.png")
iterative_solver = 1;
error1 =9e9;
error2 = 9e9;
tol = 1e-4;
%implicit jacobi_method
if iterative_solver ==1
implicit_iter = 0;
implicit_jacobi_iter =1;
while (error2 > tol)
implicit_jacobi_iter =implicit_iter + 1;
error1 = 9e9;
while(error1 > tol)
for i=2:nx-1
for j=2:ny-1
term1(i,j) = k*k1*(Told(i-1,j) + Told(i+1,j));
term2(i,j) = k*k2*(Told(i,j-1) + Told(i,j+1));
T(i,j) = k*(T_prev(i,j)) + term1(i,j) + term2(i,j);
endfor
endfor
error1=max(max(abs(Told - T)));
Told=T;
implicit_jacobi_iter = implicit_jacobi_iter + 1;
end
error2 = max(max(abs(T_prev - T)));
T_prev = T;
implicit_iter = implicit_jacobi_iter;
end
figure(2)
[C,h] = contourf(X, Y, T);
colorbar;
xlabel('x domain');
ylabel('y domain');
clabel (C, h,"fontsize", 12);
colormap(jet);
title_text=sprintf('Implicit Jacobi method for 2D unsteady heat conduction iteration number = %d', implicit_jacobi_iter);
title(title_text);
saveas(2,"temperature profile with implicit Jacobi method for 2D unsteady heat conduction.png");
end
%Implicit Guass-siedel
if iterative_solver ==2
implicit_iter = 0;
implicit_guass_siedel_iter =1;
while (error2 > tol)
implicit_guass_siedel_iter =implicit_iter + 1;
error1 = 9e9;
while(error1 > tol)
for i=2:nx-1
for j=2:ny-1
term1(i,j) = k*k1*(T(i-1,j) + T(i+1,j));
term2(i,j) = k*k2*(T(i,j-1) + T(i,j+1));
T(i,j) = k*(T_prev(i,j)) + term1(i,j) + term2(i,j);
endfor
endfor
error1=max(max(abs(Told - T)));
Told=T;
implicit_guass_siedel_iter = implicit_guass_siedel_iter + 1;
end
error2 = max(max(abs(T_prev - T)));
T_prev = T;
implicit_iter = implicit_guass_siedel_iter
end
figure(3)
[c,h] = contourf(X, Y, T);
xlabel('x domain')
ylabel('y domain')
clabel (c, h,"fontsize", 12);
colormap(jet)
colorbar
title_text=sprintf('Implicit Guass-siedel method for 2D unsteady heat conduction iteration number = %d', implicit_guass_siedel_iter);
title(title_text)
saveas(3,"temperature profile with implicit Guass-siedel method for 2D unsteady heat conduction.png")
pause(0.03)
endif
w= 2/(1+sin(pi*dx));
%Guass-siedel with successive over relaxation ( SOR)
if iterative_solver ==3
implicit_iter = 0;
implicit_guass_siedel_iter =1;
while (error2 > tol)
implicit_guass_siedel_iter =implicit_iter + 1;
error1 = 9e9;
while(error1 > tol)
for i=2:nx-1
for j=2:ny-1
term1(i,j) = k*k1*(T(i-1,j) + T(i+1,j));
term2(i,j) = k*k2*(T(i,j-1) + T(i,j+1));
T(i,j) = w*(k*(T_prev(i,j)) + term1(i,j) + term2(i,j)) + (1-w)*T_prev(i,j) ;
endfor
endfor
error1 = max(max(abs(Told - T)));
Told=T;
implicit_guass_siedel_iter = implicit_guass_siedel_iter + 1;
end
error2 = max(max(abs(T_prev - T)));
T_prev = T;
implicit_iter = implicit_guass_siedel_iter
end
figure(4)
[c,h] = contourf(X, Y, T);
xlabel('x domain')
ylabel('y domain')
clabel (c, h,"fontsize", 12);
colormap(jet)
colorbar
title_text=sprintf('Implicit Guass-siedel method with SOR for 2D unsteady heat conduction iteration number = %d', implicit_guass_siedel_iter);
title(title_text)
saveas(4,"temperature profile with implicit Guass-siedel method with SOR for 2D unsteady heat conduction.png")
pause(0.03)
endif
- Conclusions
- Steady state analysis & Transient State Analysis with 2D heat conduction equation by using the point iterative techniques were carried out.
- Octave program is written for solving the differential equation numerically.
- Steady state equations are solved explicitly.
- Transient state equations are solved explicitly and implicitly using Jacobi method, Guass-Siedel method and Guass-Siedel method with SOR.
- Faster convergence for transient state equation is observed for explicit techniques. implicit technique takes more time to provide solution.
- No. of iterations are lesser with SOR technique.
|
No. of iterations |
error |
|
|
Steady state equation for 2D heat conduction |
||
|
Explicit technique |
913 |
1e-4 |
|
Transient state equation for 2D heat conduction |
||
|
Explicit technique |
4798 |
1e-4 |
|
Implicit technique |
||
|
Jacobi method |
23105 |
1e-4 |
|
Guass-Siedel method |
20696 |
1e-4 |
|
Guass-Siedel method with SOR |
13040 |
1e-4 |