How o define two different data set of (t , T) for an ode solver that it's variables are (t , C)?

2 次查看（过去 30 天）
Farangis
Farangis 2023-1-23
评论： Davide Masiello 2023-1-23
How o define two different data set of (t , T) for an ode solver that it's variables are (t , C)?
I have an ode function that gives t (time) and C (Consentration) and there are k values that changes versus T (Temperature) and also temperature changes versus t. I have two different data set of (t , T) and I have to solve the ode equations to obtain C for both (t , T) data set in a one code.
Finally, I need to optimize k and Ea values for the two data set (which will be different from each other) and "alpha" (which should be the same for each data set).
So, I think I do not know how to define these two (t , T) data set for ode solver. I would be thankful if anyone could help me.
Here is my ode:
function dcdt = Cell_deg_60C(t,c,k0_vec,Ea_vec,alpha)
R = 8.314; % J/K.mol
%% For T=60 C
% I need to define T data in a vector like it T = [320.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15]
% and For T=80 C ---> T= [321.15 334.15 348.15 353.15 353.15 353.15 353.15 353.15 353.15 353.15]
%% here I just defined one data set of T, because I have no idea how to define the second one
if t <= 720
T = 297.15 + ((320.15-297.15)*(t/540)); % K , s
else
T = 333.15;
end
%% ------------------------------------------------------------------------------------------------
k01 = k0_vec(1);
k02 = k0_vec(2);
k03 = k0_vec(3);
Ea1 = Ea_vec(1);
Ea2 = Ea_vec(2);
Ea3 = Ea_vec(3);
k1 = (k01.*exp(-Ea1./(R.*T)));
k2 = (k02.*exp(-Ea2./(R.*T)));
k3 = (k03.*exp(-Ea3./(R.*T)));
%% D=1 G=2 M=3 O=4
dcdt = zeros(4,1);
dcdt(1) = - k1.*c(1) - k3.*c(1);
dcdt(2) = k1.*c(1);
dcdt(3) = k1.*c(1) + 2*k3.*c(1) - k2.*c(3);
dcdt(4) = alpha*k2*c(3);
dcdt = [dcdt(1);dcdt(2);dcdt(3);dcdt(4)];
% here is odesolver mfile
c0 = zeros(4,1);
c0(1) = 0.14; % mol/l
c0(2) = 0.00045;
c0(3) = 0.01529;
c0(4) = 0.00101;
k0_vec = [0.63, 0.03, 0.951];
Ea_vec = [17000, 10000, 100000];
alpha = 1.44;
tspan = 60*[
9
12
27
42
57
72
87
102
117
132
]; % second
[t,c] = ode45(@(t,c) Cell_deg_60C(t,c,k0_vec,Ea_vec,alpha),tspan,c0);
plot (t/60,c);
legend('Disac','GISA','Monosac','Otheracids')
xlabel('time (min)');
ylabel('C at T=60C (mol/L)');

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

采纳的回答

Davide Masiello
Davide Masiello 2023-1-23
Ran in:
I think your code is fine, below you can see a slightly different approach that uses interpolation to define the T(t) dependence.
Regarding optimization, you'd have to provide experimental concentration profiles to carry it out.
% Parameters and ICs
c0 = [0.14 0.00045 0.01529 0.00101]; % Initial concentrations, mol/L
k0 = [0.63, 0.03, 0.951]; % Pre-exponential factors
Ea = [17000, 10000, 100000]; % Activation energies
alpha = 1.44; % Time vector, s
% Experimental temperatures and interpolation
time = 60*[9 12 27 42 57 72 87 102 117 132];
TEMP1 = [320.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15];
TEMP2 = [321.15 334.15 348.15 353.15 353.15 353.15 353.15 353.15 353.15 353.15];
T1 = @(t)interp1(time,TEMP1,t);
T2 = @(t)interp1(time,TEMP2,t);
[t,c] = ode45(@(t,c) Cell_deg_60C(t,c,k0,Ea,alpha,T1),[time(1),time(end)],c0);
figure(1)
plot (t/60,c);
title('Kinetics with first temperature evolution')
xlabel('time (min)');
ylabel('C at T=60C (mol/L)');
legend('Disac','GISA','Monosac','Otheracids','Location','Best')
[t,c] = ode45(@(t,c) Cell_deg_60C(t,c,k0,Ea,alpha,T2),[time(1),time(end)],c0);
figure(2)
plot (t/60,c);
title('Kinetics with second temperature evolution')
xlabel('time (min)');
ylabel('C at T=80C (mol/L)');
legend('Disac','GISA','Monosac','Otheracids','Location','Best')
function dcdt = Cell_deg_60C(t,c,k0,Ea,alpha,T1)
R = 8.314; % Universal gas constantJ/K.mol
k = (k0.*exp(-Ea/(R*T1(t)))); % Rate constants
dcdt(1,1) = -k(1)*c(1)-k(3)*c(1);
dcdt(2,1) = k(1)*c(1);
dcdt(3,1) = k(1)*c(1)+2*k(3)*c(1)-k(2)*c(3);
dcdt(4,1) = alpha*k(2)*c(3);
end
  4 个评论
显示 2更早的评论

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Assembly 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by