Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Transient case:
The 2D Conduction Equation
`(delT)/(delt)=alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))`
is a parabolic 2nd order partial differential equation , we can solve this equation by using explicit and implicit method
First we descritise this equation by using finite difference method .We use forward time and central space scheme
hence our equation takes the form
`((T^(n+1))_(i,j)-(T^(n))_(i,j))/(deltat) = alpha(((T^n)_(i-1,j)-2(T^n)_(i,j)+(T^n)_(i+1,j))/(deltax^2)+((T^n)_(i,j-1)-2(T^n)_(i,j)+(T^n)_(i,j+1))/(deltay^2))`
we can re write the same equation as
`(T^(n+1))_(i,j) = (T^(n))_(i,j)+alpha*deltat(((T^n)_(i-1,j)-2(T^n)_(i,j)+(T^n)_(i+1,j))/(deltax^2)+((T^n)_(i,j-1)-2(T^n)_(i,j)+(T^n)_(i,j+1))/(deltay^2))----(1)`
we will be using this equation for solving the explicit scheme until our convergence is reached .
here `T^n` represents the present time step and `T^(n+1)` is n+1 time step, and the subscripts represent the spacial coordinates in x and y direction , for x , "i" index is used and for y "j " index is used
By observing the equation 1 its clear that if boundary conditions were given the problem can be easily solved ,but for the solution to be stable we need to satisfy the stability criteria .
`alpha.(deltat)/(deltax^2),alpha.(deltat)/(deltay^2)<=1/2`
for satisying this our time step should be extremely small .
Our 2D domain has following boundary conditions , length , height are of unit length and its made 9 equal divisions and 10 nodes,
The above solution obtained for t=14, where dt = 0.001 in explicit scheme
code for explicit scheme using octave
clear all
clc
close all
%program for Transient 2D heat conduction equation
dt = 0.001;
%since for satisfying the stability criteria a lower time step is preferred
t = 0:0.001:1.7;
%the time
x =linspace(0,1,10);
dx = x(2)-x(1);
% The x nodes
y =linspace(0,1,10);
dy = y(2)-y(1);
% The y nodes
k =10;
T =zeros(k,k);
T(k,1)=650; % northwest corner of the 2D block
T(k,k) = 850; %southeast corner of the 2D block
T(1,1) = 500; %southwest corner of the 2D block
T(1,k) = 700; %northeast corner of the 2D block
T(2:k-1,1) =400; %left side boundary
T(2:k-1,k) =800; %right side boundary
T(1,2:k-1) =600;% bottom side boundary
T(k,2:k-1) =900; % top side boundary
alpha =1.4;
Tnew=T;
k1=(alpha*dt)/((dx)^2);
k2=(alpha*dt)/((dy)^2);
for t = 1:length(t);
t
for i = 2:k-1;
for j = 2:k-1;
Tnew(i,j) = T(i,j) + (k1*((T(i+1,j)-2*T(i,j)+T(i-1,j))))+ (k2*((T(i,j+1)-2*T(i,j)+T(i,j-1))));
endfor
endfor
err = max(max(abs(Tnew - T)));
drawnow()
% the convergence plot
loglog(t,err,"-*r")
hold on
xlabel("time step")
ylabel("error ")
title("The explicit scheme for 5 seconds ")
T = Tnew;
endfor
figure("2")
contourf(x,y,T,[10],'ShowText','on')
set(gca,'YDir','reverse')
colorbar
title("explicit scheme")
SOLVING THE TRANSIENT HEAT CODUCTION EQUATION IN IMPLICIT SCHEME
Unlike the explicit scheme which is conditionally stable , we can use an implicit scheme which is unconditionally stable ,
`((T^(n+1))_(i,j)-(T^(n))_(i,j))/(deltat) = alpha(((T^(n+1))_(i-1,j)-2(T^(n+1))_(i,j)+(T^(n+1))_(i+1,j))/(deltax^2)+((T^(n+1))_(i,j-1)-2(T^(n+1))_(i,j)+(T^(n+1))_(i,j+1))/(deltay^2))`
since the RHS is on future time step that is n+1 time steps we, need to solve for all the points
This can be solved by arranging equations at each node and forming an tri diagnol system of matrix which can be solved by usimg thomas algorithm , this process is extremely difficult and time consuming
We use a simpler scheme by assuming the unkown values and using an convergence criteria and iterate it till the solution satisifes the converfence criteria
such iterative schemes are
- Jacobian iterative schemes
- Gauss seidel iterative schemes
- Point Succesive over relaxation technique
The result obtained when we use Jacobian iterative scheme after 159 iterations to converge
The result obtained when we use Gauss seidel iterative technique after 85 iterations to converge
The result obtained when we use point successive over relaxtion technique after 496 iterations to converge
in PSOR the omega value taken is 0.3 if we take omega as 1 it will converge in same number of iterations as gauss seidel method
It is observed that PSOR has very fast convergence rate when omega is greater than 1 but it can over shoot the solution if proper omega value is not taken
CODE FOR IMPLICIT METHOD
% program for implicir heat conduction equation using various iterative techniques
clear all
clc
close all
%program for Transient 2D heat conduction equation
dt = 0.1;
%since for satisfying the stability criteria a lower time step is preferred
t = 0:dt:4;
%the time
x =linspace(0,1,10);
dx = x(2)-x(1);
% The x nodes
y =linspace(0,1,10);
dy = y(2)-y(1);
% The y nodes
n =10;
T =zeros(n,n);
T(n,1)=650; % U northwest corner of the 2D block
T(n,n) = 850; %southeast corner of the 2D block
T(1,1) = 500; %Usouthwest corner of the 2D block
T(1,n) = 700; %Usoutheast corner of the 2D block
T(2:n-1,1) =400; %U left side boundary
T(2:n-1,n) =800; %U right side boundary
T(1,2:n-1) =600;% U bottom side boundary
T(n,2:n-1) =900; % top side boundary
alpha =1.4;
Tnew=T;
T_prev =T;
Told =T;
s1=(alpha*dt)/((dx)^2);
s2=(alpha*dt)/((dy)^2);
k = (1+2*s1+2*s2)^-1;
k1 = s1*k;
k2 =s2*k;
error = 9e9;
tol =1e-4;
[X,Y]=meshgrid(x,y);
method =input("enter a interative techniqu number ");
iter =0;
w=0.3;%%omega in succesive over relaxation
% for Jacobi Iterative techniques
if (method == 1)
for t = 1:length(t);
while (error > tol)
for i = 2:n-1;
for j = 2:n-1;
Tnew(i,j) = k*T_prev(i,j) + (k1*((T(i+1,j)+T(i-1,j))))+ (k2*((T(i,j+1)+T(i,j-1))));
endfor
endfor
error = max (max(abs(Tnew-T)))
%drawnow()
loglog(iter,error,"-*-r")
xlabel("no of iterations")
ylabel("error")
title(sprintf("jacobi iterative technique converged at %d iterations",iter))
hold on
T = Tnew;
iter =iter+1
endwhile
T_prev=Tnew;
endfor
endif
% for Gauss seidel Iterative techniques
if (method == 2)
for t = 1:length(t);
while (error > tol)
for i = 2:n-1;
for j = 2:n-1;
Tnew(i,j) = k*T_prev(i,j) + (k1*((Told(i+1,j)+Tnew(i-1,j))))+ (k2*((Told(i,j+1)+Tnew(i,j-1))));
endfor
endfor
error = max(max(abs(Tnew-Told)));
drawnow()
loglog(iter,error,"-*-r")
xlabel("no of iterations")
ylabel("error")
title(sprintf("gauss seidel converged at %d iterations",iter))
hold on
Told=Tnew;
iter =iter+1
endwhile
T_prev=Tnew;
endfor
endif
% for Successive Over Relaxation Iterative techniques
if (method == 3)
for t = 1:length(t);
while (error > tol)
for i = 2:n-1;
for j = 2:n-1;
Tnew(i,j) = ((1-w)*Told(i,j)) + w*((k*T_prev(i,j)) + (k1*((Told(i+1,j)+Told(i-1,j))))+ (k2*((Told(i,j+1)+Told(i,j-1)))));
endfor
endfor
error = max(max(abs(Told-Tnew)));
drawnow()
loglog(iter,error,"-*-r")
xlabel("no of iterations")
ylabel("error")
title(sprintf("point succesive over relaxation converged at %d iterations",iter))
hold on
Told = Tnew;
iter =iter+1
endwhile
T_prev=Tnew;
endfor
endif
figure("2")
contourf(X,Y,Tnew,'ShowText','on')
set(gca,'YDir','reverse')
xlabel("length")
ylabel("height")
colorbar
hold on
if (method==1)
title("point jacobi iterative scheme converged at %d iterations",iter)
elseif (method==2)
title("point gauss seidel iterative scheme converged at %d iterations",iter)
else
title("point succesive over relaxation iterative scheme converged at %d iterations",iter)
endifSTEADY STATE 2D CONDUCTION EQUATION
`(del^2T)/(delx^2)+(del^2T)/(dely^2) =0`
The steady state 2D conduction is of elliptic 2nd order partial differential equation
the discritized form of this equation
`(((T^n)_(i-1,j)-2(T^n)_(i,j)+(T^n)_(i+1,j))/(deltax^2)+((T^n)_(i,j-1)-2(T^n)_(i,j)+(T^n)_(i,j+1))/(deltay^2))=0`
for solving this equation we again need to use an tridiagnol matrix which is a time consuming prodedures .
Hence iterative techniques can be used by assuming the unknowns and using the convergence criteria
the iterative schemes used here are jacobian ,gauss seidel and point successive over relaxation
The result obtained when we use Jacobian iterative scheme after 216 iterations to converge
The result obtained when we use Gauss seidel iterative scheme after 116 iterations to converge
The result obtained when we use psor iterative scheme after 61 iterations to converge
as observed previously psor has higher convergence rate in comparison with other two iterative schemes
The code for the steady state is given by
clear all
close all
clc
% A program for solving the steady state heat conduction equation
x=linspace(0,1,10);
dx = x(2)-x(1);
% The x nodes
y=linspace(0,1,10);
dy = y(2)-y(1);
% The y nodes
k =10;
tol =1e-4;
T =zeros(k,k);
T(k,1)=650; % northwest corner of the 2D block
T(k,k) = 850; %southeast corner of the 2D block
T(1,1) = 500; %southwest corner of the 2D block
T(1,k) = 700; %northeast corner of the 2D block
T(2:k-1,1) =400; %left side boundary
T(2:k-1,k) =800; % right side boundary
T(1,2:k-1) =600; % bottom side boundary
T(k,2:k-1) =900; % top side boundary
Told=T;
Tnew =T;
k1 =(2*((dx^2) + (dy^2))/(dx^2*dy^2))^-1;
err = 9e9;
iteration =0;
w=1.3;
method = input("enter the type of iterative technique: ")
%for jacobian iterative technique
if (method == 1)
while (err > tol)
for i = 2:k-1;
for j = 2:k-1;
T(i,j) =(k1*((Told(i+1,j)+Told(i-1,j))/dx^2))+ (k1*((Told(i,j+1)+Told(i,j-1))/dy^2));
endfor
endfor
err = max(max(abs(Told-T)))
loglog(iteration,err,"-*b")
hold on
title(sprintf("Jacobi Iterative Scheme converged at %d iterations",iteration))
xlabel("Iterations")
ylabel("error")
Told = T;
iteration =iteration +1
endwhile
endif
%%%%%%%%%%%%%%%for Gauss Siedel iterative technique
if (method ==2)
while (err > tol)
for i = 2:k-1;
for j = 2:k-1;
T(i,j) = (k1*((Told(i+1,j)+T(i-1,j))/dx^2)) + (k1*((Told(i,j+1)+T(i,j-1))/dy^2));
endfor
endfor
err = max(max(abs(Told-T)));
loglog(iteration,err,"-*b")
hold on
title(sprintf("Gauss Seidel Iterative Scheme converged at %d iterations",iteration))
xlabel("Iterations")
ylabel("error")
Told = T;
iteration =iteration +1
endwhile
endif
%%%%%%%%%%%%%% for PSOR iterative technique
if (method == 3)
while (err > tol)
for i = 2:k-1;
for j = 2:k-1;
T(i,j) = ((1-w).*T(i,j)) + w.*((k1*((Told(i+1,j)+T(i-1,j))/dx^2)) + (k1*((Told(i,j+1)+T(i,j-1))/dy^2)));
endfor
endfor
err = max(max(abs(Told-T)));
loglog(iteration,err,"-*b")
title(sprintf("PSOR Iterative Scheme converged at %d iterations",iteration))
xlabel("Iterations")
ylabel("error")
hold on
Told = T;
iteration =iteration +1
endwhile
endif
figure("2")
contourf(x,y,T,[10],'ShowText','on')
set(gca,'YDir','reverse')
colorbar
title("")