import matplotlib.pyplot as plt
import numpy as np
import numpy.typing as npt
def verlet_predictor(positions: npt.NDArray[float],
velocities: npt.NDArray[float],
masses: npt.NDArray[float],
forces: npt.NDArray[float],
dt: float):
"""Étape de prédiction de Velocity-Verlet"""
def verlet_corrector(positions: npt.NDArray[float],
velocities: npt.NDArray[float],
masses: npt.NDArray[float],
forces: npt.NDArray[float],
dt: float):
"""Étape de correction de Velocity-Verlet"""7 2/ Collisions et interactions
Dans ce TP, nous allons enfin introduire une physique pour les interactions entre particules. Nous allons commencer par de simples collisions élastiques qui vont nous permettre de calculer des forces entre les particules, et discuter du potentiel d’interaction harmonique.
Objectifs
- Calculer la force à partir d’un potentiel
- Calculer les distances pour toutes les paires de particules
- Calculer les forces associées à chaque paire
- Inclure le calcul des forces dans l’intégration des équations de mouvement
- Calculer une collision avec une paroi
7.1 Prélude : organisation du code
Dans le précédent TP, nous avons implémenté les deux fonctions ci-dessous :
Vous pouvez déplacer ces fonctions dans leur propre fichier my_functions.py. Dans le script que vous écrirez pour le présent TP (et les suivants), vous pourrez écrire le code suivant au début de votre script pour pouvoir utiliser ces fonctions. Voir Section 5.1.
import sys
sys.path.append('.')
from my_functions import *
Avertissement
Il faut que le fichier my_functions.py se trouve dans le même dossier que votre fichier de TP.
7.2 Contact entre sphères élastiques
Le problème du contact entre deux sphères est, dans le cas général1, trop complexe pour être résolu à chaque collision de sphère. On va donc poser un cadre qui conservera une formulation simple, prenons deux sphères de rayons \(R\), de propriétés élastiques \((E, \nu)\), supposons un contact sans frottement et sans adhésion, et le rayon de contact (l’aire de contact est un disque) petit par rapport aux rayons des sphères. Ce problème a été résolu en 1881 par Hertz, qui a pu calculer la distance de chevauchement (voir Figure 7.1) \(\delta = 2R - d_{12}\) où \(d_{12}\) est la distance entre les centres des deux sphères :
\[ \delta = \left(\frac{9P^2}{16R{E^*}^2}\right)^\frac{1}{3},\quad \mathrm{avec}\ E^* = \frac{E}{2(1 - \nu^2)} \]
\(P\) est la force totale qui joint les deux sphères. En inversant cette relation, on peut exprimer la force en fonction de la distance \(\delta\) :
\[ P = \frac{4}{3}E^* \delta\sqrt{R\delta}\]
Notons que dans la théorie de Hertz, la quantité \(\sqrt{R\delta}\) est le rayon de l’aire de contact. Une autre façon très commune d’exprimer la force de contact est via un ressort de rigidité \(k_\alpha = \alpha E^*R\), \(\alpha\) étant sans dimension et généralement petit :
\[ P = k_\alpha\delta \]
Comparons les courbes de ces équations adimensionnalisées :
Il est intéressant de noter que la courbe de Hertz a une dérivée nulle pour \(\delta = 0\), et que la rigidité du contact \(\frac{\partial P}{\partial \delta}\) augmente alors que la modélisation du contact par un ressort suppose une rigidité constante.
7.3 Quantités conservées
Dans toute simulation, il faut s’assurer que la quantité de mouvement est bien conservée au cours du temps. Dans notre cas, puisque l’on cherche à modéliser des chocs élastiques, il faut aussi que l’énergie cinétique soit conservée durant la simulation.
7.4 Matrice de distance, calcul naïf de l’interaction
Puisqu’on doit calculer les forces qui agissent entre toutes les particules (et qui dépendent de la distance, voir Figure 7.2), on doit calculer les distances entre toutes les particules. Définissons la position de 4 particules :
Avertissement
Le code ci-dessous est donné à titre d’exemple, il ne fait pas partie des consignes du TP.
r = np.array([
[0, 1, 0, 0.5], # positions en x
[0, 0, 1, 1], # positions en y
])
On cherche donc à calculer le vecteur \(\vec{r}_{ij} = \vec{r}_j - \vec{r}_i\) pour chaque \(i\) et chaque \(j\). Une façon de procéder serait d’écrire deux boucles imbriquées.
N = r.shape[1] # le nombre de particules
D = r.shape[0] # la dimension physique du problème (ici 2D)
rij = np.zeros([D, N, N]) # on définit la matrice des vecteurs entre particules%%timeit # <- ne pas copier dans un script
for i in range(N):
for j in range(N):
rij[:, i, j] = r[:, j] - r[:, i]23.5 μs ± 59.9 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)rijarray([[[ 0. , 1. , 0. , 0.5],
[-1. , 0. , -1. , -0.5],
[ 0. , 1. , 0. , 0.5],
[-0.5, 0.5, -0.5, 0. ]],
[[ 0. , 0. , 1. , 1. ],
[ 0. , 0. , 1. , 1. ],
[-1. , -1. , 0. , 0. ],
[-1. , -1. , 0. , 0. ]]])On peut vérifier que les distances sont bien celles attendues:
dij = np.linalg.norm(rij, axis=0) # calcul des normes pour chaque colonne de rij (axis=0 <=> norme sur les colonnes)
dijarray([[0. , 1. , 1. , 1.11803399],
[1. , 0. , 1.41421356, 1.11803399],
[1. , 1.41421356, 0. , 0.5 ],
[1.11803399, 1.11803399, 0.5 , 0. ]])Cependant, les boucles for i in range(N) sont très lentes en python. Comme on utilise Numpy, on peut se servir de la fonctionnalité de “broadcast” pour éviter d’écrire des boucles:
%%timeit # <- ne pas copier dans un script
rij_ = r[:, np.newaxis, :] - r[:, :, np.newaxis]2.09 μs ± 14.7 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)On a gagné un facteur 10 sur la vitesse d’exécution ! Vérifions tout de même que le résultat est correct.
rij_ = r[:, np.newaxis, :] - r[:, :, np.newaxis]
print(np.all(rij == rij_))True
Consigne
En utilisant le code ci-dessus, implémentez la fonction ci-dessous :
def compute_dij_nij(positions: npt.NDArray[float]):
"""Calcule les distances et vecteurs normaux pour des paires de particules"""
dij = ... # distances entre paires de particules ||r_ij||
nij = ... # vecteurs normaux r_ij / ||r_ij||
return dij, nijAttention : la matrice dij calculée ci-dessus contient des zéros sur la diagonale. Le calcul de \(\vec n_{ij} = \vec r_{ij} / \|\vec r_{ij}\|\) contient donc des divisions par zéro qui produisent des NaN (Not a Number, pas un nombre). Ces NaNs sont des nombres flottants spéciaux qui se propagent quand on fait des calculs qui les contiennent. Ils sont très utiles pour déboguer un code car ils signalent une erreur probable. Dans notre cas, en revanche, on doit s’en débarrasser avant de procéder à d’autres calculs. On peut le faire avec la fonction np.nan_to_num comme ceci :
np.nan_to_num(nij, copy=False, nan=0.)Vérifiez que ce que retourne la fonction compute_dij_nij ne contient pas de NaN.
7.5 Potentiels d’interaction
On va dans un premier temps choisir la relation linéaire entre \(P\) et \(\delta\). En supposant que toutes les sphères aient le même rayon \(R\), on peut écrire l’énergie potentielle d’une paire de sphères sous la forme :
\[ U_2(d_{ij}) = \frac{1}{2}k_\alpha (2R - d_{ij})^2H(2R - d_{ij}) = \begin{cases} \frac{1}{2}k_\alpha (2R - d_{ij})^2 & \text{si } 2R \geq d_{ij}\\ 0 & \text{sinon}\end{cases} \] où \(H:\mathbb{R}\rightarrow \mathbb{R}\) est la fonction d’Heaviside :
\[ H(x) = \begin{cases} 1 & \mathrm{si}\ x \geq 0\\ 0 & \mathrm{sinon} \end{cases} \]
On utilise la fonction d’Heaviside dans l’énergie potentielle pour exprimer le fait que si deux sphères ne se touchent pas (si \(d_{ij} > 2R\)) alors l’énergie potentielle d’interaction est nulle. On implémente maintenant le calcul de la force dérivée2 au-dessus.
Consigne
def elastic_bounce_force(dij: npt.NDArray[float], R: float, ka: float):
"""Calcule l'intensité de la force de rappel élastique entre les particules, -U_2'(dij)"""
return ...Implémentez la fonction ci-dessus en utilisant la fonction \(U_2'\) dérivée plus haut. Pour intégrer la fonction d’Heaviside dans le calcul, on peut se servir du tableau de booléens suivant :
H = (2 * R - dij) > 0 # H contient p.ex. [True False ... True]Multiplier par H revient à multiplier par la fonction d’Heaviside \(H(2R - d_{ij})\).
Enfin, on peut combiner le calcul de rij et le calcul des forces en une fonction (sans oublier que \(\vec f_i = \sum_{k=1,\ k\neq i}^N \vec f_{ki}\), que l’on peut calculer avec np.sum) :
def compute_forces_collision(positions: npt.NDArray[float],
R: float, ka: float):
"""Calcule les forces sur chaque particule"""
# Étape 1: calcul de dij et nij
# Étape 2: calcul de -U_2'(dij) avec la fonction elastic_bounce_force
# Étape 3: calcul de fij = -U_2'(dij) nij
# Étape 4: somme pour calculer fi
return ...Pour finir, intégrez l’appel au bon endroit de la boucle de Velocity-Verlet.
Afin de tester votre implémentation, vous pourrez calculer à la main la force entre deux particules de rayon \(R\) = 1, raideur \(k_\alpha\) = 1, situées à une distance \(r\) = 1.5 l’une de l’autre, et comparer à votre implémentation.
7.6 Ouverture de billard
Afin de tester une situation plus complexe, nous allons simuler une ouverture (simplifiée) de billard. On donne les positions et vitesses initiales sous forme de fichiers texte :
Consigne
Utilisez la fonction np.loadtxt pour lire les fichiers3. On définit les paramètres du système suivants (attention aux unités, vérifiez les unités des fichiers des positions et vitesses):
- Rayons des particules : \(R\) = 2.8 cm
- Raideur des collisions : \(k_\alpha\) = 1 N/cm
- Masses des particules : \(m\) = 0.16 kg
Calculez le pas de temps critique \(\Delta t_\mathrm{crit} = \frac{2}{\omega}\).
Simulez l’évolution du système sans la force de gravité (suffisamment longtemps pour une collision, cf. Figure 7.3) et visualisez la trajectoire avec Ovito et la fonction write_xyz utilisée au TP précédent.
Vérifiez que la quantité de mouvement avant la collision est conservée après la collision. Vous pourrez implémenter une fonction momentum de la masse et de la vitesse qui fait le calcul de la quantité de mouvement.
Calculez l’énergie cinétique du système à la fin de chaque pas de temps, et affichez l’évolution de l’énergie cinétique au cours du temps. Cette dernière est-elle constante ? Essayez différentes valeurs du pas de temps et commentez la différence.
7.7 Contact élastique avec un mur
Le contact élastique d’une particule avec un objet immobile est relativement simple à traiter si l’on admet d’approximer l’instant et l’endroit de la collision. Puisque le choc est élastique, les énergies cinétiques avant et après choc doivent être égales. Pour une particule avec une vitesse incidente (avant choc) \(\vec v_i\), seule la composante du vecteur vitesse qui est normale à la surface rigide est affectée. On note \(\vec v_r\) la vitesse réfléchie. Dans le cas d’une surface plane, on a le schéma suivant :
Puisque la composante de vitesse parallèle au mur \(v_{i,x} = |\vec v_i|\cos(\theta_i) = v_{r,x}\) n’est pas affectée, on a \(\theta_i = \theta_r\), comme pour la réflexion d’une onde plane. La composante verticale de la vitesse réfléchie est donnée par \(v_{r,y} = -v_{i,y}\).
Consigne
En reprenant le code de l’ouverture de billard, implémentez les collisions avec un sol plat situé en contrebas de la position initiale des particules, et ajoutez la force de gravité à la force due aux chocs élastiques. Pour chaque particule, il faudra déterminer si la particule est au-dessus ou en dessous du plan du mur, puis modifier le vecteur vitesse si besoin.
Note : il faut déterminer quel est l’endroit, parmi les étapes de l’algorithme Velocity-Verlet, où il est judicieux de modifier la vitesse. Essayez les différentes possibilités et sélectionnez le résultat le plus physique.