OpenSup

SESSION 2017 PCMS006 EPREUVE SPECIFIQUE - FILIERE PC ____________________ MODELISATION DE SYSTEMES PHYSIQUES OU CHIMIQUES Jeudi 4 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, largement indépendantes. 1/15 Autour de l’équation de Poisson Ce problème s’intéresse à la résolution numérique de quelques problèmes d’électrostatique. Il se compose de deux parties. I. Étude de l’équation de Poisson et de différentes méthodes de résolution numérique. II. Deux études de cas : fil infini chargé et mouvement d’un électron entre les plaques d’un conden- sateur. Les différentes parties sont largement indépendantes. Un aide-mémoire numpy/matplotlib/pyplot est présent à la fin du sujet. Partie I - Équation de Poisson I.1 - Établissement de l’équation Q1. Rappeler l’équation de Maxwell-Gauss ainsi que la relation entre le champ ⃗E et le potentiel électrostatique V. En déduire l’équation de Poisson : ∆V + ρ ε0 = 0 . Préciser les noms et les unités usuelles de ρ et ε0. Q2. Citer plusieurs situations physiques en dehors de l’électrostatique pour lesquelles il existe une équation analogue. I.2 - Équation adimensionnée pour un problème plan On veut résoudre l’équation de Poisson dans une portion de plan P carrée de côté L. On pose : X = x/L, Y = y/L . Q3. Montrer qu’on peut écrire l’équation sous la forme suivante : ∂2V(X, Y) ∂X2 + ∂2V(X, Y) ∂Y2 + ρ′(X, Y) = 0 où ρ′(X, Y) sera exprimé en fonction de ρ, L et ε0. I.3 - Discrétisation Afin de résoudre numériquement l’équation de Poisson, on va utiliser un maillage de P, de pas h = 1/N, et on va transformer les dérivées partielles par des différences entre les valeurs de V aux différents points du maillage (on parle aussi des nœuds du maillage). La figure 1 (page suivante) représente le maillage de P pour N = 5. 2/15 Méthode de Jacobi À partir de l’équation (1), on peut exprimer : V(i, j) = 1 4(V(i + 1, j) + V(i − 1, j) + V(i, j + 1) + V(i, j − 1) + ρ′′(i, j)) . (2) La résolution s’effectue alors en deux étapes. — Initialisation a) On fixe le potentiel des points de F à la valeur imposée physiquement (bords et élec- trodes). b) On donne aux points de potentiel inconnu, donc appartenant à P\F , une valeur arbitraire V0(i, j), en général nulle. — Itérations On calcule une nouvelle valeur V1(i, j) des potentiels en appliquant l’équation (2) pour tous les points de P\F , tandis que V1(i, j) = V0(i, j) pour les points de F . Le processus est répété jusqu’à obtenir des valeurs du potentiel quasiment stables. En notant k le nombre d’itérations, on a donc pour le point de coordonnées (i, j) n’appartenant pas à la frontière : Vk+1(i, j) = 1 4(Vk(i + 1, j) + Vk(i − 1, j) + Vk(i, j + 1) + Vk(i, j − 1) + ρ′′(i, j)) . (3) La convergence de la méthode est vérifiée à l’aide du critère de convergence ek, défini par : ek = 1 N2 i, j (Vk+1(i, j) − Vk(i, j))2 . (4) Le calcul sera stoppé au bout de k itérations, quand ek deviendra inférieur à un seuil de convergence ε fixé arbitrairement. Implémentation informatique On va utiliser la bibliothèque numpy permettant une utilisation simple des tableaux de flottants à deux dimensions ; un aide-mémoire est disponible en fin de sujet. Le chargement des bibliothèques classiques est assuré par les lignes suivantes : # importation des bibliothèques import numpy as np import matplotlib.pyplot as plt import math On supposera que les tableaux numpy suivants, utilisés comme arguments dans les fonctions à définir dans les questions qui suivent, ont pour signification : — V[i,j], (i, j) ∈ ⟦0 . . . N⟧2 : tableau courant du potentiel en un point de P, — rhos[i,j], (i, j) ∈ ⟦0 . . . N⟧2 : tableau contenant la densité de charge ρ′′ en un point de P, — frontiere[i,j], (i, j) ∈ ⟦0 . . . N⟧2 : tableau de booléens indiquant si le point de coordon- nées (i, j) appartient ou non à F . En particulier, tous les points du bord du domaine seront tels que frontiere[i,j]==True. 4/15 Q6. Écrire la fonction nouveau_potentiel(V, rhos, frontiere, i, j) retournant la nou- velle valeur du potentiel au point (i, j) ∈ ⟦0 . . . N⟧2 selon l’équation (3). Q7. Montrer que pour modifier toutes les valeurs contenues dans V[i,j] pendant une itération, il est nécessaire de disposer d’une copie de ce tableau. On rappelle que l’attribut shape permet de récupérer les dimensions d’un tableau numpy. Q8. Écrire la fonction itere_J(V, rhos, frontiere) modifiant la totalité du tableau V[i,j] lors d’une seule itération et retournant l’erreur calculée conformément à l’équation (4). Q9. Écrire la fonction poisson(f_iter, V, rhos, frontiere, eps) ayant pour premier ar- gument une fonction du même type que celle définie à la question précédente, pour dernier argument eps le seuil arbitraire de convergence ε et dont le rôle est de modifier le tableau des potentiels V[i,j] jusqu’à convergence. I.5 - Améliorations Méthode de Gauss-Seidel C’est une modification de la méthode de Jacobi, pour laquelle on montre que la convergence est légèrement plus rapide. Supposons que l’on balaye le tableau des potentiels selon les indices i et j croissants : dans ces conditions, les points situés à gauche et en dessous du point courant ont déjà été calculés. On va utiliser ces nouvelles valeurs, probablement plus proches de la solution, dans la formule permettant le calcul de Vk+1(i, j). Ceci donne l’algorithme de Gauss-Seidel : Vk+1(i, j) = 1 4(Vk(i + 1, j) + Vk+1(i − 1, j) + Vk(i, j + 1) + Vk+1(i, j − 1) + ρ′′(i, j)) . (5) Q10. Montrer qu’il n’est plus nécessaire de copier le tableau V[i,j] pour la mise à jour lors d’une itération en utilisant l’équation (5). Faut-il modifier la fonction nouveau_potentiel pour passer de la méthode de Jacobi à celle de Gauss-Seidel ? Q11. Écrire la fonction itere_GS(V, rhos, frontiere) modifiant la totalité du tableau V[i,j] lors d’une seule itération et retournant l’erreur calculée conformément à l’équation (4). Méthode de Gauss-Seidel adaptative Les méthodes de Jacobi et de Gauss-Seidel n’utilisent pas la valeur de Vk(i, j) pour calculer Vk+1(i, j). La méthode de sur-relaxation (Successive Over Relaxation method) consiste à calculer la nouvelle valeur d’un nœud comme une combinaison linéaire de la valeur courante et de celle donnée par le schéma de Gauss-Seidel. En introduisant le paramètre de relaxation ω, on a alors : Vk+1(i, j) = (1 − ω)Vk(i, j) + ω 4 (Vk(i + 1, j) + Vk+1(i − 1, j) + Vk(i, j + 1) + Vk+1(i, j − 1) + ρ′′(i, j)) . (6) L’étude mathématique de cette relation permet de montrer les résultats suivants : — la méthode converge uniquement si 0 < ω < 2 et elle converge plus rapidement que la méthode de Gauss-Seidel si 1 < ω < 2, — il existe une valeur optimale de ω qui permet la convergence avec un nombre d’itérations en O(N) pour une valeur de ε fixée. 5/15 Pour la résolution de l’équation de Poisson envisagée dans ce problème (conditions de Dirichlet sur un maillage carré), on montre que la valeur optimale ωopt est : ωopt = 2 1 + π/N . (7) Q12. Écrire la fonction nouveau_potentiel_SOR(V, rhos, frontiere, i, j, omega) re- tournant la nouvelle valeur du potentiel au point (i, j) selon l’équation (6). Q13. Écrire la fonction itere_SOR(V, rhos, frontiere) optimale modifiant la totalité du ta- bleau V[i,j] lors d’une seule itération et retournant l’erreur calculée conformément à l’équa- tion (4). La résolution du problème peut alors se faire par un appel de la forme poisson(itere_SOR, V, rhos, frontiere, eps), les tableaux carrés V, rhos, frontiere étant de dimensions convenables pour représenter un maillage comportant (N + 1)2 nœuds. Q14. Quelle est la complexité temporelle de l’appel précédent quand ω = ωopt ? La figure 2 représente, pour ε = 10−4, la durée d’exécution T (en secondes) en fonction de N. Cette courbe est-elle en accord avec la complexité temporelle attendue ? Quelle serait la durée d’exécution pour N = 1 000 ? Commenter. T (gauss_seidel_SOR) = f ( N ) Figure 2 – Durée d’exécution (en secondes) de poisson(itere_SOR, V, rhos, frontiere, eps) en fonction de N 6/15 T (gauss_seidel_SOR) = f ( N ) I.6 - Détermination du champ électrique Connaissant le potentiel V[i,j], il est souvent nécessaire de calculer numériquement les compo- santes Ex et Ey du champ électrique au niveau des nœuds du maillage, qui sont alors conservées dans deux tableaux Ex et Ey dimensionnés correctement (ce sont donc des tableaux carrés de (N + 1)2 éléments). Q15. Expliquer rapidement comment il serait possible de définir la fonction calc_ExEy(Ex, Ey, V, h), permettant, à partir du tableau V et du pas du maillage h, le remplissage des deux tableaux Ex et Ey passés en arguments. Remarque : on ne demande pas d’écrire le code de la fonction, juste de décrire précisément les étapes de calcul, ainsi que les différents cas à considérer. Partie II - Deux études de cas II.1 - Fil cylindrique chargé uniformément Étude théorique On considère dans le vide un fil cylindrique infini d’axe z et de rayon R, portant une charge volumique constante ρ. Q16. En se plaçant en coordonnées cylindriques d’axe z, montrer par des considérations de symétrie et d’invariance que le champ ⃗E = E(r) ⃗ur. En déduire la forme des surfaces équipotentielles. Q17. En appliquant le théorème de Gauss, calculer le champ ⃗E dans tout l’espace. Tracer rapidement l’allure de E(r) en fonction de r. Q18. On donne : ε0 = 8,85 × 10−12 F.m−1, ρ = 1,00 × 10−5 C.m−3, R = 5,00 cm. Calculer la valeur maximale de la norme du champ électrique, ainsi que la valeur pour r = 2R. Étude numérique Pour pouvoir utiliser la méthode de Gauss-Seidel adaptative, on place le fil infini au centre d’une enceinte de longueur infinie et de section carrée (L × L), portée au potentiel nul (figure 3). Dans la suite, on prendra L = 4R = 20,0 cm. V = 0 fil chargé 2 R enceinte L 2 L L 2 L y x O Figure 3 – Fil infini dans une enceinte de section carrée, portée au potentiel nul 7/15 Le programme permettant la résolution de ce problème commence ainsi : # initialisations eps0 = 8.85e-12 # epsilon_0 L = 20.0e-2 # 20 cm N = 100 ; h = L/N # définition du maillage rho = 1.00e-5 # densité vol. de charge rho # les tableaux globaux numpy pour le cylindre chargé rhos_cyl = np.zeros((N+1,N+1)) # tableau des valeurs de rho’’ V_cyl = np.zeros((N+1,N+1)) # le potentiel Ex_cyl = np.zeros((N+1,N+1)) # la composante Ex Ey_cyl = np.zeros((N+1,N+1)) # la composante Ey # le tableau définissant la frontière est initialement # rempli entièrement par la valeur False frontiere_cyl = np.zeros((N+1,N+1), bool) Q19. Écrire la fonction dans_cylindre(x,y,xc,yc,R) retournant un résultat booléen indiquant si le point de coordonnées (x, y) est à l’intérieur ou sur le bord du cercle de centre (xc, yc) et de rayon R. Q20. Écrire la fonction initialise_rhos_cylindre(tab_rhos), initialisant le tableau rhos_cyl contenant les valeurs ρ′′(i, j) pour les nœuds du maillage. Q21. Écrire la fonction initialise_frontiere_cylindre(tab_f), mettant à True les points appartenant à la frontière, donc de potentiel fixé. La résolution numérique avec la méthode de Gauss-Seidel adaptative, utilisant les valeurs numériques précédentes et un seuil de convergence ε = 10−5, mène à la figure 4, où on a tracé un réseau de courbes équipotentielles, le potentiel V et la composante Ex du champ le long de l’axe de symétrie défini par y = L/2. En outre, la valeur calculée de Ex_cyl[50, 50] est égale à -0.0023321214257521206. Q22. Commenter le plus complètement possible ces résultats ; on veillera, en particulier, à les com- parer au modèle théorique (allure des courbes, valeurs numériques ...). Figure 4 – Équipotentielles, V(x, y) et ⃗E(x, y) · ⃗ux en fonction de x pour y = L/2 Une autre résolution est effectuée, avec une répartition de charges dans le cylindre différente de la précédente, utilisant la même valeur de la densité volumique ρ = 10−5 C.m−3. Elle mène aux courbes de la figure 5 (page suivante). 8/15 Figure 5 – Équipotentielles, V(x, y) et ⃗E(x, y) · ⃗ux en fonction de x pour la nouvelle répartition de charge Q23. Déduire de ces courbes la répartition de charges dans le cylindre dans cette deuxième situation. Calculer le champ électrique en tout point pour cette répartition de charges dans le cas d’un cylindre infini seul dans l’espace. On vérifiera que la valeur maximale du champ électrique calculée à l’aide de cette modélisation est compatible avec celle déduite de la figure 5. II.2 - Mouvement d’un électron dans un tube d’oscilloscope La figure 6 montre un tube d’oscilloscope de petite dimension, dans lequel des électrons émis par la cathode sont accélérés et déviés vers un écran luminescent. La déviation est assurée par le passage des électrons entre les plaques de deux condensateurs plans : un pour la déviation horizontale, l’autre pour la déviation verticale. L’étude qui suit ne concernera que le condensateur responsable de la déviation verticale. cathode écran Figure 6 – Petit tube d’oscilloscope, de longueur d’environ 20 cm On modélise la trajectoire d’un électron de la façon suivante (figure 7 page suivante, où la zone de déviation est grisée) : — on négligera l’effet de la pesanteur, — émis à vitesse nulle par effet thermo-électronique au niveau de la cathode portée au potentiel nul, l’électron est accéléré à l’aide d’une tension V0 > 0 afin d’acquérir à l’entrée de la zone de déviation une vitesse ⃗v0, — pendant son trajet dans la zone de déviation, il est soumis à un champ électrique ⃗E lié aux potentiels ±Vp des plaques du condensateur, de longueur ℓ et séparées par une distance d, — poursuivant son mouvement, il arrive sur la surface de l’écran à une distance ys de l’axe x, l’écran étant situé à la distance D du centre du condensateur. 9/15 Soit un point P de coordonnées (x, y) ∈ [0, L[×[0, L[. Ce point est dans la cellule (i, j), où i = ⌊x/h⌋ et j = ⌊y/h⌋. Posons rx = x − ih et ry = y − jh. Q30. Montrer alors que : Ex(x, y) ≈ Ex[i,j] + ((Ex[i+1,j] - Ex[i,j]) * rx + (Ex[i,j+1] - Ex[i,j]) * ry) / h . Écrire de même la formule permettant de calculer Ey(x, y). En déduire qu’il est possible de calculer les composantes du champ en tout point de P. On supposera dans la suite que les fonctions val_Ex(Ex,Ey,x,y,h) et val_Ey(Ex,Ey,x,y,h) sont définies et retournent les valeurs des composantes du champ électrique pour le point M de coor- données (x, y) calculées à l’aide des formules précédentes. Calcul de la trajectoire par la méthode d’Euler Q31. Montrer que les équations permettant de décrire le mouvement de l’électron par la méthode d’Euler sont les suivantes, avec δt comme petit incrément temporel et δx, δvx, δy, δvy les varia- tions pendant δt des grandeurs x, vx, y, vy : δx = vxδt δy = vyδt δvx = − e mEx(x, y)δt δvy = − e mEy(x, y)δt . Q32. Compte tenu du changement de la position d’origine du repère (imposée par la résolution nu- mérique de l’équation de Poisson, figure 8), quelles sont les conditions initiales du mouvement de l’électron ? Déterminer δt pour calculer environ 200 points successifs le long de la trajectoire. Faire l’ap- plication numérique. Q33. En tenant compte des réponses aux questions précédentes, compléter le code d’initialisation des variables de la simulation (****** dans le code suivant) : Npts = 200 o nombre de points pour le tracé de la trajectoire v0 = ****** o vitesse initiale de l’électron dt = ****** o incrément temporel o tableaux des coordonnées x et y de l’électron lx = np.zeros(Npts) ; ly = np.zeros(Npts) o tableaux des vitesses en x et en y lvx = np.zeros(Npts) ; lvy = np.zeros(Npts) oconditions initiales lx[0] = ****** ; ly[0] = ****** lvx[0] = ******* ; lvy[0] = ******* Q34. Écrire les lignes de code implémentant la boucle de remplissage des tableaux lx, ly, lvx, lvy selon la méthode d’Euler. 12/15 Comparaison théorie/simulation La figure 9 montre le résultat de la simulation précédente. On y voit le réseau de courbes équipo- tentielles, ainsi que deux trajectoires 1 et 2 , l’une étant associée au calcul théorique, l’autre à la simulation numérique. Chaque trajectoire est constituée de 200 points de calcul séparés d’une durée δt. Q35. Reproduire sommairement sur la copie la figure 9, y ajouter le tracé de quelques lignes de champ orientées dans les différentes parties de la zone de déviation. Identifier, en le justifiant, chaque trajectoire. Expliquer pourquoi la trajectoire 1 est plus courte que la trajectoire 2 . À votre avis, peut-on se contenter de l’étude théorique pour prévoir le point d’impact de l’élec- tron sur l’écran ? (On attend une réponse chiffrée.) 1 2 Figure 9 – Tracé dans la zone de déviation de quelques courbes équipotentielles, des trajectoires de l’électron (théorique et simulée numériquement) 13/15 Aide-mémoire numpy/matplotlib/pyplot Importation des bibliothèques Les bibliothèques sont importées de la façon suivante : import math import numpy as np import matplotlib.pyplot as plt Manipulation des tableaux numpy La création d’un tableau numpy à deux dimensions dont toutes les valeurs sont initialisées à 0 est faite par l’instruction np.zeros(format), format étant un doublet de la forme (n_lignes , n_colonnes) : >> t0=np.zeros((2,3)); print(t0) [[ 0. 0. 0.] [ 0. 0. 0.]] Pour avoir un tableau rempli de 1, on utilise np.ones(format) : >> t1=np.ones((2,2)); print(t1) [[ 1. 1.] [ 1. 1.]] On peut récupérer le format d’un tableau en demandant son attribut shape, ce qui retourne un doublet : >> print(t0.shape); print(t1.shape) (2, 3) (2, 2) Dans le cas d’un tableau carré, on peut donc récupérer le nombre de lignes, égal au nombre de co- lonnes, en accédant au premier élément du doublet : >> t1.shape[0] 2 On peut créer un tableau numpy de booléens en ajoutant le type bool. La valeur 0 est associée à False, la valeur 1 à True : >> np.zeros((2,3), bool) [[False False False] [False False False]]} >> np.ones((2,3),bool) [[ True True True] [ True True True]] 14/15 L’accès à un élément du tableau a (en lecture ou en modification) se fait par a[i,j], les lignes et les colonnes étant numérotées à partir de 0 : >> a=np.zeros((2,3)) ; a[0,0]=1 ; a[1,2]=2 >> a [[ 1. 0. 0.] [ 0. 0. 2.]] Une copie indépendante d’un tableau a se fait à l’aide de np.copy(a) >> b = np.copy(a) ; b[0,0]=3 ; b[1,1]=5 >> a [[ 1. 0. 0.] [ 0. 0. 2.]] >> b [[ 3. 0. 0.] [ 0. 5. 2.]] FIN 15/15 I M P R I M E R I E N A T I O N A L E – 17 1185 – D’après documents fournis