These A LECOUTRE Gautier 2018
These A LECOUTRE Gautier 2018
These A LECOUTRE Gautier 2018
Par
M. Gautier Lecoutre
Composition du Jury :
Monsieur, Ganghoffer, Jean-François Professeur à l’Université de Lorraine Président
Monsieur, Lambin, Philippe Professeur de l’Université de Namur Rapporteur
Monsieur, Yvonnet, Julien Professeur à l’université Paris-est Rapporteur
Monsieur, Devel, Michel Professeur à l’ENSMM Directeur de thèse
Monsieur, Hirsinger, Laurent Chargé de recherche CNRS Codirecteur de thèse
Monsieur, Daher, Naoum Chargé de recherche CNRS Codirecteur de thèse
Remerciement
ii
iii
Résumé
l’effet d’un gradient de déformation. Bien qu’il existe dans tous les matériaux, ce phénomène
est encore rarement utilisé car il est en général de très faible amplitude. Cependant, à l’échelle
mettant en évidence des effets flexoélectriques importants. Pour cela, nous avons choisi de
conducteur. Dans le cadre des milieux continus, nous avons tout d’abord appliqué le principe
des puissances virtuelles et la thermodynamique des milieux continus pour obtenir de façon
systématique les équations constitutives d’un matériau aux comportements couplées semi-
flexoélectrique inverse de nano-objets tels que des nanotubes de carbone décrits atome par
atome, avec des dipôles électriques permanents et induits sur chaque atome. Moyennant
quelques hypothèses d’homogénéisation, nous avons couplé ces deux approches et obtenu les
équations reliant les quantités atomistiques, calculées dans la simulation, aux quantités
continus préalablement déterminées. Les premiers résultats mettent en évidence une variation
iv
The flexoelectricity tensor of a material characterizes its ability to polarize under the
action of a deformation gradient. The phenomenon is still rarely used though it exists in every
material, because the effects are usually very weak. However, for nanoscale systems,
flexoelectricity can be largely enhanced because of a possibly much greater gradient. Thus,
the aim of this PHD thesis is to build a model that would allow us to compute the
flexoelectric effects could be used for energy conversion. For that purpose, we have studied
point of view, we have applied the principle of virtual powers and classical thermodynamics
permanent and induced dipoles to simulate the inverse flexoelectric effect on the SWCNTs.
Using homogenization hypothesis, we have coupled these two approaches by obtaining the
equations binding the atomistic quantities computed in the numerical simulations, with the
The first numerical results seem to show a notable variation of the elements of the
v
Tables des matières
Introduction ....................................................................................................................................................... 1
Chapitre 1............................................................................................................................................................ 5
1 Différences entre piézoélectricité et flexoélectricité ........................................................... 6
2 Premiers travaux sur la flexoélectricité .................................................................................... 9
3 Calcul microscopique des contributions flexoélectriques............................................... 10
4 Les différentes contributions flexoélectriques .................................................................... 12
4.1 La flexoélectricité statique ................................................................................................................ 12
4.2 La réponse dynamique ....................................................................................................................... 14
4.3 Contribution surfacique de la piézoélectricité de surface à la flexoélectricité ........... 15
5 Résultats de Calculs ab-initio ..................................................................................................... 18
6 Effets de taille ................................................................................................................................... 20
7 Méthodes expérimentales de caractérisation de la flexoélectricité ............................ 22
8 Applications de la flexoélectricité ............................................................................................ 25
9 Conclusion ......................................................................................................................................... 26
Références ....................................................................................................................................................... 28
Chapitre 2......................................................................................................................................................... 33
1 Introduction...................................................................................................................................... 34
2 Notations globales .......................................................................................................................... 36
3 Objectivité ......................................................................................................................................... 38
3.1 Définitions ............................................................................................................................................... 38
3.2 Dérivées temporelles objectives..................................................................................................... 40
3.3 Différences entre l’approche énergétique et vectorielle ...................................................... 41
4 Les champs électromagnétiques ............................................................................................... 44
4.1 Les équations de Maxwell ................................................................................................................. 44
4.2 Prise en compte explicite des différents types de charges de conduction .................... 46
4.3 Force et couple pondéromoteur ..................................................................................................... 48
5 Equations d’équilibre électro-magnéto-thermo-mécanique .......................................... 51
5.1 Construction des puissances virtuelles ....................................................................................... 51
Champs de vitesses généralisées ................................................................................................................. 51
Champs de vitesses objectives ...................................................................................................................... 52
Puissances virtuelles des forces inertielles............................................................................................. 54
vi
Puissances virtuelles des forces internes ................................................................................................ 55
Puissances virtuelles des forces volumiques extérieures ................................................................ 56
Puissance virtuelle des forces de contacts .............................................................................................. 57
Les principes ......................................................................................................................................................... 58
5.2 Équations d’équilibre locales ........................................................................................................... 61
5.3 Equations de la thermodynamique ............................................................................................... 64
6 Transformation du second principe de la thermodynamique ....................................... 65
7 Equation constitutive d’un milieu déformable .................................................................... 67
7.1 Inégalité de Clausius-Duhem par rapport à la configuration matérielle ....................... 67
7.2 Les équations constitutives dans un milieu déformable ...................................................... 69
Comparaison avec des études précédentes ............................................................................................ 72
7.3 Quelques cas spéciaux de lois de comportement dissipatives........................................... 74
8 Conclusion ......................................................................................................................................... 76
References ....................................................................................................................................................... 77
Chapitre 3......................................................................................................................................................... 79
1 Introduction...................................................................................................................................... 80
2 Les nanotubes de carbone ........................................................................................................... 80
3 Principe des simulations de mécanique moléculaire classique .................................... 82
3.1 Statique moléculaire............................................................................................................................ 83
3.2 Dynamique Moléculaire ..................................................................................................................... 85
4 Modélisation de l’énergie potentielle d’interaction entre atomes de carbone dans
un système .................................................................................................................................................. 85
4.1 Propriétés d’une « bonne » forme d’énergie potentielle ...................................................... 85
4.2 Petit historique du développement de potentiels réactifs pour le carbone ................. 86
4.3 Les forces interatomiques ................................................................................................................. 90
4.4 Prise en compte des conditions aux limites mécaniques ..................................................... 92
5 Effet d’un champ électrique extérieur sur une population d’entités polarisables . 93
5.1 Moments dipolaires électriques dans un nanotube de carbone ....................................... 93
Dipôles induits, champ local et polarisabilité ........................................................................................ 93
Dipôles permanents dans une structure carbonée ............................................................................. 94
5.2 Définition et expressions des propagateurs de champ ......................................................... 96
Introduction : potentiel et champ électrostatique créés par des dipôles ponctuels ............ 96
Régularisation des propagateurs ................................................................................................................ 99
5.3 Contribution dipolaire électrique à l’énergie du système ................................................ 101
Expression générale ....................................................................................................................................... 101
Conditions d’équilibre électrique du système .................................................................................... 102
vii
Définition du système linéaire matriciel à résoudre ....................................................................... 103
Expressions simplifiée de l’énergie induite à l’équilibre électrique ......................................... 104
5.4 Expression analytique de la force ............................................................................................... 105
6 Déroulement de la programmation ....................................................................................... 107
7 Conclusion .......................................................................................................................................109
ix
Liste des figures
fig 1. Illustration de l’effet piézoélectrique direct : une contrainte mécanique appliquée est
à l’origine d’une polarisation électrique qui produit une différence de potentiel électrique
entre les faces du cristal [1]. ....................................................................................................... 6
fig 2. Illustration de l’effet piézoélectrique indirect : une tension électrique appliquée est à
l’origine d’une déformation du matériau [1]. ............................................................................. 6
fig 3. Schéma illustrant l’effet d’une déformation uniforme (à gauche) et non uniforme (à
droite) [4].................................................................................................................................... 7
fig 4. La déformation du nuage atomique sépare le barycentre des charges négatives avec
le noyau. La déformation crée donc un champ électrique. ......................................................... 8
fig 5. Courbure de molécules polaire de cristaux liquides [5] ............................................. 8
fig 6. Schéma du déplacement des atomes dans 2 cellules dus à une déformation
macroscopique. Cas d’une déformation uniforme dans un matériau centrosymétrique (a,b,c) et
non centrosymétrique (d,e,f) et avec un gradient de déformation homogène (g,h,i). Les figure
(a,d,g) montrent leur état initial, les figures (b,e,h) montre les déplacements dans le cas d’une
déformation externe. Les figures (c,f,i) montre le déplacement réél comprenant les
déformations externes et internes. [21] .................................................................................... 11
fig 7. Schéma représentant la modélisation de la contribution surfacique de la
piézoélectricité permettant une réponse flexoélectrique [25]................................................... 15
fig 8. (a) Les orbitales des atomes de carbone dans un plan de graphène sont
symétriques et perpendiculaires au plan des orbitales . (b) Dans le nanotube de carbone,
l’écart de chaque liaison d’un atome de carbone avec ses 3 plus proches voisins et le plan
formé par ces 3 plus proches voisins est caractérisé par un angle , qui va casser les
x
fig 12. Nanocomposite utilisant des fils flexoélectriques en vue d’obtenir un comportement
analogue à un matériau piézoélectrique [73]............................................................................ 26
fig 13. Schéma du domaine avec comme contour G et comme vecteur normal n à la
surface (qui peut être discontinue) ................................................................................... 36
fig 14. Schéma des interactions entre déformations, semi conduction et thermique [19] ... 46
fig 15. Nanotube de carbone (10,0) ..................................................................................... 80
fig 16. a) représentation « boules-bâtons » des 3 types de nanotubes de carbone (n,m). b)
schéma illustrant l’origine des entiers n et m comme coordonnées du vecteur circonférence.
[4] 81
fig 17. Dimensions caractéristiques dans le modèle utilisé pour le graphène ..................... 82
fig 18. Algorithme de minimisation de l’énergie potentielle des atomes ............................ 84
fig 19. Définition de l’angle dièdre utilisé lorsque BC est une liaison double .................... 89
fig 20. Illustration des termes de torsion et de liaisons faibles rajoutés dans le potentiel
AIREBO. 89
fig 21. Définition de l’angle pyramidal entre l’atome 𝛼 et les 3 atomes formant le plan 𝜋. 94
fig 22. Algorithme de recherche de l’équilibre statique du système soumis à l’action d’un
champ électrique extérieur 𝐸0................................................................................................ 107
fig 23. Déformation du domaine par la transformation de la configuration initiale 0 à la
fig 29. Variation de la composante G211 (en Å 1 ) par rapport au champ électrique extérieur
(en V / Å ) d’un nanotube de carbone (10,0) d’une longueur de 12,39 nm ........................... 125
fig 30. Variation de la composante G211 par rapport à la composante G112 G121 lors de la
flexion d’un nanotube de carbone (10,0) d’une longueur de 12,39 nm ................................ 126
fig 31. Relation entre la déformation imposée et la contrainte 11 ................................... 130
xi
fig 32. Représentation des dipôles permanents sur un nanotube (10,0) soumis à un champ
𝐸0, 𝑥 = 𝐸0, 𝑦 = 0,1 V·Å-1. Les positions à l’équilibre des atomes sont obtenues par
optimisation, sans prise en compte des dipôles permanents dans le calcul des dipôles induits à
l’équilibre. .............................................................................................................................. 147
fig 33. Dipôle induits à l’équilibre dans la même configuration que pour la fig 32. ......... 148
fig 34. Dipôles permanents avec prise en compte des dipôles permanents pour
l’optimisation ......................................................................................................................... 149
fig 35. Dipôles induits avec prise en compte des dipôles permanents pour l’optimisation.
149
fig 36. Distribution des projections des dipôles sur l’axe horizontal, selon la position le
long de l’axe du nanotube ...................................................................................................... 150
fig 37. Distribution des projections des dipôles sur la verticale, selon la position le long de
l’axe du nanotube ................................................................................................................... 151
fig 38. Variation de l’angle de flexion d’un nanotube (8,0) de longueur 6.54 nm par rapport
au carré du champ appliqué avec un angle de 45° par rapport à l’axe du nanotube lorsqu’il
n’est pas déformé. Pour cette courbe, nous n’avons pas pris en compte les dipôles permanents
152
fig 39. Variation du gradient du second ordre de la transformation d’un nanotube (8,0) de
longueur 6.54 nm par rapport au carré de la polarisation du matériau suivant x ................... 154
fig 40. Variation de l’angle de flexion d’un nanotube (8,0) de longueur 6.54 nm par rapport
au carré du champ appliqué avec un angle de 45° par rapport à l’axe du nanotube lorsqu’il
n’est pas déformé. Pour cette courbe, nous avons pris en compte les dipôles permanents .... 155
fig 41. Variation du gradient de transformation 𝐺121 d’un nanotube (8,0) de longueur 6.54
nm par rapport au carré de la polarisation du matériau suivant x. Pour cette courbe, nous
avons pris en compte les dipôles permanents ......................................................................... 156
fig 42. Variation du coefficient f2211 d’un nanotube de carbone selon son rayon pour une
longueur constante de 20,6 nm ............................................................................................... 158
fig 43. Variation du coefficient f 2121 d’un nanotube de carbone selon son rayon pour une
longueur constante de 20.6 nm ............................................................................................... 158
fig 44. Coefficient flexoélectrique 𝑓2121 selon la longueur du nanotube (10,0) ............. 159
fig 45. Coefficient flexoélectrique 𝑓2211 selon la longueur du nanotube (10,0) ............. 160
xii
Définition des vecteurs t et t- G
+
fig 46. tangent à l’arête où le vecteur normal unitaire
extérieur n de la surface fermée est discontinue.......................................................... 165
xiii
Introduction
L’effet piézoélectrique n’existe que pour des matériaux non centrosymétriques. En
déformation homogène. Mais une déformation non homogène peut briser les symétries d’un
sous l’effet d’une déformation qui n’est pas uniforme, cet effet est appelé « la
flexoélectricité ». Ce phénomène est encore rarement utilisé bien qu’il existe dans tous les
matériaux, car il est en général d’amplitude très faible. Mais, pour des dimensions
suffisamment petites, une différence de déformations sur une faible distance donne un grand
Toutefois, au début des années 2000, Cross et ses co-auteurs ont mesuré des effets de la
flexoélectricité beaucoup plus importants que prévu pour des matériaux solides
les uns par rapport aux autres, à cause du gradient de déformation, qui provoquent une
polarisation supplémentaire du fait du mouvement relatif des barycentres des ions positifs et
négatifs dans chaque maille. Cependant, les gradients de déformation n’affectent pas
seulement la position des ions, mais affecte aussi la distribution de la densité des électrons qui
va donc affecter la polarisation totale. Les effets flexoélectriques dans les structures carbonées
grande résistance à la déformation de ces structures par rapport aux céramiques. Tout ceci est
1
détaillé dans la description sommaire de l’état de l’art du domaine qui constitue le premier
chapitre de ce manuscrit de thèse, afin de justifier notre intérêt pour la flexoélectricité dans les
structures carbonées.
des milieux continus, nous avons obtenu des contraintes sur les lois de comportement qui
peuvent être utilisées pour décrire des matériaux semi-conducteurs piézoélectriques et/ou
flexoélectriques. Pour cela, nous avons utilisé différents principes toujours valables pour la
Clausius-Duhem. C’est en utilisant cette dernière inégalité que nous avons pu retrouver les
approche de type physique des milieux continus, les premiers chercheurs à s’intéresser à ce
phénomène ont rajouté à la théorie linéaire de la piézoélectricité des termes faisant intervenir
plus possible notre travail et pour d’éventuelles recherches dans le futur, nous avons
Dans une troisième partie, nous avons repris un modèle statique d’interactions entre
dipôles électriques distribués sur chaque atome pour développer un modèle numérique
pouvant simuler la déformation d’un matériau non-métallique sous l’effet d’un champ
électrique extérieur. Nous verrons que la déformation est due à l’effet mécanique des couples
2
visant à aligner sur le champ extérieur les dipôles électriques induits par ce champ extérieur et
par les dipôles permanents créés par la courbure du feuillet de carbone pour créer le nanotube,
ainsi que par la déformation de la structure à cause du champ extérieur. Nous faisons alors
l’hypothèse que les positions d’équilibre des atomes peuvent être calculées par minimisation
itérative de l’énergie potentielle totale des liaisons entre atomes de carbone, modélisée par
potentielle décrivant les effets du champ extérieur grâce aux interactions avec les dipôles
distribués. Nous montrons alors comment le calcul de ces dipôles atomiques peut être effectué
Avant d’en arriver là, nous mettrons en relation la théorie des milieux continus et la
caractéristiques de la flexoélectricité dans l’approche milieu continu, à partir des résultats des
simulations atomistiques. Pour cela, nous avons été conduit à faire des hypothèses
d’homogénéisation pour faire le lien entre les deux approches. Pour la partie mécanique, cette
hypothèse consiste à considérer que le gradient de déformation est homogène sur un élément
déformation ou d’élasticité. Pour la partie électrique, nous avons aussi dû faire une hypothèse
faisant l’hypothèse que les dipôles supplémentaires dus à la présence du champ électrique
extérieur sont identiques dans un élément du solide considéré. A partir de cette hypothèse,
3
hypothèses d’homogénéisation, nous avons pu obtenir les équations reliant les quantités
Enfin, nous avons pu étudier la variation des éléments de l’un des tenseurs de
étudier l’impact des dipôles permanents sur la flexion des nanotubes de carbone soumis à un
4
Chapitre 1
La propriété d’un solide cristallin comme le quartz de se déformer sous une tension
électrique a été découverte par les frères Curie et a été nommé « piézoélectricité » par Hankel.
Un matériau est piézoélectrique lorsqu’il se polarise électriquement sous l’effet d’une
contrainte appliquée suivant certaines directions (voir fig 1).
fig 1. Illustration de l’effet piézoélectrique direct : une contrainte mécanique appliquée est
à l’origine d’une polarisation électrique qui produit une différence de potentiel
électrique entre les faces du cristal [1].
fig 2. Illustration de l’effet piézoélectrique indirect : une tension électrique appliquée est à
l’origine d’une déformation du matériau [1].
6
ij cijkl
E
Skl eijk Ek (1)
défini pour un champ électrique macroscopique fixé. ijS donne la susceptibilité électrique
fig 3. Schéma illustrant l’effet d’une déformation uniforme (à gauche) et non uniforme (à
droite) [4]
En revanche, une polarisation électrique peut se créer même dans des matériaux
centrosymétriques si on induit une déformation non uniforme dans le matériau. Ce mécanisme
correspond à la flexoélectricité. Une déformation non uniforme appliquée (par exemple une
flexion) est à l’origine d’une brisure de symétrie dans ces cristaux centrosymétriques. Elle fait
apparaitre une polarisation électrique (voir l’image de droite de la fig 3). C’est-à-dire une
séparation des barycentres des charges positives et négatifs qui composent le cristal.
7
Bien que dans la plupart des cas on définit la flexoélectricité comme la création d’une
polarisation induite par une dissymétrie des ions dans une maille dû à un gradient de
déformation, la polarisation peut en réalité être induite par 3 mécanismes différents :
www.MaterialsViews.com
fig 4. La déformation du nuage atomique sépare le barycentre des charges négatives avec le
noyau. La déformation crée donc un champ électrique.
8
2 Premiers travaux sur la flexoélectricité
2U j
Pi ˆijkl eijk S jk 0 ij E j (3)
xk xl
9
Les premiers calculs atomistiques des coefficients de la flexoélectricité ont été menés
par Askar et al. [13] en 1970 pour NaCl, NaI, KI et KCl, à l’aide de la théorie développée par
Mindlin. Bursian et al. ont ensuite étudié l’effet de la courbure d’une poutre de BaTiO3 sur la
polarisation ainsi créée [14]. Ils découvrirent que l’effet est d’autant plus grand que la
permittivité électrique du matériau est grande. Ce phénomène était alors en général décrit
comme étant « une piézoélectricité non locale » [15]. En 1970, Axe et al. [16] découvrirent
une manifestation de la flexoélectricité en étudiant le spectre de phonons à basse énergie dans
un matériau pérovskite ferroélectrique. Selon Tagantsev, le mot flexoélectricité fait son
apparition pour les structures cristallines pour la première fois en 1981 dans un article de
Indenbom, inspiré par la théorie des cristaux liquides de De Gennes, pour décrire des effets
similaires.
R nj
Ui
w n,i x j
dy j w nint,i (4)
x 0j
déplacement de l’atome n situé en Rn,i selon l’axe i peut s’écrire de la façon suivante :
10
fig 6. Schéma du déplacement des atomes dans 2 cellules dus à une déformation
macroscopique. Cas d’une déformation uniforme dans un matériau centrosymétrique
(a,b,c) et non centrosymétrique (d,e,f) et avec un gradient de déformation homogène
(g,h,i). Les figures (a,d,g) montrent leur état initial, les figures (b,e,h) montre les
déplacements dans le cas d’une déformation externe. Les figures (c,f,i) montre le
déplacement réél comprenant les déformations externes et internes. [21]
Supposons maintenant que l’on puisse écrire le déplacement interne w nint,i sous forme
1 Ui Uk 1
ik :
2 xk xi
1
Historiquement, Tagantsev a d’abord étudié la contribution flexoélectrique par rapport au gradient de la
déformation non symétrisé [17], [18]. Il a par la suite choisit de prendre en compte la contribution
flexoélectrique par rapport gradient de la déformation symétrisé. Par soucis de synthétiser ses travaux, nous
avons choisi d’utiliser le gradient de la déformation symétrisé dans les calculs.
11
ik
w nint,j Nnikl, j (6)
xl
ik
Pj V 1QnNnikl, j Pjext (7)
xl
Le premier terme dans (7) traduit la flexoélectricité dans le volume V puisqu’il s’agit
d’une proportionnalité entre la polarisation et le gradient de polarisation, on identifie le
tenseur de flexoélectricité:
Le dernier terme dans (7) qui est dû à la présence du premier terme dans (4) correspond à
l’effet de surface flexoélectrique. Cet effet a été introduit pour la première fois par Resta [18].
Tagantsev et al. [16] [17] identifient ainsi plusieurs types de contributions dans les
réponses de type flexoélectrique, qui n’ont pas d’analogue avec la piézoélectricité. En 2006,
Maranganti et al. élaborèrent aussi une théorie afin d’expliquer phénoménologiquement et
macroscopiquement la flexoélectricité [22]. Ces derniers ont utilisé le principe variationnel
développé par Mindlin [11] et Toupin [12] pour trouver les équations d’équilibre et les
équations constitutives de la flexoélectricité avec une approche de type dynamique des
réseaux. Notons que le choix de prendre comme variable le gradient de déformation symétrisé
ik
(dérivée par rapport à la déformation [21] ou non symétrisé (seconde dérivée par rapport
xl
2U j
au déplacement ) [23] dépend des chercheurs.
xk xl
12
kl
polarisation Pi et le gradient de déformation j kl , en l’absence de champ électrique
x j
extérieur. On introduit alors le tenseur flexoélectrique dont les coefficients sont définis par:
P
ijkl i
j kl
(9)
E 0
1 c P
0 P ² ² P f1P f2 PE
1
G (P ( x ), ( x ), x ) a P ² (10)
2 2 x x
f1 f2 ( P )
G (11)
2 x
Avec
1 c P
P ² ² P f P PE
x
(12)
2 0 2 x
13
l’intégrale soit minimale, on doit alors résoudre l’équation d’Euler-Lagrange
G d G
0 , où A représente l’une des variables du problème ( et P ). On
A dx A
x
obtient alors [21]:
P 0 E e (13)
x
P
f c P (14)
x
avec f 0 et e 0 (15)
P 0 E 0 MU 0 P (17)
x
14
²P
U c MP (18)
x 0 x ²
Dans un système fini, l’effet de surface a en général une contribution beaucoup plus
faible que l’effet du volume mais peut ne pas être négligeable dans certains cas, comme par
exemple les couches minces. Dans cette partie, nous reprendrons le calcul de Tagantsev et
Yurkov dans le calcul de la contribution surfacique de la piézoélectricité [25]. Comme le
montrent Tagantsev et al. [25], dans un matériau centrosymétrique, la symétrie est cassée à la
surface du matériau, une fine couche d’épaisseur devient alors piézoélectrique avec un
tenseur de piézoélectricité hijk .
Tagantsev et al. modélisent alors une plaque dépaisseur h ayant deux couches minces à
ses extrémités, qui sont piézoélectriques et d’épaisseur . L’épaisseur est supposée
largement inférieure à l’épaisseur h. Supposons que la plaque est uniformément déformée par
une courbure G et isotrope.
15
Si on choisit une densité volumique de potentiel thermodynamique de forme quadratique
dans la couche dépendant des déformations et de la polarisation P , on a :
F
2
P2 P h333 33 h311 11 22
c11 2
2
11 22
2
33
2
(19)
c12 11 22 11 33 22 33
où c11 c1111 , c12 c1122 , donne l’inverse de la susceptibilité électrique à déformation
constante.
Puisque la plaque est fine, on peut supposer que nous avons une relation entre la
courbure G et les déformations :
hG
11 22 (20)
2
D’après les conditions aux limites, la surface est mécaniquement libre, on a alors les
F
contraintes nulles aux limites et donc 33 0 . On a alors :
33
c12 h
33 hG 333 P (21)
c11 c11
En utilisant les équations (20) et (21), on réécrit l’équation (19) en multipliant cette
dernière équation par , on obtient une nouvelle expression de la densité d’énergie surfacique
:
1 2 c
2 P h311 12 h333 hGP 1
2 0 c11
(22)
1
h2
1
Avec 0 333
c11
et donc, donne la vraie susceptibilité diélectrique dans les conditions mécaniques actuelles.
1 donne la contribution mécanique pure. Pour connaitre le champ électrique dans la couche
, on utilise l’équation d’état E , on obtient alors :
P
16
P 0 E ehG (23)
h f 11
D e
2 f h x3
(26)
c12
Avec 0 0 et e 0 (h311 h33 )
c11
h 11
Dans le cas d’une couche mince, on a D e . On trouve alors le coefficient
x3
flexoélectrique suivant :
h
1133
eff
e (27)
eff e
Tagantsev et al. ont pu évaluer le coefficient f eff
en évaluant l’épaisseur
0 f
0.4nm avec 10 et e 1 C.m2 . A partir de ces valeurs, il trouve une valeur de
0
f eff 4V , ce qui montre que l’effet surfacique peut être très important lorsqu’on le compare
avec les effets dus au volume ( f 1 10V ).
De même, Shen et Hu [26] ont élaboré une théorie comprenant l’effet de surface en
s’inspirant des travaux de Mindlin et de Gurtin et Murdoch [27].
17
Par cette contribution, si la polarisation varie continuellement du volume à la surface du
matériau alors le gradient de polarisation crée va induire une contrainte mécanique due à
l’effet inverse de la flexoélectricité [28], [29]. Cet effet peut être aussi fort que la
flexoélectricité dans le volume.
que la dépendance entre les gradients de déformation et les champs électriques dans
l’expression de l’énergie s’explique par une invariance de Jauge. A partir de calculs ab-initio,
il démontre les effets statiques et dynamiques de la flexoélectricité, ainsi que les effets de
surface.
ainsi que des pérovskites (BaTiO3, SrTiO3..) afin de calculer la contribution flexoélectrique
statique [35].
xx
Px xxxx (28)
x
Par la suite, Stengel, Hong et Vanderbilt [36], [37] développèrent une méthode pour
Par la suite, Hong et al. [35], [37], [38] trouvèrent que les coefficients étaient de l’ordre de
0.1 1 nC / m ou alors que les coefficients f avaient une valeur allant de 10 à 20 Volts. On
notera que les coefficients flexoélectrique sont définis à champ électrique constant mais
que dans certains cas, les calculs se font à déplacement électrique constant [37]. Ce calcul D
18
est utile pour un calcul avec une constante diélectrique constante alors que le calcul de est
De même, Ponomareva et al. [39] ont utilisé une simulation ab-initio, en développant
caractérisant la flexoélectricité ont été réalisés en 2002 par Dumitrica et al. [40], mais le lien
entre leurs résultats de calculs et la flexoélectricité n’a été remarqué que 13 ans plus tard par
Kvashnin et al. [41] qui attribuèrent la création de dipôles permanents selon l’angle des
liaisons entre atomes de carbone (voir fig 8), selon le rayon du nanotube et sa chiralité, à
l’effet flexoélectrique. Kvashnin et al. purent ainsi écrire la relation entre le dipôle et l’angle
pyramidal ou la courbure c :
p m M.c (29)
fig 8. (a) Les orbitales des atomes de carbone dans un plan de graphène sont
symétriques et perpendiculaires au plan des orbitales . (b) Dans le nanotube de
carbone, l’écart de chaque liaison d’un atome de carbone avec ses 3 plus proches
voisins et le plan formé par ces 3 plus proches voisins est caractérisé par un angle
, qui va casser les symétries de distribution de charges des orbites . [40]
Kalinin et Meunier [42] ont également calculé la polarisation induite par la courbure
de la structure carbonée :
M
p (30)
R
19
où R est le rayon de courbure de la structure. Ils trouvèrent un coefficient presque identique à
celui de Dumitrica et de Kvashnin et al. Dumitrica et al. ont trouvé un dipôle de 0.82 D.Å ,
Kalinin et Meunier ont trouvé un dipôle entre 0.75 et 0.9D.Å et Kvashnin et al. ont trouvé
0.8D.Å .
fig 9. (a) Schéma d’une feuille de graphène en flexion. (b) Graphique donnant le moment
dipolaire par atome en fonction du rayon de courbure
Notons au passage que la piézoélectricité peut exister dans le graphène sous réserve
6 Effets de taille
flexoélectricité sur la constante diélectrique dans des condensateurs, ce qui leur permit
d’expliquer les phénomènes comme le « dead layer effect » pour lesquels il n’était pas
20
trouvèrent que l’effet flexoélectrique devient de plus en plus important lorsque la taille
diminue comme attendu du fait du gradient et devient non négligeable par rapport à la
piézoélectricité pour une épaisseur inférieure à 2 nm. Pour le PZT, Majdoub et al. trouvèrent
néanmoins que l’énergie récupérée par l’effet piézoélectrique est doublée lorsque son
épaisseur est de l’ordre de 20 nm, grâce à l’effet flexoélectrique [45]. Majdoub et al. [44] et
Liu et al. [46], [47] trouvèrent qu’il y a un effet de taille très important de la piézoélectricité
flexoélectricité, Liu et al. ont ainsi montré qu’il est possible d’avoir un effet piézoélectrique
21
7 Méthodes expérimentales de caractérisation de la
flexoélectricité
fig 10. Méthodes les plus utilisées pour quantifier les réponses flexoélectriques par flexion (à
gauche) ou compression (à droite) d’un matériau [4]
(voir fig 10 a). Cette méthode a pu être utilisée sur des polymères [48] ou plus récemment par
Cross afin de quantifier la réponse flexoélectrique d’une céramique de type pérovskite [49]–
22
Comme le montrent Zubko et al. [55], la déformation 11 peut induire d’autres
déformations par le couplage entre déformations causé par un coefficient de Poisson non nul.
Par exemple, pour une poutre ayant des propriétés mécaniques isotropes, on a :
12 11 1 12 (32)
comme le PMN où ils ont trouvé 12 4 106 C / m à température ambiante et ont étudié la
dépendance du tenseur selon la température et la permittivité [49], [50]. Ils se sont aperçu que
les effets flexoélectriques sont beaucoup plus importants que les estimations de Marvan[56]
(10 000 fois supérieur). Ils ont étudié le BST [53] au-dessus de sa température de Curie de
ils obtiennent une valeur de 12 de 100 106 C / m pour le BST et de seulement de
1.4 106 C / m pour le PZT [54] . Cependant, la proportionnalité entre le tenseur
flexoélectrique et la permittivité n’est plus vraie lorsque la permittivité électrique devient trop
importante et que les phénomènes non linéaires entrent en jeu. Pour le BaTiO3 , Cross et al.
ont calculé un coefficient 12 10 106 C / m [51]. D’importants effets ont été calculés par
Baskaran et al. sur des couches minces de Polyvinylidene Fluoride (PVDF) avec des valeurs
pas suffire pour connaitre tous les coefficients flexoélectrique. Récemment, Zhang et al. ont
proposé de nouvelles méthodes pour calculer différents coefficients en appliquant une torsion
sur un matériau en PVDF afin de créer une déformation inhomogène. Ils trouveront plusieurs
coefficients flexoélectriques effectifs 2312 1.1 108 C.m1 [60], 1211 7.3 1010 C.m1
Une autre méthode consiste à compresser un matériau ayant une forme pyramidale
(voir fig 10b). Les contraintes sur les faces du dessus et du dessous n’étant pas les mêmes, on
23
1011 ou 1010 C / m à partir de cette méthode [56]. Plus récemment, Ma et Cross et
Baskaran et al. ont utilisé cette méthode sur le PMN (Plomb Magnesium Niobate), BST
déformation qui en résulte sur le BST [64] ou le PMN-PT [41] [65]. Toutefois, dans ces
expériences, l’électrostriction étant le phénomène dominant, il faut alors soustraire les effets
dus à l’électrostriction des effets dus à la flexoélectricité.
fréquence des ondes acoustiques transverses excitées par des gradients de polarisation. Cette
Zubko et .al. réalisèrent aussi des expériences sur le SrTiO3 afin de calculer les
orientations du réseau. Ils trouvèrent des valeurs de l’ordre de 109 C/m à température
ambiante mais peut-être 100 fois supérieures à basse température puisque la permittivité
flexoélectricité sur les pérovskites diélectriques et des semi-conducteurs pour pouvoir les
comparer aux résultats des expérimentateurs (Zubko, Catalan, Cross). A partir des méthodes
de calcul étudiant les ondes acoustiques dans les cristaux, développées par Tagantsev, ils vont
trouver des résultats très proches des valeurs de Askar sur le NaCl et le KCl. [13]
constante diélectrique dans les couches minces. Ils supposèrent que la déformation imposée
ne se relaxe pas directement sur la surface du film mais de manière exponentielle dans le
volume de la couche mince. Ils purent ainsi mettre en évidence les effets de surface dans le
24
calcul des contributions de la flexoélectricité. Lorsque la flexoélectricité induite par cette
déformation est introduite dans l’énergie libre du système, la constante diélectrique diminue.
En faisant varier la température, ils montrèrent que les coefficients flexoélectriques dans le
8 Applications de la flexoélectricité
En travaillant la forme du matériau (voir fig 11) Fousek et al. [70] ont pu réaliser un
nanocomposite non piézoélectrique ayant une réponse globale analogue à celle d’un matériau
piézoélectrique grâce aux propriétés flexoélectriques. Cross et al. ont fait des expériences sur
un matériau de forme pyramidale avec un fort coefficient flexoélectrique (BST) et compressé
par deux coques métalliques (voir fig 11). Lorsque la compression a lieu, un gradient de
déformation se crée sur chaque pyramide et induit une polarisation électrique [63], [71]. Par la
suite, Sharma et al ont étudié ce phénomène théoriquement [72]. Ces derniers ont eu l’idée de
relier plusieurs couches de surfaces présentant des déformations non homogènes avec des
matériaux non piézoélectriques afin de simuler un matériau piézoélectrique.
fig 11. Matériau piézoélectrique effectif utilisant les effets flexoélectrique en excitant
longitudinalement le matériau [71]
25
fig 12. Nanocomposite utilisant des fils flexoélectriques en vue d’obtenir un comportement
analogue à un matériau piézoélectrique [73]
A partir des effets flexoélectriques, Gruvermann et al. [74] ont montré qu’il était
possible de changer le sens de la polarisation sur des condensateurs PZT en créant une flexion
sur le substrat en Si. De même, Lu et al. [75] ont pu inverser le sens de polarisation d’une
couche mince en BaTiO3 en créant une déformation non homogène à l’aide d’une pointe
[62] ou sur du Titanate de Strontium (STO) [76], [77]. Ils obtiennent un meilleur apport en
énergie lorsque la poutre se rapproche de l’échelle du nanomètre. Enfin, Han et al. ont montré
que l’on pouvait obtenir des nanogénérateurs en faisant croitre directement des nanoparticules
de PZT sur des nanotubes de carbone multiparois. Ils purent ainsi obtenir dans le meilleur des
cas une tension crête de 8,6 V pour un courant de 4,7 nA, lorsqu’une force de 20 N était
9 Conclusion
Nous présentons maintenant les quelques résultats les plus important pour la suite. Nous
avons vu que les lois constitutives donné par Tagantsev et al. [4], [21] caractérisant la
flexoélectricité directe et indirecte sont donnés par :
26
P 0 E e (34)
x
P
f c P (35)
x
avec f 0 et e 0 (36)
flexocouplage f ou de flexoélectricité est très différente selon qu’on les calcule par
flexoélectricité peuvent être différents selon les conditions expérimentales. Comme nous
l’avons vu dans les travaux de Hong et Vanderbilt [37], les coefficients sont différents si on
27
Références
[1] P. Dineva, D. Gross, R. Müller, et T. Rangelov, « Piezoelectric Materials », in Dynamic
Fracture of Piezoelectric Materials, vol. 212, Cham: Springer International Publishing,
2014, p. 7‑32.
[2] G. Lippmann, « Principe de la conservation de l’électricité, ou second principe de la
théorie des phénomènes électriques », J. Phys. Théorique Appliquée, vol. 10, no 1, p. 381‑
394, 1881.
[3] P. Curie et J. Curie, « Contractions et dilatations produites par des tensions électriques
dans les cristaux hémièdres à faces inclinées », Comptes Rendus L’Académie Sci., vol. 93,
p. 1137, 1881.
[4] P. Zubko, G. Catalan, et A. K. Tagantsev, « Flexoelectric effect in solids », Annu. Rev.
Mater. Res., vol. 43, p. 387–421, 2013.
[5] T. D. Nguyen, S. Mao, Y.-W. Yeh, P. K. Purohit, et M. C. McAlpine, « Nanoscale
Flexoelectricity », Adv. Mater., vol. 25, no 7, p. 946‑974, févr. 2013.
[6] K. B. Tolpygo, « Long wavelength oscillations of diamond-type crystals including long
range forces », Sov. Phys. SOLID STATE, vol. 4, no 7, p. 1297‑1301, janv. 1963.
[7] V. S. Mashkevich et K. B. Tolpygo, « Electrical, Optical and Elastic Properties of
Diamond Type Crystals 1 », J Exper Theor. Phys, vol. 32, p. 520‑525, mars 1957.
[8] V. S. Mashkevich, « Electrical, Optical, and Elastic Properties of Diamond-Type Crystals
II. Lattice Vibrations with Calculation of Atomic Dipole Moments », Sov. Phys. JETP,
vol. 5, no 4, 1957.
[9] S. M. Kogan, « Piezoelectric effect during inhomogeneous deformation and acoustic
scattering of carriers in crystals », Sov. Phys. SOLID STATE, vol. 5, no 10, p. 2069‑2070,
avr. 1964.
[10] P. Harris, « Mechanism for the Shock Polarization of Dielectrics », J. Appl. Phys., vol.
36, no 3, p. 739, 1965.
[11] R. D. Mindlin, Polarization Gradient in Elastic Dielectrics. Vienna: Springer Vienna,
1972.
[12] R. Toupin, « The Elastic Dielectric », Journal of Rational Mechanics and Analysis, p.
849‑915, 1956.
[13] A. Askar, P. C. Y. Lee, et A. S. Cakwak, « Lattice-Dynamics Approach to the Theory
of Elastic Dielectrics with Polarization Gradient », Phys. Rev. B, vol. 1, no 8, oct. 1969.
[14] E. V. Bursian et O. I. Zaikovskii, « Changes in the curvature of a ferroelectric film due
to polarization », Sov. Phys. SOLID STATE, vol. 10, no 5, p. 1121‑1124, 1968.
[15] E. V. Bursian et N. N. Trunov, « Nonlocal piezoelectric effect », Sov. Phys. SOLID
STATE, vol. 10, p. 760‑762, 1974.
28
[16] J. D. Axe, J. Harada, et G. Shirana, « Anomalous Acoustic Dispersion in
Centrosymmetric Crystals with Soft Optic Phonons », Phys. Rev. B, vol. 1, no 3, p. 1227‑
1234, févr. 1970.
[17] A. K. Tagantsev, « Theory of flexoelectric effect in crystals », Zh Eksp Teor Fiz, vol.
88, p. 2108‑2122, 1985.
[18] A. K. Tagantsev, « Piezoelectricity and flexoelectricity in crystalline dielectrics »,
Phys. Rev. B, vol. 34, no 8, p. 5883, 1986.
[19] A. K. Tagantsev, « Pyroelectric, piezoelectric, flexoelectric, and thermal polarization
effects in ionic crystals », Sov. Phys. Uspekhi, vol. 30, no 7, p. 588‑603, juill. 1987.
[20] R. M. Martin, « Piezoelectricity », Phys. Rev. B, vol. 5, no 4, p. 1607, 1972.
[21] P. V. Yudin et A. K. Tagantsev, « Fundamentals of flexoelectricity in solids »,
Nanotechnology, vol. 24, no 43, p. 432001, nov. 2013.
[22] R. Maranganti, N. D. Sharma, et P. Sharma, « Electromechanical coupling in
nonpiezoelectric materials due to nanoscale nonlocal size effects: Green’s function
solutions and embedded inclusions », Phys. Rev. B, vol. 74, no 1, juill. 2006.
[23] M. S. Majdoub, R. Maranganti, et P. Sharma, « Understanding the origins of the
intrinsic dead layer effect in nanocapacitors », Phys. Rev. B, vol. 79, no 11, mars 2009.
[24] G. Maugin, « The method of virtual power in continuum mechanics: application to
coupled fields », Acta Mech., vol. 35, no 1, p. 1–70, 1980.
[25] A. K. Tagantsev et A. S. Yurkov, « Flexoelectric effect in finite samples », J. Appl.
Phys., vol. 112, no 4, p. 044103, août 2012.
[26] S. Shen et S. Hu, « A theory of flexoelectricity with surface effect for elastic
dielectrics », J. Mech. Phys. Solids, vol. 58, no 5, p. 665‑677, mai 2010.
[27] M. E. Gurtin et A. I. Murdoch, « A continuum theory of elastic material surfaces »,
Arch. Ration. Mech. Anal., vol. 57, no 4, p. 291–323, 1975.
[28] E. A. Eliseev, A. N. Morozovska, M. D. Glinchuk, et R. Blinc, « Spontaneous
flexoelectric/flexomagnetic effect in nanoferroics », Phys. Rev. B, vol. 79, no 16, avr.
2009.
[29] A. S. Yurkov, « Elastic boundary conditions in the presence of the flexoelectric
effect », JETP Lett., vol. 94, no 6, p. 455‑458, nov. 2011.
[30] M. Stengel, « Unified ab initio formulation of flexoelectricity and strain-gradient
elasticity », Phys. Rev. B, vol. 93, no 24, juin 2016.
[31] M. Stengel, « Surface control of flexoelectricity », Phys. Rev. B, vol. 90, no 20, p.
201112, 2014.
[32] M. Stengel, « Microscopic response to inhomogeneous deformations in curvilinear
coordinates », Nat. Commun., vol. 4, nov. 2013.
[33] M. Stengel, « Flexoelectricity from density-functional perturbation theory », Phys.
Rev. B, vol. 88, no 17, p. 174106, 2013.
29
[34] R. Resta, « Towards a bulk theory of flexoelectricity », Phys. Rev. Lett., vol. 105, no
12, p. 127601, 2010.
[35] J. Hong et D. Vanderbilt, « First-principles theory of frozen-ion flexoelectricity »,
Phys. Rev. B, vol. 84, no 18, p. 180101, 2011.
[36] M. Stengel et D. Vanderbilt, « First-principles theory of flexoelectricity »,
Flexoelectricity Solids Theory Appl., p. 31, 2016.
[37] J. Hong et D. Vanderbilt, « First-principles theory and calculation of flexoelectricity »,
Phys. Rev. B, vol. 88, no 17, p. 174107, 2013.
[38] J. Hong, G. Catalan, J. F. Scott, et E. Artacho, « The flexoelectricity of barium and
strontium titanates from first principles », J. Phys. Condens. Matter, vol. 22, no 11, p.
112201, mars 2010.
[39] I. Ponomareva, A. K. Tagantsev, et L. Bellaiche, « Finite-temperature flexoelectricity
in ferroelectric thin films from first principles », Phys. Rev. B, vol. 85, no 10, mars 2012.
[40] T. Dumitrică, C. M. Landis, et B. I. Yakobson, « Curvature-induced polarization in
carbon nanoshells », Chem. Phys. Lett., vol. 360, no 1, p. 182–188, 2002.
[41] A. G. Kvashnin, P. B. Sorokin, et B. I. Yakobson, « Flexoelectricity in Carbon
Nanostructures: Nanotubes, Fullerenes, and Nanocones », J. Phys. Chem. Lett., vol. 6, no
14, p. 2740‑2744, juill. 2015.
[42] S. V. Kalinin et V. Meunier, « Electronic flexoelectricity in low-dimensional
systems », Phys. Rev. B, vol. 77, no 3, janv. 2008.
[43] S. Chandratre et P. Sharma, « Coaxing graphene to be piezoelectric », Appl. Phys.
Lett., vol. 100, no 2, p. 023114, 2012.
[44] M. S. Majdoub, P. Sharma, et T. Cagin, « Enhanced size-dependent piezoelectricity
and elasticity in nanostructures due to the flexoelectric effect », Phys. Rev. B, vol. 77, no
12, mars 2008.
[45] M. S. Majdoub, P. Sharma, et T. Çağin, « Dramatic enhancement in energy harvesting
for a narrow range of dimensions in piezoelectric nanostructures », Phys. Rev. B, vol. 78,
no 12, sept. 2008.
[46] C. Liu, S. Hu, et S. Shen, « Effect of flexoelectricity on electrostatic potential in a bent
piezoelectric nanowire », Smart Mater. Struct., vol. 21, no 11, p. 115024, nov. 2012.
[47] C. Liu, H. Wu, et J. Wang, « Giant piezoelectric response in piezoelectric/dielectric
superlattices due to flexoelectric effect », Appl. Phys. Lett., vol. 109, no 19, p. 192901,
nov. 2016.
[48] E. Fukada, G. M. Sessler, J. E. West, A. Berraissoul, et P. Günther, « Bending
piezoelectricity in monomorph polymer films », J. Appl. Phys., vol. 62, no 9, p. 3643‑
3646, nov. 1987.
[49] W. Ma et L. E. Cross, « Large flexoelectric polarization in ceramic lead magnesium
niobate », Appl. Phys. Lett., vol. 79, no 26, p. 4420, 2001.
[50] W. Ma et L. E. Cross, « Observation of the flexoelectric effect in relaxor Pb(Mg[sub
1/3]Nb[sub 2/3])O[sub 3] ceramics », Appl. Phys. Lett., vol. 78, no 19, p. 2920, 2001.
30
[51] W. Ma et L. E. Cross, « Flexoelectricity of barium titanate », Appl. Phys. Lett., vol. 88,
no 23, p. 232902, 2006.
[52] W. Ma et L. E. Cross, « Flexoelectric effect in ceramic lead zirconate titanate », Appl.
Phys. Lett., vol. 86, no 7, p. 072905, 2005.
[53] W. Ma et L. E. Cross, « Flexoelectric polarization of barium strontium titanate in the
paraelectric state », Appl. Phys. Lett., vol. 81, no 18, p. 3440, 2002.
[54] W. Ma et L. E. Cross, « Large flexoelectric polarization in ceramic lead magnesium
niobate », Appl. Phys. Lett., vol. 79, no 26, p. 4420, 2001.
[55] P. Zubko, G. Catalan, A. Buckley, P. R. L. Welche, et J. F. Scott, « Strain-Gradient-
Induced Polarization in SrTiO 3 Single Crystals », Phys. Rev. Lett., vol. 99, no 16, oct.
2007.
[56] M. Marvan et A. Havránek, « Flexoelectric effect in elastomers », in Relationships of
Polymeric Structure and Properties, Springer, 1988, p. 33–36.
[57] S. Baskaran, X. He, Q. Chen, et J. Y. Fu, « Experimental studies on the direct
flexoelectric effect in α-phase polyvinylidene fluoride films », Appl. Phys. Lett., vol. 98,
no 24, p. 242901, juin 2011.
[58] S. Baskaran et al., « Giant flexoelectricity in polyvinylidene fluoride films », Phys.
Lett. A, vol. 375, no 20, p. 2082‑2084, mai 2011.
[59] S. Baskaran, X. He, Y. Wang, et J. Y. Fu, « Strain gradient induced electric
polarization in α -phase polyvinylidene fluoride films under bending conditions », J. Appl.
Phys., vol. 111, no 1, p. 014109, janv. 2012.
[60] S. Zhang et al., « Investigation of the 2312 flexoelectric coefficient component of
polyvinylidene fluoride: Deduction, simulation, and mensuration », Sci. Rep., vol. 7, no 1,
déc. 2017.
[61] S. Zhang, M. Xu, X. Liang, et S. Shen, « Shear flexoelectric coefficient μ 1211 in
polyvinylidene fluoride », J. Appl. Phys., vol. 117, no 20, p. 204102, mai 2015.
[62] S. Zhang, X. Liang, M. Xu, B. Feng, et S. Shen, « Shear flexoelectric response along
3121 direction in polyvinylidene fluoride », Appl. Phys. Lett., vol. 107, no 14, p. 142902,
oct. 2015.
[63] L. E. Cross, « Flexoelectric effects: Charge separation in insulating solids subjected to
elastic strain gradients », J. Mater. Sci., vol. 41, no 1, p. 53‑63, janv. 2006.
[64] J. Y. Fu, W. Zhu, N. Li, et L. E. Cross, « Experimental studies of the converse
flexoelectric effect induced by inhomogeneous electric field in a barium strontium titanate
composition », J. Appl. Phys., vol. 100, no 2, p. 024112, juill. 2006.
[65] P. Hana, M. Marvan, L. Burianova, S. J. Zhang, E. Furman, et T. R. Shrout, « Study of
the Inverse Flexoelectric Phenomena in Ceramic Lead Magnesium Niobate-Lead
Titanate », Ferroelectrics, vol. 336, no 1, p. 137‑144, juill. 2006.
[66] A. K. Tagantsev, E. Courtens, et L. Arzel, « Prediction of a low-temperature
ferroelectric instability in antiphase domain boundaries of strontium titanate », Phys. Rev.
B, vol. 64, no 22, nov. 2001.
31
[67] R. Maranganti et P. Sharma, « Atomistic determination of flexoelectric properties of
crystalline dielectrics », Phys. Rev. B, vol. 80, no 5, août 2009.
[68] G. Catalan, L. J. Sinnamon, et J. M. Gregg, « The effect of flexoelectricity on the
dielectric properties of inhomogeneously strained ferroelectric thin films », J. Phys.
Condens. Matter, vol. 16, no 13, p. 2253, 2004.
[69] G. Catalan, B. Noheda, J. McAneney, L. J. Sinnamon, et J. M. Gregg, « Strain
gradients in epitaxial ferroelectrics », Phys. Rev. B, vol. 72, no 2, juill. 2005.
[70] J. Fousek, L. E. Cross, et D. B. Litvin, « Possible piezoelectric composites based on
the flexoelectric effect », Mater. Lett., vol. 39, no 5, p. 287–291, 1999.
[71] J. Y. Fu, W. Zhu, N. Li, N. B. Smith, et L. Eric Cross, « Gradient scaling phenomenon
in microsize flexoelectric piezoelectric composites », Appl. Phys. Lett., vol. 91, no 18, p.
182910, 2007.
[72] N. D. Sharma, R. Maranganti, et P. Sharma, « On the possibility of piezoelectric
nanocomposites without using piezoelectric materials », J. Mech. Phys. Solids, vol. 55, no
11, p. 2328‑2350, nov. 2007.
[73] B. Chu, W. Zhu, N. Li, et L. E. Cross, « Flexure mode flexoelectric piezoelectric
composites », J. Appl. Phys., vol. 106, no 10, p. 104109, 2009.
[74] A. Gruverman et al., « Mechanical stress effect on imprint behavior of integrated
ferroelectric capacitors », Appl. Phys. Lett., vol. 83, no 4, p. 728‑730, juill. 2003.
[75] H. Lu et al., « Mechanical Writing of Ferroelectric Polarization », Science, vol. 336, no
6077, p. 59‑61, avr. 2012.
[76] Q. Deng, M. Kammoun, A. Erturk, et P. Sharma, « Nanoscale flexoelectric energy
harvesting », Int. J. Solids Struct., vol. 51, no 18, p. 3218‑3225, sept. 2014.
[77] A. G. Moura et A. Erturk, « Electroelastodynamics of flexoelectric energy conversion
and harvesting in elastic dielectrics », J. Appl. Phys., vol. 121, no 6, p. 064110, févr. 2017.
[78] J. K. Han et al., « Nanogenerators consisting of direct-grown piezoelectrics on multi-
walled carbon nanotubes using flexoelectric effects », Sci. Rep., vol. 6, no 1, sept. 2016.
32
Chapitre 2
polarisation (voir le travail de Midlin [1] qui utilise le principe d’Hamilton) et le gradient de
déformation (voir [2], [3], [4] en Russe, traduit en anglais dans [5], [6], [7] et plus récemment,
[8] [9]). La théorie « quasi-linéaire » est en quelque sorte rudimentaire quant à sa base
conceptuelle. Elle ne permet pas de traiter de façon consistante la réponse non linéaire de la
structure, due aux grandes déformations ou encore le comportement dissipatif et couplé non
l’électromagnétisme des milieux continus peut être obtenue à partir d’une méthode
énergétique appelée « principe des puissances virtuelles » [10]–[21]. En effet, Maugin montre
que cette méthode est particulièrement bien adaptée pour le couplage entre les interactions des
ont appliqué la méthode des puissances virtuelles en prenant en compte des gradients de
méthode avec des gradients de déformations [15]. Daher et Maugin ont appliqué la méthode
aux matériaux semi-conducteurs présentant des discontinuités et des interfaces [19]. Cette
formulation théorique décrit le couplage entre la mécanique, l’électromagnétisme et la semi-
conduction des milieux continus en présence de surfaces singulières. La théorie est basée sur
la notion de dualité, qui est bien adaptée aux principes de symétrie et d’invariance tel le
principe d’objectivité (invariance par translation et rotation) des équations décrivant l’état du
système. C’est ainsi que la symétrie du tenseur des contraintes ij ji est assurée par le
quantité de mouvement.
34
Nous allons étendre les équations données dans [19] en ajoutant les gradients de
principe des puissances virtuelles mène aux équations d’équilibre dans le volume et les
charges.
section 2 de ce chapitre et dans le glossaire (voir annexe C). La section 3 est consacrée à
l’objectivité qui est primordiale dans la méthode des puissances virtuelles. Les formes
conduction décrivant les différents continua (avec des charges de différentes espèces). Les
déduites dans la section 6 après avoir utilisé le principe des puissances virtuelles, ainsi que le
formes générales sont établies dans la section 7, où une comparaison de notre travail sera faite
35
2 Notations globales
fig 13. Schéma du domaine avec comme contour G et comme vecteur normal n à la
surface (qui peut être discontinue)
vers l’extérieur de . La vitesse de chaque point, dans le référentiel Galiléen fixe RG (le
référentiel du laboratoire), est notée v . Le champ des vitesses v , le gradient de la
36
v vi
t fixed X
(2)
i
X K
F xi , K F 1 X K ,i
X K t
xi
t
(3)
1
D Dij v(i , j ) (vi , j v j ,i ) D ji (4)
2
1
ij vi , j (vi , j v j ,i ) ji
2
1
i ijk jk jk jkp p
2
xi , K X K , j ij X K ,i xi , L KL (6)
et J det(F ) 0 (7)
ìï ¶ üï ìï ¶ üï
Ñ = í ; i = 1,2,3ý ÑR = í ; K =1,2,3ý (8)
ïî ¶xi ïþ ïî ¶X K ïþ
i
G ∇F xi , KL xk , LK (9)
X K X L t
37
La divergence d’un tenseur d’ordre 2 est prise par rapport au dernier indice :
A A A (11)
Le même symbole est aussi utilisé pour la ligne G (voir annexe A, eqn (A6)).
0 J 0 ou (.v ) 0 (12)
dY Y
Y v. Y (13)
dt t
3 Objectivité
3.1 Définitions
observe doivent rester les mêmes. Les mesures effectuées dans un repère de référence sont
suffisantes pour déterminer ces propriétés physiques dans toutes les configurations en
mouvement rigide (translation ou rotation) les unes par rapport aux autres. Contrairement aux
autres types d’effort, les forces intérieures, pour lesquelles nous sommes amenés à construire
Soit R’ un repère en mouvement de corps rigide par rapport à un autre repère R. Soit
38
xi' (t ') Qij (t )x j bi (t ) t ' t t0
(14)
QikQjk QkiQkj ij det(Q) = 1
Définition 1
On dit que les mouvements x et x’ sont équivalents si et seulement si :
Définition 2
Toute quantité tensorielle est objective si pour des mouvements objectivement équivalents
elle obéit à la loi de transformation tensorielle appropriée quel que soit le temps.
S '( X , t ') S( X , t )
Sk' ( X , t ') Qkl (t )Sl ( X , t ) (15)
Skl' ( X , t ') Qkm (t )Qln (t )Smn ( X , t )
Nous allons montrer maintenant que le gradient des vitesses vK , L est objectif à partir de la
On a donc directement :
dvm dx
vk ,l ' Qkm Qkm m
dx 'l dx 'l
dvm dxn dx
Qkm Qkm m (18)
dxn dx 'l dx 'l
Qkm vm,nQln QkmQlm
On a ainsi :
39
Dkl ' QkmQln Dkl
avec Dkl v k ,l et kl vk ,l
kl ' QkmQln kl QkmQlm (19)
dvm,n
vk ,lp ' QkmQln
dx p '
dvm,n dxr
QkmQln (20)
dxr dx ' p
QkmQln Q pr vm ,nr
L’équation (17) nous montre sur l’exemple de la vitesse que la dérivée temporelle
classique n’est pas objective. Il faut donc introduire (au moins) une nouvelle dérivée par
rapport au temps.
Pour des raisons expliquées plus tard, on introduit une dérivée spécifique pour un gradient de
vecteur définie de façon similaire [12][14][15][16] :
D’autre part, la dérivée convective, notée DC , est définie de la manière suivante pour un
40
A DC A A A . v A . v (24)
.
aij DC a ij ai , j ai ,k vk , j vi ,k ak , j ai , j vi ,k ak , j (25)
Finalement, on peut écrire la relation entre les différentes dérivées du gradient de vecteurs
en utilisant (23) et (25), ce qui donne la relation :
Notons pour finir que ces dérivées ne sont pas les seules dérivées objectives. Nous les
utilisons car elles simplifieront les expressions par la suite, mais il en existe une infinité.
corps :
- L’approche vectorielle qui utilise la notion de vecteur pour représenter des efforts.
- L’approche énergétique qui détermine les forces qui s’exercent sur un système à partir
des énergies mises en jeu.
voulons connaitre le poids d’un objet, nous le soupesons, si nous voulons connaitre la rigidité
d’un ressort, nous l’étirons, etc… Dans chaque cas, nous faisons « travailler » l’objet étudié.
De plus, la méthode se révèle être d’une grande souplesse d’emploi. Selon que l’on choisira
un espace vectoriel de vitesses plus ou moins vaste, on aura par dualité une description des
efforts plus ou moins fine. Dans ce qui suit, nous allons comparer les 2 méthodes dans un cas
41
L’approche vectorielle consiste à introduire le tenseur des contraintes afin de décrire
les forces d’interactions exercées en un point du milieu. On utilise le théorème de Cauchy qui
montre qu’en tout point M d’un solide, la dépendance du vecteur contrainte T par rapport à n
est linéaire. Il existe donc un champ tensoriel du second ordre tel que pour tout point M,
on a :
Ti (n) ij n j (28)
f d T da 0
i
i (29)
x f d x T da 0
(30)
ij , j fi 0 (31)
ij ji (32)
ij v ij* ij v*i , j ij v*i , j doit se réduire à ij v*i , j ij Dij* , puisque ij vi , j ne fait pas partie
des variables objectives. On doit donc avoir [ ij ] 0 , c’est à dire ij ji . On postule ensuite
42
que la puissance virtuelle des déformations est égale à la puissance des efforts extérieurs quel
D d fi vi* dv T v
* *
ij ij i i da (33)
De plus, on a :
fi vi* d T n v da 0
*
ij , j i ij j i (35)
Cet exemple nous montre l’importance du mouvement virtuel : en supposant que l’équation
ci-dessus est vraie quelle que soit v i* , on en déduit directement que chacun des deux contenus
de parenthèses doit être nul indépendamment de l’autre, ce qui redonne l’équation d’équilibre
(31) et l’équation (28) .
classique se dégage clairement pour des milieux plus complexes qui peuvent être polarisables,
indépendantes qui ne sont pas de nature mécanique. Nous les expliciterons dans la section 5,
mais nous rappelons ici que même s’il n’existe pas de règles systématiques qui nous les
43
2) Axiome d’équiprésence : cet axiome est une mesure de précaution puisqu’il dit que
toutes les lois de comportement doivent dépendre de la même liste de variables
constitutives indépendantes. Certains critères tels que l’axiome d’admissibilité peuvent
ultérieurement éliminer la dépendance par rapport à certaines variables.
3) Axiome de déterminisme : le comportement matériel ne dépend pas des points
extérieurs au corps ni des évènements futurs. C’est un axiome d’exclusion.
4) Axiome d’objectivité : les lois de comportement doivent être indépendantes du
mouvement de l’observateur.
5) Axiome d’invariance matérielle : les orientations cristallographiques des points
matériels dans un corps donnent lieu à des symétries dans les propriétés matérielles ;
Ceci impose des restrictions sur les lois de comportement.
6) Axiome de localisation spatiale : les lois de comportement sont influencées seulement
par le voisinage proche d’une position donnée.
7) Axiome de mémoire : ce principe est l’analogue du principe de localisation spatiale
dans le domaine du temps.
8) Axiome d’admissibilité : toutes les lois de comportement ne doivent en aucune
manière violer les équations de base de la mécanique des milieux continus.
D
H J . D q f (37)
t
44
tel 0 0 c 1 (où c représente la vitesse de la lumière dans le vide). En utilisant les équations
2
q f
. J 0 (38)
t
B
P D 0 E M H (39)
0
Galiléenne, à o(v 2 ) près, on peut écrire les relations de passage entre les grandeurs dans RG et
dans RC ( x, t ) :
E vB (40)
B
v 0 E
0 0 (41)
H vD
0 (42)
J qf v (43)
P M vP (44)
De même, les équations (36) et (37) peuvent être écrites dans le référentiel en co-
mouvement RC (x,t) . Il vient alors :
45
*
B0 . B 0 (45)
*
D . D q f (46)
charge peut être considérée comme étant une vitesse. Par dualité (recherche d’une quantité
objective qui par produit scalaire avec la partie objective du champ des vitesses donnera une
puissance des efforts intérieurs invariante par tout mouvement rigidifiant), les forces internes
liées seront construites. On aura alors le schéma de la Fig 2 traduisant les différentes
interactions [16]. Les continua de charges , peuvent correspondre à des électrons, trous, ions,
impuretés, etc. donne le type de charge. Nous avons les relations suivantes:
fig 14. Schéma des interactions entre déformations, semi conduction et thermique [19]
q f qf J J (47)
J
v (48)
qf
u v v
qf (49)
46
où qf , v , v et u donnent respectivement la densité volumique des charges électriques de
type , le champ des vitesses décrivant le ième continuum, la déformation matérielle globale
et la vitesse relative. représente la densité de courant électrique des charges de type
qui diffusent dans le matériau. De même que pour la conservation de la masse (voir par
exemple [19]), on peut écrire pour chaque continuum , une équation de conservation de la
charge qui prenne en compte les possibilités de recombinaisons entre électrons et trous ou
inversement la possibilité de génération de paires :
d
dt
qf d r d (50)
r 0 (51)
et où, pour tout scalaire, vecteur ou tenseur Y, la dérivée par rapport au temps caractéristique
du mouvement du continuum est définie par :
d Y Y
v .Y (52)
dt t
L’équation (50) nous donne donc aussi la conservation locale des charges de type :
d qf
qf .v r (53)
dt
En appliquant (52) avec Y qf , on peut alors vérifier facilement que les équations (48) et
qf
. J r (54)
t
Notons qu’en sommant sur et en utilisant les eqns (47) et (51), l’équation (54) redonne
bien l’équation (38) valable pour la totalité des types de charges.
47
Notons pour finir ce paragraphe que pour les calculs ultérieurs, il est utile de réécrire
q eff
. E (55)
0
*
0 E eff
(56)
0
eff eff
où la charge effective q et le courant effectif dans RC sont donnés par:
. P
q eff qf . P q ( eff ) q (eff) qf 1 (57)
. D
*
*
. D
( eff ) 1
eff
P (eff)
qf P (58)
avec
(59)
Nous donnons ici l’expression des forces et des couples volumiques agissant sur la
matière déformable, considérée comme continue, dues aux champs électromagnétiques. Nous
notons toutefois qu’ils peuvent être calculés à partir de moyennes sur des ensembles de micro-
f em qeff eff
B P B . (60)
48
em
La densité volumique de couple pondéromoteur d’origine électromagnétique c est quant à
elle un pseudo-vecteur, puisqu’elle s’exprime comme une somme de produits vectoriels. Elle
est ainsi issue d’un tenseur d’ordre 2, dont nous verrons plus tard qu’il correspond à la partie
antisymétrique du tenseur des contraintes électromagnétiques.
em
Théorème : il existe au moins un tenseur (des contraintes) t ij et un vecteur (densité
volumique d’impulsion généralisée du champ électromagnétique) G, définis à partir des
champs électromagnétiques, tels que :
G Cijem tem
f em t em , ij (62)
t
f em L f t em (63)
G
f q eff eff
B t F (64)
t
L
1 B2
Bi B j
t 0 Ei E j
F
0 E ij
2
(66)
0 2 0
ij
tijem i Pj i Bj B ij (67)
49
F
Notons que compte tenu de leur définition respective et notamment du fait que tij est
manifestement symétrique, on a :
Cijem tem
ij
tijem (68)
em
En ce qui concerne la densité de contrainte électromagnétique de surface T , le
em
théorème de la divergence appliqué à (61), permet de définir T sur chaque surface par
[18][20]:
Ti em tijem Gi v j n j (69)
. P
1 v v eff B (70)
. D
avec
P
v eff
v (71)
.D
Compte tenu de ces définitions, on peut maintenant réécrire la force de Lorentz L f comme
[19] :
L f f avec f qf
(72)
On peut aussi noter qu’en l’absence de polarisation et d’aimantation, l’expression (70) peut
s’écrire sous la forme classique :
P 0; 0
E v B (73)
50
5 Equations d’équilibre électro-magnéto-thermo-mécanique
peuvent être déduites à partir de trois principes écrits dans leur forme globale pour un volume
[16], [20] Ces principes sont le principe des puissances virtuelles et le premier et second
principe de la thermodynamique.
Qh , , les puissances des forces virtuelles, des forces internes, des forces externes sur
des forces externes de contact et des forces « prescrites», l’énergie cinétique, l’énergie
puissances virtuelles, on construit un espace vectoriel des mouvements généralisés, puis par
dérivation par rapport au temps un espace des vitesses généralisées, puis après utilisation du
généralisées.
systématique faisant appel aux notions de dualité et d’application linéaire sur un espace
vectoriel et au principe d’objectivité. Nous commençons par le choix de l’espace vectoriel des
vitesses selon le problème considéré. Lorsque est fixé, les efforts extérieurs sont
introduits par dualité et forment un espace appelé , dual de c’est-à-dire que le produit
d’un élément de et celui d’un élément de donne un nombre, ici une densité volumique
Nous voulons étudier la déformation d’un milieu polarisable et magnétisable qui peut
51
xi , xi , i , i (74)
L’espace des champs de vitesses généralisées correspondant (en dérivant par rapport à t) est
alors défini en tout point de par :
(0)
vi , vi , i , i (75)
Pi
avec i et i i
(76)
i est donc la polarisation par unité de masse dans Rc , tandis que i est l’aimantation par
Pour une théorie du second gradient et donc pour avoir une meilleure description des
obtient :
vi , vi , j , vi , jk , vi , vi, j , i , i , j , i , i , j (77)
vi , Dij , ij , vi , jk , ui , Dij , ij , i , i , j , i , i , j (78)
avec vi , j vi , j v i , j Dij ij et vi, j vi , j v i , j Dij ij (79)
volumes mesoscopiques (petits à l’échelle du système, grand à l’échelle des atomes) des
interactions entre composants microscopiques. C’est à partir de ces forces que l’on construira
les lois de comportement qui devront être obligatoirement objectives (c’est-à-dire invariantes
52
par translation ou rotation du trièdre de référence). Les vitesses généralisées qui sont les
variables duales des forces dans l’expression des puissances doivent donc être elles aussi
objectives en vertu de l’invariance du produit scalaire. Nous construirons les forces intérieures
par dualité avec le champ de vitesses généralisées dont nous ne garderons que la partie
objective et les interprèterons avec une analyse dimensionnelle. Pour construire le champ de
vitesse objectif, on va choisir une base linéairement indépendante qui est incluse dans . On
va pour cela utiliser les dérivées de Jaumann (21) pour les champs électromagnétiques, ainsi
que les dérivées temporelles (23) pour le champ de vitesses globales et pour les vitesses des
Les dérivées de Jaumann, ̂i et ˆ i , ̂ij et ˆ ij , sont introduites de la même façon. L’espace
objectif obj est alors composé à l’aide d’un ensemble de variables objectives et
indépendantes.
obj Dij , vi , jk , ui , Dij , ˆi , ˆij , ˆi , ˆij , ˆi , ˆij , ˆi , ˆij
Cette base de obj n’est pas unique, mais permet de simplifier les expressions, notamment
pour le passage du référentiel Lagrangien au référentiel eulérien. Nous notons en effet que les
deux dérivées convectives (24) et (25) nous permettent d’introduire d’autres variables
objectives comme, P , , ou les tenseurs, et , qui seront utilisées dans la section 6.
Les champs objectifs, éléments de obj , peuvent maintenant être construits par
dualité avec les éléments correspondants de obj afin de construire les forces internes
généralisées intervenant dans le calcul de la puissance des efforts intérieurs [16] :
53
obj ij ji , ijk , L i , ij ji , L Ei , L ij , L Bi , L ij , L Ei , L ij , L Bi , L
ij
ij et ij représentent respectivement, le tenseur des contraintes intrinsèques du premier
tenseur intrinsèque des contraintes du deuxième ordre. Une analyse dimensionnelle montre
que
L
Bi et
L
Bia sont des champs d’induction magnétique qui reflètent
phénoménologiquement les interactions entre les spins et les mouvements du réseau
L
cristallin, tandis que i , L Ei et L Ei sont homogènes à des champs électriques qui reflètent
phénoménologiquement les interactions entre les multipôles provenant de la polarisation du
L
matériau et les mouvements du réseau cristallin. L’apparition des termes L ij , L
ij , ij et L
ij
voisins et ne sont pris en compte que pour les matériaux ferromagnétiques [16].
l’aimantation si le matériau est magnétique, d’où la présence de qui est pris égal à 1 si le
t
(a)
(84)
Gi v j n j vi*da
2
Le rapport gyromagnétique de l’électron est 𝛾𝑒 = 𝑔𝐿𝑎𝑛𝑑é 𝜇𝐵 /ℏ, où 𝑔𝐿𝑎𝑛𝑑é est un nombre appelé facteur de
Landé fonction des nombres quantiques correspondants aux divers opérateurs de moments cinétiques concernés
et 𝜇𝐵 = 𝑒ℏ⁄2𝑚𝑒 est appelé magnéton de Bohr.
54
vitesse angulaire virtuelle de précession du moment magnétique autour du champ
magnétique3.
Les forces internes qui correspondent aux interactions atomiques dans les réseaux
correspondante [16][20].
P
(i ) ( , V
obj ) p(i ) d (85)
avec
où, b est égal à 1 lorsque et jouent un rôle important (i.e quand le matériau a une
susceptibilité magnétique grande devant 1). ˆi* , ˆi , ˆij* , ˆij , ˆi * , ˆi , ˆij * et ˆ ij représentent
Après avoir développé la puissance virtuelle des forces internes (voir annexe B, [10][15]), on
trouve :
3
i* est telle que avec B pour le mouvement réel du moment cinétique orbital. Si on
considère alors comme un rotateur rigide, l’égalité de la résultante dynamique avec le somme des moments
(dont les couples) pour chaque particule individuelle, se traduit par le fait que cette résultante dynamique est
égale à B 1 , ce qui contribue à justifier la forme présente dans P (a) lorsque tous les moments
cinétiques des atomes du volume mésoscopique sont alignés dans le même sens (ferromagnétisme).
55
P
(i ) ( , V
obj
) ij ijk ,k L E[i j L B[i j L [i k j ,k L
[i k
j ,k vi*d ,j
L Ei L ij , j i L Bi L
ij , j
i
d
L
ij n j i L
ij n j i da
ij ijk ,k L E[ i L B[ i L
j j [i k j , k
L
[i k
j ,k n j vi*da
j
ˆ n
j
ˆ n
p p ijk nk vi*da
ijk n j nk
vi*
n
da ijk nk b j vi*ds
(87)
ij L E[i j L B[i j L [ik j ,k L
[i k
j ,k ,j
vi * qf
L *
i i
u d
E L
i
L
ij , j
i
L
B
i
L
ij , j
i
d
L
ij n j i L
ij n j i da
ij
L
E[i j L B[i j L [ik j ,k L
[i k
j ,k n j vi *da
où
Ai B j
1
2
Ai B j Aj Bi et Ci k B j
1
2
Cik B j C jk Bi (88)
fi représente les forces volumiques comme la gravitation. ij et Cij sont respectivement les
56
ij tem
ij
Cij tem
ij
(90)
Ki i Li Bi
em
où le tenseur t ij est le tenseur des contraintes électromagnétiques introduit en (62). K ij et Lij
L L
ont la même dimension que ij et ij . Comme il n’existe (pour le moment) aucun support
et où f sont des forces volumiques sur les charges faisant partie du continuum . [16].
P
(v) ( , V
) fi vi tem D* tem
ij ij
*
ij ij i i Bi i fi ui * d (92)
0 i i
i i
Qi
présence du facteur . Le terme est appelé le vecteur « traction de polarisation
0
électronique »
57
Notons que la puissance virtuelle des forces de contacts peut s’écrire d’une autre
Ti Ti em Ti Ti em ˆ j n j ˆ p n p Rij sur
où est le vecteur tangent unitaire à la ligne G , défini dans l’annexe A et eijk correspond au
tenseur de Levi-Civita (voir annexe C).
vi* Qi
da
(0)
P ( , V ) Ti vi Ri
(c)
n 0 i
(95)
0 i Ti ui
i
*
da Li vi ds
*
Les principes
système étudié, nous avons besoin de trois principes : le principe des puissances virtuelles, le
Dans le référentiel Galiléen, la puissance virtuelle des forces inertielles doit être
égale à la somme des puissances virtuelles des forces internes, externes et de contacts pour
58
Premier principe de la thermodynamique
La variation de l’énergie totale doit être égale à la somme des puissances des forces
exercées sur le système et du taux de chaleur reçue par le système Q h , [16].
particulaire de l’entropie dans le domaine doit être supérieure à l’entropie créée dans .
Mathématiquement, on a :
d
N( ) ( ) (98)
dt
1 1
K ( ) v 2 d 2 d
2 2 (99)
E ( ) d (100)
1 B2
U ( ) 0E
em 2
2 d (101)
2 0
N ( ) d (102)
Qh ( ) hd q n da (103)
59
( ) d n da (104)
où d est l’inertie de polarisation électrique, est l’énergie interne par unité massique, est
l’entropie par unité de masse, h est le taux de chaleur par unité de masse fournie par
l’extérieur à et q est le vecteur courant de chaleur représentant le taux de la chaleur
surfacique reçue, i.e. la somme des puissances des flux thermiques et du vecteur de
q q (105)
(106)
applique le principe des puissances virtuelles pour un champ de vitesses réel. On obtient:
d
K ( ) U em ( ) P( e ) ( ) P(i ) ( , obj )
dt
(107)
f . u d n da
vi Qi
P( e ) ( ) fi vi d Ti vi Ri da
n 0 i
(108)
0 i i Ti ui da Li vi ds
et P(i ) ( , obj ) donne la puissance des forces internes pour des vitesses réelles.
60
Dans les équations (86)-(92) et (108), f et T sont respectivement les forces volumiques et
étant le tenseur d’interaction des spins qui traduit les interactions entre spins magnétiques
voisins. (cf. thèse B. Collet [24])
Pour les champs virtuels v*, v *, *, * et v * n agissant sur , , et on
obtient les équations de mouvement à partir des équations (96) et (62) qui donnent les égalités
entre les mouvements et les interactions dans un milieu déformable, magnétisable, polarisable
et semi-conducteur :
v
n
61
v
tij , j fi fi em qf
i L i vi dans (110)
tij ij ijk ,k L
E[i j L B[i j L
[i k
j ,k L
[i k
j ,k (113)
v
dans (115)
avec:
(117)
dans (118)
62
sur (119)
i
eff
i L EiT (120)
L
EiT L Ei L Ei et L T
ij L
ij L ij (121)
i
L T
Bieff tel que dans (122)
ij , j
ipq L T
pj n j 0 p q 0 sur (123)
L T
ij , j Bieff 0 dans (124)
L T
ij n j 0 i sur (125)
• Si b ' = b = 0 , sinon
Avec l’induction magnétique effective B eff et le champ local d’induction du premier ordre
L
BT et du deuxième ordre L T
ij définie par:
63
Bieff Bi L BiT (127)
L
BiT L Bi L Bi et L
ij
T
L
ij L
ij (128)
Et utilisant les équations (115) et (116), on peut réécrire les équations (110) and (111) d’une
façon plus conventionnelle:
puissances virtuelles (où v * coïncide avec la vitesse réelle v ) (107), On obtient le théorème
de l’énergie :
E ( ) P( i ) ( , obj ) Q h ( ) Qem ( ) (132)
où on a :
avec qem
. f . u (134)
64
qi
h
et i (136)
À partir de l’équation (136), l’équation locale correspondant au deuxième principe de la
thermodynamique (98) est donnée par :
d
h . q . (137)
dt
libre de Helmholtz y :
(138)
On obtient ainsi l’inégalité de Clausius-Duhem dans sa forme locale à partir de (135) et (137):
d d
p(i ) qem . 0 (139)
dt dt
D’autre part, avec (86), (26), (27), (72), (49) et (115), on obtient:
p(i ) qem tijT Dij ijk vi , jk L EiT Pi L BiT i L ijT ij L
ij
T
ij
(140)
tij ui
i
,j i
a
où et t sont des tenseurs symétriques:
, ij and ij ont été définis dans (121), (128), (24) avec
TL L T L T L T
où Ei , Bi , ik , ik , Pi , i
65
Il est difficile de reconnaitre la théorie des milieux déformables semi-conducteurs à
partir des équations de la section 5 et de l’équation (139). Afin d’essayer d’en clarifier le sens,
on définit alors :
Le tenseur M ij qui a la dimension d’un potentiel chimique va être réduit comme en [19], ce
M ij ij (144)
(eff)
(146)
De plus, on doit prendre en compte les concentrations des types de charges. (voir Kiréev
particule par unité de volume dans : la fraction massique c est donnée par :
n m
c (147)
qf n q (148)
a
On peut définir la charge électrique massique pour le type , noté cl , tel que:
66
c m
c
avec (149)
q
Donc, la densité volumique des charges électriques q f donne:
On obtient ainsi l’équation de conservation de la charge pour le type de charge . (53) peut
s’écrire alors:
dc
. r (151)
dt
Compte tenu de ces nouvelles définitions, de et d’après les équations (145) et
(151), on peut réécrire l’inégalité de Clausius-Duhem (139) sous la forme suivante :
d d T
tij Dij ijk vi , jk Ei Pi
L T
dt dt
(152)
dc 1
B
L T
L T
ij L T
ij (eff)
r q . 0
i i ik ik i i
dt
(eff)
, ij et ij sont définis par les équations (141),
T L T L T L T L T
où tij , Ei , Bi , ij , ij , i , Pi , i
L’inégalité (152) qui reflète la contrainte thermodynamique (le processus doit être
thermodynamiquement admissible) est importante pour la construction des équations
constitutives (ou lois de comportement).
67
TKL J X K ,i X L, j tijT
TKLM J X K ,i X L, j X M , k ijk
PK J X K ,i Pi MK J X K ,i M i
(153)
PKL X K ,i i, L X K ,i x j , L i, j MKL X K ,i i , L X K ,i x j ,L i , j
J K J X K ,i i
QK J X K ,i qi
EKL
1
( xi , K xi , L KL ) ELK FKLM xi ,M xi ,KL FLKM
2
E K xi , K L Ei
L L
E K xi , K L Ei
L
BK xi, K L Bi L
BK xi , K L Bi
E KL J xi , K X L, j
L L
ij
L
E KL
J xi , K X L, j L
ij (154)
L
BKL J xi , K X L, j L
ij
L
BKL
J xi , K X L, j L
ij
E K (eff) xi , K i
(eff)
1
GK i xi , K
On a donc :
d d
EKL xi , K x j , L Dij FKLM xk , K x j , L xi , M vi , jk 2 xi , KL x j , M Dij
dt dt
(155)
d d
PK J X K ,i Pi MK J X K ,i i
dt dt
68
d d
PKL X K ,i x j , L ij MKL X K ,i x j , L ij
dt dt
d d dc
0
dt
dt dt
TKL A KL
dEKL
dt
dF
TMLK KLM
dt
d PK L T d MK L T d PKL T d MKL
LE KT BK E KL LBKL J K E K (eff) (156)
dt dt dt dt
R QK GK 0
où AK et
L T L
AKL
T
sont définis par (avec égal à ou ):
L
AKT LAK LAK et
L
AKL
T
LAKL LAKL
(157)
Avec :
dv
R J r , 0 J , J
dV
(158)
1
AKL 2 X K , i X M ,i FPNM TLNP 2 CKM FPNM TLNP
Pour simplifier, les variables internes ne seront pas évoquées dans ces travaux (comme
69
, c , EKL , FKLM , PK , MK , PKL , qui doit être satisfaite pour chaque
facteur. La condition nécessaire et suffisante pour que l’inégalité soit vérifiée est : (voir par
exemple [26] section 10.11):
(160)
c (161)
TKL A KL 0
EKL (162)
TMLK 0
FKLM (163)
L
E KT 0
PK (164)
L
BKT 0
MK (165)
L
E KL
T
0
PKL (166)
L
BKL
T
0
MKL (167)
E (eff)
( R D ) (168)
Q ( R D ) (169)
R
( R D ) (170)
d
J J ( R D ) (171)
dt
70
d
G G ( R D ) (172)
dt
D J K , GK , , aK , bK (174)
¶y
avec bK = -r0 (175)
¶GK
d J K d GK
J K E K (eff) R GK QK aK
dt
bK
dt
0 (176)
J K 0, 0, GK 0, aK 0, bK 0 (177)
variables utiles pour l’étude d’un problème physique. On va de plus développer y en une
série de Taylor allant jusqu’à l’ordre 2. Par exemple, pour un problème électromécanique où
0 KL K L 0 22 IJKL
1 1 1 1 1 1
IJ KL cIJKL EIJ EKL 33cKLMNOP FKLM FNOP
2 2 2 2
1
C 0 hKIJ K EIJ 12 KIJ
1
22 f IJKL EIJ 31 f HIJK FHIJ 32 f HIJKL FHIJ
2
K IJ KL K KL
2 (178)
32
cHIJKL FHIJ EKL p 0
1
K
0 K
p 0
22 1
KL
0 KL KL 0 EKL
71
On obtient alors les équations constitutives suivantes :
TKL A KL 0 cIJKL EIJ 32 cHIJKL FHIJ hKLI PI 22
f IJKL PIJ
EKL
(179)
KL 0 2KL
c 2TKL
0
TMLK 0 33cKLMNOP FNOP 31 f KLMP PP 32 f KLMNO PNO
FKLM
(180)
32 cKLMNO ENO 3 KLM 0 3KLM
c 3TKLM
0
(181)
(182)
(183)
(184)
Si l’une des quantités est grande, il faudra alors développer l’énergie y jusqu’à un ordre
supérieur.
Rappelons tout d’abord les équations d’équilibre (130) dans trouvées précédemment :
tijT, j fi fi em vi (185)
72
avec tijT tij tij
tij ij ijk ,k L
E[i j L B[i j L
[i k
j ,k L
[i k
j ,k
tij ij L
E [i j L B [i j L
[i k
j ,k L
[i k
j ,k
Ce travail généralise les travaux précédents. Par exemple, l’article [19] de Daher et Maugin ne
prend pas en compte les gradients de polarisation, d’aimantation et de déformation. Toutefois,
même avec ces termes, les équations d’équilibre (185) peuvent être réduites et avoir la même
forme que dans [19] :
tij ij L
E[i j L B[i j
(186)
tij ij L
E [i j L B [i j
De même, on obtient les mêmes équations d’équilibre que Collet [15] en négligeant la semi-
conduction et la polarisation. On a alors la forme réduite de tij et de tij :
On trouve aussi une généralisation des équations constitutives obtenus par Majdoub
et al. dans [27] pour des matériaux diélectriques. Dans leur article, ces auteurs avaient utilisés
le principe d’Hamilton introduit par Toupin & Mindlin (voir par exemple [1]) afin de prendre
Ces équations constitutives correspondent au cas particulier suivant des équations (179) et
(180):
73
TKL TKL A KL 0 cIJKL EIJ 32 cHIJKL FHIJ hKLI PI
EKL
(190)
22 f IJKL PIJ KL 0 2KL
c 2TKL
0
TMLK 0 33cKLMNOP FNOP 31 f KLMP PP 32 f KLMNO PNO
FKLM
(191)
32 cKLMNO ENO 3 KLM 0 3KLM
c 3TKLM
0
Commençons par négliger les effets de chaleur en rapport avec les effets de
aK bK GK QK 0 (171)
J E
K K
(eff)
R 0 (172)
Dans la description des processus irréversibles, on donne une relation entre les flux
E P (eff) RPQ
J Q , R A (173)
où RPQ et A dépendent des variables : R c , , EKL , FKLM , PK , MK , PKL , MKL .
Prenant en compte (173), l’inégalité des dissipations (172) prend la forme suivante :
74
J RP PQ J Q A
2
0 (174)
Pour satisfaire l’inégalité précédente, les restrictions sont imposées sur le signe de A et sur
les composantes du tenseur RPQ .
Afin de mettre en évidence les effets apparaissant dans le volume, on décompose RPQ en ses
RPQ K PQ BPQ (175)
où on a posé :
K PQ R( PQ ) KQP
BPQ RPQ BQP
(176)
E P (eff) K PQ
J Q BPQ
J Q (177)
On peut facilement remarquer que le second terme de l’équation (177) ne contribue pas à la
production d’entropie. Ce terme constitue donc un exemple de ce que l’on nomme « effet
caché ». Pour le comprendre, on va considérer le cas particulier suivant :
BPQ PQM BM c , , EKL , FKLM , PK , MK , PKL , MKL (178)
E (eff) K J J B (179)
où l’on peut reconnaitre l’effet Hall qui est bien un effet caché (ne contribuant pas à la
création d’entropie) alors que les effets tels que l’élastorésistance, la magnétorésistance, etc…
sont encore inclus dans l’équation (179).
75
8 Conclusion
Nous avons maintenant une théorie unifiée et valide pour des grandes
électromagnétiques dans le volume avec comme données les forces de volume et de surface
ayant une origine électromagnétique, imposées par l’extérieur. Le champ des vitesses
virtuelles a été construit et le principe des puissances virtuelles a été appliqué afin de trouver
production de l’entropie nous a donné l’inégalité de Clausius-Duhem, très utile pour notre
modélisation. Cette inégalité est la base sur laquelle nous avons pu construire les équations
pour des matériaux semi-conducteurs qui prennent en compte des effets usuels tel que la
piézoélectricité ou la pyroélectricité. Mais nous avons aussi pu décrire de façon très générale
(y compris pour des grandes déformations) des phénomènes plus exotiques comme la
nous avons fait tous ces développements en unités SI plutôt qu’en unités Lorentz-Heaviside
souvent utilisées dans les travaux précédents [1][20][22][29] afin d’être plus proche de toutes
les disciplines.
76
References
[1] R. D. Mindlin, « Polarization gradient in elastic dielectrics », Int. J. Solids Struct., vol. 4,
no 6, p. 637–642, 1968.
[2] A. K. Tagantsev, « Theory of flexoelectric effect in crystals », Zhurnal Eksp. Teor. Fiz.,
vol. 88, no 6, p. 2108–22, 1985.
[3] A. K. Tagantsev, « Piezoelectricity and flexoelectricity in crystalline dielectrics », Phys.
Rev. B, vol. 34, no 8, p. 5883‑5889, oct. 1986.
[4] K. B. Tolpygo, « Physical properties of the salt lattice constructed from deforming ions »,
J Exp Theor PhysUSSR, vol. 20, p. 497, 1950.
[5] K. B. Tolpygo, « Physical properties of a rock salt lattice made up of deformable ions »,
Ukr J Phys, vol. 53, p. 93–102, 2008.
[6] S. M. Kogan, « Piezoelectric effect during inhomogeneous deformation and acoustic
scattering of carriers in crystals », Sov. Phys.-Solid State, vol. 5, no 10, p. 2069–2070,
1964.
[7] R. B. Meyer, « Piezoelectric Effects in Liquid Crystals », Phys. Rev. Lett., vol. 22, no 18,
p. 918‑921, mai 1969.
[8] P. Zubko, G. Catalan, et A. K. Tagantsev, « Flexoelectric effect in solids », Annu. Rev.
Mater. Res., vol. 43, p. 387–421, 2013.
[9] P. V. Yudin et A. K. Tagantsev, « Fundamentals of flexoelectricity in solids »,
Nanotechnology, vol. 24, no 43, p. 432001, nov. 2013.
[10] P. Germain, « La méthode des puissances virtuelles en mécanique des milieux
continus. Première partie : Théorie du second gradient », J Mécanique, vol. 12, no 2, p.
236–274, Juin 1973.
[11] P. Germain, « The Method of Virtual Power in Continuum Mechanics. Part 2:
Microstructure », SIAM J. Appl. Math., vol. 25, no 3, p. 556‑575, nov. 1973.
[12] B. Collet et G. A. Maugin, « Sur l’électrodynamique des milieux continus avec
interactions », C. r. À Académie Sci., vol. t. 279, no Série B, p. 379‑382, nov. 1974.
[13] B. Collet et G. A. Maugin, « Thermodynamique des milieux continus
électromagnétiques avec interactions », C. r. À Académie Sci., vol. t. 279, no Série B, p.
439‑442, nov. 1974.
[14] G. A. Maugin, « A continuum theory of deformable ferrimagnetic bodies. I. Field
equations », J. Math. Phys., vol. 17, no 9, p. 1727‑1738, sept. 1976.
[15] B. Collet, « Higher order surface couplings in elastic ferromagnets », Int. J. Eng. Sci.,
vol. 16, no 6, p. 349‑364, 1978.
[16] G. Maugin, « The method of virtual power in continuum mechanics: application to
coupled fields », Acta Mech., vol. 35, no 1, p. 1–70, 1980.
77
[17] G. A. Maugin et J. Pouget, « Electroacoustic equations for one‐domain ferroelectric
bodies », J. Acoust. Soc. Am., vol. 68, no 2, p. 575‑587, août 1980.
[18] N. Daher et G. A. Maugin, « Virtual power and thermodynamics for electromagnetic
continua with interfaces », J. Math. Phys., vol. 27, no 12, p. 3022‑3035, déc. 1986.
[19] N. Daher et G. A. Maugin, « Deformable semiconductors with interfaces: Basic
continuum equations », Int. J. Eng. Sci., vol. 25, no 9, p. 1093‑1129, 1987.
[20] G. A. Maugin, Continuum mechanics of electromagnetic solids. North-Holland, 1988.
[21] G. A. Maugin, « The principle of virtual power: from eliminating metaphysical forces
to providing an efficient modelling tool », Contin. Mech. Thermodyn., vol. 25, no 2‑4, p.
127‑146, sept. 2011.
[22] A. C. Eringen et G. A. Maugin, Electrodynamics of Continua I: Foundations and Solid
Media. 1990.
[23] A. Askar, P. C. Y. Lee, et A. S. Cakwak, « Lattice-Dynamics Approach to the Theory
of Elastic Dielectrics with Polarization Gradient », Phys. Rev. B, vol. 1, no 8, oct. 1969.
[24] B. Collet, « Sur une théorie des premier et second gradients des milieux continus
électromagnétiques », Phd thesis, Pierre et Maire Curie, 1976.
[25] P. S. Kiréev et S. Medvedev, La physique des semiconducteurs, Moscou : Éd. Mir,
cop. 1975. 1975.
[26] A. C. Eringen, Mechanics of Continua, 2nd edition. Huntington, N.Y.: Krieger Pub
Co, 1980.
[27] M. S. Majdoub, P. Sharma, et T. Cagin, « Enhanced size-dependent piezoelectricity
and elasticity in nanostructures due to the flexoelectric effect », Phys. Rev. B, vol. 77, no
12, p. 125424, 2008.
[28] L. E. Cross, « Flexoelectric effects: Charge separation in insulating solids subjected to
elastic strain gradients », J. Mater. Sci., vol. 41, no 1, p. 53‑63, janv. 2006.
[29] G. A. Maugin et A. C. Eringen, « On the equations of the electrodynamics of
deformable bodies of finite extent », J. Mécanique, vol. 16, p. 101–147, 1977.
78
Chapitre 3
Détermination de la
configuration d’équilibre d’un
ensemble d’atomes polarisables
soumis à un champ électrique
extérieur
1 Introduction
9x) utilisé pour calculer les propriétés électromécaniques de structures carbonées comme le
graphène, les nanotubes de carbones, les fullerènes… L’objectif de ce code est de simuler la
uniforme, que nous imposons, afin de mettre en évidence la flexoélectricité inverse. Dans
cette étude numérique, nous ne prenons toutefois en compte que la déformation, le gradient de
avait également été pris en compte au chapitre 2. Nous travaillons sur des nanotubes de
carbone car, d’une part, ces matériaux ne sont pas piézoélectriques et d’autre part, la
Après avoir présenté brièvement les nanotubes de carbone et leur géométrie, nous
nulle). Nous explicitons ensuite les différentes formes d’énergie potentielle que nous utilisons
pour trouver l’équilibre du système. Un organigramme simplifié du code est présenté à la fin
de ce chapitre.
Les nanotubes de carbone identifiés et caractérisés par Ijima en 1991 [1] ont été et
sont toujours l’objet de beaucoup d’études pour leurs caractéristiques mécanique, électrique
et thermique exceptionnelles [2]. Nous nous contentons toutefois ici de rappeler uniquement
les quelques notions sur leur géométrie qui sont utiles par la suite.
80
On peut représenter un nanotube comme l’enroulement d’une feuille de graphène [3]
autour d’un certain axe. Soit 𝐶⃗ le vecteur circonférence, identifié par deux entiers n et m, tels
que 𝐶⃗ = 𝑛𝑎⃗1 + 𝑚𝑎⃗2 , où (𝑎⃗1 , 𝑎⃗2 ) est la base d’une maille élémentaire de la feuille plane de
Selon l’orientation de l’axe selon lequel on enroule le graphène par rapport aux
liaisons carbone-carbone, on peut distinguer trois types de nanotubes (cf. fig 16a) : armchair
(l’axe du tube est perpendiculaire à un tiers des liaisons carbone-carbone qui sont donc
parallèles à la circonférence), zigzag (l’axe du nanotube est parallèle à des liaisons carbone-
carbone) et chiral (aucune liaison parallèle ou perpendiculaire à l’axe du tube, celui-ci étant
alors non superposable à son image dans un miroir contrairement aux nanotubes des deux
premiers types). On peut démontrer (pour des valeurs de 𝑛 + 𝑚 suffisamment grandes pour
que l’hybridation sp2 des atomes de carbone ne soit pas trop perturbée) que lorsque 𝑛 + 𝑚 est
81
fig 17. Dimensions caractéristiques dans le modèle utilisé pour le graphène
son rayon :
√3𝑚
𝑠𝑖𝑛(𝜃) = et ‖𝐶⃗‖ = 𝐶 = 2𝜋𝑅 = 𝑎√𝑛2 + 𝑚2 + 𝑛𝑚 (1)
2√𝑛2 +𝑚2 +𝑛𝑚
avec 𝑎 = 2,46 Å. On en déduit également que la surface moyenne occupée par un atome de
Pourvu que les vitesses restent petites devant la vitesse de la lumière et que le détail
de la structure des atomes puisse être pris en compte de façon non quantique, la simulation de
la dynamique d’un système à l’échelle atomique, sous l’action de forces extérieures, peut
s’effectuer en utilisant la mécanique newtonienne avec des formes d’énergie potentielle plus
approximations peuvent se révéler judicieuses pour des systèmes contenant un grand nombre
d’atomes car, d’une part, étudier un système comportant des dizaines d’atomes de façon
complètement quantique prendrait déjà énormément de temps et, d’autre part, les masses des
noyaux étant largement plus grandes que celles des électrons, on peut étudier le mouvement
des noyaux qui se fait sur des échelles de temps beaucoup plus grandes que celui des
(traduit en anglais par S.M Blinder)) en ne considérant que les effets moyens des électrons de
manière plus ou moins empirique. En conséquence, cette approximation nous permet de
82
connaitre les énergies d’interactions entre tous les atomes à partir de la position de leur noyau
seulement. L’énergie potentielle d’interaction du système entier ne dépend donc que des
En fait, on peut même démontrer (cf. [6]) que 𝑈 peut s’exprimer par rapport aux seules
distances entre les particules, grâce au fait que l’énergie interne du système doit être
indépendante du repère utilisé pour projeter les vecteurs (critère d’objectivité, cf. chapitre 2) :
̂ ({𝑟 𝛼𝛽 }
𝑈=𝑈 ) avec 𝑟 𝛼𝛽 = ‖𝑟⃗𝛽 − 𝑟⃗ 𝛼 ‖ (4)
𝛼,𝛽=1,…,𝑁
Afin de simuler le mouvement des atomes, chaque position de noyau (donc d’atome)
𝜕 2 𝑟⃗ 𝛼 𝛼
𝐹⃗ 𝛼
∀𝛼, = 𝑎⃗ = 𝛼 (5)
𝜕𝑡 2 𝑚
où 𝐹⃗ 𝛼 , 𝑟⃗ 𝛼 , 𝑎⃗𝛼 , 𝑚𝛼 représentent respectivement la force exercée par les autres atomes sur
Pour connaitre l’équilibre statique du système, on doit trouver les positions telles
que :
∀𝛼, 𝐹⃗ 𝛼 = ⃗⃗
0 (6)
faut donc minimiser l’énergie potentielle 𝑈 en tant que fonction de 3𝑁 variables scalaires (si
on travaille en 3 dimensions) {𝑟⃗ 𝛼 }𝛼=1,…,𝑁 où 𝑁 est le nombre d’atomes.
83
Pour cela, on va utiliser l’algorithme suivant :
déplacement dans cette direction. Elles peuvent être calculées à partir de différentes méthodes
84
telles que la méthode du gradient conjugué ou celle de la plus grande pente [7], ou bien
La dynamique moléculaire (ou MD pour Molecular Dynamics) est une technique qui
permet de simuler l’évolution temporelle d’un système, avec certaines contraintes sur des
constant, potentiel chimique ou nombre de constituant constant, …). Les atomes sont alors
temps est intégrée par un algorithme appelé velocity-Verlet, modifié pour tenir compte des
contraintes macroscopiques imposées (du type température et/ou pression imposée(s)).
Lorsque le système est arrivé à un état d’équilibre dynamique, les grandeurs moyennées sur
tous les atomes (macroscopiques) fluctuent autour de valeurs moyennes (temporelles) qui en
supposant remplis les critères d’applicabilité d’une hypothèse appelée hypothèse ergodique
L’énergie d’interaction des 𝑁 atomes de carbone de notre nanotube sera supposée ici
ne dépendre que de la position de ses atomes (et pas, par exemple, de l’état de leurs électrons
de valence). L’énergie potentielle est donc une fonctionnelle des positions des noyaux des
atomes composant le système qui doit pouvoir modéliser approximativement les forces
d’interactions non seulement entre noyaux, mais aussi entre électrons et noyaux et entre les
électrons. Cette fonctionnelle doit respecter le mieux possible certains critères [10]:
85
1. Flexibilité : la forme fonctionnelle de l’énergie potentielle doit être suffisamment
flexible pour qu’elle puisse s’adapter à des systèmes aussi différents que possible.
2. Précision : l’énergie doit pouvoir se corréler précisément avec les données existantes
être utilisés pour beaucoup de systèmes différents et pas seulement ceux qui ont servi
4. Calcul rapide : le calcul de l’énergie et des forces associées doit se faire suffisamment
rapidement pour que même pour le système le plus grand envisagé, les calculs répétés
précision statistique suffisante. On peut ainsi calculer des milliards de fois cette
complexité du phénomène d’hybridation des orbitales atomiques et des forces à longue portée
associées aux attractions entre liaisons 𝜋. En 1985, Abell [11] utilisa une méthode développée
comportement universel des liaisons covalentes proposé par Rose et al. [13] peut s’écrire sous
2
la forme d’un potentiel anharmonique de type Morse 𝑉𝑀𝑜𝑟𝑠𝑒 (𝑑) = 𝑉0 ([𝑒 −𝛼(𝑑−𝑑0 ) − 1] − 1)
aux propriétés élastiques (raideur de liaison 2𝑉0 𝛼 2 ). Abell proposa alors une expression ne
comportant que des potentiels à 2 corps, donc qui ne dépendent plus que de leurs
interdistances (distance entre les atomes) et non des orientations relatives des liaisons entre
86
atomes (ce qui n’empêche pas la force sur un atome donné de dépendre de ces orientations
𝑈 𝛼𝛽 = 𝑓𝑐 (𝑟 𝛼𝛽 )[𝑉 𝑅 (𝑟 𝛼𝛽 ) − 𝑏 𝛼𝛽 𝑉 𝐴 (𝑟 𝛼𝛽 )] (7)
covalente (Bond Order) entre les atomes 𝛼 et 𝛽. L’expression de ce dernier terme suppose que
les ordres de liaisons covalentes entre atomes sont dépendants de la coordinence des atomes
en question :
𝑏 𝛼𝛽 ~ 𝑍 −𝛿 (8)
une énergie de type « 2 corps ». 𝑓𝑐 (𝑟 𝛼𝛽 ) est une fonction de coupure (variant entre 1 pour les
plus proches voisins à la distance d’équilibre et 0 pour les atomes éloignés, de façon plus ou
moins abrupte) permettant de ne considérer que les atomes qui sont assez proches pour
pouvoir effectivement établir une liaison covalente. Conjuguée avec 𝑏 𝛼𝛽 , cette fonction
permet donc de décrire des liaisons qui pourront changer d’ordre, se rompre ou s’établir au fur
et à mesure que les interdistances évoluent pendant la simulation. Ceci permet donc de rendre
En 1987, à partir des travaux d’Abell [11], Tersoff construisit un potentiel décrivant
des énergies de liaisons covalentes simple, double ou triple pouvant modéliser des liaisons
entre atomes de silicium [14], [15], de germanium [16] ou de carbone [17] ainsi que des
interactions avec des défauts [18]. Pour cela, il paramétra l’ordre 𝑏 𝛼𝛽 des liaisons covalentes,
non seulement en fonction des positions des atomes, mais aussi des angles entre les liaisons.
L’énergie 𝑈 est alors composée d’une énergie à « 2 corps » pour la répulsion et d’une énergie
87
En 1990, un potentiel appelé REBO (Reactive Empirical Bond-Order) sera paramétré
par Brenner [19] à partir du potentiel de Tersoff, en rajoutant les interactions avec les atomes
d’hydrogène ainsi que d’autres termes pour modéliser les structures ayant des radicaux libres
et des liaisons conjuguées. Cette nouvelle forme d’énergie potentielle a permis la simulation
systèmes présentant des défauts ou des dislocations. Il a ainsi pu être utilisé de nombreuses
fois avec succès pour étudier les propriétés mécaniques et thermiques de structures carbonées
comme les nanotubes de carbone [20], les fullerènes [21], du carbone amorphe [17] ou du
graphite [22]. Cependant, le potentiel REBO n’est pas approprié pour décrire n’importe quel
système hydrocarboné. Les interactions dues aux dipôles induits et permanents ainsi que les
termes de répulsion intermoléculaire ne sont pas pris en compte, alors qu’ils sont non
négligeables pour beaucoup de systèmes hydrocarbonés, dont les liquides et les systèmes en
couches comme le graphite, les oignons de fullerène ou les nanotubes de carbone à plusieurs
couches. De même pour les matériaux présentant des liaisons covalentes fortes comme le
diamant où les forces à plus longue portée deviennent importantes si on étudie les interfaces
du système. De plus, il est impossible de modéliser des collisions entre atomes à cause de
l’absence de termes divergents de répulsion lorsque les distances interatomiques sont proches
de zéro [23].
À la fin des années 90, pour résoudre les problèmes liés au potentiel REBO,
différents potentiels font leur apparition. On peut par exemple nommer le potentiel utilisé par
Che et al. [24] pour la modélisation des liaisons non covalentes ou de Pettifor et al. [25], [26].
De leur côté, Brenner et al. ont paramétré une forme un peu plus complète d’énergie sous le
dépendant de l’angle dièdre pour les liaisons doubles. Une reparamétrisation des fonctions
splines permit également d’améliorer l’accord des résultats pour des liaisons avec des
radicaux libres ou pour des liaisons conjuguées avec les résultats de calculs ab-initio. Ceci a
88
diamants [28], [29]. Brenner et al. expliquèrent en détails la paramétrisation de cette forme
fig 19. Définition de l’angle dièdre utilisé lorsque BC est une liaison double
Pendant ce temps, Stuart et al. ont élaboré une nouvelle fonctionnelle d’énergie
[31] afin de prendre en compte les effets de torsion de liaisons ainsi que les interactions de
type van der Waals modélisées par un potentiel de type Lennard-Jones et des fonctions de
lissage.
fig 20. Illustration des termes de torsion et de liaisons faibles rajoutés dans le potentiel
AIREBO.
Cette énergie a été paramétrée pour inclure les interactions entre les atomes
d’hydrogène, de carbone ou d’oxygène. Cette énergie a pu être utilisée dans les nanotubes de
carbone (voir par exemple [32]). L’énergie AIREBO s’écrit finalement:
89
𝑁 𝑁
1
𝑈 𝐴𝐼𝑅𝐸𝐵𝑂 = ∑ ∑ 𝑓𝑐 (𝑟 𝛼𝛽 ) 𝑉 𝑅 (𝑟 𝛼𝛽 ) − 𝑏 𝛼𝛽 𝑉 𝐴 (𝑟 𝛼𝛽 ) + 𝑉𝐿𝐽 (𝑟 𝛼𝛽 )
2
𝛼=1 𝛽=1
𝛽≠𝛼 [
𝑁 𝑁
(9)
𝑡𝑜𝑟
+ ∑ ∑ 𝑉𝛼𝛽𝛾𝛿
𝛾=1 𝛿=1
𝛾≠𝛼,𝛽 𝛿≠𝛼,𝛽,𝛾 ]
où les deux premiers termes représentent l’énergie REBO2, le troisième terme représente le
l’environnement des atomes. Il ne s’agit pas d’une constante mais d’une fonction compliquée
où 𝑏 𝛼𝛽−𝜎𝜋 dépend de la coordination (nombre d'atomes voisins les plus proches) des atomes
𝛼 et 𝛽, 𝑏 𝛽𝛼−𝑅𝐶 est spécifique des structures ayant des radicaux libres et/ou des liaisons
conjuguées et 𝑏 𝛽𝛼−𝐷𝐻 dépend des angles dièdres si la liaison 𝛽𝛼 est une liaison double. Les
force que va subir chaque atome. Il est aussi possible de définir une force interatomique 𝑓⃗𝛼𝛽
qui pourra être interprétée comme étant la force qu’exerce l’atome 𝛽 sur l’atome 𝛼. Dans ce
qui suit, nous allons préciser la façon de calculer la force interatomique à l’aide de dérivées de
La force interne sur un atome 𝛼 (force due aux interactions avec tous les autres atomes du
pour la forme fonctionnelle de l’énergie potentielle permet de montrer que l’on peut écrire
̂ 𝑖𝑛𝑡 ({𝑟 𝛼𝛽 }
𝑈 𝑖𝑛𝑡 = 𝑈 ) (13)
𝛼,𝛽=1,…,𝑁
On peut donc utiliser la règle de dérivation des fonctions composées pour pouvoir calculer
𝑓⃗𝑖𝑛𝑡,𝛼 en dérivant l’énergie par rapport aux distances entre chaque atome de la manière
suivante :
𝜕𝑈 𝑖𝑛𝑡 (𝑟⃗ 1 , 𝑟⃗ 2 , … , 𝑟⃗ 𝑁 ) ̂ 𝑖𝑛𝑡 ({𝑟 𝛽𝛾 }) 𝜕𝑟 𝛽𝛾
𝜕𝑈
𝑓⃗𝑖𝑛𝑡,𝛼 = − = − ∑
𝜕𝑟⃗ 𝛼 𝜕𝑟𝛽𝛾 𝜕𝑟⃗ 𝛼
𝛽𝛾
𝛽<𝛾
1 ̂ 𝑖𝑛𝑡 ({𝑟 𝛽𝛾 }) 𝜕𝑟 𝛽𝛾
𝜕𝑈 (14)
=− ∑
2 𝜕𝑟𝛽𝛾 𝜕𝑟⃗ 𝛼
𝛽𝛾
𝛽≠𝛾
Or
2
𝛾 𝛽
𝛽𝛾
𝜕𝑟
3
𝜕√∑3𝑖=1(𝑟𝑖 − 𝑟𝑖 ) 3 𝛼
1 2(𝑟𝑗 − 𝑟𝑗 )
𝛾𝛽 𝛼
1 2(𝑟𝑗 − 𝑟𝑗 )
=∑ 𝑒⃗𝑗 = ∑ [ 𝛿𝛼𝛾 − 𝛿𝛼𝛽 ] 𝑒⃗𝑗
𝜕𝑟⃗ 𝛼 𝜕𝑟𝑗𝛼 2 𝑟𝛽𝛼 2 𝑟 𝛼𝛾
𝑗=1 𝑗=1 (15)
𝛽𝛼 𝛼𝛾
𝑟⃗ 𝑟⃗
= 𝛽𝛼 𝛿𝛼𝛾 − 𝛼𝛾 𝛿𝛼𝛽
𝑟 𝑟
On peut donc calculer 𝑓⃗𝑖𝑛𝑡,𝛼 en utilisant les forces interatomiques 𝑓⃗𝛼𝛽 que les autres atomes
Remarquons que le fait que l’énergie potentielle ne dépende que des interdistances
n’empêche pas que la force s’exerçant sur un atome donné dépende des orientations de ses
différentes liaisons avec les autres atomes, à cause des termes en 𝑟⃗ 𝛼𝛽 dans la sommation.
91
4.4 Prise en compte des conditions aux limites mécaniques
Lors d’une simulation mécanique ou d’un essai mécanique, on doit toujours préciser
les conditions dans laquelle l’expérience a été faite. Par exemple, lorsque l’on veut faire un
essai de traction sur un matériau, on impose une force sur une partie du matériau et on bloque
le déplacement d’une autre partie. En simulation atomistique, on doit pouvoir imposer une
où 𝑓⃗𝛼 est une force donnée que l’on applique sur l’atome 𝛼, 𝑟⃗ 𝛼 donne la position de l’atome
𝛼 dans la configuration actuelle (déformée) et 𝑅⃗⃗ 𝛼 donne la position de ce même atome 𝛼 dans
En dérivant cette énergie potentielle par rapport à la position actuelle d’un atome 𝛽
À l’équilibre 𝑓⃗𝛽 est nulle. Nous avons donc un équilibre entre les forces internes au système
et les forces externes imposées. Dans la simulation, on trouvera donc les forces à l’équilibre
en ajoutant les forces 𝑓⃗𝛽 , imposées sur les atomes de C, aux forces intérieures calculées avec
contraint un atome à rester immobile (pour simuler un encastrement par exemple), cela
revient à considérer que la force exercée par cet encastrement 𝑓⃗𝛼 compense exactement, à
tout instant, la résultante des forces intérieures sur cet atome (− 𝜕𝑈 𝑖𝑛𝑡 ⁄𝜕𝑟⃗𝛽 ).
92
5 Effet d’un champ électrique extérieur sur une population
d’entités polarisables
Nous allons maintenant présenter le modèle adopté pour décrire la polarisation d’un
système semi-conducteur sous l’effet d’un champ extérieur, dans l’approximation où cette
polarisation est due à des dipôles induits localisés sur chaque atome, décrit par sa
polarisabilité dipolaire linéaire, ainsi qu’à des dipôles permanents créés par la courbure de la
feuille de graphène pour former un nanotube de carbone (hybridation des liaisons qui n’est
plus totalement sp2) [33], [34]. Nous serons ensuite capables de simuler l’équilibre d’un
système sous l’influence d’un champ électrique extérieur par minimisation des énergies
Sous l’effet d’un champ électrique extérieur, un ensemble d’atomes se polarise : pour
chaque atome, les forces 𝐹⃗ = 𝑞𝐸⃗⃗ sur le noyau, d’une part, et sur les électrons, d’autre part,
sont de sens opposés puisque leurs charges le sont. Les barycentres des charges positives et
négatives ne sont donc plus confondus et un dipôle (induit) apparait sur chaque atome. Dans
l’approximation de la réponse linéaire, ce dipôle 𝑝⃗𝛽 est lié linéairement au champ électrique
s’exerçant sur l’atome 𝛽 considéré, par un tenseur symétrique d’ordre 2 appelé tenseur de
polarisabilité 𝛼̿ 𝛽 :
Notons que dans l’équation ci-dessus, le champ électrique, noté 𝐸⃗⃗𝑙𝑜𝑐 (𝑟⃗𝛽 ), s’exerçant sur
l’atome 𝛽 est classiquement appelé champ local, mais qu’il ne s’agit pas du même champ
local que celui introduit au chapitre 2. Il s’agit ici du champ électrique local définit dans un
milieu où les charges électriques sont discrétisées. Ce champ est responsable de la force qui
s’exercerait sur une charge fictive de test positionnée en 𝑟⃗𝛽 , en l’absence de l’atome 𝛽. D’un
93
autre côté, le champ électrique macroscopique introduit au chapitre 2, est défini dans un
changement d’hybridation des orbitales des électrons de valence des atomes de carbone. Cet
effet a été étudié par Yakobson et ses collaborateurs. Ainsi en 2002, Dumitrika et al. [33] ont
montré par des calculs ab-initio que ce changement d’hybridation crée un dipôle permanent
sur chaque atome de carbone, fonction de l’angle Θ𝜎𝜋 entre les liaisons 𝜎 et 𝜋. Ceci fut
feuillets par Kvashnin et al. [34]. D’après ces auteurs, le dipôle permanant créé par la
où 𝜃𝑝𝛼 = Θ𝜎𝜋 − 𝜋⁄2 correspond à l’angle pyramidal entre l’atome 𝛼 et le plan 𝜋 formé par
ses 3 plus proches voisins (cf. Définition de l’angle pyramidal entre l’atome 𝛼 et les 3 atomes
formant le plan 𝜋.fig 21), 𝑛⃗⃗𝛼 correspond au vecteur unitaire perpendiculaire au plan 𝜋 et
dirigé vers l’atome 𝛼 et 𝑓𝜃 = 2,34 D/rad est une constante calculée par lissage linéaire des
valeurs calculées de 𝜇⃗ 𝛼 en fonction de 𝜃𝑝𝛼 pour des nanotubes d’indices variés, par Kvashnin
et al. [34]. Compte tenu du fait qu’un debye est environ égal à 3,33564 × 10−30 C·m ou
fig 21. Définition de l’angle pyramidal entre l’atome 𝛼 et les 3 atomes formant le plan 𝜋.
94
Pour un atome donné d’un nanotube de carbone non déformé, on peut évaluer 𝜃𝑝𝛼 facilement à
partir de l’équation (cf. [34]): 𝑠𝑖𝑛(𝜃𝑝𝛼 ) ≈ d⁄4𝑅 où 𝑑 = 1,42 Å est la distance entre l’atome 𝛼
et ses trois plus proches voisins et 𝑅 est le rayon du nanotube. Notons que pour un fullerène, il
où la courbure totale ∁ est égale à 1⁄𝑅 pour un nanotube et 2⁄𝑅 pour un fullerène. Ceci
permet à Kvashnin et al. de définir un coefficient de flexoélectricité local pour des structures
symétriques par :
𝑑
∀𝛼 = 1, … , 𝑁, 𝜇⃗ 𝛼 = 𝑓𝜃 × ∁ 𝑛⃗⃗𝛼 = 𝐹 × ∁ 𝑛⃗⃗𝛼 (22)
4
On obtient 𝐹 = 0,831 D·Å pour 0,82 dans l’article de Dumitruca et al., soit encore = 0,173 e
Å2/atome en bon accord avec les résultats compris entre 0,157 e Å2/atome et 0,187 e
Å2/atome, selon la structure, obtenus par Kalinin et Meunier [35], avec des calculs DFT. Pour
passer du moment dipolaire à la polarisation et donc diviser par le volume d’un atome dans la
feuille de graphène, ce qui pose le problème de son épaisseur ! Si, pour calculer un ordre de
grandeur, on prend comme épaisseur 0,8 Å comme Kalinin et Meunier qui citent Huang et
al.[36], il vient :
Pour des nanotubes fléchis, nous avons pour notre part calculé 𝜃𝑝𝛼 à l’aide de la
formule suivante4 :
3
1 𝑟⃗ 𝛼𝑖 − 𝑟⃗ 𝛼 𝜋
∀𝛼 = 1, … , 𝑁, 𝜃𝑝𝛼 ⁄ −1 𝛼
= 𝛩𝜎𝜋 − 𝜋 2 ≈ ∑ |𝑐𝑜𝑠 (𝑛⃗⃗ ∙ ( 𝛼 )) − | (23)
3 ‖𝑟⃗ 𝑖 − 𝑟⃗ 𝛼 ‖ 2
𝑖=1
4
(voir fig 21), ‖𝑟⃗ 𝛼𝑖 − 𝑟⃗ 𝛼 ‖ donne la distance entre deux atomes de la pyramide, et n
⃗⃗α ∙ r⃗ αi − r⃗ α donne la hauteur
𝛼
de la pyramide on obtient donc facilement l’expression de l’angle 𝜃𝑝 en faisant une moyenne sur les 3 atomes
constituant la base de la pyramide
95
Notons que pour un nanotube non-déformé, tous ces moments dipolaires se compensent à
Dans cette section, nous allons introduire une quantité tensorielle appelée
propagateur dipolaire du champ électrique en calculant le champ électrique 𝐸⃗⃗ en tout point
d’un matériau modélisé par 𝑁 dipôles électriques ponctuels {𝑝⃗𝛼 }, soumis à un champ
électrique 𝐸⃗⃗0 (𝑟⃗) dérivant d’un potentiel électrique 𝑉0 (𝑟⃗) créé par des sources de champ
matériau. Pour cela, nous allons résoudre l’équation de Poisson (∆𝑉 = − 𝜌⁄𝜀0 ) en
charges5 𝜌𝑚𝑎𝑡 (𝑟⃗) = ∑𝑁 ⃗⃗𝛿(𝑟⃗ − 𝑟⃗ 𝛼 )), associée aux dipôles du matériau, plus une
⃗𝛼 ∙ (−∇
𝛼=1 𝑝
densité de charges 𝜌𝑒𝑥𝑡 (𝑟⃗) comprenant toutes les sources de 𝑉0. En notant 𝑉(𝑟⃗) = 𝑉0 (𝑟⃗) +
𝑉1 (𝑟⃗), (avec ∆𝑉0 (𝑟⃗) = − 𝜌𝑒𝑥𝑡 (𝑟⃗)⁄𝜀0), on cherche donc un potentiel 𝑉1 solution de :
∑𝑁 ⃗𝛼 ∙ 𝛻⃗⃗ 𝛿(𝑟⃗ − 𝑟⃗ 𝛼 )
𝛼=1 𝑝
∆𝑉1 (𝑟⃗) = (24)
𝜀0
avec des conditions aux limites correspondant au vide pour 𝑉1 (condition de radiation de
Sommerfeld), ce qui correspond à des conditions aux limites pour 𝑉 telles que 𝑉 soit continu
et égal à 𝑉0 à la limite d’une distance infinie du matériau (qui est par hypothèse, de taille
5
Pour un dipôle centré en 𝑟⃗ 𝛼 :
96
∑𝑁 ⃗𝛼 ∙ 𝛻⃗⃗′𝛿(𝑟⃗′ − 𝑟⃗ 𝛼 )
𝛼=1 𝑝
𝑉1 (𝑟⃗) = ∭ 𝐺0 (𝑟⃗, 𝑟⃗′) 𝑑𝑟⃗′
𝜀0
𝑁 𝑝 ⃗𝛼 (25)
=∑ ∭ 𝐺0 (𝑟⃗, 𝑟⃗′)𝛻⃗⃗′𝛿(𝑟⃗′ − 𝑟⃗ 𝛼 )𝑑𝑟⃗′
𝜀
𝛼=1 0
où 𝐺0 est la distribution de Green associée à l’opérateur laplacien (i.e. ∆𝑟⃗ 𝐺0 (𝑟⃗, 𝑟⃗′) = 𝛿(𝑟⃗, 𝑟⃗′))
et à des conditions aux limites nulles à l’infini en tant que fonction de 𝑟⃗, soit :
−1
𝐺0 (𝑟⃗, 𝑟⃗′) = (26)
4𝜋‖𝑟⃗ − 𝑟⃗′‖
Une intégration par partie de l’intégrale (25) donne alors (compte tenu de fait que le terme
intégré est nul puisque 𝛿(𝑟⃗′ − 𝑟⃗ 𝛼 ) est toujours nul quand 𝑟⃗′ → ∞, vu que les 𝑟⃗ 𝛼 sont à
distance finie) :
𝑁 𝑝⃗𝛼 𝑁 1 𝑟⃗ − 𝑟⃗ 𝛼
𝑉1 (𝑟⃗) = − ∑ ⃗⃗ (𝑟
𝛻 ′ 𝐺0 ⃗, 𝑟⃗′) = ∑ ( )( ) ∙ 𝑝⃗𝛼 (27)
𝛼=1 𝜀0 𝛼=1 4𝜋𝜀0 ‖𝑟
⃗ − 𝑟⃗ 𝛼 ‖ 3
la dernière égalité provient du fait que ⃗⃗′(1⁄‖𝑟⃗ − 𝑟⃗′‖) = (𝑟⃗ − 𝑟⃗′)⁄‖𝑟⃗ − 𝑟⃗′‖3 =
∇
⃗⃗(1⁄‖𝑟⃗ − 𝑟⃗′‖)
−∇
Remarquons que l’on retrouve l’expression classique du potentiel créé par un dipôle en ne
Finalement :
𝑁
(1)
𝑉(𝑟⃗) = 𝑉0 (𝑟⃗) + ∑ 𝑇̅0 (𝑟⃗, 𝑟⃗ 𝛼 )𝑝⃗𝛼 (28)
𝛼=1
où
(1) 1 𝑟⃗ − 𝑟⃗′
𝑇̅0 (𝑟⃗, 𝑟⃗ 𝛼 ) = (29)
4𝜋𝜀0 ‖𝑟⃗ − 𝑟⃗′‖3
(1)
La quantité 𝑇̅0 (𝑟⃗, 𝑟⃗ 𝛼 ) est alors appelée propagateur électrique d’ordre 1 du vide. Il permet
de calculer le potentiel créé par un dipôle (ou le champ créé par une charge) en tout point de
l’espace sauf en 𝑟⃗ = 𝑟⃗ où une autre méthode est utilisé (voir le paragraphe 5.2.2).
97
ce qui permet de calculer le potentiel et le champ en tout point dès lors que les valeurs des
La quantité
(2) (1)
𝑇̿0 (𝑟⃗, 𝑟⃗′) = (−𝛻⃗⃗) ⊗ 𝑇̅0 (𝑟⃗, 𝑟⃗′)
(31)
= (1⁄4𝜋𝜀0 ) [3(𝑟⃗ − 𝑟⃗′) ⊗ (𝑟⃗ − 𝑟⃗′) − ‖𝑟⃗ − 𝑟⃗′‖2 1̿]⁄‖𝑟⃗ − 𝑟⃗′‖5
est alors appelée propagateur électrique d’ordre 2 (ou, ici, dipolaire) du vide.
Plus généralement, on peut définir les propagateurs électriques d’ordre 𝑛 + 𝑚 du
vide, par (voir le livre de A. Stone, p. 16, [37] avec toutefois une différence de signe pour les
𝑇 (𝑚,𝑛) (𝑟⃗, 𝑟⃗′) pour des situations plus compliquées que le vide, grâce à deux remplacements.
aux limites du vide est remplacée par une autre fonction de Green de cette même équation
avec les conditions aux limites appropriées pour prendre en compte les surfaces et interfaces
du matériau considéré comme un milieu continu abritant des dipôles ponctuels différents de
ses entités propres. D’autre part, (−1⁄𝜀0 ) serait remplacé par (−(𝜀̿𝑟 )−1 ⁄𝜀0 ) , avec 𝜀̿𝑟 le
tenseur de permittivité relative de ce matériau. Notons pour finir que, par exemple, lorsque le
matériau est borné par une surface plane, le calcul du potentiel électrique et de ses gradients
matériau sera tel que 𝑇 (𝑚,𝑛) (𝑟⃗, 𝑟⃗′) ≠ 𝑇 (𝑛,𝑚) (𝑟⃗, 𝑟⃗′) contrairement à ce qu’il se passe dans le
vide et que l’on ne peut donc plus définir alors de propagateur 𝑇 (𝑚+𝑛) (𝑟⃗, 𝑟⃗′) qui peut
98
Régularisation des propagateurs
localisés aux positions des noyaux des atomes pose toutefois quelques problèmes. Ainsi, si
l’on veut calculer le champ à proximité immédiate des atomes ou à l’intérieur des nuages
électroniques, l’atome ne peut plus être considéré comme un dipôle ponctuel puisque 𝑟⃗ 𝛼 ≈
𝑟⃗𝛽 . Nous ne sommes donc plus dans des conditions pour lesquelles la distance entre le point
où l’on calcule le champ et le centre du dipôle est très grande devant les dimensions
problème, nous faisons l’hypothèse que les multipôles équivalents à la distribution de charges
réelles ne sont pas ponctuels mais proviennent d’une distribution continue de charges à
Dans cette approximation, le champ électrique créé par un dipôle qui n’est plus
Ceci est en fait équivalent à utiliser un propagateur régularisé (au sens où, comme on le verra
Plutôt que d’essayer de calculer cette intégrale directement, il peut être plus
(2) (0)
𝑇̿0 (𝑟⃗, 𝑟⃗′) = (−𝛻⃗⃗) ⊗ (−𝛻⃗⃗)𝑇0 (𝑟⃗, 𝑟⃗′) (36)
99
avec 𝑇0(0) (𝑟⃗, 𝑟⃗′) = 1⁄(4𝜋𝜀0 ‖𝑟⃗ − 𝑟⃗′‖). En supposant que l’on peut intervertir intégrale sur 𝑟⃗′ et
dérivation par rapport à 𝑟⃗ (ce qui n’est pas forcément évident à cause de la divergence de la
Or ∫ 𝑇00 (⃗𝑟⃗, ⃗𝑟⃗′)𝑔𝛼 (⃗𝑟⃗′) 𝑑𝑟⃗⃗′ est calculable relativement aisément avec un logiciel de calcul
( )
Comme lim erf(𝑥) = 1, on note que l’on retrouve bien le propagateur standard du vide
𝑥→∞
lorsque 𝑎 → 0. De plus, on vérifie que le propagateur converge en 𝑟⃗ = 𝑟⃗′:
(0) 1 2 1
𝑙𝑖𝑚 𝑇̂0 (𝑟⃗, 𝑟⃗ 𝛼 ) = × × 𝛼 (39)
𝛼
‖𝑟⃗−𝑟⃗ ‖→0 4𝜋𝜀0 √𝜋 𝑅
ce qui illustre les propriétés de régularisation de la convolution d’un propagateur avec une
gaussienne.
(2)
̂
̿ ) (𝑟⃗, 𝑟⃗′), il vient :
Si on en revient au calcul de (𝑇 0
(2)
̂
̿ ) (𝑟⃗, 𝑟⃗ 𝛼 ) = (−𝛻⃗⃗) ⊗ (−𝛻⃗⃗)𝑇̂0(0) (𝑟⃗, 𝑟⃗ 𝛼 )
(𝑇 0
1 1
= ×
4𝜋𝜀0 ‖𝑟⃗ − 𝑟⃗𝛼 ‖5
‖𝑟⃗ − 𝑟⃗ 𝛼 ‖
[3(𝑟⃗ − 𝑟⃗𝛼 ) ⊗ (𝑟⃗ − 𝑟⃗ 𝛼 ) − ‖𝑟⃗ − 𝑟⃗𝛼 ‖2 1̿]𝑒𝑟𝑓 ( )
𝑅𝛼 (40)
‖𝑟⃗ − 𝑟⃗ 𝛼‖ ‖𝑟⃗ − 𝑟⃗ 𝛼‖ 3 ‖𝑟⃗ − 𝑟⃗ 𝛼‖ 2
𝛼 𝛼
6 4
× − [(𝑟⃗ − 𝑟⃗ ) ⊗ (𝑟⃗ − 𝑟⃗ ) ( (
𝑅𝛼
)+ (
𝑅𝛼
) )] 𝑒𝑥𝑝 (− (
𝑅𝛼
) )
√𝜋 √𝜋
2
2 ‖𝑟⃗ − 𝑟⃗𝛼 ‖ ‖𝑟⃗ − 𝑟⃗𝛼 ‖
+1̿ [‖𝑟⃗ − 𝑟⃗ 𝛼 ‖2 ( ( ))] 𝑒𝑥𝑝 (− ( ) )
[ √𝜋 𝑅𝛼 𝑅𝛼 ]
provenant d’une distribution de charges à symétrie sphérique avec une distribution radiale
gaussienne, en un point de l’espace où l’on ne mettrait qu’une charge fictive de test.
100
Cherchons maintenant quelle serait l’énergie d’interaction entre deux distributions de charges
à symétrie sphérique et à distribution radiale gaussienne. Pour cela, on démontre qu’il nous
faut faire la convolution du propagateur standard du vide non plus avec une mais deux
fonctions gaussiennes, et que l’on peut simplifier cela en définissant un propagateur régularisé
(0)
similaire au 𝑇̂0 de la section précédente :
𝛼𝛽 (0) (0)
𝑈𝑞−𝑞 = ∫ ∫ 𝑇0 (𝑟⃗ ′ , 𝑟⃗")𝑞 𝛼 𝑔𝛼 (𝑟⃗′)𝑞 𝛽 𝑔𝛽 (𝑟⃗") 𝑑𝑟⃗ ′ 𝑑𝑟⃗" ≡ 𝑞 𝛼 𝑇̃0 (𝑟⃗ 𝛼 , 𝑟⃗𝛽 )𝑞 𝛽 (41)
et que :
avec :
2 2
𝛼𝛽 𝛽
𝑎𝑞−𝑞 = √(𝑅𝑞𝛼 ) + (𝑅𝑞 ) (43)
comme il se doit pour que l’énergie d’interaction entre 𝛼 et 𝛽 soit la même que celle entre 𝛽 et
𝛼.
Pour le propagateur d’ordre 2, utilisé pour calculer l’énergie d’interaction entre deux
(2)
̃
̿)
dipôles « gaussiens », on définit un propagateur régularisé (𝑇 de même forme que
0
(2) 2
̂
̿ ) (𝑟⃗ 𝛼 , 𝑟⃗𝛽 ), sauf que le paramètre de régularisation est maintenant 𝑎𝑝−𝑝
(𝑇 𝛼𝛽 𝛼 2 𝛽
= √(𝑅𝑝 ) + (𝑅𝑝 ) et non
0
plus 𝑅 𝛼 .
Expression générale
Considérons une structure carbonée semi-conductrice (donc sans charge libre) dont
la réponse à un champ extérieur est calculée dans l’approximation dipolaire électrique, avec
des dipôles permanents 𝜇⃗ et/ou induits 𝑝⃗. L’expression de l’énergie due à l’interaction avec le
101
champ extérieur et à la polarisation de la structure est donnée dans Thole [38] et Mayer [39].
Elle s’exprime sous forme d’une fonction des variables, 𝑟⃗ 𝛿 , 𝑝⃗𝛿 , 𝜇⃗𝛿 et 𝐸⃗⃗0 (𝑟⃗ 𝛿 ), où 𝛿 prend toutes
avec
(2)
̃
𝐸⃗⃗𝑝⃗𝛽+𝜇⃗⃗𝛽 (𝑟⃗𝛽𝛼 ) = [(𝑇
̿ ) (𝑟⃗𝛽𝛼 )] (𝑝⃗𝛽 + 𝜇⃗ 𝛽 )
0 (46)
le champ créé en 𝑟⃗ 𝛼 par les dipôles induits et permanents situés en 𝑟⃗𝛽 et 𝑟⃗𝛽𝛼 = 𝑟⃗ 𝛼 − 𝑟⃗𝛽 .
Le premier terme de l’équation (45) correspond à l’énergie nécessaire pour polariser les
atomes (création des 𝑝⃗𝛼 ou « self-energy »). Le deuxième terme correspond à l’énergie
d’interaction entre les dipôles et le troisième à l’interaction des dipôles avec le champ
extérieur.
𝑖𝑛𝑑
(Les abréviations « dip » et « ind » dans l’énergie 𝑈𝑑𝑖𝑝 correspondent respectivement à
Notons que les dipôles permanents 𝜇⃗ 𝛼 considérés ici sont dus à la courbure de la
feuille de graphène par rapport à la situation où le nanotube est déformé. Puisque l’influence
des (𝜇⃗ 𝛼 ) est supposée être déjà prise en compte de façon phénoménologique dans le
potentiel AIREBO, on déduit les interactions entre les dipôles permanents dans le dernier
terme.
À partir de cette énergie, on peut trouver les valeurs auto-cohérentes des dipôles
𝑖𝑛𝑑
induits puisqu’elles doivent minimiser 𝑈𝑑𝑖𝑝 à l’équilibre électrique. En utilisant les équations
102
(45) et (46), ainsi que l’invariance par échange des indices d’atomes (et de coordonnées) du
(2)
̃
̿ ) , il vient :
propagateur régularisé (𝑇 0
𝑖𝑛𝑑
𝜕𝑈𝑑𝑖𝑝
∀𝛿 = 1, … , 𝑁, ⃗0⃗ = ( )
𝜕𝑝⃗𝛿 é𝑞
𝑁
(2)
−1 ̃
̿ ) (𝑟⃗𝛽𝛿 )] (𝑝⃗𝛽 + 𝜇⃗ 𝛽 ) − 𝐸⃗⃗0 (𝑟⃗ 𝛿 )
= (𝛼̿ 𝛿 ) 𝑝⃗𝛿 − ∑ [(𝑇
(47)
0
𝛽=1
𝛽≠𝛿
Les valeurs de dipôles induits à l’équilibre électrique doivent donc vérifier les équations
(auto-cohérentes) suivantes :
𝑁
(2)
∀𝛿 = 1, … , 𝑁, 𝛿 −1 ̃
(𝛼̿ ) (𝑝⃗𝛿 )𝑒𝑞 = 𝐸⃗⃗0 (𝑟⃗ 𝛿 ) + ∑ [(𝑇
̿ ) (𝑟⃗𝛽𝛿 )] ((𝑝⃗𝛽 ) + 𝜇⃗ 𝛽 )
0 𝑒𝑞 (48)
𝛽=1
𝛽≠𝛿
−1
Le champ local 𝐸⃗⃗loc (𝑟⃗ 𝛿 ) = (𝛼̿ 𝛿 ) (𝑝⃗𝛿 )𝑒𝑞 est donc la somme du champ imposé par
l’extérieur du système et des champs créés par les autres dipôles du système.
Il est en fait possible de calculer tous les dipôles induits à l’équilibre ensemble, en ré
exprimant les équations ci-dessus sous forme d’un système de 𝑁 équations linéaires entre
vecteurs que l’on met ensuite sous forme matricielle. En effet, les équations (45) peuvent se
Ceci peut donc se mettre sous la forme 𝐴̂𝑋̂ = 𝑌̂, avec 𝐴̂ une matrice symétrique à 3𝑁 lignes et
103
(2) (2)
(𝛼̿ 1 )−1 ̃
̿ ) (𝑟⃗12 )
−(𝑇 ⋯ ̃
̿ ) (𝑟⃗1𝑁 )
−(𝑇
0 0
(2) (2)
̃
̿ ) (𝑟⃗ 21 )
−(𝑇 (𝛼̿ 2 )−1 ⋯ ̃
̿ ) (𝑟⃗ 2𝑁 )
−(𝑇
𝐴̂ = 0 0 (50)
⋮ ⋮ ⋱ ⋮
(2) (2)
̃
̿ 𝑁1 ̃
̿ 𝑁2 ⋯ (𝛼̿ 𝑁 )−1
(−(𝑇)0 (𝑟⃗ ) −(𝑇)0 (𝑟⃗ ) )
𝑇
𝑋̂ = ((𝑝⃗1 )𝑒𝑞 (𝑝⃗2 )𝑒𝑞 ⋯ (𝑝⃗𝑁 )𝑒𝑞 ) (51)
et
𝑁
(2)
𝐸⃗⃗0 ⃗
(𝑟 1 ) ̃
̿ ) (𝑟⃗𝛽1 )] 𝜇⃗ 𝛽
+ ∑ [(𝑇 0
𝛽=1
𝛽≠1
𝑁
(2)
̃
𝐸⃗⃗0 (𝑟⃗ 2 ) + ∑ [(𝑇
̿ ) (𝑟⃗𝛽2 )] 𝜇⃗ 𝛽
𝑌̂ = 𝛽=1
0
(52)
𝛽≠2
⋮
𝑁
(2)
̃
𝐸⃗⃗0 (𝑟⃗ 𝑁 ) + ∑ [(𝑇
̿ ) (𝑟⃗𝛽𝑁 )] 𝜇⃗ 𝛽
0
𝛽=1
( 𝛽≠𝑁 )
Notons également que l’on peut aussi exprimer (46) sous la forme :
𝑁
(2)
∀𝛿 = 1, … , 𝑁,
−1 ̃
̿ ) (𝑟⃗𝛽𝛿 )(1 − 𝛿𝛿𝛽 )]] ((𝑝⃗𝛽 ) + 𝜇⃗ 𝛽 )
∑ [(𝛼̿ 𝛿 ) 𝛿𝛿𝛽 − [(𝑇 0 𝑒𝑞
𝛽=1 (53)
𝛿 −1 𝛿
= 𝐸⃗⃗0 (𝑟⃗ 𝛿 ) + (𝛼̿ ) 𝜇⃗
et donc que l’on peut calculer les ((𝑝⃗𝛽 )𝑒𝑞 + 𝜇⃗ 𝛽 ), puis les (𝑝⃗𝛽 )𝑒𝑞 , avec une équation
matricielle du même type (𝐴̂𝑋̂′ = 𝑌̂′) avec la même matrice, mais un second membre plus
simple à calculer.
104
𝑁
𝑖𝑛𝑑 1
𝑈𝑑𝑖𝑝,𝑒𝑞 = ∑(𝑝⃗𝛼 )𝑒𝑞 (𝛼̿ 𝛼 )−1 (𝑝⃗𝛼 )𝑒𝑞
2
𝛼=1
𝑁 𝑁
1 ̃ (2)
̿ ) (𝑟⃗𝛽𝛼 )] ((𝑝⃗𝛽 ) + 𝜇⃗ 𝛽 )
− ∑((𝑝⃗𝛼 )𝑒𝑞 + 𝜇⃗ 𝛼 ) ∙ ∑ [(𝑇
2 0 𝑒𝑞
𝛼=1 𝛽=1
𝛽≠𝛼
𝑁 𝑁 𝑁
1
− ∑((𝑝⃗ 𝑒𝑞 + 𝜇⃗ ) ∙ 𝐸⃗⃗0 (𝑟⃗ 𝛼 ) + ∑ ∑ 𝜇⃗ 𝛼 ∙ 𝐸⃗⃗𝜇⃗⃗𝛽 (𝑟⃗𝛽𝛼 )
𝛼) 𝛼
2
𝛼=1 𝛼=1 𝛽=1
𝛽≠𝛼
𝑁 𝑁
1 1
𝑖𝑛𝑑
𝑈𝑑𝑖𝑝,𝑒𝑞 = ∑(𝑝⃗𝛼 )𝑒𝑞 (𝛼̿ 𝛼 )−1 (𝑝⃗𝛼 )𝑒𝑞 − ∑((𝑝⃗𝛼 )𝑒𝑞 + 𝜇⃗ 𝛼 ) ∙ ((𝛼̿ 𝛼 )−1 (𝑝⃗𝛼 )𝑒𝑞 − 𝐸⃗⃗0 (𝑟⃗ 𝛼 ))
2 2
𝛼=1 𝛼=1
𝑁 𝑁 𝑁
1
− ∑((𝑝⃗𝛼 )𝑒𝑞 + 𝜇⃗ 𝛼 ) ∙ 𝐸⃗⃗0 (𝑟⃗ 𝛼 ) + ∑ ∑ 𝜇⃗ 𝛼 ∙ 𝐸⃗⃗𝜇⃗⃗𝛽 (𝑟⃗𝛽𝛼 )
2
𝛼=1 𝛼=1 𝛽=1
𝛽≠𝛼
𝑁 𝑁 𝑁 𝑁
1 1 1
𝑖𝑛𝑑
𝑈𝑑𝑖𝑝,𝑒𝑞 = − ∑ 𝜇⃗ 𝛼 (𝛼̿ 𝛼 )−1 (𝑝⃗𝛼 )𝑒𝑞 − ∑((𝑝⃗𝛼 )𝑒𝑞 + 𝜇⃗ 𝛼 ) ∙ 𝐸⃗⃗0 (𝑟⃗ 𝛼 ) + ∑ ∑ 𝜇⃗ 𝛼 ∙ 𝐸⃗⃗𝜇⃗⃗𝛽 (𝑟⃗𝛽𝛼 )
2 2 2 (54)
𝛼=1 𝛼=1 𝛼=1 𝛽=1
𝛽≠𝛼
Rappel de notation :
- 𝑝𝑘𝛼 , 𝜇𝑘𝛼 et 𝑟𝑘𝛼 sont les 𝑘 𝑖è𝑚𝑒 composantes, dans la base considérée, respectivement du
atome 𝛽, en se souvenant que les (𝑝⃗𝛼 )𝑒𝑞 et les 𝜇⃗ 𝛼 sont des fonctions des positions des atomes
et que l’équilibre électrique n’est pas forcément l’équilibre mécanique et donc que l’on peut
dériver aussi bien l’expression (54) que l’expression (45) pour trouver l’expression des forces
105
électriques à l’équilibre électrique (ce que nous avons vérifié explicitement après quelques
𝛽 𝛽
𝑖𝑛𝑑
𝐹𝑖 = −𝛻𝑖 𝑈𝑑𝑖𝑝,𝑒𝑞
𝑁 𝑁
1 𝛽 1 𝛽
= ∑ 𝛻𝑖 (𝜇⃗ 𝛼 )(𝛼̿ 𝛼 )−1 (𝑝⃗𝛼 )𝑒𝑞 + ∑ 𝜇⃗ 𝛼 (𝛼̿ 𝛼 )−1 𝛻𝑖 (𝑝⃗𝛼 )𝑒𝑞
2 2
𝛼=1 𝛼=1
𝑁
1 𝛽
+ ∑ 𝛻𝑖 ((𝑝⃗𝛼 )𝑒𝑞 + 𝜇⃗ 𝛼 ) ∙ 𝐸⃗⃗0 (𝑟⃗ 𝛼 )
2
𝛼=1
(55)
𝑁 𝑁 𝑁
1 𝛽 1 𝛽
+ ∑((𝑝⃗𝛼 )𝑒𝑞 + 𝜇⃗ 𝛼 ) ∙ 𝛻𝑖 𝐸⃗⃗0 (𝑟⃗ 𝛼 ) + 𝛻𝑖 ∑ ∑ 𝜇⃗ 𝛼 ∙ 𝐸⃗⃗𝜇⃗⃗𝛽 (𝑟⃗𝛽𝛼 )
2 2
𝛼=1 𝛼=1 𝛽=1
𝛽≠𝛼
( )
Compte tenu du fait que nous avons calculé les dipôles permanents en fonction des
écarts de position des atomes à la section 5.1.2, il est possible de calculer numériquement ou
𝛽 𝛽
analytiquement les 𝛻𝑖 (𝜇⃗ 𝛼 ). Il ne nous reste donc « plus qu’à » calculer les 𝛻𝑖 (𝑝⃗𝛼 )𝑒𝑞 ou, de
𝛽
façon équivalente, les 𝛻𝑖 [(𝑝⃗𝛿 )𝑒𝑞 + 𝜇⃗ 𝛿 ]. Pour cela, nous nous servons de la forme matricielle
𝛼𝛽
𝐴̂𝑋̂′ = 𝑌̂′ que nous réécrivons 𝑋̂′ = (𝐴̂−1 )𝑌̂′ et nous définissons (𝐴̂−1 ) les 𝑁 2 sous-
et donc :
𝛽
∀𝛽, 𝛿 = 1, … , 𝑁, 𝛻𝑖 [(𝑝⃗𝛿 )𝑒𝑞 + 𝜇⃗ 𝛿 ]
𝑁
𝛽 𝛿𝛼
= ∑ [𝛻𝑖 (𝐴̂−1 ) ] (𝐸⃗⃗0 (𝑟⃗ 𝛼 ) + (𝛼̿ 𝛼 )−1 𝜇⃗ 𝛼 )
𝛽=1 (57)
𝑁
𝛿𝛼 𝛽
+ ∑(𝐴̂−1 ) 𝛻𝑖 (𝐸⃗⃗0 (𝑟⃗ 𝛼 ) + (𝛼̿ 𝛼 )−1 𝜇⃗𝛼 )
𝛼=1
𝛽 𝛿𝛼
Notons que les quantités 𝛻𝑖 (𝐴̂−1 ) peuvent être calculées analytiquement grâce au fait que
𝛽 𝛽 𝛽 𝛽 𝛽
0̂ = 𝛻𝑖 (𝐴̂−1 𝐴̂) = 𝛻𝑖 (𝐴̂−1 )𝐴̂ + (𝐴̂−1 )𝛻𝑖 (𝐴̂) et donc que 𝛻𝑖 (𝐴̂−1 ) = (𝐴̂−1 )𝛻𝑖 (𝐴̂)(𝐴̂−1 ),
106
𝛽 ̃ (3) dont
puisque les éléments de 𝛻𝑖 (𝐴̂) peuvent se calculer avec des propagateurs de type (𝑇) 0
6 Déroulement de la programmation
fig 22. Algorithme de recherche de l’équilibre statique du système soumis à l’action d’un
champ électrique extérieur 𝐸⃗⃗0
107
108
7 Conclusion
Nous avons ensuite présenté un modèle atomistique (codé en langage Fortran 9x)
va nous permettre de calculer les dipôles induits résultant d’un champ extérieur que l’on va
appliquer sur la structure. C’est par l’interaction entre ces dipôles induits que notre structure
ne sera plus en équilibre mécanique. La déformation de la structure est alors calculé à partir
de l’interaction entres les dipôles via le modèle dipôles-dipôles, et par le potentiel d’AIREBO
qui va traduire toutes les liaisons entre les atomes (covalentes, Lennard-Jones..). Enfin, nous
avons utilisé les résultats de Dumitrica, Kalinin et Meunier et de Kvashnin et al. afin de
pouvoir ajouter les dipôles permanents créés par la courbure du matériau dans les interactions
les dipôles permanents et l’angle de pyramidation avec 𝑓𝜃 = 0,487 e Å/rad. C’est en rajoutant
ces dipôles créés par une hybridation des liaisons de chaque atome que nous serons
109
Références
[1] S. Iijima, « Helical microtubules of graphitic carbon », Nature, vol. 354, no 6348, p.
56‑58, nov. 1991.
[2] M. S. Dresselhaus, Éd., Carbon nanotubes: synthesis, structure, properties, and
applications. Berlin: Springer, 2001.
[3] K. S. Novoselov et al., « Two-dimensional gas of massless Dirac fermions in
graphene », ArXiv Prepr. Cond-Mat0509330, 2005.
[4] Z. Wang, « Propriétés électro-mécaniques des nanotubes de carbone », Université de
Franche-Comté, 2008.
[5] M. Born et R. Oppenheimer, « Zur Quantentheorie der Molekeln », Ann. Phys., vol.
389, no 20, p. 457‑484, 1927.
[6] N. C. Admal et E. B. Tadmor, « A Unified Interpretation of Stress in Molecular
Systems », J. Elast., vol. 100, no 1‑2, p. 63‑143, juin 2010.
[7] W. H. Press, Éd., Numerical recipes in Fortran 90: the art of parallel scientific
computing, 2. ed., Reprinted with corr. Cambridge: Cambridge Univ. Press, 1999.
[8] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, et P. Gumbsch, « Structural relaxation
made simple », Phys. Rev. Lett., vol. 97, no 17, p. 170201, 2006.
[9] D. C. Rapaport, The Art of Molecular Dynamics Simulation, 2e éd. Cambridge:
Cambridge University Press, 2004.
[10] D. W. Brenner, « The art and science of an analytic potential », Phys. Status Solidib,
vol. 217, no 1, p. 23–40, 2000.
[11] G. C. Abell, « Empirical chemical pseudopotential theory of molecular and metallic
bonding », Phys. Rev. B, vol. 31, no 10, p. 6184, 1985.
[12] P. W. Anderson, « Self-consistent pseudopotentials and ultralocalized functions for
energy bands », Phys. Rev. Lett., vol. 21, no 1, p. 13, 1968.
[13] J. H. Rose, J. R. Smith, et J. Ferrante, « Universal features of bonding in metals »,
Phys. Rev. B, vol. 28, no 4, p. 1835, 1983.
[14] J. Tersoff, « New empirical model for the structural properties of silicon », Phys. Rev.
Lett., vol. 56, no 6, p. 632, 1986.
110
[15] J. Tersoff, « New empirical approach for the structure and energy of covalent
systems », Phys. Rev. B, vol. 37, no 12, p. 6991, 1988.
[16] J. Tersoff, « Modeling solid-state chemistry: Interatomic potentials for
multicomponent systems », Phys. Rev. B, vol. 39, no 8, p. 5566, 1989.
[17] J. Tersoff, « Empirical Interatomic Potential for Carbon, with Applications to
Amorphous Carbon », Phys. Rev. Lett., vol. 61, no 25, p. 2879‑2882, déc. 1988.
[18] J. Tersoff, « Carbon defects and defect reactions in silicon », Phys. Rev. Lett., vol. 64,
no 15, p. 1757, 1990.
[19] D. W. Brenner, « Empirical potential for hydrocarbons for use in simulating the
chemical vapor deposition of diamond films », Phys. Rev. B, vol. 42, no 15, p. 9458, 1990.
[20] R. S. Ruoff, D. Qian, et W. K. Liu, « Mechanical properties of carbon nanotubes:
theoretical predictions and experimental measurements », Comptes Rendus Phys., vol. 4, no
9, p. 993‑1008, nov. 2003.
[21] D. W. Brenner, J. A. Harrison, C. T. White, et R. J. Colton, « Molecular dynamics
simulations of the nanometer-scale mechanical properties of compressed
Buckminsterfullerene », Thin Solid Films, vol. 206, no 1‑2, p. 220–223, 1991.
[22] D. H. Robertson, D. W. Brenner, et J. W. Mintmire, « Energetics of nanoscale
graphitic tubules », Phys. Rev. B, vol. 45, no 21, p. 12592, 1992.
[23] P. de Sainte Claire, K. Song, W. L. Hase, et D. W. Brenner, « Comparison of ab initio
and empirical potentials for H-atom association with diamond surfaces », J. Phys. Chem., vol.
100, no 5, p. 1761–1766, 1996.
[24] J. Che, T. Çağın, et W. A. Goddard Iii, « Generalized extended empirical bond-order
dependent force fields including nonbond interactions », Theor. Chem. Acc., vol. 102, no 1‑6,
p. 346–354, 1999.
[25] D. G. Pettifor et I. I. Oleinik, « Analytic bond-order potentials beyond Tersoff-
Brenner. I. Theory », Phys. Rev. B, vol. 59, no 13, p. 8487, 1999.
[26] I. I. Oleinik et D. G. Pettifor, « Analytic bond-order potentials beyond Tersoff-
Brenner. II. Application to the hydrocarbons », Phys. Rev. B, vol. 59, no 13, p. 8500, 1999.
[27] D. Srivastava, D. W. Brenner, J. D. Schall, K. D. Ausman, M. Yu, et R. S. Ruoff,
« Predictions of Enhanced Chemical Reactivity at Regions of Local Conformational Strain on
Carbon Nanotubes: Kinky Chemistry », J. Phys. Chem. B, vol. 103, no 21, p. 4330‑4337, mai
1999.
111
[28] A. Garg, J. Han, et S. B. Sinnott, « Interactions of carbon-nanotubule proximal probe
tips with diamond and graphene », Phys. Rev. Lett., vol. 81, no 11, p. 2260, 1998.
[29] O. A. Shenderova, D. W. Brenner, A. Omeltchenko, X. Su, et L. H. Yang, « Atomistic
modeling of the fracture of polycrystalline diamond », Phys. Rev. B, vol. 61, no 6, p. 3877,
2000.
[30] D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni, et S. B. Sinnott,
« A second-generation reactive empirical bond order (REBO) potential energy expression for
hydrocarbons », J. Phys. Condens. Matter, vol. 14, no 4, p. 783‑802, févr. 2002.
[31] S. J. Stuart, A. B. Tutein, et J. A. Harrison, « A reactive potential for hydrocarbons
with intermolecular interactions », J. Chem. Phys., vol. 112, no 14, p. 6472‑6486, avr. 2000.
[32] A. Nikitin et al., « Hydrogenation of Single-Walled Carbon Nanotubes », Phys. Rev.
Lett., vol. 95, no 22, nov. 2005.
[33] T. Dumitrică, C. M. Landis, et B. I. Yakobson, « Curvature-induced polarization in
carbon nanoshells », Chem. Phys. Lett., vol. 360, no 1, p. 182–188, 2002.
[34] A. G. Kvashnin, P. B. Sorokin, et B. I. Yakobson, « Flexoelectricity in carbon
nanostructures: Nanotubes, Fullerenes, and Nanocones », J. Phys. Chem. Lett., vol. 6, no 14,
p. 2740–2744, 2015.
[35] S. V. Kalinin et V. Meunier, « Electronic flexoelectricity in low-dimensional
systems », Phys. Rev. B, vol. 77, no 3, janv. 2008.
[36] Y. Huang, J. Wu, et K. C. Hwang, « Thickness of graphene and single-wall carbon
nanotubes », Phys. Rev. B, vol. 74, no 24, déc. 2006.
[37] A. J. Stone, The theory of intermolecular forces, 2. ed. Oxford: Oxford Univ. Press,
2013.
[38] B. T. Thole, « Molecular polarizabilities calculated with a modified dipole
interaction », Chem. Phys., vol. 59, no 3, p. 341‑350, août 1981.
[39] A. Mayer, « A monopole-dipole model to compute the polarization of metallic carbon
nanotubes », Appl. Phys. Lett., vol. 86, no 15, p. 153110, avr. 2005.
112
Chapitre 4
Dans le chapitre 2, nous avons obtenu les équations reliant les variables mécaniques,
manière générale, les variables telles que le tenseur de Cauchy du premier ou du deuxième
ordre ou les tenseurs d’élasticité sont connues expérimentalement en réalisant des tractions,
C’est à partir de la configuration d’équilibre de notre système d’atomes que nous allons
extraire les champs continus tels que l’énergie interne, les contraintes, les déformations ainsi
114
fig 23. Déformation du domaine par la transformation de la configuration initiale 0 à la
configuration finale t
Soit un référentiel RG dans lequel on a choisi un repère constitué d’une origine et de la base
cartésienne orthonormée : (e1 , e2 , e3 ) . Soit , une partie quelconque d’un système
OM 0 X X p e p OM t x x p e p (1)
M t Nt y x
(3)
(Y , t ) ( X , t )
115
Lorsque Y X tend vers 0, on peut développer par le développement de Taylor suivant :
On peut alors introduire deux tenseurs qui dépendent du point M 0 et de t. On peut les écrire
i x 2 i 2 xi
FiJ i et GiJK (6)
X J X J X J X K X J X K
On peut alors réécrire l’expression du développement limité sous forme indicielle, jusqu’à
l’ordre 2:
1
M t Nt i FiJ M 0 N0 J GiJK M 0 N0 J M 0 N0 K (7)
2
ri M ri ri M Oi OM i (8)
116
3 Représentation graphique des gradients de la transformation
Pour une meilleure explication sur le rôle de FiJ et de GiJK dans la déformation, nous
allons représenter graphiquement leur influence dans la déformation 2D d’un maillage carré
où la position est calculée par l’équation (9), avec des tenseurs F et G constants. Sans
Le titre de chaque graphique correspond aux paramètres que l’on va changer par rapport aux
paramètres initiaux. Par exemple, si le titre de la figure est « F12 a G122 b » alors on
117
fig 25. Représentation des déformations du premier ordre (G=0)
Le terme F12 donne la pente du maillage des lignes d’équation x F12 y , tandis que F11
change la taille du maillage selon x . Le graphique en bas à gauche nous montre que lorsque
F12 F21, on obtient clairement une rotation du maillage. Cette propriété peut être expliquée
par la décomposition polaire de F : lorsque F est antisymétrique alors la déformation pure U
vaut l’identité (le système n’est pas déformé) et il ne reste plus que la rotation : F=RU=R
avec U=Id.
118
fig 26. Représentation de déformations du deuxième ordre
fig 27. Représentation des grandeurs caractéristiques d’une déformation due au second
gradient de transformation G
On suppose que le maillage non déformé avait une longueur L selon l’axe x et une longueur h
selon l’axe y (qui sont égales à 2 dans la fig 24). Le tenseur G permet alors des déformations
119
non homogènes dans le maillage. Par exemple, les composantes G211 et G122 donnent la
courbure des lignes de maillage. La composante G112 G121 va changer la pente du maillage
G L
tan 112 (10)
2 2
De plus, à condition que l’angle soit petit G112L et on peut montrer que si les lignes
courbées du maillage sont perpendiculaires aux lignes représentant la section du maillage
(lignes parallèles à l’ordonnée lorsqu’il n’y a pas de déformation), alors nécessairement
G112 G211 .
extérieur. Le potentiel électrique sera donc nul et seul le potentiel AIREBO contribuera au
contenue dans une région B tel que S=A+B. On identifie les positions des atomes de chaque
partie par :
120
R A R A ,
RB R B (11)
Selon les différentes parties, on peut décomposer l’énergie potentielle des interactions
internes U par :
U AB (r A , r B )
f 0 où S (13)
r
U A (r A ) U AB (r A , r B )
f 0
r r
U int (r A ) U ext (r A , r B ) (14)
r r
f int, f ext ,
Dans cette expression, nous avons décomposé la force totale s’exerçant sur un atome de la
région A du système S isolé, en une force interne f int, et une force externe f ext , .
Principe de la méthode
121
Pour étudier l’élasticité d’un cristal, Cauchy fit l’hypothèse qu’à tout instant la
déformation est la même pour chaque atome dans le réseau. La relation entre la position d’un
atome dans l’état déformé et la position de cet atome dans l’état non déformé s’écrit alors :
r FR (15)
où le tenseur F est le même pour tous les atomes.
Plus tard, Born montra que la méthode de Cauchy est trop restrictive car elle ne peut
pas décrire la dynamique à température fixée de la plupart des cristaux. Born étendit la
méthode de Cauchy à une température non nulle en ne l’utilisant pas pour les positions
instantanées mais pour la position moyenne des atomes à l’équilibre. Born considéra aussi le
cas des cristaux avec des motifs non triviaux en rajoutant une déformation interne associée au
mouvement relatif des sous-réseaux [1], [2]. L’étude de Friesecke et Theil sur la validité de la
méthode de Cauchy-Born a montré que la méthode fonctionne tant que la déformation reste
suffisamment homogène à l’échelle de la taille des cristaux [3]. Dans le cas où la déformation
n’est plus homogène, on pourrait utiliser par exemple ce qui s’appelle « the quasicontinuum
method » développée par Tadmor et al. [4]–[6]. Cette méthode permet de coupler la méthode
continue et la méthode atomistique de façon un peu plus élaborée. Lorsque la déformation est
homogène dans un certain sous-ensemble d’atomes, il est plus rapide de ne pas considérer
tous les atomes explicitement et de supposer que tous les atomes du sous ensemble vérifient
l’hypothèse de Cauchy-Born. Cette hypothèse permettra de gagner en temps de calcul pour les
forces et les énergies entre les atomes. Il est donc plus judicieux de faire une interpolation sur
le sous-domaine en utilisant la méthode des éléments finis par exemple. Lorsque l’on se
trouve dans un domaine où l’on ne peut plus considérer que la déformation est uniforme, on
va plutôt prendre en compte tous les atomes dans les calculs et faire une simulation
atomistique avec ces atomes. Cependant, cette méthode d’homogénéisation ne fonctionne plus
pour des matériaux cristallins tels que le graphène ou les nanotubes de carbone. Lu a étudié
l’élasticité des nanotubes sans prendre en compte les effets de courbure du tube [7]. En effet,
le gradient de transformation décrit le comportement tangentiel à la surface (voir fig 28). Pour
pallier à ce problème, différentes méthodes ont été proposées par Friesecke et James qui
122
utilisent la méthode de Cosserat [8] ou par Arroyo et al. qui utilisent une méthode qu’ils
une application exponentielle [9]) ou encore par Sunyk et Steinmann qui supposèrent que le
gradient du gradient de transformation est constant [10]. Park et al. [11] utilisent la méthode
Dans notre étude, nous avons choisi d’utiliser « la méthode de Cauchy-Born étendue
et locale » développée par Sunyk et Steinmann [10] afin de faire le lien entre approches
déformation puisse être constant et suffisamment grand par rapport à l’échelle atomique pour
que les moyennes sur tous les atomes ne changent pas de façon appréciable si on ajoute ou
enlève les atomes d’une maille élémentaire. On suppose que les tenseur F et G sont constants
123
Détermination de F et G par la méthode des moindres carrés
Pour trouver les valeurs de FiJ et GiJK correspondant aux valeurs moyennes des
déformations issues de la simulation atomistique, Il faut minimiser la somme des carrés des
i 2 AX AX T X T AAX T XNX
T
(18)
avec
Le minimum de 2 (𝑋) est alors atteint lorsque X est un vecteur propre correspondant à la plus
petite des valeurs propres de A. Remarquons que N étant une matrice définie positive (puisque
N T AA ), ses valeurs propres sont forcément toutes positives.
sur l’axe x1 , auquel on impose un champ extérieur E0 dans le plan x1 , x2 , à un angle de 45°
par rapport à l’axe x1 du nanotube non déformé et bloquons les atomes du premier anneau.
Par la méthode des moindres carrés on trouve FiJ et GiJK en remarquant que le tenseur F est
assez proche d’une matrice rotation et que les composantes de G ayant les composantes les
plus grandes sont G211 et G112 G121 qui sont les composantes qui traduisent la flexion dans le
124
plan des axes x1 et x2 . Lorsque l’on fait varier la norme du champ électrique extérieur, on
obtient la fig 29 :
fig 29. Variation de la composante G211 (en Å 1 ) par rapport au champ électrique extérieur
(en V / Å ) d’un nanotube de carbone (10,0) d’une longueur de 12,39 nm
La courbure donnée par G211 est donc bien d’autant plus grande que le champ électrique
extérieur augmente. On peut aussi remarquer que les composantes G211 et G112 G121 ont une
relation de proportionnalité comme le montre la fig 30.
125
fig 30. Variation de la composante G211 par rapport à la composante G112 G121 lors de la
flexion d’un nanotube de carbone (10,0) d’une longueur de 12,39 nm
Dans le chapitre 2, nous avons obtenu les équations permettant de faire le lien entre
l’énergie interne et les variables macroscopiques telles que les tenseurs de contraintes ou le
champ local. Dans le chapitre 3, nous avons vu comment nous pouvons calculer les énergies
potentielles U entre les atomes d’un système à partir du potentiel AIREBO et du potentiel
électrique traduisant les interactions dipolaires. C’est à partir de l’énergie de Helmholtz par
126
unité de masse et l’énergie potentielle U que nous ferons le lien entre l’approche
sont constants dans le VER dont le centre est le point M (mais pas forcément les mêmes que
En utilisant l’équation (21), on fait alors le lien entre les variables atomistique et les
variables continues :
( F ) 1 U (r ) 1
U r
fi RJM RJM
3
1
PiJ 0 int k (26)
FiJ V0 FiJ 2 , A k 1 rk FiJ 2V0 , A
( F ) 1 U (r ) 1 U int rk
fi RKM RLM RKM RLM (29)
3
1
P 0
(2)
GiJK V0 GiJK 2 , A k 1 rk GiJK 4V0 , A
iJK
numériques, nous allons exprimer ces tenseurs par rapport à la force interne fi int, . En effet, si
le système contient n atomes, alors il faudra calculer n(n 1) 2 forces si on utilise les
int,
équations (26)-(29) alors que l’utilisation des forces internes fi ne nous demanderait que n
calculs.
1
2V0
f i
int,
R 2V1 f
M
J i
int,
R
M
J
(30)
A 0 A
1
V0
f i
int,
R M
J
A
128
4.4 Calcul du tenseur d’élasticité
On utilise la dérivée partielle pour trouver le tenseur d’élasticité mixte que nous
établissons à partir des équations constitutives que nous avons données durant le chapitre 2 et
ses annexes. On a :
P
D (32)
F
1
P
f i
int,
R
M
J
avec :
carbone en utilisant les équations (28) et (35). On prend dans ce cas un nanotube de carbone
(10,0) d’une longueur initiale (sans déformation) de 123,92 Å, contenant exactement 1200
atomes. L’axe du nanotube est orienté suivant l’axe 1. Dans l’expérience suivante, on va créer
129
avoir créé une traction sur le nanotube en changeant la position initiale R par la position r
E (37)
11
Le module d’Young axial vaut donc environ 0,99 TPa. Ce qui est conforme avec les tests déjà
faits dans la littérature. [12]
Notons qu’en utilisant directement les équations (33) et (35), on peut trouver le
130
1.13 0.04 0.04 0 0 0
0.04 0.45 0.14 0 0 0
0.04 0.14 0.45 0 0 0
0 0 0 0.31 0 0
0 0 0 0 0.5 0
0 0 0 0 0 0.5
On a ainsi obtenu une matrice caractéristique d’une isotropie transverse selon son axe
longitudinal de la forme :
11 c c1122 c1122 0 0 0 11
1111
22 c1122 c2222 c2233 0 0 0 22
33 c c2233 c2222 0 0 0 33
1122
2 23 0 0 0 c2222 c2233 0 0 2 23
2 31 0 0 0 0 2c1212 0 2 31
2 12 0 0 0 0 0 2c1212 212
Cette matrice nous donne donc toutes les constantes matériau usuelles comme les modules
d’Young ou les coefficients de Poisson. On a ainsi obtenu un module d’Young orthoradial de
0,45 TPa et un module longitudinal de 1,13 TPa. Nous obtenons donc un module d’Young
axial légèrement supérieur au module d’Young du test dans lequel on a imposé la
déformation. Ceci provient du fait que lorsque l’on fait une variation de la position pour
calculer les dérivées numériques, on ne réoptimise pas la géométrie de la structure. Cela
pourrait être aussi dû au fait que nous n’avions pas pris en compte les déplacements internes
du matériau. De plus nous n’avons pas pris en compte les effets de surface dans la méthode de
Cauchy-Born [11].
La façon dont on calcule la Hessienne ik lors de l’application de l’équation (33) demande
131
5 Passage de l’échelle microscopique à l’échelle macroscopique
pour la partie induite par le champ électrique extérieur
calcul avec le seul potentiel AIREBO) due à la présence d’un champ électrique extérieur
vaut :
ind
U dip
r , p , , E0 r
2 1
p p p E p r
1 N 1 1 N N
2 1 1
1
: p p
(38)
E0 r .E r
N
1 N N
p
1 2 1 1
avec : 2
E p r E p
r E r , p , T r p (39)
0
C’est à partir de cette énergie que nous allons calculer le champ local du chapitre 2, puis le
tenseur de flexoélectricité, à partir des équations constitutives.
(38) et en supposant que les n premiers atomes appartiennent à ce VER, on a l’énergie totale
ind ,VER
due aux dipôles électriques induits et restreinte au VER, notée U dip qui correspond à :
132
p p p E p r .E r
1 n 1 1 n n 1 n n
ind ,VER
U dip
2 1 2 1 1 2 1 1
(40)
n N n
p E p r p E0 r
1 n 1 1
Exprimons maintenant le fait que chaque atome du VER est à l’équilibre dans une
population donnée de dipôles permanents , soit en utilisant l’ équation (40) :
U dip
ind ,VER
. p
1 n 1 n N
T r p E0 r
1
, eq
T r p T r p
p 2 1 2 1 n 1
. p T r p E0 r 0
N
1 (41)
1
d’équilibre trouvée dans le chapitre 3, on vérifie que l’on obtient bien la même équation si on
dérive l’énergie du système entier au lieu de celle du VER, ce qui montre que les dipôles
induits sont toujours à l’équilibre quel que soit le VER que l’on choisit. Chaque dipôle induit
p du VER est donc une fonction des autres dipôles du système entier
p
, 1,..., N ; et du champ extérieur et dépend donc de l’extérieur. On
construit ensuite l’énergie interne (énergie potentielle des interactions entre constituants du
VER) et l’énergie externe (énergie potentielle des interactions avec l’extérieur) du VER. Dans
l’expression de l’énergie totale du VER, on obtient l’énergie interne en ne retenant que les
interactions entre atomes du VER. L’énergie des actions de l’extérieur du VER sur celui-ci est
alors l’ensemble des termes restants dans l’énergie totale de VER. On partitionne ainsi
l’énergie totale due aux dipôles électriques et donnée par (40) de la manière suivante :
133
p p p E p r .E r
1 n 1 1 n n 1 n n
ind ,VER
U dip
2 1 2 1 1 2 1 1
1 n 1
1
U dind
ip
,VER ,Ext
Les deux derniers termes de l’équation (42) correspondent à l’énergie des actions
mécaniques extérieures au VER sur celui-ci. Ces actions sont les actions mécaniques à
« longue » distance entre les dipôles électriques extérieurs au VER et le champ imposé par
l’extérieur du système global et les dipôles inclus dans le VER, que l’on peut réécrire de la
manière suivante :
p
n N n
, EVER r , r p E p r p E0 r
ind ,VER , Ext
U dip
1 n 1 1
n N
p E0 r E p r (43)
1 n 1
n
p EVER r
1
où EVER correspond au champ électrique créé au point par l’extérieur du VER (ou par le
VER ) c’est à dire dû aux dipôles électriques extérieurs au VER et au champ extérieur :
N
EVER r E p r E0 r (44)
n1
Ce champ électrique peut être interprété comme étant le champ électrique extérieur appliqué
aux atomes du VER de la même façon que représente le champ imposé par l’extérieur du
système total.
134
L’énergie interne dipolaire du VER
L’énergie interne due aux dipôles électriques induits du VER correspond aux deux
premiers termes de l’équation (42) que l’on peut réécrire de la manière suivante :
1 n
p p p E p r
1 1 n n
r r , p ,
ind ,VER , Int
U dip
r 2 1 2 1 1
.E r
1 n n
2 1 1
n
1 n n
T r : p p
1 1
: p
p
2 1 2 1 1
(45)
.E r
1 n n
2 1 1
n
1 1 1 n n
2 1
: p
p
2 1 1
T r : p p
p p
n n
T r :
1 1 2
Cette énergie interne due aux dipôles est fonction des variables internes objectives suivantes :
La polarisation
se calculer par :
n
1
p
P (46)
VVER 1
135
Le champ local
Dans cette partie, en plus de la notation vectorielle, on adopte aussi une notation
tensorielle avec une sommation implicite sur les indices romains. Dans la configuration
1 U dip
ind ,VER ,int
L T
VVER Pk (47)
k
F,G
Compte tenu de la séparation que nous avons faite entre l’énergie des liaisons
(induite) du champ électrique extérieur au système, seule l’énergie interne électrique dépend
L T
1 U dip
ind ,VER ,int
p , , r
(48)
VVER Pk
k
F,G
L T
Nous allons maintenant chercher à exprimer k à partir des variables microscopiques.
tenseur G était constant dans notre VER. De même, afin de calculer le champ local, nous
allons avoir besoin de faire une hypothèse d’homogénéisation sur la configuration électrique
de notre système. Nous ferons l’hypothèse que les variations de dipôles par rapport à la
configuration d’équilibre en l’absence de champ extérieur, sont constantes dans le VER, pour
donné.
L
E
T 1 U dip
ind ,VER ,int
p , , r
VVER
i
P F,G
(49)
1 U dip
ind ,VER ,int
pk U dip
ind ,VER ,int
k U dip
ind ,VER ,int
rk
VVER 1,n pk Pi 1,n k Pi 1,n 1,n rk Pi
0
136
pk k
On doit alors calculer les termes et .
Pi F,G Pi F,G
k
F et G, donc 0
Pi F,G
p V
• Montrons que k VER ki :
Pi F,G n
1 n pj j
fi
pk k
(P , p
)
Pi
Pi 1
ik 0
pk k V 1 pk k pk k V
Pi
ik
pk k VVER
(51)
n pi i
1 Pj
VVER ij (52)
D’autre part, en supposant que les dipôles sont égaux entre eux et sont aussi égaux à la valeur
moyenne macroscopique du dipôle, on pose : pk pk VER
137
Où correspond aux valeurs moyennes des vecteurs.
1
n
1 Pj Pj Pj Pj
k
VVER kj
On a alors :
p k
VVER
kj (53)
Pj n
E
L T
U dip
ind ,VER ,int
pj T r , r
n 1 n n
p j
avec
pk 1
ij
1 1
ij
Et donc :
p T r , r p j p j
1 n n n
L T 1
E ,r
n 1 1
i
et 1n ij
1
ij
(54)
1 n n 1 n
T r , r p j Eloc
n 1 1 ij n 1 , j
Notons qu’il ne semble pas anormal de trouver que le champ local (macroscopique et
continu) défini dans les articles de Maugin et ses collaborateurs soit différent de la
moyenne des champs locaux microscopiques et discontinus introduits dans les
approches atomistiques, puisqu’au chapitre 2, on a vu qu’en statique et en l’absence de
138
La constante diélectrique
L EiT
0 il
1
(55)
Pl F,G
n n
T r , r p j p j
n
1
1 1 ij
1
ij
L EiT 1 n
pk
(56)
Pl F,G n 1 pk P F,G
il n
VVER 1 n 1 n n
0 il
1
1
T r , r
n n 1 il (57)
1 1
On peut noter que cela ressemble bien à une sorte d’équation de Clausius-Mossotti
généralisée, puisque dans le cas d’un cristal (cubique par exemple) dans lequel tous les
atomes ont le même tenseur de polarisabilité6 :
𝑉𝑉𝐸𝑅 −1 𝐿̿
(𝛼̿)−1 = (𝜀0 (𝜀̿𝑟 − 1̿)) + (58)
𝑛 𝜀0
où 𝐿̿ est le tenseur des facteurs de dépolarisation de Lorentz, normalement défini par ([13]) :
1 (𝑟⃗ ′ − 𝑟⃗)
𝐿̿ = ∯ ⃗⃗⃗⃗𝑒𝑥𝑡
⊗ 𝑑𝑆′
4𝜋 𝑆𝛿(𝑟⃗) |𝑟⃗ ′ − 𝑟⃗|3 (59)
6
Dans le cas d’un cristal cubique avec des polarisabilités isotropes, on doit avoir 𝐿̿ = 1̿⁄3 pour retrouver
𝑉 (𝜀 −1)
l’équation de Clausius-Mossotti habituelle : 𝛼 = 3𝜀0 (𝜀𝑟
𝑁 𝑟 +2)
139
1 (𝑟⃗ ′ − 𝑟⃗) (𝑟⃗ ′ − 𝑟⃗)
𝐿̿ = ∯ ⃗⃗⃗⃗
⊗ 𝑑𝑆′𝑒𝑥𝑡 = ∭ ∇′ ⊗ ⃗⃗ 𝑑𝑉
4𝜋 𝑆𝛿(𝑟⃗) |𝑟⃗ ′ − 𝑟⃗|3 𝑉(𝑟⃗) 4𝜋|𝑟⃗ ′ − 𝑟⃗|3
= 𝜀0 ∭ 𝑇̿(𝑟⃗ ′ , 𝑟⃗)𝑑𝑉
𝑉(𝑟⃗)
A comparer avec ce que l’on peut obtenir par identification de (58) dans (57) :
n n
1 n n VVER
T r , r il 0 n T r , r
V
Lil 0 VER
n2
1 1
1 1
il n (61)
Tout se passe donc comme si nous calculions le facteur de dépolarisation comme moyenne
des intégrales discrétisées, calculées en prenant pour centre chaque atome du VER.
Le tenseur flexoélectrique
1 2U int L EiT
fijKL (62)
VVER Pi G jKL G jKL
1 n n
T r , r p j j
n 1
n 1 1
rk
ij
n ij
1 (63)
fijKL
1 rk G jKL
rk
Donc RKM RLM (65)
G jKL
140
Puisque le propagateur T est symétrique et le propagateur j T r , r est antisymétrique,
on obtient :
j T r , r T r , r
1
p k
fijKL j (66)
2n 1 1 1 ik ik
1
Dans cette égalité, il est utile de remarquer qu’il existe deux catégories de composantes : l’une
est proportionnelle aux dipôles induits et à leurs gradients p , p qui dépendent du champ
électrique extérieur et des dipôles permanents et l’autre est proportionnel aux dipôles
permanents et à leurs gradients qui sont fonctions de la déformation de la structure.
Les termes proportionnels à p ' correspondent aux termes suivant de l’énergie interne
exprimée en fonction des variables P et G , u Int P,G AhijKL P0hP0i G jKL qui reflète
u Int P perm ,G fijKL Pi perm G jKL qui reflète la flexoélectricité (c’est à dire une polarisation
donc les termes faisant intervenir les dipôles permanents seulement (dans une configuration
mixte) :
1 n n n
fijKL
j T r , r
2n 1 1 1
T r
k
,r R
j
k
M
K RLM (67)
ik ik
Insistons sur le fait que les dipôles induits ne peuvent pas contribuer aux effets
l’extérieur. Le coefficient ne pourrait donc pas être une constante intrinsèque du matériau. À
partir de la théorie, on remarque toutefois que l’effet flexoélectrique est dû au couplage entre
141
les dipôles induits par le champ électrique extérieur et les dipôles permanents dus à la
déformation du nanotube.
6 Conclusion
locale étendue, nous avons pu calculer les tenseurs macroscopiques purement mécaniques
comme les tenseurs des contraintes du premier ou du second ordre ou le tenseur d’élasticité, à
partir des variables microscopiques comme la position des atomes ou les forces qu’exercent
les atomes sur les autres. Pour la partie électrique, il a fallu mettre en œuvre une autre
hypothèse d’homogénéisation. Nous avons alors dû supposer que les variations de dipôles par
matériau, étaient toutes égales. C’est à partir de cette méthode que nous avons pu calculer les
1 U dip
p
ind ,VER ,int
T r , r ij p j ij p j
1 n n
n
1
L T
E ,r
VVER i (68)
i
et 1n
F ,G n 1 1 1
flexoélectricité. Nous en avons conclu que seuls les dipôles permanents créés par la courbure
peuvent contribuer qu’à des ordres supérieurs fonctions non linéaires du champ électrique. A
partir de ce chapitre, nous avons donc pu ajouter une nouvelle partie de code qui permet de
calculer chaque variable macroscopique à partir des variables utilisées durant la modélisation
moléculaire. Nous avons par exemple codé le calcul des tenseurs de contraintes, champs
142
les résultats de ce chapitre afin de caractériser les différentes constantes du matériau dans le
chapitre suivant.
143
Références
[1] J. L. Ericksen, « On the Cauchy—Born Rule », Math. Mech. Solids, vol. 13, no 3‑4, p.
199‑220, mai 2008.
[2] G. Zanzotto, « On the material symmetry group of elastic crystals and the Born rule »,
Arch. Ration. Mech. Anal., vol. 121, no 1, p. 1–36, 1992.
[3] Friesecke et Theil, « Validity and Failure of the Cauchy-Born Hypothesis in a Two-
Dimensional Mass-Spring Lattice », J. Nonlinear Sci., vol. 12, no 5, p. 445‑478, oct.
2002.
[4] E. B. Tadmor, M. Ortiz, et R. Phillips, « Quasicontinuum analysis of defects in solids »,
Philos. Mag. A, vol. 73, no 6, p. 1529‑1563, juin 1996.
[5] R. E. Miller et E. B. Tadmor, « The quasicontinuum method: Overview, applications and
current directions », J. Comput.-Aided Mater. Des., vol. 9, no 3, p. 203–239, 2002.
[6] V. B. Shenoy, R. Miller, E. B. Tadmor, D. Rodney, R. Phillips, et M. Ortiz, « An adaptive
finite element approach to atomic-scale mechanics—the quasicontinuum method », J.
Mech. Phys. Solids, vol. 47, no 3, p. 611–642, 1999.
[7] J. P. Lu, « Elastic properties of carbon nanotubes and nanoropes », Phys. Rev. Lett., vol.
79, no 7, p. 1297, 1997.
[8] G. Friesecke et R. D. James, « A scheme for the passage from atomic to continuum theory
for thin films, nanotubes and nanorods », J. Mech. Phys. Solids, vol. 48, no 6‑7, p. 1519‑
1540, juin 2000.
[9] M. Arroyo et T. Belytschko, « An atomistic-based finite deformation membrane for single
layer crystalline films », J. Mech. Phys. Solids, vol. 50, no 9, p. 1941–1977, 2002.
[10] R. Sunyk et P. Steinmann, « On higher gradients in continuum-atomistic modelling »,
Int. J. Solids Struct., vol. 40, no 24, p. 6877‑6896, déc. 2003.
[11] H. S. Park, P. A. Klein, et G. J. Wagner, « A surface Cauchy–Born model for
nanoscale materials », Int. J. Numer. Methods Eng., vol. 68, no 10, p. 1072‑1095, déc.
2006.
[12] Z. Wang, « Propriétés électro-mécaniques des nanotubes de carbone », Université de
Franche-Comté, 2008.
[13] A. D. Yaghjian, « Electric dyadic Green’s functions in the source region », Proc.
IEEE, vol. 68, no 2, p. 248‑263, 1980.
144
Chapitre 5
Dans le chapitre 3, nous avons expliqué comment calculer l’équilibre électrique d’un
système. Les dipôles d’un système d’atomes sont soit induits par un champ électrique
extérieur soit créé lors du changement d’hybridation d’orbitales causé par la courbure locale
du matériau. À partir de l’équilibre électrique, nous pouvons en déduire l’équilibre mécanique
à partir des forces entre les atomes données par la somme des contributions des interactions
« mécaniques » décrites par le potentiel AIREBO et des interactions électriques entre les
dipôles de la structure. Dans cette partie, nous simulerons la déformation d’un nanotube de
carbone semi-conducteur sous l’effet d’un champ électrique extérieur pour pouvoir calculer
les champs macroscopiques et les effets flexoélectriques. L’intérêt de choisir un nanotube de
carbone semi-conducteur est qu’il ne nécessite pas de calcul de distribution de charges
effectives comme ce serait le cas pour un nanotube métallique. C’est à partir des variables
microscopiques que nous en déduirons les champs macroscopiques ainsi que les constantes du
matériau (flexoélectriques, diélectriques, élasticité..).
atomes, d’une longueur de 12,39 nm, sous l’effet d’un champ 𝐸0,𝑥 = 𝐸0,𝑦 = 0,1 V·Å-1, après
avoir bloqué le déplacement des atomes sur le premier anneau. Les flèches représentent la
permanents et les dipôles induits. Bien que nous n’ayons pas pris en compte les dipôles
permanents dans l’optimisation, on va représenter ces dipôles en fig 32. Comme attendu, les
dipôles permanents induits par la courbure du nanotube ont toujours une direction normale à
146
la surface du nanotube et ont la même norme. En fig 33, les dipôles induits s’alignent selon le
champ extérieur 𝐸⃗⃗0 faisant un angle de 45°, avec l’axe horizontal. On remarque que les
dipôles sont d’autant plus importants qu’ils sont éloignés des bords et que les dipôles de la
partie inférieure du nanotube ont une pente moins forte par rapport aux dipôles de la partie
supérieure du nanotube.
fig 32. Représentation des dipôles permanents sur un nanotube (10,0) soumis à un champ
𝐸0,𝑥 = 𝐸0,𝑦 = 0,1 V·Å-1. Les positions à l’équilibre des atomes sont obtenues par
optimisation, sans prise en compte des dipôles permanents dans le calcul des dipôles
induits à l’équilibre.
147
fig 33. Dipôle induits à l’équilibre dans la même configuration que pour la fig 32.
Cette fois, on calcule les positions d’équilibre en prenant en compte les dipôles
permanents. Dans les figures ci-dessous, on représente les dipôles permanents, puis les
dipôles induits.
148
fig 34. Dipôles permanents avec prise en compte des dipôles permanents pour l’optimisation
fig 35. Dipôles induits avec prise en compte des dipôles permanents pour l’optimisation.
La fig 34 nous montre que les normes des dipôles permanents sont légèrement plus
importantes sur la partie supérieure du nanotube. Sur la fig 35, on peut voir que, cette fois-ci,
149
les dipôles induits pointent vers l’intérieur du nanotube puisqu’ils sont maintenant induits
aussi par les dipôles permanents.
Dans les graphiques ci-dessous, on trace, en fonction de la période, la somme sur cette
période des projections des dipôles permanents et induits selon les axes 1 (initialement
longitudinal) et 2. Pour les dipôles induits, on remarque une symétrie le long du nanotube
dans la distribution des composantes selon l’axe 1 des dipôles induits (cf. fig 36)
contrairement à ce qu’il se passe pour l’axe 2 (cf. fig 37).
Puisque les dipôles permanents sont orthogonaux à la surface du nanotube, on obtient une
contribution presque nulle en calculant la somme algébrique des projections des dipôles
permanents sur chaque période.
fig 36. Distribution des projections des dipôles sur l’axe horizontal, selon la position le long
de l’axe du nanotube
150
fig 37. Distribution des projections des dipôles sur la verticale, selon la position le long de
l’axe du nanotube
Ici, nous allons comparer les résultats de notre simulation de la déflexion d’un
obtenus par Zhao Wang durant sa thèse [1] (à cette époque, les dipôles permanents n’étaient
conducteur (8,0) soumis à un champ extérieur 𝐸⃗⃗0 faisant un angle de 45° après avoir bloqué
les premiers anneaux du nanotube pour modéliser un encastrement. En étudiant les variations
2
de en fonction du carré de la norme du champ extérieur ‖𝐸⃗⃗0 ‖ = 𝐸02 , il trouva finalement
appliqué°:
151
3.1 Simulation sans dipôles permanents
Dans un premier temps, nous faisons la même étude sans prendre en compte les
fig 38. Variation de l’angle de flexion d’un nanotube (8,0) de longueur 6.54 nm par rapport
au carré du champ appliqué avec un angle de 45° par rapport à l’axe du nanotube
lorsqu’il n’est pas déformé. Pour cette courbe, nous n’avons pas pris en compte les
dipôles permanents
Nous retrouvons bien la même loi de proportionnalité, mais avec 𝐴 = 0.0139 (V·Å-1)-2. Cette
différence s’explique probablement par le fait que nous n’avons pas pris les mêmes
paramètres de régularisation pour les atomes de bord, lors du calcul des propagateurs.
Étudions maintenant le sens physique de la relation (1). Nous pouvons démontrer que
lors du chapitre 3 :
𝑁
𝛿𝛼
∀𝛿 = 1, … , 𝑁, (𝑝⃗𝛿 )𝑒𝑞 + 𝜇⃗ 𝛿 = ∑(𝐴̂−1 ) (𝐸⃗⃗0 (𝑟⃗ 𝛼 ) + (𝛼̿ 𝛼 )−1 𝜇⃗ 𝛼 ) (2)
𝛼=1
152
où 𝐴̂ représente la matrice des propagateurs 𝑇
⃡(𝑟⃗ 𝛼𝛽 ). Or, on modélise la déformation du
nanotube dans le plan (𝑥, 𝑦) = (1,2). Donc, sans prendre en compte les dipôles permanents et
𝑁 𝑛 𝑛
1 1 𝛼𝛿
𝑃⃗⃗𝑖 = ∑ 𝑝⃗𝑖𝛼 = {∑ ∑(𝐴̂−1 )𝑖𝑗 } 𝐸0,𝑗 (3)
𝑉 𝑉
𝛼=1 𝛼=1 𝛿=1
1 𝛼𝛿
𝑃𝑖 = 𝜂𝑖 𝐸0 avec 𝜂𝑖 = {∑𝑛𝛼=1 ∑𝑛𝛿=1 ∑2𝑗=1(𝐴̂−1 )𝑖𝑗 } (4)
√2𝑉
De plus, dans le chapitre 4, on a montré que pour des angles 𝛼 suffisamment petits pour que
tan 𝛼 ≈ sin 𝛼 ≈ 𝛼, l’angle de déflexion a une relation de proportionnalité avec le gradient du
second ordre de la transformation G : 𝛼 = 𝐺112 𝐿. En utilisant les équations (1), (3) et (4), on
On alors montré que la relation (1) est similaire à l’équation constitutive traduisant la
« flexoélectrostriction », c’est-à-dire la proportionnalité entre le gradient du second ordre de
la transformation G et le carré de la polarisation électrique 𝑃⃗⃗.
Ainsi, si on trace 𝐺121 en fonction de 𝑃12 ou 𝑃22 , pour trouver 𝑎1 et 𝑎2 en faisant une
régression linéaire:
153
fig 39. Variation du gradient du second ordre de la transformation d’un nanotube (8,0) de
longueur 6.54 nm par rapport au carré de la polarisation du matériau suivant x
En faisant une régression linéaire, nous trouvons a1 2,6 108 m.C 2 . On a alors identifié un
terme qui caractérise l’induction d’un gradient de déformation due au carré d’une polarisation.
154
fig 40. Variation de l’angle de flexion d’un nanotube (8,0) de longueur 6.54 nm par rapport
au carré du champ appliqué avec un angle de 45° par rapport à l’axe du nanotube
lorsqu’il n’est pas déformé. Pour cette courbe, nous avons pris en compte les dipôles
permanents
Il est beaucoup plus difficile de converger lorsque nous prenons en compte les dipôles
permanents dans le modèle, nous avons alors réalisé une modélisation avec un champ
extérieur qui augmente beaucoup moins rapidement que lorsque nous ne prenions pas en
compte les dipôles permanents. On remarque que la pente a augmentée de 4.4% par rapport à
la pente représentée à la fig 38. L’influence des dipôles permanents n’est donc pas négligeable
dans la contribution de la déformation du nanotube. Si nous traçons le tenseur G par rapport à
la polarisation 𝑃1 , nous avons :
155
fig 41. Variation du gradient de transformation 𝐺121 d’un nanotube (8,0) de longueur 6.54
nm par rapport au carré de la polarisation du matériau suivant x. Pour cette courbe,
nous avons pris en compte les dipôles permanents
de p : 1211 1.38 107 C.m1 . Ce facteur semble être beaucoup plus important que l’ordre
de grandeur habituel que nous avons trouvé dans la littérature lors du chapitre 1 de
1010 C.m1 .
156
1 2 aireboelec 1 2 elec L Ei
fijKL (7)
V Pi G jKL V Pi G jKL G jKL
Comme nous l’avons vu précédemment, la courbure dû à la flexion est traduite par les
coefficients du tenseur du gradient du second ordre de la transformation G211 et G212 . Nous
nous intéresserons donc en particulier aux coefficients flexoélectriques f1211 , f1121 , f 2211 , f 2121
p 1 n n
ij p j
n
L
EiT
, r T r
, r
p
1
et 1n n 1 1 ij j
1 (8)
Puisque nous calculons la contribution flexoélectrique, nous calculons le champ local dû aux
dipôles permanents.
1 U dip
ind ,VER ,int
1 n n
, r T r , r j
L
EiT
et 1n V i n 1 1 ij (9)
Nous allons donc nous intéresser à la contribution des dipôles permanents sur la
flexoélectricité dans un nanotube de carbone semi-conducteur et mono-feuillet.
157
fig 42. Variation du coefficient f2211 d’un nanotube de carbone selon son rayon pour une
longueur constante de 20,6 nm
fig 43. Variation du coefficient f 2121 d’un nanotube de carbone selon son rayon pour une
longueur constante de 20.6 nm
158
Nous pouvons voir un effet de taille lorsque le rayon est inférieur à 10 Å, puis une
stabilisation lorsque le rayon dépasse 12 Å, avec une valeur d’environ 𝑓2121 = −1,9 V et de
𝑓2211 = −2,2 V.
le rayon du nanotube mais on fait varier sa longueur. Dans ce qui suit, nous calculerons le
coefficient flexoélectrique par une dérivée seconde évaluée sans déformation du nanotube,
159
fig 45. Coefficient flexoélectrique 𝑓2211 selon la longueur du nanotube (10,0)
Les coefficients commencent à être constants lorsque la longueur du nanotube atteint les
35 nm. Lorsque la longueur du nanotube est inférieure à 5 nm, l’effet de la longueur devient
important.
160
Conclusion et perspectives
Tout d’abord, nous avons utilisé la méthode des puissances virtuelles afin d’obtenir
avons appliqué les résultats généraux aux milieux semi-conducteurs déformables avec des
l’entropie, nous a donné l’inégalité très utile de Clausius-Duhem. Cette inégalité nous a
la flexoélectricité que nous pouvons trouver dans la littérature. Mais ce travail a aussi permis
d’avoir une théorie unifiée (écrite en unités SI et non en unités Lorentz-Heaviside comme
précédemment) qui peut traduire un plus grand nombre de phénomènes, ce qui permettra par
la suite d’être utilisé pour l’étude de beaucoup de phénomènes autres que la piézoélectricité
Nous avons ensuite présenté un modèle atomistique (codé en langage Fortran 9x)
l’approximation des dipôles-ponctuels distribués nous a permis de calculer les dipôles induits
résultants d’un champ extérieur appliqué sur la structure. C’est par l’interaction entre ces
dipôles induits que notre structure ne sera plus en équilibre mécanique. La déformation de la
structure est alors calculée à partir de l’interaction entre les dipôles électriques et du potentiel
AIREBO qui permet de rendre compte semi-empiriquement des liaisons entre les atomes
(covalentes, Lennard-Jones..). Enfin, nous avons utilisé les résultats de Dumitrica et al., de
161
Kalinin et Meunier et de Kvashnin et al. afin de pouvoir ajouter les dipôles permanents créés
par la courbure du matériau dans les interactions dipolaires. C’est en rajoutant ces dipôles
créés par une hybridation des liaisons 𝜎 et 𝜋 de chaque atome que nous pouvons modéliser les
Pour faire le lien entre ces 2 parties, il a fallu exprimer les variables macroscopiques
de la théorie des milieux continus en fonction des variables microscopiques données par les
modèles microscopiques en simulation atomistique. Pour cela, nous avons utilisé des
à-dire que nous avons supposé que le second gradient de la transformation était constant dans
un VER. Nous avons pu établir, à partir des données du modèle atomistique, les tenseurs
tenseur d’élasticité. Pour la partie électrique, nous avons fait une hypothèse
C’est à partir de cette méthode que nous avons pu calculer les tenseurs macroscopiques tels
que le champ électrique local ou la constante diélectrique du matériau. Enfin, à partir de nos
dipôles homogène dans un VER), nous avons pu obtenir les tenseurs de couplages
électromécaniques comme le tenseur de piézoélectricité ou de flexoélectricité. Nous en avons
conclu que seuls les dipôles permanents créés par la courbure du matériau pouvaient
contribuer à la flexoélectricité et que les dipôles induits qui dépendent du champ électrique
extérieur ne peuvent contribuer qu’à des ordres supérieurs décrivant, par exemple, des
strictions lors de déformation linéaires ou de flexion. Enfin, à partir des équations que nous
avons construites, nous avons pu illustrer numériquement comment il est possible de calculer
Le nombre d’atomes utilisés dans la simulation est limité par des contraintes de
temps d’exécution (qui varie approximativement comme le cube du nombre d’atomes) et de
162
mémoire vive disponible (la quantité nécessaire varie approximativement comme le carré du
nombre d’atomes). Pour modéliser une structure plus grande, il serait souhaitable de
modéliser un couplage entre la simulation atomistique et la méthode des éléments finis. Dans
notre étude, nous avons seulement traité le problème d’un matériau semi-conducteur dans le
vide. Il serait néanmoins possible d’élargir le modèle à des matériaux qui seraient métalliques
ou ioniques en ajoutant la distribution des charges dans notre modèle. De plus, comme nous
l’avons vu dans le chapitre 3, il serait possible de modéliser l’influence des surfaces qui
pourrait être due à l’encastrement. Dans l’état de l’art, nous avons vu que des quadripôles
étaient utilisés par Ong et Vanderbilt dans leurs calculs ab-initio des composantes des tenseurs
caractéristiques de la flexoélectricité. Dans l’un de ses articles, Maugin indique que les
modèles avec des gradients de polarisation (et d’aimantation) dipolaire sont alternatifs aux
modèles avec des multipôles d’ordre plus élevés et donnent des résultats très similaires. Il
Un autre projet ambitieux pourrait être l’ajout d’un modèle prenant en compte les
dipôles magnétiques dans le système pour pouvoir étudier des phénomènes de couplage
mécano-électromagnétiques. C’est à partir de cette étude que nous pourrions faire d’autres
recherches sur des domaines similaires à notre travail, comme l’étude du flexomagnétisme ou
de la magnétostriction.
Pour finir, notons que des expériences sur des nanotubes de carbone sont en cours au
flexoélectricité.
163
Annexes
A. Théorème de la divergence surfacique
Rappelons le théorème de la divergence surfacique donné par Collet [15] (voir aussi
[22] p. 22). On suppose que est une surface fermée traversée par une arête continue G
donnant le contour de , sur lequel on peut observer une discontinuité du vecteur unitaire
extérieur normal n . Sur chaque point G dans la fig 3, on peut introduire le vecteur unitaire
fig 46. Définition des vecteurs t et t tangent à l’arête G où le vecteur normal unitaire
+ -
Le théorème de la divergence surfacique est appliquée sur la surface et :
ˆ i qi da ni qi (ˆ j n j )da b q ds
i i
(A1)
165
ˆ i qi da ni qi (ˆ j n j )da b q ds
i i
(A2)
(
Ñ̂ j · = Pjk Ñk · = d jk - nj nk Ñk · ) (A3)
¶·
= nk Ñ k · (A4)
¶n
ˆ q da ˆ n )da
i i
ni qi ( j j
bi qi ds
(A5)
ij Dij Ei ˆi Bi ˆ i
* L * L * L
ij ˆij
* L
ij ˆ ij
(A7)
166
Puisque est un tenseur symétrique, on obtient s ij Dij* = s ij vi*, j . Donc,
P
(i ) ( , V
obj ) p(i ) d
L Eiˆi L Bi ˆ i L
ij ˆij* L ij ˆ ij d
P
(i ) ( , V
obj ) ( ij vi ), j ij , j vi ( ijk vi*, j ),k ( ijk ,k vi* ), j ijk , kj vi* d
L Ei i *ij j L Bi i *ij j (A10)
L ij i ik k , j L
,j ij
i ,j ik k , j d
Utilisant les identités suivantes :
Ai B j ij
1
2
Ai B j Aj Bi vi, j Ai B j vi, j
P
(i ) ( , V
obj ) p(i ) d
ij ijk ,k L E[i j L B[i j L [i k j ,k L [i k
j ,k ,j
vi*d
L Ei L
ij , j
i
L
Bi L
ij , j
i
d
L
ij n j i L
ij n j i da (A12)
ij ijk ,k L E[i j L B[i j L [i k j ,k L [i k
j ,k n j vi*da
ijk vi, j nk da
167
Afin de développer le dernier terme de l’équation (avec le facteur vi , j sur
ˆ P n n
avec et nk k (A14)
n
j jk k jk j k k
Le produit interne entre Tij et vi*, j donné par (A13) nous donne :
vi* ˆ ˆ T T n vi
*
ˆ v* T n
Tij vi*, j Tij j Tij vi* vi* (A15)
n n
j i ij j j ij ij j
ˆ j Tij vi* da
n jTij vi* ˆ p n p da b jTij vi* ds (A16)
- -
*ˆ vi*
Tij vi*, j da ˆ n da
n j Tij vi* p p b jTij vi* ds
T
i j ij ij j
v T n
n
da
(A17)
ˆ
ˆ *
j n j p n p Tij Tij n j n vi da b j Tij vi ds
*
168
P
(i ) ( , V
obj
) ij ijk ,k L E[i j L B[i j L [i k j ,k L
[i k
j ,k vi*d ,j
L Ei L ij , j i L Bi L
ij , j
i
d
L
ij n j i L
ij n j i da
ij ijk ,k L E[ i L B[ i L
j j [i k j , k
L
[i k
j ,k n j vi*da
j
ˆ n
j
ˆ n
p p ijk nk vi*da
ijk n j nk
vi*
n
da ijk nk b j vi*ds
(A18)
ij L E[i j L B[i j L [ik j ,k L
[i k
j ,k ,j
vi * qf
L *
i i
u d
E L
i
L
ij , j
i
L
B
i
L
ij , j
i
d
L
ij n j i L
ij n j i da
ij
L
E[i j L B[i j L [ik j ,k L
[i k
j ,k n j vi *da
p(i ) ij Dij ijk vi*, jk L Eiˆi L Bi ˆi L ijˆij* L ij ˆ ij qf L i ui *
(A19)
D E ˆ B ˆ
ij
*
ij
L
i
*
i
L
i
*
i
L
ij
*
ˆ
ij
L
ij
ˆ
ij
et d’après (134) :
qem (
. f .u ) (A20)
On a donc :
169
p(i ) qem ij Dij ijk vi , jk L Ei ( i ik k , j ) L Bi ( i ik k , j ) L
ij ( ) i ,j ik k , j
L ij ( )
i ,j ik k , j qf L
i v i
vi
ij Dij L
Ei ( i ij j ) L
Bi ( i ij j ) L
ij ( ) i ,j ik k , j L ij
( )
i ,j ik k , j (A22)
(
. f .u )
Puisque
On obtient :
p(i ) qem ij ( L E[ i j L B[ i j ) ( L [i ,k j ,k + L [i ,k
j ,k ) vi , j ijk vi , jk L Ei ( i )
L Bi ( i ) L ij ( i ), j L ij ( ) q
i ,j
L
f i v
i vi
ij ( L E[i j L B[i j ) ( L [i ,k j ,k + L
[i ,k
j ,k ) vi, j (A24)
L
Ei ( i ) L Bi ( i ) L ij ( i ), j L
ij ( )
i ,j
(
. f .u )
f qf
u v v
qf
(A26)
et donc :
f u
170
p(i ) qem ij ( L E[i j L B[i j ) ( L [i ,k j ,k L [i ,k
j ,k ) vi , j ijk vi , jk L Ei ( i )
L Bi ( i ) L ij ( i ), j L ij ( ) q
i ,j
L
f i v
i vi
L Ei ( i ) L Bi ( i ) L ij ( i ), j L
ij ( )
i ,j
(
.
)
Donc :
p(i ) qem ij ( L E[i j L B[i j ) ( L [i ,k j ,k L [i ,k
j ,k ) vi , j ijk vi , jk L Ei ( i )
L Bi ( i ) L ij ( i ), j L ij ( ) q
i ,j
L
f i vi
L Ei ( i ) L Bi ( i ) L ij ( i ), j L
ij ( )
i ,j
(
.
)
tij, j qf ( i
L i ) 0 (A30)
On obtient ainsi :
q f
L
i vi tij, j qf i
v i (A31)
t ij , j qf i
v q
i
f i
v t v t
i
ij i ,j
ij vi , j (A32)
On obtient alors :
171
p(i ) qem ij ( L E[i j L B[i j ) ( L [i ,k j ,k L [i ,k
j ,k ) tij vi , j ijk vi , jk L Ei ( i )
Bi ( i )
L L
ij ( ) ( ) q
i ,j
L
ij i ,j
f i
v t v
i
ij i ,j
L
Ei ( i ) L
Bi ( i ) L
ij ( )
i ,j
L
ij ( ) i ,j
(
.
)
(
.
)
fi vi ij i
f v ( i
i
.
)
,j
(A35)
.(
) ij i (
)
,j
.(
) .
. (
)
.( ) .
. ( eff )
.
On obtient au final :
172
pi q em tijT vi , j ijk vi , jk L Ei i L Bi i L
ij i , j L ij i , j
L
Ei i L
Bi i L
ij i, j L
ij i , j
(A36)
( eff )
.
.
On a démontré que :
Avec
(A39)
t ij ( E[i j )
T
ij
L
d d 1
tijT vi , j ijk vi , jk L Ei i q . 0
i i
dt dt
Le but est alors de calculer vi , jk et vi , j par rapport à FiJ et GiJK . Voici le détail des calculs :
173
d d d
F xi , K GiKL xi , KL
dt dt dt
vF vi , KL
v FF 1 (vi , j x j , K ), L
vi , j FiJ X J , j vi , jL x j , K vi , j x j , KL
vi , jk xk , K x j , L vi , j x j , KL
vi , jk FkK FjL vi , j G jKL
d
GiKL vi , jk FkK FjL vi ,l GlKL
dt
GiKL X Kj X Lk vi , jk vi ,l GlKL X Kj X Lk
vi , jk GiKL X Kj X Lk vi ,l GlKL X Kj X Lk
vi ,kj GiKL X Kk X Lj vi ,l GlKL X Kk X Lj
Donc :
d d T 1
dt
dt
tij vi , j ijk GiKL X Kk X Lj vi ,l GlKL X Kk X Lj Ei i
L
i i q . 0
Puisque ijk vi ,l GlKL X Kk X Lj ilk vi , j G jKL X Kk X Ll (A40)
d d 1
dt
dt
T
L
tij ilk G jKL X Kk X Ll FiJ X J , j ijk GiKL X Kj X Lk Ei i i i q . 0
On a :
d dFiJ dGiJK d i d
dt FiJ dt GiJK dt i dt dt
174
d d T
tij ilk G jKL X Kk X Ll X J , j FiJ ijk X Kj X Lk GiKL
d dt FiJ GiKL
(A41)
1
L Ei i i q . 0
i
i
La condition nécessaire et suffisante pour que cette inégalité soit vérifiée est que :
FiJ
tijT ilk G jKL X Kk X Ll X J , j ij ( L E[i j ) ilk G jKL X Kk X Ll X J , j
ijk X Kj X Lk QiKL
GiKL (A42)
L Ei
i
E. Glossaire
RC = Référentiel en co-mouvement
On rappelle qu’il existe d’autres définitions du potentiel chimique, noté m , dans d’autres
a
contextes, mais nous avons pris la même définition que Daher et al. dans [19].
Table 1: Liste des variables utilisées avec leur signification et leur dimension en unités SI
175
0 Permittivity of vacuum ( 1/ [c 2 0 ]) M 1.L3 .I 2 .T 4
1 if i = j
ij or KL Kronecker delta such that : ij 1
0 otherwise
Material domain L3
176
Gradient operator in Kt (i.e. derivative with respect to the xi
L1
variables)
q ( eff ) Effective charge density of the specie per unit volume I .T .L3
eff
i = Effective electromotive field referring to RC (also called
M .L.I 1.T 3
"Effective electromotive intensity" by Maugin)
177
i Macroscopic magnetic induction referring to RC M .I 1.T 2
i Conduction current density of charge carrier referring to RC (or
current per unit area) I .L2
(eff)
Effective conduction current of the charge carrier referring
I .L2
i
i Macroscopic electric field referring to RC for the th continuum of M .L.I 1.T 3
charge
178
(eff)
i = Effective electromotive field for the specie M .L.I 1.T 3
L
Bi Local magnetic induction reflecting spin - lattice interaction M .I 1.T 2
L
ij Second order tensor of spin spin interaction or exchange force M 2 .I 1.T 2 .L2
L
Bi Local magnetic induction of the th continuum of charge M .I 1.T 2
L
BiT Total local magnetic induction M .I 1.T 2
L
ij Second order tensor of the spin spin interaction or exchange
M 2 .I 1.T 2 .L2
force of the th continuum of charge
L
ij
T
Second order tensor of the total spin spin interaction or total
M 2 .I 1.T 2 .L2
exchange force
L
Ei = Local Electric field M .L.I 1.T 3
L
ij Second order tensor of the shell shell interaction force M 2 .I 1.L1.T 3
M .L.I 1.T 3
L
Ei Local electric field of the th continuum of charge M .L.I 1.T 3
L
EiT Total local electric field M .L.I 1.T 3
M 2 .I 1.L1.T 3
M 2 .I 1.L1.T 3
179
ij First order intrinsic symmetric stress tensor acting on the specie M .L1.T 2
fi Body force per unit volume acting on the specie M .L2 .T 2
M .T 2
Li (and Li ) Linear traction on (and the reduced one) (per unit length )
180
tijT = Total symmetric stress tensor M .L1.T 2
M .L2 .T 2
h = Radiation heat input per unit mass (i.e. power sources different from
those included in p(i ) and qem ) L2 .T 3
181
i Entropy flux vector (power per unit area and temperature) M .T 3 . 1
Pi Convective - time derivative of the polarization I .L2
(Unit of a )
aˆi DJ i Jaumann derivative of the vector
.T 1
(Unit of a )
aˆij Specific time derivative tensor of the gradient of the vector
. L1.T 1
182
In the reference configuration KR Dimension
PK Polarization I .T .L2
MK Magnetization I .L1
183
E K (eff) = Effective electric field (or electromotive intensity) for the
M .L.I 1.T 3
specie
E KL
L
Second order tensor of shell shell interaction tensor for
M 2 .I 1.L1.T 3
the specie
E KL
L T
Second order tensor of total shell shell interaction tensor M 2 .I 1.L1.T 3
L
BK Local magnetic induction M .I 1.T 2
L
BK Local magnetic induction for the specie M .I 1.T 2
L
BKT Total local magnetic induction M .I 1.T 2
L
BKL Second order tensor of spin spin interaction or the exchange
M 2 .I 1.T 2 .L2
force
L
BKL
Second order tensor of spin spin interaction or the exchange
M 2 .I 1.T 2 .L2
force for the specie
L
BKL
T
Second order tensor of total spin spin interaction or
M 2 .I 1.T 2 .L2
the total exchange force
184
12
KIJ Second order tensor of dielectric susceptibilities M 2 .L1.I 2 .T 4
22
KLIJ Polarization gradient - electric gradient coupling tensor M 2 .L4
32
cHIJKL Second order tensor of elastic moduli M .T 2
33
cKLMNOP Strain gradient strain gradient coupling tensor M .L.T 2
31
f HIJK Tensor of flexoelectric coefficients M .L2 .I 1.T 3
22
f IJKL Gradient polarization strain coupling tensor M 2 .L1.I 1.T 3
32
f HIJKL Gradient polarization strain gradient coupling tensor M 2 .I 1.T 3
1
Chemical potential polarization coupling tensor
K M 2 .L.I 2 .T 4
2
KL Chemical potential polarization gradient coupling tensor M 3 .L1.I 2 .T 4
2
Chemical potential strain coupling tensor
KL M 2 .L1.I 1.T 3
3
KLM Chemical potential strain gradient coupling tensor M 2 .I 1.T 3
M 2 .L1.I 1
Chemical potential temperature coupling factor
.T 3 . 1
185
2 0
TKL First order stress tensor at the natural state M .L1.T 2
3 0
TKLM Second order stress tensor at the natural state M .T 2
E KL0 The second order local electric field at the natural state
L
M 2 .L1.I 1.T 3
0 Helmholtz free energy per unit volume at the natural state M .L1.T 2
une conversion des équations des variables SI aux équations avec des variables données en
186
Table 2: Conversion des variables pour passer des unités CGS Lorentz-Heaviside aux unités SI
CGS Lorentz-Heaviside SI
1
e0
(d, c KL
, 12 c KLI , 22 c KLIJ ,d a , 1g aK , 2g aKL ) (d, c KL
, 12 c KLI , 22 c KLIJ ,d a , 1g aK , 2g aKL )
187
G. Lien entre G et la courbure et l’angle de déflexion
1⁄
𝑅
On trouve :
G L
tan 112
2 2
1
Soit R le rayon de courbure, avec G211 si la courbure reste perpendiculaire à sa section
R
alors :
R L ,
Pour un angle faible,
G L
tan 112
2 2 2
1
Donc G112 G211
R
188
Titre : Etude de la fléxoélectricité de nanosystèmes par le développement d'algorithmes mêlant approche
atomistique et mécanique des milieux continus.
Title : Study of flexoelectric nanosystems through the development of multi-scale algorithms mixing atomistic
approach and continuum mechanics