How o define two different data set of (t , T) for an ode solver that it's variables are (t , C)?
1 view (last 30 days)
Show older comments
Farangis
on 23 Jan 2023
Commented: Davide Masiello
on 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 Comments
Davide Masiello
on 23 Jan 2023
Do you have a vector of temperatures versus time?
Or do you have a T(t) analytical function?
Accepted Answer
Davide Masiello
on 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 Comments
Davide Masiello
on 23 Jan 2023
No problem, please consider to click on the Accept button if you think it was useful. Cheers.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!