Discussion :

[nini4487]

Bonjour à tous,
Je vais essayer de présenter mon problème. Je suis en train de faire une simulation de pollution sur matlab. On me donne mes différentes équations chimiques avec leur constante je mets tout ça dans matlab, jusque là pas trop de problème. Certaine constante sont dépendant de la luminosité solaire selon la variable "ki". Je dois ensuite résoudre mon équation différentielle en tenant compte de la variation de la constante cinétique et là ça coince. Car le ODE que j'utilise ne reconnait pas ce que je fais. Si quelqu'un pourrait m'aider ou m'indiquer comment faire... ça serait super sympa!

D'avance merci

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
clc;
clear;
 
global m
 
% Input parameters
% ----------------
 
dt = 0.01;
 
ki = 1;
 
t = (0:0.1:1440);
ki=[ ];
p=pi();
 
for i=1:length(t)
    if sin((2*p*((t(i)/60)-6))/24) > 0
        ku=sin((2*p*((t(i)/60)-6))/24);
        ki=[ku ki];
    else
        ku=0;
        ki=[ku ki];
    end
end
 
 
m.k1 = 0.533*ki;
m.k2 = 2.21*10^(-5);
m.k3 = 26.7;
m.k4a = 1.6*10^(-3*ki);
m.k4b = 2.11*10^(-3)*ki;
m.k5 = 1.62*10^4;
m.k6 = 1.22*10^4;
m.k7 = 1.62*10^4;
 
% Initial conditions
% ------------------
 
Cstart = zeros(13,1);
 
Cstart(1) = 0.01;                    % NO   --> dC(1); C(1)
Cstart(2) = 0.1;                     % NO2  --> dC(2); C(2)
Cstart(3) = 0;                       % O    --> dC(3); C(3)
Cstart(4) = 0;                       % O3   --> dC(4); C(4)
Cstart(5) = 0;                       % HO2  --> dC(5); C(5)
Cstart(6) = 0;                       % H2   --> dC(6); C(6)
Cstart(7) = 0;                       % HNO3 --> dC(7); C(7)
Cstart(8) = 0;                       % H2O  --> dC(8); C(8)
Cstart(9) = 0.1;                     % CHOH --> dC(9); C(9)
Cstart(10) = 0;                      % CO   --> dC(10); C(10)
Cstart(11) = 0;                      % OH   --> dC(11); C(11)
Cstart(12) = 0.21*10^6;              % O2   -->dC(12); C(12)
Cstart(13) = 10^6;                   % M    -->dc(13); C813
 
t = (0:dt:1440);
 
[t,C] = ode15s (@Pollution_model, t, Cstart);
 
figure(1)
plot(t/60,C(:,1:11),'-')
xlabel ('Temps [h]')
ylabel ('concentration [ppm]')
title ('Variation d espèce présente dans l air en fonction du temps')
legend ('NO', 'NO2', 'O', 'O3', 'HO2', 'H2', 'HNO3', 'H2O', 'CHOH', 'CO', 'OH', 'O2', 'M')

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
function dC = Pollution_model(t,C)
 
%1.  RNO = + r1 - r3 - r6 = + k1[NO][CO] - k3[O3][NO] - k6[NO][HO2]
%2.  RNO2 = - r1 + r3 + r6 - r7 = - k1[NO2] + k3[NO2][O2] + k6[NO2][OH] - k7[NO2][OH]
%3.  RO = + r1 - r2 = + k1[NO][O] - k2[O][O2][M]
%4.  RO3 = + r2 - r3 = + k2[O3][NO]
%5.  RHO2 = + r4a - r5 - r6 = + 2 k4a[HO2]^2 - k5[HO2][CO][H2O] - k6[HO2][NO]
%6.  RH2 = + r4b = + k4b[H2][CO]
%7.  RHNO3 = + r7 = + k7[HNO3]
%8.  RH2O = + r5 = + k5[HO2][CO][H2O]
%9.  RHCHO = - r4a - r4b - r5 = - k4a[HO2]^2[CO] - k4b[HCHO] - k5[OH][NO2]
%10. RCO = - r4a - r4b + r5 = - k4a[HCOH][OH] - k4b[H2][CO] + k5[HO2][CO][H2O]
%11. ROH = - r5 + r6 - r7 = - k5[HCOH][OH] + k6[NO2][OH] - k7[OH][NO2]
 
% NO   --> dC(1); C(1)
% NO2  --> dC(2); C(2)
% O    --> dC(3); C(3)
% O3   --> dC(4); C(4)
% HO2  --> dC(5); C(5)
% H2   --> dC(6); C(6)
% HNO3 --> dC(7); C(7)
% H2O  --> dC(8); C(8)
% CHOH --> dC(9); C(9)
% CO   --> dC(10); C(10)
% OH   --> dC(11); C(11)
% O2   --> dC(12); C(12)
% M    --> dC(13); C(13)
 
dC = zeros(13,1);
 
global m
%m.k    Constantes cinétiques
 
%Calcul des vitesses de réactions
 
 
 
r1 = m.k1* C(2);
r2 = m.k2 * C(3);
r3 = m.k3*C(4) * C(1);
r4a = m.k4a* C(9);
r4b = m.k4b* C(9);
r5 = m.k5 * C(9) * C(11);
r6 = m.k6 * C(5) * C(1);
r7 = m.k7 * C(11) * C(2);
 
%équations différentielles
 
dC(1)  = + r1 - r3 - r6;
dC(2)  = - r1 + r3 + r6 - r7; 
dC(3)  = + r1 - r2;
dC(4)  = + r2 - r3;
dC(5)  = + 2*r4a + r5 - r6;
dC(6)  = + r4b;
dC(7)  = + r7;
dC(8)  = + r5;
dC(9)  = - r4a - r4b - r5;
dC(10) = + r4b + r4a + r5;
dC(11) = - r5 + r6 - r7;
dC(12) = 0;
dC(13) = 0;
 
end

[Invité]

Bonjour,

Pourrais-tu préciser ton problème ?

Remarque :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
for i=1:length(t)
    if sin((2*p*((t(i)/60)-6))/24) > 0
        ku=sin((2*p*((t(i)/60)-6))/24);
        ki=[ku ki];
    else
        ku=0;
        ki=[ku ki];
    end
end

Attention à ce genre de boucle, dont le tableau grandi à chaque itération, qui peut prendre beaucoup de temps. Voir Qu'est-ce que la préallocation de mémoire ?
Elle construit de plus un vecteur ki inversé par rapport au vecteur t, est-ce normal ?