Bonjour !

J'ai un problème dans la résolution d'un simple système avec une ODE (dX/dt=-J).
Mon problème c'est que J est calculé dans une boucle pour chaque temps Durant un certain intervalle. Donc, je dois donner J à lsode à chaque temps, et je ne sais pas comment le faire...
Aussi, J peut être calculé de 2 manières différents, dépendant de la valeur de X (il y a une valeur limite Xc pour laquelle, si X est supérieur, il sera calculé d'une certaine manière et si X est inférieur et sera calculé avec l'autre equation).
Donc, je dois comparer la valeur de X pour chaque temps pour calculer J...
En plus, dans J, il y a plein de paramètres qui doit aussi être calculés à partir d'autres paramètres...

Est-ce que quelque a une idée de comment je peux résoudre ça ?
Merci merci !

PS: dans le code suivant, j'ai mis pas mal de choses en commentaire pour rendre le run plus rapide

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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
%Définition des paramètres globaux
Xc=1;             %Humidité des levures pour laquelle la fin de la phase de 
                  %séchage à vitesse constante est observée
p=101325;         %Pression totale (= atmosphérique) [Pa]
ke=0.12;          %Coefficient de transfert externe (déterminé exp) [m/s]
Mw=0.018;         %Masse molaire de l'eau [kg/mol]
Ma=0.02896;       %Masse molaire de l'air [kg/mol]
a=6;              %Aire de la surface extérieure des levures par unité de masse  
                  %de matière sèche [m²/kg]
Rg=8.31;          %Constante universelle des gaz parfaits [J/mol/K]
eps=0.3;          %Porosité du lit, doit être compris entre [0.3;0.6]
tau=5;            %Tortuosité, doit être comprise entre [2;10]
D=0.000025;       %Coefficient de diffusion de la vapeur d'eau dans le gaz [m²/s]
Lk=2250000;       %Chaleur latente de vaporisation de l'eau [J/kg]
cpa=1000;         %Chaleur spécifique massique de l'air sec [J/K/kg]
Tsat0=373.15;     %Température de référence pour saturation [K]
psat0=101315;     %Pression de saturation de l'eau à la température T0 [Pa]
Lm=2470000*0.018; %Chaleur latente de vaporisation (molaire)
 
%Paramètres variables
%X0=input("Entrez la teneur initiale en eau des levures:\n");  %(2.45 valeur ref)
X0=2.45;
%Xr=input("Entrez la teneur residuelle en eau des levures:\n");%(0.1 valeur ref)
Xr=0.1;
Xc=1;             %Valeur ref d'humidité critique
%R=input("Entrez le rayon moyen des grains de levure (mm):\n");
%R=R/1000;         %Conversion en m
R=0.0005;
%M=input("Entrez la masse de levure non-sechee (kg):\n");
M=1500;
Ms=M/(1+X0);      %Calcul de la masse de solide sec après séchage complet
%G=input("Entrez le debit d'air 'sec' a l'entree du lit fluidise (L/h):\n");
%G=(G/1000)*1.2;   %Conversion en kg/h (38500 valeur ref)
G=38500;
%ts=input("Entrez le temps de sechage total (s):\n");
ts=50;
 
%Conditions initiales
Y0=0.009;         %Humidité initiale de l'air (valeur ref)
T0=273.15+90;     %A CHANGER - Température initiale de l'air
 
Y=Y0;
T=T0;
 
Jres=[];         %Création de la matrice résultat pour taux d'évaporation
Yres=[Y];        %Création de la matrice résultat pour l'humidité de l'air
 
step=0.01;
temps=[0:step:ts+step]; %Création de la matrice temps - Vecteur ligne
temps=temps.';          %Transposée de la matrice temps --> Vecteur colonne
i=1;
 
for t=0:step:ts
    [X, ISTATE, MSG]=lsode("humsolide",X0,temps);
    psat=psat0*exp(-(Lm/Rg)*((1/T)-(1/T0)));
    if X(i)>=Xc
      %Calcul du taux d'évaporation
      J=a*Mw*ke*((psat/(Rg*T))-(p/(Rg*T))*(Y/((Mw/Ma)+Y)));
    elseif X(i)<=Xc
      %Calcul des paramètres
      Ri=R*(((X-Xr)/(Xc-Xr))^(1/3));
      ki=((eps*D)/(tau*((R)^2)))*(1/((1/Ri)-(1/R)));
      %Calcul du taux d'évaporation
      J=a*Mw*(1/((1/ke)+(1/ki)))*((psat/(Rg*T))-(p/(Rg*T))*(Y/((Mw/Ma)+Y))); 
  end
  Y=((Ms*J)/G)+Y0;
  T=T0-((Ms*J*Lk)/(cpa*G));
  %Ajout des résultats dans les matrices résultats
  Yres=[Yres;Y];
  Jres=[Jres;J];
  i+=1;
end
 
 
 
%plot(temps,Jres);
%hold on;
%plot(temps,Yres);
%hold on;
%plot(temps,X);
Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
 
function XDOT=humsolide(X0,t)
  XDOT=-J;
  endfunction
Céline