0% ont trouvé ce document utile (0 vote)

24 vues120 pages

Modélisation stochastique en mécanique des milieux continus

Ce document présente une modélisation stochastique du comportement mécanique de l'interphase inclusion-matrice dans un nanocomposite à partir de simulations en dynamique moléculaire. Le document décrit la préparation de simulations par dynamique moléculaire du nanocomposite et l'identification des propriétés de l'interphase. Un modèle de champ aléatoire est ensuite proposé pour modéliser les fluctuations des rigidités dans l'interphase, et ses paramètres sont identifiés par homogénéisation numérique.

Transféré par

hh1808199

Nous prenons très au sérieux les droits relatifs au contenu. Si vous pensez qu’il s’agit de votre contenu, signalez une atteinte au droit d’auteur ici.

Formats disponibles

Téléchargez aux formats PDF, TXT ou lisez en ligne sur Scribd

0% ont trouvé ce document utile (0 vote)

24 vues120 pages

Modélisation stochastique en mécanique des milieux continus

Transféré par

hh1808199

Nous prenons très au sérieux les droits relatifs au contenu. Si vous pensez qu’il s’agit de votre contenu, signalez une atteinte au droit d’auteur ici.

Formats disponibles

Téléchargez aux formats PDF, TXT ou lisez en ligne sur Scribd

Modélisation stochastique, en mécanique des milieux

continus, de l’interphase inclusion-matrice à partir de

simulations en dynamique moléculaire
Tien-Thinh Le

To cite this version:

Tien-Thinh Le. Modélisation stochastique, en mécanique des milieux continus, de l’interphase
inclusion-matrice à partir de simulations en dynamique moléculaire. Chemo-informatique. Université
Paris-Est, 2015. Français. �NNT : 2015PESC1172�. �tel-01308194�

HAL Id: tel-01308194

https://theses.hal.science/tel-01308194
Submitted on 27 Apr 2016

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
Université Paris-Est
Année 2015

THÈSE

pour obtenir le grade de

DOCTEUR DE l’UNIVERSITÉ PARIS-EST

Discipline : Mécanique

Présentée et soutenue publiquement par

Tien-Thinh LE

le 21 octobre 2015

Titre de la thèse :

Modélisation stochastique, en mécanique des milieux

continus, de l’interphase inclusion-matrice à partir
de simulations en dynamique moléculaire

Directeur de thèse : Professeur Christian Soize

Jury

Samuel Forest DR – CNRS Mines ParisTech Examinateur

Roger Ghanem PR – University of Southern California Rapporteur
Johann Guilleminot MCF HDR – Université Paris-Est Marne-la-Vallée Encadrant
Djimédo Kondo PR – Université Pierre et Marie Curie Rapporteur
Christian Soize PR – Université Paris-Est Marne-la-Vallée Directeur de thèse
UNIVERSITE PARIS-EST

Résumé
Thèse de Doctorat

Modélisation stochastique, en mécanique des milieux continus, de

l’interphase inclusion-matrice à partir de simulations en dynamique
moléculaire

par Tien-Thinh LE

Dans ce travail, nous nous intéressons à la modélisation stochastique continue et à l’iden-

tification des propriétés élastiques dans la zone d’interphase présente au voisinage des
hétérogénéités dans un nanocomposite prototypique, composé d’une matrice polymère
modèle renforcée par une nanoinclusion de silice. Des simulations par dynamique molécu-
laire (DM) sont tout d’abord conduites afin d’extraire certaines caractéristiques de confor-
mation des chaı̂nes proches de la surface de l’inclusion, ainsi que pour estimer, par des
essais mécaniques virtuels, des réalisations du tenseur apparent associé au domaine de si-
mulation. Sur la base des résultats obtenus, un modèle informationnel de champ aléatoire
est proposé afin de modéliser les fluctuations spatiales du tenseur des rigidités dans l’in-
terphase. Les paramètres du modèle probabiliste sont alors identifiés par la résolution
séquentielle de deux problèmes d’optimisation inverses (l’un déterministe et associé au
modèle moyen, l’autre stochastique et lié aux paramètres de dispersion et de corrélation
spatiale) impliquant une procédure d’homogénéisation numérique. On montre en particu-
lier que la longueur de corrélation dans la direction radiale est du même ordre de grandeur
que l’épaisseur de l’interphase, indiquant ainsi la non-séparation des échelles. Enfin, la
prise en compte, par un modèle de matrices aléatoires, du bruit intrinsèque généré par
les simulations de DM (dans la procédure de calibration) est discutée.
Remerciements
Ce travail de thèse a été réalisé au sein du Laboratoire Modélisation et Simulation
Multi-Echelle (MSME, UMR 8208 CNRS) de l’Université Paris-Est Marne-la-Vallée. Tout
d’abord, j’adresse mes remerciements à mon directeur de thèse Christian Soize et à mon
co-encadrant Johann Guilleminot pour leurs conseils, leur aide scientifique et la confiance
qu’ils m’ont témoignée tout au long de ma thèse. Je remercie chaleureusement mon co-
encadrant de thèse Johann Guilleminot : ses conseils scientifiques, son enthousiasme, sa
disponibilité et son soutien m’ont énormément apporté durant ces trois années de re-
cherche.

J’adresse mes remerciements à Monsieur Samuel Forest qui a accepté de présider le jury
de thèse. J’exprime également mes remerciements à Messieurs Djimédo Kondo et Roger
Ghanem qui m’ont fait l’honneur d’être rapporteurs de ce travail.

Je remercie également Christian Soize et Salah Naili, directeurs successifs du Laboratoire

Modélisation et Simulation Multi-Echelle, pour m’avoir accueilli au sein de leur labora-
toire et offert de bonnes conditions de travail, ainsi qu’un environnement de recherche
appréciable.

Mes remerciements s’adressent à tous les collègues du MSME, notamment ceux avec qui
j’ai eu des échanges scientifiques et amicaux.

Merci à tous mes amis, en France ou au Vietnam, pour les moments enthousiastes après
le travail : moments pendant lesquels j’ai pu me changer les idées et me confier.

Enfin, je remercie ma famille, mon grand-père, mes parents, mon frère, mes cousins et cou-
sines, pour leur soutien et encouragements depuis ma première année d’étude en France.
En particulier, je remercie Quynh Anh, pour son soutien depuis de nombreuses années.

iii
Table des matières

Résumé ii

Remerciements iii

Contenu iv

Liste des figures vii

Liste des tableaux xi

Notations xiii

1 Introduction 1
1 Les nanocomposites : présentation générale . . . . . . . . . . . . . . . . . 1
2 Caractérisation des nanocomposites . . . . . . . . . . . . . . . . . . . . . 3
3 Modélisation du comportement mécanique des nanocomposites . . . . . . 5
4 Positionnement de la recherche et objectifs de la thèse . . . . . . . . . . . 7

2 Simulation du comportement du nanocomposite par DM 11

1 Présentation générale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2 Préparation des configurations initiales pour le nanocomposite . . . . . . 18
3 Caractérisation de la zone d’interphase . . . . . . . . . . . . . . . . . . . 38
4 Essais mécaniques virtuels . . . . . . . . . . . . . . . . . . . . . . . . . . 42
5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

3 Modélisation stochastique du champ des rigidités dans l’interphase 57

1 Hypothèses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
2 Méthodologie de construction . . . . . . . . . . . . . . . . . . . . . . . . 59
3 Construction de la loi marginale d’ordre 1 . . . . . . . . . . . . . . . . . 61
4 Définition et générateur du champ aléatoire . . . . . . . . . . . . . . . . 62
5 Schéma de discrétisation . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
6 Prise en compte du bruit d’échantillonnage dans les simulations par dyna-
mique moléculaire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

v
Contenu vi

4 Identification des paramètres du modèle 71

1 Méthodologie d’identification . . . . . . . . . . . . . . . . . . . . . . . . . 71
2 Résultats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

5 Conclusions et perspectives 89
1 Conclusions générales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
2 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

Bibliographie 93
Table des figures

1.1 Visualisation par TEM des nanorenforts particulaires (diamètre moyen :

18 nm) dans une matrice polymère [28]. . . . . . . . . . . . . . . . . . . . 2
1.2 Evolution du module d’Young mesuré par nanoindentation pour différentes
fractions volumiques et tailles d’inclusion (extrait de [28]). . . . . . . . . 5
1.3 Schéma synoptique des échelles de temps et d’espaces pertinentes en fonc-
tion du cadre de modélisation retenu (adapté d’après [129]). . . . . . . . 6

2.1 Synoptique de la méthode de simulation par DM (d’après [100, 122]). . . 13

2.2 Schématisation des conditions aux limites périodiques en dimension 2 : la
cellule de base, grisée, est répliquée dans les deux directions de l’espace
définies par les vecteurs de la base canonique. . . . . . . . . . . . . . . . 14
2.3 Visualisation d’une chaı̂ne et des sites CH2 (chaque boule représente un
site). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.4 Représentation graphique des potentiels (en kcal/mol) autour des valeurs
d’équilibre : cas des potentiels de rigidité à l’élongation r 7→ Vb (r) (gauche)
et à la flexion θ 7→ Vθ (θ) (droite). . . . . . . . . . . . . . . . . . . . . . . 23
2.5 Représentation graphique des potentiels (en kcal/mol) autour des valeurs
d’équilibre : cas du potentiel de rigidité à la torsion φ 7→ Vφ (φ) (gauche)
et du potentiel de Lennard-Jones 12-6 r 7→ VLJ (r) (droite). . . . . . . . . 23
2.6 Visualisation de deux configurations initiales créées par l’algorithme SARW,
pour npc = 10 chaı̂nes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.7 Relaxation de la configuration initiale, obtenue par l’algorithm SARW, par
dynamique moléculaire (système polymère à 10 chaı̂nes). . . . . . . . . . 25
2.8 Visualisation de configurations relaxées pour npc = 10, npc = 80 et npc =
320 (de gauche à droite). . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.9 Représentation d’une configuration de silice α−quartz simulée pour npc =
10. Les atomes de silicium (Si) et d’oxygène (O) apparaissent en bleu et
rouge, respectivement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.10 De gauche à droite : transformation de la structure cristalline vers la struc-
ture amorphe à l’aide d’une relaxation par dynamique moléculaire (coupe
bi-dimensionnelle). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.11 Graph de la densité de probabilité pour la longueur de liaison Si−O. Trait
continu : modèle complet. Trait discontinu : modèle équivalent. . . . . . . 31
2.12 Graph des densités de probabilité pour les angles de valence Si−O−Si
(trait fin) et O−Si−O (trait épais). Trait continu : modèle complet. Trait
discontinu : modèle équivalent. . . . . . . . . . . . . . . . . . . . . . . . . 31

vii
Liste des ﬁgures viii

2.13 Principe de fabrication virtuelle du nanocomposite à partir d’une inclusion

de silice et de la matrice polymère à l’équilibre (reproduction d’après [18]). 33
2.14 Configuration initiale de l’inclusion de silice après traitement des atomes
de surface : : atome d’oxygène dans le cœur de l’inclusion ; : atome de
silicium dans le cœur de l’inclusion ; ◦ : atome de silicium dans la couche
externe de l’inclusion ; * : atome d’oxygène dans la couche externe de l’in-
clusion connecté avec deux atomes de silicium (contenus dans la couche ex-
terne) ; + : atome d’oxygène dans la couche externe de l’inclusion connecté
avec un atome de silicium (contenu dans la couche externe) . . . . . . . . 33
2.15 Visualisation de différentes configurations initiales des nanoinclusions :
Rp = 1.5 nm (gauche), Rp = 3 nm (centre) et Rp = 4.8 nm (droite). . . . 34
2.16 Création progressive de la porosité (de gauche à droite) par l’introduction
d’une force répulsive au centre de la boı̂te (cas npc = 10). . . . . . . . . . 35
2.17 Visualisation de deux configurations relaxées du nanocomposite pour npc =
10, avec Rp = 1.5 (à gauche) et Rp = 3 (à droite). Les fractions volumiques
sont égales à 4.8% et 28.5%, respectivement. . . . . . . . . . . . . . . . . 35
2.18 Evolution de la distance entre le centre d’inertie de l’inclusion et le centre de
la boı̂te de simulation, pour npc ∈ {10, 80} et pour une fraction volumique
de 4.8%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.19 Evolution temporelle du déplacement quadratique moyen (avec tM SD =
500 ps) pour un nanocomposite à 80 chaı̂nes, dans l’inclusion et dans la
phase polymère. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.20 Evolution temporelle du déplacement quadratique moyen (avec tM SD =
500 ps) : comparaison entre le polymère pur à 80 chaı̂nes et différents
nanocomposites (fraction volumique : 4.8%). . . . . . . . . . . . . . . . . 39
2.21 Graph de la fonction (r−Rp ) → ρn (r)/ρp pour les systèmes nanocomposites
définis par npc ∈ {10, 80, 320} (fraction volumique : 4.8%). . . . . . . . . 40
2.22 Moyenne des carrés de déplacements du polymère dans les coquilles concen-
triques du nanocomposite avec 80 chaı̂nes à 100 K, moyenné sur un inter-
valle de temps de 5 ps qui est répété 100 fois. . . . . . . . . . . . . . . . 41
2.23 Graph de la fonction (r − Rp ) → ρn (r)/ρp pour des systèmes exhibant des
longueurs de chaı̂nes différentes (Rp = 1.5 nm). . . . . . . . . . . . . . . 41
2.24 Graph de la fonction (r − Rp ) → ρn (r)/ρp pour des systèmes générés à
partir de densités différentes. . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.25 Evolution du volume spécifique en fonction de la température (cycle de
décharge) pour différents systèmes (polymère pur et nanocomposite, npc =
80). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
2.26 Courbes de déformation-contrainte pour les essais de traction sur la silice
amorphe (T = 100 K). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.27 Courbes de déformation-contrainte pour les essais de cisaillement sur la
silice amorphe (T = 100 K). . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.28 Simulation d’un essai de compression hydrostatique sur la silice amorphe
(T = 100 K). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
2.29 Configurations initiale (gauche) et déformée à 10% (droite) du nanocompo-
site, pour un essai de traction uniaxiale. Les atomes Si et O sont représentés
en bleu et en rouge, respectivement. . . . . . . . . . . . . . . . . . . . . . 49
Liste des figures ix

2.30 Configurations initiale (gauche) et déformée à 5% (droite) du nanocompo-

site, pour un essai de cisaillement. Les atomes Si et O sont représentés en
bleu et en rouge, respectivement. . . . . . . . . . . . . . . . . . . . . . . 49
2.31 Courbes de contrainte-déformation pour le polymère pur et pour le nano-
composite (à 100 K) : essai de traction uniaxiale. . . . . . . . . . . . . . 50
2.32 Courbes de contrainte-déformation pour le polymère pur et pour le nano-
composite (à 100 K) : essai de cisaillement. . . . . . . . . . . . . . . . . . 50
2.33 Réponse du polymère pur et du nanocomposite pour un essai de compres-
sion hydrostatique à 100 K. . . . . . . . . . . . . . . . . . . . . . . . . . 51

3.1 Graphique de la fonction de corrélation normalisée (xs , y s ) 7→ ρ (xs , y s ),

avec xs = (1.5, 0, π/2) et y s = (r, θ, π/2) pour 1.5 6 r 6 3.5, 0 6 θ 6 2π,
Lr = 0.2 et Lθ = 0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.2 Graphique de la fonction de corrélation normalisée (xs , y s ) 7→ ρ (xs , y s ),
avec xs = (1.5, 0, π/2) et y s = (r, θ, π/2) pour 1.5 6 r 6 3.5 et 0 6 θ 6 2π.
Figure de gauche : Lr = 2, Lθ = 0.5. Figure de droite : Lr = 0.2, Lθ = 2. 64

4.1 Représentation graphique des différents domaines considérés (inclusion de

silice, zone d’interphase et polymère pur). . . . . . . . . . . . . . . . . . 73
4.2 Visualisation du maillage pour le modèle continu (l’inclusion apparaı̂t en
magenta, l’interphase en blanc et la matrice de polymère pur en cyan). . 76
4.3 Convergence de l’algorithme d’optimisation pour le problème défini par les
Eqs. (4.10) et (4.11) : graphique de la fonction k 7→ min L1 (k). . . . . . . 79
4.4 Graphique de la fonction Niter 7→ Conv (Niter ) (% = 9, ∆t = 0.01). . . . . . 80
4.5 Graphique de la fonction τ11 7→ Decorr(τ11 ). . . . . . . . . . . . . . . . . 81
4.6 Graphique 4D de la fonction de vraisemblance h2 7→ L2 (h2 ) pour différentes
combinaisons de δ[N ] , Lr et La . La valeur de la fonction est proportionnelle
à la taille du marqueur, ainsi qu’à l’intensité de la couleur (les valeurs in-
finies de L2 sont remplacées par la valeur arbitraire -100, afin de visualiser
l’ensemble des résultats). . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4.7 Configuration instantanée des chaı̂nes polymères piégées dans des coquilles
sphériques d’une épaisseur de 5 Å dans la région de l’interphase. Sur la
figure de gauche (resp. droite), la coquille correspond à la coquille la plus
proche (resp. éloignée) de la surface de l’inclusion. . . . . . . . . . . . . . 83
4.8 Densité de probabilité de la composante (1, 1) du tenseur apparent (trait
noir). Les réalisations associées, obtenues par les simulations de DM, sont
représentées par les étoiles rouges. . . . . . . . . . . . . . . . . . . . . . . 83
4.9 Densité de probabilité de la composante (1, 2) du tenseur apparent (trait
noir). Les réalisations associées, obtenues par les simulations de DM, sont
représentées par les étoiles rouges. . . . . . . . . . . . . . . . . . . . . . . 84
4.10 Densité de probabilité de la composante (4, 4) du tenseur apparent (trait
noir). Les réalisations associées, obtenues par les simulations de DM, sont
représentées par les étoiles rouges. . . . . . . . . . . . . . . . . . . . . . . 84
int
4.11 Visualisation d’une réalisation de l’élasticité du champ aléatoire {C11 (xs ), xs ∈
ext
∂DI } (en GPa et en coordonnées sphériques) sur la surface extérieure de
la région de l’interphase (avec hyperparamètres calibrées). . . . . . . . . 85
Liste des figures x

4.12 Graph de la fonction de vraisemblance δ[N ] , δbruit 7→ L3 δ[N ] , δbruit .
Le résultat sans la modélisation du bruit correspond à la courbe δ[N ] 7→
L3 δ[N ] , 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
Liste des tableaux

1.1 Consommation industrielle de nanoparticules en Europe en 2007 [75]. . . 2

2.1 Description des systèmes considérés pour une fraction volumique en ren-
forts de 4.8%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2 Description des systèmes considérés pour une fraction volumique en ren-
forts de 28.5%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.3 Liste des paramètres des potentiels utilisés pour la modélisation du po-
lymère par dynamique moléculaire. . . . . . . . . . . . . . . . . . . . . . 22
2.4 Quelques grandeurs caractéristiques des systèmes estimées à l’équilibre. . 26
2.5 Paramètres du potentiel de Buckingham. . . . . . . . . . . . . . . . . . . 27
2.6 Paramètres géométriques définissant la cellule de base pour la silice cris-
talline de type α−quartz. . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.7 Paramètres des potentiels remplacés pour la silice. . . . . . . . . . . . . . 30
2.8 Comparaison des valeurs moyennes pour certaines propriétés structurelles
et pour la densité volumique. . . . . . . . . . . . . . . . . . . . . . . . . 30
2.9 Répartition du nombre d’atomes dans les nanoinclusions de silice (modèle
équivalent), pour Rp ∈ {1.5, 3, 4.8, 6} (nm). . . . . . . . . . . . . . . . . . 34
2.10 Paramètres des nanocomposites à l’équilibre, pour une fraction volumique
en renforts de 4.8%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.11 Paramètres des nanocomposites à l’équilibre, pour une fraction volumique
en renforts de 28.5%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.12 Estimation de Tg pour différents systèmes (polymère pur et nanocomposite). 44
2.13 Estimation des modules élastiques de la silice (à T = 100 K, en GPa) et
comparaison avec les données de la littérature (les travaux numériques et
expérimentaux sont signalés par les symboles (*) et (**) respectivement). 48
2.14 Estimation des modules élastiques (en GPa) du polymère pure avec 10
chaı̂nes à 100 K et leurs comparaisons. . . . . . . . . . . . . . . . . . . . 52
2.15 Estimation des modules élastiques du nanocomposite avec 10 chaı̂nes à 100
K et leurs comparaisons. . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

xi
Notations

Ensembles :
Mn (R) Ensemble des matrices carrées réelles n × n
MD
n (R) Ensemble des matrices carrées réelles n × n diagonales
MSn (R) Ensemble des matrices carrées symétriques réellesn × n
M+
n (R) Ensemble des matrices carrées symétriques définies-positives réelles n × n
Variables :
x Vecteur déterministe
X Vecteur aléatoire
[X] Matrice (ou tenseur d’ordre 2) déterministe
[X] Matrice (ou tenseur d’ordre 2) aléatoire
JXK Tenseur d’ordre 4 (déterministe)
Normes :
Pn
kxk Norme euclidienne dans Rn , kxk2 = x2i i=1

k[X]kF Norme de Frobenius, k[X]k2F = tr [X] [X] T

Opérateurs :
E Espérance mathématique
T
Transposition
tr ([X]) Trace matricielle
Divers :
c Constante de normalisation 1

1. Lorsqu’aucune ambiguı̈té n’est possible, la valeur de c est susceptible de changer de ligne à ligne
sans indication spécifique.

xiii
Chapitre 1

Introduction

Ce chapitre introductif est destiné à une présentation succincte des matériaux nanocom-
posites, de leur caractérisation expérimentale à leur modélisation dans différents cadres
d’étude. Une description générale de ces matériaux est tout d’abord présentée à la Sec. 1.
La caractérisation des spécificités exhibées par de telles microstructures, et notamment le
concept d’effet nano, est ensuite brièvement exposée à la Sec. 2. La Sec. 3 est dédiée aux
méthodes de modélisation analytiques et numériques pour le comportement mécanique
de milieux nanorenforcés. Les objectifs de ce travail et le positionnement de la recherche
sont enfin dégagés, sur la base des éléments précités, à la Sec. 4.

1 Les nanocomposites : présentation générale

Les nanocomposites sont des matériaux hétérogènes, en général de type matrice-inclusions,

dans lesquels les renforts sont de taille nanométrique – avec une taille caractéristique typi-
quement inférieure à 100 nm [96]. Contrairement aux matériaux composites renforcés par
des inclusions micro- et millimétriques, ces matériaux exhibent des nouveaux phénomènes
multi-physiques qui sont liés aux interactions (par exemple, entre la matière confinée et
des empilements de renforts plaquettaires) aux plus petites échelles. Ces interactions in-
duisent des modifications importantes des propriétés physiques locales, notamment de la
phase matrice (transition/modification de phase, conformation spécifique, modification
du degré de cristallinité pour une matrice polymère organique, etc.), qui sont à l’origine
d’une amélioration, en générale significative, de certaines propriétés (mécaniques, op-
tiques, électriques, etc.) macroscopiques. D’un point de vue mécanique (et de façon non

1
Introduction 2

exclusive), de nombreuses études expérimentales et numériques (voir les références ci-

dessous) ont montré que ce saut bénéfique de propriétés 1 est d’autant plus prononcé que
la taille des inclusions est petite (pour une fraction volumique en renforts suffisamment
importante) : on parle alors d’effet nano. La caractérisation et la modélisation de ce der-
nier, tout autant que son utilisation pour l’optimisation de certaines propriétés d’intérêt,
a motivé un très grand nombre de travaux à la frontière entre plusieurs communautés
scientifiques (chimistes, physico-chimistes, mécaniciens, etc.), ainsi qu’un intérêt indus-
triel croissant. En fonction de l’application considérée, la matrice est de type polymère
organique [26], métallique [8] ou céramique [120], tandis que les renforts sont de nature
particulaire, fibrillaire ou plaquettaire (voir [75] pour une synthèse dédiée aux nanocom-
posites à matrice polymère). Dans le cas des matrices polymères qui nous intéressera
dans ce travail, les inclusions particulaires les plus utilisés sont le noir de carbone et les
oxydes métalliques. A titre d’illustration, le Tab. 1.1 présente la consommation industrielle
européeenne de ce type de renforts pour l’année 2007, et la Fig. 1.1 présente une visuali-
sation par Microscopie Electronique en Transmission (TEM) de particules nanométriques
de silice dispersées dans une matrice de polysiloxane .

Type de nanoparticule Consommation (en tonnes)

Noir de carbone 2 000 000
Carbonate de calcium naturel 1 500 000
Hydroxyde d’aluminium 250 000
Silice 200 000

Table 1.1 – Consommation industrielle de nanoparticules en Europe en 2007 [75].

Figure 1.1 – Visualisation par TEM des nanorenforts particulaires (diamètre moyen :
18 nm) dans une matrice polymère [28].

Les nanocomposites à base de matrice polymère sont de fait très utilisés dans le cadre de
l’industrie, en raison de l’efficacité des procédés d’élaboration associés et du faible coût
1. On notera ici qu’il s’agit de considérations tout à fait générales et qu’a contrario, il existe des
matériaux pour lesquels l’ajout de renforts nanométriques entraı̂ne une diminution de certaines propriétés
mécaniques élastiques ou à rupture, par exemple (voir e.g. [83]).
Introduction 3

de ce type de matériau [45, 54]. La dispersion des particules dans la matrice peut s’ef-
fectuer par mélange direct du polymère soluble dans l’eau et des nanoparticules, par po-
lymérisation in situ en présence des nanoparticules ou encore par mélange de ces dernières
dans le polymère fondu [104]. Les propriétés finales du matériau dépendent de nombreux
paramètres, parmi lesquels (et de façon prépondérante) la qualité de la dispersion des
hétérogénéités dans la matrice, ou bien encore la qualité de l’adhésion entre les consti-
tuants (qui peut être améliorée par l’ajout d’agents de couplage chimiques [89]). Dans
ce cadre, on notera que l’optimisation des procédés d’élaboration est un domaine de re-
cherche très actif, en particulier pour favoriser [3] ou éviter la présence d’agrégats (qui
peuvent entraı̂ner, par l’intermédiaire de concentrations de contraintes, une diminution
des propriétés à rupture).

2 Caractérisation des nanocomposites

2.1 Caractérisation de l’interphase entre la matrice et les ren-

forts

La prédiction des propriétés physiques exhibées par les nanocomposites est un problème
central. Pour de tels systèmes, les interactions aux échelles fines (entre les charges et
la matrice, ou entre les charges elles-mêmes) ne sont plus négligeables et doivent être
interprétées et modélisées dans un cadre multi-échelle. Une difficulté majeure ici réside
dans l’interaction entre différents phénomènes et caractéristiques, comme la nature de
l’interface ou le transfert de charge [36, 75]. L’un des aspects les plus documentés dans
la littérature, tant d’un point de vue expérimental que d’un point de vue numérique, est
l’existence d’une interphase (c’est-à-dire d’un troisième milieu de propriétés inconnues a
priori, dans le cas d’un matériau matrice-inclusions) entourant les hétérogénéités.

Dans un cadre expérimental, cette présence d’une zone perturbée autour des inclusions est
mise en évidence par des mesures de Résonance Magnétique Nucléaire (RMN), utilisées
pour étudier la mobilité locale des atomes ; voir [12, 13, 34, 47, 48, 92] pour le cas de
renforts en silice, et [62] pour des résultats similaires pour des inclusions de type de noir
de carbone. Ces travaux démontrent une diminution notable de la mobilité des chaı̂nes
de la matrice au voisinage des particules, ainsi qu’une modification de la densité locale.
L’épaisseur de cette zone intermédiaire, appelée interphase, est d’autant plus faible que la
température augmente, mais demeure contante (à la précision des mesures près) lorsque
Introduction 4

la taille des renforts est modifiée (pour une même température). Par conséquent, une
diminution de la taille des renforts induit une augmentation de la surface spécifique de
ceux-ci. Cet effet, appelé effet de surface, est couramment invoqué dans le développement
de modèles multi-échelle pour expliquer l’effet nano précité. Notons enfin que l’épaisseur
de la zone d’interphase dépend des interactions entre les renforts et la matrice. Ce point
a été spécifiquement mis en évidence dans [89], où des caractérisations par RMN sont
conduites (à 350 K) pour plusieurs systèmes composés d’une même matrice, mais avec
des renforts de tailles variables et l’intégration de deux agents de couplage différents
(l’un favorisant la création de liaisons covalentes, l’autre celle de liaisons faibles). Pour
l’ensemble des configurations testées, on observe des résultats qualitatifs similaires, mais
une épaisseur de la zone perturbée très différente en fonction de l’agent de couplage
considéré.

