WikiPrépaLivrets

CCINP Modélisation PC 2018Sujet, corrigé et rapport du jury

Pas encore noté
  • Bilan enthalpique et thermique
  • Cinétique chimique
  • Équations différentielles ordinaires couplées
  • Régression linéaire matricielle
  • Méthodes numériques (Euler implicite, Newton-Raphson)

Téléchargements

Présentation du sujet

Difficulté moyenne
Modélisation de systèmes physique ou chimique PC : stabilité thermique d'un réacteur chimique exothermique
Afficher ou masquer la section

Le sujet étudie la stabilité thermique d'un réacteur parfaitement agité à double enveloppe siège d'une réaction exothermique. Une première partie établit le modèle thermique et le système d'équations différentielles couplées, une seconde partie traite numériquement des données expérimentales par régression linéaire matricielle puis intègre le système par la méthode d'Euler implicite couplée à Newton-Raphson.

  1. 1Partie 1 : modélisation du réacteurBilan enthalpique du réacteur, établissement du système d'équations différentielles couplées, étude de la stabilité thermique.
  2. 2Partie 2 : traitement numérique des donnéesRégression linéaire matricielle pour déterminer les paramètres du modèle, puis intégration numérique du système par la méthode d'Euler implicite et la méthode de Newton-Raphson.

Difficulté moyenne. Le rapport qualifie explicitement le sujet de difficulté moyenne, avec un niveau de difficulté des questions varié permettant de classer les candidats.

L'épreuve en chiffres

Moyenne 10,55 / 20 · 4 616 présents
Afficher ou masquer la section
Moyenne
10,55/ 20
Présents
4 616

Source : document officiel du concours. Notes publiées par le concours (après harmonisation le cas échéant). Courbe : estimation par une loi normale.

Ce qu'a observé le jury

5 erreurs relevées
Confusion température, énergie et puissance · Unité de l'énergie d'activation · Confusion entre adiabatique et isotherme
Afficher ou masquer la section

Le sujet fait appel à des notions transversales de chimie, physique, mathématiques et informatique, avec des parties indépendantes pour ne pas bloquer les candidats en difficulté sur l'une d'elles. La partie modélisation a été mieux traitée que la partie informatique, notamment l'écriture effective de code.

Les erreurs les plus sanctionnées

  1. 1
    Confusion température, énergie et puissanceQ1

    Des confusions entre température, énergie et puissance ont été relevées, avec un sens physique des termes mal compris et des analyses dimensionnelles peu rigoureuses.

  2. 2
    Unité de l'énergie d'activationQ9

    L'unité de l'énergie d'activation est trop souvent le joule alors qu'on attendait le J/mol.

  3. 3
    Confusion entre adiabatique et isothermeQ18

    Pour certains candidats le système adiabatique se traduisait par une différence de température nulle, pour d'autres par une dérivée temporelle nulle.

  4. 4
    Syntaxe numpy et confusion entier/réelQ21-Q26, Q40

    La syntaxe des fonctions de la bibliothèque numpy est rarement correcte, et il y a confusion entre entier et réel, la double division ne fonctionnant pas car elle renvoyait un réel.

  5. 5
    Écriture des développements limités peu rigoureuseQ27-Q34

    Il manque parfois le terme de reste o(delta t) ou le signe égal est utilisé à la place du signe de développement asymptotique.

Ce qui a été bien réussi

  • La première partie a été dans l'ensemble souvent traitée dans son intégralité.
  • Les questions 21 à 26 et 27 à 34 ont plutôt été bien traitées dans l'ensemble.
  • Quand la température finale est calculée correctement, l'analyse de la stabilité du réacteur est en général spontanée et correcte.

Conseils du jury

  • Traduire effectivement le modèle établi en code Python plutôt que de s'arrêter à l'écriture mathématique.
  • Préférer écrire un bout de code supplémentaire plutôt que d'utiliser une fonction intrinsèque dont on ne maîtrise pas la syntaxe.
  • Utiliser des noms de variables réalistes en code, par exemple B_chapeau plutôt qu'un symbole avec accent.

Synthèse rédigée par WikiPrépa à partir du rapport officiel du jury (à télécharger en PDF). Les citations sont extraites du rapport.

Ces sujets peuvent vous intéresser

Lecture du sujet en ligne

L'énoncé complet, avec les formules et les figures, sans ouvrir le PDF.
Afficher ou masquer la section

ÉPREUVE SPÉCIFIQUE - FILIÈRE PC

MODÉLISATION DE SYSTÈMES PHYSIQUES OU CHIMIQUES

Jeudi 3 mai: 8 h-12 h

