Salve,
scrivo per la prima volta perchè mi trovo davanti ad un problema per la mia tesi. Devo simulare un modello matematico (da me scritto) non molto complesso composto da una certa quantità (circa 10) variabili legate tra di loro da equazioni differenziali ordinarie e equazioni algebriche semplici. Premetto che il mio livello su Matlab non è molto avanzato ed il mio problema è il seguente:
ho scritto il codice dove vengono scritte sottoforma di matrice tutte le equazioni differenziali. Quando risolvo con il comando
ode15s alcuni parametri mi rimangono costanti per tutta la durata temporale e non c'è nessuna sorta di consumo di essi. Ora quello che non capisco è:
devo creare una sorta di loop dove ogni volta i valori calcolati vengono cambiati? Chiedo questo perchè ho dei termini cinetici che usano come variabili le soluzioni di alcune equazioni differenziali.
Forse non sono stato molto chiaro per qualunque cosa chiedete pure e vi spiego meglio.
Grazie mille. Vi supplico aiutatemi
Vi allego il codice matlab che ho scritto:
function dy = modello(t,var)
%Variabili
X1=var(1) ;
X2=var(2) ;
X3=var(3) ;
S1=var(4) ;
S2=var(5) ;
S3=var(6) ;
S4=var(7) ;
S5=var(8) ;
S3G=var(9) ;
S4G=var(10) ;
S5G=var(11) ;
%Parametri
mu_max1 = 0.212 ; %h^-1
mu_max2 = 0.0255 ; %h^-1
mu_max3 = 0.333 ; %h^-1
ks1 = 0.5 ; % g/L
ks2 = 0.537 ; % g/L
ki2 = 15.36 ; % g/L
ks3 = 3.13*10^(-7) ; % g/L
Y1 = 12.86 ; %g/gbiomassa
Y2 = 3.54 ; %g/gbiomassa
Y3 = 24.14 ; %g/gbiomassa
Y4 = 2.41 ; %g/gbiomassa
Y5 = 16.73 ; %g/gbiomassa
Y6 = 17.83 ; %g/gbiomassa
Y7 = 0.125 ; %g/gbiomassa
Y8 = 3.09 ; %g/gbiomassa
Y9 = 6.08 ; %g/gbiomassa
Y10 = 5.79 ; %g/gbiomassa
H_co2 = 7.48 ; %g/atm*L
H_h2 = 0.0016 ; %g/atm*L
H_ch4 = 0.016 ; %g/atm*L
kla1 = 0.294 ; % h^-1
kla2 = 6.6 ; % h^-1
kla3 = 0.261 ; % h^-1
Vg = 9.3 ; % L
Vl = 52.7 ; % L
T = 328 ; % K
P = 1 ; % atm
pw = 0.155 ; % atm
R = 0.0821 ; % mol*L/atm*K
q_in = 3 ; %l/h
S1_in = 41.619 ; % g/L
S2_in = 0 ; % g/L
S3G_in = 0; % g/L
S4G_in = 0.0744 ; % g/L
S5G_in = 0; % g/L
q_l = 53.34 ; % l/h
PM_co2 = 44 ; % g/mol
PM_h2 = 2 ; % g/mol
PM_ch4 = 16 ; % g/mol
D = q_l/Vl ; % h^-1
%Equazioni
% % syms X1 X2 X3 S1 S2 S3 S4 S5 S3G S4G S5G mu1 mu2 mu3 S3star S4star S5star q_out
mu1=mu_max1*S1/(S1+ks1);
mu2=mu_max2*S2/(S2+ks2+(S2^2)/ki2);
mu3=mu_max3*S4/(S4+ks3);
S3star=(S3G/(S3G+S4G+S5G))*H_co2 ;
S4star=(S4G/(S3G+S4G+S5G))*H_h2 ;
S5star=(S5G/(S3G+S4G+S5G))*H_ch4 ;
q_out=R*T/(P-pw) *VL*(((kla1/PM_co2)*(S3-S3star))-((kla2/PM_h2)*(S4-S4star))+(kla3/PM_ch4)*(S5-S5star));
dy(1,1)=(mu1-q_l/Vl-0.05*mu1)*X1 ;
dy(2,1)=(mu2-q_l/Vl-0.027*mu2)*X2 ;
dy(3,1)=(mu2-q_l/Vl-0.05*mu2)*X3 ;
dy(4,1)=(q_l/Vl)*(S1_in - var(4))-Y1*mu1*X1 ;
dy(5,1)=-(q_l/Vl)*S2 + Y2*mu1*X1-Y3*mu2*X2 ;
dy(6,1)=Y4*mu1*X1+Y5*mu2*X2-Y6*mu3*X3-kla1*(S3-(P*H_co2*S3star));
dy(7,1)=Y7*mu1*X1-Y8*mu3*X3+kla2*(P*H_h2*(S4star)-S4);
dy(8,1)=Y9*mu2*X2+Y10*mu3*X3-kla3*(S5-(P*H_ch4*S5star));
dy(9,1)=q_in*S3G_in/Vg -q_out*S3G/Vg +(Vl/Vg)*kla1*(S3-S3star);
dy(10,1)=q_in*S4G_in/Vg -q_out*S4G/Vg +(Vl/Vg)*kla2*(S4-S4star);
dy(11,1)=q_in*S5G_in/Vg -q_out*S5G/Vg +(Vl/Vg)*kla3*(S5-S5star);
end
clc ; clear all; close all
range=[0:1:720];
IC=[0.25;0.15;0.15;41.619;0;0;0;0;0;0.0744;0];
[time, sol]=ode15s(@Prova1,range,IC)