DAE ode15s coupled initial conditions / boundary condition

2 次查看（过去 30 天）
Marko
Marko 2019-11-9
评论： darova 2019-11-9
Hi guys,
I try to solve differential-algebraic equations - which is basicly describing the dissociation of some Acid in water with another Base.
The Systems looks like this:
%Initial conditions
initital_c_C0 = 0;
initital_c_C1 = 0;
initital_c_C2 = 0;
initital_c_H = 0;
initital_c_Na = 0;
initital_c_OH = 0;
initital_c_NaOH = 0;
initital_F_IA = 0.5;
initital_F_NaOH = 1;
%% Definition of differential-algebraic equation system using the so called mass matrix M(t,y)*y' = f(t,y)
N = zeros(9);
for i = 1 : (7)
N(i,i) = 1;
end
M=sparse(N);
%% Call of the ode-solver
% feeding
startvector = [initital_c_C0; initital_c_C1; initital_c_C2; initital_c_H; initital_c_Na; initital_c_OH; initital_c_NaOH; initital_F_NaOH; initital_F_IA];
disp('Start ODE15S');
t_end = 7;
time_interval = [0:1:t_end];
options = odeset('Mass', M, 'RelTol', 1e-5, 'Abstol', 1e-5, 'MStateDependence', 'weak', 'MassSingular','yes');
And my DAE
%% Startvector
c_C0 = startvector(1);
c_C1 = startvector(2);
c_C2 = startvector(3);
c_H = startvector(4);
c_OH = startvector(5);
c_Na = startvector(6);
c_NaOH = startvector(7);
F_Na = startvector(8);
F_IA = startvector(9);
%% DAE
% reaction rates
f(1,:) = f_1^2 * k_1_R * c_C1 * c_H - f_1 * k_1_H * c_C0;
f(2,:) = f_2 * f_1 * k_2_R * c_C2 * c_H - f_1 * k_2_H * c_C1 - f_1^2 * k_1_R * c_C1 * c_H + f_1 * k_1_H * c_C0;
f(3,:) = - f_2 * f_1 * k_2_R * c_C2 * c_H - f_1 * k_2_H * c_C1;
f(4,:) = - f_1^2 * k_1_R * c_C1 * c_H + f_1 * k_1_H * c_C0 - f_2 * f_1 * k_2_R * c_C2 * c_H + f_1 * k_2_H * c_C1 + k_w_H - f_1^2 * k_w_R * c_H * c_OH;
f(5,:) = k_w_H - f_1^2 * k_w_R * c_H * c_OH;
f(6,:) = k_Na_H * c_NaOH - k_Na_R * f_1^2 * c_Na * c_OH;
f(7,:) = - k_Na_H * c_NaOH + k_Na_R * f_1^2 * c_Na * c_OH;
% Feed mass balance
f(8,:) = c_Na + c_NaOH - F_Na;
f(9,:) = c_C0 + c_C1 + c_C2 - F_IA ;
So what I basicly want to implement is that:
% Feed mass balance
f(8,:) = c_Na + c_NaOH - F_Na;
f(9,:) = c_C0 + c_C1 + c_C2 - F_IA ;
is true for for each time or in other words - that the sum of the components c_C0, c_C1 and c_C2 is always F_IA (constant). Same for c_NA and c_NaOH.
How could I implement this logic in matlab?
For the sake of clarity I did not show the code where I wrote down the constants, K, k, f and so on.
Thanks!
  3 个评论
显示 1更早的评论
Marko
Marko 2019-11-9
Yes you are right.
So what I did was to use the constraints in a algebraic equation and eliminate one differential equation. This way I still have a determined system.
Anyway, the result is still wrong.
Heres my code:
%Initial conditions
initital_c_C0 = 0;
initital_c_C1 = 0;
initital_c_C2 = 0;
initital_c_H = 0;
initital_c_Na = 0;
initital_c_OH = 0;
initital_c_NaOH = 0;
%% Definition of differential-algebraic equation system using the so called mass matrix M(t,y)*y' = f(t,y)
N = zeros(9);
for i = 1 : (7)
N(i,i) = 1;
end
M=sparse(N);
%% Call of the ode-solver
% feeding
startvector = [initital_c_C0; initital_c_C1; initital_c_C2; initital_c_H; initital_c_Na; initital_c_OH; initital_c_NaOH; initital_F_NaOH; initital_F_IA];
disp('Start ODE15S');
t_end = 7;
time_interval = [0:1:t_end];
options = odeset('Mass', M, 'RelTol', 1e-5, 'Abstol', 1e-5, 'MStateDependence', 'weak', 'MassSingular','yes');
function f = equations_solver_mass_balance_pH(t, startvector)
%% Paramater
% activity coefficient - Index is z (charge number)
f_1 = 1;
f_2 = 1;
% forward reaction rate constants set to arbitrary number
k_1_H = 1;
k_2_H = 1;
k_w_H = 1;
k_Na_H = 1;
% Dissociation constants
K1 = 10^(-3.85);
K2 = 10^(-5.5);
Kw = 10^(-14);
K_Na = 10^(13);
% calculated backward reaction rate
k_1_R = k_1_H/K1;
k_2_R = k_2_H/K2;
k_w_R = k_w_H/Kw;
k_Na_R = k_Na_H/K_Na;
%Feed concentration
F_IA = 0.5;
F_Na = 1;
%% Startvector
c_C0 = startvector(1);
c_C2 = startvector(2);
c_H = startvector(3);
c_OH = startvector(4);
c_Na = startvector(5);
c_C1 = startvector(6);
c_NaOH = startvector(7);
%% DAE
% reaction rates
f(1,:) = f_1^2 * k_1_R * c_C1 * c_H - f_1 * k_1_H * c_C0;
%f(2,:) = f_2 * f_1 * k_2_R * c_C2 * c_H - f_1 * k_2_H * c_C1 - f_1^2 * k_1_R * c_C1 * c_H + f_1 * k_1_H * c_C0;
f(2,:) = - f_2 * f_1 * k_2_R * c_C2 * c_H - f_1 * k_2_H * c_C1;
f(3,:) = - f_1^2 * k_1_R * c_C1 * c_H + f_1 * k_1_H * c_C0 - f_2 * f_1 * k_2_R * c_C2 * c_H + f_1 * k_2_H * c_C1 + k_w_H - f_1^2 * k_w_R * c_H * c_OH;
f(4,:) = k_w_H - f_1^2 * k_w_R * c_H * c_OH;
f(5,:) = k_Na_H * c_NaOH - k_Na_R * f_1^2 * c_Na * c_OH;
% Feed mass balance
f(6,:) = c_C0 + c_C1 + c_C2 - F_IA ;
f(7,:) = c_Na + c_NaOH - F_Na;
end
darova
darova 2019-11-9
I think you are overcompicating the problem
%% Startvector
% ...
c_C1 = F_FIA - c_C2 - c_C0;
c_NaOH = F_Na - c_Na;
%% DAE

请先登录，再进行评论。

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

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Solver Outputs and Iterative Display 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by