# Une équation de Langevin pour la simulation à l'échelle atomistique de fluides hors-équilibre

Matthew Dobson, Frédéric Legoll, Tony Lelièvre, Gabriel Stoltz

► To cite this version:

Matthew Dobson, Frédéric Legoll, Tony Lelièvre, Gabriel Stoltz. Une équation de Langevin pour la simulation à l'échelle atomistique de fluides hors-équilibre. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01722086

HAL Id: hal-01722086

https://hal.science/hal-01722086v1

Submitted on 2 Mar 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Une équation de Langevin pour la simulation à l'échelle atomistique de fluides hors-équilibre

Matthew DOBSON 1,2,4, Frédéric LEGOLL 3,4, Tony LELIEVRE 1,4, Gabriel STOLTZ 1,4

**Résumé** — La motivation de ce travail est la simulation multiéchelle de liquides, dans laquelle des modèles de dynamique moléculaire sont couplés avec des descriptions à l'échelle du continuum (typiquement via l'équation de Navier-Stokes). Dans un tel calcul, il est nécessaire de pouvoir simuler un système moléculaire consistant avec le flot macroscopique. Typiquement, on souhaite imposer un taux de déformation (gradient symétrisé de la vitesse) et calculer numériquement, sur la base d'un modèle atomistique, le tenseur des contraintes. On obtient ainsi numériquement, sur la base de calculs à l'échelle atomistique, la loi de comportement à l'échelle macroscopique du fluide.

Notre travail consiste tout d'abord à construire, en utilisant un modèle de bain thermique (et suivant en cela les idées de D. Dürr, S. Goldstein et J. Lebowitz [3, 4]), une équation de type Langevin, permettant la simulation de liquides à température et taux de déformation fixés. On présente ensuite des simulations numériques montrant l'utilité de la dynamique obtenue pour simuler des liquides en cisaillement.

**Mots clés** — systèmes hors-équilibre, simulation multiéchelle, équation de Langevin, flot de cisaillement.

# 1 Introduction

La motivation de ce travail est la simulation multiéchelle de liquides, dans laquelle des modèles de dynamique moléculaire sont couplés avec des descriptions à l'échelle du continuum. A l'échelle macroscopique, l'évolution du fluide (considéré ici comme incompressible) est modélisée par l'équation de Navier-Stokes :

$$ \rho \left( \frac{\partial u}{\partial t} + u \cdot \nabla u \right) = f + \operatorname{div} \sigma, \quad \operatorname{div} u = 0, \quad \sigma = -p \operatorname{Id} + \tau, \qquad (1) $$

où $u$ est le champ de vitesse, $\sigma$ le champ de contrainte et $p$ le champ de pression associé à l'incompresibilité du fluide ($\rho$ est la densité, $f$ modélise les forces extérieures et $Id$ est la matrice identité).

Pour fermer le système d'équations, on a besoin d'une loi constitutive, reliant typiquement $\tau$ (ou $\sigma$) à $u$. Pour les fluides newtoniens, cette loi s'écrit

$$ \sigma = -p \operatorname{Id} + \eta (\nabla u + \nabla u^T), $$

où $\eta$ est la viscosité. Le but de ce travail est de fermer le système (1) numériquement, sur la base d'un modèle atomistique (ce qui permet de contourner la difficulté de postuler à l'échelle macroscopique une loi constitutive). A plus long terme, ce travail est relié au problème de la simulation multiéchelle de fluides, pour laquelle on va résoudre

$$ \rho \left( \frac{\partial u}{\partial t} + u \cdot \nabla u \right) = f + \operatorname{div} \sigma, \quad \operatorname{div} u = 0, $$

et où, dans chaque (ou quelques) cellule (ou élément fini) macroscopique, la relation $\nabla u \mapsto \sigma(\nabla u, T)$ est calculée sur la base d'un modèle atomistique.

<footnote n="1">¹ CERMICS, Ecole des Ponts ParisTech, {lelievre,stoltz}@cermics.enpc.fr</footnote>

<footnote n="2">² Nouvelle adresse: Department of Mathematics and Statistics, University of Massachusetts, dobson@math.umass.edu</footnote>

<footnote n="3">³ Laboratoire Navier, Ecole des Ponts ParisTech, legoll@lami.enpc.fr</footnote>

