I want Corrected code(FTCS) for nodal point formulas for 2D Transient with convection boundaries especially for Boundary and corner nodes.
I have searched for this in some other source it gives this format of code.
And it give results for 5 sec as max-295K, min-245K but im giving max-1503K ,min-787K as inputs.(Is it sounds good ?)
Please give feedback(Right/Wrong/Any modification)
%% FDMT for Convection boundary conditions
%% Units:
% All temperatures are in kelvin
%% Geomentric parameters
W=0.03175; % width of the plate(1.25 inch) (m)
H=0.2286; % Height of the plate(9 inch) (m)
%% Material properties(INCONEL 718- 1643K(M.P))
K=28.0285; % Thermal conductivity(W/m-K)
C=667.571; % Specific heat of materail(J/Kg-K)
den=7681.2857; % Density of material(Kg/m^3)
Alpha=K./(C.*den); % diffusivity of material(Inconel 718) (m/s^2)
%% Meshing
ni=116; % No.of nodes in i th direction
nj=16; % No.of nodes in j th direction
di=H./(ni-1); % Element size in i th direction (m)
dj=W./(nj-1); % Element size in j th direction (m)
dt=0.1; % Time step during take off
Fo=(Alpha.*dt)./(di.^2); % Fourier number
%% Initilizing the Matrix at time(t=0)
Tin=298.15; % Intial temp of the slab
for k=1:1
T=Tin.*ones(ni,nj,k); % Intial temp of the slab for all nodes
end
%% Intializing Boundary conditions
for t=0:4 % Total time(5sec)
%% Avg_Heat transfer coeffints
h1=[341.5929]; % Top surface
Bi1=(h1.*di)./K;
h2=[247.9693]; % Bottom surface
Bi2=(h2.*di)./K;
h3=[92.4127]; % Side surface
Bi3=(h3.*dj)./K;
Bi=max([Bi1,Bi2,Bi3]);
%% Stability criteria
p=Fo; % For interior nodes
q=Fo.*(2+Bi); % For surface nodes
q1=q;
r=Fo.*(1+Bi); % For corner nodes
r1=r;
if(p>0.25 && q1>0.5 && r1>0.25)
disp('unstable');
break
else
Th=[1503.4243];
Tl=[787.897];
%% Boundary nodes
for k=1:1
for i=1:ni
for j=1:nj
if(i==1 && j>=2 && j<=nj-1)
T(i,j,k)=(T(i+1,j,k)+(Th.*Bi1))./(1+Bi1); % Top
elseif(i==ni && j>=2 && j<=nj-1)
T(i,j,k)=(T(i-1,j,k)+(Tl.*Bi2))./(1+Bi2); % Bottom
elseif(j==nj && i>=2 && i<=ni-1)
T(i,j,k)=(T(i,j-1,k)+(Tl.*Bi3))./(1+Bi3); % Right
elseif(j==1 && i>=2 && i<=ni-1)
T(i,j,k)=(T(i,j+1,k)+(Tl.*Bi3))./(1+Bi3); % Left
end
end
end
end
for k=1:1
for i=1:ni
for j=1:nj
%% Corner nodes
if(i==1 && j==1)
T(i,j,k)=(T(i+1,j,k)+T(i,j+1))./2; % Top left
elseif(i==ni && j==1)
T(i,j,k)=(T(i-1,j,k)+T(i,j+1))./2; % Bottom left
elseif(i==ni && j==nj)
T(i,j,k)=(T(i-1,j,k)+T(i,j-1))./2; % Bottom Right
elseif(i==1 && j==nj)
T(i,j,k)=(T(i+1,j,k)+T(i,j-1))./2; % Top Right
end
end
end
end
%% Interior nodes
for k=1:1
for i=2:ni-1
for j=2:nj-1
T(i,j,k+1)=(T(i,j,k).*(1-(4.*Fo)))+((T(i+1,j,k)+T(i-1,j,k)+T(i,j+1,k)+T(i,j-1,k)).*(Fo));
end
end
end
%% Boundary nodes
for k=1:1
for i=1:ni
for j=1:nj
if(i==1 && j>=2 && j<=nj-1)
T(i,j,k+1)=(T(i+1,j,k+1)+(Th.*Bi1))./(1+Bi1); % Top
elseif(i==ni && j>=2 && j<=nj-1)
T(i,j,k+1)=(T(i-1,j,k+1)+(Tl.*Bi2))./(1+Bi2); % Bottom
elseif(j==nj && i>=2 && i<=ni-1)
T(i,j,k+1)=(T(i,j-1,k+1)+(Tl.*Bi3))./(1+Bi3); % Right
elseif(j==1 && i>=2 && i<=ni-1)
T(i,j,k+1)=(T(i,j+1,k+1)+(Tl.*Bi3))./(1+Bi3); % Left
end
end
end
end
%% Corner nodes
for k=1:1
for i=1:ni
for j=1:nj
if(i==1 && j==1)
T(i,j,k+1)=(T(i+1,j,k+1)+T(i,j+1,k+1))./2; % Top left
elseif(i==ni && j==1)
T(i,j,k+1)=(T(i-1,j,k+1)+T(i,j+1,k+1))./2; % Bottom left
elseif(i==ni && j==nj)
T(i,j,k+1)=(T(i-1,j,k+1)+T(i,j-1,k+1))./2; % Bottom Right
elseif(i==1 && j==nj)
T(i,j,k+1)=(T(i+1,j,k+1)+T(i,j-1,k+1))./2; % Top Right
end
end
end
end
%% Boundary nodes
for k=1:50
for i=1:ni
for j=1:nj
if(i==1 && j>=2 && j<=nj-1)
T(i,j,k+1)=((Fo).*((2.*T(i+1,j,k))+T(i,j+1,k)+T(i,j-1,k)+(2.*Bi1.*Th)))+((1-(4.*Fo)-(2.*Bi1)).*(T(i,j,k))); % Top
elseif(i==ni && j>=2 && j<=nj-1)
T(i,j,k+1)=((Fo).*((2.*T(i-1,j,k))+T(i,j+1,k)+T(i,j-1,k)+(2.*Bi2.*Tl)))+((1-(4.*Fo)-(2.*Bi2)).*(T(i,j,k))); % Bottom
elseif(j==nj && i>=2 && i<=ni-1)
T(i,j,k+1)=((Fo).*(T(i-1,j,k)+T(i+1,j,k)+(2.*T(i,j-1,k))+(2.*Bi3.*Tl)))+((1-(4.*Fo)-(2.*Bi3)).*(T(i,j,k))); % Right
elseif(j==1 && i>=2 && i<=ni-1)
T(i,j,k+1)=((Fo).*(T(i-1,j,k)+T(i+1,j,k)+(2.*T(i,j+1,k))+(2.*Bi3.*Tl)))+((1-(4.*Fo)-(2.*Bi3)).*(T(i,j,k))); % Left
end
end
end
%% Corner nodes
for i=1:ni
for j=1:nj
if(i==1 && j==1)
T(i,j,k+1)=((2.*Fo)*(T(i+1,j,k)+T(i,j+1,k)+(Bi1.*Th)+(Bi3.*Tl))+((1-(4.*Fo)-(Bi1.*2.*Fo)-(Bi3.*2.*Fo)).*T(i,j,k))); % Top left
elseif(i==ni && j==1)
T(i,j,k+1)=((2.*Fo).*(T(i-1,j,k)+T(i,j+1,k)+(Bi2.*Tl)+(Bi3.*Tl))+((1-(4.*Fo)-(Bi2.*2.*Fo)-(Bi3.*2.*Fo)).*T(i,j,k))); % Bottom left
elseif(i==ni && j==nj)
T(i,j,k+1)=((2.*Fo).*(T(i-1,j,k)+T(i,j-1,k)+(Bi2.*Tl)+(Bi3.*Tl))+((1-(4.*Fo)-(Bi2.*2.*Fo)-(Bi3.*2.*Fo)).*T(i,j,k))); % Bottom Right
elseif(i==1 && j==nj)
T(i,j,k+1)=((2.*Fo).*(T(i+1,j,k)+T(i,j-1,k)+(Bi1.*Tl)+(Bi3.*Tl))+((1-(4.*Fo)-(Bi1.*2.*Fo)-(Bi3.*2.*Fo)).*T(i,j,k))); % Top Right
end
end
end
%% Last but one boundary nodes
for i=1:ni
for j=1:nj
if(i==2 && j>=2 && j<=nj-1)
T(i,j,k+1)=(T(i,j,k).*(1-(4.*Fo)))+((T(i+1,j,k)+T(i-1,j,k+1)+T(i,j+1,k)+T(i,j-1,k)).*(Fo));
elseif(i==ni-1 && j>=2 && j<=nj-1)
T(i,j,k+1)=(T(i,j,k).*(1-(4.*Fo)))+((T(i+1,j,k+1)+T(i-1,j,k)+T(i,j+1,k)+T(i,j-1,k)).*(Fo));
elseif(j==nj-1 && i>=2 && i<=ni-1)
T(i,j,k+1)=(T(i,j,k).*(1-(4.*Fo)))+((T(i+1,j,k)+T(i-1,j,k)+T(i,j+1,k+1)+T(i,j-1,k)).*(Fo));
elseif(j==2&& i>=2 && i<=ni-1)
T(i,j,k+1)=(T(i,j,k).*(1-(4.*Fo)))+((T(i+1,j,k)+T(i-1,j,k)+T(i,j+1,k)+T(i,j-1,k+1)).*(Fo));
end
end
end
%% Interior nodes
for i=2:ni-1
for j=2:nj-1
T(i,j,k+1)=(T(i,j,k).*(1-(4.*Fo)))+((T(i+1,j,k)+T(i-1,j,k)+T(i,j+1,k)+T(i,j-1,k)).*(Fo));
end
end
%% Plotting
x=linspace(1,16,16);
y=linspace(1,116,116);
[X,Y]=meshgrid(x,y);
z=T(:,:,k);
z1=rot90(z,2);
h=surf(X,Y,z);
colorbar
colormap("jet(50)")
shading interp
title({'Transient state-Explicit temperature profile';['Time=',num2str(t+(k.*dt))]})
pause(0.1)
refreshdata(h)
grid on
end
end
end
