1. Équation d'état

L'objectif de cet exercice est de déterminer l'équation d'état d'un systèmes de sphères dures en 3 dimensions, en traçant la pression en fonction de la fraction volumique.

Préparation

Téléchargez les programmes ici ; l'archive est placée dans votre dossier personnel, dans le dossier Téléchargements.

Ouvrez maintenant un terminal et tapez les premières commandes :

• pwd (suivie de Entrée) permet d'afficher le chemin absolu du répertoire courant, qui est par défaut votre répertoire personnel au lancement du terminal. Notez qu'il commence par /, qui est la racine du système de fichiers.
• ls affiche le contenu du répertoire courant.
• cd [chemin] permet de changer de répertoire courant en [chemin]. Le [chemin] peut être donné de façon absolue (en commençant par /) ou, ce qui est souvent plus pratique, de façon relative (au répertoire courant). Par exemple, si vous tapez cd Téléchargements, vous vous retrouverez dans le répertoire Téléchargements qui est dans votre répertoire personnel. Vous pouvez ensuite utiliser pwd pour voir le nouveau chemin et ls pour voir son contenu ; vous devriez y trouver l'archive equation_etat.tar. Voilà d'autres chemins que vous pouvez entrer :
  • ~ : il s'agit de votre répertoire personnel. cd ~ vous emmène dans votre répertoire personnel. Vous pouvez composer ces chemins, cd ~/Téléchargements vous emmène dans votre répertoire Téléchargements.
  • .. : il s'agit du répertoire situé « au-dessus » du répertoire courant dans le système de fichiers. Exécuté dans Téléchargements, vous serez envoyé dans votre dossier personnel. Vous pouvez utiliser pwd avant et après cette commande pour comprendre son comportement.

Des raccourcis clavier facilitent l'utilisation du terminal :

• La touche tabulation permet de compléter automatiquement ce que vous êtes en train de taper. Cela permet de gagner du temps et d'éviter des fautes de frappes. Placez-vous dans votre répertoire personnel et tapez cd Tél puis tabulation : Téléchargements sera complété automatiquement. L'auto-complétion ne complète que s'il y a une seule possibilité.
• La flèche haut permet de faire défiler l'historique des commandes déjà utilisées.

Votre répertoire courant est en fait sur le réseau, et non sur la machine que vous êtes en train d'utiliser. Pour des raisons d'espace de stockage et de performance, vous lancerez les simulations localement. Vous vous placerez pour cela dans le répertoire /scratch (avec cd /scratch). Ce répertoire est commun à tous les utilisateur de cet ordinateur ; créez-y votre répertoire avec la commande mkdir [nom du binôme] (évitez les espaces dans le nom). Entrez ensuite dans votre répertoire avec la commande cd.

Copiez l'archive dans votre répertoire puis extrayez les programmes avec les commandes

mv ~/Téléchargements/equation_etat.tar ./
tar -xvf equation_etat.tar

Dans la première ligne, le point indique le répertoire courant, où on souhaite déplacer l'archive. Déplacez-vous ensuite dans le dossier créé, equation_etat.

Première expérience

Création du dossier et lancement du programme

Une bonne organisation est d'avoir un dossier pour chaque exercice, qui contiendra les programmes exécutables et le notebook Jupyter où les données seront traitées, et un sous-dossier pour chaque expérience, qui contiendra le fichier de paramètres et les fichiers de sortie.

Créez le dossier de la première expérience et copiez-y le fichier de paramètres avec la commande cp :

mkdir exp1
cp run.pot exp1/
cd exp1

