2D Heat conduction equation solver
STEADY STATE................................................................
In steady state, we take the initial guess of temperature at steady state.
But the steady state is unique and it doesn\'t depend on the initial guess.
TRANSIENT STATE(IMPLICIT)............................................
In transient state, we take the actual state of temperature at time t=0(first time step).
So, the state of temperature at any point of time before steady state depends on initial conditions. Steady state is unique in this case also.
In this method, we make use of values at (n+1)th time step which we don\'t know beforehand to calculate values at (n+1)th time step. Thus, we solve it iteratively.
TRANSIENT STATE(EXPLICIT)...........................................
In this method, we make use of values at nth time step which we know already to calculate values at (n+1)th time step. Thus, we don\'t need to solve it iteratively.
Program for solving problem:
clear all
close all
clc
n_x=51;
n_y=51;
n_t=501;
x=linspace(0,1,n_x);
y=linspace(0,1,n_y);
t=500;
dx=1/(n_x-1);
dy=1/(n_y-1);
dt=t/(n_t-1);
tol=1e-4;
T=298*ones(n_x,n_y);
%Boundary conditions
T(:,1)=400; %Left wall
T(:,end)=800; %Right wall
T(1,:)=900; %Bottom wall
T(end,:)=600; %Top wall
%Assigning average temperatures to corners of domain
T(1,1)=(900+400)/2;
T(end,end)=(800+600)/2;
T(1,end)=(800+900)/2;
T(end,1)=(400+600)/2;
%{
TRANSIENT IMPLICIT(process=1)
TRANSIENT EXPLICIT(process=2)
STEADY(process=3)
%}
%{
%JACOBI(method=1)
%GAUSS SEIDEL(method=2)
%SOR(method=3)
%}
%Define your requirement below
process=3;
method=3;
if process==1
T_new=transient_implicit(n_x,n_y,n_t,x,y,t,tol,T,dx,dy,dt,method);
elseif process==2
T_new=transient_explicit(n_x,n_y,n_t,x,y,t,T,dx,dy,dt);
elseif process==3
T_new=steady(n_x,n_y,x,y,tol,T,dx,dy,method);
end
Function for Transient Implicit solver:
function out=transient_implicit(n_x,n_y,n_t,x,y,t,tol,T,dx,dy,dt,method)
Told=T; %Initial temperature field at t=0
T_prev=T;
counter=0; %Initializing iteration counter
alpha=1e-4; %Thermal diffusivity
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
k3=1+2*k1+2*k2;
c1=1/k3;
c2=k1/k3;
c3=k2/k3;
% Time Loop
for k=1:n_t
error=1;
% Convergence Loop
while error>tol
if method==1
% Nodal Loop
for i=2:n_x-1
for j=2:n_y-1
T(i,j)=c1*(T_prev(i,j))+c2*(Told(i-1,j)+Told(i+1,j))+...
c3*(Told(i,j-1)+Told(i,j+1));
end
end
A=T-Told;
error=max(abs(A(:)));
Told=T;
counter=counter+1;
title_name=[\"Transient Implicit Jacobi method\",\'Time=\',num2str((k-1)*dt),\'seconds\'];
elseif method==2
% Nodal Loop
for i=2:n_x-1
for j=2:n_y-1
T(i,j)=c1*(T_prev(i,j))+c2*(T(i-1,j)+T(i+1,j))+...
c3*(T(i,j-1)+T(i,j+1));
end
end
A=T-Told;
error=max(abs(A(:)));
Told=T;
counter=counter+1;
title_name=[\"Transient Implicit GS method\",\'Time=\',num2str((k-1)*dt),\'seconds\'];
elseif method==3
% Nodal Loop
for i=2:n_x-1
for j=2:n_y-1
Tgs(i,j)=c1*(T_prev(i,j))+c2*(T(i-1,j)+T(i+1,j))+...
c3*(T(i,j-1)+T(i,j+1));
a=1.5;
T(i,j)=(1-a)*T(i,j)+a*Tgs(i,j);
end
end
A=T-Told;
error=max(abs(A(:)));
Told=T;
counter=counter+1;
title_name=[\"Transient Implicit SOR method\",\'Time=\',num2str((k-1)*dt),\'seconds\'];
end
end
T_prev=T;
figure(1)
surf(x,y,T)
shading interp
colorbar
title(title_name)
xlabel(\'x\',\'FontSize\',12,\'FontWeight\',\'bold\')
ylabel(\'y\',\'FontSize\',12,\'FontWeight\',\'bold\')
axis([0 1 0 1 0 900])
view(0,90)
pause(0.001)
end
out=T;
endFunction for Transient Explicit solver:
function out=transient_explicit(n_x,n_y,n_t,x,y,t,T,dx,dy,dt)
T_prev=T;
counter=0; %Initializing iteration counter
alpha=1e-4; %Thermal diffusivity
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
% Time Loop
for k=1:n_t
% Nodal Loop
for i=2:n_x-1
for j=2:n_y-1
T(i,j)=T_prev(i,j)+k1*(T_prev(i-1,j)-2*T_prev(i,j)+T_prev(i+1,j))+...
k2*(T_prev(i,j-1)-2*T_prev(i,j)+T_prev(i,j+1));
end
end
T_prev=T;
figure(1)
surf(x,y,T)
shading interp
colorbar
title_name=[\"Transient Explicit\",\'Time=\',num2str((k-1)*dt),\' seconds\'];
title(title_name)
xlabel(\'x\',\'FontSize\',12,\'FontWeight\',\'bold\')
ylabel(\'y\',\'FontSize\',12,\'FontWeight\',\'bold\')
axis([0 1 0 1 0 900])
view(0,90)
pause(0.001)
end
out=T;
end
Function for Steady state solver:
function out=steady(n_x,n_y,x,y,tol,T,dx,dy,method)
Told=T; %Initial temperature field guess
error=1; %Initializing error value
counter=0; %Initializing iteration counter
k=(2*((dx)^2+(dy)^2))/((dx^2)*(dy^2));
% Convergence Loop
while error>tol
if method==1
% Nodal Loop
for i=2:n_x-1 %Calculating value for all nodes iteratively by taking an initial guess
for j=2:n_y-1
T(i,j)=((Told(i-1,j)+Told(i+1,j))/(dx^2)+(Told(i,j-1)+Told(i,j+1))/(dy^2))/k;
end
end
A=T-Told;
error=max(abs(A(:)));
Told=T;
counter=counter+1;
title_name=[\"Steady state\",\"Jacobi method Iteration no \",num2str(counter)];
elseif method==2
% Nodal loop
for i=2:n_x-1 %Calculating value for all nodes iteratively by taking an initial guess
for j=2:n_y-1
T(i,j)=((T(i-1,j)+T(i+1,j))/(dx^2)+(T(i,j-1)+T(i,j+1))/(dy^2))/k;
end
end
A=T-Told;
error=max(abs(A(:)));
Told=T;
counter=counter+1;
title_name=[\"Steady state\",\"GS method Iteration no \",num2str(counter)];
elseif method==3
% Nodal Loop
for i=2:n_x-1 %Calculating value for all nodes iteratively by taking an initial guess
for j=2:n_y-1
Tgs(i,j)=((T(i-1,j)+T(i+1,j))/(dx^2)+(T(i,j-1)+T(i,j+1))/(dy^2))/k;
a=1.5;
T(i,j)=T(i,j)*(1-a)+a*Tgs(i,j);
end
end
A=T-Told;
error=max(abs(A(:)));
Told=T;
counter=counter+1;
title_name=[\"Steady state\",\"SOR method Iteration no \",num2str(counter)];
end
end
out=T;
surf(x,y,T)
shading interp
colorbar
xlabel(\'x\',\'FontSize\',12,\'FontWeight\',\'bold\')
ylabel(\'y\',\'FontSize\',12,\'FontWeight\',\'bold\')
title(title_name)
axis([0 1 0 1 0 900])
view(0,90)
endThe number of iterations required for convergence decreases in the order:
Jacobi>Gauss seidel>SOR
This is because GS method uses updated values in calculation while SOR method amplifies the solution by some factor.