<footnote n="4">⁴ INRIA Rocquencourt, MICMAC team-project</footnote>
La quantité $\nabla u$ étant macroscopique, on supposera dans toute la suite que, à l'échelle atomistique, elle est constante et uniforme.

La suite de ce document s'appuie sur l'article [2].

# 2 Equation de Langevin pour un fluide hors-équilibre

## 2.1 Modèle atomistique

A l'échelle atomistique, on modélise le fluide comme un système discret composé de N particules ponctuelles, de position $q_i \in \mathcal{D}$, de moment $p_i$ et de masse unité. On se donne un Hamiltonien pour ce système, $H(q,p) = \sum_{i=1}^{N} \frac{p_i^2}{2} + V(q_1, \dots, q_N)$, somme de l'énergie cinétique et de l'énergie potentielle.

Une dynamique communément admise pour simuler ce système à l'équilibre à la température T est l'équation de Langevin :

$$dq_i = p_i dt, \quad dp_i = -\nabla_i V(q) dt - \gamma p_i dt + \sqrt{2\gamma\beta^{-1}} dW_i, \quad (2)$$

où $\gamma > 0$ est un paramètre, $\beta$ est relié à la température par $\beta = \frac{1}{k_B T}$ ($k_B$ est la constante de Boltzmann) et $W_i$ est un mouvement brownien (on reconnait les équations de Newton auxquelles on a ajouté un terme de dissipation et un terme de fluctuation aléatoire).

Les propriétés d'équilibre du système sont définies comme la moyenne de certaines observables suivant la mesure de Boltzmann-Gibbs

$$\psi_G(q, p) = Z^{-1} \exp(-\beta H(q, p)). \quad (3)$$

Ces moyennes thermodynamiques s'écrivent

$$\langle A \rangle := \int_{\mathcal{D}^N \times \mathbb{R}^{3N}} A(q,p) \psi_G(q,p) dq dp. \quad (4)$$

Par exemple, l'observable correspondant à la pression est $A(q,p) = \frac{1}{3|\mathcal{D}|} \sum_{i=1}^{N} (p_i^2 - q_i \cdot \nabla_i V(q))$. Avec ce choix, $\langle A \rangle$ est la pression macroscopique d'un système de N particules, dans le domaine $\mathcal{D}$, à la température T.

La quantité (4) est une intégrale en très grande dimension. En pratique, elle est calculée comme la moyenne de l'observable A le long de la solution de (2) : pour tout $\gamma > 0$, et pour toute observable A, on a en effet

$$\lim_{T \to +\infty} \frac{1}{T} \int_0^T A(q(t), p(t)) dt = \langle A \rangle.$$

Le modèle qu'on vient de décrire, qui permet par exemple de calculer la pression en fonction de la densité et de la température, est un modèle d'équilibre. On a en particulier que la vitesse moyenne de chaque particule est nulle : $\langle p_i \rangle = 0$ pour tout i.

## 2.2 Modélisation hors-équilibre

On cherche à modifier le cadre précédent, pour que la vitesse moyenne au point q soit égale à $u(q)$, où $u$ est le champ de vitesse à l'échelle macroscopique (à divergence nulle) qu'on souhaite imposer. On se restreint ici au cas où $\nabla u$ est uniforme : $u(q) = \nabla u q$. Une difficulté est de contrôler correctement la température, d'une manière consistante avec le flot macroscopique imposé, et de telle façon à ce qu'il n'y ait pas de dérive énergétique. Cette question a déjà fait l'objet de nombreux travaux, cf. par exemple [6, 9, 10, 12].

Une idée naturelle est de remplacer l'équation de Langevin standard (2) par

$$dq_i = p_i dt, \quad dp_i = -\nabla_i V(q) dt - \gamma(p_i - \nabla u q_i) dt + \sqrt{2\gamma\beta^{-1}} dW_i. \quad (5)$$
Ceci revient à considérer une friction proportionnelle à la vitesse relative de la particule, différence entre sa vitesse $p_i$ et la vitesse macroscopique $u(q_i) = \nabla u q_i$ qu'on cherche à imposer au point $q_i$.