Ces observations expérimentales sur l’interphase dans le nanocomposite sont par ailleurs
bien corrélées par des simulations moléculaires. Ces dernières ont par ailleurs mis en
évidence une orientation préférentielle des segments de chaı̂ne, situées dans le plan tan-
gentiel au vecteur normal à la surface des nanoparticules ; voir par exemple [9, 17, 18,
30, 38, 77, 78, 83, 119] pour des simulations de dynamique moléculaire, ainsi que [125–
128, 135] pour l’utilisation de la méthode de Monte Carlo.

2.2 Caractérisation du comportement mécanique

De façon générale, l’ajout des nanorenforts permet une amélioration de plusieurs pro-
priétés physiques, parmi lesquelles les propriétés mécaniques. Ci-dessous, on restreint
essentiellement l’exposé à l’étude des propriétés élastiques linéaires qui seront utilisées
dans la suite de ce travail (voir e.g. [50] pour une analyse, par simulations de dynamique
moléculaire, du comportement non linéaire). Une illustration du caractère renforçant est
présentée dans [28], où un réseau de polysiloxane se voit renforcé par des nanoparticules
de silice. La caractérisation expérimentale est accomplie par des essais de traction uni-
axiale, pour trois valeurs (respectivement égales à 5, 10 et 15%) de la fraction volumique
en silice. La Fig. 1.2 montre l’évolution du module d’Young mesuré par nanoindenta-
tion, pour les différentes valeurs de fraction volumique (les diamètres des particules sont
également reportés ; le module d’Young de la matrice pure est d’environ 3 GPa).

On observe bien un effet de renforcement croissant lorsque la fraction volumique aug-

mente, et ce quel que soit la taille de l’inclusion. Par ailleurs, l’effet de surface est claire-
ment visualisé en comparant les résultats pour les particules de diamètres respectifs 15 et
Introduction 5

Figure 1.2 – Evolution du module d’Young mesuré par nanoindentation pour

différentes fractions volumiques et tailles d’inclusion (extrait de [28]).

30 nm, et s’avère d’autant plus prononcé que la fraction volumique augmente. Cet effet
n’est toutefois pas mis en évidence avec les particules de diamètre 60 nm, pour lesquelles
les renforts tendent à former des agrégats à l’origine de l’augmentation très significative
du module (voir [23], par exemple). D’autres résultats expérimentaux sont disponibles
dans [23, 65, 133] et les références citées.

3 Modélisation du comportement mécanique des na-

nocomposites

3.1 Modélisations moléculaires et couplages numériques

Les modélisations et simulations atomistiques constituent des outils essentiels à la compré-

hension des phénomènes physiques aux échelles nanométriques, ainsi qu’à la prédiction
de propriétés d’intérêt – voir de façon non exhaustive [2, 67, 83, 99, 129–131] en ce
qui concerne l’utilisation de simulations par dynamique moléculaire. Outre le choix et
l’identification des modèles et paramètres mis en jeu, une difficulté majeure ici réside dans
le fait que les échelles de temps et d’espace pertinentes dans ce cadre demeurent hors de
portée des architectures de calcul actuelles, et bien inférieures à celles traditionnellement
considérées en mécanique des milieux continus. Ces cadres théoriques de modélisation et
de simulation, ainsi que les différentes échelles associées, sont schématisés sur la Fig. 1.3.
Introduction 6

Méthodes de
modélisation

Chimie Modélisation Mécanique

Numérique Multi-échelle Numérique

Outils de
modélisation

Mécanique
Nanomécanique Micromécanique Macromécanique
Quantique

Structure moléculaire
Milieu continu
discrète

Echelle de
longueur (m)

10-10 10 -8 10 -6 10 0

Echelle de
temps (sec)

10-15 10-10 10 x

Figure 1.3 – Schéma synoptique des échelles de temps et d’espaces pertinentes en

fonction du cadre de modélisation retenu (adapté d’après [129]).

Afin de contourner, dans une certaine mesure, les restrictions induites par ces différences
d’échelle, de nombreuses méthodes ont été proposées afin de coupler, de façon concourante
ou séquentielle, des simulations de dynamique moléculaire et des modélisations multi-
échelles formulées dans le cadre de la mécanique des milieux continus. L’objectif principal
ici est d’extraire l’information physique (e.g. mécanique et/ou électrique) essentielle aux
plus petites échelles, par l’intermédiaire d’une description moléculaire, puis de transférer
cette information par une transformation ad hoc dans un cadre de modélisation continue.
Notons que ce passage peut par ailleurs s’effectuer de la représentation continue vers la
description moléculaire (par la prescription de conditions aux limites par exemple) et de
façon globale ou locale (i.e. sur un sous-domaine du domaine de modélisation continue, où
la physique doit être décrite de façon moléculaire), lorsque les deux types de description
sont itérativement couplés. On pourra se référer à [2, 67, 83, 86, 87] pour des cadres
méthodologiques permettant d’extraire (de façon directe ou par une approche inverse)
des propriétés mécaniques à partir de simulations de dynamique moléculaire, ainsi qu’à
[22, 24, 93, 94, 97, 98, 102, 103, 106] pour des stratégies de couplage numérique.
Introduction 7

3.2 Modélisations en mécanique des milieux continus

Du point de vue de la mécanique des milieux continus, deux types d’approche ont
été développés afin d’estimer les propriétés effectives de matériaux nanorenforcés et de
prendre en compte l’effet de surface introduit précédemment. Une première classe de
méthodes repose sur l’intégration de modèles d’interfaces dans des formulations mi-
cromécaniques classiques (e.g. de type Eshelby). La formulation repose en général sur
l’introduction d’une élasticité surfacique de type Gurtin-Murdoch [42, 43] : on obtient
alors une dépendance explicite des estimations des modules homogénéisés par rapport
aux dimensions caractéristiques des inclusions nanométriques (comme le rayon dans le cas
d’une charge sphérique). De très nombreuses applications de ce type d’approche peuvent
être trouvées dans la littérature ; voir e.g. [10, 16, 29, 51, 59, 61, 141, 142] et les références
citées. Notons que les paramètres de la loi de comportement surfacique sont inconnus a
priori et ne définissent pas un tenseur nécessairement défini-positif. Leurs valeurs peuvent
être estimées par des calculs directs ou en résolvant un problème inverse à partir de si-
mulations de dynamique moléculaire (un exemple est fourni dans [73]).

Dans un second type d’approche, l’interphase est considérée comme une phase addition-
nelle de volume fini et exhibant un tenseur d’élasticité défini-positif. Dans ce cas, les
propriétés équivalentes peuvent être déterminées par des schémas micromécaniques clas-
siques, tels que le schéma auto-cohérent généralisé [67, 83, 86, 87], ou par l’intermédiaire
d’une méthode numérique adaptée [91].

4 Positionnement de la recherche et objectifs de la

thèse

Malgré le caractère stochastique des phénomènes physiques se produisant lors de l’élabora-

tion des nanocomposites (des formes précises des inclusions aux interactions locales
de différentes natures), les différentes contributions discutées jusqu’à présent ont été
développées dans des cadres déterministes, en supposant notamment que les propriétés
dans la zone de l’interphase sont déterministes, constantes ou à gradients de propriétés.
Très récemment, plusieurs études concernant la quantification des incertitudes dans la
modélisation de matériaux nanorenforcés ont été proposées. Des analyses de sensibilité
reposant sur la propagation d’incertitudes paramétriques (portant sur les paramètres
Introduction 8

des champs de force dans les simulations de dynamique moléculaire, ou sur d’autres pa-
ramètres tels que la fraction volumique en renforts dans la formulation continue), ont
été proposées dans [137–139]. Une modélisation du bruit d’échantillonnage et d’incerti-
tudes sur les paramètres des champs de force a été également définie dans [102, 103], en
invoquant une approche d’identification Bayésienne et des représentations sur les chaos
polynomiaux [21, 39] – on trouvera une approche similaire dans [106], dans le cadre
d’un couplage concourant (concernant la transmission de conditions aux limites) entre
un modèle atomistique et un modèle continu de mécanique des fluides. Une méthode in-
verse d’identification des propriétés élastiques déterministes de l’interphase, supposée iso-
trope, est proposée dans [24]. L’approche consiste à imposer l’équivalence des propriétés
homogénéisées (au sens usuel du terme, c’est-à-dire à l’échelle du volume élémentaire
représentatif) obtenues par un calcul de dynamique moléculaire et par une formulation
continue d’homogénéisation (analytique ou numérique). Dans cette dernière, les parti-
cules sont supposées de rayon aléatoire (ce dernier étant distribué selon une loi Beta),
et distribuées de façon uniforme dans le domaine considéré. Enfin, l’analyse fiabiliste de
structures à base de matrices polymères renforcés par des nanotubes de carbone a été
étudiée dans [40], par l’intermédiaire de la propagation d’incertitudes paramétriques sur
une suite échelles pertinentes.

Ces premiers travaux constituent donc des prémices à la prise en compte de modélisations
probabilistes à l’échelle nanoscopique, ainsi qu’à leur propagation jusqu’aux échelles
d’intérêt. En particulier, on note que l’ensemble des approches proposées repose sur la
considération de lois de probabilité supposées a priori, et que l’élasticité dans la zone
d’interphase est toujours considérée comme constante d’un point de vue spatial (auquel
cas elle peut donc être modélisée par des variables aléatoires) et isotrope. L’objectif de ce
travail est donc d’apporter une contribution à ce champ de recherche, en proposant une
modélisation stochastique plus fine de l’élasticité dans l’interphase, modélisée comme un
champ aléatoire non Gaussien à valeurs tensorielles, ainsi qu’une méthodologie d’identi-
fication des hyperparamètres des modèles construits. Ces deux aspects sont abordés sur
la base de simulations de dynamique moléculaire, qui apportent un éclairage physique
sur les contraintes de modélisation mathématiques. Plus spécifiquement, ces simulations
permettent d’une part d’inférer sur certains propriétés fondamentales du champ, comme
la symétrie matérielle locale exhibée ou la structure de corrélation sous-jacente, et d’autre
part de construire une base de données de référence (sur les propriétés mécaniques macro-
scopiques, au sens de la physique statistique) dédiée à la calibration des représentations
probabilistes.
Introduction 9

Le manuscrit est bâti autour de trois chapitres. Dans le Chap. 2, on s’intéresse à la

modélisation d’un nanocomposite modèle par l’intermédiaire de simulations de dyna-
mique moléculaire. Après une brève introduction situant quelques points clés (théoriques
et pratiques) de l’approche, une présentation du système considéré, constitué d’une ma-
trice polymère renforcée par une nanoinclusion de silice, est proposée. En particulier,
la définition des champs de force est détaillée, ainsi que la procédure de génération des
configurations initiales. La caractérisation de la zone d’interphase est ensuite abordée, au
travers de l’étude de la mobilité des atomes du polymère proches de la surface de l’inclu-
sion, ainsi que de la conformation locale des segments de chaı̂ne. Les essais mécaniques
virtuels sont enfin présentés et permettent, outre la construction de la base de donnée
susmentionnée, de mettre en évidence le caractère renforçant de l’hétérogénéité.

Le Chap. 3 est dédié à la construction du modèle probabiliste pour le champ aléatoire

(non Gaussien) des rigidités dans la zone d’interphase. Celle-ci est accomplie dans le cadre
de la Théorie de l’Information et repose sur une transformation non linéaire de champs
Gaussiens sous-jacents. D’après les résultats obtenus au Chap. 2, le champ aléatoire exhibe
une symétrie isotrope transverse radiale, définie dans un système de repérage sphérique, et
voit sa structure de corrélation dépendre d’un ensemble de trois paramètres scalaires liés
aux longueurs de corrélation dans la direction radiale et les deux directions orthoradiales.
Un schéma de génération, basée sur la résolution d’une famille d’équations différentielles
stochastique d’Itô indexée en espace, est ensuite détaillé. La prise en compte du bruit
d’échantillonnage dans la représentation probabiliste est ensuite discutée et formalisée à
l’aide d’un modèle de matrices aléatoires.

L’identification des modèles stochastiques proposés fait l’objet du dernier chapitre. Les si-
mulations de dynamique moléculaire ne permettant pas une calibration directe du modèle
de champ aléatoire, un problème statistique inverse est formulé et repose sur l’équivalence
des propriétés apparentes (et donc, aléatoires) du nanocomposite estimées soit par les cal-
culs de dynamique moléculaire, soit par une méthode d’homogénéisation numérique as-
sociée à la description en mécanique des milieux continus. La méthodologie est par la suite
appliquée aux modèles probabilistes (associés au champ aléatoire des rigidités et aux ma-
trices aléatoires modélisant le bruit induit par les simulations de dynamique moléculaire)
développés au Chap. 2. On montre en particulier que les paramètres liés aux longueurs
de corrélation du champ sont du même ordre de grandeur que certains paramètres phy-
siques, tels que l’épaisseur de l’interphase, et que le bruit d’échantillonnage doit être pris
en compte afin d’estimer de façon robuste le niveau des fluctuations statistiques dans
celle-ci.
Introduction 10

Une conclusion et quelques perspectives générales viennent enfin conclure l’exposé.

Chapitre 2

Simulation du comportement du
nanocomposite par DM

Dans ce chapitre, nous nous intéressons à la modélisation du comportement mécanique

d’un nanocomposite modèle à l’aide de simulations par Dynamique Moléculaire (DM). Les
résultats issus de cette analyse seront utilisés aux chapitres 3 et 4 afin de procéder, respec-
tivement, à la construction et à l’identification inverse du champ aléatoire représentant
l’élasticité dans la zone d’interphase. Un bref rappel théorique associé à la méthode de si-
mulation est présenté à la Sec. 1. La définition des interactions entre les différentes unités,
ainsi que la méthodologie de création des configurations initiales, sont ensuite détaillées
à la Sec. 2. La caractérisation de la zone d’interphase, au travers de mesures de mobilité
et d’orientation, est abordée à la Sec. 3. Les essais mécaniques virtuels permettant l’es-
timation des propriétés apparentes du nanocomposite sont enfin définis et discutés à la
Sec. 4.

1 Présentation générale

1.1 Introduction

Les simulations par dynamique moléculaire ont pour objectif d’estimer des propriétés ma-
croscopiques et/ou d’identifier des mécanismes physiques (de déformation, de conforma-
tion, etc.) associés à un système décrit par un ensemble de particules (en général identifié

11
Simulations par dynamique moléculaire 12

à un ensemble d’atomes) Newtoniennes, satisfaisant les équations classiques de mouve-

ment. Les interactions sont alors régies par des champs de force, dont la détermination fait
toujours l’objet de discussions et de recherches actives (on notera par ailleurs l’émergence
d’approches alternatives, telles que la Dynamique Moléculaire ab initio [69] dans laquelle
les forces d’interaction sont estimées au vol à partir de calculs quantiques) – a contrario
des méthodes dites de Monte Carlo, dans lesquelles les évolutions de configuration sont
purement aléatoires et basées sur des critères d’acceptation/rejet probabilistes [72]. D’un
point de vue numérique, les premières études de ce type furent réalisées dans les années
1950, et concernaient l’analyse d’ensemble de sphères dures [4, 5]. De nombreux ouvrages
de référence sont disponibles sur ce sujet, et le lecteur intéressé pourra consulter, de façon
non-exhaustive, les références [6, 35, 64, 100, 122, 124].

1.2 Description schématique

De façon générale, la méthode de simulation par DM est basée sur quatre étapes schémati-
sées sur la Fig. 2.1. Considérons un système constitué de N particules évoluant dans un
espace de dimension 3. Les vecteurs position, vitesse et accélération de la particule i,
1 6 i 6 N , sont notés r i = (r1i , r2i , r3i ), v i = (v1i , v2i , v3i ) et ai = (ai1 , ai2 , ai3 ) respectivement.
Dans ce qui suit, pi = mi v i (avec mi la masse de la particule en question), 1 6 i 6 N ,
désigne la quantité de mouvement. On introduit par ailleurs les notations suivantes :

r := (r 1 , . . . , r N ) , p := (p1 , . . . , pN ) , a := (a1 , . . . , aN ) . (2.1)

En notant V l’énergie potentielle du système, fonction du vecteur position r, on obtient

l’expression de la force f i s’exerçant sur la i-ème particule :

∂V(r)
fi = − , 16i6N . (2.2)
∂r i

L’Hamiltonien H du système s’écrit par ailleurs

N
X kpi k2
H(r, p) = V(r) + . (2.3)
i=1
2mi

Dans un cadre d’analyse classique, des conditions aux limites périodiques sont imposées
sur les bords du domaine (borné) de simulation. Ces conditions, obtenues par une réplica-
tion virtuelle du domaine dans les trois directions de l’espace (voir la schématisation de
la Fig. 2.2 ci-dessous), se traduisent en terme de déplacements (une particule sortant du
Simulations par dynamique moléculaire 13

Etape 1: Initialisation
positions initiales
vitesses initiales
condition aux bords

Etape 2: Interaction
calcul de l'énergie potentielle
calcul des forces
calcul des accélérations

Etape 3: Intégration
intégration des équations de mouvement
nouvelles positions, vitesses
contrôle de température, pression

Non
Equilibre?

Oui

Etape 4: Collecte des résultats

répéter les étapes 2 et 3
calcul des moyennes
analyse des résultats

Figure 2.1 – Synoptique de la méthode de simulation par DM (d’après [100, 122]).

domaine voit son image la plus proche entrer dans celui-ci) et d’interactions, et permettent
de reproduire le comportement du milieu infini. Dans ce qui suit, on décrit très brièvement
chacune des étapes de la simulation par DM.

1.2.1 Etape 1 : initialisation

La première étape, correspondant à la phase d’initialisation, consiste à définir les valeurs

initiales pour les vecteurs r et p. Le choix de r est en général non trivial, et peut induire
une convergence lente vers l’état d’équilibre thermodynamique du système. Par exemple,
dans le cas de l’ensemble canonique (voir la discussion à la Sec. 1.2.5) pour lequel on sup-
pose que le nombre de particules N , le volume V occupé par le système et la température
Simulations par dynamique moléculaire 14

Figure 2.2 – Schématisation des conditions aux limites périodiques en dimension 2 :

la cellule de base, grisée, est répliquée dans les deux directions de l’espace définies par
les vecteurs de la base canonique.

T de ce dernier sont constants (l’ensemble canonique est alors désigné par comme l’en-
semble NVT), on montre que la mesure M(r, p) dr dp du système à l’équilibre est une
mesure produit, appelée mesure de Maxwell-Boltzmann, telle que

1
M(r, p) = exp{−β H(r, p)} , (2.4)
Z

où Z est la fonction de partition (jouant le rôle de la constante de normalisation) et

β est un paramètre lié à la contrainte d’énergie moyenne. On montre par ailleurs que
β = 1/(kB T ), où kB est la constante de Boltzmann (kB = 1.3807 × 10−23 J/K) et T la
température du système. Les Eqs. (2.4) et (2.3) impliquent que :
— la mesure, notée Π(r)dr, définissant le vecteur aléatoire r (à l’équilibre) est non
Gaussienne et dépend en particulier du choix de V ;
— la mesure associée au vecteur p (à l’équilibre) est une mesure Gaussienne, ce qui
suggère l’utilisation d’une telle mesure pour l’échantillonnage de la valeur initiale
pour p (notons qu’il s’agit là d’une approximation, et que cette mesure dépend
explicitement de la température d’équilibre du système).
En pratique, une initialisation pour les vecteurs position peut donc être obtenue, soit par
un échantillonnage ad hoc de la mesure Π(r)dr (à l’aide de méthodes de type Markov
Chain Monte Carlo 1 ), soit au travers d’un couplage entre un algorithme déterministe ou
stochastique et une technique de relaxation.
1. Il s’agit ici d’un problème très délicat, en général formulé dans un espace de très grande dimension.
Simulations par dynamique moléculaire 15

1.2.2 Etape 2 : calcul des interactions

La seconde étape de la simulation concerne l’évaluation des interactions et le calcul des

forces sur chaque particule, selon la formule citée précédemment :

∂V(r)
fi = − , 16i6N . (2.5)
∂r i

La définition de l’énergie potentielle V s’avère critique ici, tant d’un point de vue de
la modélisation physique (le choix limitant en effet l’ensemble des interactions et des
configurations possibles, et donc les phénomènes physiques dont la simulation peut rendre
compte) que d’un point de vue numérique. Ces aspects seront abordés de façon plus
substantielle dans la suite du mémoire.

1.2.3 Etape 3 : intégration des équations de mouvement

La résolution des équations de la dynamique requièrent la construction d’un schéma

d’intégration approprié, satisfaisant en particulier les conditions suivantes :
— la préservation de l’énergie ;
— la réversibilité en temps ;
— le caractère symplectique.
Dans ce cadre, de nombreux schémas ont été proposés, tels que le schéma de Verlet
[134] et sa variante en vitesse [121], ou le schéma prédicteur-correcteur de Gear (voir par
exemple [35, 100]). Dans la suite, nous présenterons succinctement l’algorithme de Verlet
en vitesse, qui repose sur l’utilisation de développement de Taylor au second ordre. Par
un développement du vecteur position, il vient :

1
r(t + ∆t) ≈ r(t) + v(t)∆t + a(t)(∆t)2 , (2.6)
2

De façon similaire, on a
v(t + ∆t) ≈ v(t) + a(t)∆t (2.7)

et
v(t) ≈ v(t + ∆t) − a(t + ∆t)∆t ⇔ v(t + ∆t) ≈ v(t) + a(t + ∆t)∆t , (2.8)

si bien que l’approximation moyenne suivante pour la vitesse en t + ∆t est retenue :

a(t) + a(t + ∆t)

v(t + ∆t) ≈ v(t) + ∆t . (2.9)
2
Simulations par dynamique moléculaire 16

Ci-dessus, l’estimation du vecteur accélération est directement obtenue via la seconde

équation de Newton. Les Eqs. (2.6) et (2.9) définissent donc le schéma de Verlet en
vitesse (qui demeure le schéma le plus couramment utilisé dans les codes de dynamique
moléculaire), qui diffère du schéma de Verlet en ce que le vecteur des vitesses est mis à
jour au sein de l’itération. Le schéma s’avère également plus stable numériquement, car les
ordres d’approximation sont similaires en position et en vitesse. On peut enfin montrer que
ce schéma est bien réversible en temps, symplectique, et présente des bonnes propriétés
pour la conservation de l’énergie [44], ce qui se révèle particulièrement important dans le
cas de dynamiques en temps longs.

1.2.4 Etape 4 : évaluation de quantités d’intérêt

Lorsque l’équilibre thermodynamique est atteint, et pour une quantité d’intérêt Ψ ne

dépendant que des vecteurs positions, on cherche de façon générale à évaluer la moyenne
statistique : Z
< Ψ >:= Ψ(r) Π(r) dr . (2.10)
R3N

Dans le cadre de simulations par DM, et sous réserve que le théorème ergodique s’applique,
on procède à l’évaluation de l’intégrale précédente par la relation
Z υ
1
< Ψ >= lim Ψ(r(t)) dt , (2.11)
υ→+∞ υ 0

où l’origine du temps est prise dans le régime stationnaire. Naturellement, la vitesse de
convergence de l’estimateur dans le membre de droite de l’équation ci-dessus dépend de
la mesure à l’équilibre et peut être lente, notamment lorsque celle-ci contient plusieurs
bassins d’attraction (et en fonction de la nature des barrières entre ces puits).

Un point crucial à observer ici est que la quantité d’intérêt concerne en général une
propriété macroscopique (telle que le tenseur des déformations du volume de simulation)
intégrée, par la suite, dans un formalisme de mécanique des milieux continus – on trouvera
des applications de ce type d’approches dans [2, 50, 67, 99, 109] par exemple. La définition
de champs continus locaux (représentants par exemple une densité volumique ou des
contraintes de Cauchy) à partir de simulations par dynamique moléculaire doit par ailleurs
être accomplie avec prudence, et dépend essentiellement de l’état (en équilibre ou hors
équilibre) du système. Pour un système à l’équilibre thermodynamique, les expressions des
tenseurs de contrainte usuels (premier et second tenseurs de Piola-Kirchhoff, tenseur de
Cauchy) peuvent notamment être obtenues par la dérivation, par rapport au gradient de
Simulations par dynamique moléculaire 17

déformation, d’une énergie libre de Helmhotz. Dans le cas de systèmes hors équilibre, les
champs locaux peuvent être construits comme des valeurs moyennes de fonctions définies
dans l’espace des phases, selon la procédure dite de Irving-Kirkwood-Noll (voir [55] et
[60, 79], ainsi que [1] pour un exposé synthétique) – les grandeurs macroscopiques sont
alors définies comme les moyennes spatiales des champs microscopiques. Notons enfin
qu’une approche différente, basée uniquement sur une prise de moyenne spatiale, a été
proposée dans [46, 76] (voir [1] pour une discussion).

1.2.5 Ensembles thermodynamiques et schémas de contrôle

Comme nous l’avons vu précédemment, les simulations de dynamique moléculaire sont

accomplies sous l’hypothèse que certaines variables, telles que le nombre de particules, le
volume du domaine de simulation et l’énergie du système, sont constantes (ou constantes
en moyenne, suivant la nature des variables considérées) au cours du temps. En fonction
des variables retenues (en complément du nombre de particules), les différents ensembles
de micro-états cohérents avec les contraintes considérées sont appelés ensembles thermo-
dynamiques. De façon classique, on peut notamment citer :
— l’ensemble microcanonique (NVE), pour lequel le système possède un volume et une
énergie constants en moyenne (il s’agit donc d’un système isolé avec l’extérieur, peu
représentatif des situations expérimentales) ;
— l’ensemble canonique (NVT), pour lequel le système est en équilibre thermique avec
un bain de chaleur extérieur (l’échange d’énergie s’accomplit donc sous forme de
transfert thermique) et exhibe un volume et une température constants en moyenne ;
— l’ensemble isotherme-isobare (NPT), pour lequel la pression et la température du
système sont fixées en moyenne.
La modification d’un schéma de dynamique moléculaire afin d’imposer une contrainte
de température moyenne est appelée algorithme de thermostat (ou thermostat) – voir
par exemple [53]. La définition de ce dernier nécessite la définition d’une température
instantanée, comparable avec la température cible. La température instantanée peut être
définie à l’instant t comme [122]

N
1 X i i
T (t) = m kv (t)k2 , (2.12)
kB Nddi i=1

où Nddi désigne le nombre de degrés de liberté internes au système (Nddi = 3N en l’absence
de contraintes sur les degrés de liberté), et est telle que la moyenne de T (t) soit égale
Simulations par dynamique moléculaire 18

à la température macroscopique T . De nombreuses stratégies de modification ont été

proposées et on peut citer, de façon non exhaustive, les plus classiques :
— l’utilisation d’une dynamique de Langevin ;
— le recours à une méthode de couplage stochastique, dans laquelle les vitesses de cer-
tains atomes (sélectionnés aléatoirement selon un critère probabiliste) sont échantil-
lonnées de façon instantanée selon la mesure de Maxwell-Boltzmann – il s’agit du
thermostat d’Andersen [7] ;
— la considération d’un couplage faible basé sur la dynamique de Langevin, dans la-
quelle le terme de force stochastique est négligé et le coefficient de friction déterminé
en fonction de la température instantanée, de la température cible et d’un temps
de relaxation τB (on parle alors de thermostat de Berendsen [11]) ;
— la définition d’un système étendu, dans lequel un degré de liberté additionnel est
introduit afin de modifier conjointement le champ des vitesses et la variable tem-
porelle (il s’agit de l’approche originale proposée par Nosé [81]), ou uniquement le
champ des vitesses (comme proposé par Hoover [49]) – ce thermostat déterministe
est dit de Nosé-Hoover, et dépend en pratique du choix d’un paramètre de relaxation
en temps noté τNH .
Il est important de noter que le choix des paramètres de relaxation est important, car il
conditionne notamment l’échantillonnage de tel ou tel ensemble. Les techniques employées
peuvent par ailleurs être étendues/généralisées pour permettre un contrôle de la pression
macroscopique au travers de barostats (voir par exemple [11] pour le barostat de Beredsen,
ou encore [90] par exemple). Dans ce travail, l’ensemble des simulations de dynamique
moléculaire est conduit dans le code massivement parallèle LAMMPS [95] (les calculs
présentés par la suite sont réalisés sur 48 nœuds de calcul), avec des schémas de contrôle
de Nosé-Hoover.

2 Préparation des configurations initiales pour le na-

nocomposite

Dans cette section, nous présentons l’ensemble des systèmes considérés dans les simula-
tions par DM, ainsi que les différentes modélisations mises en œuvre.
Simulations par dynamique moléculaire 19

2.1 Définition des systèmes considérés

On considère un nanocomposite modèle de type matrice-inclusion(s), dans lequel la phase

matrice est constituée d’un polymère linéaire amorphe et la phase inclusion d’une ou
plusieurs nanoinclusions de silice (la modélisation de ces deux phases sont détaillée par
la suite). Deux valeurs pour la fraction volumique en renforts sont considérées, à savoir
4.8% et 28.5%. La longueur de chaque chaı̂ne de polymère (ou de façon équivalente, le
nombre de sites par chaı̂ne) est supposée constante et connue pour l’ensemble des systèmes
considérés. Soit npc le nombre de chaı̂nes considéré pour un système donné. On note par
ailleurs Rp et nS la taille et le nombre de nanoinclusions insérées dans le domaine de
simulation, avec

