2 Potentiel d’interaction
On rappelle pour un ensemble de \(N\) particules, avec positions \(\{\vec r_i\}_{i=1}^N\), vitesses \(\{\vec v_i\}_{i=1}^N\), forces \(\{\vec f_i\}_{i=1}^N\) et masses \(\{m_i\}_{i=1}^N\), le principe fondamental de la dynamique :
\[ m_i\frac{{\mathrm{d}}\vec v_i}{{\mathrm{d}}t} = \vec f_i,\quad i\in\{1, \ldots, N\}. \]
On a vu dans le chapitre précédent comment résoudre cette équation numériquement, mais pas comment calculer la force \(\vec f_i\) qui s’exerce sur chaque particule. À cause de notre utilisation de l’algorithme Velocity-Verlet, la force de chaque particule ne peut dépendre que de la position des autres particules, et non de leur vitesse. Cela peut sembler restrictif, mais on verra que ce n’est pas une limitation importante pour simuler des solides et des fluides. En plus de cette restriction algorithmique, on va imposer que toutes les forces découlant des interactions entre particules dérivent de potentiels, c’est-à-dire qu’il existe une fonction scalaire \(U: \mathbb{R}^{3\times N} \rightarrow \mathbb{R}\) telle que :
\[ \vec f_i(\{\vec r_k\}) = -\frac{\partial U}{\partial \vec r_i}(\{\vec r_k\}), \]
où on note \(\{\vec r_k\}\) la position de toutes les particules (\(k = 1, \ldots, N\)).
Toute la physique que l’on souhaite intégrer à notre simulation est contenue dans \(U\), qui défini l’énergie potentielle totale du système (voir cette revue de potentiels inter-atomiques). La définition de \(U\) pour un système donné est une procédure délicate qui mène rarement à un résultat unique et en parfaite adéquation avec le système en question. Nous verrons dans ce cours quelles hypothèses sont raisonnables pour certains solides et fluides simples.
On remarque que \(\frac{\partial U}{\partial \vec r_i}\) est une fonction de la position de toutes les particules. Dans le cas général, on peut écrire (Müser, Sukhomlinov, et Pastewka 2023) :
\[ U(\{\vec r_k\}) = \sum_{k=1}^N U_1(\vec r_k) + \sum_{k < l}U_2(\vec r_k, \vec r_l) + \sum_{k < l < m}U_3(\vec r_k, \vec r_l, \vec r_m) + \ldots \]
\(U_1\) est l’énergie potentielle due à des forces extérieures (par exemple un champ gravitationnel uniforme, un champ magnétique, etc.). \(U_2\) définit l’énergie potentielle pour des paires de particules. Par exemple, l’énergie potentielle de gravité entre deux particules peut s’écrire \(U_2(\vec r_k, \vec r_l) = Gm_km_l / |\vec r_l - \vec r_k|\). Le système solaire peut donc être simulé uniquement en définissant l’énergie potentielle pour des paires de particules. Suivant cette logique, \(U_3\) définit l’énergie pour un triplet de particules. C’est l’énergie que l’on doit définir pour décrire une molécule comme \(\mathrm{CH}_4\) : l’énergie des atomes d’hydrogène dépend certes de la distance \(\mathrm{CH}\), mais aussi de la position des autres atomes d’hydrogène pour former un tétraèdre à l’équilibre, quand l’énergie potentielle est minimale. Dans ce cours, on se limitera dans un premier temps aux potentiels \(U_1\) et \(U_2\).
Ci-dessous quelques exemples de potentiels :
Potentiel pesanteur
Pour le cas d’un champ de gravité uniforme, on peut définir \(U_1\):
\[U_1(\vec r_k) = -m_k \vec g\cdot \vec r_k\]
Potentiel harmonique
L’énergie potentielle d’une paire de particules reliées par un ressort de rigidité \(k\) peut s’écrire :
\[U_2(\vec r_k, \vec r_l) = \frac{1}{2}k|\vec r_l - \vec r_k|^2\]
Potentiel de Coulomb
L’énergie d’interaction de deux particules chargées s’écrit comme le potentiel de gravité :
\[ U_2(\vec r_k, \vec r_l) = \frac{q_kq_l}{4\pi\varepsilon_0|\vec r_l - \vec r_k|} \]
En règle générale, l’énergie potentielle d’interaction d’une paire de particule obéit à deux relations d’invariance :
- L’énergie est invariante par translation : une paire de particule à une distance fixe a la même énergie, peu importe sa position dans l’espace
- L’énergie est invariante par rotation : une paire de particule à une distance fixe a la même énergie, peu importe l’orientation de la paire dans l’espace
Ces deux propriétés se traduisent mathématiquement par le fait que \(U_2(\vec r_k, \vec r_l)\) est seulement fonction de la distance entre les particules \(r_{kl} = |\vec r_l - \vec r_k|\). Démontrons que l’invariance par translation donne une dépendance en \(\vec r_{kl} = \vec r_l - \vec r_k\) uniquement.
L’invariance par translation veut que pour tout vecteur \(\vec v\) on ait :
\[U_2(\vec r_k + \vec v, \vec r_l + \vec v) = U_2(\vec r_k, \vec r_l).\]
Si l’on calcule la dérivée par rapport à \(\vec v\) de l’équation ci-dessus, on obtient :
\[ \frac{{\mathrm{d}}}{{\mathrm{d}}\vec v}\left(U_2(\vec r_k + \vec v, \vec r_l + \vec v)\right) = \frac{\partial U_2}{\partial \vec r_k}\mathbf{I} + \frac{\partial U_2}{\partial \vec r_l}\mathbf{I} = \vec 0\quad\Leftrightarrow\quad \frac{\partial U_2}{\partial \vec r_k} + \frac{\partial U_2}{\partial \vec r_l} = \vec 0 \tag{2.1}\]
On opère à présent le changement de variable suivant \(\vec a = \vec r_k + \vec r_l\) et \(\vec b = \vec r_k - \vec r_l\), et on définit le potentiel \(\tilde U_2\) tel que :
\[\tilde U_2(\vec a, \vec b) = U_2(\vec r_k, \vec r_l)\]
On peut donc écrire les dérivées partielles suivantes :
\[\begin{align*} \frac{\partial U_2}{\partial \vec r_k} & = \frac{\partial\tilde U_2}{\partial \vec a} + \frac{\partial\tilde U_2}{\partial \vec b},\\ \frac{\partial U_2}{\partial \vec r_l} & = \frac{\partial\tilde U_2}{\partial \vec a} - \frac{\partial\tilde U_2}{\partial \vec b}, \end{align*}\] et remplacer dans l’Équation 2.1:
\[ \frac{\partial U_2}{\partial \vec r_k} + \frac{\partial U_2}{\partial \vec r_l} = 2\frac{\partial\tilde U_2}{\partial \vec a} = \vec 0. \]
Puisque la dérivée de \(\tilde U_2\) par rapport à \(\vec a = \vec r_k + \vec r_l\) est nulle, \(\tilde U_2\) est donc constante par rapport à \(\vec a\) et dépend uniquement de \(\vec b = \vec r_k - \vec r_l\). On peut faire une démonstration similaire pour l’invariance par rotation.
Maintenant qu’une expression pour la force et pour l’énergie potentielle a été établie, essayons de calculer la force \(\vec f_i\) agissant sur une particule \(i\) pour des \(N\) particules indépendantes (\(U_2 = 0\)) et pour l’absence de forces extérieures (\(U_1 = 0\)).
2.1 Forces extérieures
Dans le cas où \(U(\{\vec r_k\}) = \sum_{k=1}^N U_1(\vec r_k)\), la force sur une particule à la position \(\vec r_i\) est simple à calculer :
\[ \vec f_i = -\frac{\partial U}{\partial \vec r_i}(\{\vec r_k\}) = -\sum_{k=1}^N \frac{\partial}{\partial \vec r_i} U_1(\vec r_k)\]
Le seul terme non-nul de la somme est pour \(k = i\), on a donc :
\[ \vec f_i = -\frac{\partial U_1}{\partial \vec r}(\vec r_i). \]
La force sur la particule \(i\) ne dépend que de \(\vec r_i\), les particules sont donc bien indépendantes.
Calculer la force pour le potentiel électrique \[U_1(\vec r_i) = \frac{q_0 q_i}{4\pi\varepsilon_0 |\vec r_i|}\]
2.2 Forces entre particules
On se place dans le cas où, par invariance, \(U_2: \mathbb{R}\rightarrow\mathbb{R}\) est fonction uniquement de la distance entre particules \(r_{kl} = |\vec r_{kl} | = |\vec r_l - \vec r_k|\) :
\[\begin{align*} U(\{\vec r_k\}) & = \sum_{k<l}U_2(\vec r_k, \vec r_l) = \sum_{k = 1}^N\sum_{l = k + 1}^N U_2(|\vec r_l - \vec r_k|)\\ \vec f_i & = -\sum_{k = 1}^N\sum_{l = k + 1}^N \frac{\partial}{\partial \vec r_i}U_2(|\vec r_l - \vec r_k|) \end{align*}\]
Calculons la dérivée par rapport à \(\vec r_i\) pour les termes de la double somme :
\[\frac{\partial}{\partial \vec r_i}U_2(|\vec r_l - \vec r_k|) = U_2'(r_{kl}) \frac{\partial r_{kl}}{\partial \vec r_i}. \]
La dérivée de \(r_{kl}\) par rapport à \(\vec r_i\) s’exprime sous la forme :
\[\begin{align*} \frac{\partial r_{kl}}{\partial \vec r_i} = \frac{\partial \sqrt{\vec r_{kl}\cdot \vec r_{kl}}}{\partial \vec r_i} = \frac{1}{2r_{kl}}\frac{\partial (\vec r_{kl}\cdot \vec r_{kl})}{\partial \vec r_i} &= \frac{1}{2r_{kl}}\left(\frac{\partial \vec r_{kl}}{\partial \vec r_i}\cdot \vec r_{kl} + \vec r_{kl} \cdot \frac{\partial \vec r_{kl}}{\partial \vec r_i}\right)\\ & = \frac{\vec r_{kl}}{r_{kl}} \frac{\partial \vec r_{kl}}{\partial \vec r_i}\\ & = \vec n_{kl} \left(\frac{\partial \vec r_{l}}{\partial \vec r_i} - \frac{\partial \vec r_{k}}{\partial \vec r_i}\right)\\ & = \vec n_{kl} (\delta_{li} - \delta_{ki}) \end{align*}\]
où \(\vec n_{kl} = \vec r_{kl} / r_{kl}\) est le vecteur \(\vec r_l - \vec r_k\) normalisé. On peut à présent exprimer la force \(\vec f_i\). Pour simplifier les notations, on pose \(\vec f_{kl} = -U_2'(r_{kl})\vec n_{kl}\) qui est la force qu’exerce la particule \(k\) sur la particule \(l\) (par la troisième loi de Newton, on a \(\vec f_{kl} = -\vec f_{lk}\)) :
\[\begin{align*} \vec f_i = -\sum_{k = 1}^N\sum_{l = k + 1}^N \frac{\partial}{\partial \vec r_i}U_2(r_{kl}) & = -\sum_{k = 1}^N\sum_{l = k + 1}^N U_2'(r_{kl})(\delta_{li} - \delta_{ki})\vec n_{kl}\\ & = \sum_{k = 1}^N\sum_{l = k + 1}^N \vec f_{kl}\delta_{li} - \sum_{k = 1}^N\sum_{l = k + 1}^N \vec f_{kl}\delta_{ki}\\ & = \sum_{k = 1}^N\sum_{l = k + 1}^N \vec f_{kl}\delta_{li} -\sum_{l = i + 1}^N\vec f_{il}\\ & = \sum_{k = 1}^{i-1}\sum_{l = k + 1}^N \vec f_{kl}\delta_{li} +\sum_{l = i + 1}^N\vec f_{li}\\ & = \sum_{k = 1}^{i-1} \vec f_{ki} +\sum_{l = i + 1}^N\vec f_{li}\\ & = \sum_{k = 1,\ k \neq i}^N\vec f_{ki} = \sum_{k = 1,\ k \neq i}^N U_2'(r_{ki})\vec n_{ik} \end{align*}\]
Logiquement, on trouve que la force totale sur la particule \(i\) est la somme des forces pour toutes les paires \((i, k)\), mais ce calcul montre la marche à suivre pour l’étape 3 de Velocity-Verlet :
- Pour chaque paire de particules \((k, l)\), calculer \(\vec f_{kl} = -U_2'(r_{kl})\vec n_{kl} = U_2'(r_{kl})\vec n_{lk}\)
- Pour chaque particule, calculer \(\vec f_i = \sum_{k=1,\ k\neq i}^N\vec f_{ki}\)
Puisqu’il faut faire une somme de \(N\) termes pour chacune des \(N\) particules, le calcul des forces tel qu’écrit ci-dessus, dit “calcul naïf”, a une complexité algorithmique de \(O(N^2)\), ce qui devient vite coûteux. On verra dans un prochain chapitre un moyen de calculer les forces en \(O(N)\) opérations.
Le TP 2 est une application de ce calcul pour des chocs entre boules de billard élastiques.