Exercice 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.

Première « expérience »

Préparation

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. On peut (doit) ensuite un sous-dossier pour chaque expérience, qui contiendra le fichier de paramètres et les fichiers de sortie.

Commencez par créer un dossier pour cet exercice :

mkdir état
cd état

Puis téléchargez les programmes ici, ouvrez l'archive en cliquant dessus, et extrayez son contenu dans le dossier de l'exercice.

Pour préparer la première expérience, créez son dossier et copiez-y le fichier de paramètres run.pot:

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:

nano run.pot

Pour enregistrer et sortir de nano, utilisez Ctrl + X.

Simulation

Vous pouvez maintenant lancer la simulation

../evmd --file run.pot

Les deux points donnent le chemin de l'exécutable qui est dans le dossier parent du dossier de l'expérience. Il prend comme argument le fichier de paramètres, ici 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, où ls -l pour avoir leur taille. La commande less permet de voir le contenu d'un fichier (utiler Q pour sortir). 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 xmovie, qui affiche les positions des centres des sphères (projetées sur un plan). Il utilise les fichiers dump :

xmovie dump*

Création du notebook Jupyter

Nous allons utiliser Jupyter pour traiter les données contenues dans le fichier stress.dat. Nous utiliserons un seul notebook pour tout l'exercice, qui sera donc dans le dossier de l'exercice (et non dans le dossier d'une expérience). Dans Jupyter, naviguez jusqu'au dossier de l'exercice (ici « état ») et créez un notebook dans ce dossier.

Dans le notebook, chargez les modules de calcul numérique et de tracé de courbes avec

%pylab inline

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

a = 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, comme vous l'avez vu dans le tutoriel,

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

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 une fonction qui, à partir du numéro de l'expérience, 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 calculer et renvoyer la pression moyenne par cette fonction.

Exercice

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 xmovie et identifiez deux régimes.

Faites trois expériences pour une fraction volumique φ = 0.5, et tracez leur pression sur le même graphe. Que se passe-t-il pour φ = 0.5 (utilisez xmovie) ? 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 à t0 (que vous pouvez ajouter comme argument optionnel). Vous devrez utiliser la commande numpy.where ; pour comprendre son fonctionnement, analysez le code suivant

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

Équation d'état

Pour chaque fraction volumique, mesurez la pression et notez la dans un fichier (que vous créerez avec l'éditeur nano). Quand la pression varie fortement au cours de l'expérience, prenez la valeur du deuxième plateau. Ordonnez les fractions volumiques dans votre fichier. Ajoutez des points dans les zones où la variation de la pression vous semble singulière.

Comparaison au développement du viriel

Vous avez rencontré le développement du viriel dans les notes de physique statistique, où le calcul du premier ordre est détaillé, et où les coefficients sont donnés jusqu'à l'ordre 10. Le diamètre des sphères dures est σ = 1, et la fraction volumique φ est reliée à la densité ρ = N ⁄ V par φ = vρ, où v = 4π(σ ⁄ 2)3 ⁄ 3 = π ⁄ 6 est le volume d'une sphère.

Recupérez les valeurs du tableau dans ce fichier, puis écrivez une fonction qui renvoie la pression calculée par le développement du viriel à l'ordre n pour une fraction volumique φ.

Comparez le développement du viriel aux différents ordres au résultat des simulations.