Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
AIM: Solving steady and unsteady 2D heat conduction problem
Theory:
- The steady 2D heat conduction equation is as given :
`frac(del^2T)(delx^2) +frac(del^2T)(dely^2) =0`
- The unsteady 2D heat conduction equation is as given:
`frac(delT)(delt)-alpha(frac(del^2T)(delx^2) +frac(del^2T)(dely^2))=0`
where `alpha` is the thermal diffusion of the system,Temperarture varies in 2 Dimension.
- Why to use iterative solver?
The big advantage of iterative solvers is that they use less memory that direct solver.For example let us take a system of linear equation `A.x=b` where A is the coefficient matrix , x is the variable matrix and b is the constant matrix. Here using direct solver will need to calculate the inverse of A which will take longer time specially when most of the elements of A is zero. However the iterative solver undersgoes only simple algebraic operations such as addition, subtraction, division, and multiplication.
The iterative solver guesses an initial intake for the system and iterates until the error is considerably low.
PROBLEM SPECIFIATIONS AND ASSUMPTIONS:
Consider a 2D square element with boundary conditions for the problem are as follows;
Top Boundary = 600 K
Bottom Boundary = 900 K
Left Boundary = 400 K
Right Boundary = 800 K
Taking whole of interior to be at 300K initially.
Assuming the thermal diffusivity `alpha`=1.2 for the system (only for unsteady system).
to be solved implementing the following methods for solving implicit equations as well as steady state
1. Jacobi
2. Gauss-seidel
3. Successive over-relaxation
SAMPLE CODE FOR STEADY STATE:
close all
clear all
clc
% solving steady 2D heat conduction equation
%(del^2 T/delx^2)+(del^2 T/dely^2) = 0)
nx=10;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
T=300*ones(nx,nx);
dx=x(2)-x(1);
dy=y(2)-y(1);
%boundary condition top=600K,left=400K,bottom=900K,right=800K
T(1,:)=600;
T(:,1)=400;
T(nx,:)=900;
T(:,nx)=800;
T(1,1)=500;
T(1,nx)=700;
T(nx,1)=650;
T(nx,nx)=850;
Told=T;
%solving the 2D steady state using jacobi
tol=1e-4;%min allowed error
error=9e9;% just to enter the loop
%T=((1/K)((Tl+Tr)/dx^2))+((1/k)((Tt+Tb)/dy^2))
K=2*((dx^2 +dy^2)/(dx^2 * dy^2));
method =1;% as per required
if (method==1)
%Convergence loop for JACOBI
jac_iter=1;
while(error>tol)
for i=2:nx-1
for j=2:ny-1
T(i,j)= (1/K)*(((Told(i-1,j)+Told(i+1,j))/dx^2)+((Told(i,j-1)+Told(i,j+1))/dy^2));
end
end
error=max(max(abs(T-Told)));
Told=T;
jac_iter=jac_iter+1;
end
contourf(x,y,T,'Showtext','on')
colormap(jet)
colorbar
set(gca,'YDIR','reverse')
pause(0.03);
end
if(method==2)
%Convergence loop for GAUSS SIDEL
gs_iter=1;
while(error>tol)
for i=2:nx-1
for j=2:ny-1
T(i,j)= (1/K)*(((T(i-1,j)+Told(i+1,j))/dx^2)+((T(i,j-1)+Told(i,j+1))/dy^2));
end
end
error=max(max(abs(T-Told)));
Told=T;
gs_iter=gs_iter+1;
end
contourf(x,y,T,'Showtext','on')
colormap(jet)
colorbar
set(gca,'YDIR','reverse')
end
if(method==3)
%Convergence loop for SOR
sor_iter=1;
w=1.3;%relaxation factor
while(error>tol)
for i=2:nx-1
for j=2:ny-1
T(i,j)= (1-w)*Told(i,j)+(w*(1/K)*(((T(i-1,j)+Told(i+1,j))/dx^2)+((T(i,j-1)+Told(i,j+1))/dy^2)));
end
end
error=max(max(abs(T-Told)));
Told=T;
sor_iter=sor_iter+1;
end
contourf(x,y,T,'Showtext','on')
colormap(jet)
colorbar
set(gca,'YDIR','reverse')
end
SAMPLE CODE FOR UNSTEADY STATE(explicit):
close all
clear all
clc
% solving steady 2D heat conduction equation
%((delT/delt)-alpha*(del^2 T/delx^2)+(del^2 T/dely^2) = 0)
nx=10;
ny=nx;
p=1.2;%p is thermal diffusivity
t=0.6; %simulation run time
nt=1000;%number of time steps
dt=t/nt;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
T=300*ones(nx,nx);
dx=x(2)-x(1);
dy=y(2)-y(1);
%boundary condition top=600K,left=400K,bottom=900K,right=800K
T(1,:)=600;
T(:,1)=400;
T(nx,:)=900;
T(:,nx)=800;
T(1,1)=500;
T(1,nx)=700;
T(nx,1)=650;
T(nx,nx)=850;
k1=p*dt /(dx^2);
k2=p*dt /(dy^2);
method=1;
Told=T;
Tprev=Told;
%solving the 2D steady state
tol=1e-4;%min allowed error
error=9e9;% just to enter the loop
if (method==1)
iter=1;
for m=1:nt
for i=2:nx-1
for j=2:ny-1
T(i,j)= Tprev(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;
Tprev=T;
iter=iter+1;
contourf(x,y,T,'Showtext','on')
colormap(jet)
colorbar
set(gca,'YDIR','reverse')
pause(0.003);
end
endSAMPLE CODE FOR UNSTEADY STATE(implicit):
close all
clear all
clc
% solving steady 2D heat conduction equation
%((delT/delt)-alpha*(del^2 T/delx^2)+(del^2 T/dely^2) = 0)
nx=10;
ny=nx;
p=1.2;%p is thermal diffusivity
t=0.6; %simulation run time
nt=1400;%number of time steps
dt=t/nt;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
T=300*ones(nx,nx);
dx=x(2)-x(1);
dy=y(2)-y(1);
k1= (p*dt)/dx^2 ;
k2= (p*dt)/dy^2 ;
%boundary condition top=600K,left=400K,bottom=900K,right=800K
T(1,:)=600;
T(:,1)=400;
T(nx,:)=900;
T(:,nx)=800;
T(1,1)=500;
T(1,nx)=700;
T(nx,1)=650;
T(nx,nx)=850;
Told=T;
Tprev=T;
term1=1/(1+ 2*k1 + 2*k2);
term2=k1*term1;
term3=k2*term1;
method=1;% as per required
% error=9e9;
tol=1e-4;
%JACOBI
if(method==1)%solving the unsteady state using jacobi as iterative solver
jac_iter=1;
%start of time loop
for m=1:nt
% The error was getting very lower value w.r.t each time step so
% that our solution ended abruptly.
% Here we reinitialized error after each time step so that we
% allowed the appropriate time for our solution to converge.
error=9e9;
%start of convergence loop
while(error>tol)
for i=2:nx-1
for j=2:ny-1
T(i,j)= term1*Tprev(i,j)+ term2*(Told(i-1,j)+Told(i+1,j))+ term3*(Told(i,j-1)+Told(i,j+1));
end
end
error=max(max(abs(T-Told)));
Told=T;
jac_iter=jac_iter+1;
end
Tprev=T;
contourf(x,y,T,'Showtext','on')
colormap(jet)
colorbar
set(gca,'YDIR','reverse')
pause(0.003);
end
end
if(method==2)%solving the unsteady state using guass sidel as iterative solver
gs_iter=1;
%start of time loop
for m=1:nt
error=9e9;
%start of convergence loop
while(error>tol)
for i=2:nx-1
for j=2:ny-1
T(i,j)= term1*Tprev(i,j)+ term2*(T(i-1,j)+Told(i+1,j))+ term3*(T(i,j-1)+Told(i,j+1));
end
end
error=max(max(abs(T-Told)));
Told=T;
gs_iter=gs_iter+1;
end
Tprev=T;
contourf(x,y,T,'Showtext','on')
colormap(jet)
colorbar
set(gca,'YDIR','reverse')
pause(0.003);
end
end
if(method==3)%solving the unsteady state using SOR as iterative solver
sor_iter=1;
w=1.3;%relaxation factor
%start of time loop
for m=1:nt
error=9e9;
%start of convergence loop
while(error>tol)
for i=2:nx-1
for j=2:ny-1
T(i,j)= ((1-w)*term1*Told(i,j))+(w*( term2*(T(i-1,j)+T(i+1,j))+ term3*(T(i,j-1)+T(i,j+1))));
end
end
error=max(max(abs(T-Told)));
Told=T;
sor_iter=sor_iter+1;
end
Tprev=T;
contourf(x,y,T,'Showtext','on')
colormap(jet)
colorbar
set(gca,'YDIR','reverse')
pause(0.003);
end
endUNDERSTANDING THE CODE:
- Steady state :
Since there is'nt a time marching solution the final steady state solution is what we need. The iterative solver however needs to run on loop inorder to calcculate the value of all nodes.Finally the maximum of the absolute of the eigen value of the matrix is calculated as error. If our obtained error is greater than the tolerence value of the system the process should be excecuted again.Hence we require a while loop comparing the error and tolerence value.
The updated value of each element of the matrix through each iteration varies with respect to the type of solver used. the following are the solvers used and the terms .
Jacobian:`T(i,j)=frac(1)(K)*(frac(T_(old)(i-1,j)+T_(old)(i+1,j))(delx^2)+ frac(T_(old)(i,j-1)+T_(old)(i,j+1))(dely^2))` where `K= frac(2*(delx^2+dely^2))(delx^2 *dely^2`
In jacobian the values of newly updated matrix is used only after the present iteration is over.
Gauss Sidel:`T(i,j)=frac(1)(K)*(frac(T(i-1,j)+T_(old)(i+1,j))(delx^2)+ frac(T(i,j-1)+T_(old(i,j+1)))(dely^2))` where `K= frac(2*(delx^2+dely^2))(delx^2 *dely^2`
In Gauss Sidel the value of updated matrix is also been use in the present iteration for other terms.
Sucessive over realaxation(SOR): `T(i,j)=(1-w)*T_(old) + w*frac(1)(K)*(frac(T(i-1,j) +T(i+1,j))(delx^2) + frac(T(i,j-1)+T(i,j+1))(dely^2))` `K= frac(2*(delx^2+dely^2))(delx^2 *dely^2`
In SOR the ralaxation factor plays major role since it's value can quickly converge or never converge the solution.
- Unsteady state (explicit):
Here we have a time marching solution but since it is in explicit form we dont need any iterative solver.
`T(i,j)=Tprev(i,j) + k1*(T_(old)(i-1,j) + 2*T_(old)(i,j) +T_(old)(i,j+1)) + k2*(T_(old)(i,j-1) +2*T_(old)(i,j) +T_(old)(i,j+1))` where `k1=frac(alpha*delt)(delx^2)` and `k2=frac(alpha*delt)(dely^2)`
Here we plot thecontour of temperature profile for each time step inorder to finnal achieve the solution.
- Unsteady state(implicit):
Here we have a marching solution and inside each time step we need to iterate the temperature of each nodes separately since the form is implicit. Since we have three types of iterative solvers we do have three similar equations.
jacobians
`T(i,j)= term1 *Tprev(i,j) + term2 *(T_(old)(i-1,j)+T_(old)(i+1,j)) + term3* (T_(old)(i,j-1) +T_(old)(i,j+1))`
Gauss sidel
`T(i,j)= term1 *Tprev(i,j) + term2 *(T(i-1,j)+T_(old)(i+1,j)) + term3* (T(i,j-1) +T_(old)(i,j+1))`
SOR
`T(i,j)= ((1-w)*term1*T_(old)(i,j) ) + (w*(term2*(T(i-1,j)+T(i+1,j)) +term3*(T(i,j-1)+T(i,j+1)))`
where Tprev is the values of the converged matrix in the previous time step and Told is the values of matrix in past iteration.
POINT to be noticed here is that when the error value should be initialized before every time step failing to which the time marching stops and yields a wrong solution which is also put up in the section of plots.
PLOTS, Results and Conclusions:
- Steady State plot:
It was observed that all the three types of iterative methods gave the same final steady state .Only the number of iterations to reach the steady state differed. We could see that the number of iterations followed a descending pattern as jacobian> gauss sidel> SOR.
- Unsteady state plot:
The above link shows the variation of temperature profile with respect to various time steps. We could see that the number of iterations followed a descending pattern as jacobian> gauss sidel> SOR.
It can be concluded that the SOR method helps us reach the final solution quickly compared to other methods. But the key to solving it using SOR lies in the selction of relaxation factor. If the relaxation factor is too low the time to solve would be greater than even jacobi and is the ralaxation factor is too high then due to over shooting it can take even longer time that jacobi.
Also the value of thermal diffusivity plays a crucial role in the system as lower values would lead to take larger time compared to greater values which actually states that the heat is conducted easily through the system.