How to smooth or clear the peak of the curve at certain time level?

19 ビュー (過去 30 日間)
Yanni
Yanni 2024 年 10 月 26 日
編集済み: John Kelly 2024 年 11 月 18 日 15:03
This is the tridaigonal system. here, i incoparate the time steps and it will stop at tolerance limit. and then i got the datas in T to plot graph in multiple time line. From begining all the curves coming downward(it has no problem). after the certain time level at 61th coloumn (in T workspace) it is coming peak. I have confused why its coming and how to clear or remove to get the smooth curve?
clc; close all; clear all;
%======================SPACEGRID====================================%
ymax=14; m=56; dy=ymax/m; y=dy:dy:ymax; %'i'th row
xmax=1; n=20; dx=xmax/n; x=dx:dx:xmax; %'j'th column
%=====================TIMEGRID======================================%
tmax=100; nt=500; dt=tmax/nt; t=dt:dt:tmax; % time at 'j'
%=====================STEADYSTATEINPUTVALUES========================%
tol=1e-2;
max_difference(1)=1; max_Iteration=1;
umax_difference(1)=1; %umax_Iteration=1;
%======================INPUTVALUES===================================%
pr=6.2;del=0;
%======================NFINPUTVALUES===============================%
PHI1=0; PHI2=0.; S=0;
rhoF=997.1; rhobetaF=20939.1; rhocpF=4166881; kF=0.613; sF=0.05;
rhoC=8933; rhobetaC=14918.11; rhocpC=3439205; kC=401; sC=59.6*10^6;
rhoZ=5600; rhobetaZ=2413.6; rhocpZ=2773120; kZ=13; sZ=1*10^-2;
KNF=kF*kC+S-1*kF-S-1*PHI1*kF-kC/kC+S-1*kF+PHI1*kF-kC;
KHNF=KNF*kZ+S-1*KNF-S-1*PHI2*KNF-kZ/kZ+S-1*KNF+PHI2*KNF-kZ;
SNF=sF*sC+S-1*sF-S-1*PHI1*sF-sC/sC+S-1*sF+PHI1*sF-sC;
SHNF=SNF*sZ+S-1*SNF-S-1*PHI2*SNF-sZ/sZ+S-1*SNF+PHI2*SNF-sZ;
RHOHNF=1-PHI1*1-PHI2+PHI1*rhoC/rhoF+PHI2*rhoZ/rhoF;
RHOCPHNF=1-PHI1*1-PHI2+PHI1*rhocpC/rhocpF+PHI2*rhocpZ/rhocpF;
RHOBETAHNF=1-PHI1*1-PHI2+PHI1*rhobetaC/rhobetaF+PHI2*rhobetaZ/rhobetaF;
E1=1/1-PHI1^2.5*1-PHI2^2.5*1/RHOHNF;
E2=RHOBETAHNF/RHOHNF;
E3=SHNF/RHOHNF;
E4=1/pr*KHNF/KNF*1/RHOCPHNF;
E5=1/RHOCPHNF;
%====================INTIALIZATION=================================%
UOLD=zeros(m,nt);
VOLD=zeros(m,nt);
TNEW=0; TOLD=TNEW*ones(m,nt); TWALL=ones(1,length(t));
A=zeros([1,m]);
B=A;
C=A;
D=A;
T=TOLD; U=UOLD; V=VOLD;
tic
%===================ENERGYEQUATION==============================%
while max_Iteration>tol
for j=1:nt
for i=1:m
if j>i
C(i)=dt*VOLD(i,j)/4*dy-dt*E4/(2*dy^2);
elseif i>j
A(i)=-dt*VOLD(i,j)/4*dy-dt*E4/(2*dy^2);
elseif i==j
B(i)=1+dt*UOLD(i,j)/2*dx+dt*E4/dy^2+dt*1/2*E5*del;
end
end
end
for j=2:nt
for i=1:m-1
if i==1
D(i)=TOLD(i,j)-dt*UOLD(i,j)*(-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1))/2*dx-dt*VOLD(i,j)*(TWALL(j)-TOLD(i-0,j))/4*dy+dt*E4*(TOLD(i-0,j)-2*TOLD(i,j)+TWALL(j)/2*dy^2)+(TOLD(i,j)/2*dt*E5*del)-dt*E4+E5*TWALL(j)/2*dy^2;
elseif i==m-1
D(i)=TOLD(i,j)-dt*UOLD(i,j)*(-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1))/2*dx-dt*VOLD(i,j)*(TOLD(i+1,j)-TNEW)/4*dy+dt*E4*(TNEW-2*TOLD(i,j)+TOLD(i+1,j)/2*dy^2)+(TOLD(i,j)/2*dt*E5*del)-dt*(E4+E5)*TNEW/2*dy^2;
else
D(i)=TOLD(i,j)-dt*UOLD(i,j)*(-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1))/2*dx-dt*VOLD(i,j)*(TOLD(i+1,j)-TOLD(i-1,j))/4*dy+dt*E4*(TOLD(i-1,j)-2*TOLD(i,j)+TOLD(i+1,j)/2*dy^2)+(TOLD(i,j)/2*dt*E5*del)-dt*TOLD(i,j)/2*dy^2;
end
end
T(:,j)=TriDiag(A,B,C,D);
dt=0.2+dt;
TOLD=T;
T(1,:)=TWALL(j);
%====================STEADY STATE=================================%
max_difference(j)=max(abs(T(:,j)-T(:,j-1))./max(1,abs(T(:,j-1))));
if max_difference(j)<tol
break
end
max_Iteration=max_difference;
end
end
function x = TriDiag(e,f,g,r)
%input
% e = subdiagonal vector
% f = diagonal vector
% g = superdiagonal vector
% r = right hand side vector
% output:
% x = solution vector
n=length(f);
for k = 2:n
factor = e(k)/f(k-1);
f(k) = f(k) - factor*g(k-1);
r(k) = r(k) - factor*r(k-1);
end
x(n) = r(n)/f(n);
for k = n-1:-1:1
x(k) = (r(k)-g(k)*x(k+1))/f(k);
end
for i = 1:length(x)
fprintf('\n%d = %f\n', i, x(i));
end
Kindly, help me in this hardest level.
Thank you.

回答 (1 件)

Karen
Karen 2024 年 10 月 28 日 11:46
編集済み: John Kelly 2024 年 11 月 18 日 15:03
Hey! It sounds like you’re seeing an unusual peak in your tridiagonal system’s graph around the 61st column. This could be due to boundary conditions or even some numerical instability as the time steps progress.

カテゴリ

Help Center および File ExchangeNetworks についてさらに検索

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by