On montre dans [2] qu'on peut dériver une version simplifiée des équations (5) à partir d'un modèle de bain thermique (en adaptant les idées de D. Dürr, S. Goldstein et J. Lebowitz [3, 4]). Plus précisément, en considérant une unique particule (de masse unité) dans un bain thermique constitué de particules de masse $m$, on peut montrer que, lorsque $m$ tend vers 0, la dynamique de la particule marquée ($q_t, p_t$) converge vers

$$dq = p dt, \quad dp = -\gamma(p - \nabla u q) dt + \sqrt{2\gamma\beta^{-1}} dW, \quad (6)$$

où le coefficient $\gamma$ est donné par une formule analytique, et dépend des propriétés géométriques de la particule marquée et de certaines caractéristiques du bain thermique (on renvoie à [2] pour plus de détails). La dérivation de (5) dans le cas de plusieurs particules marquées est une question ouverte et très difficile.

La dérivation rigoureuse, à partir d'un modèle de bain thermique, de (6), justifie ainsi d'une certaine façon la dynamique (5). Notons aussi que le modèle (5) peut être obtenu (de manière heuristique) en considérant les équations g-SLLOD, souvent utilisées pour décrire des systèmes hors-équilibre (cf. par exemple [5]), et en les complémentant par des termes de thermostat.

On se concentre dans la suite de ce document sur les résultats numériques obtenus en utilisant la dynamique hors-équilibre (5) et on renvoie à [2] pour la dérivation rigoureuse de (6).

## 2.3 Résultats numériques

En général, la mesure invariante de (5) n'est pas analytiquement connue (à l'opposé de l'équation (2), pour laquelle on sait, sous des hypothèses raisonnables sur $V$, que la mesure $\psi_G(q, p)dqdp$, où $\psi_G(q, p)$ est donnée par (3), est invariante, quelque soit la valeur de $\gamma$). On peut en fait vérifier que la densité

$$Z^{-1} \exp \left[ -\beta \left( V(q) + \sum_{i=1}^{N} \frac{(p_i - \nabla u q_i)^2}{2} \right) \right],$$

construite en analogie avec la mesure de Gibbs (3) et en utilisant la vitesse relative $p_i - \nabla u q_i$ plutôt que la vitesse absolue $p_i$, n'est pas laissée invariante par (5) !

Cette absence de résultats analytiques justifie l'exploration numérique des propriétés de (5). L'objectif est d'étudier dans quelle mesure la vitesse moyenne des particules au point $q$ est égale à la vitesse du fluide en ce point, $\nabla uq$, et si on peut effectivement calculer numériquement la contrainte dans le système atomique, $\nabla u$ et $\beta$ étant fixés. Une étude similaire, mais non identique, a été menée dans [8].

On se place dans un cas simple, celui d'un flot de cisaillement :

$$\nabla u = \begin{bmatrix} 0 & s & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \quad \text{où s est un paramètre.} \quad (7)$$

On considère un système de $N = 1000$ particules intéragissant via un potentiel de Lennard-Jones :

$$V(q) = \sum_{i=1}^{N} \sum_{j=i+1}^{N} \phi_{LJ}(\|q_i - q_j\|)$$

avec

$$\phi_{LJ}(r) = \begin{cases} 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right) + c_1 r + c_2 & \text{si } r < r_{\text{cut}}, \\ 0 & \text{si } r \ge r_{\text{cut}}. \end{cases}$$

Les constantes $c_1$ et $c_2$ sont choisies telles que $\phi_{LJ}(r)$ et $\phi'_{LJ}(r)$ soient continues en $r_{\text{cut}}$ (en pratique, on travaille avec $r_{\text{cut}} = 2.6$ et $\epsilon = 1$). Le volume $|D|$ est tel que la densité est égale à 0.7.

On applique des conditions aux limites de type Lees-Edwards [1, 5] au système. Ces conditions généralisent les conditions aux limites périodiques, de manière consistante avec le flot de cisaillement imposé.
Avec ce choix de $\nabla u$, on s'attend à ce que le flot soit uniforme dans les directions $x$ et $z$. On définit donc des "tranches" $R_k$ d'épaisseur $\Delta y$ dans la direction $y$ :