N.B. : le candidat attachera la plus grande importance à la clarté, à la précision et à la concision de la rédaction. Si un candidat est amené à repérer ce qui peut lui sembler être une erreur d'énoncé, il le signalera sur sa copie et devra poursuivre sa composition en expliquant les raisons des initiatives qu'il a été amené à prendre.
Les calculatrices sont autorisées
Le sujet est composé de deux parties (pages 1 à 14)
et d'une annexe (pages 15 à 18).

PROBLÈME

Étude d'une réaction exothermique : stabilité thermique en réacteur fermé

Présentation du problème

De nombreux procédés industriels font intervenir des réacteurs fermés pour la synthèse de molécules à haute valeur ajoutée. Pour optimiser le rendement de la synthèse, il est nécessaire de bien comprendre l'influence des paramètres physiques (comme la température de la réaction,...) sur la marche du réacteur. La maîtrise des échanges thermiques est cruciale dans le cas des réactions exothermiques car la chaleur dégagée par la réaction provoque une augmentation de la température du mélange réactionnel. Selon les conditions opératoires, cette augmentation de température peut entraîner un emballement thermique du réacteur et conduire à des dégâts irréversibles.
L'accident de Seveso le 10 juillet 1976 illustre les problèmes liés à l'emballement thermique des réacteurs. Il s'agissait d'un procédé de production de 2,4,5-trichloro-phénol à partir de 1,2,4,5-tétrachloro-benzène et de soude dans un solvant (l'éthylène glycol) à une température voisine de 150^∘C et à pression atmosphérique en réacteur fermé. La mauvaise maîtrise de la température de la réaction a entraîné le déroulement de réactions secondaires conduisant d'une part à une augmentation de la température et de la pression et d'autre part à la formation de produits secondaires toxiques : les dioxines. La rupture de la soupape de sécurité due à l'augmentation de la pression a conduit au rejet de dioxines dans l'atmosphère.
L'étude de l'influence des paramètres physiques sur la marche d'un réacteur se fait la plupart du temps à l'échelle du laboratoire dans des dispositifs de dimensions beaucoup plus petites que celles des réacteurs utilisés dans les procédés industriels. La particularité du réacteur utilisé pour la présente étude est qu'il possède une géométrie cylindrique et qu'il est équipé d'une double enveloppe externe dans laquelle circule un fluide permettant de refroidir la paroi du réacteur et d'empêcher un emballement thermique (figure 1). L'emballement thermique survient lorsque la chaleur dégagée par la réaction excède la capacité du réacteur à évacuer l'énergie.
Figure 1 - Schéma simplifié d'un réacteur fermé parfaitement agité avec une double enveloppe pour son refroidissement
Pour caractériser le comportement thermique du réacteur, on commence la plupart du temps par une étude en l'absence de réaction. Cette étude permet dans un premier temps de caractériser la capacité du réacteur à évacuer l'énergie avec la détermination du coefficient de transfert thermique à la paroi. Dans un deuxième temps, on met en œuvre dans ce réacteur une réaction exothermique R → produits. Ces études permettent de déterminer les valeurs de paramètres clefs intervenant dans les équations décrivant le comportement du réacteur (modèle théorique). L'utilisation de ce modèle théorique permet de prédire la stabilité thermique du réacteur. L'établissement du modèle théorique repose sur l'écriture de bilans de matière et de chaleur. Une fois que l'influence des conditions physiques sur la marche du réacteur est déterminée, on peut en déduire les conditions de stabilité du réacteur industriel.
Ce sujet est constitué de deux parties. La première partie porte sur la modélisation du réacteur fermé parfaitement agité avec double enveloppe. Elle permet l'établissement du système d'équations différentielles régissant les variations de la conversion du réactif et de la température en fonction du temps, ainsi que la détermination de la valeur de paramètres physico-chimiques intervenant dans ces équations. La deuxième partie porte sur le traitement numérique des données expérimentales avec la détermination des paramètres d'un modèle par régression linéaire, puis la prédiction du comportement thermique du réacteur par résolution d'un système d'équations différentielles par la méthode d'Euler implicite.

Partie I - Modélisation du réacteur fermé parfaitement agité avec double enveloppe

I. 1 - Etalonnage thermique du réacteur

Dans cette partie, on souhaite caractériser les transferts de chaleur entre un liquide contenu à l'intérieur du réacteur et la paroi en l'absence de toute réaction chimique. On supposera que la capacité thermique massique et la masse volumique de ce liquide sont constantes quelle que soit la température. Pour caractériser ces transferts de chaleur, on réalise deux expériences successives avec la température de la paroi, T_p, maintenue constante dans les deux cas.
  • La première expérience consiste à chauffer le liquide (initialement à une température identique à celle de la paroi) avec un dispositif annexe (une résistance chauffante) dissipant une puissance thermique P_(th) de 96, 0 W. On constate que la température de la phase liquide augmente, puis atteint une valeur asymptotique en régime permanent (figure 2a).
  • Après avoir atteint le régime permanent lors de la phase de chauffe, on réalise une seconde expérience en coupant le chauffage. La température du liquide décroit jusqu'à ce qu'elle tende vers la température de la paroi (figure 2b).
Figure 2 - Évolution de la température du liquide dans le réacteur
On exploite d'abord la courbe obtenue lors de la phase de chauffe (figure 2a) pour déterminer le coefficient de transfert de chaleur à la paroi, noté U (unité : W ⋅ m^(− 2) ⋅ K^(− 1) ). Ce coefficient rend compte des échanges de chaleur entre la phase réactionnelle et le fluide caloporteur dans la double enveloppe à travers la paroi du réacteur. Dans le cas présent, la température T_p étant la température de paroi du côté du fluide caloporteur, il s'agit d'un coefficient de transfert thermique global qui tient compte du transfert convectif côté réaction et du transfert par conduction dans la paroi qui sépare les deux fluides.
Pour obtenir la valeur de U, on commence par établir l'équation différentielle qui régit l'évolution de la température T en fonction du temps en réalisant un bilan d'énergie.
Le bilan d'énergie, conséquence directe du premier principe de la thermodynamique, appliqué au système constitué par la phase réactionnelle lors de la phase de chauffe, conduit à l'équation différentielle suivante (équation (1))
(ρ × V × Cp)(dT)/(dt) = P_(th) − U × A × (T − T_p),
où T est la température du fluide à l'intérieur du réacteur, ρ, V et C_p sont respectivement la masse volumique, le volume et la capacité thermique massique du fluide, P_(th) est la puissance thermique cédée par la résistance chauffante au milieu réactionnel, T_p est la température à la paroi, maintenue constante ( T_p = 320, 0 K ) et A la surface latérale du réacteur sur laquelle le fluide à l'intérieur du réacteur est en contact avec la double enveloppe.
Q1. Interpréter concrètement chacun des trois termes du bilan d'énergie en précisant leur signification physique et vérifier qu'ils sont homogènes à des puissances.
Q2. Donner l'expression de T − T_p en régime permanent. Il est rappelé que la température de la paroi, T_p, est maintenue constante tout au long des essais. Il est précisé que la température de la phase liquide à l'instant initial est égale à T_p.
Q3. D'après les résultats obtenus lors de la première expérience (figure 2a), donner la valeur de la différence de température T − T_p lorsqu'on atteint le régime permanent. Calculer la valeur du coefficient de transfert de chaleur à la paroi dans les unités SI. On donne T_p = 320, 0 K, ρ = 1000, 0 kg ⋅ m^(− 3), V = 0, 1 × 10^(− 3) m^3 et C_p = 1800, 0 J ⋅ kg^(− 1) ⋅ K^(− 1), P_(th) = 96, 0 W et A = 8, 0 × 10^(− 3) m^2.
Q4. On souhaite faire apparaître un temps caractéristique d'échange thermique τ_c du système. Montrer que le bilan d'énergie peut se mettre sous la forme suivante (équation (2)) :
(dT)/(dt) = s + (T_p − T)/(τ_c).
Donner les expressions de s et τ_c. Vérifier que τ_c est homogène à un temps.
On souhaite maintenant exploiter les résultats obtenus lors de la phase de refroidissement (figure 2b) pour confirmer la valeur du temps caractéristique d'échange thermique déterminé précédemment.
Q5. Le bilan d'énergie établi à la question Q4 est-il modifié ? Si oui, donner la nouvelle expression de (dT)/(dt).
Q6. Donner l'expression de T − T_p en fonction du temps t par résolution de l'équation différentielle. On notera T_(max) la température initiale lors de la phase de refroidissement.
Q7. Le tracé de ln(T − T_p) (avec T − T_p en K ) en fonction du temps t (figure 3) donne une droite d'équation y = − 1, 33 × 10^(− 2) × x + 3, 68 (avec x en secondes). En déduire la valeur du temps caractéristique d'échange thermique τ_c. Calculer la valeur du coefficient de transfert de chaleur à la paroi et vérifier que cette valeur correspond à celle obtenue avec la première expérience.
Figure 3 - Tracé de ln(T − T_p) en fonction de t à l'aide des points enregistrés lors de la phase de refroidissement (figure 2b)

I. 2 - Etude d'une réaction exothermique en réacteur fermé à double enveloppe

Dans cette sous-partie, on considère que l'on met en œuvre une réaction chimique exothermique R → produits ( R est dissous dans un solvant) dans le même réacteur que celui dont on a étudié les échanges thermiques dans la sous-partie précédente. Il s'agit ici de caractériser le comportement thermique du réacteur en présence d'une réaction exothermique.
Le comportement du réacteur fermé parfaitement agité avec double enveloppe peut être représenté par un système constitué de deux équations différentielles ordinaires. La première équation différentielle ordinaire représente l'évolution temporelle de la conversion du réactif R; la deuxième permet de caractériser l'évolution de la température de la réaction en fonction du temps.
Le réactif R étant dissout dans un solvant, on considère que le volume du mélange réactionnel V reste constant au cours du temps. On considère également que la réaction est homogène et qu'elle a lieu dans tout le volume du mélange réactionnel.
On commence par déterminer l'équation différentielle qui régit l'évolution de la conversion du réactif R en fonction du temps.
Q8. On considère que la réaction est d'ordre 1 par rapport au réactif R. Donner l'expression de la vitesse de la réaction (exprimée par unité de volume de mélange réactionnel) que l'on notera r (on notera C_R la concentration molaire du réactif R et k la constante cinétique). Préciser la dimension et l'unité de r dans le Système International.
Q9. Rappeler la loi d'Arrhenius donnant les variations de la constante de réaction en fonction de la température. On notera k_0 le facteur pré-exponentiel et E_a l'énergie d'activation. Préciser les dimensions et les unités SI de k, k_0 et E_a.
Q10. Écrire le bilan de matière sous la forme (dC_R)/(dt) = f(C_R, T). Préciser l'expression de f(C_R, T). Il est rappelé que le volume de la phase réactionnelle reste constant au cours du temps.
Q11. Dans le cas où le réacteur fonctionnerait en marche isotherme, résoudre l'équation différentielle et donner l'expression de la concentration de R en fonction du temps sous la forme C_R = g(t). On notera C_R^0 la concentration initiale en R.
Q12. Pour simplifier les bilans, on introduit le taux de conversion de R, noté X_R, défini par la relation suivante: X_R = (C_R^0 − C_R)/C_R^0. Exprimer l'évolution de taux de conversion X_R en fonction du temps pour le cas de la marche isotherme.
Q13. Donner l'expression de l'équation différentielle qui régit l'évolution du taux de conversion X_R en fonction du temps dans le cas général (marche quelconque, c'est-à-dire non isotherme), sans chercher à la résoudre.
L'évolution de la température en fonction du temps est régie par une seconde équation différentielle obtenue en réalisant un bilan d'énergie sur le réacteur, conséquence directe du premier principe de la thermodynamique.
La réaction chimique qui se déroule dans le réacteur produit par unité de temps une variation de l'enthalpie du système réactionnel donnée par S(t, X, T) (équation (3)) qui correspond à la chaleur dégagée par la réaction. Ce paramètre est appelé «terme source» dans la suite.
S(t, X_R, T) = − Δ_r H^0(T) × r(t, X_R, T) × V,
où V est le volume du mélange réactionnel, r est la vitesse de la réaction et Δ_r H^0(T) est l'enthalpie molaire standard de la réaction. Dans la suite, on suppose que Δ_r H^0(T) ne dépend pas de la température. On prendra Δ_r H^0(T) = Δ_r H^0(T_p) que l'on notera Δ_r H^0 pour simplifier.
Q14. Donner la dimension du terme source S(t, X_R, T).
Q15. Montrer qu'un bilan enthalpique appliqué à un système que l'on précisera avec soin permet d'aboutir à la relation suivante (équation (4))
(dT)/(dt) = J(dX_R)/(dt) − (T − T_p)/(τ_c)
où l'on exprimera J et τ_c en fonction de Δ_r H, C_R^0, ρ, C_p, V, U et A. On admettra qu'il est légitime de négliger la contribution des réactifs, des produits et des accessoires situés à l'intérieur du réacteur au travers de leur capacités thermiques devant celle du solvant.
Q16. Donner l'expression du paramètre J ainsi que sa dimension.
Q17. Calculer la valeur du paramètre J en unité SI . On donne Δ_r H^0 = − 360, 0 kJ ⋅ mol^(− 1), ρ = 1000, 0 kg ⋅ m^(− 3), C_p = 1800, 0 J ⋅ kg^(− 1) ⋅ K^(− 1) et C_R^0 = 500, 0 mol ⋅ m^(− 3).