Pour modifier les paramètres de la simulation (par exemple la fraction volumique), vous pouvez éditer le fichier de paramètres (celui qui est dans le dossier de l'expérience) avec l'éditeur en ligne de commande nano, en tapant nano run.pot. Vous pourrez ensuite utiliser Ctrl + X pour sortir en enregistrant.

Vous pouvez maintenant lancer la simulation en entrant

../evmd.linux --file run.pot

La simulation peut être interrompue après quelques minutes avec Ctrl + C dans le terminal (cette commande sert, de manière générale, à arrêter l'exécution d'un programme).

Fichiers de sortie et visualisation

Vous pouvez voir les fichiers de sortie avec la commande ls (ls -l affiche aussi leur taille). La sortie est composée de deux types de fichiers :

• Les fichiers dump….gz contiennent les positions de toutes les particules à des instants différents.
• Le fichier stress.dat contient les données suivantes, organisées en colonnes : [temps] [pression] [sigma_xx] [sigma_yy] [sigma_zz] [sigma_xy] [sigma_xz] [sigma_yz] ([sigma_ij] sont les composantes du tenseur des contraintes). La commande less permet de voir son contenu (utiliser Q pour sortir).

L'expérience peut être visualisée avec le programme VMD, qui affiche les positions des centres des sphères (le rayon affiché n'est qu'une représentation). Pour le lancer, utilisez le script vmd_dump.sh:

bash ../vmd_dump.sh

Seule la tranche du système $0 < z < 12$ est affichée. Pour afficher toutes les sphères, entrez la commande mol modselect 0 0 all dans l'invite de commande de VMD.

Lancement de Jupyter

Pour lancer l'interface Jupyter, ouvrez un nouveau terminal et lancez Jupyter avec

conda activate psa
jlab

Naviguez ensuite jusqu'à votre répertoire de travail, /scratch/[nom du binôme].

• Lancer Jupyter dans votre répertoire de travail vous permettra d'avoir accès aux trois exercices à partir de la même interface.
• Il est important d'utiliser un nouveau terminal car celui-ci est mono-tâche : vous ne pouvez pas exécuter des programmes dans la fenêtre de terminal utilisée pour lancer Jupyter. Ne fermez pas cette fenêtre !

Dans Jupyter, vous pouvez naviguer jusqu'au dossier de l'exercice et y créer un nouveau Notebook en choissant l'environnement « conda : psa ». Vous utiliserez ce Notebook pour tout l'exercice.

Dans le notebook, commencez par charger les modules Numpy et Matplotlib avec

import numpy as np
import matplotlib.pyplot as plt

Tracé de la pression en fonction du temps

Pour utiliser les données du fichier stress.dat, il faut commencer par les charger avec la commande genfromtxt :

a = np.genfromtxt('exp1/stress.dat', skip_footer=1)

Cette commande charge le contenu du fichier dont le chemin est donné par le premier argument dans la variable a. L'option skip_footer=1 permet de ne pas charger la dernière ligne du fichier, qui peut n'être que partiellement remplie à cause de notre façon brutale d'arrêter la simulation.

Pour tracer la pression (deuxième colonne) en fonction du temps (première colonne), il suffit d'utiliser

plt.xlabel('t')
plt.ylabel('pression')
plt.plot(a[:,0],a[:,1])

Équation d'état et diagramme de phase

Vous êtes maintenant prêts pour simuler un grand nombre de fractions volumiques et obtenir l'équation d'état du système, c'est à dire sa pression $P$ en fonction de sa fraction volumique $\phi$.

Fonctions de traitement en Python

Vous allez devoir traiter un grand nombre d'expériences, il est donc préférable de vous organiser en conséquence en écrivant des fonctions pour les opérations que vous allez devoir appliquer à toutes les expériences.

Commencez par écrire la fonction pression(n) qui, à partir du numéro de l'expérience n, trace la pression en fonction du temps. La chaîne de caractères associée à un entier $n$ est str(n), et les chaînes de caractères sont concaténées avec le caractère +. Regardez par exemple ce que donne

n = 12
print('exp' + str(n) + '/stress.dat')

Faites aussi calculer et renvoyer la pression moyenne par cette fonction.

Pression en fonction du temps

Tracez la pression en fonction du temps pour vos expériences (seules une ou deux courbes suffiront dans le rapport). Visualisez les expériences avec VMD et identifiez deux régimes.

Faites trois expériences pour une fraction volumique $\phi=0.5$, et tracez leur pression sur le même graphe. Que se passe-t-il pour $\phi=0.5$ (utilisez VMD) ? La pression qui nous intéresse correspond au deuxième plateau. Pour la calculer, modifiez la fonction précédente pour calculer la pression moyenne pour des temps supérieurs à $t_0$, que vous pouvez ajouter comme argument optionnel en écrivant

def pression(n, t0=0):
    ...

Vous devrez utiliser la commande numpy.where ; pour comprendre son fonctionnement, analysez le code suivant

r = np.random.rand(10)
print(r)
i = np.where(r>0.5)[0]
print(i)
print(r[i])

Équation d'état

Pour chaque fraction volumique, mesurez la pression et notez-la dans un tableau de la forme :

eqetat = np.array([[0.56, 10.5, 1],
                   [0.1, 0.2, 2],
                   ...,
                   [0.7, 70, 12]])

La dernière colone contient le numéro de l'expérience ; cela permet de retrouver facilement les résultats si un point semble aberrant.

Quand la pression varie fortement au cours de l'expérience, prenez la valeur du deuxième plateau. Ajoutez des points dans les zones où la variation de la pression vous semble singulière.

Tracez la courbe $P(\phi)$.

Comparaison au développement du viriel

Le développement du viriel est un développement de la pression en puissances de la densité $\rho$: $$\frac{P}{k_\mathrm{B}T}=\sum_{i=1}^\infty (4v)^{i-1}B_i\rho^i.$$ Le volume $v$ d'une sphère de diamètre 1 est $v=\pi/6$. La densité est reliée à la fraction volumique par $\phi=v\rho$. Les coefficients du viriel $B_i$ ont été calculés jusqu'à l'ordre 10. Vous pouvez les charger en copiant la ligne suivante (B[0] correspond au coefficient $B_1$, etc) :

B = [1, 1, 0.625, 0.2869495, 0.110252, 0.03888198, 0.01302354, 0.0041832, 0.0013094, 0.0004035]

Écrivez la fonction pression_viriel(phi, n) qui renvoie la pression calculée par le développement du viriel à l'ordre $n$ pour une fraction volumique $\phi$. Si vous n'utilisez pas de fonction trop compliquée, l'argument phi peut être un tableau.

Comparez le développement du viriel aux différents ordres au résultat des simulations. Notez que $k_\mathrm{B}T=1$ dans les simulations.

Effet de la température

Vous avez obtenu $P(\phi)$, mais l'équation d'état demande en principe aussi la dépendance en la température, $P(\phi, T)$.

Pour déterminez cette dépendance, partez de la fonction de partition de ce système, qui peut s'écrire $$Z=\frac{1}{N!}\int \prod_i (d\mathbf{r}_i\, d\mathbf{v}_i) \exp\left(-\sum_i \frac{m v_i^2}{2 k_\mathrm{B}T}\right) \prod_{i<j}\theta(|\mathbf{r}_i-\mathbf{r}_j|-1)$$ L'intégrale sur les vitesses peut être calculée, ce qui donne $$Z=\left(\frac{2\pi k_\mathrm{B}T}{m} \right)^{3N/2}\frac{1}{N!}\int \prod_i d \mathbf{r}_i \prod_{i<j}\theta(|\mathbf{r}_i-\mathbf{r}_j|-1).$$ La fonction de partition intervient dans la définition thermodynamique de la pression : $$P=k_\mathrm{B}T\frac{\partial\log(Z)}{\partial V}.$$ Utilisez cette expression pour montrer que la pression est proportionnelle à la température.

De manière générale, que se passe-t-il dans ce système quand la température augmente ? Vous pourrez adopter une vision mécanique de la pression : elle est donnée par la quantité de mouvement échangée lors de collisions par unité de temps.