DAE ode15s coupled initial conditions / boundary condition

1 次查看（过去 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 中查找有关 Programming 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by