Hello! I'm trying to model turbulent mixing with chemical reactions using engulfment model and ode45. To do this I need to discretise into sigma parts the added volume of acid to the solution. I must solve a system of ordinary differential equations for each part of the feed (each time until one of reagents is consumed). From one part to another, I need to update concentrations of reagents in the vicinity of the reaction zone (for the (k+1)th part of the feed concentrations are Ci(k), but I changed matrix indices, because I can't start from 0). I'd like to save all the final values of variables for each part. These values are essential to calculate Ci(k+1) for the next part of the feed.
The script usually stops after calculation for the 1st part. The error which constantly occurs is "index exceeds array bounds" (additionally with errors concerning used functions). I've tried to change indexing, but without success. I currently have no idea how to properly fix this, so I'd greatly appreciate any remarks and suggestions. Thank you!
This is my event function:
function [value,isterminal,direction] = eventsfunction(~,c,~)
value = c(3);
isterminal = 1;
direction = -1;
end
This is main function:
function [dcdt] = mojafunkcja(t,c,k)
g = 9.80665;
Dm = 0.052;
T = 273.15+20;
ro = 998.2063;
M1 = -52.843;
M2 = 3703.6;
M3 = 5.866;
M4 = -5.879*10^(-29);
M5 = 10;
mi = exp(M1+M2/T+M3*log(T)+M4*T^M5);
ni = mi/ro;
Nm = 500/60;
Rem = (Nm*Dm^2)/ni;
Frm = (Dm*Nm^2)/g;
Eum = 1.1*Frm^((1-log10(Rem))/40);
Vo = 0.2;
epsilon = (Eum*Nm^3*Dm^5)/(Vo*10^(-3));
Vbo = 0.001;
E = 0.05776*(epsilon/ni)^(1/2);
Io = 0.10740;
k2 = 10^(2.1766*Io+8.9484);
k3 = 5.6*10^9;
k33 = 7.5*10^6;
tf = 4.5;
sigma = 100;
Q = (Vbo*0.001/sigma)/tf;
u = pi*Dm*Nm;
Lambdac = (Q/(pi*u))^(1/2);
ts = 2*(Lambdac^2/epsilon)^(1/3);
C1(1) = 0.01179;
C2(1) = 0.00234;
C3(1) = 0.0933;
C4(1) = 0.0898;
C5(1) = 0;
C6(1) = 0;
dcdt = zeros(8,1);
dcdt(1) = E*(C1(k)-c(1))*(1-(c(7)*exp(-t/ts))/(Vbo/sigma))-5*k2*(c(1)^2)*c(2)*(c(3)^2)-k3*c(5)*c(1)+k33*c(6);
dcdt(2) = E*(C2(k)-c(2))*(1-(c(7)*exp(-t/ts))/(Vbo/sigma))-k2*(c(1)^2)*c(2)*(c(3)^2);
dcdt(3) = -E*C3(k)*(1-(c(7)*exp(-t/ts))/(Vbo/sigma))-6*k2*(c(1)^2)*c(2)*(c(3)^2);
dcdt(4) = E*C3(k)*(1-(c(7)*exp(-t/ts))/(Vbo/sigma))+E*(C4(k)-c(4))*(1-(c(7)*exp(-t/ts)));
dcdt(5) = E*(C5(k)-c(5))*(1-(c(7)*exp(-t/ts))/(Vbo/sigma))+3*k2*(c(1)^2)*c(2)*(c(3)^2)-k3*c(5)*c(1)+k33*c(6);
dcdt(6) = E*(C6(k)-c(6))*(1-(c(7)*exp(-t/ts))/(Vbo/sigma))+k3*c(5)*c(1)-k33*c(6);
dcdt(7) = E*c(7)*(1-(c(7)*exp(-t/ts))/(Vbo/sigma));
dcdt(8) = E*c(8)*(1-(c(7)*exp(-t/ts))/(Vbo/sigma));
end
This is my script:
clc;
Vo = 0.2;
Vbo = 0.001;
sigma = 100;
h = 10^(-6);
t0 = 0;
deltat = 1;
tk = t0+deltat;
c10 = 0;
c20 = 0;
c30 = 0.41273;
c40 = 0;
c50 = 0;
c60 = 0;
c70 = Vbo/sigma;
c80 = (Vbo/sigma)/(Vo+Vbo/sigma);
c0 = [c10 c20 c30 c40 c50 c60 c70 c80];
tspan = t0:h:tk;
for k = 1:(sigma)
C1(1) = 0.01179;
C2(1) = 0.00234;
C3(1) = 0.0933;
C4(1) = 0.0898;
C5(1) = 0;
C6(1) = 0;
options = odeset('Events',@eventsfunction);
[t,c,tk,ck] = ode45(@(t,c) mojafunkcja(t,c,k),tspan,c0,options);
tkk(k) = tk;
c1k(k) = ck(1,1);
c2k(k) = ck(1,2);
c3k(k) = ck(1,3);
c4k(k) = ck(1,4);
c5k(k) = ck(1,5);
c6k(k) = ck(1,6);
c7k(k) = ck(1,7);
c8k(k) = ck(1,8);
C1(k+1) = C1(k)*(1-c8k(k))+c1k(k)*c8k(k);
C2(k+1) = C2(k)*(1-c8k(k))+c2k(k)*c8k(k);
C3(k+1) = C3(k)*(1-c8k(k));
C4(k+1) = C4(k)*(1-c8k(k))+c4k(k)*c8k(k);
C5(k+1) = C5(k)*(1-c8k(k))+c5k(k)*c8k(k);
C6(k+1) = C6(k)*(1-c8k(k))+c6k(k)*c8k(k);
Vc(k) = Vo+(Vbo*k)/(sigma);
Xs(k) = ((C6(k+1)+C5(k+1))*Vc(k))/(c30*(Vbo*k/(sigma)))*(2+C3(1)/(3*C2(1)));
end