Rp ∈ {1.5, 3, 4.8, 6} , npc ∈ {10, 80, 320, 640} , (2.13)

en fonction des configurations (ci-dessus, les rayons sont exprimés en nm). Afin d’identifier
les résultats associés à un système donné, les conventions suivantes sont introduites :

Notation/Composition npc nS Rp (nm)

fv4.8–10c–1/R1.5 10 1 1.5
fv4.8–80c–1/R3 80 1 3
fv4.8–80c–8/R1.5 80 8 1.5
fv4.8–320c–1/R4.8 320 1 4.8
fv4.8–640c–1/R6 640 1 6

Table 2.1 – Description des systèmes considérés pour une fraction volumique en ren-
forts de 4.8%.

Notation/Composition npc nS Rp (nm)

fv28.5–10c–1/R3 10 1 3
fv28.5–10c–8/R1.5 10 8 1.5
fv28.5–80c–1/R6 80 1 6
fv28.5–80c–8/R3 80 8 3

Table 2.2 – Description des systèmes considérés pour une fraction volumique en ren-
forts de 28.5%.
Simulations par dynamique moléculaire 20

2.2 Modélisation de la matrice polymère

2.2.1 Description générale du système

Dans ce travail, la phase matrice est constituée d’un polymère prototypique linéaire
amorphe (cette dernière caractéristique conférant un caractère isotrope au matériau). Ce
dernier est constitué de npc chaı̂nes, chacune étant composée de 1000 sites CH2 représentés
par un modèle équivalent (approche de type “coarse-graining”). La masse molaire de
chaque unité CH2 est de 14.0273 g/mol ([18]). Il convient de noter que si le modèle intro-
duit est parfois utilisé dans la littérature pour la modélisation du polyéthylène, ce dernier
exhibe a contrario une structure semi-cristalline [14]. Une représentation d’une chaı̂ne et
des unités CH2 est fournie sur la Fig. 2.3.

Figure 2.3 – Visualisation d’une chaı̂ne et des sites CH2 (chaque boule représente un
site).

2.2.2 Définition des champs de force

Dans le cas de la modélisation d’une chaı̂ne polymère, l’expression de l’énergie potentielle

s’écrit classiquement comme suit :

V(r) = VB (r) + VUB (r) , (2.14)

où le terme VB représente l’ensemble des interactions entre m sites voisins le long de la
chaı̂ne (2 6 m 6 4), et VUB définit les interactions à distance. La notation i − j est
Simulations par dynamique moléculaire 21

introduite pour désigner ci-après, de façon symbolique, la liaison entre les atomes (ou
sites) i et j. Le premier terme VB est défini par
X X X
VB (r) = Vb (rij ) + Vθ (θijk ) + Vφ (φijk` ) , (2.15)
{i,j} ∈ Bb {i,j,k} ∈ Bθ {i,j,k,`} ∈ Bφ

avec r ij := r j − r i et rij = kr ij k. Dans l’Eq. (2.15), θijk représente l’angle (r\

ji , r jk ) formé

par les liaisons i − j et j − k, tandis que φijk` désigne l’angle entre les plans engendrés
par les couples (r ij , r kj ) et (r jk , r `k ). Les ensembles Bb , Bθ et Bφ sont les ensembles
de paires, triplets et quadruplets de sites successifs le long d’une chaı̂ne, toutes chaı̂nes
considérées par ailleurs. Il s’en suit que card(Bb ) = npc (nP − 1), card(Bθ ) = npc (nP − 2) et
card(Bφ ) = npc (nP − 3), avec nP le degré de polymérisation. La forme algébrique définie
par l’Eq. (2.15) comporte trois types de contribution, à savoir :
— un terme de rigidité à l’élongation, défini par

1
Vb (r) = Kb (r − r0 )2 , (2.16)
2

où r0 est une longueur de liaison d’équilibre et Kb un paramètre du modèle ;

— un terme de rigidité flexionnelle donné par

1
Vθ (θ) = Kθ (cos(θ) − cos(θ0 ))2 , (2.17)
2

avec θ0 l’angle de flexion à l’équilibre et Kθ un paramètre ;

— un terme de rigidité à la torsion, tel que

5
X
Vφ (φ) = Am cosm−1 (φ) , (2.18)
i=1

où {Am }5i=1 désigne un ensemble de paramètres et φ ∈ [−π, π].

En pratique, l’ensemble des paramètres ci-dessus est soit déterminé par une méthode
inverse formulée à partir de données expérimentales, soit estimé sur la base de simulations
ab initio. Le second terme VUB modélisant les interactions à distance est défini par la
relation suivante :
N X
X
VUB (r) = VLJ (rij ) , (2.19)
i=1 j>i
Simulations par dynamique moléculaire 22

où VLJ est le potentiel de Lennard-Jones 12-6 tronqué, dont l’expression s’écrit
σ 6
σ 12
VLJ (r) = 4ε − , r 6 rc , (2.20)
r r
= 0, r > rc . (2.21)

Ci-dessus, ε est la constante d’énergie du potentiel décrivant la profondeur du puit

d’énergie à son minimum (traduisant la force des interactions), et σ est la constante
de distance pour laquelle les forces d’attraction et de répulsion s’équilibrent (le potentiel
est donc nul pour r = σ). La grandeur rc représente le rayon de coupure au delà du-
quel les deux sites considérés ont une interaction nulle (ici, rc = 1.4 nm). Les valeurs des
différents paramètres des potentiels utilisées dans ce travail sont extraites de la littérature
(voir par exemple [18, 70]) et reportées dans le Tab. 2.3. Il convient de noter ici la liai-

Potentiel Paramètre Valeur Unité

VLJ ε 0.1133 kcal/mol
σ 0.43 nm
Vb Kb 70 000 kcal/(mol.nm2 )
r0 0.15 nm
Vθ Kθ 124.28 kcal/mol
θ0 112.81 deg
Vφ A1 2.1109 kcal/mol
A2 4.3229 kcal/mol
A3 1.1665 kcal/mol
A4 -7.6004 kcal/mol
A5 0 kcal/mol

Table 2.3 – Liste des paramètres des potentiels utilisés pour la modélisation du po-
lymère par dynamique moléculaire.

son quasi-rigide pour l’élongation entre deux sites adjacents (la longueur de liaison entre
ces derniers est de 0.153 nm, avec un coefficient de variation inférieur à 1%), ce qui im-
plique que la longueur des chaı̂nes ne varie pas au cours des simulations. Cette longueur
totale de chaı̂ne, égale à 153 nm environ, permet une estimation qualitativement satis-
faisante des différentes propriétés physiques, et des discussions concernant l’influence de
ce paramètre sur les prédictions macroscopiques peuvent être trouvées dans [17, 37, 77].
Les représentations graphiques des différents potentiels ainsi définis sont reportés sur les
Figs. 2.4 et 2.5.
Simulations par dynamique moléculaire 23

800 8

700 7

600 6

500 5

V θ (θ)
V b (r)

400 4

300 3

200 2

100 1

0 0
0 0.05 0.1 0.15 0.2 0.25 0.3 95 100 105 110 115 120 125 130
r θ

Figure 2.4 – Représentation graphique des potentiels (en kcal/mol) autour des valeurs
d’équilibre : cas des potentiels de rigidité à l’élongation r 7→ Vb (r) (gauche) et à la flexion
θ 7→ Vθ (θ) (droite).

7 200

6
150
5

100
V LJ(r)

4
V φ (φ)

3
50

2
0
1

0 ¹ 50
¹ 150 ¹ 100 ¹ 50 0 50 100 150 0.3 0.35 0.4 0.45 0.5 0.55 0.6
φ r

Figure 2.5 – Représentation graphique des potentiels (en kcal/mol) autour des valeurs
d’équilibre : cas du potentiel de rigidité à la torsion φ 7→ Vφ (φ) (gauche) et du potentiel
de Lennard-Jones 12-6 r 7→ VLJ (r) (droite).

2.2.3 Génération des configurations initiales

Il existe de nombreux algorithmes permettant de générer des configurations initiales du

système polymère à l’équilibre. Parmi celles-ci, on distingue principalement :
— la méthode de Monte Carlo et ses variantes, pour lesquelles la loi de distribution
des positions est échantillonnée et les nouveaux états acceptés selon un critère de
type Metropolis ;
— les méthodes de type marche aléatoire, qui produisent en général des configurations
hors équilibre, et qui sont donc par la suite combinées avec des techniques de relaxa-
tion (par dynamique moléculaire, par exemple) afin d’obtenir des configurations à
l’équilibre ;
Simulations par dynamique moléculaire 24

— des méthodes hybrides, dans lesquelles une méthode de Monte Carlo est utilisée
afin d’explorer les grandes modifications de conformation puis couplée de façon
séquentielle avec un algorithme de relaxation par dynamique moléculaire.
On notera que d’un point de vue théorique, le premier type d’approche est le plus perti-
nent, car il repose notamment sur un échantillonnage direct de la mesure invariante du
système. En pratique, cet échantillonnage peut toutefois s’avérer délicat, surtout pour des
systèmes de grande dimension. Dans ce travail, et par souci de simplicité, un ensemble
de configurations hors équilibre (pour un système donné, défini par le nombre de chaı̂nes
npc , avec npc ∈ {10, 80, 320, 640}) est tout d’abord généré par un algorithme de marche
aléatoire avec une condition de non recouvrement, appelé “Self Avoiding Random Walk”
(SARW) – voir [14, 50]. Deux exemples de configurations initiales obtenues par cette
technique sont représentées sur la Fig. 2.6, pour le système à 10 chaı̂nes.

Figure 2.6 – Visualisation de deux configurations initiales créées par l’algorithme

SARW, pour npc = 10 chaı̂nes.

Une relaxation par dynamique moléculaire est ensuite appliquée à chaque configuration
générée par SARW, par l’intermédiaire de schémas de Nosé–Hoover [49, 81] (avec des
temps de relaxation égaux à 2 ps et 1 ps pour la pression et la température, respective-
ment) et sous conditions de type NVT ou NPT. Pour l’ensemble des cas étudiés, le pas
de temps dans les simulations est fixé à 2 fs. La température cible pour l’équilibre final
est de 100 K, ce qui se situe bien en dessous de la température de transition vitreuse
du polymère (qui peut être estimée à environ 270 K ; voir le Tab. 2.12). Ce choix est
notamment motivé par le fait que dans ces conditions, le polymère exhibe un comporte-
ment élastique linéaire, et que les simulations de dynamique moléculaire permettent une
estimation satisfaisante des propriétés de volume.

La phase relaxation est accomplie comme suit.

— Dans un premier temps, la configuration générée par la méthode SARW est équilibrée
sous condition NVT à une température de 500 K, pendant une durée de 220 ps.
Simulations par dynamique moléculaire 25

Cette température induit une grande mobilité des chaı̂nes et favorise donc une ex-
ploration rapide des configurations à l’équilibre. Cette évolution est illustrée sur la
Fig. 2.7 ci-dessous.

Figure 2.7 – Relaxation de la configuration initiale, obtenue par l’algorithm SARW,

par dynamique moléculaire (système polymère à 10 chaı̂nes).

— Dans un second temps, une diminution de la température vers la valeur cible de

100 K est appliquée, toujours sous condition NVT. La vitesse de refroidissement est
égale à 1 K/ps.
— Enfin, un traitement NPT de 2000 ps est appliqué, avec des valeurs cibles de
température et de pression (isotrope) égales à 100 K et 0 bar, respectivement.
L’équilibre est alors obtenu après 500 ps, les 1500 ps suivantes étant destinées à
l’évaluation, par une estimation ergodique, de certaines grandeurs caractérisant no-
tamment le domaine de simulation.
Des configurations à l’équilibre obtenues par la procédure ci-dessus sont représentées
sur la Fig. 2.8, pour npc ∈ {10, 80, 320}. Dans le cas des systèmes à 10 et 80 chaı̂nes,
des ensembles de 20 configurations initiales indépendantes ont été générés. En raison
des temps calculs importants, une unique configuration à l’équilibre a été simulée pour
npc = 320 et npc = 640. Notons que ces dernières ne seront pas utilisées par la suite pour
l’identification des représentations stochastiques et servent de résultats de référence (en
terme de convergence vis-à-vis de la taille du système modélisé et simulé par dynamique
moléculaire) pour certaines propriétés estimées, notamment en terme de conformation. Le
tableau suivant synthétise les valeurs moyennes estimées de la longueur caractéristique
moyenne de la boı̂te de simulation à l’équilibre (la moyenne étant prise sur les trois
directions du repère cartésien et le cas échéant, sur l’ensemble des réalisations) et de la
densité du polymère.
Simulations par dynamique moléculaire 26

Figure 2.8 – Visualisation de configurations relaxées pour npc = 10, npc = 80 et

npc = 320 (de gauche à droite).

npc 10 80 320 640

Longueur caractéristique moyenne de la boı̂te (nm) 6.60 13.20 20.94 26.4
Densité volumique moyenne (kg/m3 ) 812.75 811.48 811.74 811.96

Table 2.4 – Quelques grandeurs caractéristiques des systèmes estimées à l’équilibre.

2.3 Simulation des inclusions de silice

L’objectif de cette section est de développer le modèle associé aux renforts nanoscopiques
de silice pour les simulations de dynamique moléculaire. Dans ce cadre, un modèle de
silice amorphe est retenu et permet d’obtenir un comportement macroscopique (i.e. à
l’échelle du domaine de simulation par DM) isotrope, même pour les faibles diamètres
d’inclusions (par exemple, pour npc = 10). La méthodologie générale est bâtie autour des
deux étapes suivantes [17–19, 58, 123, 132, 136] :
— tout d’abord, une structure cristalline de silice, de type α−quartz, est générée puis
rendue amorphe par une phase de relaxation par DM ;
— ensuite, un modèle élastique équivalent, dans lequel les interactions à distance sont
remplacées par des interactions de liaisons (faisant intervenir des potentiels har-
moniques dont les paramètres sont calibrés afin de garantir certaines propriétés de
conformation), est construit afin de réduire de façon substantielle les coûts calculs.
Ces étapes sont présentées dans la suite de cette section.
Simulations par dynamique moléculaire 27

2.3.1 Modélisation des interactions à distance : définition des champs de

force

La silice de type α−quartz est composée d’atomes de silicium (Si) et d’oxygène (O) (de
masses molaires respectives 28.0860 g/mol et 15.9994 g/mol), interagissant au travers
d’interactions à distance. Ces dernières font intervenir des interactions de type Van der
Waals et des interactions électriques. L’interaction de Van der Waals est modélisée par
un potentiel de Buckingham [132] :

r B3
VBCK (r) = B1 exp − − 6 , (2.22)
B2 r

où {Bi }3i=1 est un ensemble de paramètres dont les valeurs, extraites de la littérature,
sont fournies dans le Tab. 2.5 ci-dessous. Les interactions électriques sont modélisées par

Liaison B1 (eV) B2 (Å) B3 (eV /Å6 )

Si−Si 0 0.0657 0
Si−O 18,003.7572 0.2052 133.5381
O−O 1,388.7730 0.3622 175

Table 2.5 – Paramètres du potentiel de Buckingham.

la loi de Coulomb, formulée par le potentiel

q i qj
VC (rij ) = , (2.23)
4π0 rij

où qi et qj les charges électriques de l’atome i et j respectivement (qSi = 2.4e et qO =

−1.2e, avec e la charge élémentaire du proton), et 0 est la permittivité du vide.

2.3.2 Génération de la silice amorphe : prise en compte des interactions à

distance

La silice α−quartz est obtenue par la réplication, dans les trois directions de l’espace,
d’une cellule de base définie par des longueurs et des angles notés (a, b, c) et (α, β, γ)
respectivement [107]. Les valeurs de ces paramètres géométriques sont reportées ci-dessous
dans le Tab. 2.6.

D’un point de vue numérique, la réplication de cette cellule de base est réalisée à l’aide
de l’outil Pizza.py, disponible dans le logiciel LAMMPS. Cet outil est composé d’un en-
semble de scripts Python permettant des opérations de pré- et post-traitements pour les
Simulations par dynamique moléculaire 28

a (Å) b (Å) c (Å) α (deg) β (deg) γ (deg)

4.9137 4.9137 5.4047 90 90 120

Table 2.6 – Paramètres géométriques définissant la cellule de base pour la silice cris-
talline de type α−quartz.

simulations par DM. Notons ici que le nombre de réplications nécessaires pour les simula-
tions des systèmes nanocomposites (et donc, le coût calcul) dépend du nombre de chaı̂nes
considéré (puisque la fraction volumique en renforts est considérée comme constante).
Dans le cas du système à 10 (respectivement 80 et 320) chaı̂nes, 7 (respectivement 14 et
30) réplications sont nécessaires dans chaque direction. Un exemple de silice cristalline
simulée pour npc = 10 est représenté sur la Fig. 2.9. La dimension caractéristique finale

Figure 2.9 – Représentation d’une configuration de silice α−quartz simulée pour

npc = 10. Les atomes de silicium (Si) et d’oxygène (O) apparaissent en bleu et rouge,
respectivement.

du volume ainsi généré est de 3.6 nm, ce qui permet d’extraire une inclusion de rayon
1.5 nm (en conformité avec la fraction volumique cible de 4.8%). Notons qu’en raison du
temps calcul nécessaire pour générer des configurations cristallines de grande taille, ces
dernières peuvent être utilisées pour générer une population d’inclusions de plus faible
diamètre.

Afin d’obtenir des configurations amorphes, plusieurs étapes de relaxation par DM sont
ensuite appliquées (voir [17–19, 136]). Une première phase d’équilibrage sous condition
NPT à 300 K et pour une pression (isotrope) cible à 0 bar est accomplie. Pour ce faire,
des schémas de Nosé–Hoover sont utilisés comme algorithmes de contrôle. Dans un second
temps, la température cible est augmentée jusqu’à 10 000 K, assurant ainsi la transition de
la structure cristalline vers la structure amorphe. Un refroidissement vers la température
Simulations par dynamique moléculaire 29

de 100 K est ensuite imposé, suivi d’un traitement NPT conduisant à l’équilibre final à 100
K et pour une pression isotrope de 0 bar (notons que ces différentes étapes s’accompagnent
de variations dans la forme de la boı̂te de simulation, et que ces variations sont corrigées
jusqu’à l’obtention d’un domaine de simulation de forme cubique). La transition entre les
structures cristalline et amorphe peut être visualisée sur la Fig. 2.10.

Figure 2.10 – De gauche à droite : transformation de la structure cristalline vers

la structure amorphe à l’aide d’une relaxation par dynamique moléculaire (coupe bi-
dimensionnelle).

2.3.3 Construction d’un modèle de silice amorphe équivalent

En pratique, la simulation des inclusions de silice en intégrant les interactions à distance

(modélisées par les potentiels de Buckingham et de Coulomb) s’avère très coûteuse. Une
méthode de construction d’un modèle équivalent, ne faisant intervenir que des liaisons
élastiques, a été donc proposée dans [18] afin de réduire significativement les temps de
calcul. La stratégie consiste à insérer des liaisons élastiques fictives pour maintenir la
liaison Si−O et les angles de valence O−Si−O (principalement) et Si−O−Si de façon
équivalente au modèle intégrant les interactions à distance (voir [18] pour les détails
techniques). De façon qualitative, cette approche permet de conserver certaines propriétés
structurelles de la silice amorphe, ainsi que certaines propriétés physiques telles que la
densité volumique, l’énergie et les propriétés élastiques macroscopiques (ce point sera
illustré à la section 4.3). Une méthode alternative basée sur l’utilisation, dans le modèle
équivalent, d’un potentiel de Lennard-Jones peut être trouvée dans [9, 68]. Ici, les liaisons
élastiques sont modélisées par l’intermédiaire de potentiels de rigidité à l’élongation (défini
par l’Eq. (2.16)) et à la flexion (voir l’Eq. (2.17)). Les paramètres de ces potentiels sont
déterminés par une méthode inverse, sous la contrainte d’équivalence précédemment citée.
Les valeurs ainsi obtenues sont reportées, pour l’ensemble des liaisons en jeu, dans le
Simulations par dynamique moléculaire 30

Potentiel LiaisonParamètre Valeur Unité

Vb Si−O Kb 57573.01 kcal/(mol.nm2 )
r0 0.1625 nm
Vθ Si−O−Si Kθ 95.6022 kcal/mol
θ0 145 deg
O−Si−O Kθ 143.4034 kcal/mol
θ0 109.5 deg

Table 2.7 – Paramètres des potentiels remplacés pour la silice.

Tab. 2.7 ci-dessous. Une fois la substitution effectuée, les domaines de silice amorphe sont
ensuite soumis à une nouvelle relaxation par DM, à une température de 100 K et pour
une pression de 0 bar. Les configurations initiales alors obtenues sont ensuite utilisées afin
d’une part de déterminer (par des essais mécaniques virtuels, qui seront détaillés dans la
suite du mémoire) le comportement mécanique de la silice, et d’autre part d’extraire les
nanoinclusions qui seront insérées dans la matrice polymère.

2.3.4 Confrontation entre les deux modèles de silice

Dans cette section, on s’intéresse à la comparaison, en terme de longueur de liaisons

et d’angles de valence, entre le modèle avec interactions à distance (qualifié de modèle
complet) et le modèle avec liaisons élastiques (dit modèle équivalent). La comparison sur
la base des propriétés mécaniques sera présentée ultérieurement (voir la Sec. 4.3). Les
Figs. 2.11 et 2.12 présentent les estimations des densités de probabilité de la longueur
de liaison Si−O et des angles de valence Si−O−Si et O−Si−O (ce dernier étant le plus
critique), pour les deux modèles considérés (et associés au plus petit système, npc = 10).

Les valeurs moyennes extraites de ces données (la moyenne étant prise en espace et en
temps) pour les deux modèles, ainsi que celles de la densité, sont reportées sur le Tab. 2.8.

Propriété Modèle complet Modèle équivalent Unité

Longueur Si−O 0.1615 0.1625 nm
Angle O−Si−O 109.4223 109.3335 deg
Angle Si−O−Si 149.5454 141.1060 deg
Densité volumique 2318.13 2452.42 kg/m3

Table 2.8 – Comparaison des valeurs moyennes pour certaines propriétés structurelles
et pour la densité volumique.
Simulations par dynamique moléculaire 31

Densité de probabilité
40

0
0.155 0.16 0.165 0.17
Longueur de la liaison Si¹ O (nm)

Figure 2.11 – Graph de la densité de probabilité pour la longueur de liaison Si−O.

Trait continu : modèle complet. Trait discontinu : modèle équivalent.

0.1

0.09

0.08
Densité de probabilité

0.07

0.06 Liaison O¹ Si¹ O (modèle complet)

Liaison O¹ Si¹ O (modèle équivalent)
0.05
Liaison Si¹ O¹ Si (modèle complet)
0.04 Liaison Si¹ O¹ Si (modèle équivalent)

0.03

0.02

0.01

0
80 100 120 140 160 180 200
Angle de valence (deg)

Figure 2.12 – Graph des densités de probabilité pour les angles de valence Si−O−Si
(trait fin) et O−Si−O (trait épais). Trait continu : modèle complet. Trait discontinu :
modèle équivalent.

En ce qui concerne la longueur de la liaison Si−O, on observe donc des valeurs moyennes
similaires et une augmentation sensible de la variance lors du passage du modèle complet
au modèle équivalent. Dans le cas de l’angle O−Si−O (dont l’influence sur le compor-
tement macroscopique du système est la plus significative, car le nombre de liaisons de
ce type est important), on observe par ailleurs une bonne préservation de la moyenne et
de la variance. Enfin, on observe pour l’angle Si−O−Si des augmentations sensibles de
moyenne et d’écart-type.
Simulations par dynamique moléculaire 32

Pour le système faisant intervenir 7 (respectivement 14 et 30) réplications dans chaque

direction, la densité volumique moyenne est estimée à 2452.42 (respectivement 2441.93
et 2461.78) kg/m3 . Cette valeur est sensiblement inférieure à celle obtenue dans [18],
où la densité est estimée numériquement à 2518.6 kg/m3 . Enfin, et d’un point de vue
expérimental, la silice amorphe à 300 K et 1 bar de pression présente une densité volu-
mique de 2203 kg/m3 [107] (pour la structure quartz naturel, la densité volumique est
comprise entre 2635–2660 kg/m3 , selon [63]). On constate donc un accord qualitatif rai-
sonnable entre les propriétés calculées à partir du modèle complet et celles déterminées
à partir du modèle équivalent.

2.4 Simulation du nanocomposite

2.4.1 Principe d’assemblage

Afin d’obtenir des configurations initiales pour le nanocomposite modèle, la procédure

d’assemblage proposée dans [9, 17, 18] est suivie (voir [25, 77] pour des méthodologies
alternatives). Cette stratégie, représentée de façon schématique sur la Fig. 2.13, repose
sur trois étapes, à savoir :
— l’extraction, à partir du domaine cubique de silice amorphe, d’une inclusion sphérique
de rayon donné ;
— la gonflement du domaine occupé par la matrice polymère, puis l’insertion d’un
espace vide en son centre ;
— l’assemblage des deux systèmes précédents, suivi d’une phase de relaxation finale.

L’extraction d’une nanoinclusion de silice est effectuée suivant la méthodologie détaillée

dans [18, 19, 68] (voir également [9, 58, 77, 78] pour des discussions concernant la
modélisation de telles particules), dans laquelle l’existence de groupes silanols Si−O−H à
l’interface matrice-inclusion n’est pas prise en compte. Cette approche nécessite en parti-
culier la définition d’atomes d’oxygène et de silicium dans le cœur de l’inclusion, ainsi que
dans une zone externe couramment appelée écorce, pour lesquels des tables de connecti-
vité sont établies afin de préserver au maximum le caractère sphérique de l’inclusion. Ces
différents atomes sont schématisés sur la Fig. 2.14.

La répartition finale moyenne (i.e. en considérant l’ensemble des configurations initiales)

entre les atomes de silicium et d’oxygène obtenue dans les nanoinclusions de silice de
rayons variables figure dans le Tab. 2.9.
Simulations par dynamique moléculaire 33

Figure 2.13 – Principe de fabrication virtuelle du nanocomposite à partir d’une in-

clusion de silice et de la matrice polymère à l’équilibre (reproduction d’après [18]).

Figure 2.14 – Conﬁguration initiale de l’inclusion de silice après traitement des atomes
de surface : : atome d’oxygène dans le cœur de l’inclusion ; : atome de silicium
dans le cœur de l’inclusion ; ◦ : atome de silicium dans la couche externe de l’inclusion ;
* : atome d’oxygène dans la couche externe de l’inclusion connecté avec deux atomes de
silicium (contenus dans la couche externe) ; + : atome d’oxygène dans la couche externe
de l’inclusion connecté avec un atome de silicium (contenu dans la couche externe)

Des exemples de conﬁgurations obtenues sont représentées sur la Fig. 2.15 pour Rp ∈
{1.5, 3, 4.8} (nm).

Il convient de noter que la rugosité de la surface des inclusions est d’autant plus grande
que le rayon de celle-ci est faible (et dépend du diamètre du cœur utilisé dans la création
Simulations par dynamique moléculaire 34

Rp (nm) 1.5 3 4.8 6

Nombre d’atomes Si 272 2 460 10 978 21 776
Nombre d’atomes O 638 5 331 23 096 45 316
Nombre total d’atomes 910 7 791 34 074 67 092

Table 2.9 – Répartition du nombre d’atomes dans les nanoinclusions de silice (modèle
équivalent), pour Rp ∈ {1.5, 3, 4.8, 6} (nm).

Figure 2.15 – Visualisation de différentes configurations initiales des nanoinclusions :

Rp = 1.5 nm (gauche), Rp = 3 nm (centre) et Rp = 4.8 nm (droite).

des inclusions), ce qui a pour conséquence de générer des incertitudes plus importantes
pour des systèmes renforcés par des inclusions de petits rayons (par exemple, dans le cas
Rp = 1.5 nm).

Afin de pouvoir insérer la nanoinclusion dans la matrice polymère, la volume occupé par
celle-ci est tout d’abord dilaté. Pour ce faire, une augmentation de volume sensiblement
supérieure au volume de la sphère circonscrivant l’inclusion est appliquée. Par la suite, une
porosité de forme sphérique et de rayon Rpore > Rp est créée par l’intermédiaire d’une force
répulsive appliquée au centre de la boı̂te de simulation (en pratique, Rpore = Rp + 0.2).
Celle-ci prend la forme suivante :

Fpore (r) = −Kpore (r − Rpore )2 , r 6 Rpore (2.24)

= 0, r > Rpore , (2.25)

