How o define two different data set of (t , T) for an ode solver that it's variables are (t , C)?
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Farangis
le 23 Jan 2023
Commenté : Davide Masiello
le 23 Jan 2023
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)');
2 commentaires
Davide Masiello
le 23 Jan 2023
Do you have a vector of temperatures versus time?
Or do you have a T(t) analytical function?
Réponse acceptée
Davide Masiello
le 23 Jan 2023
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 commentaires
Davide Masiello
le 23 Jan 2023
No problem, please consider to click on the Accept button if you think it was useful. Cheers.
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Assembly dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!