41 views (last 30 days)

Show older comments

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
on 26 Feb 2021

Edited: David Goodmanson
on 26 Feb 2021

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.

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

Start Hunting!