où r désigne la distance entre le site CH2 considéré et le centre de la boı̂te et Kpore est un
paramètre du modèle traduisant la raideur de la répulsion. L’introduction de cette force
a donc pour effet de repousser les chaı̂nes du polymère en hors de la porosité, comme
illustré sur la Fig. 2.16.

Lorsque la configuration du système polymère poreux est à l’équilibre, l’inclusion est

insérée : on obtient alors la configuration géométrique initiale. Cette dernière subit alors
une dernière phase de relaxation par dynamique moléculaire (contrôlée par des schémas
Simulations par dynamique moléculaire 35

Figure 2.16 – Création progressive de la porosité (de gauche à droite) par l’introduc-
tion d’une force répulsive au centre de la boı̂te (cas npc = 10).

de Nosé-Hoover), dans laquelle seules les interactions de Van der Walls entre les sites
CH2 et les atomes d’oxygène sont prises en compte par un potentiel de Lennard-Jones.
Les paramètres du potentiel sont estimés par la formule de Lorentz–Berthelot [6], et sont
donnés par σ = 0.31 nm et ε = 0.1288 kcal/mol respectivement. Le nanocomposite est
tout d’abord relaxé à 100 K pendant 220 ps. Une augmentation de température jusqu’à
500 K est ensuite appliquée, avec une vitesse d’échauffement de 1 K/ps. Le système est
alors maintenu à 500 K pendant 200 ps, puis refroidi à 100 K (avec une vitesse 1 K/ps).
Enfin, une dernière étape sous condition NPT, pour une température de 100 K et pour
une pression cible de 0 bar, est utilisée pendant 2000 ps afin d’obtenir une configura-
tion à l’équilibre, à partir de laquelle les caractérisations structurelle et physique seront
conduites. La Fig. 2.17 suivante présente deux configurations à l’équilibre du nanocom-
posite, avec npc = 10, Rp = 1.5 et Rp = 3 nm.

Figure 2.17 – Visualisation de deux configurations relaxées du nanocomposite pour

npc = 10, avec Rp = 1.5 (à gauche) et Rp = 3 (à droite). Les fractions volumiques sont
égales à 4.8% et 28.5%, respectivement.
Simulations par dynamique moléculaire 36

Il est intéressant de souligner que les boı̂tes de simulation sont cubiques à l’équilibre
(avec des déviations angulaires inférieures à 1 deg), ce qui traduit le caractère isotrope
des systèmes (polymères et nanocomposites).

2.4.2 Première caractérisation des nanocomposites

Lorsque chaque configuration est à l’état d’équilibre, les fractions massique et volumique
sont contrôlées afin d’obtenir des systèmes équivalents en terme d’effets de renforcement.
Notons ici que l’évaluation de la fraction volumique pour un système discret est délicate
(d’autant plus que les inclusions de faible diamètre sont rugueuses par ailleurs), et re-
pose en général sur une discrétisation du domaine de simulation à l’aide de tessellations
de Voronoi (voir notamment [32, 71, 84, 105]). Dans ce travail, une technique alternative
d’insertion de pores virtuels est utilisée [18], et permet une estimation qualitativement sa-
tisfaisante de la fraction volumique, pour un coût calcul raisonnable. Les caractéristiques
moyennes (suivant les configurations et le cas échéant, les directions de l’espace) ainsi ob-
tenues sont reportées sur les Tab. 2.10 et Tab. 2.11, où figurent également les longueurs
caractéristiques (c’est-à-dire les longueurs des mailles) des boı̂tes cubiques de simulation
par dynamique moléculaire.

npc 10 80 320 640

Rp (nm) 1.5 3 4.8 6
Longueur caractéristique de la boı̂te (nm) 6.67 13.37 21.34 26.81
Densité volumique (kg/m3 ) 884.16 888.43 882.52 881.41
Fraction massique 0.1128 0.1209 0.1312 0.1296

Table 2.10 – Paramètres des nanocomposites à l’équilibre, pour une fraction volu-
mique en renforts de 4.8%.

npc 10 80
Rp (nm) 3 6
Longueur caractéristique de la boı̂te (nm) 7.34 14.75
Densité volumique (kg/m3 ) 1231.67 1248.78
Fraction massique 0.5236 0.5347

Table 2.11 – Paramètres des nanocomposites à l’équilibre, pour une fraction volu-
mique en renforts de 28.5%.

Afin de détecter une éventuelle perte de symétrie du système lors des phases de création et
d’assemblage du nanocomposite, la position du centre d’inertie de l’inclusion par rapport
au centre du domaine de simulation est caractérisée. La Fig. 2.18 représente l’évolution
Simulations par dynamique moléculaire 37

de cette quantité en fonction du temps (sur un intervalle de 500 ps), pour les systèmes
définis par npc = {10, 80} et pour une fraction volumique de 4.8%.

0.2

0.18

0.16

0.14
distance (nm)

0.12
n pc = 10
0.1
n pc = 80
0.08

0.06

0.04

0.02

0
0 100 200 300 400 500
t (ps)

Figure 2.18 – Evolution de la distance entre le centre d’inertie de l’inclusion et le

centre de la boı̂te de simulation, pour npc ∈ {10, 80} et pour une fraction volumique de
4.8%.

On observe une distance moyenne du centre d’inertie au centre du domaine égale à 0.0527
nm pour le système à 10 chaı̂nes, et égale à 0.1512 nm pour le système à 80 chaı̂nes. Ces
distances sont très inférieures aux tailles des domaines de simulation (voir le Tab. 2.10),
ce qui confirme le centrage satisfaisant de l’inclusion dans les boı̂tes de simulation (rap-
pelons de plus ici que les inclusions ne sont pas parfaitement sphériques). Notons enfin
qu’en raison de la taille des systèmes considérés, le cas npc = 80 exhibe des fluctuations
statistiques sensiblement plus faibles.

2.4.3 Caractérisation de la mobilité atomique

La caractérisation de la mobilité de certains atomes et/ou sites peut être classiquement

accomplie en étudiant le déplacement quadratique moyen (MSD, suivant l’acronyme an-
glais courant), défini comme [14, 100]

NA Z tM SD
1 X
MSD(t0 , tM SD ) := kr A(i) (t + t0 ) − r A(i) (t0 )k2 dt , (2.26)
NA × tM SD i=1 0

où t0 est une origine arbitraire de temps (dans le régime stationnaire), tM SD désigne
l’amplitude de l’intervalle de temps sur lequel la moyenne est évaluée et A est la liste des
Simulations par dynamique moléculaire 38

indices des atomes et/ou des sites à étudier (de cardinal NA ). On peut donc, par exemple,
étudier la mobilité des sites CH2 dans la phase polymère, ou celle des atomes de silicium
et d’oxygène dans l’inclusion de silice – le coefficient de diffusion étant proportionnel à
la dérivée du MSD par rapport au temps, une pente plus grande traduit des phénomènes
de diffusion plus importants. L’évolution temporelle du déplacement quadratique moyen
(à 100 K) pour les atomes constituant l’inclusion et pour les sites de la matrice polymère
est représentée sur la Fig. 2.19.

0.5

0.45

0.4 nanopar t ic ule

polymè r e
0.35
MSD ( Å 2 )

0.3

0.25

0.2

0.15

0.1

0.05
0 0.5 1 1.5 2 2.5
pas de temps 5
x 10

Figure 2.19 – Evolution temporelle du déplacement quadratique moyen (avec tM SD =

500 ps) pour un nanocomposite à 80 chaı̂nes, dans l’inclusion et dans la phase polymère.

On constate que les atomes de silicium et d’oxygène sont moins mobiles que les sites
CH2 , ce qui s’explique par une valeur plus grande du coefficient de diffusion dans la phase
polymère [14]. La comparaison entre le MSD dans différents nanocomposites et pour le
polymère pur (avec npc = 80) est illustrée sur la Fig. 2.20.

On observe à nouveau une faible mobilité au sein des systèmes considérés, ce qui est une
conséquence de la faible température du système (100 K).

3 Caractérisation de la zone d’interphase

Dans cette partie, nous étudions plus spécifiquement certains propriétés locales des chaı̂nes
de polymère, afin d’obtenir des caractérisations géométrique et structurelle des segments
dans le voisinage de l’inclusion.
Simulations par dynamique moléculaire 39

0.5

0.45

0.4

0.35

MSD ( Å 2 )
0.3

0.25
p oly mè r e pur ave c 80 chaı̂ne s
0.2 nanoc omp os it e ave c n p c = 10
nanoc omp os it e ave c n p c = 80
0.15 nanoc omp os it e ave c n p c = 320
nanoc omp os it e ave c n p c = 640
0.1

0.05
0 0.5 1 1.5 2 2.5
pas de temps 5
x 10

Figure 2.20 – Evolution temporelle du déplacement quadratique moyen (avec tM SD =

500 ps) : comparaison entre le polymère pur à 80 chaı̂nes et différents nanocomposites
(fraction volumique : 4.8%).

3.1 Densité volumique locale

Afin de caractériser l’évolution de la densité locale du polymère proche de la surface de

l’inclusion, on considère une suite de coquilles sphériques, concentriques avec l’inclusion
et de volume constant noté Vρ . Chaque coquille est définie par un rayon intérieur Ri (le
rayon extérieur étant alors déterminé par le couple Vρ et Ri ), et repérée par le rayon
moyen Rm = (Ri + Re )/2. Notons que par ce procédé, l’épaisseur des coquilles diminue
lorsque Rm augmente. La densité radiale à l’instant t (dans la coquille de rayon moyen
Rm ) est estimée par la formule suivante :

Mp N (Rm , t)
ρen (Rm , t) = , (2.27)
NA Vρ

où Mp est la masse molaire du polymère, N (Rm , t) est le nombre d’unités CH2 dans la
coquille considérée à l’instant t et NA désigne le nombre d’Avogadro (NA = 6.022 ×
1023 mol−1 ). La densité radiale est enfin estimée par l’estimateur ergodique :
Z ∆
m 1
ρn (R ) = ρen (Rm , t) dt , (2.28)
∆ 0

avec ∆ un intervalle de temps donné. En pratique, des analyses de convergence doivent

être conduites sur les paramètres Vρ et ∆. Le graphique de la densité normalisée par celle
du polymère pur (notée ρp ) est présenté ci-dessous sur la Fig. 2.21.
Simulations par dynamique moléculaire 40

1.5

ρ n (r )/ρ p
nanoc omp osit e ave c n p c = 10
nanoc omp osit e ave c n p c = 80
0.5 nanoc omp osit e ave c n p c = 320

0
0 1 2 3 4 5 6
r − R p (nm)

Figure 2.21 – Graph de la fonction (r − Rp ) 7→ ρn (r)/ρp pour les systèmes nanocom-

posites définis par npc ∈ {10, 80, 320} (fraction volumique : 4.8%).

On observe donc des fluctuations du champ radial de densité volumique au voisinage

de l’inclusion, alors que le champ tend vers la valeur correspondant au polymère pur
lorsque r > Rp + 2 (nm), indépendamment du système considéré. Il convient également
de souligner que cette caractéristique est également exhibée pour une fraction volumique
de 28.5%.

3.2 Mobilité atomique au voisinage de la nanoinclusion

On s’intéresse ici à la mobilité des sites CH2 dans les couches précédemment introduites.
L’ensemble A contient, pour r fixé, l’ensemble des indices des atomes contenus dans la
coquille de rayon moyen r (le cardinal de A fluctue donc sensiblement avec le temps). Le
résultat est présenté dans la Fig. 2.22 suivante. On observe donc une perte de mobilité (de
l’ordre de 25%) des sites CH2 les plus proches de la surface de l’inclusion, et ce jusqu’à
une distance de la surface d’environ 2 nm.

3.3 Etude paramétrique

Dans cette section, on s’intéresse à l’influence de certains paramètres des simulations

(par dynamique moléculaire) sur les propriétés précédemment étudiées. Dans un premier
temps, l’influence de la longueur des chaı̂nes est caractérisée en observant l’évolution
radiale du champ de densité volumique. Le nombre total de sites CH2 est supposé constant
Simulations par dynamique moléculaire 41

0.3
0.29
0.28
0.27

MSD ( Å 2 )
0.26
0.25
0.24
0.23
0.22
0.21
0.2
0 5 10 15 20 25 30 35
r − R p ( Å)

Figure 2.22 – Moyenne des carrés de déplacements du polymère dans les coquilles
concentriques du nanocomposite avec 80 chaı̂nes à 100 K, moyenné sur un intervalle de
temps de 5 ps qui est répété 100 fois.

et égal à 10 000, et deux longueurs de chaı̂ne sont considérées, à savoir 50 et 1000. Dans le
premier cas, la matrice polymère contient donc 200 chaı̂nes, alors que le second système est
composé de 10 chaı̂nes. La représentation graphique de la densité volumique normalisée
est représentée sur la Fig. 2.23.

1.5

1
ρ n (r )/ρ p

0.5 200 chaı̂ne s ( 50 s it e s CH 2 par chaı̂ne )

10 chaı̂ne s ( 1000 s it e s CH 2 par chaı̂ne )

0
0 0.5 1 1.5 2
r − R p (nm)

Figure 2.23 – Graph de la fonction (r − Rp ) 7→ ρn (r)/ρp pour des systèmes exhibant

des longueurs de chaı̂nes différentes (Rp = 1.5 nm).

Dans les deux cas, on observe des fluctuations (suivant la direction radiale) de la den-
sité analogues à celles observées précédemment (bien que la fréquence des oscillations
Simulations par dynamique moléculaire 42

semble plus chaotique, eu égard à la taille des systèmes considérés), ce qui confirme la
faible incidence du paramètre de longueur de chaı̂ne (notons que cette influence dépend
naturellement de la quantité d’intérêt).

