Provide Corrected nodal point formulas for 2D Transient with convection boundaries

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

Answer 1

[Anagha Mittal]

Hi Viswa,

The code mostly seems correct and should work as per the desired result. However, please note the following modifications that can be made to the above code for improving it:

1. You may revise the stability criteria condition as follows:

if(p>0.25 && q1>0.5 && r1>0.25)

This line will only trigger if all three conditions are true, but you want to trigger instability if any of these conditions are true. This should be modified to:

if(p>0.25 || q1>0.5 || r1>0.25)    

2. You may consider revising the equations for the corner nodes to better incorporate the boundary conditions.The script contains a large block of code for updating temperatures at the last iteration (50). It seems to be somewhat different from the previous temperature updates. Verify if these equations are indeed what you intend, particularly with the mixed use of k and k+1 time levels within the equations.

3. Also, a generic suggestion is to remove duplicated code and possibly create functions for repetitive tasks, particularly for updating the temperatures of nodes here. You could consolidate the repeated sections into a single section or function to avoid duplication, which would make your script more maintainable and less prone to errors.

Hope this helps!