# Les simulations de dynamique moléculaire comme outil pour une meilleure compréhension de la physique locale du contact

Mathieu Solar, Christian Gauthier, Hendrik Meyer, Olivier Benzerara, Hervé

Pelletier, Christophe Fond, Jörg Baschnagel, Robert Schirrer

► To cite this version:

Mathieu Solar, Christian Gauthier, Hendrik Meyer, Olivier Benzerara, Hervé Pelletier, et al. Les simulations de dynamique moléculaire comme outil pour une meilleure compréhension de la physique locale du contact. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01412113

HAL Id: hal-01412113

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

Submitted on 8 Dec 2016

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 CC BY 4.0 - Attribution - International License
# Les simulations de dynamique moléculaire comme outil pour une meilleure compréhension de la physique locale du contact

M. Solar, C. Gauthier, H. Meyer, O. Benzerara, H. Pelletier, C. Fond, J. Baschnagel, R. Schirrer

ICS, CNRS - UPR 22, 23 rue du Loess, BP 84047, F-67034 Strasbourg Cedex 2
{solar, hmeyer, gauthier}@ics.u-strasbg.fr

## Résumé —

Les simulations de dynamique moléculaire se révèlent intéressantes dans l'étude de la physique locale du contact. En effet, cette technique de modélisation offre d'une part, une modélisation plus réaliste de la matière, et d'autre part, elle utilise une formulation plus physique que phénoménologique, donc plus thermodynamique. L'utilisation d'une telle technique de modélisation se révèle donc complémentaire à l'approche proposée par la MMC. Nous proposons une première étude de l'orientation des chaînes de polymères sous la pointe, lors d'un contact normal.

**Mots clés** — film, polymère, contact normal, dynamique moléculaire, orientation, chaîne.

## 1 Présentation générale

L'amélioration de la résistance à la rayure des surfaces de polymères (massifs ou en couches minces), passe par la réduction du coefficient de frottement, par la recherche d'un minimum de déformation plastique et d'une recouvrance après rayure. Il se révèle actuellement important de différencier le comportement de surface de celui du volume en étudiant la zone confinée sous le contact, à l'interface entre l'indenteur et le polymère, en indentation comme en rayure [1, 9, 3, 4]. Du fait également de leurs propriétés viscoélastiques et de leur comportement en cisaillement confiné (frottement) complexe, l'amélioration du comportement surfacique des polymères nécessite une meilleure compréhension de la physique locale du processus d'indentation et de rayure. Le besoin d'étudier la physique locale du contact reste donc entier, et, le recours aux simulations de dynamique moléculaire DM trouve ici pleinement son sens. En effet, cette technique de simulation permet de modéliser la matière avec des détails moléculaires. De part sa formulation thermodynamique "implicite", cette dernière permet d'observer d'éventuels changements fondamentaux dans les propriétés de la matière. Cependant, cette technique reste très consommatrice en temps de calcul, et, les échelles de temps et d'espace simulées restent encore très petites par rapport à celles de la MMC. Les premiers travaux réalisés ont montré des similitudes entre expérimentations (tests de nanoindentation), modélisation par EF et modélisation par DM. Ce papier présente une première étude sur l'orientation des chaînes macromoléculaires sous l'indenteur pendant un contact normal.

## 2 Modélisation du contact par simulations de DM

