Problem using ODE15s
显示 更早的评论
Hey everyone. I'm trying to model this chemical reaction:
This is a combination of equilibrium reactions (forward and backward reactions) and reactions in series.
I have written a code that does not work really well, the CO2 concentration seems ok but the other concentrations are just step functions and it is not supposed to look like that. My main guess is that something is wrong with the rates of reaction but I cannot see what is wrong. My second guess is that I have chosen wrong coding strategy.
clc
close all
clear all
c_0 = (46./1e6).*1000; %mol/m^3
H2CO3_0=0; %mol/m^3
HCO3_0=6437e-6*1000; %mol/m^3
H_0=0.00000001*1000; %mol/m^3
tspan= [0 1000];
u0=[c_0;H2CO3_0;HCO3_0;H_0]; %initial conditions for ODE's
%diff equations----------------------------
%u(1)=concentration CO2
%u(2) concentration H2CO3
%u(3) concentration HCO3-
%u(4) concentration H+
k_a=18; %s^-1 Backward constant
k_b=0.04; %s^-1 Forward constant
k_21=1e7; %s^-1 Forward constant
k_12=5e10/1000; %m^3/(mol*s) Backward constant
pak1=@(t,u) [(k_a.*u(2))-(k_b.*u(1)); %dCO2/dt
((k_12.*u(3).*u(4))-k_21.*u(2)); %dH2CO3/dt
-((k_12.*u(3).*u(4))-k_21.*u(2)); %dHCO3/dt
-((k_12.*u(3).*u(4))-k_21.*u(2))]; %dH/dt
opt=odeset('AbsTol',1e-10,'RelTol',1e-12);
[t u] = ode15s(pak1, [tspan], u0, opt);
subplot(2,2,1)
plot(t./3600,u(:,1).*(1e6./1000),'-')
set(gca, 'Nextplot','add','box','on','FontSize',14)
xlabel('Time (hr)','Interpreter','LaTeX');
ylabel('CO$_2$ concentration ($\mu$mol L$^{-1}$)','Interpreter','LaTeX');
subplot(2,2,2)
plot(t./3600,u(:,2).*(1e6./1000),'-')
set(gca, 'Nextplot','add','box','on','FontSize',14)
xlabel('Time (hr)','Interpreter','LaTeX');
ylabel('H$_2$CO$_3$ concentration ($\mu$mol L$^{-1}$)','Interpreter','LaTeX');
subplot(2,2,3)
plot(t./3600,u(:,3).*(1e6./1000),'-')
set(gca, 'Nextplot','add','box','on','FontSize',14)
xlabel('Time (hr)','Interpreter','LaTeX');
ylabel('HCO$_3^{-1}$ concentration ($\mu$mol L$^{-1}$)','Interpreter','LaTeX');
subplot(2,2,4)
plot(t./3600,-log10(u(:,4)./1000),'-')
set(gca, 'Nextplot','add','box','on','FontSize',14)
xlabel('Time (hr)','Interpreter','LaTeX');
ylabel('pH','Interpreter','LaTeX');
2 个评论
Walter Roberson
2021-2-4
What is that 1000 constant multiplier? moles/m^3 hints to me that the constant should be 1000^3 if you converting mm to m. On the other hand the more common base unit for volume is cm in which case 100^3 would seem appropriate.
Jose Bolivar
2021-2-4
编辑：Jose Bolivar
2021-2-4
采纳的回答
Bjorn Gustavsson
2021-2-4
编辑：Bjorn Gustavsson
2021-2-4
You have an error in your continuity-equation, you have forgotten the loss of H2CO3 due to the first backwards reaction H2CO3 -> CO2 + H20: k_a, and you haven't included the source due to the CO2 + H2O: k_b either. I think it should look like this:
pak1=@(t,u) [(k_a.*u(2))-(k_b.*u(1)); %dCO2/dt
((k_12.*u(3).*u(4))+k_b.*u(1)-(k_21 + k_a).*u(2)); %dH2CO3/dt
-((k_12.*u(3).*u(4))-k_21.*u(2)); %dHCO3/dt
-((k_12.*u(3).*u(4))-k_21.*u(2))]; %dH/dt
I haven't checked for more. Your problem is painfully familiar to me, in order to get these right I typically have to print out the code, and explicitly mark the source and loss-terms of each reaction making sure that they are properly included in the reactions for each reactant and product (multicoloured pens help.)
HTH
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Chemistry 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
