Help with odes and two different reactions

10 次查看（过去 30 天）

显示更早的评论

Tom Goodland 2022-1-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1628960-help-with-odes-and-two-different-reactions

评论： Torsten 2022-1-18

𝑑𝐶𝐴/𝑑𝑡 = 𝑟1𝐴 = −𝑘𝐴𝐶𝐴𝐶𝐵 ( 9 ) The first two have the same rate constant

𝑑𝐶𝐵/𝑑𝑡 = 𝑟1𝐵 = 𝑟1𝐴 = −𝑘𝐴𝐶𝐴𝐶𝐵 ( 10 )

d𝐶𝑆/𝑑𝑡 = 𝑟2𝑆 = −𝑘𝑆𝐶𝑆 ( 11 )

A + B → C + ½ D (gas)

S → 3 D (gas) + misc liquids and solids

These are the two reactions, do I have to do two seperate odes because of the different reactions (and rate constants)? Or is there a way to work out all the differentials in the same ode with initial concentrations known?

If any of my code is needed to answer the question let me know.

Thanks

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Torsten 2022-1-14

Look at the second example (Solve Stiff ODE) to see how you can solve all three ODEs in one call to the integrator.

https://de.mathworks.com/help/matlab/ref/ode15s.html

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

HWIK 2022-1-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1628960-help-with-odes-and-two-different-reactions#answer_874755

在 MATLAB Online 中打开

Just write a function for the ODEs and feed it to the solver:

function out = name_your_function(tspan, inputs)
%Get your inputs and assign them to variables:
%(keep the same order as in your initial conditions!)
ca = inputs(1);
cb = inputs(2);
cs = inputs(3);
%Define your relationships:
%choose your values of k1 & k2:
k1 = ;
k2 = ;
%Reaction rates:
r1 = k1*ca*cb;
r2 = k2*cs;
%Define your ODEs:
dca = -r1;
dcb = -r1;
dcs = -r2;
%Assign your ourput:
out = [dca; dcb; dcs];
end

Solve over a given time span and for whatever initial conditions you have, using whatever solver seems to suit your system best.

33 个评论
显示 31更早的评论隐藏 31更早的评论

Tom Goodland 2022-1-15

I've tried to add the head space pressure and temperature to the ode so I've added stuff to both the function and script. however I've been getting errors:

tspan = [0 200]; %just as an example

initial_conditions = [4.3 5.1 3 0 4.4 422]; %its the values of the variables you are solving for at t=0

%there must be as many inputs to the initial conditions as variables solved

%Solve your ODE function (I picked ode45 in this case,could be another):

[t,y] = ode15s(@(t,y)thorsten(t,y(1),y(2),y(3),y(4),y(5),y(6)), tspan, initial_conditions);

%output t is your time range, and y an array with all your variables where:

ca = y(:,1);

cb = y(:,2);

cs = y(:,3);

cd = y(:,4);

P = y(:,5);

T_2 = y(:,6);

plot(t,[ca,cb,cs,cd,P,T_2]);

Here's the script:

function out = thorsten(tspan,ca,cb,cs,cd,P,T_2)

%Get your inputs and assign them to variables:

%(keep the same order as in your initial conditions!)

%Define your relationships:

k_0_1 = 4E14;

k_0_2 = 1E84;

Ea1 = 128000;

Ea2 = 800000;

R_gas = 8.314;

V0 = 4000;% L

VH = 5000;

Cv1 = 3360; % mol / h * atm

Cv2 = 53600; % mol / h * atm

Hr1 = -45400; % J / mol

Hr2 = -3.2E5;

T_2 = 422; % K

Ta = 373;

UA = 0;

NCp = 1.26E7; % J / K

%choose your values of k1 & k2:

k1 = k_0_1*exp(-Ea1/(R_gas*T_2));

k2 = k_0_2*exp(-Ea2/(R_gas*T_2));

%Reaction rates:

r1 = k1*ca*cb;

r2 = k2*cs;

%Fd

Fd = ((-0.5*r1) - (3*r2))*V0;

% Fvent

Fvent = Fd;

if Fd < 11,400; % mol / h

Fvent = (P - 1)*Cv1;

if P < 28.2 % atm

end

Fvent = (P - 1)*(Cv1 + Cv2);

if P > 28.2

end

%Define your ODEs:

dca = -r1;

dcb = -r1;

dcs = -r2;

dcd = r1/2 + 3*r2;

dp = (Fd - Fvent)*((R_gas*T_2)/VH);

dT = (V0*(r1*Hr1 + r2*Hr2) - UA(T_2 - Ta))/NCp;

%Assign your ourput:

out = [dca; dcb; dcs; dcd,dp,dT];

end

I've named the temperature T_2 because I've used T previously in different work. dT is representing the change in T_2.

The errors im getting are:

Index exceeds the number of array elements. Index must not exceed 1.

Error in thorsten (line 43)

dT = (V0*(r1*Hr1 + r2*Hr2) - UA(T_2 - Ta))/NCp;

