[Matlab] Résoudre un système linéaire itératif

a marqué ce sujet comme résolu.
Auteur du sujet

Bonjour à tous,

La semaine dernière je vous demandais comment résoudre un système non-linéaire relativement simple, cette semaine c’est un système linéaire itératif (car il est implicite). Pour faire simple, je dois résoudre un système assez gros (j’ai pas compté mais ça doit faire une quarantaines d’équations) mais de façon itérative… C’est-à-dire que je fais une estimation d’un paramètre (ici, la température) et je recalcule tous les paramètres jusqu’à avoir convergence (mon critère est: abs(Told - Tnew) < 0.1).

J’ai pas encore tenté d’écrire le code mais j’ai déjà toutes mes équations (je pense) et donc je me demandais quelle était l’approche la plus simple et plus efficace pour le faire ? Comme vous le savez, j’y connais très peu en programmation / Matlab donc ça m’aiderait d’avoir une explication de comment le coder!

Je mettrai mon code après parce qu’il risque de toute façon d’avoir des pépins :p

Merci d’avance!

Edit:

D’un point de vue de la structure du code ça devrait je pense ressembler à ça (conceptuellement)

BLOC AVEC LE SYSTEME D’EQ LINEAIRE (44) Ax = B

À LA FIN DU SYSTÈME, CALCULER DES VARIABLES AVEC LES SOLUTIONS DE CE SYSTEME

CALCULER DE NOUVELLES VARIABLES (POUR OBTENIR MA TEMPERATURE)

CONDITION "IF" T_ANCIEN - T_NOUVEAU < 0.1 => OK, LES SOLUTIONS SONT LES BONNES

SINON ("ELSE") -> RECALCULER

(J’ai mis en majuscule pour montrer que je vois ça sous forme de différents blocs / étapes).

Compliqué à mettre en oeuvre? :p

Édité par sotibio

+0 -0

Salut,

je dois résoudre un système assez gros (j’ai pas compté mais ça doit faire une quarantaines d’équations) mais de façon itérative

C’est ridiculement petit, résoudre ça de façon directe serait beaucoup plus malin. Il t’est imposé d’utiliser une méthode itérative ? Est-ce que ta matrice est particulière (symétrique et/ou creuse par exemple) ?

Sinon, les méthodes itératives sont déjà implémentées, voir lsqr, gmres, minres pour ne citer que les plus générales (hésite pas à piocher dans "See Also" à la fin des docs de ces fonctions). Remarque qu’elles signalent toutes qu’elles sont pensées pour des matrices grandes (40, c’est pas grand) et creuses (avec plein de zéros).

I don’t mind that you think slowly, but I do mind that you are publishing faster. — W. Pauli

+1 -0
Auteur du sujet

Ça m’est imposé car c’est un problème d’ingénierie et je dois estimer une température dans une colonne de distillation sauf que je sais pas si la bonne mais comme je connais le flux de produit que je veux sortir, je dois faire en sorte d’avoir les bonnes températures.

J’ai pas complètement écrit la matrice la matrice mais elle est assez creuse (pas symétrique par contre). Je pourrais la copier ici mais je sais pas si c’est une bonne idée au vu de la taille de celle-ci :p

Je vais regarder les fonctions que tu proposes!

+0 -0

Ça m’est imposé car c’est un problème d’ingénierie et je dois estimer une température dans une colonne de distillation sauf que je sais pas si la bonne mais comme je connais le flux de produit que je veux sortir, je dois faire en sorte d’avoir les bonnes températures.

Je suis en train de me dire qu’on ne parle peut être pas la même langue, alors je vais essayer de mieux comprendre ton problème.

Est-ce que tu dois résoudre un système du genre MT=RMT=R avec MM une matrice constante, TT ton vecteur d’inconnues et RR un vecteur quelconque constant ? Si oui, je ne vois pas en quoi le problème concret que tu tentes de résoudre va imposer le choix d’une méthode itérative.

Ou est-ce que MM et RR dépendent de TT auquel cas ton problème n’est plus linéaire et le choix d’une méthode itérative serait plus compréhensible.

Plutôt que copier la matrice, t’aurais pas une forme générale pour les équations vérifiées par la température ?

I don’t mind that you think slowly, but I do mind that you are publishing faster. — W. Pauli

+1 -0
Auteur du sujet

Je vais expliquer le problème d’abord pour être plus clair sinon ça reste très flou pour vous. Je dois résoudre un problème de chimie industrielle pour la production d’un polymère (styrène). J’ai deux réactions qui ont lieu dans le réacteur et mes trois derniers éléments (Flash, Dist 1 et Dist 2) sont des unités de séparations où je ne connais pas la température. Cependant je connais la production de styrène que je veux obtenir à la fin. Vu que j’ai cette infos (et d’autres comme la pression etc.), je peux écrire pleins d’équations de conservation de masse etc.

Synthèse industrielle du styrène
Synthèse industrielle du styrène

Voici quelques équations (je les mets pas toutes mais elles se ressemble concernant la forme):

μ0EB=5μ0B\mu _0^{EB} = 5\mu _0^B

iμ0iM(i)=iμ3.1iM(i)+iμ4.2iM(i)+iμ5.1iM(i)\sum\limits_i {\mu _0^iM(i) = } \sum\limits_i {\mu _{3.1}^iM(i)} + \sum\limits_i {\mu _{4.2}^iM(i) + \sum\limits_i {\mu _{5.1}^iM(i)} }

(pour 5 "i" différents car j’ai 5 espèces au total, i = EB, CY, H, ST, H ; M(i) c’est des constantes [masse molaire]).

μ2EB=μ1EB(1XA)\mu _2^{EB} = \mu _1^{EB}(1 - {X_A})

μ2H=μ1H+μ1EBXA3μ1HXB\mu _2^H = \mu _1^H + \mu _1^{EB}{X_A} - 3\mu _1^H{X_B}

(XA et XB sont des constantes connues)

Puis, dans la première unité de séparation [3. FLASH sur le schéma] on enlève "H". Donc, faut imaginer qu’après 3. on a plus que quatre espèces, ie.

μ4.1H=μ4.2H=0\mu _{4.1}^H = \mu _{4.2}^H = 0 (pareil pour le flux 5)

Finalement, on a les équations qui vont faire intervenir "indirectement" la température au niveau des unités de séparation (3, 4 et 5). Voici un exemple et c’est pareil pour toutes les espèces:

μ4.1i=μ3.2iξ4i(T4D,T4B)\mu _{4.1}^i = \mu _{3.2}^i\xi _{4}^i(T_{4}^D,T_{4}^B)

ξ(...)\xi (...) est une valeur qui est comprise entre 0 et 1 (c’est simplement le pourcentage de ce qui va se retrouver en haut de l’unité de séparation pour l’espèce ii). Ce ξ(...)\xi (...) dépend de la volatilé relative et donc de la température. Ces expressions ne sont plus linéaires par contre.