Le second paramètre vis-à-vis duquel la sensibilité des résultats doit être étudié est la
densité initiale de la matrice polymère, fixée jusqu’à présent à 500 kg/m3 . Ci-dessous,
la Fig. 2.24 présente l’évolution radiale de la densité volumique pour différents systèmes
(npc ∈ {10, 80} et pour différentes densités initiales.

1.6

1.4

1.2
ρ n (r)/ρ p

0.8

0.6 n pc = 10, de ns ité init iale de 500 kg/m 3

n pc = 10, de ns ité init iale de 700 kg/m 3
0.4
n pc = 80, de ns ité init iale de 500 kg/m 3
0.2 n pc = 80, de ns ité init iale de 700 kg/m 3

0
0 0.5 1 1.5 2 2.5 3 3.5
r − R p (nm)

Figure 2.24 – Graph de la fonction (r − Rp ) 7→ ρn (r)/ρp pour des systèmes générés à

partir de densités différentes.

On constate que la densité initiale n’a pas d’impact significatif sur la fréquence des oscil-
lations de la densité, confirmant ainsi la valeur de l’épaisseur de la zone d’interphase.

4 Essais mécaniques virtuels

Les sections précédentes ont été consacrées à la génération de configurations initiales,

à l’équilibre, pour le nanocomposite. Dans ce qui suit, nous nous intéressons à la ca-
ractérisation du comportement mécanique du nanocomposite, ainsi que de ses consti-
tuants (à savoir la matrice polymère modèle et l’inclusion de silice), au travers de simula-
tions de dynamique moléculaire. Les propriétés du polymère étant fortement sensibles à
la température, l’estimation de la température de transition vitreuse est considérée dans
un premier temps, à la section 4.1. La méthodologie d’estimation des propriétés élastiques
est ensuite présentée, ainsi que les résultats pour chaque système considéré.
Simulations par dynamique moléculaire 43

4.1 Calcul de la température de transition vitreuse

La température de transition vitreuse (notée Tg ) correspond à la température de passage

de l’état caoutchouteux à l’état vitreux, et se traduit notamment par un phénomène
de relaxation et une modification importante de la mobilité des unités (atomes et sites)
composant le système [33, 50]. Par conséquent, le passage en-dessous ou au-delà de Tg
se traduit par des évolutions importantes des propriétés physiques, telles que le volume
spécifique, le coefficient de dilatation, la viscosité ou bien encore les propriétés élastiques et
thermiques. De fait, l’identification précise et la modélisation des phénomènes physiques à
l’origine de cette transition d’état reste un sujet actif de recherche, et dépend en particulier
de la structure sous-jacente (amorphe, semi-cristalline, etc.) du polymère. Dans le cadre de
simulations de DM, l’estimation de Tg est en général réalisée en caractérisant l’évolution
du volume spécifique dans un cycle de charge-décharge en température. En effet, lorsque
le système est chauffé au delà de la température de transition vitreuse, la modification de
structure induite par la transition d’état implique un comportement différent lors que la
décharge : la température pour laquelle la rupture de pente apparaı̂t est alors considérée
comme une estimation raisonnable de Tg . Suivant la procédure détaillée dans [66], le
système est tout d’abord chauffé de 100 K à 500 K, puis subit une phase de relaxation
à 500 K pendant 500 ps. Le domaine est ensuite refroidi jusqu’à la température finale
de 100 K, avec une vitesse de −0.2 K/ps suffisamment faible pour ne pas perturber
l’équilibre thermodynamique à chaque décrément. Cette procédure est illustrée, dans le
cas du polymère pur et du nanocomposite (avec 80 chaı̂nes), sur la Fig. 2.25. Le Tab. 2.12

¹3
x 10
1.4
Volume spécifique (m 3 /kg)

1.3

1.2

1.1
100 150 200 250 300 350 400 450 500
Temp erature (K)

Figure 2.25 – Evolution du volume spécifique en fonction de la température (cycle de

décharge) pour différents systèmes (polymère pur et nanocomposite, npc = 80).
Simulations par dynamique moléculaire 44

suivant réunit les estimations de la température de transition vitreuse Tg pour plusieurs

systèmes.

Système Tg (K)
Polymère pur, npc = 10 264
Polymère pur, npc = 80 277
fv4.8–10c–1/R1.5 268
fv4.8–80c–1/R3 281

Table 2.12 – Estimation de Tg pour différents systèmes (polymère pur et nanocom-

posite).

Si l’on observe une variation sensible de Tg entre les systèmes polymère et nanocomposite
(notons que Tg ≈ 300 K pour le polyéthylène [50]), cette modification ne semble pas
significative (contrairement aux résultats obtenus dans [66], où une variation de 50 K est
constatée). Cette différence peut être expliquée par le fait que le procédé de caractérisation
est très sensible au bruit (et donc, au procédé de lissage des résultats numériques ou
au calcul ergodique). Dans la section suivante, on considère la simulation numérique
du comportement mécanique des systèmes en dessous de la température de transition
vitreuse.

4.2 Principe

Afin de réaliser des essais mécaniques virtuels, chaque configuration initiale testée est tout
d’abord relaxée par DM pendant 2000 ps. Les moyennes des caractéristiques géométriques
de la boı̂te de simulation (i.e. les trois longueurs et les trois angles formés par les vecteurs
de base) sont ensuite estimées sur un intervalle de temps de 1500 ps. Ces valeurs servent de
références pour la caractérisation des déformations macroscopiques subies par le domaine
de DM. Chaque configuration est alors soumise à des chargements mécaniques virtuels
sous conditions NPT [2, 67, 99, 109], modélisant ainsi les essais traditionnels considérés
en mécanique des milieux continus. Notons ici que plusieurs essais sont simulés de façon
indépendante à partir de la même configuration initiale. Le tenseur de pression cible est
écrit comme :
pαβ (t)
P αβ (t)ij = (δiα δjβ + δiβ δjα ) , (2.29)
2
où 1 6 α, β 6 3 sont des entiers indexant le chargement (le cas α = β = 1 correspondant
à un essai de traction suivant la direction définie par l’axe x1 par exemple), et t 7→ pαβ (t)
est une fonction décroissante monotone correspondant à une discrétisation quasi-statique
Simulations par dynamique moléculaire 45

d’une fonction de chargement linéaire. Dans le cas du polymère pur (respectivement de

la silice), cette dernière correspond à un décrément de charge de 10−4 GPa par 2 ps
(respectivement 10−3 GPa par 2 ps). La vitesse de chargement ainsi sélectionnée per-
met de ne pas modifier de façon significative l’équilibre thermodynamique du système,
et permet également de limiter les effets de viscosité. Lorsque d’un échelon de pres-
sion a été imposé, la convergence vers l’équilibre du système est contrôlée, puis l’état
de déformation macroscopique est estimé par le théorème ergodique. Cette procédure,
analogue dans son principe à la méthode d’homogénéisation numérique en mécanique
des milieux continus (voir [2, 18, 50, 88, 99, 101] pour d’autres approches, par exemple),
permet de reconstruire la réalisation du tenseur d’élasticité apparent associé à la confi-
guration initiale. Notons que cette dernière est sensiblement anisotrope, étant donné le
caractère non-sphérique de l’inclusion et la forme de la boı̂te de simulation. En introdui-
sant une hypothèse supplémentaire d’isotropie, l’estimation du module d’Young peut être
accomplie en réalisant un unique essai de traction suivant une direction définie par un
vecteur de base. En général, cette approche fournit des résultats qui différent sensible-
ment suivant la direction de sollicitation : le module final est alors défini, par convention,
comme la moyenne des différents modules obtenus. Dans ce cas, l’estimation du module
de compressibilité s’avère en général plus robuste, et peut être accomplie en mesurant
le changement de volume et son évolution par rapport à la pression appliquée [67] (voir
également [2]). Dans ce travail, deux types de résultats sont extraits de ces simulations,
à savoir
— les réalisations des tenseurs apparents, sans hypothèse de symétrie matérielle ;
— la projection de ces derniers sur l’ensemble des tenseurs isotropes.
Les résultats concernant les différentes phases, ainsi que le système nanocomposite, sont
présentés dans les sections suivantes.

4.3 Comportement mécanique de la silice

Considérons dans un premier temps la caractérisation de la silice, pour laquelle l’état de

pression macroscopique est choisi de telle façon que la déformation maximale soit égale
à 0.15%. D’un point de vue numérique, les essais sont réalisés sur une configuration de
silice amorphe, obtenue par relaxation de la structure cristalline à 7 réplications suivant
chaque direction. Les courbes de déformation-contrainte estimées pour les trois essais de
tractions et pour les trois essais de cisaillement sont représentées sur les Figs. 2.26 et
Simulations par dynamique moléculaire 46

0.12
t r ac t ion suivant x 1
t r ac t ion suivant x 2
t r ac t ion suivant x 3
0.1

0.08
σ (GPa)
0.06

0.04

0.02

0
0 0.5 1 1.5
ϵ ¹3
x 10

Figure 2.26 – Courbes de déformation-contrainte pour les essais de traction sur la

silice amorphe (T = 100 K).

0.1
c isaille me nt dans le plan ( x 2 , x 3 )
c isaille me nt dans le plan ( x 1 , x 3 )
0.08 c isaille me nt dans le plan ( x 1 , x 2 )
σ (GPa)

0.06

0.04

0.02

0
0 0.5 1 1.5
ϵ ¹3
x 10

Figure 2.27 – Courbes de déformation-contrainte pour les essais de cisaillement sur

la silice amorphe (T = 100 K).

2.27 suivantes (on rappelle que chaque point expérimental est calculé par le théorème
ergodique).

On observe clairement le caractère isotrope de la silice, ainsi que sa grande rigidité.

La réponse obtenue pour un essai de compression hydrostatique est représentée sur la
Fig. 2.28.

En réalisant cette séquence d’essais virtuels sur les 20 configurations initiales, le tenseur
apparent moyen (écrit suivant la convention de Kelvin-Voigt [27]) peut être estimé et vaut
Simulations par dynamique moléculaire 47

0.12

−(σ 11 + σ 22 + σ 33)/3 (GPa)

0.1

0.08

0.06

0.04

0.02

0
0 0.5 1 1.5 2 2.5 3
−(ϵ 11 + ϵ 22 + ϵ 33) ¹3
x 10

Figure 2.28 – Simulation d’un essai de compression hydrostatique sur la silice amorphe
(T = 100 K).

(en GPa) :
 
79.30 16.67 17.36 0.29 1.14 0.61
 

 77.17 14.80 1.53 0.66
1.07 
 
silice  77.48 −1.08 −0.23 1.14 
[C
e
DM ] = 
  . (2.30)
 61.86 0.71 −0.23 
 

 sym 59.68 0.32 
59.95

La projection du tenseur ci-dessus sur l’espace des tenseurs isotropes s’écrit alors (en
GPa) :  
77.48 16.51 16.51 0 0 0
 

 77.48 16.51 0 0 0 

 
iso e silice
 77.48 0 0 0 
P ([C DM ]) =   . (2.31)

 60.97 0 0 

 

 sym 60.97 0 

60.97

Des estimations des différents modules élastiques sont alors extraites, et sont résumées
avec plusieurs résultats d’autres études (numériques et expérimentales) dans le Tab. 2.13
suivant.

Notons que le résultat expérimental obtenu dans [107] concerne une silice amorphe à
une température (ambiante) de 300 K et à une pression de 1 bar. On observe donc
Simulations par dynamique moléculaire 48

Eq. (2.31) Réf. [66] (*) Réf. [67] (*) Réf. [107] (**)
Mod. de compressibilité 36.84 31.8 28 37
Mod. de cisaillement 30.49 27 26.3 31
Mod. d’Young 71.69 60 54 73
Coeff. de Poisson 0.18 0.19 0.223 0.17

Table 2.13 – Estimation des modules élastiques de la silice (à T = 100 K, en GPa) et
comparaison avec les données de la littérature (les travaux numériques et expérimentaux
sont signalés par les symboles (*) et (**) respectivement).

une bonne adéquation entre les estimations obtenues par les simulations de dynamique
moléculaire avec le modèle équivalent de silice et les travaux expérimentaux ou développés
avec des modèles intégrant des interactions à distance. Il convient enfin de noter que les
résultats sur les propriétés macroscopiques exhibent peu de bruit, et qu’il n’est donc
pas nécessaire d’estimer les propriétés élastiques de la silice pour des systèmes de plus
grandes dimensions (c’est-à-dire avec un nombre de réplications plus grand dans chaque
direction).

4.4 Comportement mécanique du polymère pur et du nanocom-

posite

Dans cette section, nous procédons à une étude similaire des propriétés mécaniques
élastiques du polymère pur et du système nanocomposite (pour npc = 10). Rappelons
que dans chaque cas, les systèmes subissent des traitements en température et en pression
identiques, afin d’imposer des traitements analogues, et que les propriétés sont évaluées
à une faible température (100 K) – température pour laquelle les matériaux exhibent
des comportements élastiques linéaires (rigides). Dans ce qui suit, les estimations sont
conduites pour une amplitude de déformation inférieure à 1%, avec un taux de charge-
ment de 10−4 GPa par 2 ps. Des exemples d’une configuration initiale (à l’équilibre) et
de la configuration déformée associée sont représentées sur les Figs. 2.29 et 2.30, pour un
essai de traction uniaxiale et pour un essai de cisaillement respectivement. Notons que
les états de déformation représentés sur ces figures ne correspondent pas à ceux utilisés
dans les estimations des propriétés élastiques, et sont suffisamment importantes pour per-
mettre une visualisation graphique des résultats. Les courbes de contrainte-déformation
correspondantes, ainsi que celles associées au système polymère pur, sont représentées
ci-dessous sur les Figs. 2.31 et 2.32. L’effet de renforcement de la nanoinclusion est
clairement mis en évidence, et s’avère plus prononcé dans le cas de la sollicitation en
Simulations par dynamique moléculaire 49

Figure 2.29 – Configurations initiale (gauche) et déformée à 10% (droite) du nano-

composite, pour un essai de traction uniaxiale. Les atomes Si et O sont représentés en
bleu et en rouge, respectivement.

Figure 2.30 – Configurations initiale (gauche) et déformée à 5% (droite) du nano-

composite, pour un essai de cisaillement. Les atomes Si et O sont représentés en bleu
et en rouge, respectivement.

cisaillement. Cet effet de renforcement est également bien observé au travers de l’essai de
compression hydrostatique, pour lequel les résultats sont représentés sur la Fig. 2.33.

Le tenseur d’élasticité apparent moyen du polymère pur, estimé à partir des 20 réalisations
indépendantes, est estimé (en GPa) à
 
6.78 4.66 4.68 0 −0.05 −0.05
 

 6.24 4.38 −0.07 −0.01 −0.03 

 
pp  6.23 −0.04 −0.04 0.03 
[C
e ]=
DM
 . (2.32)

 1.84 0.02 −0.01 
 

 sym 1.84 −0.04 
1.89
Simulations par dynamique moléculaire 50

0.12

0.1

σ 11 (GPa)
0.08

0.06

polymèr e pur
0.04
nanocompos ite

0.02

0
0 0.01 0.02 0.03 0.04 0.05
ϵ 11

Figure 2.31 – Courbes de contrainte-déformation pour le polymère pur et pour le

nanocomposite (à 100 K) : essai de traction uniaxiale.

0.1

0.08
σ 12 (GPa)

0.06

0.04

0.02
polymèr e pur
nanocompos ite
0
0 0.01 0.02 0.03 0.04 0.05
ϵ 12

Figure 2.32 – Courbes de contrainte-déformation pour le polymère pur et pour le

nanocomposite (à 100 K) : essai de cisaillement.

Les estimations des coefficients de variation pour les composantes a priori non nulles sont
données ci-dessous
 
0.20 0.23 0.21
 

 0.15 0.19 

 
 0.14 
[δ] =   , (2.33)

 0.07 

 

 sym 0.10 

0.07
Simulations par dynamique moléculaire 51

0.18
polymèr e pur

−(σ 11 + σ 22 + σ 33)/3 (GPa)

0.16 nanocompos ite
0.14

0.12

0.1

0.08

0.06

0.04

0.02

0
0 0.005 0.01 0.015 0.02 0.025 0.03
−(ϵ 11 + ϵ 22 + ϵ 33)

Figure 2.33 – Réponse du polymère pur et du nanocomposite pour un essai de com-

pression hydrostatique à 100 K.

et le coefficient de variation du tenseur apparent est égal à δ[Ce pp

DM ]
= 0.22. Il est important
de noter que le nombre de réalisations indépendants est ici bien trop faible pour permettre
d’obtenir des estimateurs statistiques avec un niveau de convergence satisfaisant. Ces
résultats permettent toutefois de dégager des tendances qualitatives, notamment pour la
comparaison entre les différents systèmes considérés (matrice pure et nanocomposite).

La projection d’un tenseur d’élasticité triclinique [C tri ] sur l’ensemble des tenseurs iso-
tropes est définie, au sens de la distance euclidienne, par les modules de compressibilité
k iso et de cisaillement µiso tels que [31] (voir également [80] ; voir [74] et les références
incluses pour une discussion générale)

1 tri
k iso = tri tri tri tri tri

C11 + C22 + C33 + 2(C23 + C31 + C12 ) (2.34)
9

1
µiso = tri tri tri tri tri tri tri tri tri

2(C11 + C22 + C33 − C23 − C31 − C12 ) + 3(C44 + C55 + C66 ) , (2.35)
30
Simulations par dynamique moléculaire 52

epp
avec {Cijtri }16i,j63 les composantes de [C tri ]. La projection du tenseur [C DM ] sur l’ensemble
des tenseurs isotropes s’écrit donc, en GPa :
 
6.42 4.57 4.57 0 0 0
 

 6.42 4.57 0 0 0 

 
iso e pp
 6.42 0 0 0 
P ([C DM ]) =   . (2.36)

 1.85 0 0 

 

 sym 1.85 0 

1.85

Les estimations des modules élastiques pour la matrice polymère sont listées dans le
Tab. 2.14. On constate un écart non négligeable entre les estimations obtenues dans ce

Eq. (2.36) Réf. [66] Réf. [67]

Mod. de compressibilité 5.19 4.20 4.40
Mod. de cisaillement 0.93 0.60 0.62
Mod. d’Young 2.62 1.70 1.86
Coeff. de Poisson 0.42 0.42 0.43

Table 2.14 – Estimation des modules élastiques (en GPa) du polymère pure avec 10
chaı̂nes à 100 K et leurs comparaisons.

travail et celles disponibles dans la littérature. Afin d’expliquer ces différences, on peut
souligner que :
— les codes de résolution et certains algorithmes de contrôle sont différents ;
— dans ce travail, les modules élastiques sont obtenus par identification sur le ten-
seur apparent moyen projeté sur l’ensemble des tenseurs isotropes, et non par une
régression directe sur les différents essais mécaniques ;
— l’algorithme de génération des configurations est différent par rapport aux travaux
présentés dans [66, 67], dans lesquels une méthode hybride basée sur une combinai-
son de Pivot de Monte Carlo et de relaxation par dynamique moléculaire est utilisée
(on rappelle qu’ici, les configurations sont obtenues par un algorithme de marche
aléatoire couplé à une relaxation par DM) ;
— le potentiel de LJ définissant les interactions de Van der Waals dans le polymère
est complet, c’est-à-dire qu’il contient des parties attractives et répulsives (contrai-
rement aux travaux de Brown et de ses collaborateurs, dans lesquels seule la partie
répulsive est considérée, avec une pression macroscopique de 5000 bar compensant
l’absence de la partie attractive).
Simulations par dynamique moléculaire 53

Enfin, il convient de noter que les modules estimés par une procédure de régression
violent en général les relations algébriques de l’élasticité, c’est-à-dire que la relation E =
9kµ/(3k + µ) (avec E le module d’Young) n’est pas vérifiée par exemple. A contrario, ces
relations sont satisfaites, par construction, dans la procédure présentée dans ce mémoire.

Pour le nanocomposite, le tenseur apparent moyen est (en GPa)

 
6.38 3.95 3.87 −0.09 −0.05 −0.03
 

 6.28 3.84 −0.1 −0.07
0.03 
 
ncp  6.13 −0.05 −0.09 −0.01 
[C DM ] = 
e  , (2.37)

 2.34 0.01 −0.05 
 

 sym 2.39 −0.01 
2.39

et la matrice des coefficients de variation pour les composantes est donnée par
 
0.0777 0.0892 0.0989
 

 0.0631 0.0909 

 
 0.0686 
[δ] =   . (2.38)

 0.0414 

 

 sym 0.0442 

0.0461

ncp
Le coefficient de variation du tenseur apparent [C
e ] est estimé à δ e ncp = 0.10. On
DM [C DM ]
constate donc que le tenseur d’élasticité apparent du nanocomposite présente un niveau
de fluctuation plus faible que le tenseur apparent associé à la matrice pure, ce qui peut
être expliqué par les faibles fluctuations exhibées par la silice. Le tenseur apparent moyen
projeté (sur l’ensemble des tenseurs isotropes) est (en GPa) :
 
6.26 3.89 3.89 0 0 0
 

 6.26 3.89 0 0 0 

 
iso e ncp
 6.26 0 0 0 
P ([C DM ]) =   , (2.39)

 2.38 0 0 

 

 sym 2.38 0 

2.38

ce qui permet de déduire des estimations des modules élastiques. Celles-ci sont résumées
dans le Tab. 2.15 suivant. En comparant les Tabs. 2.14 et 2.15, on note donc une aug-
Simulations par dynamique moléculaire 54

Eq. (2.39) Réf. [66] Réf. [67]

Mod. de compressibilité 4.68 4.20 4.34
Mod. de cisaillement 1.19 0.70 0.66
Mod. d’Young 3.29 1.90 1.90
Coeff. de Poisson 0.38 0.42 0.43

Table 2.15 – Estimation des modules élastiques du nanocomposite avec 10 chaı̂nes à

100 K et leurs comparaisons.

mentation significative du module de cisaillement pour le nanocomposite, tandis que

l’estimation du module de compressibilité d’avère plus faible. Dans le cas du module de
cisaillement, le résultat peut s’expliquer par la perte de mobilité locale (autour du renfor-
cement) des chaı̂nes. Cette dernière, associée à l’orientation préférentielle des segments
des chaı̂nes, entraı̂ne une augmentation de la rigidité locale dans la zone correspondant à
l’interphase, avec un effet significatif sur les propriétés en cisaillement. De façon similaire
aux résultats associés à la matrice polymère, les différences quantitatives observées entre
les résultats présentés ici et ceux proposés dans [66, 67] sont également imputables à la
stratégie de calcul des propriétés mécaniques, ainsi qu’au choix des champs de force et
de leurs paramètres (notons par ailleurs que les densités volumiques finales à l’équilibre
sont sensiblement différentes). Dans le cas du module de compressibilité, on notera que la
procédure de projection peut générer un biais lorsque les réalisations sont ne sont pas par-
faitement isotrope. Afin de clarifier ce point, le degré d’anisotropie peut être caractérisé
en évaluant l’écart entre les réalisations et leurs projections sur l’espace des tenseurs
isotropes. En introduisant la variable aléatoire dncp
iso mesurant cet écart dans le cas du
nanocomposite par exemple,

encp (ωi ) − P iso (C

kC encp (ωi ))k
dncp
iso (ωi ) := DM DM
, (2.40)
encp (ωi ))k
kP iso (C DM

e ncp , on note que les

encp (ωi ) désigne la i-ème réalisation du tenseur apparent C
où C DM DM
moyennes de cet écart à l’isotropie sont estimées à 0.129 et 0.078 pour le polymère pur et le
nanocomposite, respectivement. Par conséquent, le comportement du polymère présente
des fluctuations tricliniques plus importantes que le nanocomposite, ce qui peut être à
l’origine de l’estimation plus faible du module de compressibilité. Il est intéressant de sou-
ligner que le système polymère pur avec npc = 80 présente des fluctuations tricliniques
moins importantes, avec une distance moyenne à l’isotropie diminuée de moitié.
Simulations par dynamique moléculaire 55

5 Conclusion

Ce chapitre a été dédié à l’analyse, par des simulations de dynamique moléculaire, du

comportement mécanique (élastique) d’un nanocomposite modèle – constitué d’une ma-
trice polymère prototypique renforcée par une inclusion nanométrique de silice. Après
avoir détaillé les différentes interactions et processus de génération des configurations
initiales (suivant les travaux présentés dans [18]), nous nous sommes attachés dans un
premier temps à la caractérisation de la zone d’interphase. La faible mobilité des sites
CH2 au voisinage de l’inclusion, ainsi que l’orientation préférentielle des chaı̂nes, suggèrent
une modélisation géométrique de l’interphase comme une coquille sphérique d’épaisseur
constante, dans laquelle le tenseur des rigidités présente des fluctuations spatiales iso-
tropes transverses dans le système de repérage sphérique. Dans un second temps, des
essais mécaniques virtuels ont été conduits afin d’obtenir des réalisations indépendantes
du tenseur d’élasticité apparent associé au domaine de simulation – notons ici que l’ob-
jectif de ces calculs n’est pas de mettre en évidence un quelconque effet nano et que la
fraction volumique considérée dans cette étude est en général trop faible pour qu’un tel
effet soit constaté. Le caractère aléatoire exhibé par ce tenseur apparent provient d’une
part de l’aléa de microstructure (dans la mesure où les réalisations de l’inclusion, de faible
diamètre et non sphérique, différent d’une simulation à l’autre) et d’autre part, du bruit
intrinsèque généré par les simulations de DM. La modélisation probabiliste de ces objets,
à savoir du champ aléatoire des rigidités dans l’interphase et du bruit d’échantillonnage,
font précisément l’objet du chapitre suivant.
Chapitre 3

Modélisation stochastique du champ

des rigidités dans l’interphase

Dans ce chapitre, nous considérons la construction d’un modèle probabiliste pour le

champ aléatoire non Gaussien, à valeurs tensorielles, représentant l’élasticité (non iso-
trope) dans la zone d’interphase. D’un point de vue géométrique, et d’après les résultats du
Chap 2, cette dernière est modélisée comme une coquille sphérique d’épaisseur constante
dans laquelle les segments de chaı̂ne s’orientent préférentiellement de façon tangentielle
par rapport au vecteur de base radial (en coordonnées sphériques). Les hypothèses de
modélisation sont tout d’abord présentées à la Sec. 1. La méthodologie de construction,
basée en partie sur la Théorie de l’Information, est ensuite détaillée à la Sec. 2. La
définition de la loi marginale d’ordre un du champ est spécifiquement discutée à la Sec. 3,
tandis que le champ et son générateur de réalisations indépendantes sont introduits à
la Sec. 4. La stratégie de génération reposant sur la définition d’une famille d’équations
différentielles stochastiques d’Itô indexée en espace, un schéma d’intégration adapté est
par la suite discuté à la Sec. 5. Enfin, la modélisation du bruit induit par les simulations
de DM à l’aide d’un modèle de modèle de matrice aléatoire (à fluctuations statistiques
tricliniques) est présentée à la Sec. 6.

1 Hypothèses

Dans ce qui suit, on note (O, er , eθ , eϕ ) la base du système de coordonnées sphériques,

dans lequel un point générique de R3 est noté xs = (r, θ, ϕ). Soit x = (x1 , x2 , x3 ) un point

57
Modélisation stochastique 58

quelconque de R3 , exprimé dans le système de coordonnées cartésiennes, avec ∀x ∈ R3 :

x1 = r sin(ϕ) cos(θ) (3.1)

x2 = r sin(ϕ) sin(θ) (3.2)
x3 = r cos(ϕ) (3.3)

où r ∈ [0, +∞[, θ ∈ [0, 2π] et ϕ ∈ [0, π]. Sur la base des résultats issus des simulations par
dynamique moléculaire, les hypothèses de modélisations suivantes peuvent être formulées.
— D’un point de vue géométrique, la région de l’interphase peut être modélisée comme
une coquille sphérique d’épaisseur constante eI = 2 nm :

DI := {xs = (r, θ, ϕ)|r ∈ [Rp , Rp + eI ], θ ∈ [0, 2π], ϕ ∈ [0, π]} , (3.4)

avec Rp le rayon de la nanoparticule associée. Rappelons que Rp dépend du système

considéré et en particulier, du nombre de chaı̂nes npc inclues dans la boı̂te (puisque
les différentes configurations ont une fraction volumique en renfort constante) et de
la température.
— Le matériau dans l’interphase exhibe un comportement local isotrope transverse,
défini au point xs par le vecteur normal normé n(xs ) = er (xs ).
— La valeur moyenne du champ aléatoire du tenseur d’élasticité dans l’interphase
est indépendante du point considéré en coordonnées sphériques (en coordonnées
cartésiennes, la valeur moyenne dépend du point x et le milieu associé apparaı̂t donc
comme hétérogène) – s’il s’agit ici d’une hypothèse de modélisation simplificatrice,
il convient de souligner par exemple que l’orientation tangentielle des segments de
chaı̂ne ne varie que de façon très modérée dans la direction radiale.
— Les propriétés probabilistes du champ aléatoire du tenseur d’élasticité dans l’in-
terphase sont identiques dans les directions définies par les vecteurs de base eθ et
eϕ .
Il convient de noter que l’hypothèse concernant la valeur moyenne peut être modifiée sans
difficulté, en ajoutant par exemple une dépendance dans la direction radiale - le milieu
présenterait alors un gradient de propriétés. Toutefois, une telle modélisation nécessiterait
des informations expérimentales très riches en vue de la calibration du modèle, ce qui n’est
pas le cadre d’étude retenu dans ce travail.
Modélisation stochastique 59

2 Méthodologie de construction

Soit E sym (n) le sous-ensemble de M+

6 (R) correspondant à l’ensemble des tenseurs du
second-ordre isotropes transverses par rapport au vecteur normal normé n. On note
{[C s (xs )], xs ∈ DI } le champ aléatoire à valeurs dans E sym (n) (avec n = (1, 0, 0))
représentant le champ des rigidités dans l’interphase. Dans un premier temps, et pour
tout xs de DI , on considère la décomposition algébrique suivante :

∀xs ∈ DI , [C s (xs )] = η[I6 ] + [Ms (xs )] , (3.5)

où η 1 est un paramètre introduit afin de garantir la propriété d’ellipticité uniforme

pour le problème aux limites stochastique associé, [I6 ] est la matrice identité (6 × 6)
et {[Ms (xs )], xs ∈ DI } est le champ aléatoire à valeur dans E sym (n). En introduisant
ensuite la valeur moyenne Ms = E{[C s (xs )]} − η[I6 ], définie-positive et indépendante de
xs dans le repérage sphérique, on considère la normalisation suivante :

∀xs ∈ DI , [Ms (xs )] = [Ms ]1/2 [N s (xs )] [Ms ]1/2 , (3.6)

où {[N s (xs )], xs ∈ DI } est un nouveau champ aléatoire à valeurs dans E sym (n) et tel que
E{[N s (xs )]} = [I6 ], ∀xs ∈ DI . Clairement, la matrice déterministe [Ms ]1/2 est inversible,
ce qui implique que la matrice aléatoire

[N s (xs )] = [Ms ]−1/2 [Ms (xs )] [Ms ]−1/2 (3.7)

est à valeurs dans E sym (n). En introduisant la représentation matricielle {[E i ]}5i=1 de la
base de Walpole [140] pour l’ensemble E sym (n), définie par les tenseurs d’ordre quatre
suivants :

JE1 K := [p] ⊗ [p] , JE2 K := 12 [q] ⊗ [q] , JE3 K := √1 ([p]

2
⊗ [q] + [q] ⊗ [p]) ,
(3.8)
JE4 K := [q]⊗[q] − JE2 K , JE5 K := JIK − JE1 K − JE2 K − JE4 K ,

avec [p] := n ⊗ n, [q] := [I3 ] − [p] et JIKijkl := (δik δjl + δil δjk )/2, il vient

5
X
s s
[N (x )] = Nis (xs )[E i ] , (3.9)
i=1

où le champ aléatoire {N s (xs ) = (N1s (xs ), . . . , N5s (xs )), xs ∈ DI } est à valeurs dans le
sous-ensemble Vsym ⊂ R5 (rappelons que la classe considérée ici est l’isotropie transverse)
Modélisation stochastique 60

tel que

Vsym := {v ∈ R5 | v1 > 0 , v2 > 0 , v4 > 0 , v5 > 0 , v1 v2 − v32 > 0} . (3.10)

La génération de réalisations indépendantes du champ {N s (xs ), xs ∈ DI } nécessite donc

la prise en compte de la contrainte de support définie par l’Eq. (3.10), ce qui peut s’avérer
problématique avec des techniques classiques de type Metropolis-Hastings (où le taux
rejet peut être très important dans le cas de variances élevées). Afin de contourner cette
difficulté, on introduit le champ aléatoire {[G s (xs )], xs ∈ DI } tel que

∀xs ∈ DI , [N s (xs )] = exp{[G s (xs )]} (3.11)

où [G s (xs )] admet la décomposition algébrique

5
X
s s
[G (x )] = Gsi (xs )[E i ] . (3.12)
i=1

On montre que ce champ aléatoire est défini de façon unique, et que {Gsi (xs ), xs ∈ DI },
i ∈ {1, . . . , 5}, est un champ aléatoire scalaire non-Gaussien à valeurs dans R [41]. La
transformation exponentielle permet donc de relaxer la contrainte de support associée
au caractère défini-positif de [N s (xs )], ce qui constitue un ingrédient clé en vue de la
construction d’un générateur de réalisations pour les champs introduits. Introduisons
enfin le champ aléatoire {Gs (xs ), xs ∈ DI } à valeurs dans R5 tel que

∀xs ∈ DI , Gs (xs ) = (Gs1 (xs ), . . . , Gs5 (xs )) . (3.13)

Dans ce qui suit, on s’intéresse donc à la construction d’un modèle probabiliste pour le
champ aléatoire {G(xs ), xs ∈ DI }. Pour ce faire, et suivant [41, 114], la méthodologie (à
deux étapes) suivante est adoptée.
— Dans un premier temps, on s’intéresse à la construction de la famille de lois mar-
ginales d’ordre un pour le champ aléatoire {G(xs ), xs ∈ DI }. Cette étape consiste
à définir la densité de probabilité de la variable aléatoire G(xs ) à valeurs dans R5 ,
pour tout xs de DI . Cette construction est accomplie dans le cadre de la Théorie
de l’Information et plus spécifiquement, en invoquant le principe du Maximum de
l’Entropie [56, 108].
— Dans un second temps, le champ aléatoire {G(xs ), xs ∈ DI } est défini par l’in-
termédiaire d’une transformation non-linéaire locale d’un champ Gaussien centré
sous-jacent, à valeurs dans R5 , telle que {G(xs ), xs ∈ DI } exhibe comme famille de
Modélisation stochastique 61

lois marginales d’ordre un la famille définie lors de la première étape. Notons dès à
présent que le champ non-Gaussien ainsi défini possède une structure de corrélation
héritée de la transformation du champ Gaussien sous-jacent.

3 Construction de la loi marginale d’ordre 1

3.1 Définition des contraintes

Soit xs fixé dans DI . On considère donc la construction explicite de la densité de pro-

babilité pGs (xs ) de la variable aléatoire Gs (xs ). Comme mentionné précédemment, cette
construction repose sur le principe du Maximum d’Entropie, pour lequel les contraintes
algébriques (sous forme d’espérances mathématiques) traduisent l’information objective
disponible sur la variable considérée. Dans ce travail, les deux contraintes suivantes sont
considérées [41] : ( !)
5
X
E exp Gsi (xs )[E i ] = [I6 ] , (3.14)
i=1

5
X
E {Gsi (xs )} tr([E i ]) = ν , |ν| < +∞ , (3.15)
i=1

où tr(·) désigne l’opérateur de trace et ν un paramètre supposé indépendant de xs .

La contrainte définie par l’Eq. (3.14) est induite par la normalisation introduite dans
l’Eq. (3.6), et signifie que la valeur moyenne du champ est supposée connue. La seconde
contrainte définie par l’Eq. (3.15) traduit la finitude des moments d’ordre deux pour la
matrice aléatoire [N s (xs )] et son inverse [N s (xs )]−1 [113] :

E{k[N s (xs )]k2 } < +∞ , E{k[N s (xs )]−1 k2 } < +∞ . (3.16)

3.2 Solution explicite

On peut alors montrer qu’une approximation de la densité de probabilité recherchée s’écrit

[118]
∀g ∈ R5 , pGs (xs ) (g) = c exp (−Φ(g)) , (3.17)
Modélisation stochastique 62

avec c la constante de normalisation et Φ la fonction potentiel définie comme suit :

5
! 5
!
X X
−2
∀g ∈ R5 , Φ(g) = α × δ[N ] × tr exp gi [Ei ] − gi [Ei ] . (3.18)
i=1 i=1

Dans l’Eq. (3.18), α est un paramètre qui dépend de la symétrie matérielle considérée
(α ≈ 0.82 pour l’ensemble E sym (n) des tenseurs isotropes transverses, quel que soit le
vecteur n), et δ[N ] est un paramètre scalaire indépendant de xs permettant de contrôler
le niveau des fluctuations statistiques exhibées par [N s (xs )]. Dans les sections suivantes,
le champ aléatoire {G(xs ), xs ∈ DI } et son générateur de réalisations sont définies par
l’intermédiaire de la définition d’une famille de processus de diffusion indexée par xs ∈ DI .

4 Définition et générateur du champ aléatoire

4.1 Définition d’une famille de processus de Wiener normalisés

Soit {Ξ(xs ) = (ξ 1 (xs ), . . . , ξ 5 (xs )) , xs ∈ DI } le champ aléatoire Gaussien centré, à va-

leurs dans R5 et de composantes statistiquement indépendantes, définies comme suit. Soit
(xs , y s ) 7→ [RΞ (xs , y s )] la fonction de corrélation continue du champ {Ξ(xs ), xs ∈ DI },
à valeurs dans l’ensemble des matrices réelles diagonales et telle que

[RΞ (xs , y s )]ii = E{ξ i (y s )ξ i (xs )} =: ρ (xs , y s ) (3.19)

avec
[RΞ (xs , xs )]ii = 1 , 16i65 (3.20)

pour tout couple (xs , y s ) dans DI × DI . Par construction, la structure de corrélation

du champ aléatoire est complètement définie par la fonction de corrélation normalisée
ρ. Cette dernière est supposée exhibant une structure séparable, qui s’écrit ici pour tous
couples xs = (r, θ, ϕ) et y s = (r0 , θ0 , ϕ0 ) dans DI :

∀(xs , y s ) ∈ DI × DI , ρ (xs , y s ) := ρr (τr ) × ρθ (r, r0 , τθ ) × ρϕ (r, r0 , τϕ ) , (3.21)

où τr := r − r0 , τθ := θ − θ0 et τϕ := ϕ − ϕ0 sont les variables des translations radiale

et orthoradiales. Ici, les fonctions de corrélation ci-dessus sont définies au travers des
Modélisation stochastique 63

paramétrages suivants :
 2
2Lr 2 πτr
ρr (τr ) := sin ,





 πτr 2L r !
0 2

1 r + r τ

θ

ρθ (r, r0 , τθ ) := exp − sin2 , (3.22)

 2 L θ 2
!
0 2

1 r + r τ

ϕ

0
 ρϕ (r, r , τϕ ) := exp − 2 sin2 .



Lϕ 2

Ci-dessus, Lr représente la longueur de corrélation spatiale selon la direction radiale,

tandis que Lθ et Lϕ est un ensemble de paramètres qui contrôlent les taux de décorrélation
selon les directions définies par eθ et eϕ . Pour des valeurs données de r et r0 dans [0, Rp ],
les fonctions ρθ et ρϕ sont 2π-périodiques et satisfont les propriétés suivantes, liées à la
symétrie sphérique du problème :
— ρθ (r, r0 , −τθ ) = ρθ (r, r0 , τθ ) pour tout τθ dans [−2π, 2π] ;
— ρϕ (r, r0 , −τϕ ) = ρϕ (r, r0 , τϕ ) pour tout τϕ dans [−π, π].
Il convient de noter que d’autres formes algébriques, cohérentes avec la symétrie du
problème, peuvent être utilisées pour les fonctions de corrélation ρr , ρθ et ρϕ . De même, les
fonctions de corrélation suivant eθ et eϕ peuvent être représentées par des développements
en séries de Fourier. Bien que plus générales, ces représentations dépendent toutefois d’un
nombre important de paramètres (suivant la qualité des approximations recherchées), ce
qui s’avère peu adapté au cadre retenu d’une identification inverse à partir de données
limitées. Des graphiques de la fonction de corrélation normalisée sont proposés sur les
Figs. 3.1 et 3.2, pour différentes valeurs des paramètres Lr et Lθ (pour cette application,
τϕ = 0). Ces illustrations permettent de visualiser clairement l’incidence de chaque pa-
ramètre sur la structure de corrélation du champ Gaussien considéré. Introduisons enfin
le champ aléatoire Gaussien W = {W (t, xs ), t > 0, xs ∈ DI }, centré et à valeurs dans
R5 , tel que :
— ∀xs ∈ DI , W (0, xs ) = 0 presque-sûrement (p.s.) ;
— la dérivée généralisée Dt W de W (par rapport au paramètre t) est le bruit blanc
Gaussien cylindrique noté B, pour lequel la fonction de covariance [CB ] est donnée,
pour ∀(xs , y s ) ∈ DI × DI et pour tout τ ∈ R, par la relation suivante :

[CB (xs , y s , t + τ, t)] := δ0 (τ )[RΞ (xs , y s )] , (3.23)

avec δ0 la distribution de Dirac au point 0.

Par construction, W = {W (t, xs ), t > 0, xs ∈ DI } est tel que :
Modélisation stochastique 64

2
2
0
0
¹2 ¹2

Figure 3.1 – Graphique de la fonction de corrélation normalisée (xs , y s ) 7→ ρ (xs , y s ),

avec xs = (1.5, 0, π/2) et y s = (r, θ, π/2) pour 1.5 6 r 6 3.5, 0 6 θ 6 2π, Lr = 0.2 et
Lθ = 0.5.

1 1

2 2
2 2
0 0
0 0
¹2 ¹2 ¹2 ¹2

Figure 3.2 – Graphique de la fonction de corrélation normalisée (xs , y s ) 7→ ρ (xs , y s ),

avec xs = (1.5, 0, π/2) et y s = (r, θ, π/2) pour 1.5 6 r 6 3.5 et 0 6 θ 6 2π. Figure de
gauche : Lr = 2, Lθ = 0.5. Figure de droite : Lr = 0.2, Lθ = 2.

— pour xs fixé dans DI , W = {W (t, xs ), t > 0} est un processus de Wiener normalisé

(appelé encore mouvement Brownien standard) à valeurs dans R5 – notons ici que
la variable t ne doit pas être confondue avec une variable de temps et en particulier,
avec la variable de temps utilisée dans les simulations par dynamique moléculaire ;
— pour t fixé dans [0, +∞[, W = {W (t, xs ), xs ∈ DI } est un champ Gaussien coloré.
Dans ce qui suit, la famille de processus de Wiener indexée par DI va être utilisée afin
de construire un générateur de réalisations robuste, basé sur la définition d’une famille
de processus de diffusion. Ce point fait l’objet de la prochaine section.
Modélisation stochastique 65

4.2 Générateur de réalisations

Soit xs ∈ DI un point fixé dans la zone d’interphase. On introduit le processus de Markov

{(U (t, xs ), V (t, xs )), t > 0}, à valeurs dans R5 × R5 , satisfaisant l’équation différentielle
stochastique d’Itô [41] :

 dU (t, xs ) = V (t, xs ) dt
% √ , (3.24)
 dV (t, xs ) = −∇u Φ(U (t, xs )) − V (t, xs ) dt + η dW (t, xs )
2

pour tout t ∈ R+ , où Φ est la fonction potentiel définie par l’équation Eq. (3.18),
% ∈]0, +∞[ est un paramètre numérique à ajuster (sur la base d’analyse de vitesses
de convergence, par exemple) et {W (t, xs ), t > 0} est le processus de Wiener défini à
la section précédente. Les conditions initiales sont données par U (0, xs ) = U 0 (xs ) et
V (0, xs ) = V 0 (xs ) presque sûrement, où la distribution de probabilité de la variable
aléatoire (U 0 (xs ), V 0 (xs )) est supposée connue. On montre par ailleurs que Φ est telle
que la fonction u 7→ k∇u Φ(u)k :
— est localement bornée sur R5 ;
— vérifie la propriété : E{k∇u Φ(u)k} < +∞.
On montre alors le résultat de convergence suivant [115] :

loi
lim U (t, xs ) = Gs (xs ) . (3.25)
t → +∞

L’intégration simultanée des éléments de la famille d’équations différentielles stochas-

tiques d’Itô, dont chaque élément est défini à xs fixé par l’Eq. (3.24), permet donc de
définir un algorithme d’échantillonnage pour le champ aléatoire {Gs (xs ), xs ∈ DI }. Un
schéma numérique d’intégration associé est proposé à la section suivante.

5 Schéma de discrétisation

Dans ce travail, l’équation différentielle stochastique définie par l’équation Eq (3.24) est
discrétisée par un schéma de Störmer-Verlet [20]. Soit ∆t le pas d’intégration, et soit
tk = (k − 1)∆t, k > 1, une discrétisation régulière de l’intervalle d’intégration. Pour
Modélisation stochastique 66

k > 1, le schéma de Störmer-Verlet s’écrit

 ∆t (k) s
 U (k+1/2) (xs ) = U (k) (xs ) + V (x )
2


 √

(k+1) 1 − ζ (k) s ∆t (k+1/2) %
V s
(x ) = V (x ) − ∇u Φ(U s
(x )) + ∆W (k+1) (xs ) ,
 1 + ζ 1 + ζ 1 + ζ
 U (k+1) (xs ) = U (k+1/2) (xs ) + ∆t V (k+1) (xs )



2
(3.26)
avec U (k) (xs ) := U (tk , xs ), V (k) (xs ) := V (tk , xs ), ζ := %∆t/4 et
√
∆W (k+1) (xs ) = W (tk+1 , xs ) − W (tk , xs ) := ∆t Ξ(xs , θk+1 ) . (3.27)

Ci-dessus, xs 7→ Ξ(xs , θk+1 ) est la (k + 1)-ième réalisation indépendante du champ

aléatoire Gaussien vectoriel {Ξ(xs ), xs ∈ DI } introduit à la section 4.1. Les conditions
initiales sont données par U (1) (xs ) = u(0) et V (1) (xs ) = v (0) , avec u(0) et v (0) deux
vecteurs arbitraires dans R5 . La condition de convergence définie par l’Eq. (3.25) s’écrit
alors :
s s s
G (x ) = lim lim U (tk , x ) . (3.28)
∆r ↓ 0 k → +∞

Notons, pour conclure cette partie algorithmique, que plusieurs stratégies peuvent être
mises en œuvre pour l’échantillonnage du champ Gaussien {Ξ(xs ), xs ∈ DI }, parmi
lesquelles la méthode spectrale (développée dans [110, 111] dans le cas de champs ho-
mogènes ; voir [117] pour des développements complémentaires, notamment dans le cas
non homogène) ou une méthode de factorisation. Dans ce travail, la méthode basée sur
la factorisation de Cholesky est utilisée en raison du faible temps calcul associé (la
décomposition de Cholesky de la matrice de covariance étant calculée dans une phase
de prétraitement).

6 Prise en compte du bruit d’échantillonnage dans

les simulations par dynamique moléculaire

Dans cette section, nous nous intéressons à la prise en compte, dans la modélisation de
mécanique des milieux continus, du bruit d’échantillonnage présent dans les simulations
par Dynamique Moléculaire – ce dernier pouvant biaiser les résultats de l’identification
Modélisation stochastique 67

inverse. Plus spécifiquement, nous supposerons ici que le bruit ne modifie pas les lon-
gueurs de corrélation spatiales (qui sont induites par les modèles d’interaction entre en-
tités atomiques), mais peut avoir un impact sur l’estimation du niveau des fluctuations
statistiques. Pour ce faire, nous procédons de la façon suivante :
— dans le cas des zones de silice et de polymère pur (i.e. hors zone d’interphase),
les propriétés mécaniques élastiques sont modélisées par des matrices aléatoires
centrées sur des moyennes données (celles considérées sans la modélisation du bruit)
et dépendante d’un même paramètre mesurant le niveau du bruit ; dans ce cas, les
propriétés des phases sont aléatoires mais ne dépendent pas du point de l’espace
considéré ;
— dans la région de l’interphase, un terme additionnel indépendant du point, centré
sur la matrice identité et générant des fluctuations tricliniques est introduit dans
la décomposition algébrique du champ (exprimée en coordonnées sphériques) ; les
fluctuations de ce terme sont contrôlées par le paramètre de dispersion introduit
précédemment pour les zones de silice et de polymère pur.
Le modèle de matrices aléatoires utilisé est défini ci-dessous.

6.1 Modèle probabiliste pour des matrices aléatoires à fluctua-

tions tricliniques

Afin de modéliser la matrice aléatoire à fluctuations tricliniques, l’ensemble SE+ de ma-

trices aléatoires construit dans [112] pour la modélisation des incertitudes en dynamique
des structures est considéré. Soit [Y ] la matrice aléatoire à valeurs dans M+
6 (R), satisfai-
sant les équations de contrainte suivantes :

E{[Y ]} = [I6 ] , (3.29)

E{ln (det([Y ]))} = ν[Y ] , |ν[Y ] | < +∞ . (3.30)

On montre alors que la densité de probabilité de la matrice aléatoire [Y ] construite par

le principe du maximum d’entropie s’écrit

5 + 2`
p[Y ] ([Y ]) = IM+6 (R) ([Y ]) k det([Y ])`−1 exp{− tr([Y ])} , (3.31)
2
Modélisation stochastique 68

avec [Y ] 7→ IM+6 (R) ([Y ]) la fonction indicatrice de l’ensemble M+

6 (R), k la constante de
normalisation (dont l’expression explicite peut être trouvée dans [112]) et

(1 − δ[Y ] 2 ) 1 + δ[Y ] 2
`= n + . (3.32)
2δ[Y ] 2 2δ[Y ] 2

Ci-dessus, δ[Y ] est le paramètre mesurant le niveau des fluctuations statistiques de [Y ],

avec r
n+1
0 < δ[Y ] < . (3.33)
n+5
Dans ce qui suit, la notation [Y ] ∼ SE+ (δ[Y ] ) signifie que la matrice aléatoire [Y ] est
définie par la densité de probabilité donnée par l’Eq. (3.31) et dépendante du paramètre
de dispersion δ[Y ] .

6.2 Modélisation du bruit dans les différentes phases

Dans les zones associées à la silice et au polymère pur, les propriétés mécaniques sont
modélisées par les matrices aléatoires d’élasticité notées [C silice ] et [C pp ], respectivement.
Ces dernières sont définies à partir des résultats estimés par les simulations de dynamique
moléculaire selon

[C silice ] := [Lsilice ]T [Z silice ][Lsilice ] , [C pp ] := [Lpp ]T [Z pp ][Lpp ] , (3.34)

où les matrices [Lsilice ] et [Lpp ] sont issues des décompositions de Cholesky des valeurs
moyennes définies au chapitre précédent,

[C silice ] = [Lsilice ]T [Lsilice ] , [C pp ] = [Lpp ]T [Lpp ] , (3.35)

et
[Z silice ] ∼ SE+ (δb ) , [Z pp ] ∼ SE+ (δb ) . (3.36)

Par construction, les matrices aléatoires [C silice ] et [C pp ] exhibent donc des fluctuations
tricliniques autour des modèles moyens estimés par les simulations de DM. Dans l’in-
terphase, on introduit le nouveau champ aléatoire {[Cb s (xs )], xs ∈ DI }, défini par une
modification de l’expression algébrique de l’Eq. (3.5) :

s
∀xs ∈ DI , b (xs )] = η[I6 ] + [Ms (xs )]1/2 [Z D ] [Ms (xs )]1/2 ,
[C (3.37)
I
Modélisation stochastique 69

où {[Ms (xs )] , xs ∈ DI } est le champ aléatoire défini aux sections précédentes, et [Z DI ] ∼
SE+ (δb ). Le modèle de champ avec bruit dépend alors
— des paramètres δ[N ] , Lr , Lθ et Lϕ définissant le modèle du champ {[Ms (xs )] , xs ∈
DI } à valeurs dans E sym (n) ;
— du paramètre de dispersion mesurant le niveau de bruit, δb .

7 Conclusion

Dans ce chapitre, nous avons proposé la construction d’un modèle probabiliste pour
le champ aléatoire non Gaussien du tenseur d’élasticité dans la zone d’interphase. La
méthodologie générale consiste à écrire le champ comme une transformation non linéaire
de champs Gaussiens sous-jacents, la transformation étant telle que le champ ainsi défini
exhibe une famille cible de lois marginales d’ordre un. Cette dernière est spécifiée dans
le cadre de la Théorie de l’Information et plus précisément, en invoquant le principe du
Maximum d’Entropie. Le problème d’optimisation associé est formulé en intégrant des
contraintes mathématiques liées à la valeur moyenne, ainsi qu’à la finitude des moments
d’ordre deux pour le tenseur des rigidités et son inverse. La symétrie matérielle est no-
tamment inférée à partir des résultats des simulations de DM présentés au Chap. 2, qui
permettent également de définir une structure de corrélation ad hoc pour les champs Gaus-
siens sous-jacents. Le générateur de réalisations associé est ensuite présenté et fait interve-
nir une famille d’équations différentielles stochastiques d’Itô discrétisées, en chaque point
de l’espace, par un schéma de Störmer-Verlet. La modélisation du bruit d’échantillonnage
induit par les simulations de DM est enfin abordée. Pour ce faire, un modèle de matrices
aléatoires à fluctuations tricliniques est utilisé afin d’introduire un bruit anisotrope dans
la zone d’interphase, ainsi que pour la modélisation des propriétés dans le polymère pur et
la silice. Les modèles ainsi construits offrent l’avantage de ne dépendre que sur un nombre
limité de paramètres (typiquement, un paramètre de dispersion et des paramètres liés à la
corrélation spatiale pour le modèle de champ), ce qui s’avère particulièrement intéressant
en vue de leur identification. Ce point est l’objectif du chapitre suivant.
Chapitre 4

Identification des paramètres du

modèle

Ce chapitre est dédié à la calibration du modèle probabiliste associé au champ aléatoire

non-Gaussien des rigidités dans l’interphase (voir le Chap. 3). Pour ce faire, un problème
statistique inverse est formulé, et repose sur l’utilisation des résultats des essais mécaniques
virtuels simulés par dynamique moléculaire, décrits au Chap. 2. La stratégie générale
d’identification est tout d’abord présentée à la Sec. 1. La méthode est basée sur l’équivalen-
ce (postulée) des propriétés apparentes définies par les descriptions moléculaires et conti-
nues, et fait intervenir un problème d’homogénéisation sans séparation des échelles, for-
malisé à la Sec. 1.3. La définition des fonctions coût est ensuite détaillée à la Sec. 1.4. Les
résultats numériques sont alors présentés dans la Sec. 2. En particulier, l’incidence de la
modélisation du bruit induit par les simulations de dynamique moléculaire est discuté.

1 Méthodologie d’identification

1.1 Principe

D’après le Chap. 3, le modèle probabiliste pour le champ aléatoire {[C s (xs )], xs ∈ DI }
représentant le champ des rigidités dans la zone d’interphase DI (sans prise en compte
du bruit induit par les simulations de DM) s’écrit :

∀xs ∈ DI , [C s (xs )] = η[I6 ] + [Ms ]1/2 exp{[G s (xs )]} [Ms ]1/2 . (4.1)

71
Identification du modèle probabiliste 72

Par conséquent, le modèle stochastique dépend de deux vecteurs déterministes, notés

h1 et h2 respectivement, définis comme suit (on rappelle que η désigne un paramètre
arbitraire de régularisation, dont la valeur est fixée a priori ) :
— le vecteur d’hyperparamètres h1 := (m1 , . . . , m5 ) ∈ S1 ⊂ R5 rassemble les va-
leurs des coordonnées de la matrice moyenne [Ms ] décomposée, en coordonnées
sphériques, sur la base de Walpole {[E i ]}5i=1 (la définition de S1 , induite en parti-
culier par la définie-positivité de [Ms ] est précisée par la suite) ;
— le vecteur h2 = (δ[N ] , Lr , La ) ∈ S2 est composé des hyperparamètres liés au niveau
de fluctuation de la matrice aléatoire [N s (xs )] = exp{[G s (xs )]} au point xs (rappe-
lons que ce paramètre est indépendant de xs , de par les hypothèses considérées lors
de la construction du modèle), ainsi qu’aux paramètres des fonctions de corrélation
des champs Gaussiens sous-jacents – on a supposé ici, sur la base des résultats des
simulations de dynamique moléculaire, que Lθ = Lϕ =: La .
Dans ce travail, les valeurs optimales des vecteurs ci-dessus, désignées par hopt opt
1 et h2 res-
pectivement, sont identifiées en postulant l’équivalence des propriétés apparentes obtenues
à partir des simulations de dynamique moléculaire (avec npc = 10) et celles déterminées
dans une formulation de Mécanique des Milieux Continus (MMC), en prenant en compte
le champ aléatoire des propriétés mécaniques élastiques dans l’interphase. Cette stratégie
de calibration est motivée par le fait des réalisations du champ tensoriel de rigidité locale
ne peuvent pas être extraites de simulations par DM et transférées de manière licite dans
la formulation continue MMC, pour au moins deux raisons principales. Premièrement, la
définition et l’interprétation de propriétés locales de MMC basées sur celles extraites à
partir des résultats DM (voir e.g. [88, 101]) est encore une question ouverte (sur laquelle
quelques éléments ont été reportés au Chap. 2), notamment lorsque la matière est confinée
[67]. Deuxièmement, les conditions de chargement appliquées dans les essais virtuels de
dynamique moléculaire (dans lesquels les propriétés locales sont de plus tricliniques par
rapport au repère de référence de la boı̂te de simulation) sont très différentes de celles qui
sont considérées dans le cadre de la MMC, de sorte qu’imposer une correspondance, dans
un certain sens probabiliste, entre les deux champs locaux serait discutable. Nous propo-
sons donc de développer une méthode de calibration basée sur la résolution d’un problème
inverse statistique ne faisant intervenir que les propriétés macroscopiques (dans la termi-
nologie usuelle de la physique statistique) et mésoscopiques (au sens de la mécanique des
milieux continus).
Identification du modèle probabiliste 73

1.2 Description du modèle en mécanique des milieux continus

Dans cette section, la modélisation équivalente en MMC du domaine de simulation par

DM (avec une matrice polymère à 10 chaı̂nes) est présentée. Soit D le domaine ouvert
borné de R3 , de bord ∂D, tel que D = (]−LD /2, LD /2[)3 , avec LD = 6.68 (nm). Ci-dessous,
BO (r) désigne la boule de rayon r (en nm), centrée à l’origine O du repère cartésien
(O, e1 , e2 , e3 ). D’un point de vue géométrique, l’inclusion de silice est modélisée par la
boule Dsilice := BO (1.5), tandis que la région de l’interphase est représentée par le domaine
DI := D ∩ (BO (3.5) \ BO (1.5)). Si le diamètre de la boule BO (3.5) est légèrement plus
grand que la longueur du côté du domaine D, cette caractéristique n’a pas d’effet notable
sur les résultats numériques présentés ci-dessous. Le domaine occupé par le polymère pur
est donc DP := D \ (Dsilice ∪ DI ). Ces différents domaines sont représentés sur la Fig. 4.1.

Figure 4.1 – Représentation graphique des différents domaines considérés (inclusion

de silice, zone d’interphase et polymère pur).

Les propriétés déterministes de chaque phase isotrope (silice et polymère pur) sont ex-
traites des simulations de dynamique moléculaire, présentées au Chap. 2, et sont rappelées
ci-dessous (k et µ désignent le modules de compressibilité et de cisaillement, respective-
ment) en GPa :
— pour l’inclusion de silice, on a k silice = 36.8, µsilice = 30.5 ;
— dans le cas du polymère pur, k pp = 5.19 et µpp = 0.93.
Dans la région de l’interphase, les propriétés élastiques sont modélisés par le modèle
de champ aléatoire introduit au Chap. 3 et pour lequel la décomposition algébrique est
rappelée dans l’Eq. (4.1). Afin d’alléger les notations, on considère à partir de cette
section que le champ est indexé et exprimé dans le système de coordonnées cartésiennes
(on omettra donc l’exposant “s”).
Identification du modèle probabiliste 74

1.3 Problème d’homogénéisation pour le modèle continu

1.3.1 Formulation théorique

Afin de formuler le problème inverse statistique, introduisons à présent la procédure d’ho-

mogénéisation linéaire dans le cadre de la mécanique des milieux continus. D’un point de
vue théorique, l’exposé ci-dessous est largement emprunté à l’ouvrage de référence [15].
Notons toutefois que contrairement aux hypothèses classiques énoncées dans la théorie
de l’homogénéisation des milieux aléatoires, le domaine D n’est pas supposé être un vo-
lume élémentaire représentatif, de sorte que les propriétés homogénéisées sont qualifiées
de propriétés apparentes et dépendent des conditions aux limites appliquées (par opposi-
tion aux propriétés effectives) [52] (voir [85] pour un exposé détaillé, ainsi que [57] pour
une mise en œuvre numérique). Ces propriétés apparents présentent donc des fluctua-
tions statistiques non négligeables, qui seront précisément utilisées lors de la phase de
calibration du modèle. Afin d’assurer la cohérence avec les simulations de DM, dans les-
quelles la variable de contrôle macroscopique est le tenseur des pressions, le problème
d’homogénéisation est formulé sous des conditions de contraintes homogènes au contour.
Pour les valeurs sélectionnées de h1 ∈ S1 et h2 ∈ S2 et pour une réalisation donné du
champ aléatoire dans l’interphase, on introduit donc le problème aux limites suivant :

−div [σ(x)] = 0 , ∀x ∈ D , (4.2)

[σ(x)] n(x) = [Σ] n(x) , ∀x ∈ ∂D , (4.3)

où x 7→ [σ(x)] est le champ de contrainte locale, n(x) est le vecteur unitaire normal à
∂D au point x et [Σ] est le tenseur des contraintes macroscopiques. La loi constitutive
s’écrit :
[σ(x)] = JC(x)K : [(x)] , ∀x ∈ D , (4.4)

où JC(x)K est la représentation tensorielle d’ordre quatre du tenseur d’ordre deux [C(x)].
Ce dernier est défini par morceaux, déterministe pour x ∈ Dsilice ou x ∈ DP et stochas-
tique pour x ∈ DI . Afin de favoriser un exposé synthétique, la même notation est utilisée
pour [C(x)], indépendamment de sa nature (aléatoire ou déterministe). On montre qu’il
existe un champ de tenseurs de concentration des contraintes x 7→ JB(x)K, à valeurs dans
l’ensemble des tenseurs de quatrième ordre (avec symétries mineures), tel que

[σ(x)] = JB(x)K : [Σ] , ∀x ∈ D (4.5)

Identification du modèle probabiliste 75

et Z
1
JB(x)K dx = JIK , (4.6)
|D| D

où l’on rappelle que JIK désigne le tenseur identité symétrique du quatrième ordre. Par
construction, on a
Bijk` (x) = σijk` (x) , 1 6 i, j, k, ` 6 3 , (4.7)

où x 7→ [σ k` (x)] est le champ des tenseurs des contraintes solution du problème aux
limites précédents, dans lequel

1
Σk`
ij := (δik δj` + δi` δjk ) . (4.8)
2

Il s’en suit que le champ x 7→ JB(x)K peut être complétement reconstruit en résolvant
le problème aux limites définis par les Eqs. (4.2) et (4.5), pour six valeurs différentes
du tenseur des contraintes macroscopiques [Σk` ]. La réalisation du tenseur de souplesse
apparent JSK,
e associée à la réalisation du champ aléatoire dans l’interphase, est alors
définie comme Z
e = 1
JSK JS(x)K : JB(x)K dx , (4.9)
|D| D

où x 7→ JS(x)K est le champ du tenseur de souplesse local de quatrième ordre (JS(x)K =
JC(x)K−1 pour tout x dans D). Soit JCK e = JSK e −1 est le tenseur d’élasticité apparent
associé, et notons [Ce MMC (h1 , h2 )] sa représentation de tenseur du second ordre (où
la dépendance sur h1 et h2 est rendue explicite pour la suite de l’exposé). De même,
[C
eMMC (h1 )] désigne le tenseur apparent (dans les mêmes conditions aux limites) obtenu
en considérant une interphase avec des propriétés déterministes définies par le vecteur h1
en coordonnées sphériques.

1.3.2 Résolution numérique

Dans ce travail, le problème aux limites défini par les Eqs. (4.2) et (4.5) est résolu, pour
toute réalisation du champ aléatoire dans la région de l’interphase, par la méthode des
éléments finis. Le domaine D est discrétisé avec des éléments tétraèdre à 4 nœuds (avec
un point d’intégration de Gauss par élément), comme illustré sur la Fig. 4.2.

Le modèle numérique final est constitué de 190 310 éléments (correspondant à un total de
102 561 degrés de liberté). Il convient de noter qu’une telle densité du maillage fournit au
moins quatre points d’intégration par longueur de corrélation (en moyenne), quelle que
Identification du modèle probabiliste 76

Figure 4.2 – Visualisation du maillage pour le modèle continu (l’inclusion apparaı̂t

en magenta, l’interphase en blanc et la matrice de polymère pur en cyan).

soit la direction ou la configuration testée – assurant ainsi un échantillonnage satisfaisant

de la structure de corrélation du champ aléatoire.

1.4 Définition des fonctions coût

Afin de réduire le coût calcul associé à l’exploration des espaces admissibles S1 et S2 ,

nous proposons ci-après de procéder de façon séquentielle, en deux étapes :
— dans un premier temps, la valeur de hopt
1 est recherchée, sans considérer de fluctua-
tions statistiques dans la zone d’interphase – la distance retenue dans la fonction
coût sera donc déterministe ;
— dans un second temps, la valeur optimale du second hyperparamètre h2 à va-
p
leurs dans S2 =]0, 7/11[×[0, +∞[×[0, +∞[ est déterminée à l’aide du principe du
maximum de vraisemblance, en considérant la valeur h1 = hopt
1 trouvée à l’étape
précédente.
La définition précise de chaque fonction coût est détaillée dans les sections suivantes. Du
point de vue des notations, on rappelle ici que le tenseur apparent des rigidités défini au
sens des essais virtuels par dynamique moléculaire (voir la Sec. 4) est noté [C
e DM ], et que
eMD (ωi )]}20 .
l’ensemble des réalisations obtenues numériquement est désigné par {[C i=1

1.4.1 Calibration du modèle moyen

Suivant les notations pour le tenseur apparent définies à la Sec. 1.3.1, le vecteur hopt
1 des
hyperparamètres optimaux pour la valeur moyenne [Ms ] du champ aléatoire {[C s (xs )], xs ∈
Identification du modèle probabiliste 77

DI } (exprimé en coordonnées sphériques ici) est définie comme

hopt
1 := arg min L1 (h1 ) , (4.10)
h1 ∈ S1

où L1 est la fonction coût

kP iso (E{[C
e DM ]}) − P iso ([C
eMMC (h1 )])kF
L1 (h1 ) := , (4.11)
kP iso (E{[C
e DM ]})kF

et le sous-ensemble S1 de R5 est tel que la matrice moyenne est définie-positive, c’est-à-

dire :

S1 = {m ∈ R5 | m1 > 0 , m2 > 0, m4 > 0, m5 > 0, m1 m2 − m23 > 0} . (4.12)

Dans l’Eq. (4.11), P iso est l’opérateur de projection sur l’ensemble des tenseurs isotropes,
introduit à la Sec. 4.4 du Chap. 2.

1.4.2 Calibration du niveau de fluctuation statistique et des paramètres des

fonctions de corrélation

Comme mentionné précédemment, l’identification du vecteur h2 est accomplie en fixant

la valeur de h1 à la valeur hopt
1 définie par l’Eq. (4.10), et en invoquant le principe du
maximum de vraisemblance. Il vient donc :

hopt
2 := arg max L2 (h2 ) , (4.13)
h2 ∈ S2

p
où S2 =]0, 7/11[×[0, +∞[×[0, +∞[, L2 est la fonction de vraisemblance donnée par

20
Y
L2 (h2 ) := p[Ce MMC (hopt ,h2 )] [C
eDM (ωi )] (4.14)
1
i=1

e MMC (hopt , h2 )]. En pratique,

et p[Ce MMC (hopt ,h2 )] est la densité de probabilité estimée de [C
1
1
le nombre de réalisations calculées pour le tenseur apparent ne permet pas d’obtenir un
estimateur par noyau multidimensionnel convergé. Pour cette raison, on substitue à la
définition précédente de L2 la forme suivante (la notation n’est pas modifiée ici, par souci
Identification du modèle probabiliste 78

de concision), basée sur l’estimation de densités marginales d’ordre un :

20
!
Y Y
L2 (h2 ) := pCeMMC (hopt ,h2 )jk C
eDM (ωi )jk . (4.15)
1
i=1 16j6k66

Notons enfin que d’un point de vue numérique, la recherche de maximum est réalisée de
façon équivalente sur la fonction h2 7→ ln (L2 (h2 )).

1.4.3 Remarque sur l’échantillonnage du champ aléatoire par la méthode

ergodique

Soit t∗ := (k ∗ − 1)∆t le premier instant dans le régime stationnaire, avec ∆t le pas de

discrétisation du schéma utilisé lors de la résolution des équations différentielles stochas-
tiques (voir la Sec. 5 du Chap. 3). Pour 1 6 i 6 j 6 5, un indice de décorrélation τijε (xs )
associé aux processus stochastiques {Ui (t, xs ), t > t∗ } et {Uj (t, xs ), t > t∗ } (où xs est un
point quelconque de DI ) est introduit et tel que

|E Uj ((k + τij (xs ))∆t, xs )Ui (k∆t, xs ) | 6 ε

(4.16)

pour un certain seuil de décorrélation arbitraire 0 < ε 1, avec k > k ∗ . Un indice de

décorrélation global τ pour le processus stochastique {U (t, xs ), t > t∗ } peut donc être
défini comme

τ = sup max τijε (xs ) . (4.17)
xs ∈ DI 16i6j 65

En pratique, pour une valeur suffisamment petite de ε 1, ce résultat signifie que

{U ((k + τ ε )∆t, xs ), xs ∈ DI } et {U (k∆t, xs ), xs ∈ DI }, k > k ∗ sont des réalisations fai-
blement corrélées du champ aléatoire {Gs (xs ), xs ∈ DI }. Cette stratégie d’échantillonnage
est particulièrement pertinente lorsque la vitesse de convergence vers la solution station-
naire (pour tous les points de DI ) est lente et nécessite un nombre d’itérations beau-
coup plus grand que τ ε (de sorte que l’échantillonnage par la méthode de Monte Carlo
serait moins efficace en terme de temps calcul). La valeur du paramètre ε doit donc
être soigneusement sélectionnée afin d’obtenir un compromis acceptable entre la qualité
d’échantillonnage et le coût de calcul.
Identification du modèle probabiliste 79

2 Résultats

2.1 Identification du modèle moyen

Le problème d’optimisation défini par les Eqs. (4.10–4.11) est résolu en utilisant la fonction
fmincon de Matlab, sous les contraintes d’inégalité issue de la définition de l’espace ad-
missible S1 (voir l’Eq. (4.12)). La convergence de l’algorithme de recherche est représentée
sur la Fig. 4.3, où minL1 (k) désigne (avec un abus de notation) la valeur minimale de la
fonction coût L1 déterminée à l’itération k.

0
10

¹2
10
min L 1 (k )

¹4
10

¹6
10

¹8
10
0 5 10 15 20 25 30 35 40
Itération k

Figure 4.3 – Convergence de l’algorithme d’optimisation pour le problème défini par

les Eqs. (4.10) et (4.11) : graphique de la fonction k 7→ min L1 (k).

On observe une convergence rapide de l’algorithme, avec une valeur inférieure de la fonc-
tion coût égale à 3.3 × 10−8 après 40 itérations. En notant que hopt opt opt
1 = (m1 , . . . , m5 ) et
que [Ms ] = 5i=1 mopt i
P
i [E ] (en coordonnées sphériques), il vient (en GPa) :

 
5.51 4.02 4.02 0 0 0
 

 5.19 1.16 0 0 0 

 
s
 5.19 0 0 0 
[M ] =   . (4.18)

 4.03 0 0 

 

 Sym. 4.25 0 

4.25
Identification du modèle probabiliste 80

Un effet de renforcement en cisaillement est clairement observé pour le modèle moyen dans
la région de l’interphase, ce qui est en accord avec la perte de la mobilité du polymère
dans le voisinage de l’inclusion, ainsi qu’avec résultats numériques obtenus par ailleurs
(voir par exemple [99]).

2.2 Calibration du niveau de fluctuation statistique et de la

structure de corrélation (sans bruit)

2.2.1 Choix des paramètres pour l’intégration des équations différentielles

stochastiques

Pour x fixé dans DI , l’équation différentielle stochastique d’Itô définie par l’Eq. (3.24) est
intégrée avec % = 9 et en utilisant un pas ∆t égal à 0.01 – ces valeurs ayant été estimées à
partir d’une analyse de convergence. Plus spécifiquement, la convergence vers la solution
stationnaire (pour x ∈ DI fixé arbitrairement) est caractérisée par la convergence la
fonction Niter 7→ Conv (Niter ) telle que :

Niter
1 X 2
Conv (Niter ) := Uk (x) , Niter > 1 . (4.19)
Niter k=1

Le graphique de cette fonction est représenté sur la Fig. 4.4.

1.2
1.1
1
0.9
Conv(Niter)

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0 2 4 6 8 10
Niter 4
x 10

Figure 4.4 – Graphique de la fonction Niter 7→ Conv (Niter ) (% = 9, ∆t = 0.01).

Identification du modèle probabiliste 81

On observe une convergence raisonnable après 20 000 itérations. On pose donc k ∗ =

20 000. Par ailleurs, le paramètre ε impliqué dans l’échantillonnage ergodique est choisi
égal à 0.1. Le graphique de la fonction

τ11 7→ Decorr(τ11 ) := |E{U1 ((k + τ11 )∆t, x)U1 (k∆t, x)}|

(avec k > k ∗ et x fixé arbitrairement ; la dépendance de τ11 par rapport à x n’est pas
reportée afin de simplifier les notations), où l’espérance mathématique est estimée de
manière ergodique, est représenté sur la Fig. 4.5.

0.9

0.8

0.7
Decorr(τ 11)

0.6

0.5

0.4

0.3

0.2

0.1

0
0 200 400 600 800 1000
τ 11

Figure 4.5 – Graphique de la fonction τ11 7→ Decorr(τ11 ).

ε
On constate donc que τ11 = 200 pour ε = 0.1. En pratique, cette valeur correspond au
maximum τ ε recherché (voir l’Eq. (4.17)), c’est-à-dire que dans le régime stationnaire (k >
20 000), une réalisation indépendante du champ aléatoire {G(x), x ∈ DI } est obtenue
toutes les 200 pas de temps. Enfin, un ensemble de 300 réalisations indépendantes est
généré et utilisé par la suite afin d’estimer, pour des valeurs données des hyperparamètres,
les différentes densités de probabilité considérées.

2.2.2 Résolution

Le problème d’optimisation non convexe défini par les Eqs. (4.13) et (4.15) est résolu par
un échantillonnage régulier des hyperparamètres dans leurs espaces admissibles respectifs.
Afin d’interpréter plus facilement les résultats d’un point de vue physique, la longueur
de corrélation dans la direction radiale est normalisée par l’épaisseur de l’interphase eI ,
Identification du modèle probabiliste 82

tandis que les paramètres liés à la (dé)corrélation angulaire sont normalisées par le terme
2πRext , avec Rext := Rp + eI . Le graph de la fonction de coût L2 est représenté sur la
Fig. 4.6 pour différentes combinaisons des paramètres. Les valeurs optimales sont données

−30
0.5
−40
0.4 −50
δ[N ]

0.3 −60
0.2 −70
0.1 −80
20 4
15 −90
10 3
5 2 −100
2πRext /La 1 eI /Lr

Figure 4.6 – Graphique 4D de la fonction de vraisemblance h2 7→ L2 (h2 ) pour

différentes combinaisons de δ[N ] , Lr et La . La valeur de la fonction est proportion-
nelle à la taille du marqueur, ainsi qu’à l’intensité de la couleur (les valeurs infinies
de L2 sont remplacées par la valeur arbitraire -100, afin de visualiser l’ensemble des
résultats).

ci-dessous :

opt
δ[N ] = 0.3 , Lopt
r ≈ 0.9eI = 1.8 (nm) , Lopt
a ≈ 3.85 rad × nm . (4.20)

opt
La valeur du paramètre δ[N ] indique que les fluctuations statistiques dans l’interphase
sont importantes et ne peuvent pas être totalement annihilées par une procédure d’ho-
mogénéisation (puisque la condition Lopt
r eI n’est pas satisfaite, par exemple). La
valeur optimale du paramètre La montre en outre que, pour une position radiale donnée,
le champ aléatoire d’élasticité est corrélé sur une demi-sphère et donc, par symétrie, sur
la sphère. D’un point de vue qualitatif, ce résultat est cohérent avec les résultats obtenus
par les simulations de dynamique moléculaire et notamment, avec la conformation locale
des chaı̂nes autour de l’inclusion de silice – voir l’illustration sur la Fig. 4.7. Les densités
de probabilité de certaines composantes, ainsi que les réalisations obtenues par les simula-
tions de DM, sont représentées sur les Figs. 4.8, 4.9 et 4.10. Enfin, et à titre d’illustration,
int
une réalisation du champ aléatoire {C11 (xs ), xs ∈ ∂DIext } identifié (où ∂DIext représente
la surface extérieure de DI ) est représentée sur la Fig. 4.11.
Identification du modèle probabiliste 83

Figure 4.7 – Configuration instantanée des chaı̂nes polymères piégées dans des co-
quilles sphériques d’une épaisseur de 5 Å dans la région de l’interphase. Sur la figure de
gauche (resp. droite), la coquille correspond à la coquille la plus proche (resp. éloignée)
de la surface de l’inclusion.

1
Simulations DM
0.9
Simulations MMC
0.8
Densité de probabilité

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0
4 4.5 5 5.5 6 6.5 7 7.5 8
C 11 (GPa)

Figure 4.8 – Densité de probabilité de la composante (1, 1) du tenseur apparent (trait

noir). Les réalisations associées, obtenues par les simulations de DM, sont représentées
par les étoiles rouges.

2.3 Calibration avec bruit d’échantillonnage de DM

Dans cette section, l’effet du bruit d’échantillonnage induit par les simulations de dy-
namique moléculairesur les résultats d’identification est étudié. D’un point de vue phy-
sique, il est raisonnable de supposer qu’un tel bruit n’a pas d’impact sur la structure de
corrélation du champ aléatoire d’élasticité dans la région interphase, mais peut conduire à
Identification du modèle probabiliste 84

1
Simulations DM
0.9
Simulations MMC
0.8

Densité de probabilité
0.7

0.6

0.5

0.4

0.3

0.2

0.1

0
1.5 2 2.5 3 3.5 4 4.5 5 5.5
C 12 (GPa)

Figure 4.9 – Densité de probabilité de la composante (1, 2) du tenseur apparent (trait

noir). Les réalisations associées, obtenues par les simulations de DM, sont représentées
par les étoiles rouges.

5.5
Simulations DM
5
Simulations MMC
4.5
4
Densité de probabilité

3.5
3
2.5
2
1.5
1
0.5
0
1.8 2 2.2 2.4 2.6 2.8
C 44 (GPa)

Figure 4.10 – Densité de probabilité de la composante (4, 4) du tenseur apparent (trait

noir). Les réalisations associées, obtenues par les simulations de DM, sont représentées
par les étoiles rouges.

une surestimation du δ[N ] mesurant le niveau des fluctuations statistiques. Par conséquent,
ce bruit est pris en compte ci-dessous :
— en modélisant les propriétés (constantes en espace) du polymère pur et de la silice
par des matrices aléatoires ;
— en considérant un bruit triclinique, indépendant de xs , dans l’interphase où le champ
aléatoire d’élasticité est écrit en coordonnées sphériques.
Identification du modèle probabiliste 85

Figure 4.11 – Visualisation d’une réalisation de l’élasticité du champ aléatoire

int (xs ), xs ∈ ∂D ext } (en GPa et en coordonnées sphériques) sur la surface extérieure
{C11 I
de la région de l’interphase (avec hyperparamètres calibrées).

Pour ce faire, les modèles probabilistes introduits à la Sec. 6.2 du Chap. 3 sont uti-
lisés. Les propriétés élastiques dans la silice et dans le polymère pur sont donc définies
par les Eqs. (3.34), (3.35) et (3.36), dans lesquelles on suppose que le paramètre de
dispersion est identique – et noté δbruit . Le nouveau champ aléatoire dans l’interphase,
présentant des fluctuations statistiques tricliniques autour de la classe isotrope trans-
verse (en coordonnées sphériques), est noté {[Cb s (xs )], xs ∈ DI }. Ce champ est défini
par l’Eq. (3.37), dans laquelle on suppose également que [Z DI ] ∼ SE+ (δbruit ). En notant
[C
e MMC (δ[N ] , δbruit )] la variable aléatoire modélisant le tenseur de rigidité apparent ob-
tenu par l’homogénéisation de la microstructure aléatoire ainsi introduite (en considérant
les valeurs de hopt
1 et des paramètres Lopt
r et Lopt
a identifiées à la Sec. 2.2.2), les valeurs
optimales pour les paramètres δ[N ] and δbruit sont déterminées en maximisant la fonction

20
!
Y Y
L3 δ[N ] , δbruit = pCeMMC (δ[N ] ,δbruit )jk C
eDM (ωi )jk , (4.21)
i=1 16j6k66

p p
sur l’espace [0, 7/11[×[0, 7/11[. Le graph de la fonction de vraisemblance δ[N ] , δbruit 7→

L3 δ[N ] , δbruit est représenté sur la Fig. 4.12.
opt opt
Les valeurs optimales sont égales à δ[N ] = 0.2 et δbruit = 0.1 respectivement, ce qui montre
que le fait de négliger le bruit d’échantillonnage pour les résultats obtenus par DM induit
une surestimation de δ[N ] .
Identification du modèle probabiliste 86

"
L 3 δ [N ], δ bruit
¹ 50

¹ 100
!

¹ 150
0.5
0.4 0
0.05
0.1
0.3 0.15
0.2
0.2 0.25
0.3
0.35
0.1 0.4 δ bruit
δ [N ]

Figure 4.12 – Graph de la fonction de vraisemblance δ[N ] , δbruit 7→ L3 δ[N ] , δbruit .
Le résultat sans la modélisation du bruit correspond à la courbe δ[N ] 7→ L3 δ[N ] , 0 .

3 Conclusion

Dans ce dernier chapitre, nous avons présenté une méthodologie d’identification des
modèles probabilistes construits au Chap. 3. En ce qui concerne le champ aléatoire du
tenseur des rigidités dans l’interphase, l’identification du modèle moyen, du paramètre
de dispersion et des deux paramètres contrôlant la corrélation du champ dans la di-
rection radiale et dans les directions orthoradiales est accomplie de façon séquentielle,
par l’intermédiaire de deux problèmes inverses. De façon générale, l’idée consiste à mi-
nimiser l’écart observé entre les résultats estimés par les essais virtuels en dynamique
moléculaire (voir le Chap. 2) et ceux obtenus par une procédure d’homogénéisation
numérique en mécanique des milieux continus (sous des conditions de contraintes ho-
mogènes au contour). Le modèle moyen est donc calibré dans un premier temps en im-
posant, dans un cadre déterministe, l’équivalence du tenseur homogénéisé du domaine
dans lequel les propriétés de l’interphase correspondent au modèle moyen recherché, et la
moyenne des réalisations obtenues par dynamique moléculaire. Dans un second temps, le
paramètre de dispersion et les paramètres liés à la corrélation spatiale sont identifiés par
l’intermédiaire du principe du maximum de vraisemblance, formulé sur les composantes
du tenseur apparent. On montre que la longueur de corrélation dans la direction radiale
est proche de l’épaisseur de l’interphase, tandis que le champ est corrélé, dans les direc-
tions orthoradiales, sur une hémisphère – en accord qualitatif avec la morphologie locale
des chaı̂nes proches de la surface de l’inclusion. Ces résultats montrent en particulier que
Identification du modèle probabiliste 87

les propriétés dans la zone d’interphase ne peuvent pas être homogénéisées au sens clas-
sique du terme, puisque la condition de séparation des échelles n’est pas satisfaite dans la
direction radiale, par exemple. Ces résultats statistiques sont enfin raffinés par l’identifi-
cation du paramètre de dispersion associé au bruit d’échantillonnage, ce qui permet une
meilleure estimation du paramètre contrôlant le niveau des fluctuations du champ dans
l’interphase.
Chapitre 5

Conclusions et perspectives

1 Conclusions générales

Ce travail a été consacré à la modélisation probabiliste et à l’identification, dans un cadre

de mécanique des milieux continus, des propriétés d’élasticité aléatoires dans la zone
d’interphase présente au voisinage des hétérogénéités dans un nanocomposite modèle,
constitué d’une matrice polymère renforcée par une inclusion de silice.

Dans un premier temps, des simulations par dynamique moléculaire ont été conduites,
dans la lignée des travaux présentés dans [18], afin de caractériser l’interphase, tant
d’un point de vue géométrique que du comportement mécanique élastique. En accord
avec la littérature, on constate que la zone d’interphase peut être modélisée par une
coquille sphérique d’épaisseur constante, indépendante de la taille de l’inclusion (mais
dépendante de la température). Les simulations moléculaires montrent par ailleurs l’orien-
tation préférentielle des chaı̂nes du polymère proches de la surface de l’inclusion, si bien
que le tenseur d’élasticité dans la zone perturbée peut être modélisé comme un tenseur ex-
hibant une isotropie transverse radiale en coordonnées sphériques. Des essais mécaniques
virtuels ont par la suite été réalisés et permettent d’estimer des réalisations indépendantes
du tenseur d’élasticité apparent associé au domaine de simulation moléculaire. Ce tenseur
présente des fluctuations statistiques non négligeables, induites d’une part par les varia-
tions de rugosité des réalisations de la particule de silice (de faible diamètre), et d’autre
part par le bruit intrinsèque généré par les calculs de DM.

Dans une seconde partie de ce travail, la construction d’un modèle probabiliste pour le
champ aléatoire des rigidités dans l’interphase a été abordée dans le cadre de la Théorie

89
Conclusions et perspectives 90

de l’Information. La méthodologie générale consiste à écrire le champ comme une trans-

formation non linéaire de champs Gaussiens sous-jacents, la transformation étant telle
que le champ ainsi défini exhibe une loi marginale d’ordre 1 cible, indépendante de la
position (d’après les résultats obtenus par DM). La définition de cette dernière repose
sur une décomposition algébrique ad hoc, ainsi que sur l’utilisation d’une transforma-
tion exponentielle, ce qui permet de formuler le problème de maximisation d’entropie
dans un espace non contraint (en l’occurrence, R5 ). Le générateur de réalisations pour
ce modèle non standard a ensuite été présenté. La stratégie consiste à construire une
famille d’équations différentielles stochastiques d’Itô indexée en espace (et dont chaque
élément est discrétisé par un schéma de Störmer-Verlet), telle que la famille de mesures
invariantes associées coincident par construction avec la famille de lois marginales d’ordre
1 à prescrire. La modélisation du bruit d’échantillonnage induit par les simulations de DM
a ensuite été accomplie par l’intermédiaire d’un modèle de matrices aléatoires à fluctua-
tions tricliniques. Notons ici que ce choix de modélisation est introduit, dans le cadre de
ces travaux, afin de pouvoir dissocier les fluctuations intrinsèques du système de celles
générées par les simulations de DM – l’objectif ici étant de mieux estimer le paramètre
de dispersion du champ dans l’interphase.

La dernière partie du mémoire a été dédiée au développement d’une méthodologie de

calibration pour les modélisations probabilistes introduites. La stratégie proposée dans ce
cadre repose sur la résolution séquentielle de deux problèmes inverses, l’un lié à l’identi-
fication du modèle, l’autre à celle des paramètres de dispersion et de corrélation spatiale.
Conceptuellement, l’idée consiste à minimiser (dans un certain sens) l’écart observé entre
les résultats estimés par les essais virtuels en dynamique moléculaire et ceux obtenus
par une procédure d’homogénéisation numérique en mécanique des milieux continus –
sous des conditions de contraintes homogènes au contour. On montre que la longueur
de corrélation calibrée dans la direction radiale est proche de l’épaisseur de l’interphase,
tandis que le champ d’élasticité est corrélé, dans les directions orthoradiales, sur une
hémisphère (ce qui se révèle en accord avec la morphologie locale des chaı̂nes proches de
la surface de l’inclusion, d’après les simulations de DM). D’un point de vue multi-échelle,
ces résultats indiquent par ailleurs que pour la microstructure considérée, les propriétés
dans la zone d’interphase ne peuvent pas être homogénéisées au sens classique du terme
– puisque la condition de séparation des échelles n’est pas satisfaite dans la direction
radiale. En combinant la représentation de champ avec un modèle de matrices aléatoires
associé au bruit intrinsèque des simulations de DM, on observe une diminution de la va-
leur du paramètre de dispersion dans la zone d’interphase, soulignant ainsi l’importance
de la prise de ce bruit dans la méthodologie d’identification inverse.
Conclusions et perspectives 91

2 Perspectives

Ce travail a concerné le développement et l’identification, sur la base de simulations de

DM, d’une représentation probabiliste en champ pour le tenseur d’élasticité dans la zone
d’interphase. Plusieurs pistes de développements futurs peuvent en être dégagées. Tout
d’abord, il serait intéressant de revisiter les résultats d’identification
— sur la base d’une description plus riche du comportement du polymère en général,
et dans l’interphase en particulier – notamment pour une température proche de
l’ambiante, pour laquelle les phénomènes physiques (e.g. la viscosité) peuvent être
plus complexes à simuler, notamment au travers de calculs de DM ;
— en intégrant la modélisation probabiliste proposée dans une démarche plus générale,
dans laquelle la représentation est étendue par l’intermédiaire d’un développement
sur les chaos polynomiaux (voir [82, 116]).
Une modification de la température induisant par ailleurs un changement de compor-
tement des chaı̂nes du polymère dans l’interphase (dont l’épaisseur varie alors), la ca-
ractérisation de la symétrie matérielle exhibée et de l’évolution des paramètres de corréla-
tion lors de ces changements de température est une question ouverte que les méthodologies
développées pourraient permettre d’aborder. Une problématique sous-jacente ici serait
alors de clarifier dans quel cadre les propriétés de l’interphase peuvent être homogénéisées
(c’est-à-dire telles que le résidu statistique après le changement d’échelle soit négligeable).
Le développement de méthodologies permettant l’extraction directe de propriétés méso-
scopiques continues à partir de simulations de DM et pour des conditions aux limites
cohérentes pour le mécanicien des milieux continus constituerait une contribution intéres-
sante pour le couplage entre les deux types de représentation. Enfin, l’étude de la condi-
tion de séparation des échelles (dans le cadre de la mécanique des milieux continus)
dans le cas d’une microstructure contenant plusieurs inclusions entourées d’interphases
modélisées par des champs stochastiques, ainsi que l’identification inverse d’un modèle
d’interface équivalent, serait une contribution intéressante pour l’analyse micromécanique
des matériaux nanorenforcés.
Bibliographie

[1] N.C. Admal and E. B. Tadmor. A unified interpretation of stress in molecular

systems. Journal of Elasticity, 100 :63–143, 2010.

[2] A. Adnan, C. T. Sun, and H. Mahfuz. A molecular dynamics simulation study to

investigate the effect of filler size on elastic properties of polymer nanocomposites.
Composites Science and Technology, 67(3-4) :348–356, 2007.

[3] N. D. Albérola, K. Benzarti, C. Bas, and Y. Bomal. Interface effects in elastomers

reinforced by modified precipitated silica. Polymer Composites, 22(2) :312–325,
2001.

[4] B. J. Alder and T. E. Wainwright. Phase transition for a hard sphere system.
Journal of Chemical Physics, 27 :1208–1209, 1957.

[5] B. J. Alder and T. E. Wainwright. Studies in molecular dynamics. i. general method.

Journal of Chemical Physics, 31 :459–466, 1959.

[6] M. P. Allen and D. J. Tildesley. Computer Simulation of Liquids. Clarendon Press

Oxford, 1989.

[7] H.C. Andersen. Molecular dynamics simulations at constant pressure and/or tem-
perature. The Journal of Chemical Physics, 72 :2384–2393, 1980.

[8] S.R. Bakshi, D. Lahiri, and A. Argawal. Carbon nanotube reinforced metal matrix
composites - a review. International Materials Reviews, 55, 2010.

[9] D. Barbier, D. Brown, A. C. Grillet, and S. Neyertz. Interface between end-

functionalized peo oligomers and a silica nanoparticle studied by molecular dy-
namics simulations. Macromolecules, 37(12) :4695–4710, 2004.

[10] Y. Benveniste. A general interface model for a three-dimensional curved thin ani-
sotropic interphase bandween two anisotropic media. Journal of the Mechanics and
Physics of Solids, 54 :708–734, 2006.
93
Bibliographie 94

[11] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R.

Haak. Molecular-dynamics with coupling to an external bath. Journal of Chemical
Physics, 81 (8) :3684–3690, 1984.

[12] J. Berriot, F. Lequeux, L. Monnerie, H. Montes, D. Long, and P. Sotta. Filler-

elastomer interaction in model filled rubbers, a h-1 nmr study. Journal of Non-
Crystalline Solids, 307 :719–724, 2002.

[13] J. Berriot, F. Martin, H. Montes, L. Monnerie, and P. Sotta. Reinforcement of model

filled elastomers : characterization of the crosslinking density at the filler-elastomer
interface by h-1 nmr measurements. Polymer, 44(5) :1437–1447, 2003.

[14] K. Binder. Monte Carlo and Molecular Dynamics Simulations in Polymer Science.
Oxford University Press, 1995.

[15] M. Bornet, T. Bretheau, and P. Gilormini. Homogénéisation en mécanique des

matériaux 1 : matériaux aléatoires élastiques et milieux périodiques. Hermes
Science, 2001.

[16] S. Brisard, L. Dormieux, and D. Kondo. Hashin-shtrikman bounds on the bulk

modulus of a nanocomposite with spherical inclusions and interface effects. Com-
putational Materials Science, 48 :589 – 596, 2010.

[17] D. Brown, P. Mélé, S. Marceau, and N.D. Albérola. A molecular dynamics study
of a model nanoparticle embedded in a polymer matrix. Macromolecules, 36(4) :
1395–1406, 2003.

[18] D. Brown, V. Marcadon, P. Mélé, and N. D. Albérola. Effect of filler particle size on
the properties of model nanocomposites. Macromolecules, 41(4) :1499–1511, 2008.

[19] A. Bródka and T. W. Zerda. Properties of liquid acetone in silica pores : Molecular
dynamics simulation. The Journal of Chemical Physics, 104 :6319, 1996.

[20] K. Burrage, I. Lenane, and G. Lythe. Numerical methods for second-order stochastic
differential equations. Siam Journal on Scientific Computing, 29(1) :245–264, 2007.

[21] R.H. Cameron and W.T. Martin. The orthogonal development of non-linear func-
tionals in series of fourier-hermite functionals. Annals of Mathematics, 48 :385–392,
1947.

[22] V. Carnevale, C. Micheletti, F. Pontiggia, and R. Potestio. Multiscale Approaches

to Protein Modeling. Springer New York, 2011.
Bibliographie 95

[23] E. Chabert. Propriétés mécaniques de nanocomposites à matrice polymère : ap-

proche expérimentale et modélisation. PhD thesis, INSA de Lyon, 2002.

[24] S. Chang, S. Yang, H. Shin, and M. Cho. Multiscale homogenization model for ther-
moelastic behavior of epoxy-based composites with polydisperse sic nanoparticles.
Composite Structures, 128 :342–353, 2015.

[25] J. Cho and C. T. Sun. A molecular dynamics simulation study of inclusion size
effect on polymeric nanocomposites. Computational Materials Science, 41(1) :54–
62, 2007.

[26] J.M.G. Cowie and V. Arrighi. Polymers : Chemistry and Physics of Modern Mate-
rials. CRC Press, 2007.

[27] S. C. Cowin and M. M. Mehrabadi. The structure of the linear anisotropic elastic
symmetries. Journal of the Mechanics and Physics of Solids, 40(7) :1459–1471,
1992.

[28] J. Douce, J.-P. Boilot, J. Biteau, L. Scodellaro, and A. Jimenez. Effect of filler size
and surface condition of nano-sized silica particles in polysiloxane coatings. Thin
Solid Films, 1-2 :114–122, 2004.

[29] H.L Duan, J Wang, Z.P Huang, and B.L Karihaloo. Eshelby formalism for nano-
inhomogeneities. Proceedings of the Royal Society of London A : Mathematical,
Physical and Engineering Sciences, 461(2062) :3335–3353, 2005.

[30] Qi Chen et al. Nanoscale and effective mechanical behavior and fracture of silica
nanocomposites. Composites Science and Technology, 68 :3137–3144, 2008.

[31] F.I. Fedorov. Theory of elastic waves in crystals. Plenum Press, 1968.

[32] J.T. Fern, D.J. Keffer, and W.V. Steele. Vapor-liquid equilibrium of ethanol by
molecular dynamics simulation and voronoi tessellation. J. Phys. Chem. B, 111(46) :
13278–13286, 2007.

[33] Capaldi. F.M., M.C. Boyce, and G.C. Rutledge. Molecular response of a glassy
polymer to active deformation. Polymer, 45(4) :1391–1399, 2003.

[34] D. Fragiadakis. Glass transition and molecular dynamics in

poly(dimethylsiloxane)/silica nanocomposites. Polymer, 46 :6001–6008, 2005.

[35] D. Frenkel. Understanding Molecular Simulation, Second Edition : From Algorithms

to Applications (Computational Science). Academic Press, 2002.
Bibliographie 96

[36] S. Y. Fu, X. Q. Feng, B. Lauke, and Y. W. Mai. Effects of particle size, particle ma-
trix interface adhesion and particle loading on mechanical properties of particulate
polymer composites. Composites Part B Engineering, 39(6) :933–961, 2008.

[37] A. Ghanbari, M. C. Böhm, and F. Müller-Plathe. A simple reverse mapping pro-

cedure for coarse-grained polymer models with rigid side groups. Macromolecules,
44(13) :5520–5526, 2011.

[38] A. Ghanbari, T. V. M. Ndoro, F. Leroy, M. Rahimi, M. C. Böhm, and F. Müller-

Plathe. Interphase structure in silica-polystyrene nanocomposites : A coarse-grained
molecular dynamics study. Macromolecules, 45(1) :572–584, 2012.

[39] R. Ghanem and P.D. Spanos. Stochastic Finite Elements : A Spectral Approach.
Springer-Verlag, 1991.

[40] H. Ghasemi, R. Rafiee, X. Y. Zhuang, J. Muthu, and T. Rabczuk. Uncertainties

propagation in metamodel-based probabilistic optimization of cnt/polymer com-
posite structure using stochastic multi-scale modeling. Computational Materials
Science, 85 :295–305, 2014.

[41] J. Guilleminot and C. Soize. Stochastic model and generator for random fields with
symmetry properties : Application to the mesoscopic modeling of elastic random
media. Multiscale Modeling & Simulation, 11(3) :840–870, 2013.

[42] M. Gurtin and A. Murdoch. A continuum theory of elastic material surfaces. Archive
for Rational Mechanics and Analysis, 57 :291–323, 1975.

[43] M. Gurtin, J. Weissmuller, and F. Larche. A general theory of curved deformable

interfaces in solids at equilibrium. Philosophical Magazine A, 78 :1093–1109, 1998.

[44] E. Hairer, C. Lubich, and G. Wanner. Geometric Numerical Integration. Springer-

Verlag Berlin Heidelberg, 2006.

[45] T. Hanemann. Polymer-nanoparticle composites : From synthesis to modern appli-

cations. Material, 3 :3468–3517, 2010.

[46] R. J. Hardy. Formulas for determining local properties in molecular dynamics

simulations : Shock waves. The Journal of Chemical Physics, 76 :622, 1982.

[47] S. E. Harton, S. K. Kumar, H. C. Yang, T. Koga, K. Hicks, E. Lee, J. Mijovic,

M. Liu, R. S. Vallery, and D. W. Gidley. Immobilized polymer layers on spherical
nanoparticles. Macromolecules, 43(7) :3415–3421, 2010.
Bibliographie 97

[48] A. P. Holt, P. J. Griffin, V. Bocharova, A. L. Agapov, A. E. Imel, M. D. Dadmun,

J. R. Sangoro, and A. P. Sokolov. Dynamics at the polymer/nanoparticle interface
in poly(2-vinylpyridine)/silica nanocomposites. Macromolecules, 47(5) :1837–1843,
2014.

[49] W. G. Hoover. Canonical dynamics - equilibrium phase-space distributions. Physical

Review A, 31(3) :1695–1697, 1985.

[50] D. Hossain, M. A. Tschopp, D. K. Ward, J. L. Bouvard, P. Wang, and M. F. Hors-

temeyer. Molecular dynamics simulations of deformation mechanisms of amorphous
polyethylene. Polymer, 51(25) :6071–6083, 2010.

[51] Z.P. Huang, J. Wang, H.L. Duan, and B.L. Karihaloo. Size-dependent effective
elastic constants of solids containing nano-inhomogeneities with interface stress.
Journal of the Mechanics and Physics of Solids, 53 :1574–1596, 2005.

[52] C. Huet. Application of variational concepts to size effects in elastic heterogeneous

bodies. Journal of the Mechanics and Physics of Solids, 38(6) :813–841, 1990.

[53] P. H. Hünenberger. Thermostat algorithms for molecular dynamics simulations.

Advances in Polymer Science, 173 :105–149, 2005.

[54] F. Hussain, M. Hojjati, M. Okamoto, and R.E. Gorga. Polymer-matrix nanocom-

posites, processing, manufacturing, and application : An overview. Journal of Com-
posite Materials, 40 :17, 2006.

[55] J. H. Irving and J. G. Kirkwood. The statistical mechanical theory of transport

processes. iv. the equations of hydrodynamics. The Journal of Chemical Physics,
18 :817, 1950.

[56] E.T. Jaynes. Information theory and statistical mechanics. Physical Review,
106(4)/108(2) :620–630/171–190, 1957.

[57] T. Kanit, S. Forest, I. Galliet, V. Mounoury, and D. Jeulin. Determination of the

size of the representative volume element for random composites : statistical and
numerical approach. International Journal of Solids and Structures, 40 :3647–3679,
2003.

[58] J. Lane, D. Matthew, A.E. Ismail, M. Chandross, C.D. Lorenz, and G.S. Grest.
Forces between functionalized silica nanoparticles in solution. Phys. Rev. E, 79,
2009.
Bibliographie 98

[59] H. Le Quang and Q.C. He. Size-dependent effective thermoelastic properties of

nanocomposites with spherically anisotropic phases. Journal of the Mechanics and
Physics of Solids, 55 :1889–1921, 2007.

[60] R.B. Lehoucq and A. Von Lilienfeld-Toal. Translation of walter noll’s “derivation
of the fundamental equations of continuum thermodynamics from statistical me-
chanics”. Journal of Elasticity, 100 :5–24, 2010.

[61] H. Lequang and Q.C. He. Variational principles and bounds for elastic inhomoge-
neous materials with coherent imperfect interfaces. Mechanics of Materials, 40 :
865–884, 2008.

[62] G. Leu. Nmr characterization of elastomer-carbon black interactions. Macromole-

cules, 37 :6883–6891, 2004.

[63] D. R. Lide. Handbook of Chemistry and Physics. CRC Press, 1993.

[64] Wing Kam Liu and Eduard G. Karpov. Nano Mechanics and Materials : Theory,
Multiscale Methods and Applications. Wiley, 2006.

[65] Y.Q. Liu, H.T. Cong, W. Wang, C.H. Sun, and H.M. Cheng. Aln nanoparticle-
reinforced nanocrystalline al matrix composites : Fabrication and mechanical pro-
perties. Materials Science and Engineering : A, 505 :151–156, 2009.

[66] V. Marcadon. Effets de taille et d’interphase sur le comportement mécanique de

nanocomposites particulaires. PhD thesis, Ecole Polytechnique, France, 2005.

[67] V. Marcadon, D. Brown, E. Herve, P. Mele, N. D. Alberola, and A. Zaoui. Confron-

tation between molecular dynamics and micromechanical approaches to investigate
particle size effects on the mechanical behaviour of polymer nanocomposites. Com-
putational Materials Science, 79 :495–505, 2013.

[68] S. Marceau. Architecture multiéchelle de nanocomposites. PhD thesis, Université

de Savoie, France, 2003.

[69] D. Marx and J. Hutter. Ab Initio Molecular Dynamics : Basic Theory and Advanced
Methods. Cambrigde University Press, 2009.

[70] S. L. Mayo, B. D. Olafson, and W. A. Goddard. Dreiding - a generic force-field for

molecular simulations. Journal of Physical Chemistry, 94(26) :8897–8909, 1990.
Bibliographie 99

[71] N. N. Medvedev, A. Geiger, and W. Brostow. Distinguishing liquids from amor-

phous solids : Percolation analys is on the voronoi network. J. Chem. Phys., 93 :
8337–8342, 1990.

[72] N. Metropolis. Equation of state calculations by fast computing machines. Journal

of Chemical Physics, 21(6) :1087–1092, 1953.

[73] R. E. Miller and V. B. Shenoy. Size-dependent elastic properties of nanosized

structural elements. Nanotechnology, 11(3) :139–147, 2000.

[74] M. Moakher. On the averaging of symmetric positive-definite tensors. Journal of

Elasticity, 82 :273–296, 2006.

[75] J. Moczo and B. Pukanszky. Polymer micro and nanocomposites : Structure, in-
teractions, properties. Journal of Industrial and Engineering Chemistry, 14(5) :
535–563, 2008.

[76] A. I. Murdoch and D. Bedeaux. Continuum equations of balance via weighted

averages of microscopic quantities. Proceedings of the Royal Society of London.
Series A : Mathematical and Physical Sciences, 445(1923) :157–179, 1994.

[77] T. V. M. Ndoro, E. Voyiatzis, A. Ghanbari, D. N. Theodorou, M. C. Böhm, and

F. Müller-Plathe. Interface of grafted and ungrafted silica nanoparticles with a
polystyrene matrix : Atomistic molecular dynamics simulations. Macromolecules,
44(7) :2316–2327, 2011.

[78] T. V. M. Ndoro, M. C. Böhm, and F. Müller-Plathe. Interface and interphase

dynamics of polystyrene chains near grafted and ungrafted silica nanoparticles.
Macromolecules, 45(1) :171–179, 2012.

[79] W. Noll. Die herleitung der grundgleichungen der therm omechanik der kontinua aus
der statistichen mechanik. Journal of Rational Mechanics and Analysis, 445(1923) :
627–646, 1955.

[80] A. N. Norris. The isotropic material closest to a given anisotropic material. Journal
of Mechanics of Materials and Structures, 1 :223–238, 2006.

[81] S. Nose. A unified formulation of the constant temperature molecular-dynamics

methods. Journal of Chemical Physics, 81(1) :511–519, 1984.
Bibliographie 100

[82] A. Nouy and C. Soize. Random fields representations for stochastic elliptic boun-
dary value problems and statistical inverse problems. European Journal of Applied
Mathematics, 25 :339–373, 2014.

[83] G. M. Odegard, T. C. Clancy, and T. S. Gates. Modeling of the mechanical pro-

perties of nanoparticle/polymer composites. Polymer, 46(2) :553–562, 2005.

[84] A. Okabe, B. Boots, K. Sugihara, and S.N. Chiu. Spatial tessellations : concepts
and applications of Voronoi diagrams. John Wiley & Sons, Inc.,, 2000.

[85] M. Ostoja-Starzewski. Microstructural Randomness and Scaling in Mechanics of

Materials. Chapman & Hall/CRC/Taylor & Francis, 2008.

[86] B. Paliwal and M. Cherkaoui. Estimation of anisotropic elastic properties of nano-

composites using atomistic-continuum interphase model. International Journal of
Solids and Structures, 49(18) :2424–2438, 2012.

[87] B. Paliwal, M. Cherkaoui, and O. Fassi-Fehri. Effective elastic properties of nano-

composites using a novel atomistic-continuum interphase model. Comptes Rendus
Mecanique, 340(4-5) :296–306, 2012.

[88] G. J. Papakonstantopoulos, M. Doxastakis, P. F. Nealey, J. L. Barrat, and J. J.

de Pablo. Calculation of local mechanical properties of filled polymers. Physical
Review E, 75(3), 2007.

[89] A. Papon. Low-field nmr investigations of nanocomposites : Polymer dynamics and

network effects. Macromolecules, 44 :913–922, 2011.

[90] M. Parrinello and A. Rahman. Polymorphic transitions in single crystals : A new

molecular dynamics method. Journal of Applied Physics, 52 :7182–7190, 1981.

[91] R. D. Peng, H. W. Zhou, H. W. Wang, and L. Mishnaevsky. Modeling of nano-

reinforced polymer composites : Microstructure effect on young’s modulus. Com-
putational Materials Science, 60 :19–31, 2012.

[92] Z. S. Petrović, Y. J. Cho, I. Javni, S. Magonov, N. Yerina, D. W. Schaefer, J. Ilavsky,

and A. Waddon. Effect of silica nanoparticles on morphology of segmented poly-
urethanes. Polymer, 45(12) :4285–4295, 2004.

[93] S. Pfaller, G. Possart, P. Steinmann, M. Rahimi, F. Müller-Plathe, and M.C. Böhm.

A comparison of staggered solution schemes for coupled particle-continuum systems
modeled with the arlequin method. Computational Mechanics, 49(5) :565–579, 2012.
Bibliographie 101

[94] S. Pfaller, M. Rahimi, G. Possart, P. Steinmann, F. Muller-Plathe, and M. C.

Bohm. An arlequin-based method to couple molecular dynamics and finite element
simulations of amorphous polymers and nanocomposites. Computer Methods in
Applied Mechanics and Engineering, 260 :109–129, 2013.

[95] S. Plimpton. Fast parallel algorithms for short-range molecular-dynamics. Journal

of Computational Physics, 117(1) :1–19, 1995.

[96] L.S. Schadler P.M. Ajayan and P.V. Braun. Nanocomposite science and technology.
Wiley, 2003.

[97] R. Potestio, C. Peter, and K. Kremer. Computer simulations of soft matter : Linking
the scales. Entropy, 16 :4199–4245, 2014.

[98] M. Rahimi, H.A. Karimi-Varzaneh, M.C. Böhm, F. Müller-Plathe, S. Pfaller,

G. Possart, and P. Steinmann. Nonperiodic stochastic boundary conditions for
molecular dynamics simulations of materials embedded into a continuum mecha-
nics domain. The Journal of Chemical Physics, 134 :154108, 2011.

[99] M. Rahimi, I. Iriarte-Carretero, A. Ghanbari, M. C. Bohm, and F. Muller-Plathe.

Mechanical behavior and interphase structure in a silica-polystyrene nanocomposite
under uniaxial deformation. Nanotechnology, 23(30), 2012.

[100] D. C. Rapaport. The Art of Molecular Dynamics Simulation. Cambrigde University

Press, 2004.

[101] E. Riccardi, M. C. Bohm, and F. Muller-Plathe. Molecular dynamics method to

locally resolve poisson’s ratio : Mechanical description of the solid-soft-matter in-
terphase. Physical Review E, 86(3), 2012.

[102] F. Rizzi, H. N. Najm, B. J. Debusschere, K. Sargsyan, M. Salloum, H. Adalsteinsson,

and O. M. Knio. Uncertainty quantification in md simulations. part i : Forward
propagation. Multiscale Modeling & Simulation, 10(4) :1428–1459, 2012.

[103] F. Rizzi, H. N. Najm, B. J. Debusschere, K. Sargsyan, M. Salloum, H. Adalsteinsson,

and O. M. Knio. Uncertainty quantification in md simulations. part ii : Bayesian
inference of force-field parameters. Multiscale Modeling & Simulation, 10(4) :1460–
1492, 2012.

[104] B.A. Rozenberg and R. Tenne. Polymer-assisted fabrication of nanoparticles and

nanocomposites. Progress in Polymer Science, 33 :40–112, 2008.
Bibliographie 102

[105] C.H. Rycroft. Multiscale modeling in granular flow. PhD thesis, Massachusetts
Institute of Technology, 2007.

[106] M. Salloum, K. Sargsyan, R. Jones, B. Debusschere, H. N. Najm, and H. Adal-

steinsson. A stochastic multiscale coupling scheme to account for sampling noise
in atomistic-to-continuum simulations. Multiscale Modeling & Simulation, 10(2) :
550–584, 2012.

[107] H. Scholze. Glass : nature, structure, and properties. Springer-Verlag, 1991.

[108] C.E. Shannon. A mathematical theory of communication. Bell System Tech. J.,
27 :379–423 and 623–659, 1948.

[109] W. Shinoda, M. Shiga, and M. Mikami. Rapid estimation of elastic constants by

molecular dynamics simulation under constant stress. Physical Review B, 69(13),
2004.

[110] M. Shinozuka. Simulations of multivariate and multidimensional random processes.

The Journal of the Acoustical Society of America, 39 :357–367, 1971.

[111] M. Shinozuka and C.M. Jan. Digital simulation of random processes and its appli-
cations. Journal of Sound and Vibration, 25 :111–128, 1972.

[112] C. Soize. A nonparametric model of random uncertainties for reduced matrix models
in structural dynamics. Probabilistic Engineering Mechanics, 15 :277–294, 2000.

[113] C. Soize. A nonparametric model of random uncertainties for reduced matrix models
in structural dynamics. Probabilistic Engineering Mechanics, 15(3) :277–294, 2000.

[114] C. Soize. Non-gaussian positive-definite matrix-valued random fields for elliptic

stochastic partial differential operators. Computer Methods in Applied Mechanics
and Engineering, 195 :26–64, 2006.

[115] C. Soize. Tensor-valued random fields for meso-scale stochastic model of anisotropic
elastic microstructure and probabilistic analysis of representative volume element
size. Probabilistic Engineering Mechanics, 23 :307–323, 2008.

[116] C. Soize. A computational inverse method for identification of non-gaussian random

fields using the bayesian approach in very high dimension. Computer Methods in
Applied Mechanics and Engineering, 200 :3083–3099, 2011.
Bibliographie 103

[117] C. Soize and F. Poirion. Numerical methods and mathematical aspects for simula-
tion of homogeneous and non homogeneous Gaussian vector fields. In P. Krée and
W. Wedig, editors, Probabilistic Methods in Applied Physics, pages 17–53. Springer-
Verlag, Berlin, 1995.

[118] B. Staber and J. Guilleminot. Approximate solutions of lagrange multipliers for

information-theoretic random field models. SIAM/ASA Journal on Uncertainty
Quantification, accepté pour publication le 26 Mai 2015.

[119] F. W. Starr, T. B. Schroder, and S. C. Glotzer. Molecular dynamics simulation

of a polymer melt with a nanoscopic particle. Macromolecules, 35(11) :4481–4492,
2002.

[120] M. Sternitzke. Structural ceramic nanocomposites. Journal of the European Cera-

mic Society, 17 :1061–1082, 1997.

[121] W. C. Swope, H. C. Andersen, P. H. Berens, and K. R. Wilson. A computer-

simulation method for the calculation of equilibrium-constants for the formation
of physical clusters of molecules - application to small water clusters. Journal of
Chemical Physics, 76(1) :637–649, 1982.

[122] E.B. Tadmor and R.E. Miller. Modeling Materials : Continuum, Atomistic and
Multiscale Techniques. Cambridge University Press, 2011.

[123] S. Tsuneyuki. Molecular dynamics simulation of silica with a first-principles inter-

atomic potential. Molecular Engineering, 6 :157–182, 1996.

[124] M. E. Tuckerman. Statistical Mechanics : Theory and Molecular Simulation. Oxford

University Press, 2010.

[125] M. Vacatello. Monte carlo simulations of polymer melts filled with solid nanopar-
ticles. Macromolecules, 34(6) :1946–1952, 2001.

[126] M. Vacatello. Monte carlo simulations of the interface between polymer melts and
solids. ? chain stiffness. Macromolecular Theory and Simulations, 10(3) :187–195,
2001.

[127] M. Vacatello. Molecular arrangements in polymer-based nanocomposites. Macro-

molecular Theory and Simulations, 11(7) :757–765, 2002.

[128] M. Vacatello. Predicting the molecular arrangements in polymer-based nanocom-

posites. Macromolecular Theory and Simulations, 12(1) :86–91, 2003.
Bibliographie 104

[129] P. K. Valavala and G. M. Odegard. Modeling techniques for determination of

mechanical properties of polymer nanocomposites. Reviews on Advanced Materials
Science, 9(1) :34–44, 2005.

[130] P. K. Valavala, T. C. Clancy, G. M. Odegard, T. S. Gates, and E. C. Aifantis.

Multiscale modeling of polymer materials using a statistics-based micromechanics
approach. Acta Materialia, 57(2) :525–532, 2009.

[131] P. K. Valavala, G. M. Odegard, and E. C. Aifantis. Influence of representative

volume element size on predicted elastic properties of polymer materials. Modelling
and Simulation in Materials Science and Engineering, 17(4), 2009.

[132] B. W. H. van Beest, G. J. Kramer, and R. A. van Santen. Force fields for silicas
and aluminophosphates based on ab initio calculations. Physical Review Letters,
64 :1955, 1990.

[133] A. Vassiliou, D. Bikiaris, K. Chrissafis, K.M. Paraskevopoulos, S.Y. Stavrev, and

A. Docoslis. Nanocomposites of isotactic polypropylene with carbon nanoparticles
exhibiting enhanced stiffness, thermal stability and gas barrier properties. Compo-
sites Science and Technology, 68 :933–943, 2008.

[134] L. Verlet. Computer “experiments” on classical fluids. i. thermodynamical proper-

ties of lennard-jones molecules. Physical Review, 159 :98–103, 1967.

[135] G. G. Vogiatzis, E. Voyiatzis, and D. N. Theodorou. Monte carlo simulations of a

coarse grained model for an athermal all-polystyrene nanocomposite system. Eu-
ropean Polymer Journal, 47(4) :699–712, 2011.

[136] K. Vollmayr, W. Kob, and K. Binder. Cooling-rate effects in amorphous silica : A

computer-simulation study. Physical Review B, 54 :15808, 1996.

[137] N. Vu-Bac, T. Lahmer, H. Keitel, J. Zhao, X. Zhuang, and T. Rabczuk. Stochas-

tic predictions of bulk properties of amorphous polyethylene based on molecular
dynamics simulations. Mechanics of Materials, 68 :70–84, 2014.

[138] N. Vu-Bac, T. Lahmer, Y. Zhang, X. Zhuang, and T. Rabczuk. Stochastic predic-

tions of interfacial characteristic of polymeric nanocomposites (pncs). Composites
Part B-Engineering, 59 :80–95, 2014.

[139] N. Vu-Bac, R. Rafiee, X. Zhuang, T. Lahmer, and T. Rabczuk. Uncertainty quan-

tification for multiscale modeling of polymer nanocomposites with correlated para-
meters. Composites Part B-Engineering, 68 :446–464, 2015.
Bibliographie 105

[140] L. J. Walpole. 4th-rank tensors of the 32 crystal classes - multiplication tables.

Proceedings of the Royal Society of London Series a-Mathematical Physical and
Engineering Sciences, 391(1800) :149–179, 1984.

[141] J. Wang, H.L. Duan, Z.P. Huang, and Z.Y. Luo. Stress concentration tensors of
inhomogeneities with interface effects. Mechanics of Materials, 37 :723–736, 2005.

[142] J. Wang, H.L. Duan, Z. Zhang, and Z.P. Huang. An anti-interpenetration model
and connections bandween interphase and interface models in particle-reinforced
composites. International Journal of Mechanical Sciences, 47 :701–708, 2005.

Vous aimerez peut-être aussi

Machine thermique nano-électro-mécanique
Pas encore d'évaluation
Machine thermique nano-électro-mécanique
233 pages
Transport de Phonons en Nanostructures
Pas encore d'évaluation
Transport de Phonons en Nanostructures
171 pages
Comportement dynamique des structures non-linéaires
Pas encore d'évaluation
Comportement dynamique des structures non-linéaires
180 pages
Bacle P. 2017 Organisation Et Dynamique D'espèces Chargées (Électrolyte Argile)
Pas encore d'évaluation
Bacle P. 2017 Organisation Et Dynamique D'espèces Chargées (Électrolyte Argile)
251 pages
These Anne Lenoble
Pas encore d'évaluation
These Anne Lenoble
173 pages
Simulation des mémoires PCM par champs de phase
Pas encore d'évaluation
Simulation des mémoires PCM par champs de phase
121 pages
Modélisation des écoulements diphasiques
Pas encore d'évaluation
Modélisation des écoulements diphasiques
142 pages
Table These Collette 2007
Pas encore d'évaluation
Table These Collette 2007
4 pages
Dura 0299000505Z TH
Pas encore d'évaluation
Dura 0299000505Z TH
172 pages
Thèse sur la plasmonique à Strasbourg
Pas encore d'évaluation
Thèse sur la plasmonique à Strasbourg
220 pages
Propriétés Optiques des Nanostructures Métalliques
Pas encore d'évaluation
Propriétés Optiques des Nanostructures Métalliques
194 pages
HDR HLaurent
Pas encore d'évaluation
HDR HLaurent
147 pages
Fatigue des pièces Ti-6Al-4V : rugosité SLM/EBM
Pas encore d'évaluation
Fatigue des pièces Ti-6Al-4V : rugosité SLM/EBM
219 pages
Propriétés mécaniques des aérogels de silice
Pas encore d'évaluation
Propriétés mécaniques des aérogels de silice
141 pages
Étude des résines époxy pour rotomoulage
Pas encore d'évaluation
Étude des résines époxy pour rotomoulage
151 pages
Méthodes numériques en dynamique moléculaire
Pas encore d'évaluation
Méthodes numériques en dynamique moléculaire
26 pages
Mécanique Des Milieux Continus: To Cite This Version
Pas encore d'évaluation
Mécanique Des Milieux Continus: To Cite This Version
173 pages
Propriétés des Poudres Nanostructurées de Fer
Pas encore d'évaluation
Propriétés des Poudres Nanostructurées de Fer
144 pages
Proprietés Physique Meca Des Materiaux
100% (1)
Proprietés Physique Meca Des Materiaux
206 pages
Modélisation du transport thermique en thermoélectricité
Pas encore d'évaluation
Modélisation du transport thermique en thermoélectricité
189 pages
Modélisation Hyperélastique Innovante
Pas encore d'évaluation
Modélisation Hyperélastique Innovante
168 pages
Mécanique des milieux continus expliquée
Pas encore d'évaluation
Mécanique des milieux continus expliquée
173 pages
Homogénéisation des matériaux viscoélastiques
Pas encore d'évaluation
Homogénéisation des matériaux viscoélastiques
149 pages
Mécanique des Milieux Continus
100% (1)
Mécanique des Milieux Continus
484 pages
Trinite
Pas encore d'évaluation
Trinite
217 pages
Simulation de mise en forme à chaud
Pas encore d'évaluation
Simulation de mise en forme à chaud
164 pages
Simulation du Modèle d'Ising Dilué en Physique
Pas encore d'évaluation
Simulation du Modèle d'Ising Dilué en Physique
145 pages
Soufflage de membranes non-linéaires
Pas encore d'évaluation
Soufflage de membranes non-linéaires
161 pages
Modélisation du changement de phase interfacial
Pas encore d'évaluation
Modélisation du changement de phase interfacial
175 pages
Mécanique des milieux continus 2009-2010
Pas encore d'évaluation
Mécanique des milieux continus 2009-2010
484 pages
These - 2012 07 02
Pas encore d'évaluation
These - 2012 07 02
140 pages
Résistance des Matériaux: Cours ENA Oujda
Pas encore d'évaluation
Résistance des Matériaux: Cours ENA Oujda
64 pages
Modélisation des composants optoélectroniques
Pas encore d'évaluation
Modélisation des composants optoélectroniques
155 pages
Mécaniquedesmilieuxcontinus: François Sidoroff
Pas encore d'évaluation
Mécaniquedesmilieuxcontinus: François Sidoroff
172 pages
Mécanique des Milieux Continus
Pas encore d'évaluation
Mécanique des Milieux Continus
507 pages
Synthèse des phases MAX et MXènes
Pas encore d'évaluation
Synthèse des phases MAX et MXènes
153 pages
Modèle multi-bulles en cavitation
Pas encore d'évaluation
Modèle multi-bulles en cavitation
106 pages
Fatigue des Matériaux à Mémoire de Forme
Pas encore d'évaluation
Fatigue des Matériaux à Mémoire de Forme
42 pages
Mécanique des milieux continus expliquée
Pas encore d'évaluation
Mécanique des milieux continus expliquée
172 pages
Mécaniquedesmilieuxcontinus: François Sidoroff
Pas encore d'évaluation
Mécaniquedesmilieuxcontinus: François Sidoroff
172 pages
Projet d'étude sur la dalle d'essai CTE
100% (1)
Projet d'étude sur la dalle d'essai CTE
83 pages
Mécanique des milieux continus : Concepts clés
67% (3)
Mécanique des milieux continus : Concepts clés
172 pages
Gestion des bandes de localisation en composites
Pas encore d'évaluation
Gestion des bandes de localisation en composites
157 pages
Propriétés dynamiques des molécules adsorbées
Pas encore d'évaluation
Propriétés dynamiques des molécules adsorbées
162 pages
Dynamique des structures et méthodes d'analyse
Pas encore d'évaluation
Dynamique des structures et méthodes d'analyse
81 pages
Propagation Acoustique Dans Les Milieux Poreux Hétérogènes: To Cite This Version
Pas encore d'évaluation
Propagation Acoustique Dans Les Milieux Poreux Hétérogènes: To Cite This Version
134 pages
Caractérisation par ultrasons des FGM
Pas encore d'évaluation
Caractérisation par ultrasons des FGM
175 pages
Dissertation 1
Pas encore d'évaluation
Dissertation 1
5 pages
Chapitre
100% (1)
Chapitre
14 pages
Moniteur de Courant AC Par LM358
Pas encore d'évaluation
Moniteur de Courant AC Par LM358
14 pages
DTL 3eme
Pas encore d'évaluation
DTL 3eme
1 page
Contrôle de Réception des Marchandises
Pas encore d'évaluation
Contrôle de Réception des Marchandises
8 pages
Se mentir à soi-même : vérité ou illusion ?
0% (1)
Se mentir à soi-même : vérité ou illusion ?
3 pages
Convention Collective Enseignement Privé Hors Contrat
Pas encore d'évaluation
Convention Collective Enseignement Privé Hors Contrat
116 pages
La Fiscalité Du Cabinet Médical Au Maroc Fes Nov 2023 VF-2
Pas encore d'évaluation
La Fiscalité Du Cabinet Médical Au Maroc Fes Nov 2023 VF-2
52 pages
Étude des Frontières et Divisions Politiques
Pas encore d'évaluation
Étude des Frontières et Divisions Politiques
4 pages
Entretien Avec Un Chrétien À Qui Appartient L'argent Hugues Et Rodrigue
Pas encore d'évaluation
Entretien Avec Un Chrétien À Qui Appartient L'argent Hugues Et Rodrigue
37 pages
Resume Mondialisation
Pas encore d'évaluation
Resume Mondialisation
3 pages
Ô Jésus Christ, Transfiguré - Eric DASSIE
Pas encore d'évaluation
Ô Jésus Christ, Transfiguré - Eric DASSIE
2 pages
Rite, Idéologie, Fonction (P. Henriet)
Pas encore d'évaluation
Rite, Idéologie, Fonction (P. Henriet)
15 pages
Établissements avec option cinéma audiovisuel
Pas encore d'évaluation
Établissements avec option cinéma audiovisuel
18 pages
5 Chap1
Pas encore d'évaluation
5 Chap1
4 pages
Équipement Électrique: Caractéristiques
Pas encore d'évaluation
Équipement Électrique: Caractéristiques
26 pages
Combat Proterre : TTA 150 Édition 2018
Pas encore d'évaluation
Combat Proterre : TTA 150 Édition 2018
244 pages
Morphologie des Fleurs des Angiospermes
Pas encore d'évaluation
Morphologie des Fleurs des Angiospermes
43 pages
La Salsa Dura
100% (1)
La Salsa Dura
9 pages
Épreuve de mathématiques Baccalauréat 2023
Pas encore d'évaluation
Épreuve de mathématiques Baccalauréat 2023
4 pages
Analyse des chaînes alimentaires et décomposeurs
Pas encore d'évaluation
Analyse des chaînes alimentaires et décomposeurs
15 pages
TP Word 1 A 8 Vo2010 v1
Pas encore d'évaluation
TP Word 1 A 8 Vo2010 v1
34 pages
Escaliers PDF
Pas encore d'évaluation
Escaliers PDF
21 pages
Corrigé brevet des collèges 2020
Pas encore d'évaluation
Corrigé brevet des collèges 2020
3 pages
Routage Avec QoS Dans L'internet (QoS Routing in The Internet)
100% (1)
Routage Avec QoS Dans L'internet (QoS Routing in The Internet)
35 pages
Activités de compréhension en FLE
Pas encore d'évaluation
Activités de compréhension en FLE
1 page
CHAPITRE 2 (La Production À CT)
Pas encore d'évaluation
CHAPITRE 2 (La Production À CT)
4 pages
Anatomie Du Perinee
Pas encore d'évaluation
Anatomie Du Perinee
26 pages
Évaluation Diagnostique 5ème Sciences
Pas encore d'évaluation
Évaluation Diagnostique 5ème Sciences
2 pages
Blanc 2020 Assabila
Pas encore d'évaluation
Blanc 2020 Assabila
5 pages