Les simulations de DM, réalisées ici sous MD Spherical, permettent de simuler l'évolution au cours du temps de la matière, sous d'éventuelles sollicitations mécaniques, par intégration numé-
rique des équations du principe fondamental de la dynamique de Newton (algorithme de Verlet). Une description plus complète de cette technique de simulation est faite dans les ouvrages suivants [10] [11]. La technique de nano-indentation étant le plus souvent à la limite du VER d'un polymère, l'utilisation de la DM est justifiée. Quelques travaux précurseurs existent sur l'étude de l'essai de nano-indentation [5]. Dans notre cas, un polymère est modélisé par un modèle mathématique générique type "coarse-grained" (ou granulation de la matière). Une particule désigne dans ce cas un ou plusieurs groupe(s) d'atome(s) d'une macrochaîne, et est considérée comme étant une sphère quasi-impénétrable, de diamètre $\sigma$ et de masse $m$. Dans ce travail, l'indentation avec un cône parfait a été simulée sur des films de polymères amorphes linéaires, donc sans ramifications de nature chimique entre les chaînes macromoléculaires. Les noeuds sont de nature physique (ou topologique) uniquement. En simulations de DM, les résultats sont présentés en unités adimensionnées de Lennard-Jones (Unités LJ, symboles avec exposants *), et les plus utilisées sont rassemblées dans le tableau 1, où sont également donnés des ordres de grandeurs pour un PMMA<fnref n="1" />. Pour alléger les notations, l'exposant * sera volontairement omis dans cette partie, sauf si une éventuelle ambiguïté subsiste.

## 2.1 Le modèle mathématique de chaînes : modèle de chaînes flexibles de polymères

Les films de polymère simulés sont constitués dans notre cas de 192 chaînes de 64 particules, reposant sur un substrat infiniment rigide (cf. Fig. 1). Une caractérisation complète de ce film- modèle a été réalisée par Peter et al. [6]. Les indentations simulées utilisent un cône parfait de demi-angle au sommet $\theta$. La taille de la boîte de simulation utilisée est de $(30\sigma \times 30\sigma \times 30\sigma)$. La taille du film est de $(30\sigma \times 30\sigma \times \approx 15\sigma)$. Les conditions aux limites périodiques (incluant la condition d'image minimale) sont appliquées au système suivant les axes $\vec{X}$ et $\vec{Z}$. La simulation se fait dans le système thermodynamique NVT : le nombre $N$ de particules dans la boîte reste constant ; le volume $V$ de la boîte reste constant, même si le film peut se dilater (son épaisseur augmenté avec la température) ; la température $T$ moyenne de la boîte reste constante (la régulation en température se fait avec un thermostat DPD (Dissipative Particule Dynamics) [7]). La surface supérieure du film est laissée libre, de sorte que le film s'arrange librement en surface (pression $P_h = 0$ en surface). Dans ce travail, le modèle de chaînes flexibles est utilisé comme modèle générique grain - ressort pour représenter les chaînes de polymère linéaire. La connectivité des chaînes est assurée par un potentiel d'élongation harmonique, noté $U_{bond}$, (1) ($k_{bond} = 1111 \cdot \varepsilon/\sigma^2$ est la raideur de la liaison, $R_i$ est la longueur de la liaison et $r_0 = 0.967\sigma$ est sa longueur à l'équilibre). Les particules d'une même chaîne, qui ne sont pas directement liées, et les particules appartenant à des chaînes différentes interagissent via un potentiel de Lennard-Jones 12–6, noté $U_{LJ}$, (1) tronqué ($r_c = 2.3\sigma$ est le rayon de coupure) et décalé ($C = 0.02684\sigma$ pour assurer la continuité du potentiel). La pointe est modélisée par un potentiel $U_{LJ}$ (zone interdite par l'indenteur), où $r$ devient la distance normale entre la surface de la pointe et une particule du film de polymère. Les valeurs suivantes ont été prises $r_c = 1.12\sigma$ et $C = 0.99982\sigma$, si bien que l'interaction est purement répulsive. Ceci associé au fait que la pointe est parfaitement lisse, le contact modélisé est parfaitement glissant.

$$ U_{bond}(R_i) = \frac{k_{bond}}{2} \cdot (R_i - r_0)^2 \quad \text{et} \quad U_{LJ}(r) = \begin{cases} 4\varepsilon \cdot [(\sigma/r)^{12} - (\sigma/r)^6] + C & \text{si } r < r_c \\ 0 & \text{si } r \ge r_c \end{cases} \quad (1) $$

<footnote n="1">¹Les valeurs suivantes ont été prises : $m = 1.660 \cdot 10^{-25} [kg]$, ce qui correspond à une masse molaire de $M = 100 [g \cdot mol^{-1}]$ [2]. À noter également que $k_b$ est la constante de Boltzmann : $k_b \approx 1.380 \cdot 10^{-23} [J \cdot K^{-1}]$, $\tau$ est l'unité de temps avec $\tau = \sqrt{(m \cdot \sigma^2)/\varepsilon}$, et enfin que $\sigma$ and $\varepsilon$ sont aussi souvent notés $\varepsilon_{LJ}$ and $\sigma_{LJ}$.</footnote>
## 2.2 Les outils pour l'étude de l'orientation des chaînes

Lorsque la pointe s'enfonce dans le film, l'orientation des chaînes est modifiée à plus ou moins longue distance. L'influence de l'enfoncement de l'indenteur sur l'orientation des liaisons est évaluée en utilisant le polynôme de Legendre d'ordre 2 (il recherche une symétrie à 180° sans distinction du sens entre les vecteurs ; les polynômes de Legendre servent au développement d'une fonction en harmoniques sphériques). Pour quantifier l'influence de l'enfoncement de la pointe, une moyenne 2 est réalisée sur toutes les valeurs prises par $P_2(X)$ dans un certain domaine du film, ce qui exige de discrétiser le film. Quant à la variable $X$, elle représente le produit scalaire entre un vecteur de liaison $\vec{u}_i$ (normalisé implicitement dans le calcul), qui lie la particule $i$ à la particule $i+1$ ($\vec{u}_i = i(i+1)$) et un certain vecteur de référence. Ainsi, la valeur de $\langle P_2(X) \rangle$ varie de -0.5 (tous les vecteurs sont perpendiculaires au vecteur de référence) à 1 (tous les vecteurs sont parallèles au vecteur de référence), en passant par 0 (pas d'orientation préférentielle, ce qui correspond à une structure amorphe).

**Choix de la variable X**
L'étude se situe dans la base locale à l'indenteur parfaitement conique (à base circulaire). Le point $A_i$ (3) est la projection orthogonale du point $X_i$ (2) sur le cône 3. Un vecteur de liaison est repéré par son vecteur position $\vec{OX}_i = \vec{X}_i$.

$$ \vec{OX}_i = X[i] \cdot \vec{X} + Y[i] \cdot \vec{Y} + Z[i] \cdot \vec{Z} \quad (2) $$

$$ \vec{OA}_i = X[A_i] \cdot \vec{X} + Y[A_i] \cdot \vec{Y} + Z[A_i] \cdot \vec{Z} \quad (3) $$

$$ k = \frac{Y[i] \cdot \tan \theta - \tan \theta \cdot \delta_p}{r[i] \cdot (1 + \tan^2 \theta)} + \frac{\tan^2 \theta}{1 + \tan^2 \theta} \quad (4) $$

Les vecteurs de référence sont ceux de la base locale 4, à savoir $\vec{U}_n$ (5), $\vec{U}_g$ (6) et $\vec{U}_t$ (7) (cf. Fig 2). L'étude de l'orientation des liaisons sous la pointe sera jugée par les produits scalaires suivants : $X = \cos \alpha^{gen} = \vec{u}_i \cdot \vec{U}_g / \|\vec{u}_i\|$, $X = \cos \alpha^{nor} = \vec{u}_i \cdot \vec{U}_n / \|\vec{u}_i\|$ et $X = \cos \alpha^{tan} = \vec{u}_i \cdot \vec{U}_t / \|\vec{u}_i\|$.

$$ \vec{U}_n = \left\{ X[i] \cdot \vec{X} - \tan \theta \cdot r[i] \cdot \vec{Y} + Z[i] \cdot \vec{Z} \right\} / \left\{ r[i] \cdot \sqrt{1 + \tan^2 \theta} \right\} \quad (5) $$

$$ \vec{U}_g = \tan \theta \cdot \left\{ X[i] \cdot \vec{X} + (r[i] / \tan \theta) \cdot \vec{Y} + Z[i] \cdot \vec{Z} \right\} / \left\{ r[i] \cdot \sqrt{1 + \tan^2 \theta} \right\} \quad (6) $$

$$ \vec{U}_t = \vec{U}_n \wedge \vec{U}_g = \left\{ -Z[i] \cdot \vec{X} + X[i] \cdot \vec{Z} \right\} / r[i] \quad (7) $$

**Discrétisation du film**
Une moyenne doit être réalisée sur toutes les valeurs prises par $P_2(X)$ dans un certain domaine du film. Le film est à ce moment là discrétisé en sous domaines, pour chacun desquels il sera calculé la valeur de $\langle P_2(X) \rangle$. Il est important de remarquer que la base locale n'est pas définie à la pointe du cône, car une singularité y est présente. Une "zone morte" est alors définissable, et l'étude sur l'orientation des liaisons y a un traitement spécifique, qui ne sera pas présenté ici 5 (cf. Fig 2). Le film est alors scindé en deux zones (cf. Fig 3) :

1. dans la zone morte 6, où la base locale n'est pas définie.

2. en dehors de la zone morte 7, où la base locale est définie.

<footnote n="2">2. La moyenne d'une variable aléatoire $X$ sera notée $\langle X \rangle$.</footnote>

<footnote n="3">3. Le rayon position vaut : $r[i] = \sqrt{(X[i])^2 + (Z[i])^2}$; $A_i$ et $X_i$ sont dans le même plan radial, avec $X[A_i] = k \cdot X[i]$ et $Z[A_i] = k \cdot Z[i]$ (4). Il vient donc que : $Z[A_i]/X[A_i] = Z[i]/X[i]$.</footnote>

<footnote n="4">4. La base locale est composée de deux vecteurs $\vec{U}_g$ et $\vec{U}_t$, constituants le plan tangent, et d'un troisième vecteur $\vec{U}_n$ normal à la surface. Cette base se calcule avec les coordonnées cylindriques [8].</footnote>

<footnote n="5">5. Les résultats ont montrés que la structure du film y reste globalement amorphe, même si la pointe est enfoncée.</footnote>

<footnote n="6">6. Elle est définie par l'équation $Y[i] < -\tan \theta \cdot r[i] + \delta_p$.</footnote>

<footnote n="7">7. Elle est définie par l'équation $Y[i] \geq -\tan \theta \cdot r[i] + \delta_p$.</footnote>
L'orientation des liaisons sera évaluée par les grandeurs définies précédemment. Le film est découpé en hauteur suivant $Y_{bl} = ||\vec{O_1A_i}||$ et radialement suivant $R_{bl} = ||\vec{A_iX_i}||$ (8)<fnref n="8" />. La résolution suivant $Y_{bl}$ ne se fera pas ($\Delta Y_{bl}$ est pris suffisamment grand pour englober toute la hauteur du film), et celle suivant $R_{bl}$ se fera par pas de $\Delta R_{bl} = 1\sigma$ (cf. Fig 3).

$$R_{bl} = |1 - k| \cdot r[i] \cdot \sqrt{1 + \tan^2 \theta} \quad \text{et} \quad Y_{bl} = |k| \cdot r[i] \cdot \frac{\sqrt{1 + \tan^2 \theta}}{\tan \theta} \quad (8)$$

### 3 Résultats sur l'orientation des liaisons

L'étude présentée ici portera sur l'indentation avec des cônes d'angle $\theta$ valant 30°, 45° et 60°. Le pilotage est en déplacement, avec une vitesse constante de chargement / déchargement ($d\delta/dt = V = 10^{-3}$). Les températures d'expérience (numérique) seront $T = 0.2$ ($< T_g$) et $T = 0.4$ ($\approx T_g$).

A l'observation de la figure 4, il en ressort que la structure du film est bien amorphe en dehors de la zone morte, lorsque la pointe n'est pas présente. Ceci est confirmé par le fait que $\langle P_2(X) \rangle$ fluctue autour de 0.

Lorsque la pointe est enfoncée, l'orientation des liaisons tend à être normale au vecteur $\vec{U}_n$. En effet, $\langle P_2(\cos \alpha^{nor}) \rangle$ tend vers -0.5, ce que montre la figure 5. Les vecteurs de liaison ont donc tendance à être dans le plan tangent au cône, qui est le plan $(\vec{U}_g; \vec{U}_t)$ (figure 5), et ceci toujours en dehors de la zone morte. Par contre, il ne semble pas y avoir d'orientation préférentielle pour l'un ou l'autre des vecteurs du plan tangent au cône. En effet, $\langle P_2(\cos \alpha^{tan}) \rangle$ tend vers 0.3 et $\langle P_2(\cos \alpha^{gen}) \rangle$ tend vers 0.2, ceci révélant une position intermédiaire des vecteurs de liaison entre les deux vecteurs, $\vec{U}_g$ et $\vec{U}_t$, du plan tangent. Par ailleurs, cette zone perturbée par l'enfoncement de la pointe s'étend sur $R_{bl} \approx 2\sigma$ (figure 5), et ceci pour toutes les géométries de cône étudiées, et pour les deux températures, ce que montre la figure 6.

### 4 Conclusion et discussion

Les simulations de DM se révèlent comme étant un outil intéressant dans l'étude de la physique locale du contact. L'étude sur l'orientation des liaisons lors du contact normal, avec un indenteur parfaitement conique, a révélée que les liaisons sont plaquées dans le plan tangent local au cône, et ceci dans une zone de environ $R_{bl} = 2\sigma$, le long du cône. Une étude (non présentée ici) a également révélée que cette orientation se conserve dans le temps. En effet, des observations réalisées sur une durée de $\approx 80\,000\tau$, ont révélées que les vecteurs de liaisons restaient dans le plan tangent : aucun retour vers une structure amorphe n'a été observé, lorsque la pointe est enfoncée. A la pointe du cône (singularité, où la base locale n'est pas définie), une zone morte est définie, et, à l'intérieur de cette dernière, les liaisons conservent une orientation globalement amorphe. Cette analyse est alors un premier pas vers l'étude des contraintes virielles sous la pointe (tenseur viriel) et vers l'étude de la couche cisaillée sous l'indenteur, lors d'un contact tangentiel.

<table><caption>Tableau 1 – Principales unités adimensionnées LJ, en comparaison avec les unités SI.</caption><thead><tr><th>symbole LJ</th><th>Désignation</th><th>symbole SI</th><th>Ordre de grandeur</th></tr></thead><tbody><tr><td>$r^* = (1/\sigma) \cdot r$</td><td>longueur</td><td>$r$ [m]</td><td>$\sigma \approx 0.5 \cdot 10^{-09}$ [m]</td></tr><tr><td>$T^* = (k_b/\varepsilon) \cdot T$</td><td>température</td><td>$T$ [K]</td><td>$\varepsilon/k_b \approx 8.5 \cdot 10^{+02}$ [K]</td></tr><tr><td>$t^* = (1/\tau) \cdot t$</td><td>temps</td><td>$t$ [s]</td><td>$\tau \approx 1.9 \cdot 10^{+00}$ [ps]</td></tr><tr><td>$F^* = (\sigma/\varepsilon) \cdot F$</td><td>force</td><td>$F$ [N]</td><td>$\varepsilon/\sigma \approx 2.3 \cdot 10^{+01}$ [pN]</td></tr></tbody></table>

<footnote n="8">⁸$R_{bl}$ est le rayon dans le plan $(\vec{U}_n, \vec{U}_t)$ à partir du point $A_i$, et, $Y_{bl}$ la norme de $\vec{O_1A_i} \propto \vec{U}_g$.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1 – (a) : schéma montrant l’échelle spatiale servant de cadre à cette étude ; (b) : illustration de la modélisation du contact en simulations de *DM*. La pointe est modélisée avec des particules pour une meilleure clareté, mais, elle est en fait parfaitement lisse.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2 – L’étude de l’orientation des liaisons se fait par rapport aux vecteurs de la base locale au cône : $R_{bl}$ et $Y_{bl}$ repèrent un vecteur liaison. Il est rappelé que $P_2(X) = 1.5 \cdot X^2 - 0.5$, et donc que $\langle P_2(X) \rangle = \langle 1.5 \cdot X^2 - 0.5 \rangle = 1.5 \cdot \langle X^2 \rangle - 0.5$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 4 – Variations de $\langle P_2(X) \rangle (R_{bl})$ autour de 0 dans le cas où la pointe n’est pas présente ( = structure amorphe du film). Pour le calcul des orientations, il est pris $\delta_p = 8\sigma$.</figcaption></figure>

<caption>Figure 5 – Variations de $\langle P_2(X) \rangle (R_{bl})$ dans le cas où la pointe est enfoncée. Pour le calcul des orientations, il est pris $\delta_p = 8\sigma$, ce qui correspond à $\delta = 6.02\sigma$.</caption>
<figure><img src="image_5.png" /></figure>

Figure 6 – Variations de $\langle P_2(\cos \alpha^{\text{nor}}) \rangle (R_{bl})$ dans le cas où la pointe est enfoncée, pour toutes les géométries de cône étudiées, et les deux températures. Pour le calcul des orientations, il est pris $\delta_p = 8\sigma$, ce qui correspond à $\delta = 6.02\sigma$ pour $T = 0.2$, et, à $\delta = 6.50\sigma$ pour $T = 0.4$.

Références

[1] E. Charrault. *Physique des polymères*, Strasbourg : thèse Université Louis Pasteur, 206 p, 2007. http://eprints-scd-ulp.u-strasbg.fr:8080/921/.

[2] B. Schnell. *Physique des polymères*, Strasbourg : thèse Université Louis Pasteur, 206 p, 2006. http://eprints-scd-ulp.u-strasbg.fr:8080/545/.

[3] N. Chen, N. Maeda, M. Tirrell, J. Israelachvili. Adhesion and Friction of Polymer Surfaces : The Effect of Chain Ends. *Macromolecules* **38**, 3491-3503, 2005.

[4] L. Léger, H. Hervet, L. Bureau. Friction mechanisms at polymer-solid interfaces. *C. R. Chimie* **9**, 80-89, 2006.

[5] G. Pätzold, A. Linke, T. Hapke, D. W. Heermann. Computer simulation of nanoindentation into polymer films. *Z. Phys. B 104*, 513-521, 1997.

[6] S. Peter, H. Meyer, J. Baschnagel. Thickness - dependent reduction of the glass-transition temperature in thin polymer films with a free surface. *Journal of Polymer Science* **44**, 2951-2967, 2006.

[7] T. Soddemann, B. Dünweg, K. Kremer K. Dissipative particle dynamics : A useful thermosat for equilibrium and nonequilibrium molecular dynamics simulations. *Physical Review E* **68**, 046702-1 - 046702-8, 2003.

[8] P. Steimann. On boundary potential energies in deformational and configurational mechanics. *J. Mech. Phys. Solids* (2007), doi: 10.1016/j.jmps.2007.07.001.

[9] B. J. Briscoe. *K. Friedrich, Composite materials series - Volume 1 : Friction and wear of polymer composites*, Elsevier, Amsterdam, 1986.

[10] M. P. Allen, D. J. Tildesley. *Computer Simulation of Liquids*, Oxford Science Publications, USA, 1987.

[11] D. Frenkel, B. Smit. *Understanding Molecular Simulation from Algorithms to Applications*, Academic Press, UK, 1996.