On définit une volatilité relative

αij=αijbasαijhaut=Ps,i(Tbas)Ps,i(Thaut)Ps,j(Tbas)Ps,j(Thaut){\overline \alpha _{ij}} = \sqrt {\alpha _{ij}^{bas}\alpha _{ij}^{haut}} = \sqrt {\frac{{{P_{s,i}}({T^{bas}}){P_{s,i}}({T^{haut}})}}{{{P_{s,j}}({T^{bas}}){P_{s,j}}({T^{haut}})}}}

Ce que je peux encore écrire en fonction des températures sachant que

logPs=ABT+C\log {P_s} = A - \frac{B}{{T + C}}

AA, BB et CC sont des constantes que je connais.

Là, on imagine qu’on a fait une itération… Puis je peux calculer mes compositions (fractions molaires car je connais tous les flux, c’est bêtement relié) et calculer une nouvelle température Tnew=ixiPs,i{T_{new}} = \sum\limits_i {{x_i}{P_{s,i}}} .

L’exercice me donne un critère de convergence T5.1,oldT5.1,new<0.1[K]\left| {{T_{5.1,old}} - {T_{5.1,new}}} \right| < 0.1[K]

Et mes inconnues sont tous les flux μ\mu (sauf pour ST car je connais la quantité que je veux produire).

J’espère que c’est un peu plus clair. Je sais que c’est assez long et pourtant j’ai réduit au maximum le nombre d’équations… Ça m’a pris pas mal de temps à fixer toutes ces équations, maintenant reste à l’implémenter dans Matlab ce qui a l’air plus complexe que je ne le pensais :p

Édité par sotibio

+0 -0

Je vais être honnête, je comprends pas grand chose à ces équations ni à tes notations. Mais si j’essaye de résumer, tu as des équations liant les μ\mu (et non les TT comme je croyais au début…) avec des coefficients qui sont des fonctions des TT, elles même fonctions des μ\mu. Donc ton système n’est pas linéaire, en fait, mais de la forme M(μ)μ=RM(\mu)\mu = R.

Si je comprends bien, c’est dans un cadre scolaire puisque tu parles d’exercice. Vu la façon dont tu présentes les choses, j’ai l’impression que ton système se comporte de façon sympathique et la solution μ\mu est un point fixe de la suite définit par récurrence par M(μn)μn+1=RnM(\mu_n)\mu_{n+1} = R_n. Donc probablement que ce qui est attendu est juste de construire la suite des μn\mu_n à coup de linsolve jusqu’à ce que T(μn+1)T(μn)T(\mu_{n+1}) - T(\mu_n) converge à la précision souhaitée.

I don’t mind that you think slowly, but I do mind that you are publishing faster. — W. Pauli

+2 -0
Auteur du sujet

Exactement! Désolé si c’était très ambigu. Je vais tenter de faire un linslove avec la documentation MATLAB et revenir ici si j’ai des soucis. Mais du coup, si c’est un linsolve c’est bien un système linéaire, non?

Est-ce qu’un fsolve ferait l’affaire aussi ? C’est simplement pour éviter de mettre les équations sous forme de matrice car ça risque d’être compliqué vu le nombre d’équations et de paramètres.

Édité par sotibio

+0 -0

Cette réponse a aidé l’auteur du sujet

Mais du coup, si c’est un linsolve c’est bien un système linéaire, non?

Dans l’expression M(μn)μn+1=RnM(\mu_n)\mu_{n+1}=R_n, μn+1\mu_{n+1} est solution d’un système linéaire. Ça n’empêche que le μ\mu que tu cherches lui, est solution d’une équation non linéaire.

Est-ce qu’un fsolve ferait l’affaire aussi ? C’est simplement pour éviter de mettre les équations sous forme de matrice car ça risque d’être compliqué vu le nombre d’équations et de paramètres.

Quitte à utiliser fsolve, t’aurait autant de mettre ton équation sur μ\mu plutôt que te t’embêter à passer par la suite des μn\mu_n. Il serait par ailleurs absurde de passer par fsolve pour résoudre un système linéaire. Ce sera pas plus compliqué de passer par une matrice de toute façon.

I don’t mind that you think slowly, but I do mind that you are publishing faster. — W. Pauli

+1 -0
Auteur du sujet

Merci @adri1.

Je voyais pas trop comment créer une matrice avec toutes ces variables mais j’ai écrit le code avec fsolvefsolve. Il manque encore des données (les pressions de saturation que je dois chercher avec les coefficients d’Antoine) et le calcul à la fin de mes fractions molaires. Cependant, je vois pas du tout comment faire une boucle / résoudre de manière itérative ce système.

Je laisse quand même toutes les équations pour les curieux même si elles sont pas vraiment utile pour répondre à la question (je suppose). J’ai 45 équations et 45 variables à trouver. Mes paramètres qui dépendent de TT sont les ξ\xi (enfin, indirectement car ça dépend de la volatilité relative qui elle dépend de la pression saturante donc de TT). J’ai plusieurs températures (4 car deux colonnes et 2 températures par colonnes) mais l’exercice me demande de mettre comme critère de convergence uniquement T5.1,oldT5.1,new<0.1\left| {{T_{5.1,old}} - {T_{5.1,new}}} \right| < 0.1 ce qui me perturbe aussi…

J’espère que c’est pas trop bordélique mais tout commentaire est bon à prendre pour que je m’améliore :-)

function fval = root11d(x)
%B = Benzene, ST = Styrene, EB = Ethylbenzene, H = Hydrogen, CY =
%Cyclohexane

%Molar masses in kg/mol
MW_B = (78.11)/1000; 
MW_ST = (104.15)/1000;
MW_EB = (106.17)/1000;
MW_H = (2.01)/1000;
MW_CY = (84.16)/1000;

%Conversion rates of rxn 

X_A = 0.5;
X_B = 0.95;

%Splitting ratio xi 

