Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Please see the attached word document for the output plots and explainination.
clear all
close all
clc
L=1;%Length of the domain
omega=1.25;% SOR factor
nx=100;%grid size
ny=nx;
x=linspace(0,L,nx);
dx=x(2)-x(1);
y=linspace(0,L,ny);
dy=y(2)-y(1);
T=300*ones(nx,ny);
T(:,1)=400;
T(1,:)=900;
T(:,end)=800;
T(end,:)=600;
T(1,1)=(900+400)/2;
T(1,end)=(900+800)/2;
T(end,1)=(400+600)/2;
T(end,end)=(800+600)/2;
Told=T;
k1=2*(dx^2+dy^2)/(dx^2*dy^2);
tol=1e-4;
error=9e9;
solver=input('Enter the number for the solver method: n Jacobi=1 n Gauss Seidel=2 n SOR=3 n Solver Type: ');
if solver == 1
tic
jacobi_iter=1;
while (error>tol)
for j=2:ny-1
for i=2:nx-1
H=((Told(i-1,j)+Told(i+1,j))/(k1*dx^2));
V=((Told(i,j-1)+Told(i,j+1))/(k1*dy^2));
T(i,j)=H+V;
end
end
error=max(max(abs(Told-T)))
Told=T;
jacobi_iter=jacobi_iter+1
timecounter=toc;
figure(1)
[C,h]=contourf(x,y,T,'ShowText','on');
colorbar
colormap(jet)
xlabel('xx');
ylabel('yy');
clabel(C,h)
title(sprintf('2d Steady State Heat Conduction n Method:Jacobi n No of iterations=%d n Time:%f s',jacobi_iter,timecounter))
end
end
if solver == 2
tic
gs_iter=1;
while (error>tol)
for j=2:ny-1
for i=2:nx-1
H=((T(i-1,j)+Told(i+1,j))/(k1*dx^2));
V=((T(i,j-1)+Told(i,j+1))/(k1*dy^2));
T(i,j)=H+V;
end
end
error=max(max(abs(Told-T)))
Told=T;
gs_iter=gs_iter+1
timecounter_1=toc;
figure(1)
[C,h]=contourf(x,y,T,'ShowText','on');
colorbar
colormap(jet)
xlabel('xx');
ylabel('yy');
clabel(C,h)
title(sprintf('2d Steady State Heat Conduction n Method:Gauss Seidel n No of iterations=%d n Time:%f s',gs_iter,timecounter_1))
end
end
if solver == 3
tic
sor_iter=1;
while (error>tol)
for j=2:ny-1
for i=2:nx-1
H=((T(i-1,j)+Told(i+1,j))/(k1*dx^2));
V=((T(i,j-1)+Told(i,j+1))/(k1*dy^2));
T(i,j)=(Told(i,j)*(1-omega))+(omega*(H+V));
end
end
error=max(max(abs(Told-T)))
Told=T;
sor_iter=sor_iter+1
timecounter_2=toc;
figure(1)
[C,h]=contourf(x,y,T,'ShowText','on');
colorbar
colormap(jet)
xlabel('xx');
ylabel('yy');
clabel(C,h)
title(sprintf('2d Steady State Heat Conduction Method:SOR No of iterations=%d Time:%f s',sor_iter,timecounter_2))
end
end
clear all
close all
clc
%% intializing variables
Nx = 51;
Ny = 51;
x = linspace (0,1,Nx);
y = linspace (0,1,Ny);
dx = 1/(Nx-1);
dy = 1/(Ny-1);
alpha = 0.0001;
dt = 1;
tol = 1e-4;
%% Temperature gride
T = 298*ones(Nx,Ny);
T(1,:) = 900; %bottom temp
T(end,:) = 600; %Top temp
T(:,1) = 400; %Left temp
T(:,end) = 800; %right temp
%edge refining
T(1,1) = (900+400)/2;
T(1,end) = (900+800)/2;
T(end,1) = (400+600)/2;
T(end,end) = (800+600)/2;
%% CFL number calculation and adaptive Time-step control
CFL_number = (alpha*dt)*((1/dx^2)+(1/dy^2));
%% Calculation of temperature distribution explicitly
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[xx,yy] = meshgrid(x,y);
Told = T;
T_prev_dt = T;
tic;
counter = 1;
for z = 1:2500
% Explicit Iteration Scheme
for i = 2:Nx-1
for j = 2:Ny-1
term1 = (1-2*k1-2*k2)*Told(i,j);
term2 = k1*(Told(i-1,j)+Told(i+1,j));
term3 = k2*(Told(i,j-1)+Told(i,j+1));
T(i,j) = term1+term2+term3;
end
end
Told = T;
counter = counter+1;
time_counter = toc;
figure (1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf('2D Unstedy state heat conduction Method:Explicit Method Number of iterations: %d(Time: %0.2fs) Computation Time: %0.4f CFL Number: %0.4f', counter-1,z*dt,time_counter,CFL_number));
xlabel('xx');
ylabel('yy');
end
clear all
close all
clc
%% Initializing Variable
Nx = 51;
Ny = 51;
x = linspace (0,1,Nx);
y = linspace (0,1,Ny);
dx = 1/(Nx-1);
dy = 1/(Ny-1);
alpha = 0.0001;
dt = 1;
tol = 1e-4;
%% Temperature gride
T = 298*ones(Nx,Ny);
T(1,:) = 900; %bottom temp
T(end,:) = 600; %Top temp
T(:,1) = 400; %Left temp
T(:,end) = 800; %right temp
%edge refining
T(1,1) = (900+400)/2;
T(1,end) = (900+800)/2;
T(end,1) = (400+600)/2;
T(end,end) = (800+600)/2;
%% CFL number calculation and adaptive Time-step control
CFL_number = (alpha*dt)*((1/dx^2)+(1/dy^2));
%% Calculation of temperature distribution using implicit method
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[xx,yy] = meshgrid(x,y);
Told = T;
T_prev_dt = T;
term1 = 1/(1+2*k1+2*k2);
term2 = k1*term1;
term3 = k2*term1;
tic;
counter = 1;
solver = input('Enter the number for the solver method: n Jacobi=1 n Gauss Seidal=2 n SOR=3 n Solver type: ');
% Jacobin Method
if solver == 1
for k = 1:2500
error = 9e9;
while(tol<error)
for j = 2:Ny-1
for i = 2:Nx-1
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;
counter = counter+1;
end
T_prev_dt = T;
figure (1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
time_counter = toc;
title(sprintf('2D Unsteady state heat conduction Method:Implicit Jacobi Number of iterations: %d (Time: %0.2fs) Computation Time: %0.4f CFL Number: %0.4f',counter-1,k*dt,time_counter,CFL_number));
xlabel('xx');
ylabel('yy');
end
end
% Gauss Seidel Method
if solver == 2
for k = 1:2500
error = 9e9;
while (tol<error)
for j = 2:Ny-1
for i = 2:Nx-1
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;
counter = counter+1;
end
T_prev_dt = T;
figure (1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
time_counter = toc;
title(sprintf('2D Unsteady state heat conduction Method:Implicit Gauss Seidel Number of iterations: %d (Time: %0.2fs) Computation Time: %0.4f CFL Number: %0.4f',counter-1,k*dt,time_counter,CFL_number));
xlabel('xx');
ylabel('yy');
end
end
%SOR method
if solver == 3
for k=1:2500
error = 9e9;
while(tol<error)
for j = 2:Ny-1
for i = 2:Nx-1
H= (T(i-1,j)+Told(i+1,j));
V= (T(i,j-1)+Told(i,j+1));
T(i,j) =(1-1.3)*(Told(i,j))+(1.3*((T_prev_dt(i,j)*term1)+(H*term2)+(V*term3)));
end
end
error=max(max(abs(Told-T)));
Told = T;
counter = counter+1;
end
T_prev_dt = T;
figure(1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
time_counter = toc;
title(sprintf('2D Unsteady state heat conduction n Method:Implicit SOR n Number of iterations: %d (Time: %0.2fs) n Computation Time: %0.4f n CFL number: %0.4f',counter-1,k*dt,time_counter,CFL_number));
xlabel('xx');
ylabel('yy');
end
end