Error using Odearguments - vector must return a column
Show older comments
I kept on trying to reword and rearrange my code however it asked me to do it so that the line return a column vector. I tried reshape(), [],1, (does not work). Also tried .' at the end which also doesnt work. I'm not sure where I can continue with this?
%r_H2= -(-Rds - Rwf - Rh);
%Rd = kds*(CH20*((theta_O2-X)/(1+eps*X)))*(CH20*((1-X)/(1+eps*X)))^2
%Rwf = kwf* (CH20*((theta_O2-(b/a)*X)/(1+wf_eps*X)))^.5*CH20.*((1-X)./(1+wf_eps*X))
%Rh = kh* (CH2O2*(theta_H2-(d_b/d_a)*X)) *CH2O2*(1-X)
wspan = [0 1000];
dxdw = @(W,X)(kds.*(CH20.*((theta_O2-X)./(1+eps.*X))).*(CH20.*((1-X)./(1+eps.*X))).^2 + ...
kwf.* (CH20.*((theta_O2-(b/a).*X)./(1+wf_eps.*X))).^.5*CH20.*((1-X)./(1+wf_eps.*X))+ kh.* (CH2O2.*(theta_H2-(d_b/d_a).*X)) .*CH2O2.*(1-X));
IC = [0 0 0 0 0 0 0 0 0 0 0];
[W,X] = ode45(dxdw,wspan,IC);
plot(W,X)
title ('Catalyst v.s Conversion')
xlabel('Weight of Catalyst (kg)')
ylabel('Conversion of H2')
Accepted Answer
You only have one differential equation, so only one initial condition is necessary.
%r_H2= -(-Rds - Rwf - Rh);
%Rd = kds*(CH20*((theta_O2-X)/(1+eps*X)))*(CH20*((1-X)/(1+eps*X)))^2
%Rwf = kwf* (CH20*((theta_O2-(b/a)*X)/(1+wf_eps*X)))^.5*CH20.*((1-X)./(1+wf_eps*X))
%Rh = kh* (CH2O2*(theta_H2-(d_b/d_a)*X)) *CH2O2*(1-X)
wspan = [0 1000];
dxdw = @(W,X)(kds.*(CH20.*((theta_O2-X)./(1+eps.*X))).*(CH20.*((1-X)./(1+eps.*X))).^2 + ...
kwf.* (CH20.*((theta_O2-(b/a).*X)./(1+wf_eps.*X))).^.5*CH20.*((1-X)./(1+wf_eps.*X))+ kh.* (CH2O2.*(theta_H2-(d_b/d_a).*X)) .*CH2O2.*(1-X));
% IC = [0 0 0 0 0 0 0 0 0 0 0];
IC = eps;
[W,X] = ode45(dxdw,wspan,IC);
plot(W,X)
title ('Catalyst v.s Conversion')
xlabel('Weight of Catalyst (kg)')
ylabel('Conversion of H2')
All the missing constants prevent me from running this.
.
4 Comments
Anh Bui
on 17 Apr 2024
Your ‘dxdw’ function produces an (11x11) matrix of identical (matching) columns in each iteration. That was the reason for the ‘column vector’ error. I wrote a ‘helper’ function, ‘dxdwc’ that returns only the first column of the matrix, and changed the solver to ode15s because the system takes a while to integrate, and like every kinetic equation I’ve ever encountered, is ‘stiff’ because of the wide range of parameter magnitudes, so this is appropriate.
There may still be problems. (I have no idea what the result is supposed to look like, however it makes sense that the more catalyst added the more quickly the system reaches equilibrium.)
Try this —
%% Temperature range
T = 313; % Kelvin
%T = 40+273; %K
%% Pressure
% P = [20 30 50]; % Bar
% P = P.*(10^5); % Pa
P = 5.1*(10^6); %Pa
%% Constants (A, Ea, and R)
A_ds = 4.5*(10^7);
A_wf = 6.5*(10^7);
A_d = 1.005*(10^7);
A_h = 280000;
Ea_ds = 6.9*(10^3); % J/mol
Ea_wf = 20*(10^3); % J/mol
Ea_d = 53*(10^3); % J/mol
Ea_h = 44.3*(10^3); % J/mol
R = 8.314; % J/molK or Pa*m3/molK
%% Solve for kinetic constants (ki)
kds = A_ds.*exp((-Ea_ds)./(R.*T));
kwf = A_wf.*exp((-Ea_wf)./(R.*T));
kd = A_d.*exp((-Ea_d)./(R.*T));
kh = A_h.*exp((-Ea_h)./(R.*T));
%initial molar flow mol/s
F_H2_in = 15;
F_O2_in = 87.85;
F_CO2_in = 278.82;
yH2o = .04;
theta_O2= F_O2_in/F_H2_in;
theta_CO2 = F_CO2_in/F_H2_in;
theta_H2O2 = 0;
delta = -1;
eps = yH2o*delta;
%Temperature change
T0T = 1;
%Pressure Drop
PP0 = 1;
%initial mass flow g/s
M_H2 = F_H2_in/1.02;
M_O2 = F_O2_in/16;
M_CO2 = F_CO2_in/44;
%density of gas mixtures g/L
rho_mix = 77.555;
%volumetric flow L/s
v_H2 = M_H2/rho_mix;
v_O2 = M_O2/rho_mix;
v_CO2 = M_CO2/rho_mix;
v0 = v_CO2+v_O2+v_H2;
%initial concentration of H2 mol/L
CH20 = F_H2_in/v0;
%mass balance for direct synthesis
X = [0:0.1:1];
CH2 = CH20.*((1-X)./(1+eps.*X))*T0T*PP0;
CO2 = CH20.*((theta_O2-X)./(1+eps*X))*T0T*PP0;
CH2O2 = CH20.*((theta_H2O2+X)./(1+eps*X))*T0T*PP0;
Rds = kds.*CO2.*CH2.^2;
plot (X,CH2,X,CO2,X,CH2O2)
legend('C_H2','C_O2','C_H2O2',Location='east')
title('Direct Synthesis Conversion')
xlabel('Conversion')
ylabel('Concentration (mol/L)')
%rate of direct synthesis
%2H2+O2->2H2O
%coefficient of the chemical equation
a = 2; %H2 coeff
b = 1; %O2 coeff
c = 2;%water coeff
wf_yH20 = .04;
wf_delta = -a-(b/a)+(c/a);
wf_eps = wf_yH20*delta;
wf_C_H2 = CH20.*((1-X)./(1+wf_eps*X))*T0T*PP0;
wf_C_O2 = CH20.*((theta_O2-(b/a).*X)./(1+wf_eps.*X))*T0T*PP0;
wf_C_H2O = CH20.*((0+(c/a).*X)/(1+wf_eps.*X))*T0T*PP0;
wf_C_CO2 =((CH20*theta_CO2)./(1+eps*X))*T0T*PP0;
Rwf = kwf.*(wf_C_O2.^.5).*wf_C_H2;
plot (X,wf_C_H2,X,wf_C_CO2,X,wf_C_H2O,X,wf_C_CO2)
legend('C_H2','C_O2','C_H2O2','C_CO2',Location='east')
title('Water Formation Conversion')
xlabel('Conversion')
ylabel('Concentration (kmol/L)')
%H2O2(l)-> .5 O2 (g) + H2O(l)
%coeff of the chemical equation
d_a = 1;%H2O2 coeff
d_b = .5;%O2 coeff
d_c = 1; %water coeff
d_C_H2O2= CH2O2.*(1-X);
d_C_O2 = CH2O2.*(0+(d_b/d_a).*X);
d_C_H2O = CH2O2.*(0+(d_c/d_a).*X);
Rd= kd.*(d_C_H2O2);
plot (X,d_C_H2O2,X,d_C_O2,X,d_C_H2O)
legend('C_H2O2','C_O2','C_H2O',Location='northwest')
title('Decomposition Conversion')
xlabel('Conversion')
ylabel('Concentration (kmol/L)')
%H2O2(l)+H2(g)-> 2H2O(l)
%coeff of the chemical equation
h_a = 1;%coeff of h2o2
h_b = 1; %coeff of h2
h_c = 2; %coeff of water
%molar flow rate kmol/hr
F_H2O2 = 33.95;
F_H2 = 10.72;
h_yH20 = .04;
h_delta = 0;
h_eps = h_yH20*h_delta;
theta_H2 = F_H2/F_H2O2;
h_C_H2O2 = CH2O2.*(1-X);
h_C_H2 = CH2O2.*(theta_H2-(h_b/h_a).*X);
h_C_H2O = CH2O2.*(0+(h_c/h_a).*X);
Rh = kh.*h_C_H2.*h_C_H2O2;
plot (X,h_C_H2O2,'-o',X,h_C_H2,'ob',X,h_C_H2O,'-r')
legend('C_H2O2','C_H2','C_H2O',Location='northwest')
title('Hydrogenation Conversion')
xlabel('Conversion')
ylabel('Concentration (kmol/L)')
r_H2O2 = Rds-Rd-Rh;
r_H2O= Rwf+Rd+2.*Rh;
r_O2 = -(-Rds-0.5.*Rwf-0.5.*Rh);
r_H2= -(-Rds - Rwf - Rh);
plot(X,r_H2,'-o',X,r_H2O,'ob',X,r_O2,'b',X,r_H2O2,'r')
legend('Rate of H2','Rate of water','Rate of O2','Rate of H2O2')
title('Rate vs Conversion')
xlabel('Conversion')
ylabel('Rate of Components')
%r_H2= -(-Rds - Rwf - Rh);
%Rd = kds*(CH20*((theta_O2-X)/(1+eps*X)))*(CH20*((1-X)/(1+eps*X)))^2
%Rwf = kwf* (CH20*((theta_O2-(b/a)*X)/(1+wf_eps*X)))^.5*CH20.*((1-X)./(1+wf_eps*X))
%Rh = kh* (CH2O2*(theta_H2-(d_b/d_a)*X)) *CH2O2*(1-X)
wspan = [0 1000]/250;
% dxdw = @(W,X)(kds.*(CH20.*((theta_O2-X)./(1+eps.*X))).*(CH20.*((1-X)./(1+eps.*X))).^2 + kwf.* (CH20.*((theta_O2-(b/a).*X)./(1+wf_eps.*X))).^.5*CH20.*((1-X)./(1+wf_eps.*X))+ kh.* (CH2O2.*(theta_H2-(d_b/d_a).*X)) .*CH2O2.*(1-X));
IC = ones(size(X))*1E-8;
[W,X] = ode15s(@(W,X)dxdwc(W,X, kds,CH20,theta_O2,eps,kwf,b,a,wf_eps,kh,CH2O2,theta_H2,d_b,d_a),wspan,IC);
Out = [W X]
figure
plot(W,X)
Ax = gca;
Ax.XScale = 'log';
grid
title ('Catalyst v.s Conversion')
xlabel('Weight of Catalyst (kg)')
ylabel('Conversion of H2')
function d = dxdwc(W,X, kds,CH20,theta_O2,eps,kwf,b,a,wf_eps,kh,CH2O2,theta_H2,d_b,d_a)
dxdw = @(W,X)(kds.*(CH20.*((theta_O2-X)./(1+eps.*X))).*(CH20.*((1-X)./(1+eps.*X))).^2 + kwf.* (CH20.*((theta_O2-(b/a).*X)./(1+wf_eps.*X))).^.5*CH20.*((1-X)./(1+wf_eps.*X))+ kh.* (CH2O2.*(theta_H2-(d_b/d_a).*X)) .*CH2O2.*(1-X));
% Q = dxdw(W,X)
for k = 1:numel(X)
dv(k,:) = dxdw(W,X(k));
end
d = dv(:,1);
end
The ‘helper’ function ‘dxdwc’ extracts the first column from the matrix created by ‘dxdw’ (all the columns appear to be identical) and returns it to ode15s. It does not otherwise change the code in any way.
There may still be problems with the differential equation. I defer to you to address any problems that remain.
.
Anh Bui
on 17 Apr 2024
Star Strider
on 17 Apr 2024
As always, my pleasure!
More Answers (1)
This works - although I'm not sure if it is what you want:
wspan = [0 1000];
CH2O2 = @(X)CH20.*((theta_H2O2+X)./(1+eps*X))*T0T*PP0;
dxdw = @(W,X)(kds.*(CH20.*((theta_O2-X)./(1+eps.*X))).*(CH20.*((1-X)./(1+eps.*X))).^2 + kwf.* (CH20.*((theta_O2-(b/a).*X)./(1+wf_eps.*X))).^.5*CH20.*((1-X)./(1+wf_eps.*X))+ kh.* (CH2O2(X).*(theta_H2-(d_b/d_a).*X)) .*CH2O2(X).*(1-X));
IC = eps*ones(11,1);
[W,X] = ode15s(dxdw,wspan,IC);
And better use a variable name different from eps. It usually is a predefined constant is MATLAB:
eps
Categories
Find more on Chemistry in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