%Unit 4
xi_4_lk = 0.99; %EB
xi_4_hk = 0.01; %ST 
P_s_4_lk_B = 10^(-;%Antoine Eq for saturation pressure (@ bottom and top of the col)
P_s_4_lk_V = ;
P_s_4_hk_B = ;
P_s_4_hk_V = ;
P_s_4_B_B = ; 
P_s_4_B_V = ; 
P_s_4_CY_B = ;
P_s_4_CY_V =;
P_4 = 0.05; %pressure of the col 4 [bar]

alpha_4_lkhk = ((P_s_4_lk_B*P_s_4_lk_V)/(P_s_4_hk_B*P_s_4_hk_V))^(1/2) % relative volatility between lk and hk (geometric mean). B = Bottom of the col , V = Vapour (ie top of the col)
N_4_min = log(xi_4_lk*(1-xi_4_hk)/xi_4_hk*(1-xi_4_lk))/log(alpha_4_lkhk); %Min number of stages (to calculate splitting ratios of the other comp.)
alpha_4_Bhk = ((P_s_4_B_B*P_s_4_B_V)/(P_s_4_hk_B*P_s_4_hk_V))^(1/2);% rel. vol. between B and ST (geom. mean)
alpha_4_CYhk = ((P_s_4_CY_B*P_s_4_CY_V)/(P_s_4_hk_B*P_s_4_hk_V))^(1/2):% rel. vol. between CY and ST (geom. mean)
xi_4_B = ((alpha_4_Bhk)^(N_4_min)*xi_4_hk)/(1+(alpha_4_Bhk-1)*xi_4_hk);
xi_4_ST = 0.01
xi_4_EB = 0.99
xi_4_CY = ((alpha_4_CYhk)^(N_4_min)*xi_4_hk)/(1+(alpha_4_CYhk-1)*xi_4_hk);

%Unit 5

xi_5_lk = 0.99; %CY
xi_5_hk = 0.01; %EB
P_s_5_lk_B = ;%Antoine Eq for saturation pressure (@ bottom and top of the col)
P_s_5_lk_V = ;
P_s_5_hk_B = ;
P_s_5_hk_V = ;
P_s_5_B_B = ; 
P_s_5_B_V = ; 
P_s_5_ST_B = ;
P_s_5_ST_V =;
P_5 = 1; %pressure of the col 5 [bar]

%Variables Def / Notation: u_flow_spcies, eg. u_0_B is the flow rate "0" of B (benzene)
u_0_B = x(1); 
u_0_ST = x(2); 
u_0_EB = x(3); 
u_0_H = x(4); 
u_0_CY = x(5);

u_1_B = x(6); 
u_1_ST = x(7); 
u_1_EB = x(8); 
u_1_H = x(9); 
u_1_CY = x(10);

u_2_B = x(11); 
u_2_ST = x(12); 
u_2_EB = x(13); 
u_2_H = x(14); 
u_2_CY = x(15);

u_31_B = x(16); 
u_31_ST = x(17); 
u_31_EB = x(18); 
u_31_H = x(19); 
u_31_CY = x(20);

u_32_B = x(21); 
u_32_ST = x(22); 
u_32_EB = x(23); 
u_32_H = x(24); 
u_32_CY = x(25);

u_41_B = x(26); 
u_41_ST = x(27); 
u_41_EB = x(28); 
u_41_H = x(29); 
u_41_CY = x(30);

u_42_B = x(31); 
u_42_ST = x(32); 
u_42_EB = x(33); 
u_42_H = x(34); 
u_42_CY = x(35);

u_51_B = x(36); 
u_51_ST = x(37); 
u_51_EB = x(38); 
u_51_H = x(39); 
u_51_CY = x(40);

u_52_B = x(41); 
u_52_ST = x(42); 
u_52_EB = x(43); 
u_52_H = x(44); 
u_52_CY = x(45);


%Define sys of eq 

fval(1,1) = 5*u_0_B - mu_0_EB;
fval(2,1) = u_42_B*MW_B + u_32_B*MW_B + u_51_B*MW_B - u_0_B*MW_B ...
    + u_42_ST*MW_ST + u_32_ST*MW_ST + u_51_ST*MW_ST - u_0_ST*MW_ST ...
    + u_42_EB*MW_EB + u_32_EB*MW_EB + u_51_EB*MW_EB - u_0_EB*MW_EB ...
    + u_42_H*MW_H + u_32_H*MW_H + u_51_H*MW_H - u_0_H*MW_H ...
    + u_42_CY*MW_CY + u_32_CY*MW_CY + u_51_CY*MW_CY - u_0_CY*MW_CY;
fval(3, 1) = u_0_B + u_52_B - u_1_B; 
fval(4, 1) = u_0_ST + u_52_ST - u_1_ST;
fval(5, 1) = u_0_EB + u_52_EB - u_1_EB;
fval(6, 1) = u_0_H + u_52_H - u_1_H;
fval(7, 1) = u_0_CY + u_52_CY - u_1_CY;
fval(8, 1) = u_1_EB*(1-X_A)-u_2_EB;
fval(9, 1) = u_1_ST + u_1_EB*X_A -u_2_ST;
fval(10, 1) = u_1_B*X_B-u_2_CY+u_1_H;
fval(11, 1) = u_1_EB*X_A -3*u_1_B*X_B -u_2_H + u_1_H; 
%Flash Eq 
fval(12, 1) = u_31_B; 
fval(13, 1) = u_31_ST;
fval(14, 1) = u_31_EB; 
fval(15, 1) = u_31_CY;
fval(16, 1) = u_2_H - u_31_H;
fval(17, 1) = u_2_B - u_3_B; 
fval(18, 1) = u_2_ST - u_3_ST;
fval(19, 1) = u_2_EB - u_3_EB; 
fval(20, 1) = u_2_CY - u_3_CY;
%Distilation Column 1
fval(21, 1) = u_32_B*xi_4_B - u_41_B; %xi to be defined!!
fval(22, 1) = u_32_ST*xi_4_ST - u_41_ST;
fval(23, 1) = u_32_EB*xi_4_EB - u_41_EB;
fval(24, 1) = u_32_CY*xi_4_CY - u_41_CY;
fval(25, 1) = u_41_H; 
fval(26, 1) = u_42_H; 
fval(27, 1) = u_32_B*(1 - xi_4_B) - u_42_B; 
fval(28, 1) = u_32_ST*(1 - xi_4_ST) - u_42_ST - 2000/MW_ST;
fval(29, 1) = u_32_EB*(1 - xi_4_EB) - u_42_EB;
fval(30, 1) = u_32_CY*(1 - xi_4_CY) - u_42_CY;
fval(31, 1) = u_51_H;
fval(32, 1) = u_52_H;
fval(33, 1) = u_41_B*xi_5_B - u_51_B;
fval(34, 1) = u_41_ST*xi_5_ST - u_51_ST;
fval(35, 1) = u_41_EB*xi_5_EB - u_51_EB;
fval(36, 1) = u_41_CY*xi_5_CY - u_51_CY;
fval(37, 1) = u_41_B*(1 - xi_5_B) - u_52_B;
fval(38, 1) = u_41_ST*(1 - xi_5_ST) - u_52_ST;
fval(39, 1) = u_41_EB*(1 - xi_5_EB) - u_52_EB;
fval(40, 1) = u_41_CY*(1 - xi_5_CY) - u_52_CY;
%Feed constraints 
fval(41, 1) = u_O_ST; 
fval(42, 1) = u_0_CY;
fval(43, 1) = u_0_H;
fval(44, 1) = u_0_B + u_52_B - u_1_B;
fval(45, 1) = u_0_EB + u_52_EB - u_1_EB;

%Mole fractions
x_blabla = u_jrrjrk / ejjeor; 
%Calculate new bubble temperature T_51
T_51 = x_blabla*P_s_blabla + ..;
% If old(T_51_new - T_51_old) < 0.1 -> solution is OK
% Else -> Solve again the system of equations with T_51_new

PS: Considérez les équations comme correctes évidemment ^^

Édité par sotibio

+0 -0

Cette réponse a aidé l’auteur du sujet

Je voyais pas trop comment créer une matrice avec toutes ces variables mais j’ai écrit le code avec fsolvefsolve.

Encore une fois, c’est complètement aberrant. Tu utilises un solveur hyper-général pensé pour des équations non linéaires, et tu injectes un problème linéaire dedans. C’est inefficace possible, ça ne simplifie pas le code, et en plus ça rend ton intention beaucoup plus difficile à lire. Personne ne va s’attendre à ce que tu utilises fsolve pour résoudre un système linéaire, et encore moins pour résoudre nn fois un système linéaire purement accessoire pour résoudre un système non linéaire qui aurait pu être injecté directement dans fsolve. C’est un peu comme si tu voulais détruire un immeuble infesté de mouche, et que tu utilises un bulldozer pour tuer les mouches une par une au lieu de détruire directement l’immeuble.

Il manque encore des données (les pressions de saturation que je dois chercher avec les coefficients d’Antoine) et le calcul à la fin de mes fractions molaires. Cependant, je vois pas du tout comment faire une boucle / résoudre de manière itérative ce système.

Les valeurs successives du vecteur μ\mu sont notées μn\mu^n. Le cœur de ton programme va ressembler à ça :

  • tant que T51(μn)T51(μn+1)<0.1|T_{51}(\mu^n)-T_{51}(\mu^{n+1})|<0.1

    • μnμn+1\mu^n\leftarrow \mu^{n+1}
    • calcul de M(μn)M(\mu^n) et R(μn)R(\mu^n)
    • calcul de μn+1\mu^{n+1} solution de M(μn)μn+1=R(μn)M(\mu^n)\mu^{n+1}=R(\mu^n)

Pour le calcul de M(μn)M(\mu^n) et R(μn)R(\mu^n), ça demande un peu d’organisation. Je vais noter avec un indice les différences composantes du vecteur inconnu μ\mu. Disons pour simplifier que tu as trois inconnues, μ1\mu_1, μ2\mu_2 et μ3\mu_3. Tu as trois équations, mettons

μ1+α(μ)μ2=β(μ)μ2+3μ3=γ(μ)μ1δ(μ)μ3=0\begin{aligned} \mu_1+\alpha(\mu)\mu_2 &= \beta(\mu)\\ \mu_2 + 3\mu_3 &= \gamma(\mu)\\ \mu_1 - \delta(\mu)\mu _3 &= 0 \end{aligned}

avec α(μ)\alpha(\mu), β(μ)\beta(\mu), γ(μ)\gamma(\mu) et δ(μ)\delta(\mu) des coefficients qui dépendent de μ\mu.

Pour écrire ça sous forme matricielle, il va falloir mettre une équation par ligne dans la matrice. Si je reprends le système d’au dessus et que je rajoute les exposants dénotant les itérations, ça va ressembler à ça :

(1α(μn)001310δ(μn))μn+1=(β(μn)γ(μn)0)\begin{pmatrix} 1 & \alpha(\mu^n) & 0\\ 0 & 1 & 3 \\ 1 & 0 & -\delta(\mu^n) \end{pmatrix}\mu^{n+1}= \begin{pmatrix} \beta(\mu^n)\\ \gamma(\mu^n)\\ 0 \end{pmatrix}

Et sous forme de code, ça va ressembler à ça (en supposant que MM et RR existent déjà) :

M = 0;
R = 0;

% First equation
M(1, 1) = 1;
M(1, 2) = alpha(mu_old);
R(1) = beta(mu_old);

% Second equation
M(2, 2) = 1;
M(2, 3) = 3;
R(2) = gamma(mu_old);

% Third equation
M(3, 1) = 1;
M(3, 3) = -delta(mu_old);

% Solve for mu_new
mu_new = linsolve(M, R);

Tout ça sera bien sûr dans un while qui va répéter la manœuvre jusqu’à convergence.

EDIT : inversion des indices lignes/colonnes dans le code Matlab…

Édité par adri1

I don’t mind that you think slowly, but I do mind that you are publishing faster. — W. Pauli

+1 -0
Auteur du sujet

Salut,

Je comprends tout à faire le concept mais pas tellement comment l’écrire… J’ai transcrit tout ça dans la matrice qui est donc une 45×4545 \times 45. J’ai pas encore mis les température T51T_{51} et T52T_{52} mais bon pour ce qui en est des équations on en est presque là (ça sera pareil que pour T41T_{41} et T42T_{42}). À part ça, c’est sûrement plus optimal d’utiliser linsolve mais ça m’a pris du temps vu la taille de la matrice :p

Ce que je comprends pas c’est quand est-ce que je dois intégrer "mu_old" dans mes équations… Est-ce que c’est au-moment où j’ai une dépendence sur les température ? (donc où j’ai des "xi")

Aussi, je ne vois pas super bien où je pourrais intégrer la boucle whilewhile. Est-ce que je met tout le système ? J’aurai aussi des termes non-linéaires vu que la température apparaît dans les équations d’Antoine pour le calcul des pressions saturantes (P_s dans le code).

Merci pour ton aide!

%B = Benzene, ST = Styrene, EB = Ethylbenzene, H = Hydrogen, CY =
%Cyclohexane


%Molar masses in kg/mol
MW_B = (78.11)/1000; 
MW_ST = (104.15)/1000;
MW_EB = (106.17)/1000;
MW_H = (2.01)/1000;
MW_CY = (84.16)/1000;

%Conversion rates of rxn 

X_A = 0.5;
X_B = 0.95;

%Antoine parameters 

A_B = 4.73;
B_B = 1316.55;
C_B = - 35.58;
A_EB = 4.08;
B_EB = 1419.30;
C_EB = -60.54;
A_ST = 4.22; 
B_ST = 1525.06; 
C_ST = -56.40;
A_CY = 4.14;
B_CY = 1316.55;
C_CY = - 35.58;

%Splitting ratio xi

%Unit 4
xi_4_lk = 0.99; %EB
xi_4_hk = 0.01; %ST 
P_s_4_lk_B = 10^(-(A_EB - B_EB/(T_42 + C_EB))) ;%Antoine Eq for saturation pressure (@ bottom and top of the col)
P_s_4_lk_V = 10^(-(A_EB - B_EB/(T_41 + C_EB)));
P_s_4_hk_B = 10^(-(A_ST - B_ST/(T_42 + C_ST)));
P_s_4_hk_V = 10^(-(A_ST - B_ST/(T_41 + C_ST)));
P_s_4_B_B = 10^(-(A_B - B_B/(T_42 + C_B))); 
P_s_4_B_V = 10^(-(A_B - B_B/(T_41 + C_B))); 
P_s_4_CY_B = 10^(-(A_CY - B_CY/(T_42 + C_CY)));
P_s_4_CY_V = 10^(-(A_CY - B_CY/(T_41 + C_CY)));
P_4 = 0.05; %pressure of the col 4 [bar]

%Define K-value (ie. P_sat / P_column)
K_4_lk_B = P_s_4_lk_B / P_4;
K_4_lk_V = P_s_4_lk_V / P_4;
K_4_hk_B = P_s_4_hk_B / P_4;
K_4_hk_V = P_s_4_hk_V / P_4;
K_4_B_B = P_s_4_B_B / P_4;
K_4_B_V = P_s_4_B_V / P_4;
K_4_CY_B = P_s_4_B_CY_B / P_4;
K_4_CY_V = P_s_4_CY_V / P_4;


alpha_4_lkhk = (K_4_lk_B*K_4_lk_V/(K_4_hk_B*K_4_hk_V))^(1/2) % relative volatility between lk and hk (geometric mean). B = Bottom of the col , V = Vapour (ie top of the col)
N_4_min = log(xi_4_lk*(1-xi_4_hk)/xi_4_hk*(1-xi_4_lk))/log(alpha_4_lkhk); %Min number of stages (to calculate splitting ratios of the other comp.)
alpha_4_Bhk = ((K_4_B_B*K_4_B_V)/(K_4_hk_B*K_4_hk_V))^(1/2);% rel. vol. between B and ST (geom. mean)
alpha_4_CYhk = ((K_4_CY_B*K_4_CY_V)/(K_4_hk_B*K_4_hk_V))^(1/2):% rel. vol. between CY and ST (geom. mean)
xi_4_B = ((alpha_4_Bhk)^(N_4_min)*xi_4_hk)/(1+(alpha_4_Bhk-1)*xi_4_hk);
xi_4_ST = 0.01
xi_4_EB = 0.99
xi_4_CY = ((alpha_4_CYhk)^(N_4_min)*xi_4_hk)/(1+(alpha_4_CYhk-1)*xi_4_hk);


M = 0;
R = 0;

% First equation
M(1, 1) = 5;
M(1, 2) = -1;

% Second equation
M(2, 1) = -MW_B;
M(2, 2) = -MW_EB;
M(2, 3) = -MW_ST;
M(2, 4) = -MW_CY;
M(2, 5) = -MW_H;
M(2, 21) = MW_B;
M(2, 22) = MW_EB;
M(2, 23) = MW_ST;
M(2, 24) = MW_CY;
M(2, 25) = MW_H;
M(2, 31) = MW_B;
M(2, 32) = MW_EB;
M(2, 33) = MW_ST;
M(2, 34) = MW_CY;
M(2, 35) = MW_H;
M(2, 36) = MW_B;
M(2, 37) = MW_EB;
M(2, 38) = MW_ST;
M(2, 39) = MW_CY;
M(2, 40) = MW_H;

% Third equation
M(3, 1) = 1;
M(3, 6) = -1;
M(3, 36) = 1; 

% Fourth Eq

M(4, 2) = 1;
M(4, 7) = -1;
M(4, 37) = 1; 

% Fifth Eq

M(5, 3) = 1;
M(5, 8) = -1;
M(5, 38) = 1; 

% Eq 6

M(6, 4) = 1;
M(6, 9) = -1;
M(6, 39) = 1;

% Eq 7

M(7, 5) = 1;
M(7, 10) = -1;
M(7, 40) = 1;

% Eq 8 

M(8, 7) = (1-X_A); 
M(8, 12) = -1; 

% Eq 9

M(9, 7) = X_A;
M(9, 8) = 1; 
M(9, 13) = -1;

% Eq 10

M(10, 6) = X_B;
M(10, 10) = 1;
M(10, 14) = -1;

% Eq 11

M(11, 6) = -3*X_B;
M(11, 7) = X_A; 
M(11, 10) = 1;
M(11, 15) = -1;

% Eq 12

M(12, 16) = -1; 

% Eq 13

M(13, 17) = -1;

% Eq 14

M(14, 18) = -1;

% Eq 15

M(15, 19) = -1:

% Eq 16

M(16, 15) = 1;
M(16, 20) = -1;

% Eq 17

M(17, 11) = 1; 
M(17, 21) = -1;

% Eq 18

M(18, 12) = 1; 
M(18, 22) = -1;

% Eq 19

M(19, 13) = 1; 
M(19, 23) = -1;


% Eq 20

M(20, 14) = 1; 
M(20, 24) = -1;

% Eq 21

M(21, 21) = xi_4_B;
M(21, 26) = -1;

% Eq 22

M(22, 22) = xi_4_EB;
M(22, 27) = -1;


% Eq 23

M(23, 23) = xi_4_ST;
M(23, 28) = -1;


% Eq 24

M(24, 24) = xi_4_CY;
M(24, 29) = -1;

% Eq 25

M(25, 30) = 1; 

% Eq 26 

M(26, 35) = 1;

% Eq 27

M(27, 21) = (1 - xi_4_B);
M(27, 31) = -1; 

% Eq 28

M(28, 22) = (1 - xi_4_EB);
M(28, 32) = -1; 

% Eq 29

M(29, 23) = (1 - xi_4_ST);
M(29, 33) = -1;
R(29) = 2000/MW_ST;

% Eq 30

M(30, 24) = (1 - xi_4_CY);
M(30, 34) = -1;

% Eq 31 & 32  

M(31, 40) = 1; 
M(32, 45) = 1;

% Eq 33

M(33, 26) = xi_5_B;
M(33, 36) = -1;

% Eq 34

M(34, 27) = xi_5_EB;
M(34, 37) = -1;

% Eq 35

M(35, 28) = xi_5_ST;
M(35, 38) = -1;

% Eq 36

M(36, 29) = xi_5_CY;
M(36, 39) = -1;

% Eq 37

M(37, 26) = (1 - xi_5_B);
M(37, 41) = -1;

% Eq 38

M(38, 27) = (1 - xi_5_EB);
M(38, 42) = -1;

% Eq 39

M(39, 28) = (1 - xi_5_ST);
M(39, 43) = -1;


% Eq 40

M(40, 29) = (1 - xi_5_CY);
M(40, 44) = -1;

% Eq 41, 42, 43, 44 & 45 - Feed constraints 

M(41, 3) = 1;
M(42, 4) = 1; 
M(43, 5) = 1;
M(44, 1) = 1;
M(44, 6) = -1;
M(44, 41) = 1;
M(45, 2) = 1;
M(45, 7) = -1;
M(45, 42) = 1;

% Solve for mu_new
mu_new = linsolve(M, R);
+0 -0

Je comprends tout à faire le concept […] Ce que je comprends pas c’est quand est-ce que je dois intégrer "mu_old" dans mes équations… Est-ce que c’est au-moment où j’ai une dépendence sur les température ? (donc où j’ai des "xi")

Aussi, je ne vois pas super bien où je pourrais intégrer la boucle while. Est-ce que je met tout le système ? J’aurai aussi des termes non-linéaires vu que la température apparaît dans les équations d’Antoine pour le calcul des pressions saturantes (P_s dans le code).

Ces questions me font dire que tu n’as pas compris l’idée derrière l’algorithme. mu_old va en effet apparaître à travers certains des coefficients (les xi qui dépendent des T qui eux même dépendent de mu).

Pour la boucle while, je l’ai écrit dans ce message. Le but, c’est d’itérer tant que les températures n’ont pas convergées. Et oui tu vas avoir des termes non linéaires en mu_old dans la matrice, c’est le but de toute la méthode : transformer une équation non linéaire en une limite d’un processus itératif linéaire.

I don’t mind that you think slowly, but I do mind that you are publishing faster. — W. Pauli

+1 -0
Auteur du sujet

Si je reprends le code dans mon message précédent, j’aurai donc un while partir de la ligne 33 jusqu’à la ligne 300? Boucle de la forme:

while abs(T51(mu_old) - T51(mu_new)) <0.1

end

Et par exemple, à la ligne 194 où j’ai une dépendance sur T, j’aurai xi_5_EB(mu_old) ?

+0 -0

Cette réponse a aidé l’auteur du sujet

Si je reprends le code dans mon message précédent, j’aurai donc un while partir de la ligne 33 jusqu’à la ligne 300? Boucle de la forme:

Jusqu’à la toute fin, même, sinon ça sert pas à grand chose si mu_new et mu_old ne sont pas mis à jour au fur et à mesure.

Et par exemple, à la ligne 194 où j’ai une dépendance sur T, j’aurai xi_5_EB(mu_old) ?

Oui. Enfin, ce sera peut être pas écrit comme ça exactement mais c’est l’idée : la valeur de xi_5_EB va dépendre de mu_old d’une manière ou d’une autre.

I don’t mind that you think slowly, but I do mind that you are publishing faster. — W. Pauli

+1 -0
Auteur du sujet

D’accord. Par contre, je vois pas comment est-ce que je peux récupérer une solution.. Est-ce que c’est par exemple munew(2) qui va me renvoyer ${\mu {0,EB}}$ ?

Je demande car j’en ai besoin pour calculer les fractions molaire que je calculerais donc juste après la ligne mu_new = linsolve(M, R). Et du coup, ma boucle ira jusqu’après le calcul des nouvelles températures (ie. juste après les fractions molaires…)

Par contre, j’avoue que je vois pas du tout comment je pourrais écrire la valeur de xi_5_EB en fonction de mu_old avec du code…

+0 -0

Est-ce que c’est par exemple munew(2) qui va me renvoyer μ0,EB{\mu {0,EB}} ?

Ben ça je peux pas le deviner à ta place, je ne sais pas à quelle position dans le vecteur μ\mu tu as mis μ0,EB{\mu {0,EB}}

Je demande car j’en ai besoin pour calculer les fractions molaire que je calculerais donc juste après la ligne mu_new = linsolve(M, R). Et du coup, ma boucle ira jusqu’après le calcul des nouvelles températures (ie. juste après les fractions molaires…)

Hmmm, t’aurais pas intérêt à mettre ces calculs là au début de la boucle plutôt qu’à la fin ? Ça t’évite d’avoir à les initialiser.

Par contre, j’avoue que je vois pas du tout comment je pourrais écrire la valeur de xi_5_EB en fonction de mu_old avec du code…

Ben… Ça aussi c’est pas un truc que je peux deviner, c’est toi qui est censé savoir comment ces coefficients se calculent. C’est toi qui m’avait dit que pour une valeur de μ\mu, tu pouvais calculer les TT et en déduire les ξ\xi. Ben il suffit d’écrire juste ça.

Au fait, je viens de voir mais là :

while abs(T51(mu_old) - T51(mu_new)) <0.1

ta condition est dans le mauvais sens. while, ça veut dire tant que, pas jusqu’à ce que.

Édité par adri1

I don’t mind that you think slowly, but I do mind that you are publishing faster. — W. Pauli

+0 -0
Auteur du sujet

Est-ce que c’est par exemple munew(2) qui va me renvoyer μ0,EB{\mu {0,EB}} ?

Ben ça je peux pas le deviner à ta place, je ne sais pas à quelle position dans le vecteur μ\mu tu as mis μ0,EB{\mu {0,EB}}

Oui oui, je veux dire du point de vue conceptuel. En Matlab, quand c’est un vecteur à une dimension il suffit de faire mu_new(2) et pas mu_new(2,1) ou pas?

Hmmm, t’aurais pas intérêt à mettre ces calculs là au début de la boucle plutôt qu’à la fin ? Ça t’évite d’avoir à les initialiser.

Sûrement :p Le soucis si je les mets au début c’est que la première fois que la boucle va tourner elle sera perdue car elle aura pas encore de valeurs pour mes mu_new ?

Ben… Ça aussi c’est pas un truc que je peux deviner, c’est toi qui est censé savoir comment ces coefficients se calculent. C’est toi qui m’avait dit que pour une valeur de μ\mu, tu pouvais calculer les TT et en déduire les ξ\xi. Ben il suffit d’écrire juste ça.

Ah, ça je l’ai fait avant de lister toutes les équations (lignes 35–60 pour l’unité 4 avant d’initialiser M et R à zéro):

%Splitting ratio xi


%Unit 4
xi_4_lk = 0.99; %EB
xi_4_hk = 0.01; %ST 
P_s_4_lk_B = 10^(-(A_EB - B_EB/(T_42 + C_EB))) ;%Antoine Eq for saturation pressure (@ bottom and top of the col)
P_s_4_lk_V = 10^(-(A_EB - B_EB/(T_41 + C_EB)));
P_s_4_hk_B = 10^(-(A_ST - B_ST/(T_42 + C_ST)));
P_s_4_hk_V = 10^(-(A_ST - B_ST/(T_41 + C_ST)));
P_s_4_B_B = 10^(-(A_B - B_B/(T_42 + C_B))); 
P_s_4_B_V = 10^(-(A_B - B_B/(T_41 + C_B))); 
P_s_4_CY_B = 10^(-(A_CY - B_CY/(T_42 + C_CY)));
P_s_4_CY_V = 10^(-(A_CY - B_CY/(T_41 + C_CY)));
P_4 = 0.05; %pressure of the col 4 [bar]

%Define K-value (ie. P_sat / P_column)
K_4_lk_B = P_s_4_lk_B / P_4;
K_4_lk_V = P_s_4_lk_V / P_4;
K_4_hk_B = P_s_4_hk_B / P_4;
K_4_hk_V = P_s_4_hk_V / P_4;
K_4_B_B = P_s_4_B_B / P_4;
K_4_B_V = P_s_4_B_V / P_4;
K_4_CY_B = P_s_4_CY_B / P_4;
K_4_CY_V = P_s_4_CY_V / P_4;


alpha_4_lkhk = (K_4_lk_B*K_4_lk_V/(K_4_hk_B*K_4_hk_V))^(1/2) % relative volatility between lk and hk (geometric mean). B = Bottom of the col , V = Vapour (ie top of the col)
N_4_min = log(xi_4_lk*(1-xi_4_hk)/xi_4_hk*(1-xi_4_lk))/log(alpha_4_lkhk); %Min number of stages (to calculate splitting ratios of the other comp.)
alpha_4_Bhk = ((K_4_B_B*K_4_B_V)/(K_4_hk_B*K_4_hk_V))^(1/2);% rel. vol. between B and ST (geom. mean)
alpha_4_CYhk = ((K_4_CY_B*K_4_CY_V)/(K_4_hk_B*K_4_hk_V))^(1/2):% rel. vol. between CY and ST (geom. mean)
xi_4_B = ((alpha_4_Bhk)^(N_4_min)*xi_4_hk)/(1+(alpha_4_Bhk-1)*xi_4_hk);
xi_4_ST = 0.01
xi_4_EB = 0.99
xi_4_CY = ((alpha_4_CYhk)^(N_4_min)*xi_4_hk)/(1+(alpha_4_CYhk-1)*xi_4_hk);


%Unit 5
xi_5_lk = 0.99; %CY
xi_5_hk = 0.01; %EB 
P_s_5_lk_B = 10^(-(A_CY - B_CY/(T_52 + C_CY))) ;%Antoine Eq for saturation pressure (@ bottom and top of the col)
P_s_5_lk_V = 10^(-(A_CY - B_CY/(T_51 + C_CY)));
P_s_5_hk_B = 10^(-(A_EB - B_EB/(T_52 + C_EB)));
P_s_5_hk_V = 10^(-(A_EB - B_EB/(T_51 + C_EB)));
P_s_5_B_B = 10^(-(A_B - B_B/(T_52 + C_B))); 
P_s_5_B_V = 10^(-(A_B - B_B/(T_51 + C_B))); 
P_s_5_ST_B = 10^(-(A_ST - B_ST/(T_52 + C_ST)));
P_s_5_ST_V = 10^(-(A_ST - B_ST/(T_51 + C_ST)));
P_5 = 1; %pressure of the col 4 [bar] 

%Define K-value (ie. P_sat / P_column)
K_5_lk_B = P_s_5_lk_B / P_5;
K_5_lk_V = P_s_5_lk_V / P_5;
K_5_hk_B = P_s_5_hk_B / P_5;
K_5_hk_V = P_s_5_hk_V / P_5;
K_5_B_B = P_s_5_B_B / P_5;
K_5_B_V = P_s_5_B_V / P_5;
K_5_ST_B = P_s_5_ST_B / P_5;
K_5_ST_V = P_s_5_ST_V / P_5;


alpha_5_lkhk = (K_5_lk_B*K_5_lk_V/(K_5_hk_B*K_5_hk_V))^(1/2) % relative volatility between lk and hk (geometric mean). B = Bottom of the col , V = Vapour (ie top of the col)
N_5_min = log(xi_5_lk*(1-xi_5_hk)/xi_5_hk*(1-xi_5_lk))/log(alpha_5_lkhk); %Min number of stages (to calculate splitting ratios of the other comp.)
alpha_4_Bhk = ((K_5_B_B*K_5_B_V)/(K_5_hk_B*K_5_hk_V))^(1/2);% rel. vol. between B and ST (geom. mean)
alpha_5_SThk = ((K_5_ST_B*K_5_ST_V)/(K_5_hk_B*K_5_hk_V))^(1/2):% rel. vol. between CY and ST (geom. mean)
xi_5_B = ((alpha_5_Bhk)^(N_5_min)*xi_5_hk)/(1+(alpha_5_Bhk -1)*xi_5_hk);
xi_5_EB = 0.01
xi_5_CY = 0.99
xi_5_ST = ((alpha_5_SThk)^(N_5_min)*xi_5_hk)/(1+(alpha_5_SThk -1)*xi_5_hk);

Là, j’ai déjà défini ce qu’était ces ξ\xi. Mon but juste avant ces quelques lignes est de donner des températures (notées T_52, T_51, … dans le code) approximatives pour les faire converger par la suite vers les températures réelles. Mais dans tout ça j’ai jamais dans mon code un "mu_old" ce qui me perturbe énormément pour l’écriture de la boucle.

Au fait, je viens de voir mais là :

while abs(T51(mu_old) - T51(mu_new)) <0.1

ta condition est dans le mauvais sens. while, ça veut dire tant que, pas jusqu’à ce que.

adri1

Ah merci! Je dois donc changer le sens de l’inégalité :-)