$$R_k = \mathcal{D} \cap \{k\Delta y \le y \le (k+1)\Delta y\}. \quad (8)$$

Pour intégrer la dynamique (5), on utilise un schéma de splitting, en intégrant séparément la partie Hamiltonienne (grace au schéma de Verlet) et les autres termes (intégrés analytiquement). Le schéma s'écrit

$$ 
\begin{align*}
p_i^{n+1/2} &= p_i^n - \frac{\Delta t}{2} \nabla_i V(q^n), \\
q_i^{n+1} &= q_i^n + \Delta t p_i^{n+1/2}, \\
p_i^* &= p_i^{n+1/2} - \frac{\Delta t}{2} \nabla_i V(q^{n+1}), \\
p_i^{n+1} &= \alpha p_i^* + (1-\alpha) \nabla u q_i^{n+1} + \left(\frac{1-\alpha^2}{\beta}\right)^{1/2} G_i^n,
\end{align*}
$$

où $\alpha = \exp(-\gamma\Delta t)$ et $G_i^n \sim \mathcal{N}(0, \text{Id})$. En pratique, on travaille avec $\Delta t = 0.005$, $\beta = 1$ et $\gamma = 1$.

Pour tout temps $t_\ell = \ell\Delta t$, on définit la valeur moyenne de la vitesse dans la tranche $R_k$ par

$$\mathring{V}(t_\ell, k) = \frac{\sum_{n=0}^\ell \sum_{i=1}^N p_i^n \mathbf{1}_{R_k}(q_i^n)}{\sum_{n=0}^\ell \sum_{i=1}^N \mathbf{1}_{R_k}(q_i^n)}. \quad (9)$$

Sur la Figure 1, on montre la valeur de cette vitesse moyenne au temps final de la simulation $t_\ell = T = 500$ : on constate qu'elle est bien en accord avec le flot $u(q) = \nabla uq$ imposé. On montre sur la même figure la distance entre la vitesse moyenne calculée (au temps final) et le flot souhaité :

$$\mathrm{dist}(t_\ell, k) = \left\| \mathring{V}(t_\ell, k) - u_{\mathrm{bkgrd}}(R_k) \right\| = \left[ (\mathring{V}_x(t_\ell, k) - s y_{R_k})^2 + (\mathring{V}_y(t_\ell, k))^2 + (\mathring{V}_z(t_\ell, k))^2 \right]^{1/2}. \quad (10)$$

On vérifie bien que la vitesse moyenne obtenue est proche de celle souhaitée.

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Gauche : vitesse moyenne (9) – selon $x$ – des particules dans les tranches (8) pour un flot de cisaillement (7) avec $s = 0.05$. Droite : distance (10) en fonction de $y$.</figcaption></figure>