I. 3 - Stabilité thermique du réacteur

Une première condition de stabilité, valable pour le cas d'une marche adiabatique, est que la température finale T_f de la réaction soit inférieure à une température maximale T_(max), telle que T_(max) = 1, 25 × T_p.
Q18. Exprimer l'équation différentielle (équation (4)) dans le cas d'un fonctionnement adiabatique.
Q19. Déterminer alors l'expression de la température T en fonction du taux de conversion X_R, du paramètre J et de T_0 la température initiale du mélange réactionnel.
Q20. Donner l'expression de la température T_f atteinte en fin de réaction dans le cas d'une marche adiabatique sachant que T_0 = T_p = 320, 0 K. Conclure quant à la stabilité du réacteur dans le cas de cette étude. Donner la signification physique du paramètre J.

Partie II - Traitement numérique des données expérimentales

Les algorithmes demandés au candidat devront être réalisés dans le langage Python. On supposera les bibliothèques «numpy » et «matplotlib.pyplot » chargées. Une annexe présentant les fonctions usuelles de Python est disponible pages 15 à 18. Les commentaires suffisants à la compréhension du programme devront être apportés et des noms de variables explicites devront être utilisés lorsque ceux-ci ne sont pas imposés.

II. 1 - Détermination des paramètres d'un modèle par régression linéaire