+0 -0

Le soucis si je les mets au début c’est que la première fois que la boucle va tourner elle sera perdue car elle aura pas encore de valeurs pour mes mu_new ?

Rien compris, lors du premier tour de boucle, mu_new n’a pas encore été calculé une seule fois, comment diable compte tu te servir de sa valeur ? Au deuxième tour de boucle, le mu_old sera le mu_new du tour précédent et sera utilisé à ce moment là.

Ah, ça je l’ai fait avant de lister toutes les équations (lignes 35–60 pour l’unité 4 avant d’initialiser M et R à zéro):

Elles sortent d’où T_41 et T_42 ? Et les valeurs des différents xi ? De ce que j’avais compris, c’était ça qui dépendait de μ\mu (et donc dans ta boucle, de mu_old). Si ça dépend pas de μ\mu, alors j’ai absolument rien compris à ton problème depuis le départ et j’ai joyeusement écrit des M(μ)M(\mu) alors que ta matrice ne dépend pas de μ\mu:-°

I don’t mind that you think slowly, but I do mind that you are publishing faster. — W. Pauli

+0 -0
Auteur du sujet

Rien compris, lors du premier tour de boucle, mu_new n’a pas encore été calculé une seule fois, comment diable compte tu te servir de sa valeur ? Au deuxième tour de boucle, le mu_old sera le mu_new du tour précédent et sera utilisé à ce moment là.

