Mid term project - Solving the steady and unsteady 2D heat conduction problem
This mid-term project was a part of the course - Introduction to CFD using MATLAB and OpenFOAM
STEADY STATE:
Discretization of the equation-
For Successive Over-relaxation method, the following formula is used-
Here, h^2,S(i,j) is the source term. But, as we are not considering the source term in our problem, this term is neglected.
%Steady State 2D conduction
clear all
close all
clc
for a=1:3
n=11; %number of cells
p=linspace(0,1,n); %Creating grid points in X-direction
q=linspace(0,1,n); %Creating grid points in Y-direction
[X Y]=meshgrid(p,q); %Creating a 2D grid
dx=1/(n-1); %Step size in X and Y directions
dy=dx;
beta=1.5; %Overrelaxation factor
%Initializing tolerance and error
tolerance=1e-4;
error=9e9;
%Setting initial BCs
T=ones(n,n);
T(1,:)=900; %Temperature at Bottom boundary
T(end,:)=600; %Temperature at Top boundary
T(:,1)=400; %Temperature at Left boundary
T(:,end)=800; %Temperature at Right boundary
Told=T; %Setting Told=T will later be used to calculate T at higher order points
%a=1;
iter=1;
%Jacobi Method
while error>tolerance
for k=2:n-1 %X-Direction
for l=2:n-1 %Y-Direction
if a==1 %Jacobi
T(k,l)=(Told(k-1,l)+Told(k+1,l)+Told(k,l-1)+Told(k,l+1))/4;
elseif a==2 %Gauss Seidal
T(k,l)=(T(k-1,l)+T(k+1,l)+T(k,l-1)+T(k,l+1))/4;
elseif a==3 %Successive over-relaxation
T(k,l)=((T(k-1,l)+T(k+1,l)+T(k,l+1)+T(k,l-1))*(beta/4))+((1-beta)*Told(k,l));
end
end
end
error1=max(abs(T-Told)); %Calculating error after each time step
error=max(error1); %Calculating maximum error
Told=T;
iter=iter+1;
%Plotting graphs
figure(a)
[C,h]=contourf(X,Y,T);
set(gca,'XTick',(0:dx:1));
set(gca,'YTick',(0:dy:1));
xlabel('X-Axis');
ylabel('Y-Axis');
clabel(C,h,'Fontsize',9);
title_text=sprintf('Iteration Number = %d',iter);
title(title_text);
pause(0.03);
end
%a=a+1;
endRESULTS:
1. Jacobi Method:
2. Gauss-Seidal Method:
3. Successive Over-relaxation Method:
TRANSIENT STATE:
Discretization:
A) Explicit Method:
%Transient state 2D Heat conduction equation using explicit method
clear all
close all
clc
%Initializing the grid
L=1;
n=21;
%Creating the grid
nx=linspace(0,L,n);
ny=linspace(0,L,n);
[X Y]=meshgrid(nx,ny);
dx=L/(n-1);
dy=dx; %Since grid is unoform in X and Y directions
%Initializing domain conditions
alpha_x=1e-2; %Small value taken for simulation purposes
alpha_y=1e-2; %With alpha=1e-4, unable to get solution even after 20,000 iterations
%beta=1.5; %Succesive overrelaxation factor
t=800;
dt=0.01; %dt taken small for simulation purposes
T=298*ones(n,n);
T(1,:)=900; %Temperature at Bottom boundary
T(end,:)=600; %Temperature at Top boundary
T(:,1)=400; %Temperature at Left boundary
T(:,end)=800; %Temperature at Right boundary
Told=T;
K1=(alpha_x*dt)/(dx^2); %For simplification
K2=(alpha_y*dt)/(dy^2);
%Initializing tolerance, error and iteration
tolerance=1e-2; %tolerance taken small to get solution early. For 1e-4, unable to get solution even after 20,000 iterations
error=9e9;
iter=1;
while error>tolerance
%for i=1:t
%X Direction Loop
for j=2:n-1
for k=2:n-1 %Y Direction Loop
Tx=K1*(Told(j-1,k)-(2*Told(j,k))+Told(j+1,k));
Ty=K2*(Told(j,k-1)-(2*Told(j,k))+Told(j,k+1));
T(j,k)=Told(j,k)+Tx+Ty;
%T(i,j)=Told(i,j)+(k1*(Told(i-1,j)-(2*Told(i,j))+Told(i+1,j)))+(k2*(Told(i,j-1)-(2*Told(i,j))+Told(i,j+1)));
end
end
error1=max(abs(T-Told)); %Calculating errors after each iteration
error=max(error1); %Clculating the maximum error
Told=T;
iter=iter+1;
%Plotting the grid
figure(1)
[C,h]=contourf(X,Y,T);
%set(gca,'XTick',(0:dx:1));
%set(gca,'YTick',(0:dy:1));
xlabel('X-Axis');
ylabel('Y-Axis');
clabel(C,h,'Fontsize',9);
title_text=sprintf('Iteration Number = %d',iter);
title(title_text);
pause(0.003);
%end
endRESULT:
B) Implicit Method:
The 2D Heat conduction equation is solved implicitly using the Thomas Tridiagonal Matrix method.
The Tridiagonal system is of the form:
Equations in x1,x2,...,xn are found by converting the matrix into a lower triangular matrix with diagonal elements 1.
The values of x1,x2,...,xn are calculated by back-substitution starting from xn.
%2D Head conduction using implicit method
clear all
close all
clc
L=1; %Initializing domain length
n=21;
dt=0.1; %Taken large for simulation purpose
t=1000; %No of iterations
alpha_x=1e-2; %For simulation purpose, value of alpha is taken large
alpha_y=1e-2;
%Creating the grid
x=linspace(0,L,n);
y=linspace(0,L,n);
[X Y]=meshgrid(x,y);
dx=x(2)-x(1);
dy=dx;
Kx=(alpha_x*dt)/(dx^2); %For simplification
Ky=(alpha_y*dt)/(dy^2);
%Initial temperatures
T=298*ones(n,n);
T(n,1:n-1)=900;
T(1,2:n)=600;
T(1:n-1,1)=400;
T(2:n,n)=800;
%Told=T
k1=Kx/2;
k2=Ky/2;
A=zeros(n*n,n*n); %Initializing TDM
B=zeros(n*n,1); %Initializing RHS matrix
%Initializing the Tridiagonal Matrix
for i=2:n-1;
for j=2:n-1;
A(j+n*(i-1),j+n*(i-1))=1+(2*k1)+(2*k2);
A(j+n*(i-1),j+n*(i-1)-1)=-k1;
A(j+n*(i-1),j+n*(i-1)+1)=-k1;
A(j+n*(i-1),j+n*(i-2))=-k2;
A(j+n*(i-1),j+n*(i))=-k2;
B(j+n*(i-1))=T(i,j);
end
end
%Setting boundary conditions in the vertical directions and intermediate
%values as 1
for j=2:n
A(j*n,j*n)=1;
B(j*n)=800;
A(n*(j-2)+1,n*(j-2)+1)=1;
B(n*(j-2)+1)=400;
end
%Setting boundary conditions in the horizontal directions and intermediate
%values as 1
for k=2:n
A(k,k)=1;
A((n*n)-k+1,(n*n)-k+1)=1;
B(k)=600;
B((n*n)-k+1)=900;
end
iter=1;
%Time loop
for l=1:t
X1=inv(A)*B; %Solving the TDA matrix by taking the inverse
B=X1; %Setting B=X1 to update values in B
for p=2:n-1
for q=2:n-1
T(p,q)=B(q+n*(p-1)); %Back substitution to get values of Temperature from matrix B
end
end
%Plotting the graph
figure(1)
[C,h]=contourf(Y,X,T);
colorbar;
colormap(jet);
%set(gca,'XTick',(0:dx:1));
%set(gca,'YTick',(0:dy:1));
clabel(C,h);
t_text=sprintf('Iteration Number = %d',iter);
title('2D transient state heat conduction');
title(t_text);
xlabel('X-Axis');
ylabel('Y-Axis');
iter=iter+1;
endRESULT: