4  Smoothed Particle Hydrodynamics

Figure 4.1: Écoulement d’un fluide autour d’un obstacle cylindrique à \(\mathrm{Re} = 100\), simulé avec la méthode Smoothed Particle Hydrodynamics.

Dans le TP4, nous avons vu comment simuler un système de particules en phase liquide ou gazeuse en imposant une température suffisamment élevée. Cette méthode de simulation des fluides n’est en revanche pas très adaptée à la simulation de fluides aux échelles macroscopiques (ou continues) :

C’est à cause de ces limitations que d’autres méthodes de simulations que la dynamique moléculaire sont préférées aux échelles continues. La méthode Smoothed Particle Hydrodynamics (hydrodynamique à particules lissées) ou SPH, permet la simulation de fluides aux échelles macroscopiques, tout en conservant des propriétés très similaires à la dynamique moléculaire et aux approches que l’on a développées dans les TP jusqu’à présent, c’est-à-dire :

En revanche, le point de départ de la méthode SPH n’est pas l’équation de Newton (même si elle finit par apparaître dans le développement de la méthode), mais celle de Navier-Stokes :

\[ \rho \frac{D \vec v}{D t} = -\nabla p + \mu \nabla^2 \vec{v} + \vec f_\mathrm{ext} \tag{4.1}\]

\(\vec v\) est la vitesse, \(p\) est la pression, \(\mu\) la viscosité dynamique, \(\rho\) la densité et \(\vec f_\mathrm{ext}\) les forces extérieures (par exemple la gravité). \(\nabla p\) est le gradient de \(p\) et \(\nabla^2\vec v = \mathrm{div}(\nabla \vec v)\) est le laplacien vecteur de \(v\) (la divergence de la matrice jacobienne de \(v\), en notation indicielle \((\nabla^2\vec v)_{i} = \vec v_{i,jj}\)).

La SPH est donc une technique pour discrétiser l’équation de Navier-Stokes, c’est-à-dire transformer l’équation différentielle aux dérivées partielles (4.1) en un ensemble d’équations algébriques (des équations ne comportant que des opérations de base, addition, multiplication, etc.).

Une excellente introduction aux principes de la SPH est donnée par Koschier et al. (2019) : ce chapitre théorique s’appuie sur leur travail de synthèse.

4.1 Ensemble de particules et approximation

Le principe de base de la SPH est de représenter un fluide continu comme un ensemble de particules \(\{\vec x_j\}_{j=1}^N\) de masses \(\{m_j\}_{j=1}^N\). Comme les particules sont ponctuelles, et que le fluide est continu, il faut pouvoir exprimer des champs, comme la vitesse ou la pression, entre les particules. Supposons que l’on veuille calculer la pression à un point \(\vec x\), et que chaque particule a une valeur de pression \(p_j\), alors la pression en \(\vec x\) est calculée en faisant une moyenne pondérée des pressions sur les particules qui se trouvent dans un rayon \(h\) autour de \(\vec x\).

Figure 4.2: Cercle de rayon \(h\) autour de \(\vec x\)

On interpole donc la pression1 \(p(\vec x)\) de cette façon :

\[ p(x) = \sum_{j = 1}^N p_j \frac{m_j}{\rho(\vec x_j)} W(|\vec x - \vec x_j|, h) \tag{4.2}\]

On remarque que la pondération associée à chaque \(p_k\) dépend :

  • de la masse de la particule \(m_j\)
  • de la densité du fluide, \(\rho(\vec x_j)\), évaluée à la position \(\vec x_j\), que l’on note simplement \(\rho_j\)
  • d’une fonction \(W(d, h)\), qui est une fonction de pondération, et qui dépend de la distance \(|\vec x - \vec x_j|\) et du rayon \(h\).
Note

Pour que l’interpolation fonctionne, il faut que \(W\) satisfasse certaines conditions :

  • que \(W\) soit normalisée, c’est-à-dire \(\int_{\mathbb R^n}W(|\vec x|, h)\,\mathrm{d}x = 1\)
  • que \(W\) tende vers la distribution de Dirac quand \(h\) tend vers zéro : \(\lim_{h\rightarrow 0}W(|\vec x|, h) = \delta(\vec x)\)
  • que \(W\) soit positive \(W(|\vec x|, h) \geq 0\)
  • que \(W\) soit nulle pour des distances plus grandes que \(h\), \(W(|\vec x - \vec y|, h) = 0\) pour \(|\vec x - \vec y| \geq h\) (c’est la propriété illustrée sur la Figure 4.2)

Il est commun de choisir \(W\) polynomiale, typiquement, en deux dimensions :

\[W(d, h) = \frac{40}{7\pi h^2} \times \begin{cases} 6\left(\left(\frac{d}{h}\right)^3 - \left(\frac{d}{h}\right)^2\right) + 1 & 0 \leq \frac{d}{h} \leq \frac{1}{2} \\ 2\left(1 - \frac{d}{h}\right)^3 & \frac{1}{2} < \frac{d}{h} \leq 1 \\ 0 & \text{sinon}\end{cases} \]

On remarquera que l’on a besoin de \(\rho_j\) pour l’interpolation. On peut donc réutiliser l’Équation 4.2 pour interpoler la densité :

\[ \rho(\vec x) = \sum_{j = 1}^N \rho_j \frac{m_j}{\rho_j}W(|\vec x - \vec x_j|, h) = \sum_{j=1}^N m_jW(|\vec x - \vec x_j|, h)\]

La densité est donc simplement la moyenne pondérée des masses dans un cercle de rayon \(h\) autour de \(\vec x\). On peut évaluer l’équation ci-dessus en \(\vec x_i\), la position de la particule \(i\), pour connaître \(\rho_i\):

\[ \rho_i = \sum_{j = 1}^N m_jW(|\vec x_i - \vec x_j|, h) \tag{4.3}\]

4.2 Calcul des forces de pression

Le premier terme qui nous intéresse dans l’Équation 4.1 est le gradient de la pression \(\nabla p\). On peut reprendre l’Équation 4.2 et calculer le gradient :

\[ \nabla p = \frac{\partial}{\partial \vec x} p = \sum_{j = 1}^N p_j \frac{m_j}{\rho_j} \frac{\partial}{\partial \vec x} W(|\vec x - \vec x_j|, h) \tag{4.4}\]

Parallèle avec la dynamique moléculaire

Bien que les quantités manipulées soient différentes, on remarque que l’Équation 4.4 a une structure familière : si l’on reprend l’équation des forces du chapitre sur les interactions, on peut écrire la force totale \(\vec f_i\) que subit la particule \(i\), comme :

\[ \vec f_i = \sum_{j = 1, j\neq i}^N -U_2'(|\vec x_i - \vec x_j|)\vec n_{ij} = \sum_{j=1, j\neq i}-\frac{\partial}{\partial \vec x_i} U_2(|\vec x_i - \vec x_j|) \tag{4.5}\]

La structure de la somme de l’Équation 4.5 est donc presque identique à l’Équation 4.4 : on peut interpréter la fonction \(W\) comme étant similaire à l’énergie d’une paire de particules \(U_2\). Pour évaluer le gradient de pression pour une particule \(\vec x_i\), on peut donc écrire :

\[\nabla p_i = \sum_{j = 1}^N p_j \frac{m_j}{\rho_j} W'(|\vec x_i - \vec x_j|, h) \vec n_{ij} \]

Important

Notez que l’Équation 4.4 n’exclue pas le terme \(j = i\) de la somme, contrairement à l’Équation 4.5.

L’équation ci-dessous, bien que correcte, souffre d’un défaut de précision et n’est pas souvent utilisée en pratique. On lui préférera l’approximation suivante (Price 2012), qui a presque la même structure :

\[ \nabla p_i = \sum_{j = 1}^N m_j \left(\frac{p_i}{\rho_i^2}+\frac{p_j}{\rho_j^2}\right) W'(|x_i - x_j|, h) \vec n_{ij} \tag{4.6}\]

Équation d’état

Il nous reste à exprimer la pression pour une particule, \(p_i\), en fonction de la densité du fluide \(\rho_i\). Cela se fait avec une équation d’état, \(p_i = p(\rho_i)\), typiquement :

\[ p(\rho) = k\left(\frac{\rho}{\rho_0} - 1\right)\]

\(k\) est une pression (qui joue le rôle de constante de raideur) et \(\rho_0\) est la densité du fluide au repos.

Important

On note immédiatement que cette relation entre pression et densité implique que le fluide est compressible. Les formulations incompressibles de la SPH sont bien plus complexes à établir, il est donc important de garder ce fait en tête. La constante \(k\) contrôle la compressibilité du fluide, plus \(k\) est élevé, plus le fluide approche la limite incompressible.

Étapes de calcul

Pour calculer les forces de pression, il faut donc procéder aux étapes suivantes :

  1. Calcul des listes de voisins, de \(d_{ij}\) et \(\vec n_{ij}\) avec une longueur de troncation \(h\)
  2. Calcul de la densité \(\rho_i\) pour chaque particule avec \[ \rho_i = \sum_{j} m_j W(d_{ij}, h)\]
  3. Calcul de la pression \(p_i\) pour chaque particule avec l’équation d’état
  4. Calcul du gradient de pression pour chaque particule \(\nabla p_i\) avec \[ \nabla p_i = \sum_j m_j \left(\frac{p_i}{\rho_i^2}+\frac{p_j}{\rho_j^2}\right) W'(d_{ij}, h) \vec n_{ij} \]

Densité des particules et forces de pression pour une distribution aléatoire, \(h = 6R\)\(R\) est le rayon d’une particule. Plot réalisé avec les fonctions plt.scatter et plt.quiver de matplotlib.

4.3 Calcul des forces de viscosité

Le second terme de l’équation de Navier-Stokes est issu de la viscosité d’un fluide newtonien. Ici, le laplacien vecteur est approximé par l’équation suivante (Koschier et al. 2019) :

\[ \nabla^2 \vec v_i = 12\sum_{j=1}^N \frac{m_j}{\rho_j} \frac{(\vec v_i - \vec v_j)\cdot (\vec x_i - \vec x_j)}{d_{ij}^2} W'(d_{ij}, h)\vec n_{ij} \tag{4.7}\]

De nouveau, cette équation a la même structure que pour le calcul des forces.

Astuce

En Python, le produit scalaire \((\vec v_i - \vec v_j)\cdot(\vec x_i - \vec x_j)\) peut se calculer avec la fonction np.einsum :

# i, j sont les listes de voisins (comme au TP5)
xij = positions[:, i] - positions[:, j]
vij = velocities[:, i] - velocities[:, j]
result = np.einsum('ij,ij->j', vij, xij)

Étapes de calcul

Le calcul des forces de viscosité suit une structure similaire aux forces de pression :

  1. Calcul des listes de voisins, de \(d_{ij}\) et \(\vec n_{ij}\) avec une longueur de troncation \(h\)
  2. Calcul de la densité \(\rho_i\) pour chaque particule avec \[ \rho_i = \sum_{j} m_j W(d_{ij}, h)\]
  3. Calcul des forces de viscosité avec l’Équation 4.7 et \(\mu\)

4.4 Calcul des sommes avec listes de voisins

La fonction assemble_forces introduite au TP5 peut être utilisée pour calculer les sommes dont on a besoin pour les forces de pression (Équation 4.6) et de viscosité (Équation 4.7).

En revanche, pour la densité, on devra utiliser np.bincount pour faire la somme :

# i, j sont les listes de voisins (comme au TP5)
Wij = W(dij, h)
rho = np.zeros_like(masses)
rho += 0.5 * np.bincount(j, weights=Wij * masses[j], minlength=natoms)
rho += 0.5 * np.bincount(i, weights=Wij * masses[i], minlength=natoms)
rho += masses * W(0, h)  # terme pour i == j

4.5 Conditions aux limites

Pour simuler un environnement infini, dont on a parfois besoin, par exemple pour la simulation d’un écoulement autour d’un obstacle (voir Figure 4.1), on a recours à des conditions périodiques aux limites, illustrées ci-dessous.

Conditions périodiques aux limites

Avec des conditions périodiques aux limites, le domaine d’intérêt est répété dans deux directions. Concrètement, cela implique que les particules qui “sortent” du domaine par un côté doivent “re-rentrer” du côté opposé, en conservant leur vitesse.

Les distances entre particules sont affectées par les conditions périodiques, puisqu’une particule proche d’un bord du domaine peut interagir avec les particules de l’autre côté du domaine.

Astuce

La fonction get_neighbour_info introduite au TP5 permet de calculer correctement les distances entre particules. Voyons un exemple d’un domaine de taille 100 \(\times\) 100 avec conditions périodiques :

domain = np.array([100., 100.])
periodic_boundary = np.array([True, True])

i, j, dij, rij = get_neighbour_info(positions, cutoff, 'ijdD',
                                    domain, periodic_boundary)

Les tableaux dij et rij contiennent les distances et vecteurs corrigés pour prendre la condition périodique en compte.

4.6 Intégration en temps

L’intégration en temps peut se faire avec Velocity-Verlet (voir le chapitre sur l’intégration en temps), avec cependant une modification due à la viscosité : dans Velocity-Verlet, la vitesse pendant l’itération est la vitesse au demi pas de temps, ce qui cause soucis puisque la vitesse contribue aux forces par l’Équation 4.7. Il faut donc utiliser une vitesse extrapolée \(\vec v'\) au pas de temps complet (voir étape 3 ci-dessous) uniquement pour le calcul de la force de viscosité (étape 6). Le schéma est donc le suivant :

  1. \(\vec v_{i+\frac{1}{2}} \gets \vec v_i + \frac{\Delta t}{2}\left(\frac{\vec f_i}{m}\right)\)
  2. \(\vec r_{i+1} \gets \vec r_i + \Delta t \vec v_{i+\frac{1}{2}}\)
  3. \(\vec v'_{i+1} \gets \vec v_{i+\frac{1}{2}} + \frac{\Delta t}{2}\left(\frac{\vec f_i}{m}\right)\) (la vitesse au pas de temps complet)
  4. Calcul des listes de voisins et de \(\rho_{i+1}\), \(p_{i+1}\)
  5. Calcul de \(\vec f^\text{pression}_{i+1}\)
  6. Calcul de \(\vec f^\text{viscosité}_{i+1}\) avec \(\vec v'_{i+1}\)
  7. \(\vec f_{i+1} \gets \vec f^\text{pression}_{i+1} + \vec f^\text{viscosité}_{i+1} + \vec f^\text{ext}\)
  8. \(\vec v_{i+1} \gets \vec v_{i+\frac{1}{2}} + \frac{\Delta t}{2}\left(\frac{\vec f_{i+1}}{m}\right)\)

4.7 Affichage de la vitesse

La fonction suivante permet la visualisation de la vitesse des particules dans Ovito.

Code
def write_exyz(fname, positions, velocities):
    import ase
    from ase.io import write

    positions = np.vstack([positions, np.zeros(positions.shape[1])])
    velocities = np.vstack([velocities, np.zeros(positions.shape[1])])

    atoms = ase.Atoms('H'*positions.shape[1],
                      positions=positions.T,
                      velocities=velocities.T)
    write(fname, atoms)
Koschier, Dan, Jan Bender, Barbara Solenthaler, et Matthias Teschner. 2019. « Smoothed Particle Hydrodynamics Techniques for the Physics Based Simulation of Fluids and Solids ». Eurographics 2019 - Tutorials, 41 pages. https://doi.org/10.2312/egt.20191035.
Price, Daniel J. 2012. « Smoothed Particle Hydrodynamics and Magnetohydrodynamics ». Journal of Computational Physics, Special Issue: Computational Plasma Physics, 231 (3): 759‑94. https://doi.org/10.1016/j.jcp.2010.12.011.

  1. La pression, ou tout autre champ scalaire, vectoriel ou tensoriel (par exemple la vitesse).↩︎