Mais… alors, faut bien que je le mette à la fin vu que j’ai aucune idée des "mu_old"… J’aurai une idée de leur valeur uniquement après la première itération.

Elles sortent d’où T_41 et T_42 ? Et les valeurs des différents xi ? De ce que j’avais compris, c’était ça qui dépendait de μ\mu (et donc dans ta boucle, de mu_old). Si ça dépend pas de μ\mu, alors j’ai absolument rien compris à ton problème depuis le départ et j’ai joyeusement écrit des M(μ)M(\mu) alors que ta matrice ne dépend pas de μ\mu:-°

adri1

T41 et T_42 (idem pour T_51 et T_52) sont des températures de ces fameuses unités de séparation. C’est pour ces paramètres que je "devine" les valeurs… Les ξi\xi _i (xi) ce sont des formules de mon cours pour ces unités de séparation. Ces ξi\xi _i ne dépendent pas des flux μ\mu mais uniquement du paramètre $N{min}$ (calculé juste avant par une autre formule) et de la volatilité relative (ie. de la température -> T_41 et T_42 pour l’unité 4)

Pour répéter, tous les μ\mu sont mes inconnues mais les quatre températures TT le sont aussi (même si je les devine au début, à la fin elles auront des valeurs différents en fonction de la convergence).

