Why i cant solve the 3 odes

I dont know why it is not working:

%% Clean the stage

clear;

close all;

clc;

%% Identification the global variables

global k2 Ka k3 ndA0 ndB0 alpha R P0 Lambda CpA CpB T0 eps v0

%% Declaration of the variables

P0 = 150; % Inlet pressure, atm

alpha = 0.6004; % Pressure drop parameter, 1/(kg catalyst)

k2=4.6*10^8;

Ka= 0.010985321;

k3=2.96*10^-6;

Lambda=-109;

R=8.314*1000;

CpA= 0.02919;

CpB=0.02922;

eps=-0.5;

T0=703;

ndA0=1798.94;

ndB0=3*ndA0;

v0=1000;

%% Declaration of the weight span and the initial conditions

W = [0 60]; % Catalyst weight span, kg catalyst

xi = [0 1 703]; % xi(1)=eps at the inlet, xi(2)=P/P0 at inlet xi(3)=Ti

%% Solution of ODE system

[W,x]=ode113(@funisopbr, W, xi);

%% Extraction and plotting of the results

% Calculation of molar flowrates

ndA = ndA0-ndA0*x(:,1); % Molar flowrate of A

ndB = ndB0-ndA0*x(:,1); % Molar flowrate of B

ndC = ndA0*x(:,1); % Molar flowrate of C

% Plot of molar flowrates vs. weight of catalyst

figure(1)

plot(W,ndA,'r', W,ndB,'k',W,ndC, 'b','linewidth', 2);

title('Molar flowrates vs. weight of catalyst');

xlabel('Weight of catalyst, kg');

ylabel('Molar flowrate per tube, mol/s');

legend('nd_A','nd_B','nd_C');

% Plot of pressure vs. weight of catalyst

figure(2)

plot(W,x(:,2)*10,'m','linewidth',2);

title('Pressure vs. weight of catalyst');

xlabel('Weight of catalyst, kg');

ylabel('Pressure, atm');

% Plot of temperature vs. volume

figure(2)

plot(W,x(:,3),'linewidth',2);

title('Temperature vs. volume');

xlabel('Volume, L');

ylabel ('Temperature, K');

[21:23, 29.5.2019] Mayg: function dxdw=funisopbr(w,x)

%% Identification of global variables

global k2 Ka k3 ndA0 ndB0 alpha R Lambda CpA CpB Pa Pb Pc ra T0 eps v0

%% Decleration of ODE system

% x(1)=conv, x(2)=y x(3)=T

Pa=ndA0/v0*(1-x(1)/(1+eps*x(1)))*x(2)*T0/x(3)*R*x(3);

Pb=ndA0/v0*(3-3*x(1))/(1+eps*x(1))*x(2)*T0/x(3)*R/x(3);

Pc=ndA0/v0*(2*x(1)/(1+eps*x(1)))*x(2)*T0/x(3)*R/x(3);

ra =(k2*(1/3*Pa*Ka^2-Pc^2/Pb)/(1+k3*Pc/Pb^1.5)^0.75);

dxdw=[-ra/ndA0;

-alpha*x(3)/T0/(2*y)*(1+eps*x(1));

ra*Lambda/(ndA0*CpA+1/3*ndB0*CpB)];

end

[Star Strider]

This version of your code (without the global calls) runs. It just takes forever, so after 45 minutes on a fast computer without results (MATLAB is still ‘busy’), I stopped it:

function dxdw=funisopbr(w,x, k2, Ka, k3, ndA0, ndB0, alpha, R, Lambda, CpA, CpB, T0, eps, v0)
%% Decleration of ODE system
% x(1)=conv, x(2)=y x(3)=T
Pa=ndA0/v0*(1-x(1)/(1+eps*x(1)))*x(2)*T0/x(3)*R*x(3);
Pb=ndA0/v0*(3-3*x(1))/(1+eps*x(1))*x(2)*T0/x(3)*R/x(3);
Pc=ndA0/v0*(2*x(1)/(1+eps*x(1)))*x(2)*T0/x(3)*R/x(3);
ra =(k2*(1/3*Pa*Ka^2-Pc^2/Pb)/(1+k3*Pc/Pb^1.5)^0.75);
dxdw=[-ra/ndA0;
    -alpha*x(3)/T0/(2*x(2))*(1+eps*x(1));
    ra*Lambda/(ndA0*CpA+1/3*ndB0*CpB)];