Error in untitled2>@(t,y)thorsten(t,y(1),y(2),y(3),y(4),y(5),y(6)) (line 5)

[t,y] = ode15s(@(t,y)thorsten(t,y(1),y(2),y(3),y(4),y(5),y(6)), tspan, initial_conditions);

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode15s (line 152)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in untitled2 (line 5)

[t,y] = ode15s(@(t,y)thorsten(t,y(1),y(2),y(3),y(4),y(5),y(6)), tspan, initial_conditions);

Do you know why im getting these errors, any help would be appreciated.

Thanks

Torsten 2022-1-15

在 MATLAB Online 中打开

function main
tspan = [0 200]; %just as an example
initial_conditions = [4.3 5.1 3 0 4.4 422]; %its the values of the variables you are solving for at t=0
%there must be as many inputs to the initial conditions as variables solved
%Solve your ODE function (I picked ode45 in this case,could be another):
[t,y] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6)), tspan, initial_conditions);
%output t is your time range, and y an array with all your variables where:
ca = y(:,1);
cb = y(:,2);
cs = y(:,3);
cd = y(:,4);
P = y(:,5);
T_2 = y(:,6);
figure(1)
plot(t,[ca,cb,cs,cd])
figure(2)
plot(t,P)
figure(3)
plot(t,T_2)
end
function out = fun(tspan,ca,cb,cs,cd,P,T_2)
%Get your inputs and assign them to variables:
%(keep the same order as in your initial conditions!)
%Define your relationships:      
k_0_1 = 4E14;
k_0_2 = 1E84;
Ea1 = 128000;
Ea2 = 800000;
R_gas = 8.314;
V0 = 4000;% L
VH = 5000;
Cv1 = 3360; % mol / h * atm
Cv2 = 53600; % mol / h * atm
Hr1 = -45400; % J / mol
Hr2 = -3.2E5;
Ta = 373;
UA = 0;
NCp = 1.26E7; % J / K
%choose your values of k1 & k2:
k1 = k_0_1*exp(-Ea1/(R_gas*T_2));
k2 = k_0_2*exp(-Ea2/(R_gas*T_2));
%Reaction rates:
r1 = k1*ca*cb;
r2 = k2*cs;
%Fd
Fd = (-0.5*r1-3*r2)*V0;
%
if Fd < 11.4 % mol / h
    Fvent = Fd;
else 
    if P < 28.2 % atm
        Fvent = (P - 1)*Cv1;
    else
        Fvent = (P - 1)*(Cv1 + Cv2);
    end
end
%Define your ODEs:
dca = -r1;
dcb = -r1;
dcs = -r2;
dcd = r1/2 + 3*r2;
dp = (Fd - Fvent)*((R_gas*T_2)/VH);
dT = (V0*(r1*Hr1 + r2*Hr2) - UA*(T_2 - Ta))/NCp;
%Assign your ourput:
out = [dca; dcb; dcs; dcd;dp;dT];
end

Check whether Fvent is calculated correctly.

Torsten 2022-1-17

在 MATLAB Online 中打开

This is a rough guideline.

If errors pop up, first make an attempt to correct them on your own, please.

function main
  tspan = [0 200]; %just as an example
  initial_conditions = [4.3 5.1 3 0 4.4 422]; %its the values of the variables you are solving for at t=0
  options = odeset('Events',@myEventsFcn);
  iflag = 1;
  [t1,y1,te,ye,ie] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),iflag), tspan, initial_conditions,options);
  iflag = 2;
  [t2,y2] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),iflag), [te tspan(end)], ye);
  t = [t1;te;t2];
  y = [y1;ye;y2];
  %output t is your time range, and y an array with all your variables where:
  ca = y(:,1);
  cb = y(:,2);
  cs = y(:,3);
  cd = y(:,4);
  P = y(:,5);
  T_2 = y(:,6);
  figure(1)
  plot(t,[ca,cb,cs,cd])
  figure(2)
  plot(t,P)
  figure(3)
  plot(t,T_2)
end
function [value,isterminal,direction] = myEventsFcn(t,y)
  value = y(6) - 455; % The value that we want to be zero
  isterminal = 1;  % Halt integration 
  direction = 0;   % The zero can be approached from either direction
end
function out = fun(tspan,ca,cb,cs,cd,P,T_2,iflag)
%Get your inputs and assign them to variables:
%(keep the same order as in your initial conditions!)
%Define your relationships:      
k_0_1 = 4E14;
k_0_2 = 1E84;
Ea1 = 128000;
Ea2 = 800000;
R_gas = 8.314;
V0 = 4000;% L
VH = 5000;
Cv1 = 3360; % mol / h * atm
Cv2 = 53600; % mol / h * atm
Hr1 = -45400; % J / mol
Hr2 = -3.2E5;
Ta = 373;
if iflag == 1
  UA = 0;
elseif iflag == 2
  UA = 2.77E6;
