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.

Answer 1

[David Goodmanson]

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.

[Yong Khee Foo]

Thank you so much, this really helped a lot!

And yes, Xgf --> 1 is desirable.