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');

[Walter Roberson]

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]

Hi Walter

I multiply with 1000 so all my units are in mol/m^3. I think the code should work maybe I have defined wrong the rates of reactions, what do you think?

Answer 1

[Bjorn Gustavsson]

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

[Jose Bolivar]

Yes! I had the feeling something was off with the equilibrium reactions, thank you so much for helping me again! :D

[Bjorn Gustavsson]

My pleasure!