J’espère que c’est plus clair…

+0 -0

Mais… alors, faut bien que je le mette à la fin vu que j’ai aucune idée des "mu_old"… J’aurai une idée de leur valeur uniquement après la première itération.

Ben… Non, tu te donnes un mu_old pour commencer l’itération. Enfin, vu ce que tu dis ensuite, j’ai probablement toujours pas compris correctement le problème.

T_41 et T_42 (idem pour T_51 et T_52) sont des températures de ces fameuses unités de séparation. C’est pour ces paramètres que je "devine" les valeurs… Les ξi\xi _i (xi) ce sont des formules de mon cours pour ces unités de séparation.

Je sais pas ce que sont ces "fameuses unités de séparation". Encore une fois, je comprends pas ton problème chimique, mais si tu me donnes un problème mathématiquement bien posé on arrivera à avancer. Quand tu parles de "formules" dans tous les sens, ça m’avance pas parce que ça me dit pas qui dépend de qui. J’ai juste besoin de savoir qui est fonction de quoi dans ton problème, et si ces fonctions sont linéaires ou non.

Ces ξi\xi _i ne dépendent pas des flux μ\mu mais uniquement du paramètre NminN_{min} (calculé juste avant par une autre formule) et de la volatilité relative (ie. de la température -> T_41 et T_42 pour l’unité 4)

