Runge-Kutta - Plot Graphic - Help

Question

0 个投票

Hey guys! I wrote a code to study the dynamic of some parameters in a reactor. The problem is that I applied RK, but one of the variables is an explicit equation and, therefore, I'm having trouble to figure it out how to plot the graphic of Tcw2. Can somebody give me an idea? I need Tcw2 versus t (time).....

Thank you! Best regards!

function dx=EstudoDeCaso_RK4(t,x)
%Legenda
%x(1) ---> V
%x(2) ---> h
%x(3) ---> Ca
%x(4) ---> Cb
%x(5) ---> Cc
%x(6) ---> Cd
%x(7) ---> Cw
%x(8) ---> Tr
%Parâmetros
Vu=4;                   %m³
d=(2*Vu/pi)^(1/3);      %m
H=2*d;                  %m  
A=pi*d^2/4;             %m²
Vo=0.2;                 %m³
ve=0.004;               %m³/s
Cwe=55;                 %kmol/m³
Cbe=1;                  %kmol/m³
Fbe=Cbe*ve*1000;        %mol/s
MMa=88/1000;            %kg/mol
MMw=18/1000;            %kg/mol
MMb=40/1000;            %kg/mol
MMmist=(Cbe/(Cbe+Cwe)*MMb+Cwe/(Cbe+Cwe)*MMw); %kg/mol
Cp=75.24;               %J/(kmol.K)
Cpa=170.7;              %J/(kmol.K)
theta_w=Cwe/Cbe;        %(-)
UA=3000;                %J/(s.K)
To=300;                 %K
Tref=273;               %K
Tcw1=285;               %K
mcw=100;                %kg/s
dH=-7.9076E+07;         %J/kmol
ko=0.39175;             %m³/(kmol.s)
E=5472.7;               %K^-1  ----> Ea/R
%Equações Explícitas
Na=x(3)*x(1);
Nb=x(4)*x(1);
Nc=x(5)*x(1);
Nd=x(6)*x(1);
Nw=x(7)*x(1);
NCp=Cp*(Nb+Nc+Nd+Nw)+Cpa*Na;
kd=ko*exp(E*(1/273-1/x(8)));
Kc=10^(3885.44/x(8));
ra=kd*(x(5)*x(6)/Kc-x(3)*x(4));
rb=ra;
rc=kd*(x(3)*x(4)-x(5)*x(6)/Kc);
rd=rc;
Tcw2=x(8)-(x(8)-Tcw1)*exp(-UA/(mcw*Cp));
Qdot=mcw*Cp*(Tcw1-x(8))*(1-exp(-UA/(mcw*Cp)));
%Equações Diferenciais
dx(1)=ve;
dx(2)=ve/A;
dx(3)=ra-ve/x(1)*x(3);
dx(4)=rb+ve/x(1)*(Cbe-x(4));
dx(5)=rc-ve/x(1)*x(5);
dx(6)=rd-ve/x(1)*x(6);
dx(7)=ve/x(1)*(Cwe-x(7));
dx(8)=(Qdot-Fbe*Cp*(1+theta_w)*(x(8)-To)-ra*x(1)*dH)/NCp;
%Vetor linha em vetor coluna
dx=dx';
end