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.
7.1 Prélude : organisation du code
Dans le précédent TP, nous avons implémenté les deux fonctions ci-dessous :
Important
Notez l’apparition du dt dans les paramètres des fonctions. Vérifiez votre code pour que ce dt apparaisse bien et soit correctement utilisé dans le corps de la fonction.
Vous pouvez copier ou déplacer ces fonctions dans leur propre fichier verlet.py. Dans le script que vous écrirez pour le présent TP (et les suivants), vous pourrez écrire from verlet import * pour pouvoir utiliser ces fonctions.
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_1\) et \(R_2\), de propriétés élastiques \((E_1, \nu_1)\) et \((E_2, \nu_2)\), 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 = R_1 + R_2 - 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}\ \frac{1}{R} = \frac{1}{R_1} + \frac{1}{R_2},\quad \mathrm{et}\ \frac{1}{E^*} = \frac{1 - \nu_1^2}{E_1} + \frac{1 - \nu_2^2}{E_2} \]
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. Quel que soit le choix de modélisation, il faut s’assurer que l’énergie totale est conservée lors de la collision (puisque l’on suppose un choc élastique).
7.3 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 (qui forment un tétraèdre) :
Avertissement
Le code ci-dessous est donné à titre d’exemple, réfléchissez avant de copier-coller dans votre script.
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]
D = r.shape[0]
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]26.9 µs ± 72.4 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
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.48 µs ± 52.6 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.4 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 ont 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}), \] où \(H:\mathbb{R}\rightarrow \mathbb{R}\) est la fonction de Heaviside :
\[ H(x) = \begin{cases} 1 & \mathrm{si}\ x \geq 0\\ 0 & \mathrm{sinon} \end{cases} \]
On utilise la fonction de 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ée 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"""
return ...Implémentez la fonction ci-dessus en utilisant la fonction \(U_2'\) dérivée plus haut. On pourra se servir de la fonction np.heaviside pour distinguer les cas où \(2R - d_{ij} \geq 0\) et \(2R - d_{ij} < 0\), ou encore 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 de 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 à compute_forces dans la boucle de Velocity-Verlet du TP précédent.
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.5 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 fichiers. On définit les paramètres du système suivants (attention aux unités):
- 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.
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.6 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.
qui prend en compte la rugosité, le frottement, l’adhésion, etc.↩︎