Week 5 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Given condition: 2D HEAT CONDUCTION EQUATION
`(delT)/(delt)=alpha((del^2T)/(del^2x)+(del^2T)/(del^2y))`
Boundary condition:
Top Boundary(T_top) = 600 K
Bottom Boundary(T_bottom) = 900 K
Left Boundary(T_left) = 400 K
Right Boundary(T_right) = 800 K
Absolute error criteria (error_tol)= 1e-4
OBJECTIVE:Solve Steady state and Unsteady state 2D heat conduction equation by using the iterative techniques and write code for both implicit and explicit schemes .
EXPLICIT & IMPLICIT METHODS:
The explicit method and implicit method are numerical analysis methods used to solve a time-dependent differential equation.
Explicit method :This method,calculates the system status at a future time from the currently known system status.
In explicit method, iterative methods are not used,because all values are known values .Just substituting known values in the equation will get solutions.
Implicit method: This method calculates the system status at a future time from the system statuses at present and future times.
In implicit method, iterative methods are used,because all unknown values are presented.
Iterative methods:
Iterative method is a mathematical procedure that uses an initial guess to generate a sequence of improving approximate solutions for a class of problems, in which the n-th approximation is derived from the previous ones.
Jacobi method:
Jacobi method is a iterative method at where all updated values are used for the computations only after all the variables have been updated.
Gauss-Seidel method:
Gauss-Seidel method is a iterative method,where the updated variables are used in the computations as
soon as they are updated.
Successive over-Relaxation (SOR):
SOR is a variant of the Gauss-Seidel method for solving a system of linear equations, resulting in faster convergence. It is also used for slowly converging iterative process.
If relaxation factor <1 then method is under-relaxed.
relaxation factor >1 then method is over-relaxed.
STEADY STATE CONDITION:
Steady state conduction is the form of conduction that happens when the temperature differences driving the conduction are constant,so that after an equilibrium time,the spatial distribution of temperature in conducting object does not change any furthur.
so,2D steady state heat conduction equation as below
`((del^2T)/(del^2x)+(del^2T)/(del^2y))=0`
Therefore,`(delT)/(delt)=0`
Steady state conditon is time independent differential equation.Hence both implicit and explicit methods are not applicable.
CONCEPT:
Use TAYLOR TABLE method to derive required approximation linear equation to given differential equation.
Update all initial and Boundary conditions
Create initial sample array of "T " based on boundary conditions.Assume remaining unknown values .
Create copy array "T" as "Told" to use it in iterative methods
Use Iterative methods to solve unknown values
- While loop used to control "error criteria".
- In while loop,2 "for loops" are used to calculate "T" values along grid points.Hence it is 2D problem, 2 indicies namely "i" & "j" was used along "x" axis and "y" axis respectively.
- Based on Iteration method, "T" equation in for loop will vary
- After "for loop", error value is calculated.
- Update "Told "as" T"
- Update "number of iteration value".
- Finally check error criteria
- while loop continues until error criteria was fails.Once error criteria fails,while loop exit.
DERIVATION OF APPROXIMATE LINEAR EQUATION TO DIFFERENTIAL EQUATION:
`(del^2T)/(del^2x)=(T_(i+1,j)-2*T_(i,j)+T_(i-1,j))/(dx^2)`
`(del^2T)/(del^2y)=(T_(i,j+1)-2*T_(i,j)+T_(i,j-1))/(dy^2)`
`(del^2T)/(del^2x)+(del^2T)/(del^2y)=0`
`(T_(i+1,j)-2*T_(i,j)+T_(i-1,j))/(dx^2)+(T_(i,j+1)-2*T_(i,j)+T_(i,j-1))/(dy^2)=0`
`2*T_(i,j)/(dx^2)+2*T_(i,j)/(dy^2)=(T_(i+1,j)+T_(i-1,j))/(dx^2)+(T_(i,j+1)+T_(i,j-1))/(dy^2)`
`T_(i,j)*(2/(dx^2)+2/(dy^2))=(T_(i+1,j)+T_(i-1,j))/(dx^2)+(T_(i,j+1)+T_(i,j-1))/(dy^2)`
After solving,finally we get
`T_(i,j)=c_1*(T_(i+1,j)+T_(i-1,j))+c_2*(T_(i,j+1)+T_(i,j-1))`
where
`c_1=dy^2/(2*(dx^2+dy^2))`
`c_2=dx^2/(2*(dx^2+dy^2))`
MAIN PROGRAM:
[In OCTAVE,indices are reversed,because octave make calculations based on rows and columns.Hence here rows are "j" index and columns are "i" index]
%% Solving steady state 2D heat conduction
%d^2T/dx^2+d^2T/dy^2=0
clear all
close all
clc
%%Initial condition for x axis(coloumn as i)
%start point of 'x'=sx
sx=0;
%Length of 'x'=Lx
Lx=1;
%number of grid points in 'x axis'=nx
nx=10;
% array of 'x'
x=linspace(sx,Lx,nx);
%%Initial condition for y axis(row as j)
%start point of 'y'=sy
sy=0;
%Length of 'y'=Ly
Ly=1;
%number of grid points in 'y axis'=ny
ny=10;
% array of 'y'
y=linspace(sy,Ly,ny);
%%Input
T_left=400;
T_right=800;
T_top=600;
T_bottom=900;
%%calculation of 'dx'&'dy' values
dx=Lx/(nx-1);
dy=Ly/(ny-1);
%%Initial array
T=273*ones(ny,nx);
%%Apply BC
T(:,1)=T_left;
T(:,end)=T_right;
T(1,:)=T_top;
T(end,:)=T_bottom;
T(1,1)=(T_left+T_top)/2;
T(1,end)=(T_top+T_right)/2;
T(end,1)=(T_left+T_bottom)/2;
T(end,end)=(T_bottom+T_right)/2;
%%Copy of array
Told=T;
Tinitial=T;
%tolerance values
error_tol=1e-4;
%%constant
c1=dy^2/(2*(dx^2+dy^2));
c2=dx^2/(2*(dx^2+dy^2));
%initial iteration
jac_iteration=0;
gs_iteration=0;
sor_iteration=0;
%%method loop
%no.of method = nm;
nm=3;
for m=1:nm
method=m;
if (method==1)
%% JACOBI METHOD
error_jac=9e9;
while(error_jac>error_tol)
% loop to calculate T(j,i)
for j = 2:ny-1
for i= 2:nx-1
T(j,i)=c1*(Told(j,i+1)+Told(j,i-1))+c2*(Told(j+1,i)+Told(j-1,i));
end
end
error_jac = max(max(abs(T-Told)));
Told=T;
jac_iteration=jac_iteration+1;
figure 1
[c,h]=contourf(x,y,T);
clabel(c,h)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
title_text=sprintf('jacobi_iteration=%d',jac_iteration);
title(title_text)
pause(0.01)
end
end
figure 2
subplot(2,1,1)
surf(x,y,T)
title_text=sprintf('jacobi_iteration=%d',jac_iteration);
title(title_text)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
colorbar;
subplot(2,1,2)
pcolor(x,y,T)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
if(method==2)
%% GAUSS SEDAL METHOD
T=Tinitial;
Told=Tinitial;
error_gs=9e9;
while(error_gs>error_tol)
% loop to calculate T(j,i)
for j = 2:ny-1
for i= 2:nx-1
T(j,i)=c1*(T(j,i+1)+T(j,i-1))+c2*(T(j+1,i)+T(j-1,i));
end
end
error_gs = max(max(abs(T-Told)));
Told=T;
gs_iteration=gs_iteration+1;
figure 3
[c,h]=contourf(x,y,T);
clabel(c,h)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
title_text=sprintf('gauss sedial_iteration=%d',gs_iteration);
title(title_text)
pause(0.01)
end
end
figure 4
subplot(2,1,1)
surf(x,y,T)
title_text=sprintf('gauss sedial_iteration=%d',gs_iteration);
title(title_text)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
colorbar;
subplot(2,1,2)
pcolor(x,y,T)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
if(method==3)
%% SOR METHOD
T=Tinitial;
Told=Tinitial;
error_sor=9e9;
w=1;
while(error_sor>error_tol)
% loop to calculate T(j,i)
for j = 2:ny-1
for i= 2:nx-1
T(j,i)=[(1-w)*Told(j,i)]+[w*(c1*(T(j,i+1)+T(j,i-1))+c2*(T(j+1,i)+T(j-1,i)))];
end
end
error_sor = max(max(abs(T-Told)));
Told=T;
sor_iteration=sor_iteration+1;
figure 5
[c,h]=contourf(x,y,T);
clabel(c,h)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
title_text=sprintf('SOR_iteration=%d',sor_iteration);
title(title_text)
pause(0.01)
end
end
figure 6
subplot(2,1,1)
surf(x,y,T)
title_text=sprintf('SOR_iteration=%d',sor_iteration);
title(title_text)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
colorbar;
subplot(2,1,2)
pcolor(x,y,T)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
end
error_jac
jac_iteration
error_gs
gs_iteration
error_sor
sor_iteration
RESULTS:
OBSERVATION FROM RESULT:
- Gauss sedial method has high convergance rate while compare to Jacobi method
- Gauss sedial method and SOR method has same convergence rate because we assume relaxation factor w=1.
UNSTEADY OR TRANSIENT STATE CONDITION:
Unsteady state conduction is the form of conduction that happens when the temperature differences driving the conduction are varying wrt to time.But after an particular time it was expected to attain steady state.Practically where problem, initially unsteady state only.
so,2D steady Unstate heat conduction equation as below
`(delT)/(delt)=alpha((del^2T)/(del^2x)+(del^2T)/(del^2y))`
Unsteady state conditon is time dependent differential equation.Hence both implicit and explicit methods are applicable.
DERIVATION OF APPROXIMATE LINEAR EQUATION TO DIFFERENTIAL EQUATION:
`(del^2T)/(del^2x)=(T_(i+1,j)-2*T_(i,j)+T_(i-1,j))/(dx^2)`
`(del^2T)/(del^2y)=(T_(i,j+1)-2*T_(i,j)+T_(i,j-1))/(dy^2)`
`(delT)/(delt)=(T_(i,j)^(n+1)-T_(i,j)^(n))/(dt)`
`(delT)/(delt)=alpha((del^2T)/(del^2x)+(del^2T)/(del^2y))`
`(T_(i,j)^(n+1)-T_(i,j)^(n))/(dt)=alpha((T_(i+1,j)-2*T_(i,j)+T_(i-1,j))/(dx^2)+(T_(i,j+1)-2*T_(i,j)+T_(i,j-1))/(dy^2))^(n or n+1)`
IMPLICIT METHOD:
This method calculates the system status at a future time from the system statuses at present and future times.
CONCEPT:
Use TAYLOR TABLE method to derive required approximation linear equation to given differential equation.
Update all initial and Boundary conditions
Create initial sample array of "T " based on boundary conditions.Assume remaining unknown values .
Create copy array "T" as "Told" & "Tprev"&"Tinitial" to use it in iterative methods
- Told="T" values of previous iteration
- "Tprev" = "T" values of Previous time step
Use Iterative methods to solve unknown values
- "for loop" is used to control Time loop
- While loop used to control "error criteria".
- In while loop,2 "for loops" (space loop)are used to calculate "T" values along grid points.Hence it is 2D problem, 2 indicies namely "i" & "j" was used along "x" axis and "y" axis respectively.
- Based on Iteration method, "T" equation in for loop will vary
- After "for loop", error value is calculated.
- Update "Told "as" T"
- Update "number of iteration value".
- Finally check error criteria
- while loop continues until error criteria was fails.Once error criteria fails,while loop exit.
- After 'while loop', "Tprev" was updated
- Graph plot will update
- Process repeat until final time step.
DERIVATION OF APPROXIMATE LINEAR EQUATION TO DIFFERENTIAL EQUATION:
`(T_(i,j)^(n+1)-T_(i,j)^(n))/(dt)=alpha((T_(i+1,j)-2*T_(i,j)+T_(i-1,j))/(dx^2)+(T_(i,j+1)-2*T_(i,j)+T_(i,j-1))/(dy^2))^(n+1)`
After solving,finally we get
`T_(i,j]^(n+1)=(T_(i,j)^(n)+(c_1*[T_(i+1,j)+T_(i-1,j)]^(n+1))+(c_2*[T_(i,j+1)+T_(i,j-1)]^(n+1)))/(1+(2*c_1)+(2*c_2))`
where
`c_1=(alpha*dt)/(dx^2)`
`c_2=(alpha*dt)/(dy^2)`
MAIN PROGRAM:
%% Solving Transient state 2D heat conduction by implicit method
%Transient state 2D heat conduction equation:dT/dt=alpha*(d^2T/dx^2+d^2T/dy^2)
clear all
close all
clc
%%Initial condition for x axis(coloumn as i)
%start point of 'x'=sx
sx=0;
%Length of 'x'=Lx
Lx=1;
%number of grid points in 'x axis'=nx
nx=15;
% array of 'x'
x=linspace(sx,Lx,nx);
%%Initial condition for y axis(row as j)
%start point of 'y'=sy
sy=0;
%Length of 'y'=Ly
Ly=1;
%number of grid points in 'y axis'=ny
ny=15;
% array of 'y'
y=linspace(sy,Ly,ny);
%%Initial condition for t axis(row as k)
%start point of 't'=st
st=1;
%Length of 't'=Lt
Lt=1800;
%number of grid points in 't axis'=nt
nt=180;
% array of 'y'
t=linspace(st,Lt,nt);
%%Input
T_left=400;
T_right=800;
T_top=600;
T_bottom=900;
%%calculation of 'dx'&'dy' values
dx=Lx/(nx-1);
dy=Ly/(ny-1);
dt=Lt/nt;
%%Initial array
T=273*ones(ny,nx);
%%Apply BC
T(:,1)=T_left;
T(:,end)=T_right;
T(1,:)=T_top;
T(end,:)=T_bottom;
T(1,1)=(T_left+T_top)/2;
T(1,end)=(T_top+T_right)/2;
T(end,1)=(T_left+T_bottom)/2;
T(end,end)=(T_bottom+T_right)/2;
%%Copy of array
Tinitial=T;
%tolerance values
error_tol=1e-4;
%CONSTANTS
alpha=1.27e-4;
c1=(alpha*dt)/dx^2;
c2=(alpha*dt)/dy^2;
c=c1+c2;
%Iteration values
jac_iteration=0;
gs_iteration=0;
sor_iteration=0;
%%IMPLICIT METHOD:
%%method loop
%no.of method = nm;
nm=3;
for m=1:nm
method=m;
if (method==1)
%% Jacobi Method
%%array initialization for time loop
Tprev=Tinitial;
Told=Tinitial;
%time loop
for n=1:nt
error_jac=9e9;
while(error_jac>error_tol)
% space loop to calculate T(j,i)
for j=2:ny-1
for i=2:nx-1
T(j,i)=(Tprev(j,i)+c1*[Told(j,i+1)+Told(j,i-1)]+c2*[Told(j+1,i)+Told(j-1,i)])/(1+(2*c1)+(2*c2));
end
end
error_jac = max(max(abs(T-Told)));
Told=T;
jac_iteration=jac_iteration+1;
end
Tprev=T;
figure 1
[c,h]=contourf(x,y,T);
clabel(c,h)
title_text=sprintf('Time loop =%d,Total no.of Jacobi iteration =%d',n,jac_iteration);
title(title_text)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
pause(0.01)
end
jac_iteration_final_loop=jac_iteration
T_jac_final_loop=T;
figure 2
subplot(2,1,1)
surf(x,y,T_jac_final_loop)
title_text=sprintf('Time loop =%d,Total no.of Jacobi iteration =%d',n,jac_iteration);
title(title_text)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
colorbar;
subplot(2,1,2)
pcolor(x,y,T_jac_final_loop)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
end
if (method==2)
%% Gauss Sedial Method
%%array initialization for time loop
Tprev=Tinitial;
Told=Tinitial;
%time loop
for n=1:nt
error_gs=9e9;
while(error_gs>error_tol)
% space loop to calculate T(j,i)
for j=2:ny-1
for i=2:nx-1
T(j,i)=(Tprev(j,i)+c1*[T(j,i+1)+T(j,i-1)]+c2*[T(j+1,i)+T(j-1,i)])/(1+(2*c1)+(2*c2));
end
end
error_gs = max(max(abs(T-Told)));
Told=T;
gs_iteration=gs_iteration+1;
end
Tprev=T;
figure 3
[c,h]=contourf(x,y,T);
clabel(c,h)
title_text=sprintf('Time loop =%d,Total no.of Gauss sedial iteration =%d',n,gs_iteration);
title(title_text)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
pause(0.01)
end
gs_iteration_final_loop=gs_iteration
T_gs_final_loop=T;
figure 4
subplot(2,1,1)
surf(x,y,T_gs_final_loop)
title_text=sprintf('Time loop =%d,Total no.of Gauss sedial iteration =%d',n,gs_iteration);
title(title_text)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
colorbar;
subplot(2,1,2)
pcolor(x,y,T_gs_final_loop)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
end
if (method==3)
%% SOR Method
%%array initialization for time loop
Tprev=Tinitial;
Told=Tinitial;
w=1;
%time loop
for n=1:nt
error_sor=9e9;
while(error_sor>error_tol)
% space loop to calculate T(j,i)
for j=2:ny-1
for i=2:nx-1
T(j,i)=[Told(i,j)*(1-w)]+w*[(Tprev(j,i)+c1*[T(j,i+1)+T(j,i-1)]+c2*[T(j+1,i)+T(j-1,i)])/(1+(2*c1)+(2*c2))];
end
end
error_sor = max(max(abs(T-Told)));
Told=T;
sor_iteration=sor_iteration+1;
end
Tprev=T;
figure 5
[c,h]=contourf(x,y,T);
clabel(c,h)
title_text=sprintf('Time loop =%d,Total no.of SOR iteration =%d',n,sor_iteration);
title(title_text)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
pause(0.01)
end
sor_iteration_final_loop=sor_iteration
T_sor_final_loop=T;
figure 6
subplot(2,1,1)
surf(x,y,T_sor_final_loop)
title_text=sprintf('Time loop =%d,Total no.of SOR iteration =%d',n,sor_iteration);
title(title_text)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
colorbar;
subplot(2,1,2)
pcolor(x,y,T_sor_final_loop)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
end
end
RESULTS:
OBSERVATION FROM RESULT:
- Total no.of iteration is totally depends on dx,dy,dt,`alpha`,c1,c2,c values.
- Gauss sedial method has high convergance rate while compare to Jacobi method
- Gauss sedial method and SOR method has same convergence rate because we assume relaxation factor w=1.
EXPLICIT METHOD:
This method,calculates the system status at a future time from the currently known system status.Hence it has all known values in RHS side,iteration methods are not needed.
DERIVATION OF APPROXIMATE LINEAR EQUATION TO DIFFERENTIAL EQUATION:
`(T_(i,j)^(n+1)-T_(i,j)^(n))/(dt)=alpha((T_(i+1,j)-2*T_(i,j)+T_(i-1,j))/(dx^2)+(T_(i,j+1)-2*T_(i,j)+T_(i,j-1))/(dy^2))^(n)`
After solving,finally we get
`T_(i,j)^(n+1)=T_(i,j)^(n)+(c_1*[T_(i+1,j)-2*T_(i,j)+T_(i-1,j)]^n)+(c_2*[T_(i,j+1)-2*T_(i,j)+T_(i,j-1)]^(n))`
where
`c_1=(alpha*dt)/(dx^2)`
`c_2=(alpha*dt)/(dy^2)`
CONCEPT:
Use TAYLOR TABLE method to derive required approximation linear equation to given differential equation.
Update all initial and Boundary conditions
Create initial sample array of "T " based on boundary conditions.Assume remaining unknown values .
Create copy array "T" as "Tprev"&"Tinitial" to use it in iterative methods.
- "Tprev" = "T" values of Previous time step
Use Iterative methods to solve unknown values
- "for loop" is used to control Time loop
- It 2 "for loops" (space loop)are used to calculate "T" values along grid points.Hence it is 2D problem, 2 indicies namely "i" & "j" was used along "x" axis and "y" axis respectively.
- for loop continues until indicies exist.
- After 'for loop', "Tprev" was updated
- Graph plot will update
- Process repeat until final time step.
MAIN PROGRAM:
%% Solving Transient state 2D heat conduction by explicit method
%Transient state 2D heat conduction equation:dT/dt=alpha*(d^2T/dx^2+d^2T/dy^2)
clear all
close all
clc
%%Initial condition for x axis(coloumn as i)
%start point of 'x'=sx
sx=0;
%Length of 'x'=Lx
Lx=1;
%number of grid points in 'x axis'=nx
nx=15;
% array of 'x'
x=linspace(sx,Lx,nx);
%%Initial condition for y axis(row as j)
%start point of 'y'=sy
sy=0;
%Length of 'y'=Ly
Ly=1;
%number of grid points in 'y axis'=ny
ny=15;
% array of 'y'
y=linspace(sy,Ly,ny);
%%Initial condition for t axis(row as k)
%start point of 't'=st
st=1;
%Length of 't'=Lt
Lt=1800;
%number of time loop =nt
nt=180;
% array of 'y'
t=linspace(st,Lt,nt);
%%Input
T_left=400;
T_right=800;
T_top=600;
T_bottom=900;
%%calculation of 'dx'&'dy' values
dx=Lx/(nx-1);
dy=Ly/(ny-1);
dt=Lt/nt;
%%Initial array
T=273*ones(ny,nx);
%%Apply BC
T(:,1)=T_left;
T(:,end)=T_right;
T(1,:)=T_top;
T(end,:)=T_bottom;
T(1,1)=(T_left+T_top)/2;
T(1,end)=(T_top+T_right)/2;
T(end,1)=(T_left+T_bottom)/2;
T(end,end)=(T_bottom+T_right)/2;
%%Copy of array
Tinitial=T;
Tprev=T;
%tolerance values
error_tol=1e-4;
%CONSTANTS
alpha=1.27e-4;
c1=(alpha*dt)/dx^2;
c2=(alpha*dt)/dy^2;
c=c1+c2;
time=0;
%%EXPLICIT METHOD:
%time loop
for n=1:nt
% space loop to calculate T(j,i)
for j=2:nx-1
for i=2:nx-1
T(j,i)=Tprev(j,i)+c1*[Tprev(j,i+1)-2*Tprev(j,i)+Tprev(j,i-1)]+c2*[Tprev(j+1,i)-2*Tprev(j,i)+Tprev(j-1,i)];
end
end
Tprev=T;
time=time+1;
figure 1
[c,h]=contourf(x,y,T);
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
clabel(c,h)
title_text=sprintf('Time loop=%d',time);
title(title_text)
pause(0.01)
end
figure 2
subplot(2,1,1)
surf(x,y,T)
title_text=sprintf('Time loop=%d',time);
title(title_text)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
colorbar;
subplot(2,1,2)
pcolor(x,y,T)
xlabel('x axis')
ylabel('y axis')
zlabel('Temperature')
RESULT:
OBSERVATION FROM RESULT:
- Iteration methods not required
- Results from both implicit and explict methods are almost same.