As you can see above, the vector C of initial conditions is 4x2, but the time derivative vector dCsdt in function "concentrations" is 4x1.
Have successfully used ode45 once in this code, but with the same size parameters am getting an error about initial conditions vector being longer than the return vector.
4 views (last 30 days)
Show older comments
%% Problem 1
% Constants throughout the calculation
l = 2/100; %[m]
d = 1.1/1000; %[m]
e = 0.65; %[-]
u = 0.21; %[m/s]
Sh = 3.0; %[-]
p = 0.98; %[Pa]
Ru = 8.3145; %[J/mol*K]
Ga = (4*e)/d;
rho = 31.0175; %[kg/m3]
T = 380;
Cn2 = 27.7606;
Cco = 0.1551;
Co2 = 1.5509;
Cco2 = 1.5509;
% Diffusion Rates
Dn2 = 0.220/100/100; %[m2/s]
Dco = 0.212/100/100; %[m2/s]
Do2 = 0.176/100/100; %[m2/s]
Dco2 = 0.165/100/100; %[m2/s]
% Arrhenius Variables
k1 = 1E13;
E1 = 68E3;
Aco = (2E-2);
dHco = -50E3;
Ao2 = (5E-11);
dHo2 = -60E3;
% Equilibrium Constants
Kco = Aco*exp(-dHco/(Ru*T));
Ko2 = Ao2*exp(-dHo2/(Ru*T));
K = [Kco Ko2];
% Mass Transfer to the Surface
kn2 = (Sh*Dn2)/d;
kco = (Sh*Dco)/d;
ko2 = (Sh*Do2)/d;
kco2 = (Sh*Dco2)/d;
k = [kn2 kco ko2 kco2];
% Number of points
N = 30;
% Length of Catalyst
l = 2/100; %[m]
% Length of Segments
dx = (N-1)/l;
% Ideal Gas Constant
Ru = 8.3145; %[J/mol*K]
% Temperatures
T = 380; %[K]
% Pressure
p = 0.98*100000; %[Pa]
% Inlet Mole Fractions
Xn2 = 0.895;
Xco = 0.005;
Xo2 = 0.050;
Xco2 = 0.050;
% Density
rho = p./(Ru*(T));
% Inlet Partial Pressures
ppn2 = Xn2*p;
ppco = Xco*p;
ppo2 = Xo2*p;
ppco2 = Xco2*p;
pp = [ppn2; ppco; ppo2; ppco2];
dt = 0.01;
% Step 1: All initial concentrations
C = [Xn2; Xco; Xo2; Xco2]*rho;
Cn2 = C(1);
Cco = C(2);
Co2 = C(3);
Cco2 = C(4);
Cs = [0; 0; 0; 0];
tspan = linspace(0,1,2);
% Step 2: Solve species on the surface at the boundary
[t,Cs] = ode45(@(t,C) concentrations(t,C,pp,k1,E1,Ru,T,K,k,Ga,e),tspan,C);
Csout = [Cs(1,1); Cs(end,2); Cs(end,3); Cs(end,4)];
% Step 3a, 3b, 3c: Propagate the bulk gas species (i+1) by solving for the
% next node via given equations solved for euler implicit
Csout(:,2) = Csout(:,1);
Xn2 = Csout(1,2)/(Csout(1,2)+Csout(2,2)+Csout(3,2)+Csout(4,2));
Xco = Csout(2,2)/(Csout(1,2)+Csout(2,2)+Csout(3,2)+Csout(4,2));
Xo2 = Csout(3,2)/(Csout(1,2)+Csout(2,2)+Csout(3,2)+Csout(4,2));
Xco2 = Csout(4,2)/(Csout(1,2)+Csout(2,2)+Csout(3,2)+Csout(4,2));
pp = [Xn2*p Xco*p Xo2*p Xco2*p];
for j = 1:4
C(j,2) = (((dx*k(j)*Ga*Csout(j,2))/(u*e))+C(j,1))/(1+((dx...
*k(j)*Ga)/(u*e)));
end
size(C)
[t,Cs] = ode45(@(t,C) concentrations(t,C,pp,k1,E1,Ru,T,K,k,Ga,e),tspan,C);
% Csout(:,2) = [Cs(1,1); Cs(end,2); Cs(end,3); Cs(end,4)];
% for i = 1:N
% Rn2 = 0;
% Rco = (k1*exp(-E1/(Ru*T))*Kco*Ko2*ppco*(ppo2^(1/2)))/((1+(Kco*ppco)+(Ko2*(ppo2^1/2)))^2);
% Ro2 = (1/2)*Rco;
% Rco2 = -Rco;
% R(:,i) = [Rn2; Rco; Ro2; Rco2];
% % Setting up each node propagation
% C(:,i+1) = C(:,i);
% Cs(:,i+1) = Cs(:,i);
% Cg(:,i) = [0 0 0 0];
% Cgs(:,i) = [0 0 0 0];
% while abs(C(2,i+1)-Cg(2,i))>=1E-5 && abs(Cs(2,i+1)-Cgs(2,i))>=1E-5
%
% Cg(:,i) = C(:,i+1);
% Cgs(:,i) = Cs(:,i+1);
% for j = 1:4
% C(j,i+1) = (((dx*k(j)*Ga*Cgs(j,i))/(u*e))+C(j,i))/(1+((dx...
% *k(j)*Ga)/(u*e)));
% end
% [t,Cs] = ode45(@(t,C) concentrations(t,C,pp,k1,E1,Ru,T,K,k,Ga,e),tspan,C);
% Cs = [Cs(1,1); Cs(end,2); Cs(end,3); Cs(end,4)];
%
% % Mole Fractions
% Xn2 = Cs(1,i+1)/(Cs(1,i+1) + Cs(2,i+1) + Cs(3,i+1) + Cs(4,i+1));
% Xco = Cs(2,i+1)/(Cs(1,i+1) + Cs(2,i+1) + Cs(3,i+1) + Cs(4,i+1));
% Xo2 = Cs(3,i+1)/(Cs(1,i+1) + Cs(2,i+1) + Cs(3,i+1) + Cs(4,i+1));
% Xco2 = Cs(4,i+1)/(Cs(1,i+1) + Cs(2,i+1) + Cs(3,i+1) + Cs(4,i+1));
% % Partial Pressures
% ppn2 = p*Xn2;
% ppco = p*Xco;
% ppo2 = p*Xo2;
% ppco2 = p*Xco2;
%
% end
% end
%% Functions Used (Write first w/ initial conditions)
function dCsdt = concentrations(t,C,pp,k1,E1,Ru,T,K,k,Ga,e)
% Equilibrium Constants
Kco = K(1);
Ko2 = K(2);
% Mass Transfer to the Surface
kn2 = k(1);
kco = k(2);
ko2 = k(3);
kco2 = k(4);
% Surface Concentration Guesses
Csn2 = C(1);
Csco = C(2);
Cso2 = C(3);
Csco2 = C(4);
% Bulk Gas Concentrations
Cn2 = C(1);
Cco = C(2);
Co2 = C(3);
Cco2 = C(4);
% Partial Pressure Updates
ppn2 = pp(1);
ppco = pp(2);
ppo2 = pp(3);
ppco2 = pp(4);
% Reaction rate expressions to update each iteration
Rn2 = 0;
Rco = (k1*exp(-E1/(Ru*T))*Kco*Ko2*ppco*(ppo2^(1/2)))/((1+(Kco*ppco)+(Ko2*(ppo2^1/2)))^2);
Ro2 = (1/2)*Rco;
Rco2 = -Rco;
% Surface Concentrations of the Species
dCsdt = [(kn2*Ga*(Cn2-Csn2)-Rn2)/(1-e);(kco*Ga*(Cco-Csco)-Rco)/(1-e);...
(ko2*Ga*(Co2-Cso2)-Ro2)/(1-e);(kco2*Ga*(Cco2-Csco2)-Rco2)/(1-e)];
size(dCsdt)
end
3 Comments
Torsten
on 11 Dec 2024
Edited: Torsten
on 11 Dec 2024
I don't know what you do in your code, but if you solve for n concentrations C1,...,Cn, you need n time derivatives dC1/dt,...,dCn/dt that reflect their temporal evolution.
The initial values for C1,...,Cn at time t=0 have to be passed in your C-vector, the time derivatives dC1/dt,...,dCn/dt have to be passed in dCsdt.
The vector C passed to "concentrations" is the vector of concentrations at time t and its size is equal to the size of the vector of initial conditions.
function dCsdt = concentrations(t,C,pp,k1,E1,Ru,T,K,k,Ga,e)
Answers (1)
Sam Chak
on 11 Dec 2024
Hi @Raoul
There is no error in organizing the code structure clearly with constants and states (variables that change over time). You should avoid unnecessary duplication by passing confusing vectors around. The number of states should be consistent with the number of ordinary differential equations (ODEs) represented by dCsdt. However, it seems that you intend to update some 'constant variables.' If so, you will need to modify the code to include the necessary ODEs that describe the rate of change of those 'constant variables.'
In the simulation below, it appears that the 'states' are in equilibrium. Please verify this as I'm not as expert in concentration dynamics.
tspan = linspace(0, 100, 10001);
Xn2 = 0.895;
Xco = 0.005;
Xo2 = 0.050;
Xco2 = 0.050;
rho = 31.0175; % [kg/m3]
C0 = [Xn2; Xco; Xo2; Xco2]*rho; % Initial states are re-labeled as C0
[t, C] = ode45(@concentrations, tspan, C0);
plot(t, C), grid on
function dCsdt = concentrations(t, C)
% Constants throughout the calculation
l = 2/100; % [m]
d = 1.1/1000; % [m]
e = 0.65; % [-]
u = 0.21; % [m/s]
Sh = 3.0; % [-]
p = 0.98; % [Pa]
Ru = 8.3145; % [J/mol*K]
Ga = (4*e)/d;
rho = 31.0175; % [kg/m3]
T = 380;
% Bulk Gas Concentrations
Cn2 = 27.7606;
Cco = 0.1551;
Co2 = 1.5509;
Cco2 = 1.5509;
% % Bulk Gas Concentrations % Unnecessary to pass fixed constants
% Cn2 = C(1);
% Cco = C(2);
% Co2 = C(3);
% Cco2 = C(4);
% Diffusion Rates
Dn2 = 0.220/100/100; %[m2/s]
Dco = 0.212/100/100; %[m2/s]
Do2 = 0.176/100/100; %[m2/s]
Dco2 = 0.165/100/100; %[m2/s]
% Arrhenius Variables
k1 = 1E13;
E1 = 68E3;
Aco = (2E-2);
dHco = -50E3;
Ao2 = (5E-11);
dHo2 = -60E3;
% Equilibrium Constants
Kco = Aco*exp(-dHco/(Ru*T));
Ko2 = Ao2*exp(-dHo2/(Ru*T));
% K = [Kco Ko2]; % Unnecessary to pass fixed constants
% % Equilibrium Constants, K
% Kco = K(1);
% Ko2 = K(2);
% Mass Transfer to the Surface
kn2 = (Sh*Dn2)/d;
kco = (Sh*Dco)/d;
ko2 = (Sh*Do2)/d;
kco2 = (Sh*Dco2)/d;
% k = [kn2 kco ko2 kco2]; % Unnecessary to pass fixed constants
% % Mass Transfer to the Surface, k
% kn2 = k(1);
% kco = k(2);
% ko2 = k(3);
% kco2 = k(4);
% Number of points
N = 30;
% Length of Catalyst
l = 2/100; %[m]
% Length of Segments
dx = (N - 1)/l;
% Ideal Gas Constant
Ru = 8.3145; %[J/mol*K]
% Temperatures
T = 380; %[K]
% Pressure
p = 0.98*100000; %[Pa]
% Inlet Mole Fractions
Xn2 = 0.895;
Xco = 0.005;
Xo2 = 0.050;
Xco2 = 0.050;
% Density
rho = p./(Ru*(T));
% Inlet Partial Pressures
ppn2 = Xn2 *p;
ppco = Xco *p;
ppo2 = Xo2 *p;
ppco2 = Xco2*p;
% pp = [ppn2; ppco; ppo2; ppco2]; % Unnecessary to pass fixed constants
% % Partial Pressure Updates, pp
% ppn2 = pp(1);
% ppco = pp(2);
% ppo2 = pp(3);
% ppco2 = pp(4);
% dt = 0.01;
% Step 1: All initial concentrations
% C = [Xn2; Xco; Xo2; Xco2]*rho; % "C" is reserved for the "ODE states"
% Cn2 = C(1);
% Cco = C(2);
% Co2 = C(3);
% Cco2 = C(4);
% Cs = [0; 0; 0; 0];
% Reaction rate expressions to update each iteration
Rn2 = 0;
Rco = (k1*exp(-E1/(Ru*T))*Kco*Ko2*ppco*(ppo2^(1/2)))/((1+(Kco*ppco)+(Ko2*(ppo2^1/2)))^2);
Ro2 = (1/2)*Rco;
Rco2 = -Rco;
% Surface Concentration Guesses, C
Csn2 = C(1);
Csco = C(2);
Cso2 = C(3);
Csco2 = C(4);
% Surface Concentrations of the Species
dCsdt = [(kn2 *Ga*(Cn2 - Csn2) - Rn2) /(1 - e);
(kco *Ga*(Cco - Csco) - Rco) /(1 - e);
(ko2 *Ga*(Co2 - Cso2) - Ro2) /(1 - e);
(kco2*Ga*(Cco2 - Csco2) - Rco2)/(1 - e)];
% size(dCsdt)
end
0 Comments
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)