Steady and Transient state Analysis of Two Dimensional Temperature Heat conduction Equation
Title : Steady state analysis & Transient State Analysis
Objective: 1.To Solve the 2D heat conduction equation by using the point iterative techniques for 1. Jacobi, 2. Gauss-seidel, 3. Successive over-relaxation.
2.code for both implicit and explicit schemes for the Transient Analysis.
3. Code for steady state.
Theory:
`(∂T)/(∂t)= α((∂^2T)/(∂x^2)+(∂^2T)/(∂y^2))`
The above equation is 2D tempreature heat conduction equation, to solve this equation there are two approch.
1. Transient equation.
2. Steady state.
1. Transient equation.
In the transient case there are two methods are use to solve eqaution.
1. Explicit method
`T^(n+1)(i,j)- T^(n)(i,j)-α((∂^2T)/(∂x^2)+(∂^2T)/(∂y^2))=0`
After discritization of grid on x & y-axis we get follwing equation
Use of point technique, LHS use of forword diffrence on RHS use of Central diffrence as there is second order derivative. on RHS use ^n power term. which known by previous time step.
`T^(n+1)(i,j)- T^(n)(i,j)-αΔt((T(i+1,j)-2T(i,j)+T(i-1,j))/(Δx^2)+(T(j+1,i)-2T(i,j)+T(j-1,i))/(Δy^2))^n=0`
simplify the above term.
`T^(n+1)(i,j)=T^(n)(i,j)+αΔt((T_R-2T_P+T_L)/(Δx^2)+(T_T-2T_P+T_B)/(Δy^2))^n`
After taking future (n+1) term of Tempreature on LHS
`T^(n+1)(i,j)= T_P^(n)(i,j)*(1-2K_1+2K_2)+K_1(T_R+T_L)^n+K_2(T_T+T_B)^n`
where, K1 & k2 are the courant friedrichs lewy no. which use for stable the solution for explicit. genrally it uses between 0 to 0.75 are uses for stable the solution.
k=(k1+k2), where k=courant friedrichs lewy no
`K1=(αΔT)/(Δx^2) & K2=(αΔT)/(Δy^2)`
Whrere H & V are the
`H=((Told_(i-1,j) +Told_(i+1,j))`
`H=((Told_(j-1,i) +Told_(j+1,i))`
`T(i,j)= T_Old^(n)(i,j)*(1-2K_1+2K_2)+K_1*(H)+K_2*(V)` ....Equation (1)
The above equation on RHS side there is all known term by previous time step. so, we get all uncoupled linear equation. so, we are not use any iterative solvers technique. so we are solving the equation (1) by time marching method. just simply taken no.of time step. so, above equation (1). is the final eqaution for code in MATLAB.
Implicit Method:
In case of implcit we are take the RHS side ^ (n+1) , we are finding future time step value. so, we get some couple eqaution . this eqaution is very difficult solve by simple matrix. we need to use some iterative solver Technique. following are the tecnique
1. Jacobi solver
2. Gauss seidel Solver
3. Sucessive over Relaxation
After discritization of grid on x & y-axis we get follwing equation
`T^(n+1)(i,j)=T^(n)(i,j)+αΔt((T_R-2T_P+T_L)/(Δx^2)+(T_T-2T_P+T_B)/(Δy^2))^(n+1)`
After rearrange the tem we get it
`T^(n+1)(i,j)= T^n_Prev(i,j)+([K_1(T_R+T_L)^(n+1)+K_2(T_T+T_B)^(n+1))/(1+2K_1+2K_2)`
for use above equation in MATLAB we need to do some simplification
where k1 & k2,
in case of implicit the courant friedrichs lewy no. which use for stable the solution
`K1=(αΔT)/(Δx^2) & K2=(αΔT)/(Δy^2)`
term1= (1+2*K1+2*K2)
term2= k1*term1
term3=k2*term1
`H=((Told_(i-1,j) +Told_(i+1,j))`
`V=((Told_(j-1,i) +Told_(j+1,i))`
So, final eqaution obtained
`T(i,j)= T_prev(i,j)*term1+term2*(H)+term3*(V)``------equation(2)
the above eqaution (2) use for solve the Jacobi solver in matlab code
Gauss seidel Solver
in case of gauss seidal we use previous geneartion calculated value of T(i-1), T(j-1) are uses. this is the diffrence between Gauss seidel Solver & Jacobi solver.
same equation (2) we are use. only change is that
`H=((T_(i-1,j) +Told_(i+1,j))`
`V=((T_(j-1,i) +Told_(j+1,i))`
`T(i,j)= T_prev(i,j)*term1+term2*(H)+term3*(V)`------equation(3)
final above eqaution(3) for MATLAB code.
Sucessive over Relaxation
In this method there use of following equation
`H=((T_(i-1,j) +Told_(i+1,j))`
`V=((T_(j-1,i) +Told_(j+1,i))`
`Tgauss= T_prev(i,j)*term1+term2*(H)+term3*(V)``
`T_(i,j)=Told(i,j)*(1-omega)+omega*(Tgauss)`-----equatiion (4`
where omega >1= overrelaxed
omega<1= under realaxed
above equation (4) uses for MATLAB code
2. Steady state
`α((∂^2T)/(∂x^2)+(∂^2T)/(∂y^2))=0`
where on RHS side there is time derivative term is zero , so we are not use of no.of time step to solve the code.
After discritization of grid on x & y-axis we get follwing equation
`((T_R-2T_P+T_L)/(Δx^2)+(T_T-2T_P+T_B)/(Δy^2))=0`
After rearrange the term,
`(2*(Δx^2+Δy^2)/(Δx^2*Δy^2))=K`
`T_P=1/K((T_R+T_L)/(Δx^2)+(T_T+T_B)/(Δy^2))`
As, dx=dy so,
where H & V
`H=((Told_(i-1,j) +Told_(i+1,j))`
``V=((Told_(j-1,i) +Told_(j+1,i))``
`T_(i,j)=1/4((H)+(V))`---------equation (5)
the above eqaution (5) are uses in MATLAB code for Jacobi solver
Gauss seidel Solver
`H=((T_(i-1,j) +Told_(i+1,j))`
``V=(T_(j-1,i) +Told_(j+1,i))`
`T_(i,j)=1/4((H)+(V))`---------equation (6)`
the above eqaution (6) are uses in MATLAB code for Gauss seidel
Sucessive over Relaxation
`H=((T_(i-1,j) +Told_(i+1,j))`
`V=(T_(j-1,i) +Told_(j+1,i))`
`T_(i,j)=1/4((H)+(V))`
`T_(i,j)=Told(i,j)*(1-omega)+omega*(T(i,j))`-----Equation (7)
the above eqaution (6) are uses in MATLAB code for Sucessive over Relaxation.
CODE For Transient Equation
1. CODE for Implicit method.
Main CODE
Note: For run code of diffrent iterative solver use 1- Jacobi, 2- Guass siedal, 3- SOR
clear all
close all
clc
% TRANSIANT IMPLICIT METHOD
%code for solve a second order 2D Temprature heat conduction equation
% dT/dt=aplha*(d^2T/dx^2 + d^2T/dy^2) =0
%input of given variable
L=1; % One unit square
nx=11; % no.of grid point along x & y direction are same
ny=nx;
nt=400; % no.of time step taken to analysis the temp. distribution.
x=linspace(0,L,nx); % no.of grid point is discretize along x-axis
y=linspace(0,L,ny); % no.of grid point is discretize along y-axis
dx=L/(nx-1); % distance between two grid point along x-axis
dy=L/(ny-1); % distance between two g666666rid point along y-axis
t=4; % simulation time is asummed
dt=0.01;
alpha=1; % therma diffusivity should be between 1 to 2
nt = t/dt; % for determine th no.of time step or dt.
% consider that temprature is one at all row & colomn
T=ones(nx,ny);
% the boundry condition as per given in problem
T_L=400; % Temprature at left side
T_R=800; % Temprature at right side
T_T=600; % Temprature at top side
T_B=900; % Temprature at bottom side
% locate the boundry condition
T(:,1)=T_L; % 1st value in row
T(:,end)=T_R; % last value in row
T(1,:)=T_T; % 1st colomn value
T(end,:)=T_B; % last colomn value
% define the temprature at corner or edges
T(end,1)=0.5*(T_L+T_B); % bottom left
T(end,end)=0.5*(T_R+T_B); % bottom right
T(1,1)=0.5*(T_T+T_L); % Top left
T(1,end)=0.5*(T_T+T_R); %Top right
% define K1 & K2 It\'s called CFL
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
%simplify terms
term1=(1+(2*k1)+(2*k2))^-1;
term2=(k1*term1);
term3=(k2*term1);
omega=1.3; % for SOR Solvers it decide solver methods is under & over relaxation
iterative_solver=1; %(Iterative solver,1=Jacobi, 2=Gauss seidel, 3=SOR)
% error value assume(it shoul be greater than tolearance value)
Told=T;
T_Prev =T;
%enter in to the while loop
error=1e+8;
%tolerance uses to check the error value reduces or not
tol=1e-4;
%jacobi method
if iterative_solver==1
tic; % for start stopwatch
jacobi_iter=1;
for k=1:nt
error=1e8; % start Time loop
while(error > tol) % start convergence loop
for i=2:nx-1 %nodal loop
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)=(T_Prev(i,j)*term1)+(H*term2)+(V*term3);
end
end
error=max(max(abs(Told-T)));
%update old velocity for end the while loop
Told=T;
jacobi_iter=jacobi_iter+1;
end
T_Prev=T;
end
% End of time loop
toc;
F=toc;
[P,Q]=contourf(x,y,T); % creating field contour plot
clabel(P,Q);
colormap(jet)
colorbar
title_1=sprintf(\'iteration number=%d\',jacobi_iter);
title_2=sprintf(\'solving time =%d\',F);
title({\'Transient Implicit Jacobi solver\',title_1,title_2});
xlabel(\'x-axis\')
ylabel(\'y-axis\')
end
%GAUSS SIEDAL METOD
if iterative_solver==2
tic; % for start stopwatch
Gauss_iter=1;
for k=1:nt
error=1e8; % start Time loop
while(error > tol) % start convergence loop
for i=2:nx-1 %nodal loop
for j=2:ny-1
H=(T(i-1,j)+Told(i+1,j));
V=(T(i,j-1)+Told(i,j+1));
T(i,j)=(T_Prev(i,j)*term1)+(H*term2)+(V*term3);
end
end
error=max(max(abs(Told-T)));
%update old velocity for end the while loop
Told=T;
Gauss_iter=Gauss_iter+1;
end
T_Prev=T;
end
% End of time loop
toc;
F=toc;
[P,Q]=contourf(x,y,T); % creating field contour plot
clabel(P,Q);
colormap(jet)
colorbar
title_1=sprintf(\'iteration number=%d\',Gauss_iter);
title_2=sprintf(\'solving time =%d\',F);
title({\'Transient Implicit Gauss solver\',title_1,title_2});
xlabel(\'x-axis\')
ylabel(\'y-axis\')
end
% SOR Iterative solver
if iterative_solver==3
tic; % for start stopwatch
sor_iter=1;
for k=1:nt
error=1e8; % start Time loop
while(error > tol) % start convergence loop
for i=2:nx-1 %nodal loop
for j=2:ny-1
H=(T(i-1,j)+Told(i+1,j));
V=(T(i,j-1)+Told(i,j+1));
Tgs=(T_Prev(i,j)*term1)+(H*term2)+(V*term3);
T(i,j)=(Told(i,j)*(1-omega))+(omega*(Tgs));
end
end
error=max(max(abs(Told-T)));
%update old velocity for end the while loop
Told=T;
sor_iter=sor_iter+1;
end
T_Prev=T;
end
% End of time loop
toc;
F=toc;
[P,Q]=contourf(x,y,T); % creating field contour plot
clabel(P,Q);
colormap(jet)
colorbar
title_1=sprintf(\'iteration number=%d\',sor_iter);
title_2=sprintf(\'solving time =%d\',F);
title({\'Trasient Implicit SOR Solver\',title_1,title_2});
xlabel(\'x-axis\')
ylabel(\'y-axis\')
end
Output Result of Implicit method of Transient Equation
a.) Jacobi Solver
b) Guass Seidal Solver
c) SOR Solver
Grpah conclusion of Implict method
1.from the above graph it\'s found that SOR solver take less no. of iteration=1085 at omega=1.3 & also simulation time 0.014 sec take to solve the eqaution.
2. For omega>1.3 for SOR method take more no. of iteration & take more simultion time.
3. Same for omega <1.3 SOR methhod take more no.of iteration.
4. Guass siedal solver take less no. of iteration=1485 & simulation time=0.015 sec to solve the eqaution comapre with the Jacobi method.
5. For less than 10 no.of time step, solution is start to unstable.
2. CODE for Explicit method.
Main code
close all
clear all
clc
% EXPLICIT METHOD
%input of given variable
L=1; % One unit square
nx=11; % no.of grid point along x & y direction are same
ny=nx;
nt=400; % no.of time step taken to analysis the temp. distribution.
x=linspace(0,L,nx); % no.of grid point is discretize along x-axis
y=linspace(0,L,ny); % no.of grid point is discretize along y-axis
dx=L/(nx-1); % distance between two grid point along x-axis
dy=L/(ny-1); % distance between two g666666rid point along y-axis
t=4; % simulation time is asummed
dt=0.01;
nt=t/dt; % for determine th no.of time step or dt.
alpha=1; % therma diffusivity should be between 1 to 2
% consider that temprature is one at all row & colomn
T=ones(nx,ny);
% the boundry condition as per given in problem
T_L=400; % Temprature at left side
T_R=800; % Temprature at right side
T_T=600; % Temprature at top side
T_B=900; % Temprature at bottom side
% locate the boundry condition
T(:,1)=T_L; % 1st value in row
T(:,end)=T_R; % last value in row
T(1,:)=T_T; % 1st colomn value
T(end,:)=T_B; % last colomn value
% define the temprature at corner or edges
T(end,1)=0.5*(T_L+T_B); % bottom left
T(end,end)=0.5*(T_R+T_B); % bottom right
T(1,1)=0.5*(T_T+T_L); % Top left
T(1,end)=0.5*(T_T+T_R); %Top right
% error value assume(it shoul be greater than tolearance value)
Told=T;
% define K1 & K2 It\'s called CFL
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
% C.F.L No
k=0.2555; % k is C.F.L no it is uses for stable the solution
% if we go beyond the 0.25 the solution become unstalble &
% solution become blow-up
counter=0;
tic;
for p=1:nt
% start Time loop
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)=Told(i,j)*(1-4*k)+k*((H)+(V));
end
end
error=max(max(abs(Told-T)));
Told=T;
counter=counter+1;
end
toc;
F=toc;
%plot graph
figure
[P,Q]=contourf(x,y,T); % creating field contour plot
clabel(P,Q);
colorbar
colormap(jet);
title_1=sprintf(\'iteration number=%d\',counter);
title_2=sprintf(\'solving time =%d\',F);
title_3=sprintf(\'CFL=%d\',k);
title({\'Transient-Explicit\',title_1;title_2,title_3});
xlabel(\'x-axis\');
ylabel(\'y-axis\') ;
Output Result of Explicit method of Transient Equation
Grpah conclusion of Explict method
1.In case of Explict method CFL no . have much imporatnce to stablize the solution. at CFL=2.55 it is found that solution is stable & if we go above or below the same CFL no. solution become unstable.
2. For less than 40 no.of time step, solution is start to unstable.
Overall conclusion of Implcit & Explicit
1. Even though implcit method is very difficult , but it give us a stable solution by diffrent solvers method. SOR method take very less time & n.o of iteration to stable solution & solve eqaution.
2. Implicit method is not that much affected by CFL No. as there is in case of explicit CFL have much importance to stable the solution.
3. No. of time step that much affected to Implicit & Explicit. for less than 40 no.of time step, solution is start to unstable.
3.Steady State code for solve 2D Temprature heat conduction equation.
Code for steady state
Note: For run the code of diffrent iterative solver use 1- Jacobi, 2- Guass siedal, 3- SOR
clear all
close all
clc
% STEADY STATE METHOD
%code for solve a second order 2D Temprature heat conduction equation
% dT/dt=aplha*(d^2T/dx^2 + d^2T/dy^2) =0
%input of given variable
L=1; % One unit square
nx=11; % no.of grid point along x & y direction are same
ny=nx;
nt=400; % no.of time step taken to analysis the temp. distribution.
x=linspace(0,L,nx); % no.of grid point is discretize along x-axis
y=linspace(0,L,ny); % no.of grid point is discretize along y-axis
dx=L/(nx-1); % distance between two grid point along x-axis
dy=L/(ny-1); % distance between two g666666rid point along y-axis
t=4; % simulation time is asummed
dt=0.01;
alpha=1; % therma diffusivity should be between 1 to 2
nt = t/dt; % for determine th no.of time step or dt.
% consider that temprature is one at all row & colomn
T=ones(nx,ny);
% the boundry condition as per given in problem
T_L=400; % Temprature at left side
T_R=800; % Temprature at right side
T_T=600; % Temprature at top side
T_B=900; % Temprature at bottom side
% locate the boundry condition
T(:,1)=T_L; % 1st value in row
T(:,end)=T_R; % last value in row
T(1,:)=T_T; % 1st colomn value
T(end,:)=T_B; % last colomn value
% define the temprature at corner or edges
T(end,1)=0.5*(T_L+T_B); % bottom left
T(end,end)=0.5*(T_R+T_B); % bottom right
T(1,1)=0.5*(T_T+T_L); % Top left
T(1,end)=0.5*(T_T+T_R); %Top right
omega=1.5; %for SOR Solvers it decide solver methods is under & over relaxation
iterative_solver=1; %(Iterative solver,1=Jacobi, 2=Gauss seidel, 3=SOR)
% error value assume(it shoul be greater than tolearance value)
Told=T;
%enter in to the while loop
error=1e+8;
%tolerance uses to check the error value reduces or not
tol=1e-4;
%jacobi method
if iterative_solver==1
tic; % for start stopwatch
jacobi_iter=1;
error=1e8; % start Time loop
while(error > tol) % start convergence loop
for i=2:nx-1 %nodal loop
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)=0.25*(H+V);
end
end
error=max(max(abs(Told-T)));
%update old velocity for end the while loop
Told=T;
jacobi_iter=jacobi_iter+1;
end
toc;
F=toc;
[P,Q]=contourf(x,y,T); % creating field contour plot
clabel(P,Q);
colormap(jet)
colorbar
title_1=sprintf(\'iteration number=%d\',jacobi_iter);
title_2=sprintf(\'solving time =%d\',F);
title({\'Steady state Jacobisolver\',title_1,title_2});
xlabel(\'x-axis\')
ylabel(\'y-axis\')
end
%GAUSS SIEDAL METOD
if iterative_solver==2
tic; % for start stopwatch
Gauss_iter=1;
error=1e8; % start Time loop
while(error > tol) % start convergence loop
for i=2:nx-1 %nodal loop
for j=2:ny-1
H=(T(i-1,j)+Told(i+1,j));
V=(T(i,j-1)+Told(i,j+1));
T(i,j)=0.25*(H+V);
end
end
error=max(max(abs(Told-T)));
%update old velocity for end the while loop
Told=T;
Gauss_iter=Gauss_iter+1;
end
toc;
F=toc;
[P,Q]=contourf(x,y,T); % creating field contour plot
clabel(P,Q);
colormap(jet)
colorbar
title_1=sprintf(\'iteration number=%d\',Gauss_iter);
title_2=sprintf(\'solving time =%d\',F);
title({\'Steady state Gauss solver\',title_1,title_2});
xlabel(\'x-axis\')
ylabel(\'y-axis\')
end
% SOR Iterative solver
if iterative_solver==3
tic; % for start stopwatch
sor_iter=1;
error=1e8; % start Time loop
while(error > tol) % start convergence loop
for i=2:nx-1 %nodal loop
for j=2:ny-1
H=T(i-1,j)+Told(i+1,j);
V=T(i,j-1)+Told(i,j+1);
T(i,j)=0.25*(H+V);
T(i,j)=(Told(i,j)*(1-omega))+omega*T(i,j);
end
end
error=max(max(abs(Told-T)));
%update old velocity for end the while loop
Told=T;
sor_iter=sor_iter+1;
end
toc;
F=toc;
[P,Q]=contourf(x,y,T); % creating field contour plot
clabel(P,Q);
colormap(jet)
colorbar
title_1=sprintf(\'iteration number=%d\',sor_iter);
title_2=sprintf(\'solving time =%d\',F);
title({\'Steady state SOR Solver\',title_1,title_2});
xlabel(\'x-axis\')
ylabel(\'y-axis\')
end
Output Result of Steady State method
A)Jacobi Solver
2.Gauss Seidal Solver
3. Sucessive over Relaxation solver
Graph conclusion of Steady state method
1. SOR method take very less no.of iteration =42 & solving time=0.0047sec comapre with the other two method Jacobi & Gauss Seidal.
2. For SOR method at omega=1.5 gives a less no.of less iteration & simulation time.
3. For beyond & less than the omega=1.5 , SOR solver take more no. of iteartion to solve the equation.
Overall Conclusion of Project
1. steady state method take very less no.of iteration & solving time to solve eqaution & to stable the solution compare with other two method. beacuse steady state method is not depends upon the time derivative term
2. Explicit method is very much depends upon the CFL no. for stable the unstable solution. for more & less than CFL=2.55 found that solution is start to unstable.
3. Eventhough Implcit method is very difficult for SOR method it take very less no. of iteration & simulation time to solve eqaution & to stable the solution.
4. For SOR method omega value much affect on solution to reduce the no.of iteration & simulation time.