Méthode de Résolution

Dans cette partie on présentera  la discrétisation des équations qu'on cherche à résoudre ainsi que les schémas numériques utilisés pour la résolution, par des schémas on illustreras la dépendance des codes et leur implémentation ainsi que les méthodes de calcul utilisées et la gestion des différentes routines et leur communication.

En effet, les équations obtenues précédemment sont non linéaires, les domaines seront discrétisés en une multitude de tranche, il convient alors de faire le choix judicieux du type de la discrétisation ainsi qu'au schémas de résolution.

Partie dense

Le code simulateur de la partie dense est effectuée par le biais d'un code impléménté en Fortran 90. La résolution de cette partie réside dans la résolution d'un système de trois équations non linéaires en débits molaires :

\begin{equation}
\left\lbrace
\begin{array}{ccc}
F_{O_{2_e}}-k_c \rho_g \frac{Y_{O_2}S_{pc}\alpha_c V_{rac}}{M_{O_2}}-F_{O_{2_{guess}}}=0\\
F_{C_e}-k_c \rho_g \frac{Y_{O_2}S_{pc}\alpha_c V_{rac}}{M_{O_2}}-F_{C_{guess}}=0\\
\sum_{1}^{5} F_{ie}H_{ie} - F_{C_{O_2}}\Delta H^0_{C_{O_2}} - \sum_{1}^{5} F_{is}H_{is}=0
\end{array}
\right.
\end{equation}

A partir de ces équations qu'on considére comme composantes d'un vecteur on construit un vecteur dont   la norme Euclidienne ce est  à minimiser, qui est équivalent à trouver où il s'annule. L'idée est de partir d'un vecteur contenant les guess : des grandeurs qu'on aura à estimer et  puis refaire les itérations partant à chaque fois d'un nouveau guess, le point le plus difficile dans cette approche est de calculer à chaque fois le guess pour se rapprocher le plus vite possible de la solution.

L'algorithme utilisé dans ce cas est celui de Levenberg-Marquardt, vu que c'est un algorithme plus stable que celui de Gauss Newton (celui utilisé avant), plus rapide et arrive à trouver dans la plupart des cas la solution.

L'algorithme de Levenberg-Marquardt, souvenet appelé algorithme LM, permet la résolution d'un système non linéaire de plusieurs paramètres à travers la minimisation de la fonction qui aux inconnues du système associe les équations formant le système. Il repose sur l'interpolation de l'algorithme de Gauss Newton, et de l'algorithme du gradient. Cet interpolation des deux dernières algorithmes assure une stabilité relativement inconditionnelle par rapport à celle associé à l'algorithme de Gauss newton, puisque en se soucie moins des estimations grossière de la solution,car on arrive à avoir des solutions avec une convergence rapidement atteinte.

La procédure de l'algorithme est itérative. On part d'un paramètre initial, que l'on suppose « assez proche » d'un minimum, et qui constituera le vecteur $\beta$ de départ. Dans beaucoup de cas, un paramètre de départ « standard », tel que $\beta^T$=(1,1,…,1) fonctionnera sans problème. Dans certains cas, il n'y a convergence que si le vecteur de départ n'est pas trop éloigné de la solution.

À chaque itération, on remplace p par une nouvelle estimation $\beta + \delta$. Afin de déterminer q, les fonctions $ fi(\beta + \delta) $sont approchées en étant linéarisées :

    f$(\beta+\delta) ≈ f(\beta) + J \delta $

où on a noté J la jacobienne de f en $\beta$.

À un minimum de la somme des carrés S, on a $\nabla_{q} S=0$. En dérivant le carré de l'expression de droite, qui s'annule donc, on obtient :

    $(J^T J)\beta = J^T [y − f(\beta)] $

d'où l'on tire aisément $\delta$ en inversant $J^T J$.
Dans l'inversion matricielle, tout va dépendre du rapport entre la valeur propre la plus "grande" en norme, et la valeur propre la plus petite ; ce rapport, appelé "conditionnement de matrice", va concrètement refléter la robustesse de l'algorithme face au bruit.
Le point essentiel de l'algorithme de Levenberg-Marquardt est d'approcher cette équation, en l' « amortissant » un peu.
(On parle alors de "chargement de la diagonale" afin de contourner le mauvais conditionnement le cas échéant, problème que l'on retrouve avec l'algorithme de Capon et que l'on peut résoudre par une décomposition en valeurs singulières) :

   $ (J^T J + λ.diag(J^T J))q = J^T[y − f(\beta)].$

Le facteur d'amortissement positif λ est ajusté à chaque nouvelle itération. Si la diminution de S est rapide, on peut utiliser une valeur plus faible - ce qui rapproche l'algorithme de celui de Gauss-Newton. Si en revanche une itération est peu efficace, on peut augmenter λ - ce qui rapproche cette fois l'algorithme de celui du gradient. Un tel facteur d'amortissement est utilisé par exemple dans la régularisation de Tikhonov, utilisée pour résoudre certains problèmes linéaires.

Si on a effectué plus d'un certain nombre d'itérations, ou bien que l'on s'est approché suffisamment d'un minimum, la procédure se termine et renvoie le paramètre p comme estimation de la solution.

Le code est modélisé par le schéma suivant :​ 

 


 

Partie transporté

dans cette partie , on s'interesse à la simulation de la partie transportée, les équations explicitées dans la partie précédente sont résolues grâce à un algorithme de Merson d'ordre 4 codé en Fortran 90, ici le schéma utilisé est un celui de Runge Kutta, ce schéma repose sur le principe d'itération, c'est-à-dire qu'une première estimation de la solution est utilisée pour calculer une seconde estimation, plus précise, et ainsi de suite, ici cette procédure est répétée quatre fois, voici un schéma qui explique en détail la procédure suivie pour la résolution ainsi que les communication entre les différents bouts du code.