end
NCp = 1.26E7; % J / K
%choose your values of k1 & k2:
k1 = k_0_1*exp(-Ea1/(R_gas*T_2));
k2 = k_0_2*exp(-Ea2/(R_gas*T_2));
%Reaction rates:
r1 = k1*ca*cb;
r2 = k2*cs;
%Fd
Fd = (-0.5*r1-3*r2)*V0;
%
if Fd < 11.4 % mol / h
    Fvent = Fd;
else 
    if P < 28.2 % atm
        Fvent = (P - 1)*Cv1;
    else
        Fvent = (P - 1)*(Cv1 + Cv2);
    end
end
%Define your ODEs:
dca = -r1;
dcb = -r1;
dcs = -r2;
dcd = r1/2 + 3*r2;
dp = (Fd - Fvent)*((R_gas*T_2)/VH);
dT = (V0*(r1*Hr1 + r2*Hr2) - UA*(T_2 - Ta))/NCp;
%Assign your ourput:
out = [dca; dcb; dcs; dcd;dp;dT];
end

Tom Goodland 2022-1-17

I am thinking some type of loop would be best for this challenge but I'm not sure what, does anyone know what loop code would be best for the problem I described in the 2 previous comments on this thread?

Thanks

Torsten 2022-1-18

在 MATLAB Online 中打开

You mean something like this ?

function main
tspan = [0 10]; %just as an example
initial_conditions = [4.3 5.1 3 0 4.4 422]; %its the values of the variables you are solving for at t=0
%there must be as many inputs to the initial conditions as variables solved
%Solve your ODE function (I picked ode45 in this case,could be another):
Tstop_left = 450.0;
Tstop_right = 495.0;
Tstop = (Tstop_left + Tstop_right)/2.0;
while Tstop_right - Tstop_left >= 0.1
    options = odeset('RelTol',1e-8,'AbsTol',1e-8,'Events',@(t,y)myEventsFcn(t,y,Tstop));
    iflag = 1;
    [t1,y1,te,ye,ie] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),iflag), tspan, initial_conditions,options);
    iflag = 2;
    options = odeset('RelTol',1e-8,'AbsTol',1e-8,'InitialStep',1e-8);
    [t2,y2] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),iflag), [te tspan(end)], ye, options);
    Tmax = max(y2(:,6))  
    if Tmax > 500
        Tstop_left = Tstop_left;
        Tstop_right = Tstop;
        Tstop = (Tstop_left + Tstop_right)/2;
    else
        Tstop_left = Tstop;
        Tstop_right = Tstop_right;
        Tstop = (Tstop_left + Tstop_right)/2;
    end 
end
Tstop
t = [t1;te;t2];
y = [y1;ye;y2];
%output t is your time range, and y an array with all your variables where:
ca = y(:,1);
cb = y(:,2);
cs = y(:,3);
cd = y(:,4);
P = y(:,5);
T_2 = y(:,6);
figure(1)
plot(t,[ca,cb,cs,cd])
figure(2)
plot(t,P)
figure(3)
plot(t,T_2)
end
function [value,isterminal,direction] = myEventsFcn(t,y,Tstop)
value = y(6) - Tstop; % The value that we want to be zero
isterminal = 1;  % Halt integration 
direction = 0;   % The zero can be approached from either direction
end
function out = fun(tspan,ca,cb,cs,cd,P,T_2,iflag)
%Get your inputs and assign them to variables:
%(keep the same order as in your initial conditions!)
%Define your relationships:      
k_0_1 = 4E14;
k_0_2 = 1E84;
Ea1 = 128000;
Ea2 = 800000;
R_gas = 8.314;
V0 = 4000;% L
VH = 5000;
Cv1 = 3360; % mol / h * atm
Cv2 = 53600; % mol / h * atm
Hr1 = -45400; % J / mol
Hr2 = -3.2E5;
Ta = 373;
if iflag == 1
    UA = 0;
elseif iflag == 2
    UA = 2.77E6;
end
NCp = 1.26E7; % J / K
%choose your values of k1 & k2:
k1 = k_0_1*exp(-Ea1/(R_gas*T_2));
k2 = k_0_2*exp(-Ea2/(R_gas*T_2));
%Reaction rates:
r1 = k1*ca*cb;
r2 = k2*cs;
%Fd
Fd = (-0.5*r1-3*r2)*V0;
%
if Fd < 11.4 % mol / h
    Fvent = Fd;
else 
    if P < 28.2 % atm
        Fvent = (P - 1)*Cv1;
    else
        Fvent = (P - 1)*(Cv1 + Cv2);
    end
end
%Define your ODEs:
dca = -r1;
dcb = -r1;
dcs = -r2;
dcd = r1/2 + 3*r2;
dp = (Fd - Fvent)*((R_gas*T_2)/VH);
dT = (V0*(r1*-Hr1 + r2*-Hr2) - UA*(T_2 - Ta))/NCp;
%Assign your ourput:
out = [dca; dcb; dcs; dcd;dp;dT];
end