2d heat conduction
code-1
% MATLAB Program to generate temperature profile of 2D Transient State Heat Conduction Equation explicitly:
clear all;
close all;
clc;
% Initializing Variables:
Nx = 51;
Ny = Nx;
%unit square
x = linspace(0,1,Nx);
y = linspace(0,1,Ny);
dx = 1/(Nx-1);
dy = 1/(Ny-1);
alpha = 1e-4;
dt = 0.01;
error = 6e9;
tol = 1e-4;
% Initializing the Temperature grid:
T = 298*ones(Nx,Ny);
T(1,:) = 900; % Bottom Temperature
T(end,:) = 600; % Top Temperature
T(:,1) = 400; % Left Temperature
T(:,end) = 800; % Right Temperature
% 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/(dx^2);
% Calculation of Temperature Distribution explicitly using iterative solvers:
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[xx,yy] = meshgrid(x,y);
Told = T;
counter = 1;
tic
for k = 1:20000
% Gauss Jacobi explicit Iteration Scheme
for i = 2:Nx-1
for j = 2:Ny-1
term1 = Told(i,j);
term2 = k1*(Told(i-1,j)-2*Told(i,j)+Told(i+1,j));
term3 = k2*(Told(i,j-1)-2*Told(i,j)+Told(i,j+1));
T(i,j) = term1+term2+term3;
end
end
Told = T;
counter = counter+1;
end
exe_time=toc;
%plotting of the results
figure(1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf(\'2D Unsteady state heat conduction (Explicit Method) / Number of iterations: %d. / Computation Time: %0.2f./ CFL Number: %0.4f./ execution time: %0.4f.\',counter-1,k*dt,CFL_Number,exe_time));
xlabel(\'xx\');
ylabel(\'yy\');
code-2
% Unsteady state 2D heat conduction by implicit method.3
clear all
close all
clc
%Initialization and grid formation
nx=51;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=1/(nx-1);
dy=1/(ny-1);
alpha = 1e-4;%Thermal Diffusibility
t_start = 0;
t_end = 1/alpha;
dt = 0.01;
error = 1e9; %Initializing error value
tol = 1e-4; %tolerence
CFL_Number = alpha*dt/(dx^2);%courant number
Nt = t_start:dt*50:t_end; %Time spacing
%BOUNDARY CONDITIONS
T=298*ones(nx,ny);
T(1,:) = 900; % Bottom Temperature
T(end,:) = 600; % Top Temperature
T(:,1) = 400; % Left Temperature
T(:,end) = 800; % Right Temperature
% Refinementing mesh at the edges
T(1,1) = (900+400)/2;
T(1,end) = (900+800)/2;
T(end,1) = (400+600)/2;
T(end,end) = (800+600)/2;
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[xx,yy] = meshgrid(x,y);
Told = T;
T_prev_dt = Told;
counter = 1;
solver = input(\'Enter Number: \')
tic;
for k = 1:20000
%GAUSS JACOBI SOLVER
if solver == 1
term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
error = 1e9;
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
H = (Told(i-1,j)+Told(i+1,j));
V = (Told(i,j-1)+Told(i,j+1));
T(i,j) = term1*Told(i,j)+term2*H+term3*V;
end
end
error = max(max(abs(Told-T)));
Told = T;
counter = counter + 1;
end
end
%GAUSS SIEDEL SOLVER
if solver == 2
term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
error = 1e9;
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
H = (T(i-1,j)+Told(i+1,j));
V = (T(i,j-1)+Told(i,j+1));
T(i,j) = term1*T_prev_dt(i,j)+term2*H+term3*V;
end
end
error = max(max(abs(Told-T)));
Told = T;
counter = counter + 1;
end
end
%SOR SOLVER
if solver == 3
a = 0.9; %Relaxation Parameter
term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
error = 1e9;
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
H=((1-a)*T(i-1,j)+a*(T(i-1,j)+Told(i+1,j)));
V=((1-a)*T(i,j-1)+a*(T(i,j-1)+Told(i,j+1)));
T(i,j)=(T_prev_dt(i,j))*term1+(term2*H)+(term3*V);
end
end
error = max(max(abs(Told-T)));
Told = T;
counter = counter + 1;
end
end
T_prev_dt = T;
end
execution_time = toc;
%plotting the graphs;
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf(\'2D Unsteady state heat conduction (implicit Method) / Number of iterations: %d / Computation Time: %0.2f/ CFL Number: %0.4f\',counter-1,k*dt,CFL_Number));
xlabel(\'xx\');
ylabel(\'yy\');
code-3
%Code for steady state 2D heat equation (SOLVING IMPLICITLY)
clear all
close all
clc
%Initialization and grid formation
nx=51;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=1/(nx-1);
dy=1/(ny-1);
k = 2*(1/(dx^2)+1/(dy^2));
acc = 1e-4;
error = 1e8;%Initializing error value
%BOUNDARY CONDITIONS
T=298*ones(nx,ny); % Setting standard temperature at all nodes initially
T(1,:) = 900; % Bottom Temperature
T(end,:) = 600; % Top Temperature
T(:,1) = 400; % Left Temperature
T(:,end) = 800; % Right Temperature
% Refinement of the created mesh at the edges
T(1,1) = (900+400)/2;
T(1,end) = (900+800)/2;
T(end,1) = (400+600)/2;
T(end,end) = (800+600)/2;
%Design of the iterative solver
x=y;
[xx,yy] = meshgrid(x,y);
Told = T;
SOLVER = input(\'Enter Number: \')
tic;
%GAUSS JACOBI SOLVER
if SOLVER == 1
iteration_count_jacobi = 1;
while (error > acc)
for i = 2:nx-1
for j = 2:ny-1
%Discretized Terms of the Equation
term1 = (Told(i-1,j)+Told(i+1,j))/(k*(dx^2));
term2 = (Told(i,j-1)+Told(i,j+1))/(k*(dy^2));
T(i,j) = term1+term2;
end
end
error = max(max(abs(Told-T)));
Told = T;
iteration_count_jacobi = iteration_count_jacobi + 1;
end
end
%GAUSS SIEDEL SOLVER
if SOLVER == 2
iteration_count_siedel = 1;
while (error > acc)
for i = 2:nx-1
for j = 2:ny-1
%Discretized Terms of the Equation
term1 = (T(i-1,j)+T(i+1,j))/(k*(dx^2));
term2 = (T(i,j-1)+T(i,j+1))/(k*(dy^2));
T(i,j) = term1+term2;
end
end
error = max(max(abs(Told-T)));
Told = T;
iteration_count_siedel = iteration_count_siedel + 1;
end
end
%SUCCESSIVE OVER RELAXATION SOLVER
if SOLVER == 3
a = 2/(1+sin(pi*dx)); %Relaxation Parameter
iteration_count_SOR = 1;
while (error > acc)
for i = 2:nx-1
for j = 2:ny-1
%Discretized Terms of the Equation
term1 = (T(i-1,j)+Told(i+1,j))/(k*(dx^2));
term2 = (T(i,j-1)+Told(i,j+1))/(k*(dy^2));
T(i,j) = a*(term1+term2)-(a-1)*Told(i,j);
end
end
error = max(max(abs(Told-T)));
Told = T;
iteration_count_SOR = iteration_count_SOR + 1;
end
end
time_elapsed = toc;
figure(1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
if (SOLVER == 1)
title(sprintf(\'Temperature profile of steady state heat conduction using GAUSS JACOBI METHOD. \\n No. of Iterations: %d. \\n Computation Time: %0.3f s.\',iteration_count_jacobi,time_elapsed));
end
if (SOLVER == 2)
title(sprintf(\'Temperature profile of steady state heat conduction using GAUSS SIEDEL METHOD. \\n No. of Iterations: %d. \\n Computation Time: %0.3f s.\',iteration_count_siedel,time_elapsed));
end
if (SOLVER == 3)
title(sprintf(\'Temperature profile of steady state heat conduction using SOR method. \\n No. of Iterations: %d. \\n Computation Time: %0.3f s. \\n Relaxation Parameter: %f\',iteration_count_SOR,time_elapsed,real(a)));
end
xlabel(\'xx\');
ylabel(\'yy\');Transient state or unsteady state by Explicit method
Jacobi method:-
Explicit method of the Jacobi method has taken very less time to converge than that of the implicit method. iteration are within in the specified range.
Transient state or unsteady state by implicit method
1.Jacobi method.
Iterations have crossed the specified range of 20k and reached around 265k steps to converge. This is because of the of many unknown values where are as in the explicit methid the values most of the variables are known and it is easy to calculate.
2. Gauss seidel
This method has taken less time to converge among other methods in implicit form but still the iteration count is high about 50k steps.
3. Successive over relaxation:-
For this method I have experimented various values of omega \'a\' from 0.2 to 2 but the solution seems to converge only when the value of the \'a\' is between 0.2 to 0.9.
faster convergence was observed at the 0.9 of omega value. This is the slowest of all the methods. Finding the right omega value is the challenge in this method.
Steady state by implicit method
1.Jacobi method:-
In steady state the convergence of the solution is very fast as there is no time time involved in it.
2. Gauss seidel:-
In this method the convergence taking less time and steps in all the cases.
3. Successive over relaxation:-
SOR is performing very fast in the steady state. By this we can understand that time poses challenge to this method for converging.