Provide Corrected nodal point formulas for 2D Transient with convection boundaries
2 次查看(过去 30 天)
显示 更早的评论
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
0 个评论
回答(1 个)
Anagha Mittal
2023-9-20
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!
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Multiobjective Optimization 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!