Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Aim: To solve the 2 Dimensional heat conduction equation for steady and unsteady flow by using iterative methods like jacobi,gauss sidel and succesive over relaxation method by using mat lab.For unstaedy flow we are using both implicit and explicit methods.
Introduction and mathemetical equations:
Consider general three dimension heat conduction equation with out internal heat generation if we are writing for cartesian coordinates,its given by
`(delT)/(delt)=alpha((del^(2)T)/(delx^(2))+(del^(2)T)/(dely^(2))+(del^(2)T)/(delz^(2)))`
where T=temperature in Kelvin
alpha=thermal diffusivity in m2/s
t=time in seconds
For two dimension and steady flow the equation becomes
`(del^(2)T)/(delx^(2))+(del^(2)T)/(dely^(2))=0`
if we discretize above equation by cenral difference
`(T_(i-1,j)-2T_(i,j)+T_(i+1,j))/(trianglex)^(2)`+`(T_(i,j+1)-2T_(i,j)+T_(i,j-1))/(triangley^(2))=0`
If we solve the above equation by using iterative methods
1.Jacobi method
2.Gauss-siedel method
3.Successive over relaxation method
Jacobi method:In this iteration technique we assume the initial values at all nodes to calculate the new values which updated once we find the all the nodes values and this process will goes on until we will arrive at an convergence solution.
Gauss siedel method:In this iteration technique we assume the values obtained in jacobi iteration for next equation immediately with out waiting for completion of first cycle.
Succesive over relaxation method:This method is an updated version of above two methods.In this we use relaxation factor to ease the process
Given data in problem
Top temperature:600K
Bottom temperature:900K
Left side temperature:400K
Right side temperature:800K
Assumptions:
1.Domain is unit square
2.number of grid points along x direction is equal to number of grid points along y direction,nx=ny
Assumed values for given problem
Tol(Tolerance)=1e-3
Omega(Relaxation factor)=1.2
Thermal diffusivity(alpha)=1.4
Error=100
MATLAB CODE FOR STEADY TWO DIMENSIONAL CONDUCTION
clear all
close all
clc
%solving steady state two dimension conduction equation
%input parameters
nx=10;
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=1/(nx-1);
dy=1/(ny-1);
%initializing temp in x and y directions
T=300*ones(nx,ny);
%Assigning boundary conditions
T(1,:)=600;
T(nx,:)=900;
T(:,1)=400;
T(:,ny)=800;
omega=1.4;
%corner temperatures
T(1,1)=(600+400)/2;
T(1,nx)=(600+800)/2;
T(ny,1)=(900+400)/2;
T(nx,ny)=(800+900)/2;
%Assigning error and tolerance
k=2*(dx^(2)+dy^(2))/(dx^(2)*dy^(2));
error=100;
tol=1e-3;
iterative_solver=3;
% calculation of steady two dim equation by jacobi method
Told=T;
if iterative_solver==1
jacobi_iter=1;
while(error>tol)
for i=2:nx-1
for j=2:ny-1
term1=(Told(i-1,j)+Told(i+1,j))/((dx)^(2));
term2=(Told(i,j-1)+Told(i,j+1))/((dy)^(2));
T(i,j)= (term1+term2)/k;
end
end
error=max(max(abs(Told-T)));
Told =T;
jacobi_iter=jacobi_iter+1;
end
colormap(jet)
A=contourf(x,y,T);
clabel(A)
set(gca,'YDir','reverse')
xlabel('X-AXIS')
ylabel('Y-AXIS')
title(sprintf('No of jacobi iterations=%d',jacobi_iter))
end
%Gauss siedel mrthod
if iterative_solver==2
gs_iter=1;
while(error>tol)
for i=2:nx-1
for j=2:ny-1
term1=(T(i-1,j)+Told(i+1,j))/((dx)^(2));
term2=(T(i,j-1)+Told(i,j+1))/((dy)^(2));
T(i,j)= (term1+term2)/k;
end
end
error=max(max(abs(Told-T)));
Told =T;
gs_iter=gs_iter+1;
end
colormap(jet)
A=contourf(x,y,T);
clabel(A)
set(gca,'YDir','reverse')
xlabel('X-AXIS')
ylabel('Y-AXIS')
title(sprintf('No of gauss_siedel_ iterations=%d',gs_iter))
end
% SOR method
if iterative_solver==3
sor_iter=1;
while(error>tol)
for i=2:nx-1
for j=2:ny-1
term1=(T(i-1,j)+Told(i+1,j))/((dx)^(2));
term2=(T(i,j-1)+Told(i,j+1))/((dy)^(2));
T(i,j)= (1-omega)*Told(i,j)+(omega*(term1+term2))/k;
end
end
error=max(max(abs(Told-T)));
Told =T;
sor_iter=sor_iter+1;
end
colormap(jet)
A=contourf(x,y,T);
clabel(A)
set(gca,'YDir','reverse')
xlabel('X-AXIS')
ylabel('Y-AXIS')
title(sprintf('No of SOR iterations=%d',sor_iter))
endResults for steady state two dimension conduction equation
JACOBI METHOD
GAUSS-SIEDEL MRTHOD
SUCCESSIVE OVER RELAXATION METHOD
CONCLUSION
Jacobi method converges after 171 iterations,gauss-siedel method converges after 94 iterations where as SOR method converges after 40 iterations.Based on these results SOR method converges quickly compared to jacobi and gauss-sidel methods.It is because of using relaxation factor omega in successive over relaxation.
.......................................................................................................................................
Now we will solve the unsteady(transient) state two dimensional conduction equation by explicit and implicit method.For explicit method we solve directly and for implicit method we will use again iterative schemes jacobi,gauss-siedel,successive over relaxation methods.
`(delT)/(delt)=alpha((del^(2)T)/(delx^(2))+(del^2T)/(dely^(2)))`
time derivative discretized by forward differencing scheme,spatial terms is discretized by central differencing scheme we will get final equation
`T_(i,j)^(n+1)=T_(i,j)^(n)+k1(T_(i-1,j)-2T_(i,j)+T_(i+1,j))``+k2(T_(i,j+1)-2T_(i,j)+T_(i,j-1))`
where `k1=alpha(triangleT)/(trianglex)^(2)`
`k2=alpha(triangleT)/(triangley)^(2)`
MATLAB PROGRAMME FOR TRANSIENT EQUATION USING EXPLICIT METHOD
clear all
close all
clc
%solving steady state two dimension conduction equation
%input parameters
nx=10;
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=1/(nx-1);
dy=1/(ny-1);
alpha=1.4;
dt=0.001;
%initializing temp in x and y directions
T=300*ones(nx,ny);
%Assigning boundary conditions
T(1,:)=600;
T(nx,:)=900;
T(:,1)=400;
T(:,ny)=800;
%corner temperatures
T(1,1)=(600+400)/2;
T(1,nx)=(600+800)/2;
T(ny,1)=(900+400)/2;
T(nx,ny)=(800+900)/2;
%Assigning error and tolerance
k1=alpha*(dt/dx^(2));
k2=alpha*(dt/dy^(2));
cfl=k1+k2;
% calculation of steady two dim equation by jacobi method
Told=T;
nt=1400;
for w=1:nt
for i=2:nx-1
for j=2:ny-1
term1=Told(i-1,j)-2*Told(i,j)+Told(i+1,j);
term2=Told(i,j+1)-2*Told(i,j)+Told(i,j-1);
T(i,j)= Told(i,j)+k1*(term1)+k2*(term2);
end
end
Told =T;
colormap(jet)
A=contourf(x,y,T);
clabel(A)
set(gca,'YDir','reverse')
xlabel('X-AXIS')
ylabel('Y-AXIS')
title(sprintf('EXPLICIT METHOD'))
end
Results for explicit method
........................................................................................................................................
Now we will solve by transient 2D heat conduction equation by implicit method
time derivative discretized by forward differencing scheme,spatial terms is discretized by central differencing scheme we will get final equation
`T_(i,j)^(n+1)=1/(1+2k1+2k2)(T_(i,j)^(n)+k1(T_(i-1,j)+T_(i+1,j))^(n+1)``+k2(T_(i,j+1)+T_(i,j-1))^(n+1))`
MATLAB CODE FOR IMPLICIT METHOD
%clear all
close all
clc
%solving steady state two dimension conduction equation
%input parameters
nx=10;
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=1/(nx-1);
dy=1/(ny-1);
alpha=1.4;
dt=0.001;
%initializing temp in x and y directions
T=300*ones(nx,ny);
%Assigning boundary conditions
T(1,:)=600;
T(nx,:)=900;
T(:,1)=400;
T(:,ny)=800;
omega=1.2;
%corner temperatures
T(1,1)=(600+400)/2;
T(1,nx)=(600+800)/2;
T(ny,1)=(900+400)/2;
T(nx,ny)=(800+900)/2;
%Assigning error and tolerance
k1=alpha*(dt/dx^(2));
k2=alpha*(dt/dy^(2));
tol=1e-3;
iterative_solver=3;
% calculation of steady two dim equation by jacobi method
Told=T;
T_prev_dt=T;
nt=1400;
if iterative_solver==1
jacobi_iter=1;
for w=1:nt
error=10e4;
while(error>tol)
for i=2:nx-1
for j=2:ny-1
term1=(1+(2*k1)+(2*k2))^-1;
term2=k1*term1;
term3=k2*term1;
H=(Told(i-1,j)+Told(i+1,j));
V=(Told(i,j-1)+Told(i,j+1));
T(i,j)= (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
end
end
error=max(max(abs(Told-T)));
Told =T;
jacobi_iter=jacobi_iter+1;
end
T_prev_dt=T;
end
colormap(jet)
A=contourf(x,y,T);
clabel(A)
set(gca,'YDir','reverse')
xlabel('X-AXIS')
ylabel('Y-AXIS')
title(sprintf('No of jacobi iterations=%d',jacobi_iter))
end
%Gauss siedel mrthod
if iterative_solver==2
gs_iter=1;
for w=1:nt
error=10e4;
while(error>tol)
for i=2:nx-1
for j=2:ny-1
term1=(1+(2*k1)+(2*k2))^-1;
term2=k1*term1;
term3=k2*term1;
H=(T(i-1,j)+Told(i+1,j));
V=(T(i,j-1)+Told(i,j+1));
T(i,j)= (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
end
end
error=max(max(abs(Told-T)));
Told =T;
gs_iter=gs_iter+1;
end
T_prev_dt=T;
end
colormap(jet)
A=contourf(x,y,T);
clabel(A)
set(gca,'YDir','reverse')
xlabel('X-AXIS')
ylabel('Y-AXIS')
title(sprintf('No of gauss_siedel_ iterations=%d',gs_iter))
end
% SOR method
if iterative_solver==3
sor_iter=1;
for w=1:nt
error=100;
while(error>tol)
for i=2:nx-1
for j=2:ny-1
term1=(1+(2*k1)+(2*k2))^(-1);
term2=k1*term1;
term3=k2*term1;
H=(T(i-1,j)+Told(i+1,j));
V=(T(i,j-1)+Told(i,j+1));
Tgs= (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
T(i,j)= (1-omega)*Told(i,j)+omega*Tgs;
end
end
error=max(max(abs(Told-T)));
Told =T;
sor_iter=sor_iter+1;
end
T_prev_dt=T;
end
colormap(jet)
A=contourf(x,y,T);
clabel(A)
set(gca,'YDir','reverse')
xlabel('X-AXIS')
ylabel('Y-AXIS')
title(sprintf('No of SOR iterations=%d',sor_iter))
endPLOTS AND RESULTS FOR IMPLICIT METHOD
JACOBI METHOD
GAUSS-SIEDEL METHOD
SUCCESSIVE OVER RELAXATION METHOD
Matlab code steps above all programmes
1.First the number of grid points along x-axis and y-axis nx=10 and ny=10 points defined
2.calculating dx and dy using linspace command
3.The constants like error,tolerance,omega value,temperature condition at boundaries and edges are defined
4.First if loop we are defining,under if loop,we are using while loop(convergence loop) is exectuted in which T(i,j) is calculated under for loop i=2:nx-1 and again for loop by parameter j=2:ny-1,T(i,j) is calculated,finally error value calculated and after iteration number updated and temperatures are updated
5.While loop ends by when error falls behind tolerance value and T(i,j) is plotted using contourf and color map command
6.The above steps are repeated for Jacobi,Gs-sidel,SOR in which T(i,j) values are altered for each methods for steady state 2D conduction.
7. In explicit method iterative solver(convergence loop,while) is not used.Time loop used(nt=1:1400),then tempearture T(i,j) is found after which oldtemperatures and errors are updated.
8.In implicit scheme,under if loop time loop defined,then error and tolerance values are assigned and in this case convergence loop used T(i,j) calculated under for loop i=2:nx-1 and again for loop by parameter j=2:ny-1 the error value is found out and Told,iteration numbers are updated,while loop end by after Tprevious is updated for each time loop.
Results and Conclusion
| S.no | Method | No of iterations |
| 1 | Unsteady (explicit method) | Not applicable |
| 2 | Unsteady(implicit scheme,jacobi) | 2898 |
| 3 | Unsteady(implicit scheme,Gauss sidel method) | 2550 |
| 4 | Unsteady(implicit scheme,successive over relaxation method) | 2443 |
1.From the above results,we observe explicit scheme does not need any iterations because it does not need any guess values and convergence for explicit solution is based on the previous state values which is already computed
2.Implicit method requires guess values and iteration for convergence depends on future values which are not known at the time of convergence
3.from the above results SOR has fastest convergence rate,followed by gauss-siedel and then jacobi.This is because SOR method is optimised form og gauss sidel method,it scales up convergencerate of gauss-siedel by factor relaxation factor omega,in our problem we assumed omega=1.2