Pour calculer la valeur du temps caractéristique d'échange thermique du réacteur à la question Q7, un modèle de régression linéaire simple a été estimé à partir des points expérimentaux enregistrés lors de la phase de refroidissement.
Le modèle de régression linéaire simple (fonction affine) est un modèle de régression d'une variable expliquée ( ln(T − T_p) dans notre cas) sur une variable explicative (le temps t dans notre cas) dans lequel on fait l'hypothèse que la fonction qui relie la variable explicative à la variable expliquée est linéaire dans ses paramètres.
Soit n le nombre de points expérimentaux. Le modèle linéaire simple s'écrit de la manière suivante pour un point i(1 ≤ i ≤ n)
y_i = β_1 ˆ × x_i + β_0 ˆ
où β_0 ˆ et β_1 ˆ sont les paramètres du modèle, y_i est la variable expliquée et x_i est la variable explicative.
On propose de déterminer les paramètres du modèle par deux méthodes directes.
La méthode consiste à écrire le modèle (équation (5)) sous la forme matricielle Y = L × Bˆ.Bˆ est un vecteur colonne contenant les paramètres du modèle β_0 ˆ et β_1 ˆ, Y est un vecteur colonne contenant les n valeurs y_i et L une matrice à n lignes et 2 colonnes, telle que L(i, 1) = (∂y_i)/(∂β_0 ˆ) et L(i, 2) = (∂y_i)/(∂β_1 ˆ). Rappelons que Bˆ = (L^t × L)^(− 1) × L^t Y où L^t est la matrice transposée de L.
Q21. Donner les expressions de (∂y_i)/(∂β_0 ˆ) et (∂y_i)/(∂β_1 ˆ). En déduire la valeur des coefficients de la matrice L.
Q22. On suppose que les vecteurs colonnes Y et X, qui contiennent les valeurs y_i et x_i(1 ≤ i ≤ n), sont déjà créés. Donner le code permettant de créer la matrice L.
Q23. Donner le code permettant de déterminer les coefficients de la matrice P = L^t × L. Préciser les dimensions de la matrice P.
Q24. Donner le code permettant de déterminer les coefficients de la matrice Q = L^t × Y. Préciser les dimensions de la matrice Q.
Q25. Donner le code permettant de créer une fonction inv_mat(M) qui renvoie la matrice inverse de la matrice M de dimension ( 2 × 2 ) donnée comme argument d'entrée.
Q26. On note N la matrice inverse de M. Donner le code permettant de déterminer les coefficients β_0 ˆ et β_1 ˆ de la matrice Bˆ.