On montre sur la Figure 2 la variance des vitesses dans chaque tranche $R_k$ : la variance de la vitesse dans les trois directions est proche de $\beta^{-1}$, comme souhaité (la température s'interprète encore comme la variance de la vitesse).

Nous constatons donc que la vitesse des particules évoluant suivant (5) est proche (au moins en terme de moyenne et de variance) de celle qu'on cherche à imposer. On s'intéresse maintenant au calcul du tenseur des contraintes. On utilise pour cela la formule de Irving-Kirkwood [7] :

$$\sigma = \frac{1}{|\mathcal{D}|} \sum_{i=1}^{N} \left( (p_i - \nabla u q_i) \otimes (p_i - \nabla u q_i) + \frac{1}{2} \sum_{j=1, j \neq i}^{N} (q_i - q_j) \otimes f^{(ij)} \right) \quad (11)$$
<figure><img src="image_3.png" /><figcaption>Fig. 2 – Variance de la vitesse en fonction de y (on a moyenné en x, z et t) : vitesse selon x (gauche), vitesse selon y (droite). Les résultats pour la vitesse selon z sont similaires.</figcaption></figure>

où

$$f^{(ij)} = - \frac{\phi'_{LJ}(\|q_i - q_j\|)(q_i - q_j)}{\|q_i - q_j\|}$$

est la force exercée par l'atome j sur l'atome i. Le premier terme dans (11) est le terme de pression cinétique, qui permet de retrouver, dans le cas des gaz parfaits ($V \equiv 0$), la loi de Mariotte. Le second terme permet de prendre en compte les forces entre les particules. De manière standard, nous avons soustrait la vitesse moyenne locale dans le terme de pression cinétique (qui dépend donc de la vitesse relative et non absolue, cf. [7]).

On montre sur la Figure 3 le coefficient $\sigma_{12}$ du tenseur des contraintes (i.e. la contrainte de cisaillement), obtenue en moyennant l'observable $\sigma$ sur l'ensemble de la trajectoire. Plus précisément, nous avons considéré 10 trajectoires indépendantes pour chaque valeur $s = 0, 0.01, 0.02, \dots, 0.07$ dans (7) et pour chaque valeur du paramètre $\gamma = 0.1, 1.0$ et 10. Nous avons ensuite moyenné les résultats obtenus sur ces 10 réalisations, afin d'obtenir de petits intervalles de confiance.

On observe sur la Figure 3 une relation linéaire entre la contrainte de cisaillement $\sigma_{12}$ et le taux de cisaillement $s$, ce qui permet de définir une viscosité de cisaillement $\eta = -\frac{\sigma_{12}}{s}$. Pour $\gamma \le 1$, on obtient le résultat $\eta \approx 1.2175$, ce qui est consistant avec les valeurs données dans la littérature (cf. par exemple [11, Table IV]). On remarque aussi que les résultats dépendent a priori de $\gamma$ (même si, pour $\gamma = 0.1$ et $\gamma = 1$, ils sont très proches). Ceci n'est pas surprenant : on s'attend en effet à ce que la mesure invariante de (5) dépende de $\gamma$.

<figure><img src="image_4.png" /><figcaption>Fig. 3 – Contrainte de cisaillement $\sigma_{12}$ en fonction du taux de cisaillement $s$, pour plusieurs valeurs de $\gamma$.</figcaption></figure>

Les résultats numériques résumés ici montrent que la dynamique (5) donne des résultats tout à fait satisfaisants.
Références

[1] M. P. Allen et D. J. Tildesley, *Computer simulation of liquids*, Clarendon Press, New York, 1989.

[2] M. Dobson, F. Legoll, T. Lelièvre et G. Stoltz. *Derivation of Langevin dynamics in a nonzero background flow field*, preprint arXiv 1203.3773 et HAL 00680278.

[3] D. Dürr, S. Goldstein et J. L. Lebowitz. *A mechanical model of Brownian motion*, Comm. Math. Phys., 78(4):507–530, 1981.

[4] D. Dürr, S. Goldstein et J. Lebowitz. *A mechanical model for the Brownian motion of a convex body*, Z. Wahrscheinlichkeit, 62(4):427–448, 1983.

[5] D. J. Evans et G. P. Morriss. *Statistical mechanics of nonequilibrium liquids*, ANU E Press, Canberra, 2007.

[6] N. G. Hadjiconstantinou. *Discussion of recent developments in hybrid atomistic-continuum methods for multiscale hydrodynamics*, Bull. Pol. Acad. Sci-Te., 53(4):335–342, 2005.

[7] J. H. Irving et J. G. Kirkwood, *The statistical mechanical theory of transport processes. IV. The equations of hydrodynamics*, J. Chem. Phys., 18:817–829, 1950.

[8] R. Joubaud et G. Stoltz, *Nonequilibrium shear viscosity computations with Langevin dynamics*, Multiscale Model. Sim., 10:191–216, 2012.

[9] S. T. O'Connell et P. A. Thompson, *Molecular dynamics–continuum hybrid computations: A tool for studying complex fluid flows*, Phys. Rev. E, 52(6):R5792–R5795, 1995.

[10] W. Q. Ren et W. E. *Heterogeneous multiscale method for the modeling of complex fluids and micro-fluidics*, J. Comp. Phys., 204(1):1–26, 2005.

[11] R. Rowley et M. Painter, *Diffusion and viscosity equations of state for a Lennard-Jones fluid obtained from molecular dynamics simulations*, Int. J. Thermophys., 18(5):1109–1121, 1997.

[12] T. Werder, J. H. Walther et P. Koumoutsakos, *Hybrid atomistic-continuum method for the simulation of dense fluid flows*, J. Comp. Phys., 205(1):373–390, 2005.