Pareil, tu parles de formule mais sans me dire qui dépend de quoi je suis dans le flou le plus complet.

Pour répéter, tous les μ\mu sont mes inconnues mais les quatre températures TT le sont aussi (même si je les devine au début, à la fin elles auront des valeurs différents en fonction de la convergence).

Comme t’as dit trois trucs différents au moins, je me suis un peu perdu en route.

Si j’essaye de re-résumer, tu as des flux (un vecteur de μ\mu) inconnus reliés par plein d’équations qui peuvent se mettre sous forme matricielle. Les coefficients de cette matrice dépendent des coefficients ξ\xi qui je comprends bien sont eux-même des fonctions des températures qui sont aussi des inconnues du problème. Jusque là j’ai bon, ou il me manque un truc ?

Maintenant, ce qui me pose question, c’est comment tu calcules ces températures, de quelles autres variables elles dépendent. Sans lien entre les températures et les flux je vois mal comment tu comptes t’en sortir.

I don’t mind that you think slowly, but I do mind that you are publishing faster. — W. Pauli

+0 -0
Auteur du sujet

Ben… Non, tu te donnes un mu_old pour commencer l’itération. Enfin, vu ce que tu dis ensuite, j’ai probablement toujours pas compris correctement le problème.

Je peux mettre n’importe quoi ? Je peux pas initialiser à zéro car je vais devoir faire un ratio pour obtenir une fraction (molaire).

