Imaginary numbers when solving ode (ode15s)
显示 更早的评论
I tried solving 4 first-order differential equations simultaneously using ode15s. 2 of them are mass balance equations in a counter-current moving bed reactor, while the other 2 accounts for the heat balance.
1 function file is used containing the 4 ODEs and relevant values. 1 script file was used to run the function file to solve the ODEs.
Function file:
function dvar= movingbed_scaleupFR (z, var)
Xgf= var(1) ; Xsf= var(2); Tgf= var(3); Tsf= var(4);
%Xgf= CH4 conversion in the FR; Xsf= Solid conversion in the FR
%Tgf= CH4 temp in the FR; Tsf= Solid temp in the FR
Cgf= Cgofuel*(1-Xgf)*Tgofuel/((1+e*Xgf)*Tgf); %methane concentration in the FR
n= 0.681; %order of reaction
%solve for Sm
Np= rhobulk/mp; %number of particles per unit volume
Sm= Np*Ap; %surface area per unit volume
%%%%%% Sm solved %%%%%%%
%Cp equations
%solving for the ODEs
dXgfdz=(A/Fgofuel)*(rhos/b)*kf*Cgf^n;
dXsfdz=(A/Fsofuel)*(-1*rhos*kf)*Cgf^n;
dvar=[dXgfdz; dXsfdz; dTgfdz; dTsfdz];
end
Script:
function to_run_movingbed_scaleupFR
var0= [0, 1, 1223, 1223]; %Initial conditions for [Xgf, Xsf, Tgf, Tsf]
bedheight= [0 20]; %setting bed height
[z,y]= ode15s(@movingbed_scaleupFR, bedheight, var0);
%Graphing plotting
yyaxis left
hold on
plot(z,(y(:,1)),'-r','displayname','Xgf')
plot(z,(y(:,2)),'-b','displayname','Xsf')
legend;
ylabel('conversion');
xlabel('Length/m');
hold off
end
Running the script, i get imaginary parts as shown below.
Warning: Imaginary parts of complex X and/or Y arguments ignored.
> In to_run_movingbed_scaleupFR (line 11)
Warning: Imaginary parts of complex X and/or Y arguments ignored.
> In to_run_movingbed_scaleupFR (line 12)
Warning: Imaginary parts of complex X and/or Y arguments ignored.
> In to_run_movingbed_scaleupFR (line 20)
Warning: Imaginary parts of complex X and/or Y arguments ignored.
> In to_run_movingbed_scaleupFR (line 21)
Final bed height whereby complete conversion occurs seemed wonky too, I was expecting bed height with a much longer length. So I assume the imaginary numbers along the way affected the output thereafter.
I'm not sure where I did wrong, I've checked the equations and values multiple times and can't find any mistakes.
采纳的回答
David Goodmanson
2021-2-26
编辑：David Goodmanson
2021-2-26
Hello YF,
whenever unexpected complex quantities start appearing, it's time to look for the usual suspects, which tend to be either negative quantities raised to fractional powers, or inverse trig functions with arguments outside the domain +-1. In this case there are three such lines
h= 0.33*Nre^(1/3)*kg/dp; %heat transfer coefficient, W/(m^2.K)
Xgfdz=(A/Fgofuel)*(rhos/b)*kf*Cgf^n;
dXsfdz=(A/Fsofuel)*(-1*rhos*kf)*Cgf^n;
Nre appears to only directly involve positive quanities, so by far the most suspicious line is
Cgf= Cgofuel*(1-Xgf)*Tgofuel/((1+e*Xgf)*Tgf); %methane concentration in the FR
n is not an integer, so if Xgf increases to become greater than 1, things go complex. And that is what happens. After running the code this can be seen with
figure(2); grid on
Xgf = y(:,1);
n = numel(Xgf);
plot(1:n,real(Xgf),1:n,1e6*imag(Xgf))
% show numbers
format long
real(Xgf)-1
(if you run the code more than once in succession, adding 'figure(1)' before your set of plots helps in keeping that figure from intruding into figure 2).
Is it true that Xgf --> 1 is physically correct? If so, you could omit all the points after Xgf first passes 1 (on the grounds of numerical imprecision) and call it good. If Xgf --> 1 is not correct then it looks like there will be more work to do.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