end

%% Declaration of the variables
P0    = 150;                % Inlet pressure, atm
alpha = 0.6004;             % Pressure drop parameter, 1/(kg catalyst)
k2=4.6*10^8;
Ka= 0.010985321;
k3=2.96*10^-6;
Lambda=-109;
R=8.314*1000;
CpA= 0.02919;
CpB=0.02922;
eps=-0.5;
T0=703;
ndA0=1798.94;
ndB0=3*ndA0;
v0=1000;
%% Declaration of the weight span and the initial conditions
W  = [0 60];  % Catalyst weight span, kg catalyst
xi = [0 1 703];   % xi(1)=eps at the inlet, xi(2)=P/P0 at inlet xi(3)=Ti

%% Solution of ODE system
[W,x]=ode113(@(w,x)funisopbr(w,x,k2, Ka, k3, ndA0, ndB0, alpha, R, Lambda, CpA, CpB, T0, eps, v0),W, xi );

%% Extraction and plotting of the results
% Calculation of molar flowrates
ndA = ndA0-ndA0*x(:,1);    % Molar flowrate of A
ndB = ndB0-ndA0*x(:,1);  % Molar flowrate of B
ndC = ndA0*x(:,1);         % Molar flowrate of C

% Plot of molar flowrates vs. weight of catalyst
figure(1)
plot(W,ndA,'r', W,ndB,'k',W,ndC, 'b','linewidth', 2);
title('Molar flowrates vs. weight of catalyst');
xlabel('Weight of catalyst, kg');
ylabel('Molar flowrate per tube, mol/s');
legend('nd_A','nd_B','nd_C');

% Plot of pressure vs. weight of catalyst
figure(2)
plot(W,x(:,2)*10,'m','linewidth',2);
title('Pressure vs. weight of catalyst');
xlabel('Weight of catalyst, kg');
ylabel('Pressure, atm');
 

% Plot of temperature vs. volume
figure(2)
plot(W,x(:,3),'linewidth',2);
title('Temperature vs. volume');
xlabel('Volume, L');
ylabel ('Temperature, K');

Since I didn’t wait for it to finish, I have no idea if there are other problems, so I’m not posting this as an Answer.

[Muge Savas]

Thank you for your time anyway.

[Star Strider]

My pleasure.

Answer 1

[Star Strider]

Your ODE system is apparently ‘stiff’. If you change the solver to ode15s:

%% Declaration of the weight span and the initial conditions
W  = [0 60]+eps;  % Catalyst weight span, kg catalyst
xi = [0 1 703];   % xi(1)=eps at the inlet, xi(2)=P/P0 at inlet xi(3)=Ti
%% Solution of ODE system
[W,x]=ode15s(@(w,x)funisopbr(w,x,k2, Ka, k3, ndA0, ndB0, alpha, R, Lambda, CpA, CpB, T0, eps, v0),W, xi );

and use semilogx for the plots:

figure(1)
semilogx(W,ndA,'r', W,ndB,'k',W,ndC, 'b','linewidth', 2);
title('Molar flowrates vs. weight of catalyst');
xlabel('Weight of catalyst, kg');
ylabel('Molar flowrate per tube, mol/s');
legend('nd_A','nd_B','nd_C')
grid
% Plot of pressure vs. weight of catalyst
figure(2)
semilogx(W,x(:,2)*10,'m','linewidth',2);
title('Pressure vs. weight of catalyst');
xlabel('Weight of catalyst, kg');
ylabel('Pressure, atm')
grid
 
% Plot of temperature vs. volume
figure(3)
semilogx(W,x(:,3),'linewidth',2);
title('Temperature vs. volume');
xlabel('Volume, L');
ylabel ('Temperature, K')
grid

it solves in a few seconds (on my computer), and produces decent plots.

My apologies for not thinking of thie earlier.