II. 2 - Prédiction du comportement thermique du réacteur

Dans cette sous-partie, on souhaite utiliser le modèle constitué du système d'équations différentielles établies dans la Partie I qui décrit l'évolution du taux de conversion du réactif X_R et de la température de réaction T en fonction du temps pour prédire le comportement thermique du réacteur en présence d'une réaction exothermique.
Pour simplifier les notations, on met le système d'équations différentielles sous la forme suivante (équation (6)) :
{(dX_R)/(dt) = f_1(t, X_R, T); (dT)/(dt) = f_2(t, X_R, T)
La méthode de résolution proposée pour résoudre le système d'équations différentielles est la méthode d'Euler implicite à pas fixe. Cette méthode est préférée car elle donne de meilleurs résultats que la méthode explicite pour les systèmes dits raides (un système raide est un système qui est caractérisé par une évolution rapide des phénomènes en fonction du temps, ce qui est le cas ici pour la température de réaction).
Soit une variable y qui dépend du temps t. Comme la méthode d'Euler explicite, la méthode d'Euler implicite consiste à évaluer la valeur de y(t + Δt) à partir de celle de y(t) et de la dérivée (∂y)/(∂t) (figure 4, page suivante). La différence entre les deux méthodes réside dans le choix de l'abscisse à laquelle est évaluée la dérivée (∂y)/(∂t). Dans le cas de la méthode explicite, elle est évaluée en t : (∂y)/(∂t)(t) comme le montre le schéma de la figure 4a, page suivante. Pour la méthode implicite, elle est évaluée en t + Δt : (∂y)/(∂t)(t + Δt) (figure 4b, page suivante).
Figure 4 - Approximation de y(t + Δt) par la méthode d'Euler
Dans le cas d'un schéma implicite (figure 4b), l'expression de y(t + Δt) en fonction de y(t) et de la dérivée (∂y)/(∂t)(t + Δt) évaluée en t + Δt est obtenue en réalisant un développement limité dit rétrograde : y(t) = y(t + Δt) − Δt × (∂y)/(∂t)(t + Δt) + o(Δt).
Q27. À l'aide d'un développement limité rétrograde de la fonction X_R(t), donner l'expression de X_R(t + Δt) à l'ordre 1 en fonction de X_R(t) et de sa dérivée partielle par rapport à t, (dX_R)/(dt)(t + Δt) évaluée en t + Δt.
Q28. En déduire une valeur approchée de (dX_R)/(dt)(t + Δt) à l'ordre 0 en fonction de X_R(t), X_R(t + Δt) et Δt.
Q29. À l'aide d'un développement limité rétrograde de la fonction T(t), donner l'expression de T(t + Δt) à l'ordre 1 en fonction de T(t) et de sa dérivée partielle par rapport à t, (dT)/(dt)(t + Δt).
Q30. En déduire une valeur approchée de (dT)/(dt)(t + Δt) à l'ordre 0 en fonction de T(t), T(t + Δt) et Δt.
On procède à la discrétisation des équations. On note X_(Ri) la conversion évaluée au temps t_i, X_(Ri + 1) la conversion évaluée au temps t_(i + 1) et (dX_R)/(dt)|_(t_(i + 1)) la dérivée de X_R évaluée à l'instant t_(i + 1). De même, on note T_i la température évaluée au temps t_i, T_(i + 1) la température évaluée au temps t_(i + 1) et (dT)/(dt)|_(t_(i + 1)) la dérivée de T évaluée à l'instant t_(i + 1).
Q31. Donner l'expression de (dX_R)/(dt)|_(t_(i + 1)) en fonction de X_(Ri), X_(Ri + 1) et Δt.
Q32. Donner l'expression approchée de X_(Ri + 1) en fonction de X_(Ri), Δt et de la fonction f_1(t_(i + 1), X_(Ri + 1), T_(i + 1)), évaluée en t_(i + 1).
Q33. Donner l'expression de (dT)/(dt)|_(t_(i + 1)) en fonction de T_i, T_(i + 1) et Δt.
Q34. Donner l'expression approchée de T_(i + 1) en fonction de T_i, Δt et de la fonction f_2(t_(i + 1), X_(Ri + 1), T_(i + 1)), évaluée en t_(i + 1).
On constate que les expressions obtenues aux questions Q32 et Q34 constituent un système non linéaire dont les inconnues sont X_(Ri + 1) et T_(i + 1). On propose d'utiliser la méthode de Newton-Raphson pour trouver les valeurs de X_(Ri + 1) et T_(i + 1) à chaque itération de la méthode d'Euler.
La méthode de Newton-Raphson pour la résolution d'un système de n équations non linéaires à n inconnues x = (x_1, …, x_n), mis sous la forme de l'équation (7) suivante,
g(x) = [g_1(x_1, …, x_n); …; g_n(x_1, …, x_n)] = 0
est une extension de la méthode de Newton permettant de trouver la racine d'une fonction d'une variable.
On peut démontrer la formule de récurrence suivante (équation (8)) :
x^(j + 1) = x^j − (Dg(x^j))^(− 1)g(x^j),
où x^(j + 1) est la valeur du vecteur x à l'itération j + 1, x^j est la valeur du vecteur x à l'itération j, g(x^j) est la valeur de g(x) à l'itération j et Dg(x^j) est la matrice Jacobienne évaluée en x^j (équation (9)) :
Dg(x^j) = [(∂g_1)/(∂x_1), …, (∂g_1)/(∂x_n); …, …; (∂g_2)/(∂x_1), …, (∂g_2)/(∂x_n)]_(x = x^j).
La formule de récurrence s'accompagne du choix d'une valeur initiale, notée x^0, et d'un critère d'arrêt, par exemple ‖x^(j + 1) − x^j‖ ≤ ε.
Q35. Transformer les expressions obtenues aux questions Q32 et Q34 pour les mettre sous la forme g_1(X_(Ri + 1), T_(i + 1)) = 0 et g_2(X_(Ri + 1), T_(i + 1)) = 0.
Q36. Donner les expressions de (∂g_1)/(∂X_(Ri + 1)), (∂g_1)/(∂T_(i + 1)), (∂g_2)/(∂X_(Ri + 1)) et (∂g_2)/(∂T_(i + 1)) permettant de construire la matrice Jacobienne Dg(X_(Ri + 1), T_(i + 1)).
Q37. Écrire une fonction mat_Dg(x) qui a pour argument d'entrée un vecteur x contenant les valeurs de X_(Ri + 1) et T_(i + 1) à l'itération j et qui retourne la matrice Jacobienne Dg(x^j). On supposera que les paramètres suivants ont été au préalable déclarés comme variables globales : Δt = 0, 01 s, k_0 = 5, 0 s^(− 1), E_a = 20000, 0 J.mol^(− 1), J = 100, 0 K, τ_c = 75 s et R = 8, 314 J.K^(− 1).mol^(− 1).
Q38. Écrire le code permettant de calculer les valeurs du vecteur x à l'itération j + 1 à l'aide de l'équation (8) lors d'une boucle de l'algorithme de Newton-Raphson. On notera x_old la valeur de x a l'itération j et x _new la valeur de x à l'itération j + 1. De même, on notera x_euleri et T_euleri les vecteurs contenant les valeurs de X_(Ri) et T_i à l'itération i de la méthode d'Euler implicite. Pour l'inversion de matrice, on utilisera la fonction inv_mat(M) écrite à la question Q25.
Q39. Pour obtenir la valeur de x_new par la méthode de Newton-Raphson, on souhaite créer une boucle itérative avec condition. La condition d'arrêt porte sur la valeur absolue de la différence des températures T_(i + 1)^j et T_(i + 1)^(j + 1) que l'on souhaite inférieure à 10^(− 5) K. Pour les valeurs initiales de X_(ri + 1)^0 et T_(i + 1)^0 on prendra respectivement 0,5 et T_p + J/2. Écrire le code correspondant. Il est inutile de recopier l'intégralité du code écrit à la question précédente ; on indiquera néanmoins sa place dans le code de cette question.
Maintenant que le code permettant de trouver les valeurs de X_(Ri + 1) et T_(i + 1) par la méthode de Newton-Raphson lors d'une itération de la méthode d'Euler implicite a été établi, on souhaite calculer les valeurs pour l'ensemble des itérations de la méthode d'Euler implicite. On rappelle que X_R(t = 0) = 0 et T(t = 0) = T_p = 320, 0 K. En plus de noter x_−euleri et T_euleri les vecteurs contenant les valeurs de X_(Ri) et T_i pour chaque itération i de la méthode d'Euler implicite, on notera t_euleri le vecteur contenant les valeurs de t_i à chaque itération. L'intégration sera réalisée sur l'intervalle [t_0, t_f] avec t_0 = 0 s et t_f = 1000 s. Le pas de temps Δt est de 0, 01 s.
Q40. Donner le code permettant de calculer le nombre m d'intervalles Δt compris dans l'intervalle [t_0, t_f] ( m est un entier).
Q41. Écrire le code permettant de calculer les valeurs des éléments des vecteurs x_euleri, т_euleriet t_euleri à chaque itération de la méthode d'Euler implicite.
Q42. Donner le code permettant de tracer les graphes de la figure 5 montrant l'évolution de la conversion et de la température en fonction du temps que l'on obtiendrait en réalisant la résolution numérique du système d'équations différentielles (simulation réalisée avec k_0 = 5, 0 s^(− 1), E_a = 20000, 0 J ⋅ mol^(− 1), J = 100, 0 K, τ_c = 75 s).
Figure 5 - Température et conversion calculées à partir du modèle constitué du système d'équations différentielles déterminées dans la Partie I

FIN

Annexe : Fonctions de Python

1. Bibliothèque numpy de Python

Dans les exemples ci-dessous, la bibliothèque numpy a préalablement été importée à l'aide de la commande : import numpy as np
On peut alors utiliser les fonctions de la bibliothèque, dont voici quelques exemples :

np.array(liste)

Description : fonction permettant de créer une matrice (de type tableau) à partir d'une liste. Argument d'entrée : une liste définissant un tableau à 1 dimension (vecteur) ou 2 dimensions (matrice).
Argument de sortie : un tableau (matrice).
Exemples: np.array([4,3,2])
⇒ [4, 3, 2]
np.array([[5],[7],[1]])
=> [[5]
    [7]
    [1]]
np.array([[3,4,10],[1,8,7]])
=> [[3 4 10]
    [1 8 7]]
A[i, j].
Description : fonction qui retourne l'élément ( i + 1, j + 1 ) de la matrice A. Pour accéder à l'intégralité de la ligne i + 1 de la matrice A, on écrit A[i, : ]. De même, pour obtenir toute la colonne j + 1 de la matrice A, on utilise la syntaxe A[ :, j].
Arguments d'entrée : une liste contenant les coordonnées de l'élément dans le tableau A.
Argument de sortie : l'élément ( i + 1, j + 1 ) de la matrice A.
ATTENTION: en langage Python, les lignes d'un tableau A de dimension n x m sont
numérotées de 0 à n-1 et les colonnes sont numérotées de 0 à m-1
Exemple : A = np.array([[3, 4, 10], [1, 8, 7]])
A[0,2]
=> 10
A[1,:]
⇒ [1, 8, 7]
A[:,2]
d> [10 7]

np.zeros((n,m))

Description : fonction créant une matrice (tableau) de dimensions n × m dont tous les éléments sont nuls.
Arguments d'entrée : un tuple de deux entiers correspondant aux dimensions de la matrice à créer.
Argument de sortie : un tableau (matrice) d'éléments nuls.
Exemple: np.zeros((3,4))
⇒ [[0, 0, 0, 0]
[0 0 0 0 0 ]
[0 0 0 0 0 0 ]]

np.linspace(Min,Max,nbElements)

Description: fonction créant un vecteur (tableau) de nbElements nombres espacés régulièrement entre Min et Max. Le 1^(er) élément est égal à Min, le dernier est égal à Max et les éléments sont espacés de (Max-Min)/(nbElements - 1) :
Arguments d'entrée : un tuple de 3 entiers.
Argument de sortie : un tableau (vecteur).
Exemple: np.linspace(3,25,5)
⇒ [3, 8.5, 14, 19.5, 25]

np.loadtxt('nom_fichier',delimiter='string',usecols=[n])

Description : fonction permettant de lire les données sous forme de matrice dans un fichier texte et de les stocker sous forme de vecteurs.
Arguments d'entrée : le nom du fichier qui contient les données à charger, le type de caractère utilisé dans ce fichier pour séparer les données (par exemple un espace ou une virgule) et le numéro de la colonne à charger (ATTENTION, la première colonne porte le numéro 0 ). Argument de sortie : un tableau.
Exemple: data=np.loadtxt('fichier.txt',delimiter=' ',usecols=[0]) #dans cet exemple data est un vecteur qui correspond à la première #colonne de la matrice contenue dans le fichier 'fichier.txt'.

2. Bibliothèque matplotlib.pyplot de Python

Cette bibliothèque permet de tracer des graphiques. Dans les exemples ci-dessous, la bibliothèque matplotlib.pyplot a préalablement été importée à l'aide de la commande :
import matplotlib.pyplot as plt
On peut alors utiliser les fonctions de la bibliothèque, dont voici quelques exemples :

plt.plot(x,y)

Description : fonction permettant de tracer un graphique de n points dont les abscisses sont contenues dans le vecteur x et les ordonnées dans le vecteur y. Cette fonction doit être suivie de la fonction plt.show() pour que le graphique soit affiché.
Arguments d'entrée : un vecteur d'abscisses x (tableau de dimension n ) et un vecteur d'ordonnées y (tableau de dimension n ).
Argument de sortie : un graphique.
Exemple: x = np.linspace (3, 25, 5)
y = sin(x); plt.plot (x, y); plt.title('titre_g raphique'); plt.xlabel('x'); plt.ylabel('y'); plt.show()

plt.title('titre')

Description : fonction permettant d'afficher le titre d'un graphique.
Argument d'entrée : une chaîne de caractères.

plt.xlabel('nom')

Description : fonction permettant d'afficher le contenu de nom en abscisse d'un graphique.
Argument d'entrée : une chaîne de caractères.

plt.ylabel('nom')

Description : fonction permettant d'afficher le contenu de nom en ordonnée d'un graphique.
Argument d'entrée : une chaîne de caractères.

plt.show()

Description : fonction réalisant l'affichage d'un graphe préalablement créé par la commande plt.plot( x, y ). Elle doit être appelée après la fonction plt.plot et après les fonctions plt.xlabel et plt.ylabel.

3. Fonction intrinsèque de Python

sum(x)

Description : fonction permettant de faire la somme des éléments d'un vecteur ou tableau.
Arguments d'entrée : un vecteur ou un tableau de réel, entier ou complexe.
Argument de sortie : un scalaire y qui est la somme des éléments de x .
Exemple : y = sum(x)
//y retourne la somme des éléments de x .

Questions fréquentes

3 questions
Sur quels chapitres porte le sujet de modélisation PC 2018 CCINP ?
Afficher ou masquer la section

Sur quels chapitres porte le sujet de modélisation PC 2018 CCINP ?

Il porte sur le bilan thermique et la cinétique d'un réacteur chimique exothermique, puis sur la régression linéaire matricielle et les méthodes numériques d'intégration (Euler implicite, Newton-Raphson) en Python.

Quelles erreurs le jury a-t-il le plus relevées ?

Des confusions entre température, énergie et puissance, une unité d'énergie d'activation fausse, une confusion entre adiabatique et isotherme, et une syntaxe numpy souvent incorrecte.

Le sujet de modélisation PC 2018 est-il difficile ?

Le rapport le qualifie de difficulté moyenne, avec des questions de niveau varié qui permettent de classer les candidats.

Pas de description pour le moment