Week 5 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Steady and transient analysis of 2D heat conduction
2D heat equation.
For steady state analysis the time derivative is zero as the function is varying only with space coordination and not with time ,therefore the time derivative term is zero in the equation .
equation is the most important for solving using iterative techniques.
We can note that the equation is elliptical form hence employing central difference scheme for the problem statement.
We solve for T(i,j)
Final equation for steady state is given above .
Steady state analysis of 2D heat conduction equation .
We use three iterative techniques to solve the above equation which are as follows
- Jacobi method
- Gauss seidel
- Successive over relaxation
These three method vary depending on the value it picks to solve the equation .
We define a boundary of 2D as unit square of one units on x and y and we discretize the domain for 10 equal spaced points and we define the boundary conditions bottom 900k,left 400k,right 800k,top 600k .and we assign average values to the corner points and take a copy of the temperature array as Told .
We have 3 different iterations for three different methods and controlled by conditional statement (IF).assigning 1 for Jacobi ,2 for gauss and 3 for sor method, the main difference is injacobi we use the old values for computing while in gauss method the previously calcuted value for the neighbor cells are used and sor method we assign a value to omega to make a small jump towards the next iteration hence it converges fast compared to other methods . the number of iterations for Jacobi method > gauss siedal >successive order relaxation .
We assign tolerance term to be 1e-4 for convergence when error term falls below this value the loop breaks and plots the final x axis y axis and temperature.this is the procedure for iterative method for solving the 2D heat equation.
***************************************
%defining the domian
nx=10;
ny=10;
l=1;
x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
%initiualizing boundary condition
t=400*ones(nx,ny);
t(:,1)=400;
t(:,end)=800;
t(1,:)=600;
t(end,:)=900;
% initializing corner values
t(1,1)=(600+400)/2;
t(nx,ny)=(900+800)/2;
t(1,end)=(600+800)/2;
t(end,1)=(900+400)/2;
% copying value of temperature
told=t;
%defining loop values
tol=1e-4;
error=9e9
k=(2*(dx^2+dy^2))/(dx^2*dy^2);
w=1.25;%omega value
iteration = 3;%1-jacobi,2for gauss,3 for sor method
%programing
if iteration==1%conditional loop
iter_count=1;
while error > tol%convergence loop
for i=2:nx-1 %space loop
for j=2:ny-1 %space loop
term_h=(told(i-1,j)+told(i+1,j))/(dx^2);
term_v=(told(i,j-1)+told(i,j+1))/(dy^2);
t(i,j)=(term_h+term_v)/k;
end
end
error_x=max(abs(told-t));%error
error=max(error_x);
told=t;%updating previous values
colormap(jet)%plotting
contourf(x,y,t)
iter_count=iter_count+1;
title_text=sprintf('iteration count for steadystate jacobi method=%d',iter_count);
xlabel('x axis');
ylabel('y axis');
set(gca,'ydir','reverse')
title(title_text)
colorbar
pause(0.03)
end
end
if iteration==2
iter_count=1;
while error > tol
for i=2:nx-1
for j=2:ny-1
term_h=(t(i-1,j)+told(i+1,j))/(dx^2);
term_v=(t(i,j-1)+told(i,j+1))/(dy^2);
t(i,j)=(term_h+term_v)/k;
end
end
error_x=max(abs(told-t));
error=max(error_x);
told=t;
colormap(jet)
contourf(x,y,t)
iter_count=iter_count+1;
title_text=sprintf('iteration count staedy state for gauss seidal=%d',iter_count);
title(title_text)
xlabel('x axis');
ylabel('y axis');
set(gca,'ydir','reverse')
colorbar
pause(0.03)
end
end
if iteration==3
iter_count=1;
while error > tol
for i=2:nx-1
for j=2:ny-1
term_h=(t(i-1,j)+told(i+1,j))/(dx^2);
term_v=(t(i,j-1)+told(i,j+1))/(dy^2);
t(i,j)=(1-w)*told(i,j)+w*(term_h+term_v)/k;
end
end
error_x=max(abs(told-t));
error=max(error_x);
told=t;
colormap(jet)
contourf(x,y,t)
iter_count=iter_count+1;
title_text=sprintf('iteration count for staedy state sor method=%d',iter_count);
title(title_text)
xlabel('x axis');
ylabel('y axis');
set(gca,'ydir','reverse')
colorbar
pause(0.03)
end
end
The code has 3 condition statement for three iterative techniques first we have to define the temperature domain top left right and bottom with the given temperature and have we have assigned a guess value for the inner nodes 400k to all the nodes and the average values between orthogonal walls by taking the average value.
we have a convergence loop after the conditional statement which ensures our value to the assigned node is below the error term defined. once this condition is satisfied the loop breaks.and we plot the graph with x axis and y axis and temperature .
we get an animation of the plot with iteration count to know when the system has reached stable condition for temperature .
RESULTS AND DISCUSSIONS
CASE 1 - JACOBI METHOD
we get the plot for the axis and temperature distribution along the domain. and we have the plots per iteration this method uses the old values for calculating the present cell values hence it takes more iterations to solve the function .
CASE-2 GAUSS SEIDEL METHOD
We can see that number of iteration for GS method is less compared to jacobi method thats because it uses the updated values from the previous cell to compute values of the current cell.hence the convergence is faster than jacobi method.
CASE -3 SOR METHOD
This is an interesting method here we use a term omega which is a parameter that lets the solution converge faster depending on the the value of (w).
w<1(under relaxation)
w>1(over relaxation)
in this program we have value for (W)as 1.25 and we can clearly see that the number of iteration for jacobi is greater than gass seidel and sor method . the sor method has the least value of iteration and converges faster and it pirely depends on omega value .
NUMBER OF ITERATIONS
JACOBI METHOD>GAUSS SEIDEL>SOR METHOD .
TRANSIENT ANALYSIS
where n is the time step.
we have two cases one is explicit and other is implicit where in explicit we calculate for the present time step .
For the same problem statement we are applying transient condition hence we have the temperature and other specifications remaining same .
In case of transient solution since time is a varying parameter we use time loop ,convergence and space loop implicit case as the computation requires an update for Tprevious after every time step as they remain constant for convergence part .
Explicit condition we don’t use the convergence loop as they don’t employ these methods hence we have one time and two space loop we give enough time step to solve for steady state solution to arrive . We calculate for n time step for explicit condition and for implicit we calculate for the next time step hence we have an convergence loop an error term derived when less than the tolerance the iteration breaks and the nest time step starts .
The code is given below
Explicit case we need to make sure that CFL value along x and y combined should less than 0.5
CFL<0.5thats the stability criteria for explicit condition to maintain cfl to that value we use the value of dt and dx according .
alpha=1.3
dt=0.002
CFLx=alpha*dt/(dx^2);
CFLy=alpha*dt/(dy^2);
We have the total CFL to be 0.42. Once we have defined the parameters we call in the loop for time and define time period time period = nt/dt ,(100/0.002)=time steps .in explicit methos CFL number is most important we have to maintain it below unit value or the solution may diverge . we do not use a convergence loop in explicit method as we are playing with time steps and grid points .we employ three FOR loops (one time and two space loop) and update the old values after the end of space loop.
clc
clear all
close all
nx=10;
ny=nx;
x=linspace(0,1,nx);
y=x;
dx=1/(nx-1);
dy=1/(ny-1);
dy=y(2)-y(1);
%initiualizing boundary condition
t=400*ones(nx,ny);
t(:,1)=400;
t(:,end)=800;
t(1,:)=600;
t(end,:)=900;
t(1,1)=(600+400)/2;
t(nx,ny)=(900+800)/2;
t(1,end)=(600+800)/2;
t(end,1)=(900+400)/2;
told=t;
%convergence critiria
tol=1e-4;
error=9e9;
%loop specification
iteration=1;%1for jacobi,2-gs,3-sor
alpha=1.3
dt=0.002
k1=alpha*dt/(dx^2);
k2=alpha*dt/(dy^2);
if iteration==1;
iter_count=1;
for it=1:100
for i=2:nx-1 %space loop
for j=2:ny-1
term1=(told(i,j))*(1-2*k1-2*k2)
term2=(k1*(told(i-1,j)+told(i+1,j)))
term3=k2*(told(i,j-1)+told(i,j+1));
t(i,j)=term1+term2+term3;
end
end
told=t;
%plotting x,yaxis and temperature
colormap(jet)
contourf(x,y,t)
iter_count=iter_count+1;
title_text=sprintf(' transient explicit jacobi method');
xlabel('x axis');
ylabel('y axis');
set(gca,'ydir','reverse')
title(title_text)
colorbar
pause(0.03)
end
end
if iteration==2;
iter_count=1;
for it=1:75
for i=2:nx-1
for j=2:ny-1
term1=(told(i,j))*(1-2*k1-2*k2)
term2=(k1*(t(i-1,j)+told(i+1,j)))
term3=k2*(t(i,j-1)+told(i,j+1))
t(i,j)=term1+term2+term3;
end
end
error_x=max(abs(told-t));
error=max(error_x);
told=t;
colormap(jet)
contourf(x,y,t)
iter_count=iter_count+1;
title_text=sprintf('transient explicit gauss seidel method');
xlabel('x axis');
ylabel('y axis');
set(gca,'ydir','reverse')
title(title_text)
colorbar
pause(0.03)
end
end
if iteration==3;
iter_count=1;
for it =1:75
for i=2:nx-1
for j=2:ny-1
w=1.1
term1=(told(i,j))*(1-2*k1-2*k2)
term2=(k1*(t(i-1,j)+told(i+1,j)))
term3=k2*(t(i,j-1)+told(i,j+1))
t(i,j)=(1-w)*told(i,j) + w*(term1+term2+term3);
end
end
error_x=max(abs(told-t));
error=max(error_x);
told=t;
colormap(jet)
contourf(x,y,t)
iter_count=iter_count+1;
title_text=sprintf(' transient explicit sor method');
xlabel('x axis');
ylabel('y axis');
set(gca,'ydir','reverse')
title(title_text)
colorbar
pause(0.03)
end
endThe above code has the main program for all three techniques controled by a conditional staement IF (iteration) and plots the animation depending on the selection of the method .has the same code structure for all the three .we have given enough time steps for the solution to converge .
RESULTS AND DISCUSSION
CASE 1 JACOBI METHOD
CASE -2 GAUSS SEIDEL METHOD
CASE-3 SOR METHOD
We can see that all three plots are the same and we have converged over time we can see the animation from the program . we can also conclude when the time steps are less the solution does not converge and we do not get the right solution the iteration just halts after the given time steps this is one important parameter to understand .the explicit method is easy to understand and simple way of solving a given problem .
IMPLICIT METHOD OF SOLVING
This method is little complex as it computes the value from the next time step as mentioned in the formulae above we n+1 th time step to evalute temperature at the current node .in implicit method we dont need to have CFL condition we can have larger time steps to solve the problem .we have two space and a time loop .first time loop and then before enterning the we have a convergence loop that determines the error term for each and every time step hence we have more number of iterations .as the count keeps growing for each and every time step .
clear all
close all
clc
%defining space cordinates
nx=10;
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=1/(nx-1);
dy=y(2)-y(1);
tol=1e-4;
dt=10
alpha=1.4;
nt=140;
t=400*ones(nx,ny);
t(:,1)=400;
t(:,end)=800;
t(1,:)=600;
t(end,:)=900;
told=t;
tprev=told;
k1=alpha*dt/(dx^2);
k2=alpha*dt/(dy^2);
iteration=3;
if iteration==1
iter_count=1;
for k=1:nt
error=1e9
while error>tol
for i=2:nx-1
for j=2:ny-1
term1=(1+2*(k1+k2));
term2=(told(i-1,j)+told(i+1,j))*k1;
term3=(told(i,j-1)+told(i,j+1))*k2;
t(i,j)=(tprev(i,j)+term2 +term3)/term1;
end
end
error=max(max(abs(told-t)));
told=t;
iter_count=iter_count+1;
end
tprev=t;
colormap(jet)
contourf(x,y,t)
title_text=sprintf('iteration count for transient implicit jacobi method=%d',iter_count);
xlabel('x axis');
ylabel('y axis');
set(gca,'ydir','reverse')
title(title_text)
colorbar
pause(0.03)
end
end
if iteration==2
iter_count=1;
for k=1:nt
error=1e9
while error>tol
for i=2:nx-1
for j=2:ny-1
term1=(1+2*(k1+k2));
term2=(t(i-1,j)+told(i+1,j))*k1;
term3=(t(i,j-1)+told(i,j+1))*k2;
t(i,j)=(tprev(i,j)+term2 +term3)/term1;
end
end
error=max(max(abs(told-t)));
told=t;
iter_count=iter_count+1;
end
tprev=t;
colormap(jet)
contourf(x,y,t)
title_text=sprintf('iteration count fortransient implicit gauss seidel method=%d',iter_count);
xlabel('x axis');
ylabel('y axis');
set(gca,'ydir','reverse')
title(title_text)
colorbar
pause(0.03)
end
end
if iteration==3
iter_count=1;
for k=1:nt
error=1e9
while error>tol
for i=2:nx-1
for j=2:ny-1
w=1.25
term1=(1+2*(k1+k2));
term2=(t(i-1,j)+told(i+1,j))*k1;
term3=(t(i,j-1)+told(i,j+1))*k2;
t(i,j)=(1-w)*told(i,j) +(w*(tprev(i,j)+term2 +term3))/(term1);
end
end
error=max(max(abs(told-t)));
told=t;
iter_count=iter_count+1;
end
tprev=t;
colormap(jet)
contourf(x,y,t)
title_text=sprintf('iteration count for transient implicitSOR method=%d',iter_count);
xlabel('x axis');
ylabel('y axis');
set(gca,'ydir','reverse')
title(title_text)
colorbar
pause(0.003)
end
end
we have the code for the implicit method of solving the problem we have all three methods and defined by a conditional statement IF.we can see from the results that the solution has converged and we have obtained an exact solution for all three case with the iteration count as expected with jacobi having more iteration .
CASE -1 JACOBI METHOD
CASE 2 GAUSS SEIDEL METHOD
CASE 3 SOR METHOD
RESULTS AND DISCUSSIONS .
it is important to note that with increase in time steps the number of iteration incraeses we can play arount with the term 'nt'in the problem to see the number of iteration it takes to solve for the entire time loop and the time required for solving also can be printed in the same plot . We can see that all three cases lead to the same results but the way to go ahead in a problem statement is more important and undewrstanding what parameters play a role in converging a solution and stability parameters like CFL number and time steps are terms that can lead to divergence and we may end up in a completely wrong solution . therefore first we should clearly have knowledge in the steps that the problem statement defines . we can use SOR method for analysis and the main challenge is in selecting the right value for omwga the parameter that makes the solution converge faster .