Je sais pas ce que sont ces "fameuses unités de séparation". Encore une fois, je comprends pas ton problème chimique, mais si tu me donnes un problème mathématiquement bien posé on arrivera à avancer. Quand tu parles de "formules" dans tous les sens, ça m’avance pas parce que ça me dit pas qui dépend de qui. J’ai juste besoin de savoir qui est fonction de quoi dans ton problème, et si ces fonctions sont linéaires ou non.

OK. Désolé. Du coup, sans parler chimie mais simplement mathématiquement voici:

  • Toutes les équations que j’ai écris dans le code M(.., ..) sont linéaires et sont simplement des conservations de masse;

  • Cependant, dans ces équation on retrouve un paramètre qu’il faut évaluer avec des équations non-linéaire. C’est les fameux ξ\xi (y en a plusieurs mais ils sont évalués avec les mêmes équations - simplement les valeurs changent).

  • Ces fameuses équations pour trouver ξ\xi sont les suivantes:

Étape 1: on cherche un paramètre appelé N_min

Nmin=ln[ξlk(1ξhk)ξhk(1ξlk)]/lnαlk,hk{N_{\min }} = \ln \left[ {\frac{{{\xi _{lk}}(1 - {\xi _{hk}})}}{{{\xi _{hk}}(1 - {\xi _{lk}})}}} \right]/\ln {\alpha _{lk,hk}}

αij=αijBαijV=Ps,i(TB)Ps,i(TV)Ps,j(TB)Ps,j(TV){\overline \alpha _{ij}} = \sqrt {\alpha _{ij}^B\alpha _{ij}^V} = \sqrt {\frac{{{P_{s,i}}({T^B}){P_{s,i}}({T^V})}}{{{P_{s,j}}({T^B}){P_{s,j}}({T^V})}}}.

On évalue les Ps,iP_{s,i} avec une autre équation non-linéaire qui est:

Ps=10(ABC+T){P_s} = {10^{ - (A - \frac{B}{{C + T}})}}

où la température apparaît!

Étape 2: on calcule le ξ\xi correspondant avec l’équation non-linéaire suivante:

ξj=αj,hkNminξhk/[1+(αj,hkNmin1)ξhk]{\xi _j} = \alpha _{j,hk}^{{N_{\min }}}{\xi _{hk}}/\left[ {1 + (\alpha _{j,hk}^{{N_{\min }}} - 1){\xi _{hk}}} \right]

ξhk\xi_{hk} est une constante (écrit dans le code)

Maintenant, ce qui me pose question, c’est comment tu calcules ces températures, de quelles autres variables elles dépendent. Sans lien entre les températures et les flux je vois mal comment tu comptes t’en sortir.

adri1

C’est la que j’ai besoin des fractions molaires (donc des mu). Chaque nouvelle température sera calculée à partir de là avec une relation (que j’ai pas écrite car j’en suis pas encore sûr.. je dois demander vérification).

J’espère que c’est (finalement) clair!

Édité par sotibio

+0 -0
Vous devez être connecté pour pouvoir poster un message.
Connexion

Pas encore inscrit ?

Créez un compte en une minute pour profiter pleinement de toutes les fonctionnalités de Zeste de Savoir. Ici, tout est gratuit et sans publicité.
Créer un compte