# De l'atome à la pièce: une promenade numérique au travers des différentes échelles de la rhéologie

Francisco Chinesta

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

## To cite this version:

Francisco Chinesta. De l'atome à la pièce: une promenade numérique au travers des différentes échelles de la rhéologie. Rhéologie, 2009, 15, pp.1-30. hal-01006717

HAL Id: hal-01006717

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

Submitted on 14 Dec 2016

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# De l'atome à la pièce : une promenade numérique au travers des différentes échelles de la rhéologie

F. Chinesta

EADS Corporate Foundation International Chair, GeM : CNRS - Centrale de Nantes
1 rue de la Noe, BP 92101, F-44321 Nantes cedex 3

Francisco.Chinesta@ec-nantes.fr

**Résumé** : La modélisation des fluides complexes nécessite une approche multi-échelle, multi-physique et multi-modèle, ainsi que la définition de passerelles permettant le passage entre les différentes échelles afin de pouvoir décrire le comportement rhéologique aux échelles nanométrique, micrométrique ou encore macroscopique, ainsi que de pouvoir définir le comportement à chaque échelle à partir de ceux intervenants aux échelles les plus fines. Cet article passe en revue les différents types de description ainsi que les passerelles reliant les différentes échelles. Afin de pouvoir élargir la portée de ce travail, j'ai volontairement décidé d'essayer de rendre ce document accessible aux non-spécialistes. Ce document retrace mes principaux travaux de recherche dans le domaine de la rhéologie numérique. La plupart d'entre eux, si ce ne sont pas tous, sont faits en collaboration avec mon collègue et ami Amine Ammar du Laboratoire de Rhéologie de Grenoble, à qui j'adresse ma gratitude dès le début de ce document pour lequel je le considère comme co-auteur.

Ce travail est dédié à la mémoire de Maurice Couette, Professeur, chercheur et homme de grande culture, comme l'ont brillamment défini Jean-Michel et Monique Piau (Piau et Piau, *Rhéologie*, 8, 1-4, 2005), avec tout l'honneur que me confère le fait d'avoir été le récipiendaire d'un prix qui porte son nom.

**Mots-clé** : Rhéologie numérique, Echelle quantique, Dynamique moléculaire, Théorie cinétique, Suspensions, Polymères, Fluides complexes, Fluides micro et nano-structurés, Fokker-Planck, Schrödinger.

[Abridged English version on last page]

## 1. Introduction

Mon collègue Malcolm Mackley de l'Université de Cambridge répond à ses étudiants quand ceux-ci l'interrogent à propos du contenu de leurs conférences : « Il faut simplement raconter une belle histoire ! ». Voilà le défi : raconter une belle histoire, simple, mais tout un défi !

Dans cet article, je vous propose un voyage qui nous amènera de l'échelle quantique, là où les liaisons moléculaires sont présentes, jusqu'à l'échelle macroscopique, où se situent la plupart des écoulements industriels. Je n'oublierai pas également de prévoir quelques arrêts au niveau nanométrique et micrométrique : les micro et nanosciences, ainsi que les micro- et nano-technologies, nous y obligent !

Toute histoire nécessite un début et l'établir est loin d'être un choix simple. J'ai donc décidé de débuter mon histoire dans l'ancien laboratoire Cavendish, au centre ville de Cambridge, il y a un peu plus d'un siècle. Comme l'annonce une plaque fixée sur l'un des murs de l'ancien laboratoire (Fig. 1), Joseph John Thomson y a découvert l'électron en 1897. Ce

fait marquant est à l'origine d'une nouvelle physique ayant d'importantes conséquences dans le domaine de la rhéologie, comme nous allons le voir par la suite.

<figure><img src="image_3.png" /><figcaption>Figure 1. La découverte de l'électron en 1897 marque la naissance de l'ère atomique</figcaption></figure>
## 2. L'échelle atomique

### 2.1 Introduction

Avec la découverte de l'électron, il a fallu établir un modèle atomique. Le plus primitif, celui de Thomson, représente l'atome comme une « pastèque », où la chair devait être positive pour équilibrer la charge négative des « pépins » (les électrons). La découverte d'autres particules, cette fois-ci positives, lors de désintégrations radioactives (voir les travaux de Becquerel, Pierre et Marie Curie, Rutherford ...) a conduit à des modèles atomiques plus sophistiqués. L'un des pionniers était le modèle de Rutherford, dans lequel les électrons tournaient autour du noyau massif et positif, comme le font les planètes autour du soleil avec des trajectoires parfaitement décrites par les équations d'Isaac Newton. Cependant, un problème surgit du fait que toutes charges électriques accélérées (comme les électrons autour du noyau) émettent un rayonnement et perdent par conséquence de l'énergie. Ainsi, la seule possibilité qui s'offrait à ces électrons tournant dans une espèce de système « planétaire » était de se diriger inexorablement vers le noyau.

D'autres modèles ont vu le jour. Avec les travaux de Max Planck qui a quantifié l'énergie (qui ne sera plus une variable continue), les gains ou pertes d'énergie doivent nécessairement être un multiple du paquet fondamental d'énergie, le quanta, de valeur $h$ (la constante de Planck). Ce résultat a permis à Niels Bohr de proposer un modèle atomique stable, puisque seule l'émission ou l'absorption d'une certaine quantité précise d'énergie permet à l'électron de changer d'orbite. Bien évidemment, l'histoire est longue, et de nombreuses améliorations de ce modèle furent proposées, mais nous en resterons là.

En 1905, Albert Einstein publie trois articles dont il est difficile de savoir lequel des trois était le plus décisif!

Dans l'un de ces articles, traitant l'effet photoélectrique, Einstein tranche une ancienne dispute entre les défenseurs du caractère particulaire de la lumière (dont Isaac Newton était le porte-drapeau !) et les partisans de sa nature ondulatoire (Christian Huygens, Thomas Young, Augustin Fresnel...). Définitivement, la lumière doit être à jamais conçue avec une nature ondulatoire ou corpusculaire. Ces deux descriptions sont pertinentes et décrivent une réalité qui possède les deux facettes : l'effet photoélectrique prouve son caractère corpusculaire (la lumière est composée de photons) et les

phénomènes de diffraction et d'interférence sont le constat de sa nature ondulatoire.

L'éminent français Louis de Broglie va étendre cette double nature de la lumière à la matière, en sentant les bases de la théorie ondulatoire de la matière. Ainsi, toutes les particules élémentaires, dont l'électron fait partie, n'ont pas seulement une nature corpusculaire, mais elles ont également un caractère ondulatoire qui devrait générer des phénomènes d'interférence. Il en est convaincu, mais la preuve vient encore de l'autre côté de la Manche, du fils de Joseph John Thomson, George Thomson, qui prouve de façon incontestable la nature ondulatoire de la matière. Lui, le fils, tout comme son père, ont reçu le prix Nobel de physique. Le père, Joseph John, l'a eu pour avoir prouvé la nature corpusculaire des électrons, tandis que le fils, Georges, l'a eu pour avoir prouvé le contraire, sa nature ondulatoire. Voici une dispute entre père et fils qui se solde par deux prix Nobels !

Le deuxième article d'Albert Einstein de 1905 porte sur le mouvement brownien (nous y reviendrons quand nous nous intéresserons à la théorie cinétique). Décisif comme son premier, celui-ci prouvait la nature atomique de la matière. Certes, cette nature atomique de la matière était entrevue d'après les travaux en chimie de Boyle, Lavoisier, Dalton, Gay Lussac, Avogadro... Les approches les plus mécanistes de cette théorie atomique ont pris du terrain grâce aux travaux de Maxwell et Boltzmann. Cependant, les travaux de ces deux génies n'ont pas fait à l'époque l'unanimité, loin de là. Supposer cette nature atomique posait pas mal de problèmes. Premièrement, la physique macroscopique semblait fortement non-linéaire et surtout irréversible : les travaux de Sadi Carnot en étaient la preuve et ont permis d'aboutir au deuxième principe de la thermodynamique! Par contre, à l'échelle atomique tout semblait trop simple et, ce qui est encore pire, absolument réversible. On retrouve cette réversibilité dans les simulations de dynamique moléculaire. Enfin, l'article d'Einstein s'est imposé trop tard puisque Boltzmann s'est suicidé en 1906, suite au discrédit que sa théorie cinétique lui a rapporté. Cependant, sur sa tombe, il a fait graver la définition microscopique de l'entropie ! L'histoire lui donnera raison !

Finalement, le troisième article d'Einstein de 1905 introduit la relativité restreinte, écartant l'existence de référentiels absolus et établissant la vitesse de la lumière comme limite supérieure des vitesses et son indépendance du référentiel considéré. Les transfor-
rmations de Lorentz déplacent définitivement celles
de Galilée et les lois de l'électromagnétisme de
Maxwell retrouvent finalement la paix, après avoir
eu de sérieux soucis avec les transformations de
référentiels.

Dix ans plus tard, Albert Einstein revient en scène
avec sa plus grande création, la théorie générale de
la relativité, mais nous en resterons là.

## 2.2 L'équation de Schrödinger

En mécanique classique, l'évolution d'un système nécessite la connaissance de la position et de la vitesse, toutes les deux étant supposées accessibles avec une précision illimitée. En mécanique quantique, les choses s'annoncent différentes à cause du principe d'incertitude de Werner Heisenberg qui a établi une limite à la détermination simultanée de la vitesse et de la position.

Malgré l'opposition de scientifiques de renoms dont
fait partie Einstein, le principe d'incertitude s'est
finalement imposé : l'expérience d'Alain Aspect en
1975 est sans appel !

En bref, toute particule atomique doit être distribuée
dans l'espace physique (comme le serait une onde)
et sa vitesse ainsi que sa position ne sont pas
définies avant d'effectuer une mesure. C'est à
l'instant de la mesure que la fonction d'onde se
matérialise comme une particule dont la vitesse et la
position ne peuvent pas être déterminées au-delà
d'une précision limite.

Ainsi, au lieu de parler de la position d'une
particule, nous devrions plutôt parler de la
probabilité de la trouver, lors d'une mesure, dans
une position donnée **x** au temps **t**. Cette probabilité
peut être notée par :

$$P_x(\mathbf{x}, t) \tag{1}$$

et de même en ce qui concerne les vitesses :

$$P_p(\mathbf{p}, t) \tag{2}$$

où l'on préfère habituellement utiliser la quantité de
mouvement **p** (**p** = **m****v**, **m** étant la masse) plutôt que
la vitesse **v** elle-même. Etant positives, ces deux
probabilités peuvent être réécrites comme :

$$P_x(\mathbf{x}, t) = |\Psi_x(\mathbf{x}, t)|^2 \tag{3}$$

$$\text{et } P_p(\mathbf{p}, t) = |\Psi_p(\mathbf{p}, t)|^2 \tag{4}$$

où $\Psi_x(\mathbf{x}, t)$ et $\Psi_p(\mathbf{p}, t)$ sont respectivement les
fonctions d'onde associées à la position et à la
quantité de mouvement.

Cependant, ces deux fonctions d'onde ne sont pas
indépendantes, le principe d'incertitude de Heisen-
berg en est la cause. On peut prouver que l'une est la
transformée de Fourier de l'autre [1]. Il est bien
connu que la transformée de Fourier d'une fonction
parfaitement localisée en espace fait intervenir
toutes les fréquences et, réciproquement, une
fonction qui ne fait intervenir qu'une seule
fréquence sera localisée partout dans l'espace. Cela
pose un « petit » problème puisqu'en mécanique
classique le système est déterminé par deux
variables (position et vitesse), tandis qu'en
mécanique quantique nous n'avons qu'une fonction
d'onde indépendante, soit $\Psi_x(\mathbf{x}, t)$, soit $\Psi_p(\mathbf{p}, t)$. Ce
problème peut être résolu si l'on suppose que la
fonction d'onde contient une partie réelle et une
partie imaginaire. Par la suite, nous n'allons
considérer que $\Psi_x(\mathbf{x}, t)$ avec :

$$P_x(\mathbf{x}, t) = |\Psi_x(\mathbf{x}, t)|^2 = \bar{\Psi}_x(\mathbf{x}, t) \cdot \Psi_x(\mathbf{x}, t) \tag{5}$$

Si l'on suppose connue la fonction d'onde, on peut
calculer la valeur moyenne de la position. En nous
plaçant, par simplicité, dans un espace uni-
dimensionnel, on aura pour la position :

$$\langle x \rangle = \int_{-\infty}^{+\infty} x P_x(x,t) dx = \int_{-\infty}^{+\infty} x \bar{\Psi}_x(x,t) \cdot \Psi_x(x,t) dx \tag{6}$$

et en appliquant une propriété de la transformation
de Fourier et le théorème de Parseval, on obtient
pour la quantité de mouvement [1] :

$$\langle p \rangle = \int_{-\infty}^{+\infty} p P_p(x,t) dp = \int_{-\infty}^{+\infty} \bar{\Psi}_x(x,t) \frac{\hbar}{i} \frac{\partial \Psi_x(x,t)}{\partial x} dx \tag{7}$$

(où $\hbar$ est égal à $h/2\pi$, $h$ étant la constant de Planck)
ou encore pour l'énergie :

$$H(\mathbf{x}, t) = \frac{\mathbf{p}^2}{2m} + V(\mathbf{x}) \tag{8}$$

(où $V$ est l'énergie potentielle), dont la valeur
moyenne est donnée par :

$$\langle H \rangle = \int_{-\infty}^{+\infty} \bar{\Psi}_x \left( -\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x} + V(x) \right) \Psi_x dx \tag{9}$$

Ces valeurs moyennes sont reliées par le théorème
d'Ehrenfest :

$$\langle p \rangle = m \frac{d \langle x \rangle}{dt} \tag{10}$$

qui a une forte ressemblance avec l'équation
correspondante de la mécanique classique.
Si le système est composé de $N_e$ électrons, la fonction d'onde s'écrit :

$$ \Psi(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_{N_e}, t) \qquad (11) $$

où l'indice $x$ précédemment associé à la fonction d'onde est volontairement omis pour alléger la notation. Une fois que cette fonction sera connue, on peut facilement déduire la distribution spatiale de l'électron $\rho_e(\mathbf{x}, t)$, par intégration :

$$ \rho_e(\mathbf{x}, t) = \int_{\mathbb{R}^{3(N_e-1)}} |\Psi|^2 d\mathbf{x}_1 \cdots d\mathbf{x}_{e-1} d\mathbf{x}_{e+1} \cdots d\mathbf{x}_{N_e} \qquad (12) $$

Ainsi, la distribution électronique totale qui en résulte s'écrit :

$$ \rho(\mathbf{x}, t) = \sum_{e=1}^{e=N_e} \rho_e(\mathbf{x}, t) \qquad (13) $$

Cette expression est d'une importance capitale dans la description des propriétés des matériaux. Cependant, il nous manque une brique essentielle, l'équation permettant le calcul de la fonction d'onde. Erwin Schrödinger proposa en 1926 une équation permettant de la calculer. Pour un système composé d'un seul électron et d'un noyau, ce dernier étant fixé dans $\mathbf{X}$ (modèle de Born-Oppenheimer), la forme adimensionnelle de l'équation de Schrödinger s'écrit :

$$ i \frac{\partial \Psi(\mathbf{x}, t; \mathbf{X})}{\partial t} = -\frac{1}{2m} \nabla^2 \Psi(\mathbf{x}, t; \mathbf{X}) + V_{\text{en}} \Psi(\mathbf{x}, t; \mathbf{X}) \qquad (14) $$

où $V_{\text{en}}$ est le potentiel d'interaction coulombien électron-noyau :

$$ V_{\text{en}} = \frac{q_e \cdot q_n}{\|\mathbf{x} - \mathbf{x}_0\|} \qquad (15) $$

où la charge de l'électron est $q_e = -1$, tandis que celle du noyau correspond à son numéro atomique : $q_n = +Z$. Dans (14), $m = 1$ pour les électrons.

Cette équation est, certes, loin d'être évidente mais elle est censée être un premier principe, comme l'avait été deux siècles auparavant l'équation fondamentale de la dynamique de Newton $\mathbf{F} = m\mathbf{a}$. Tout comme celle de Newton, l'équation de Schrödinger ne dérive d'aucune autre équation et sa validité ne peut se prouver que par un accord entre ses prédictions et des mesures expérimentales. Nous n'avons donc pas le droit de nous poser la question suivante : et comment a-t-elle été obtenue ? Cette équation est le produit du raisonnement d'une tête bien éclairée, celle d'Erwin Schrödinger, et sa proposition lui a bien valu le prix Nobel !

Quand on s'intéresse à des distributions indépendantes du temps, il est licite d'écrire :

$$ \Psi = \Psi_x(\mathbf{x}) \cdot \Psi_t(t) \qquad (16) $$

ce qui conduit, après substitution dans l'équation de Schrödinger, à :

$$ \frac{H(\Psi_x)}{\Psi_x} - \frac{T(\Psi_t)}{\Psi_t} = 0 \qquad (17) $$

$$ \text{où } T(\Psi) = -i \frac{\partial \Psi}{\partial t} \qquad (18) $$

et l'opérateur hamiltonien est donné par :

$$ H = -\frac{1}{2m} \nabla^2 + V_{\text{en}} \qquad (19) $$

L'équation (17) implique que :

$$ \left\{ \begin{aligned} -\frac{T(\Psi_t)}{\Psi_t} &= E \\ \frac{H(\Psi_x)}{\Psi_x} &= E \end{aligned} \right. \qquad (20) $$

où $E$ est une constante dont on verra qu'elle correspond à l'énergie. D'où :

$$ \Psi_t(t) = A \exp(-iEt) \qquad (21) $$

et la partie spatiale vérifie le problème aux valeurs propres :

$$ H(\Psi_x(\mathbf{x})) = E \Psi_x(\mathbf{x}) \qquad (22) $$

La solution fondamentale (le « ground state ») est associée à la plus petite valeur propre (énergie). Si l'on note par $\Psi_x^n(\mathbf{x})$ la solution du problème :

$$ H(\Psi_x^n(\mathbf{x})) = E_n \Psi_x^n(\mathbf{x}) \qquad (23) $$

impliquant par ailleurs que :

$$ \Psi_t^n(t) = A \exp(-iE_n t) \qquad (24) $$

il est immédiat de prouver l'indépendance temporelle de la fonction d'onde associée :

$$ |\Psi^n(\mathbf{x}, t)|^2 = |\Psi_t^n(t) \cdot \Psi_x^n(\mathbf{x})|^2 = |\Psi_x^n(\mathbf{x})|^2 \qquad (25) $$

Pour un système multiélectronique, et toujours dans le cadre du modèle de Born-Oppenheimer (dépendance paramétrique de la fonction d'onde avec les coordonnées des noyaux) et dans le cas de la recherche des solutions indépendantes du temps, la fonction d'onde peut s'écrire :
$$ \Psi = \Psi(\mathbf{x}_1, \dots, \mathbf{x}_{N_c}; \mathbf{X}_1, \dots, \mathbf{X}_{N_n}) \qquad (26) $$

où $N_n$ fait référence au nombre de noyaux. La fonction d'onde est solution de l'équation de Schrödinger indépendante du temps :

$$ \begin{aligned} H(\Psi(\mathbf{x}_1, \dots, \mathbf{x}_{N_c}; \mathbf{X}_1, \dots, \mathbf{X}_{N_n})) &= E\Psi(\mathbf{x}_1, \dots, \mathbf{x}_{N_c}; \mathbf{X}_1, \dots, \mathbf{X}_{N_n}) \end{aligned} \qquad (27) $$

où le hamiltonien s'écrit :

$$ H(\Psi) = - \sum_{e=1}^{e=N_c} \frac{\nabla_e^2 \Psi}{2m} + \sum_{e=1}^{e=N_p} \sum_{n=1}^{n=N_n} V_{en} \Psi + \sum_{e=1}^{e=N_c} \sum_{e'=e+1}^{e'=N_c} V_{ee'} \Psi \qquad (28) $$

avec les potentiels coulombiens électron-électron et électron-noyau donnés par :

$$ V_{ee'} = \frac{1}{\|\mathbf{x}_e - \mathbf{x}_{e'}\|} \qquad (29) $$

$$ V_{en} = - \frac{Z}{\|\mathbf{x}_e - \mathbf{X}_n\|} \qquad (30) $$

La principale difficulté associée à la résolution de l'équation de Schrödinger pour un système multiélectronique est son caractère multidimensionnel. Ainsi, on aura :

$$ \Psi = \Psi(\mathbf{x}_1, \dots, \mathbf{x}_{N_c}; \mathbf{X}_1, \dots, \mathbf{X}_{N_n}) : \mathbb{R}^{3 \times N_c} \to \mathbb{C} \qquad (31) $$

Quand on veut calculer la solution d'un problème et qu'on ne peut pas la déterminer analytiquement, on a toujours recours à une technique de discrétisation appropriée. Ainsi, devant l'impossibilité de calculer la solution en tout point (il y en a toute une infinité !), on se limite à calculer une solution approchée pour un certain nombre de points (finis et souvent relativement réduits). Une fois que la solution est obtenue en ces points, elle pourra s'étendre à tout point via l'emploi d'une interpolation adéquate. Jusqu'ici, tout semble simple, mais la difficulté frappe à la porte ! Si nous sommes en dimension 1 et si nous imaginons que la solution en tout point peut s'écrire à partir de la solution connue dans $N$ points, dans le cas bidimensionnel, nous aurons besoin de calculer la solution dans $N \times N$ points, et dans le cas d'un modèle défini dans un espace de dimension $D$, nous aurons besoin de la calculer dans $N^D$ points. Dans le cas de l'équation d'onde multiélectronique (31) avec une trentaine d'électrons (nombre ridicule quand on le compare au nombre d'électrons intervenant dans la plupart des molécules) et un maillage très grossier $N \approx 10$, on aura besoin de calculer la solution dans

$10^{90}$ points, à partir desquels la solution pourrait être définie via une interpolation dans tous les points du domaine $\Omega$ où la fonction d'onde conjointe des trente électrons est définie, $\Omega \subset \mathbb{R}^{90}$. Cependant, cette résolution est aujourd'hui et restera à jamais inabordable parce que l'on estime que $10^{80}$ est le nombre de particules élémentaires dans notre univers ! Donc ce n'est pas la peine de chercher un ordinateur plus puissant !

Cette difficulté est connue en mathématiques sous le nom de malédiction de la dimensionnalité. A ce propos, je me permets d'inclure un paragraphe extrait d'une conférence de R.B. Laughlin (prix Nobel de physique en 1998) à l'Académie de Sciences des Etats-Unis en 2000 [2] :

> « ... the theory of Everything is not even remotely a theory of every thing. We know **this equation is correct** because it has been solved accurately for small number of particles –isolated atoms and small molecules- and **found to agree in minute detail with experiments**. However it cannot be solved accurately when the number of particles exceeds about 10. No computer existing, or that will ever exist, can break this barrier because it is a catastrophe of dimension. If the amount of computer memory required to represent the quantum wavefunction of one particle is $N$ then the amount required to represent the wavefunction of $k$ particles is $N^k$. It is possible to perform approximate calculations ... But the schemes for approximating are not first-principles deductions but rather art keyed to experiment, ...»

Cet extrait prouve la difficulté de résoudre une telle équation et donc d'accéder à la vraie distribution électronique qui devrait, à son tour, nous permettre une compréhension de la structure et de la mécanique des matériaux à l'échelle la plus fine.

Laughlin porte la responsabilité de la difficulté de réussir une telle résolution, même pour un système composé d'une dizaine d'électrons, à la dimensionnalité du modèle. Comme nous l'avons vu, la complexité (taille) du système discret évolue exponentiellement avec la dimension de l'espace.

## 2.3 Résolution de modèles multidimensionnels

Nous avons proposé récemment une technique capable de « by-passer », dans le cas de nombreux modèles, la malédiction de la dimensionnalité [3]. Nous allons supposer un modèle défini dans le domaine $\Omega = (-L, L)^D = \mathbb{P}^D$ où l'on a distribué $N$ nœuds dans chacune des $D$ directions de l'espace. Ainsi, si $N \approx 1000$ et $D \approx 100$, nous aurions eu
besoin d'avoir recours à une technique de discrétisation standard (éléments finis, différences finies, volumes finis, ...), de $N^D \approx 1000^{100} = 10^{300}$ dégrées de liberté, chiffre astronomique au delà de l'imaginable !

Les techniques de Monte-Carlo ont permis de s'affranchir légèrement de la malédiction de la dimensionnalité mais le bruit statistique inhérent à ces techniques n'a pas permis une simulation précise des modèles hautement multidimensionnels. Du côté des techniques déterministes, la méthode connue comme « sparse grids » [4] a rendu possible la résolution de modèles de dimension élevée. Mais certains spécialistes [5] ont estimé à une vingtaine au maximum les dimensions des modèles capables d'être traités par la méthode « sparse grids ».

Nous allons décrire la technique que nous avons proposée et que nous avons appliquée avec succès à la résolution de modèles issus de la chimie quantique mais aussi à la résolution des modèles cinétiques des polymères, comme nous le verrons ultérieurement.

Soit $u(\mathbf{x}) = u(x_1, x_2, \dots, x_D)$ la solution d'une équation aux dérivées partielles. Une approximation de cette fonction pourrait s'écrire sous la forme :

$$u(\mathbf{x}) \approx \sum_{i=1}^{i=Q} F_1^i(x_1) \times \dots \times F_D^i(x_D) = \sum_{i=1}^{i=Q} \prod_{k=1}^{k=D} F_k^i(x_k) \quad (32)$$

Malgré la multi-dimensionnalité du problème, l'approximation (32) entraîne $Q \times N \times D$ inconnues, à la place des $N^D$ nécessaires dans le cadre d'autres techniques de discrétisations habituelles.

La question qui se pose immédiatement, et de façon naturelle, est comment calculer toutes ces fonctions que l'approximation (32) entraîne.

Pour décrire ce calcul, nous allons considérer le problème à résoudre sous la forme :

$$\mathfrak{I}(u) = \mathfrak{A} \quad (33)$$

dont la forme faible s'écrit :

$$\int_{\mathbb{J}^D} u^*(\mathfrak{I}(u) - \mathfrak{A}) dx_1 \cdots dx_D = 0 \quad (34)$$

Maintenant, si l'on suppose qu'à l'itération $n$, les $n$ premiers produits de fonctions ont été calculés :

$$u^n(\mathbf{x}) \approx \sum_{i=1}^{i=n} F_1^i(x_1) \times \dots \times F_D^i(x_D) \quad (35)$$

alors, à cette itération, on cherche à calculer le nouveau produit de fonctions :

$$u^{n+1}(\mathbf{x}) \approx u^n(\mathbf{x}) + F_1^{n+1}(x_1) \times \dots \times F_D^{n+1}(x_D) \quad (36)$$

afin de vérifier au mieux la forme faible (34). La fonction test $u^*$ dans (34) est exprimée par :

$$\begin{aligned} u^*(\mathbf{x}) = & F_1^*(x_1) \times \dots \times F_D^{n+1}(x_D) + \dots \\ & + F_1^{n+1}(x_1) \times \dots \times F_D^*(x_D) \end{aligned} \quad (37)$$

Le problème étant non linéaire, la résolution ne peut se faire qu'à l'aide d'un schéma de linéarisation itératif. Si l'on considère une stratégie de directions alternées, chaque fonction $F_k^{n+1}(x_k)$ est alors obtenue à l'itération $r$ du solveur non-linéaire, $[F_k^{n+1}(x_k)]^{(r)}$, en supposant toutes les autres fonctions connues à l'itération précédente $r-1$. Ainsi, la procédure d'enrichissement peut être facilement implémentée dans des plateformes de calcul parallèle.

Pour le calcul de $[F_k^{n+1}(x_k)]^{(r)}$, la fonction test qui en résulte est :

$$\begin{aligned} u^*(\mathbf{x}) = & [F_1^{n+1}(x_1)]^{(r-1)} \times \dots \times [F_{k-1}^{n+1}(x_{k-1})]^{(r-1)} \times \\ & F_k^*(x_k) \times [F_{k+1}^{n+1}(x_{k+1})]^{(r-1)} \times \dots \times [F_D^{n+1}(x_D)]^{(r-1)} \end{aligned} \quad (38)$$

avec

$$\begin{aligned} u^{n+1}(\mathbf{x}) &\approx u^n(\mathbf{x}) + \\ & [F_1^{n+1}(x_1)]^{(r-1)} \times \dots \times [F_{k-1}^{n+1}(x_{k-1})]^{(r-1)} \times \\ & [F_k^{n+1}(x_k)]^{(r)} \times [F_{k+1}^{n+1}(x_{k+1})]^{(r-1)} \times \dots \times [F_D^{n+1}(x_D)]^{(r-1)} \end{aligned} \quad (39)$$

En introduisant les équations (38) et (39) dans la forme faible (34) et en intégrant par rapport à $x_1, \dots, x_{k-1}, x_{k+1}, \dots, x_D$, on obtient finalement un problème unidimensionnel permettant de calculer

$$[F_k^{n+1}(x_k)]^{(r)}.$$

$$\text{Quand } \max_{k \in \{1, \dots, D\}} \left\| [F_k^{n+1}(x_k)]^{(r)} - [F_k^{n+1}(x_k)]^{(r-1)} \right\| < \varepsilon, \quad (40)$$

on peut arrêter le solveur non-linéaire dédié à l'enrichissement, en écrivant :

$$F_k^{n+1}(x_k) = [F_k^{n+1}(x_k)]^{(r)}; \quad k=1, \dots, D \quad (41)$$
Finalement, le critère d'arrêt pour l'algorithme global de construction de la solution approximée du problème (33) est donné par :

$$ \|u^{n+1}(\mathbf{x}) - u^n(\mathbf{x})\| < \varepsilon \qquad (42) $$

qui permet finalement d'écrire :

$$ u(\mathbf{x}) \approx u^{n+1}(\mathbf{x}) \qquad (43) $$

Une technique similaire fut utilisée par P. Ladeveze au sein de la méthode LATIN pour séparer les fonctions spatio-temporelles. Il appela cette approximation « approximation radiale » [6].

## 2.4 Résolvons l'équation de Schrödinger !

Nous avons appliqué la technique que nous venons de décrire à la résolution de l'équation de Schrödinger, introduite dans la Section 2.2., pour un système composé de 1, 2, 3, 5 et 10 électrons. La Figure 2 représente la distribution électronique. Celle-ci est calculée en appliquant les expressions (12) et (13) à l'état fondamental (« ground state ») de la fonction d'onde solution de l'équation de Schrödinger indépendante du temps (27).

<figure><img src="image_4.png" /><figcaption>Figure 2. Distribution électronique d'un système composé (courbes du bas vers le haut) de 1, 2, 3, 5 et 10 électrons.</figcaption></figure>

Sur la Figure 2, nous pouvons constater que la densité électronique associée à $N_e$ électrons n'est que $N_e$ fois celle qui correspond à un seul électron, c'est-à-dire :

$$ \rho(\mathbf{x}; N_e) = N_e \cdot \rho(\mathbf{x}; N_e = 1) \qquad (44) $$

ce qui implique que tous les électrons se trouvent sur la même orbitale. La forme de cette orbitale est sphérique (Fig. 2), centrée autour du noyau ($\mathbf{x} = \mathbf{0}$) et avec une densité électronique décroissante avec la distance par rapport au noyau.

Nous savons que cette situation est interdite à cause du principe d'exclusion de Pauli, qui affirme que,

pour les fermions (tels que les électrons), on ne peut pas trouver deux particules élémentaires avec des nombres quantiques identiques sur la même orbitale. Cependant, nous voyons déjà que l'équation de Schrödinger que nous avons présentée précédemment n'incorpore pas le nombre quantique dit de spin, et rien dans cette équation n'interdit que tous les électrons se disposent sur l'orbitale d'énergie minimale. Bien entendu l'équation de Schrödinger a précédé dans le temps le principe d'exclusion de Pauli. Ainsi, il a fallu introduire de façon ad-hoc ce principe d'exclusion dans la résolution de l'équation de Schrödinger, comme nous verrons plus loin.

Dans nos travaux, nous avons voulu résoudre l'équation de Schrödinger, sans se soucier de « ce petit caprice dû à Pauli », pour une population assez large d'électrons. Dans [7-8], nous avons présenté des résultats même pour un système composé de 50 électrons et aucune difficulté numérique n'a été recensée. Bien entendu, nous sommes allés plus loin que la limite que le Professeur Laughlin avait prédite. Pourtant, nous sommes en mesure d'affirmer que la vraie difficulté dans la résolution de l'équation de Schrödinger pour des systèmes moléculaires n'est pas dans son caractère multidimensionnel. Celui-ci peut être détourné si, à la place d'une technique de discrétisation standard, on utilise une technique basée sur la séparation de variables, comme celle que nous avons introduite dans la Section 2.3 et qui, par ailleurs, avec une formulation légèrement différente, a été amplement utilisée en chimie quantique, étant à la base des approches dites de Hartree-Fock ou encore de post-Hartree-Fock [9-11].

Cependant, malgré cette touche d'optimisme, il ne faut pas oublier que le système que nous avons résolu n'a aucun sens physique puisqu'il considère les électrons comme des bosons (c'est-à-dire des particules qui n'obéissent pas au principe d'exclusion de Pauli) et ceci est simplement faux. Mais nous avons volontairement voulu séparer les deux difficultés, celles provenant de la multidimensionnalité de celles dues à l'imposition du principe d'exclusion dont nous allons nous occuper dans la prochaine section.

## 2.5 Le principe d'exclusion de Pauli

La vérification de l'exclusion nécessite d'empêcher que deux électrons avec les mêmes nombres quantiques se situent sur la même orbitale. Pour cela, il a fallu inclure dans les coordonnées de chaque électron son spin, c'est-à-dire, dans le reste de cette section $\mathbf{x} = (x, y, z, s)$, où la nouvelle coordonnée $s$ est discrète et prend les valeurs $s = \pm 1/2$.
Maintenant, la fonction d'onde multi-électronique est écrite sous une forme séparée :

$$ \Psi(\mathbf{x}_1, \dots, \mathbf{x}_{N_e}) = \sum_n \Psi^n(\mathbf{x}_1, \dots, \mathbf{x}_{N_e}) \quad (45) $$

avec

$$ \Psi^n(\mathbf{x}_1, \dots, \mathbf{x}_{N_e}) = \Phi_1^n(\mathbf{x}_1) \times \dots \times \Phi_{N_e}^n(\mathbf{x}_{N_e}) \quad (46) $$

Nous devons imposer l'antisymétrie de la fonction d'onde pour garantir le principe d'exclusion. La façon la plus simple (et mathématiquement la plus optimale) de le faire consiste à définir les déterminants de Slater :

$$ D^n(\mathbf{x}_1, \dots, \mathbf{x}_{N_e}) = \begin{vmatrix} \Phi_1^n(\mathbf{x}_1) & \dots & \Phi_1^n(\mathbf{x}_{N_e}) \\ \vdots & \ddots & \vdots \\ \Phi_{N_e}^n(\mathbf{x}_1) & \dots & \Phi_{N_e}^n(\mathbf{x}_{N_e}) \end{vmatrix} \quad (47) $$

Ainsi, la fonction d'onde peut s'exprimer par :

$$ \Psi(\mathbf{x}_1, \dots, \mathbf{x}_{N_e}) = \sum_n D^n(\mathbf{x}_1, \dots, \mathbf{x}_{N_e}) \quad (48) $$

Le déterminant s'annulant si deux lignes sont identiques, l'approximation (48) permet de garantir le principe d'exclusion.

## 2.6 Simulation de quelques molécules simples

Nous nous rendons compte que l'approximation (48) va rapidement mettre en échec toute sorte de technique d'approximation. En effet, le nombre de produits de fonctions résultant du développement du déterminant (47) est $N_e!$, et nous savons bien que le factoriel explose rapidement. Ainsi, si nous considérons 10 électrons dans le système atomique, nous aurons $10! = 3\,628\,800$ produits de fonctions dans chaque terme de la somme (48). Cette valeur est astronomique ! Nous l'avons compris, la difficulté n'est pas dans le caractère multi-dimensionnel de l'équation de Schrödinger, mais dans l'imposition de l'antisymétrie que le principe de Pauli exige. Pour cette raison, nous ne pourrons résoudre que des systèmes contenant au plus une dizaine de « vrais » électrons ! Finalement, les prévisions de Laughlin restent valables !

Nous allons donc nous contenter de résoudre le même système que précédemment avec une population d'électrons évoluant de 1 à 5, ceux-ci « vivant » dans un espace unidimensionnel. La Figure 3 montre les densités électroniques pour les différentes populations électroniques : $N_e = 1, 2, 3, 4, 5$ (courbes du bas vers le haut).

<figure><img src="image_5.png" /><figcaption>Figure 3. Distribution électronique pour un système composé (courbes du bas vers le haut) de 1, 2, 3, 4 et 5 électrons.</figcaption></figure>

Nous constatons que la relation (44) n'est plus vérifiée, certainement à cause des effets du principe d'exclusion. Pour voir plus clair, nous allons représenter $\rho_e - \rho_{e-1}$ pour $e = 2, \dots, 5$ dans les Figures 4 à 7. Nous constatons que le deuxième électron va se placer sur la même orbitale que le premier (ayant des spins contraires, rien ne l'empêche), celle-ci étant sphérique (dite de type s) et centrée autour du noyau. Par contre, le troisième doit chercher une nouvelle orbitale qui ressemble fortement à la forme connue des orbitales p. Le quatrième électron peut se situer dans la même orbitale que le troisième (son spin étant l'opposé) et, finalement, le cinquième devra encore chercher une nouvelle orbitale.

<figure><img src="image_6.png" /><figcaption>Figure 4. $\rho_2(x) - \rho_1(x)$</figcaption></figure>

On pourrait également considérer des petites molécules ne comportant que quelques électrons. Nous avons imaginé deux molécules diatomiques et nous avons calculé la distribution électronique associée au « ground state » pour lequel l'énergie est minimale. Bien entendu, compte tenu de l'utilisation du modèle de Born-Oppenheimer, cette énergie est fonction de la distance entre les deux noyaux.
<figure><img src="image_7.png" /><figcaption>Figure 5. $\rho_3(x) - \rho_2(x)$</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. $\rho_4(x) - \rho_3(x)$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7. $\rho_5(x) - \rho_4(x)$</figcaption></figure>

La Figure 8 montre l'évolution de cette énergie pour deux molécules (une hypothétique, celle du He₂, l'autre réelle, celle du LiH). Dans cette figure, nous constatons une forme du potentiel interatomique bien connue par tout le monde, avec une distance interatomique qui minimise l'énergie, et donc définit la distance interatomique au zéro absolu de température. L'asymétrie du potentiel rend compte du coefficient de dilatation thermique et, finalement,

le gradient du potentiel (pente des courbes) donne les forces agissant sur les noyaux quand on s'éloigne de la position d'équilibre (celle-ci est définie par le minimum de la courbe de potentiel où la pente, et donc la force interatomique, s'annule).

<figure><img src="image_10.png" /><figcaption>Figure 8. Potentiels interatomiques pour deux systèmes diatomiques.</figcaption></figure>

Bien évidemment, tous les résultats que nous venons d'obtenir auraient pu être obtenus dans un espace physique à trois dimensions. Les Figures 9 et 10 illustrent la forme des orbitales attendues et obtenues, respectivement, pour un atome isolé.

<figure><img src="image_11.png" /><figcaption>Figure 9. Orbitales atomiques de type s et p.</figcaption></figure>

Le potentiel interatomique est le point d'entrée des vraies simulations de dynamique moléculaire, comme nous le verrons par la suite. Cependant, la difficulté est désormais évidente : nous ne pouvons pas résoudre la vraie équation de Schrödinger quand le système fait intervenir de nombreux électrons (disons une dizaine !). Des approximations existent : Hartree-Fock, post-Hartree-Fock, DFT – Density Functional Theory, ... [9-11], mais elles ne sont que des approximations, indispensables certes, mais des approximations quand même !
<figure><img src="image_12.png" /><figcaption>Figure 10. Orbitales atomiques calculées de type s et p.</figcaption></figure>

## 2.7 Difficultés et perspectives

Nous venons de souligner les difficultés associées à la résolution de l'équation de Schrödinger, essentiellement à cause de la contrainte d'antisymétrie imposée par le principe d'exclusion de Pauli.

Quelqu'un pourrait se dire que le problème pourrait être détourné en considérant une version plus évoluée de l'équation « fondamentale », comme c'est le cas de l'équation de Dirac [11].

L'équation de Dirac, au contraire de celle de Schrödinger, est en parfait accord avec la relativité restreinte, en particulier avec la transformation de Lorentz. Cette équation est a priori plus appropriée pour décrire les atomes lourds où les effets relativistes jouent un rôle décisif.

L'équation de Dirac fait intervenir le « spinneur » à quatre composantes. Deux pour la particule, l'électron en l'occurrence, et deux pour son antiparticule, le positron. Les deux composantes de chacun prennent en compte les deux valeurs du spin. Malgré toute cette richesse, et contrairement à ce que certains ont affirmé, l'équation de Dirac ne permet pas de se passer de l'antisymétrie. Celle-ci doit encore entre imposée en faisant appel aux déterminants de Slater. Nous voici de retour à la case de départ !

Nous n'avons pu résoudre que certains systèmes atomiques très réduits et, dans ce cas, les résultats

obtenus ont été en parfait accord avec ceux obtenus par résolution de l'équation non-relativiste de Schrödinger.

## 2.8 Des techniques « ab initio » à la dynamique moléculaire

Enfin, si l'on avait pu réussir à résoudre l'équation de Schrödinger, chose impossible pour l'instant comme nous venons de le justifier, on aurait pu avoir accès à la distribution électronique d'un système moléculaire : $\rho(\mathbf{x}; \mathbf{X}_1, \dots, \mathbf{X}_{N_n})$, qui dépendrait paramétriquement des positions des $N_n$ noyaux, notées ci-dessus par $\mathbf{X}_1, \dots, \mathbf{X}_{N_n}$ (nous nous plaçons toujours dans le modèle de Born-Oppenheimer).

Avec la densité électronique obtenue, nous pourrions calculer facilement la force que cette distribution continue de charge négative exerce sur chacun des noyaux qui, grâce au théorème de Hellmann-Feynman [1], n'est qu'une force de type coulombienne. Ainsi, chaque noyau verra l'action répulsive de tous les autres noyaux et l'effet de la distribution électronique négative. Voici la force résultante pour le noyau $k$ :

$$ \mathbf{F}_k = -\sum_{\substack{j=1 \\ j \neq k}}^{j=N_n} \frac{Z_j \cdot Z_k}{|\mathbf{X}_j - \mathbf{X}_k|^3} (\mathbf{X}_j - \mathbf{X}_k) + \int_{\mathbb{R}^3} \frac{Z_k \cdot \rho(\mathbf{x})}{|\mathbf{x} - \mathbf{X}_k|^3} (\mathbf{x} - \mathbf{X}_k) \, d\mathbf{x} \quad (49) $$

Connaissant la force agissant sur chacun des noyaux, on peut calculer leurs accélérations en appliquant la loi fondamentale de la dynamique de Newton :

$$ \mathbf{a}_k = \frac{\mathbf{F}_k}{m}, \quad k = 1, 2, \dots, N_n \quad (50) $$

qui permet à son tour de recalculer leurs vitesses et positions. Le schéma le plus simple consisterait à effectuer :

$$ \left. \begin{array}{l} \mathbf{v}_k \leftarrow \mathbf{v}_k + \mathbf{a}_k \Delta t \\ \mathbf{x}_k \leftarrow \mathbf{x}_k + \mathbf{v}_k \Delta t \end{array} \right\} \quad k = 1, 2, \dots, N_n \quad (51) $$

Le schéma (51), bien qu'étant le plus simple, est le moins précis. Dans les simulations « ab initio », qui procèdent en calculant pour une position donnée des noyaux, (i) la densité électronique (via la résolution de l'équation de Schrödinger ou d'une de ses approximations) ; (ii) la force agissant sur chaque noyau (Eq. (49)) et (iii) les accélérations de tous les noyaux (Eq. (50)) pour finir avec le repositionnement de tous les noyaux et ainsi revenir au point (i), on préfère à la place de (51), utiliser un schéma capable de conserver l'énergie. Ces schémas sont
connus sous le nom de schémas symplectiques, et celui de Verlet est l'un des plus populaires dans le cadre des simulations « ab initio ».

Les simulations « ab initio » sont très gourmandes du point de vue calcul. En plus, on devrait se conformer à la résolution d'une des approximations de l'équation de Schrödinger ou encore de Dirac, la résolution directe n'étant plus possible. Ainsi, des approches plus pragmatiques se sont imposées. Elles sont connues sous le nom de simulations de **dynamique moléculaire**, où la seule différence par rapport à la procédure que nous venons de décrire plus haut, se trouve dans le calcul des forces agissant sur chacun des noyaux. En dynamique moléculaire, on stipule un potentiel (parfois empirique, d'autres fois inspiré des approches quantiques) à deux, trois... noyaux et on calcule la force résultante sur chaque noyau.

La Figure 11 illustre une forme typique d'un potentiel interatomique à deux particules, décrit par le potentiel dit de Lennard-Jones. D'autres alternatives, très nombreuses, existent [12].

<figure><img src="image_13.png" /><figcaption>Figure 11. Potentiel interatomique de Lennard-Jones.</figcaption></figure>

Les paramètres de la loi de Lennard-Jones doivent être identifiés expérimentalement, les prédictions effectuées seront subsidiaires de la précision avec laquelle cette identification sera faite.

### 3. Dynamique moléculaire

Comme nous venons de le voir, la dynamique moléculaire n'est qu'une approche « coarse grained » des simulations « ab initio ».

Une simulation de dynamique moléculaire nécessite de définir les positions et vitesses initiales de toutes les particules (atomes) faisant partie du modèle $X_1^0, \dots, X_{N_n}^0$ et $v_1^0, \dots, v_{N_n}^0$, où l'exposant fait référence au pas de temps considéré.

Il faut également spécifier tous les potentiels agissant entre les différentes particules, puisqu'ils

dépendent de la nature du système (le type de liaison existant entre elles, ...). Avec les positions et les potentiels définis, la force agissant sur chaque noyau peut être calculée, d'où les accélérations, les vitesses et les positions peuvent être recalculées. Bien évidemment, si d'autres forces agissent sur chaque particule, leur prise en compte ne pose pas de problèmes majeurs.

Ce type de simulation est conceptuellement très simple, mais sa validité est assez limitée comme nous le verrons ultérieurement.

Dans les Figures 12 et 13, nous représentons la molécule du poly-(n-propyl isocianate) PPIC dans la modélisation multi-échelle des cristaux liquides polymères [13], avec ses angles de flexion et torsion, respectivement (coordonnées généralisées), sur lesquels portent les potentiels considérés.

<figure><img src="image_14.png" /><figcaption>Figure 12. Molécule de PPIC (chaîne principale et groupe latéral): définition des coordonnées généralisées de flexion dans le plan.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13. Molécule de PPIC (chaîne principale et groupe latéral) : définition des coordonnées généralisées en torsion.</figcaption></figure>

Cette molécule est décrite en dynamique moléculaire par l'hélice représentée à la Figure 14, avant et après la prise en compte des molécules de solvant.

Le principal intérêt de ce type de simulation est l'élucidation des effets corporatifs des molécules des
cristaux liquides. Cependant, la vraie simulation atomistique est encore hors de portée côté calcul. Un « coarse graining » additionnel s'avère nécessaire [13].

<figure><img src="image_16.png" /><figcaption>Figure 14. Représentation DM (dynamique moléculaire) de la molécule de PPIC en absence de solvant (gauche) et recouverte de solvant (droite).</figcaption></figure>

En effet, les simulations atomistiques ont d'importantes limitations. D'un côté, les différents potentiels sont plus ou moins empiriques (dans le meilleur des cas, ils sont inspirés des modélisations quantiques). Ainsi, la validité de ces potentiels est souvent prise en défaut. De l'autre côté, de nombreuses difficultés persistent du point de vue numérique. Le nombre d'atomes explose rapidement avec la dimension du système. A présent, des volumes aux alentours du micron cubique définissent les limites de la capacité de calcul. Une autre difficulté est associée aux pas de temps qu'on doit considérer afin de garantir la stabilité des schémas d'intégration. Souvent ces pas de temps sont de l'ordre de la femto-seconde, impliquant rarement des temps de simulation au delà de la nanoseconde. Dans la plupart des applications le transitoire simulé est de l'ordre de la picoseconde.

D'autres difficultés, plus conceptuelles, persistent également. Parmi elles, on peut citer [8] :

- les difficultés liées à la prise en compte de la température qui induit des fluctuations rapides entraînant des pas de temps extrêmement petits ;

- la difficulté de coupler des discrétisations atomistiques avec des descriptions issues de la mécanique des milieux continus ;

- la difficulté de la prise en compte des conditions aux limites ou encore

- le traitement des régimes transitoires.

## 4. Dynamique brownienne

Comme nous venons de le voir, les simulations atomistiques ne permettent pas de traiter des systèmes au delà de quelques micromètres. Ainsi, de telles simulations trouvent toute leur place dans la simulation des échelles nanométriques et micrométriques que l'on retrouve dans les nanosciences et les nanotechnologies. D'un autre côté, ces simulations permettent de dériver des approches plus grossières (« coarse graining ») à fort contenu physique.

Cependant, on s'intéresse parfois aux échelles au-delà de celles capables d'être traitées avec des modélisations totalement atomistiques. Pour illustrer notre discours, nous allons imaginer un système composé d'une particule en suspension dans un solvant, de masse similaire à celle du solvant. Bien évidemment, une modélisation qui prendrait en compte toutes les particules du solvant et la particule en question est trop riche pour les moyens actuels de calcul, d'où la nécessité de simplifier la modélisation. Une possibilité est de focaliser l'analyse sur la particule à laquelle on s'intéresse, en oubliant l'individualité des particules de solvant qui ne seront retenues dans le modèle que par leurs effets sur la particule qu'on analyse.

Ainsi, l'équation du mouvement de la particule dans le cas 1D (dans le cas 3D, elle est formellement équivalente) s'écrit :

$$m \frac{d^2 x}{dt^2} + \xi \left( V(x) - \frac{dx}{dt} \right) + G = R \quad (52)$$

Elle est connue comme l'équation de Langevin. Les effets du solvant sont présents au travers du terme de friction $\xi(V(x) - dx/dt)$ qui établit la proportionnalité existante entre l'effort de friction et la différence entre la vitesse du fluide $V(x)$ à l'endroit où se trouve la particule et celle de la particule $dx/dt$. L'autre terme dû à la présence du solvant, celui du membre de droite $R$, traduit le bombardement que la particule est en train de subir par les particules du solvant. Cette contribution est aléatoire. Finalement, $G$ représente toutes les autres forces agissant sur la particule.

Dans un milieu isotrope, $R$ aura une moyenne nulle. En ce qui concerne son écart type, nous devrons affiner l'analyse. Pour cela, on définit [1] :

$$B_{\Delta t} = \int_{t=0}^{t=\Delta T} \frac{R(t)}{m} dt \quad (53)$$
En supposant que la particule subit $p$ impacts dans l'intervalle $\Delta t$, cette équation peut se réécrire :

$$B_{\Delta t} = \sum_{i=1}^{i=p} \frac{R(t_i)}{m} \delta t \qquad (54)$$

où $\delta t = \Delta t/p$ et $t_i = (i-1) \delta t$. Le théorème de la limite centrale nous permet de supposer une distribution normale pour $B_{\Delta t}$ :

$$B_{\Delta t} = N(0, q\Delta t) \qquad (55)$$

où $q$ est une valeur à présent inconnue. La dépendance de $\Delta t$ s'explique par le fait que le nombre de variables aléatoires $R(t_i)$ est proportionnel à la taille du pas de temps considéré $\Delta t$.

Pour calculer la valeur de $q$, il suffit d'intégrer l'équation de Langevin pour obtenir la distribution de vitesses à l'équilibre qui devra coïncider avec la distribution de Maxwell-Boltzmann, d'où l'on obtient [1] :

$$q = \frac{2\xi K_B T}{m^2} \qquad (56)$$

où $K_B$ est la constante de Boltzmann.

L'équation (56) reflète la dépendance de la dispersion avec la température du solvant (intensité des impacts) et la dépendance inverse avec la masse de la particule.

Dans la suite, nous allons considérer différents types de systèmes.

## 4.1 Le cas d'une fibre courte au sein d'un écoulement

Dans cette section, nous allons dériver les équations qui gouvernent l'évolution de l'orientation d'une fibre élancée, rigide et de mase négligeable, en suspension dans un écoulement.

La dimension de la fibre étant micro ou nanométrique, le gradient de vitesses du fluide qui l'entoure sera homogène à l'échelle de la fibre. Ainsi, dans la suite, nous allons considérer une fibre de longueur $2L$ dont le centre de gravité se trouve à l'origine du système de coordonnées. La fibre est considérée infiniment élancée et rigide, donc sa conformation est complètement décrite par son orientation : celle-ci est définie par le vecteur unitaire $\mathbf{p}$. Lors de l'évolution, l'origine du système de coordonnées restera solidaire avec le centre de gravité de la fibre. La fibre sera modélisée par deux « perles », où s'appliquent les différentes forces, reliées par un connecteur rigide (Fig. 15).

<figure><img src="image_17.png" /><figcaption>Figure 15. Modélisation d'une fibre au sein d'un écoulement.</figcaption></figure>

Pour l'instant, nous allons négliger les effets browniens qui seront introduits ultérieurement, pour nous focaliser sur l'action de la traîné produite par le fluide. Ainsi, chaque perle est soumise à une force hydrodynamique $\mathbf{F}_H$ :

$$\mathbf{F}_H = \xi (\mathbf{V} - L\dot{\mathbf{p}}) \qquad (57)$$

où $\mathbf{V}$ est la vitesse du fluide dans le point occupé par la « perle », celle-ci (la vitesse du fluide) étant supposée non perturbée par la présence de la fibre, et $L\dot{\mathbf{p}}$ représente la vitesse de la « perle ». L'équation (57) peut s'écrire sous la forme :

$$\mathbf{F}_H = \xi L (\nabla \mathbf{V} \cdot \mathbf{p} - \dot{\mathbf{p}}) \qquad (58)$$

L'équilibre des forces (imposé par l'annulation des forces d'inertie) est automatiquement vérifié puisque, si l'on a $\mathbf{F}_H$ à un bout de la fibre, on aura $-\mathbf{F}_H$ à l'autre bout. Maintenant, avec un moment d'inertie également nul, la résultante des moments devra également s'annuler :

$$\mathbf{M} = 2L\mathbf{p} \times \mathbf{F}_H = \mathbf{0} \qquad (59)$$

Ainsi, la force devra être colinéaire à la fibre, c'est à dire :

$$\mathbf{F}_H = \alpha \mathbf{p} \qquad (60)$$

d'où

$$\xi L (\nabla \mathbf{V} \cdot \mathbf{p} - \dot{\mathbf{p}}) = \alpha \mathbf{p} \qquad (61)$$

Si l'on multiplie les deux termes par $\mathbf{p}$ et si l'on tient compte de l'inextensibilité de la fibre, qui se traduit par $\mathbf{p} \cdot \mathbf{p} = 1$ et donc $\dot{\mathbf{p}} \cdot \mathbf{p} = 0$, on aura :

$$\alpha = \xi L (\mathbf{p} \cdot \nabla \mathbf{V} \cdot \mathbf{p}) \qquad (62)$$

Introduisant l'expression (62) dans (60), on obtient pour la force :

$$\mathbf{F}_H = \xi L (\mathbf{p} \cdot \nabla \mathbf{V} \cdot \mathbf{p}) \mathbf{p} \qquad (63)$$

qui est bien alignée avec la direction de la fibre.
En ce qui concerne la cinématique, il suffit d'introduire (62) dans (61), ce qui conduit à :

$$ \dot{\mathbf{p}} = \nabla \mathbf{V} \cdot \mathbf{p} - (\mathbf{p} \cdot \nabla \mathbf{V} \cdot \mathbf{p}) \mathbf{p} \quad (64) $$

qui n'est rien d'autre que la particularisation de la solution de Jeffery [14] à des fibres infiniment élancées.

## 4.2 Suspension diluée non-brownienne de fibres courtes

Dans le cas d'une suspension diluée de fibres courtes, nous pouvons négliger les interactions entre elles et donc les équations (63) et (64) resteront valables pour chacune d'entre elles.

La passerelle entre l'échelle microscopique (celle où l'évolution de l'orientation des fibres a lieu) et l'échelle macroscopique (où a lieu l'écoulement) peut être définie à travers le calcul du tenseur des contraintes prenant en compte la présence des fibres.

Le tenseur des contraintes peut être calculé en rajoutant à la contribution visqueuse du solvant (celui-ci considéré newtonien pour l'instant), l'extra-contrainte $\tau$ traduisant la contribution des fibres. Cette dernière est calculée avec l'expression du Viriel (connue dans le domaine de la rhéologie sous le nom d'expression de Kramer) :

$$ \sigma = -p\mathbf{I} + 2\eta\mathbf{D} + \tau \quad (65) $$

avec

$$ \tau = \sum_k \mathbf{p}_k \otimes \mathbf{F}_H = \sum_k \lambda \mathbf{D} : (\mathbf{p}_k \otimes \mathbf{p}_k \otimes \mathbf{p}_k \otimes \mathbf{p}_k) \quad (66) $$

où $\mathbf{D} = \frac{\nabla \mathbf{V} + (\nabla \mathbf{V})^T}{2}$, $\mathbf{p}_k$ définit l'orientation de la kième fibre et $\lambda$ est un paramètre qui dépend de la nature du solvant (à travers le coefficient de friction), de la concentration des fibres et de leur longueur.

L'expression (66) nous permet de constater la symétrie du tenseur des contraintes résultant.

## 4.3 Fibre soumise aux effets browniens

Dans cette section, nous allons considérer une fibre isolée, en suspension dans un fluide au repos et soumis au bombardement des molécules de solvant ou aux interactions avec d'autres fibres. Par souci de simplicité, nous allons contraindre l'orientation des fibres dans le plan.

On pourrait sans doute imaginer des forces aléatoires dirigées dans la direction de la fibre (Fig. 16).

<figure><img src="image_18.png" /><figcaption>Figure 16. Fibre soumise à des impacts dirigés dans sa propre direction.</figcaption></figure>

Dans la situation illustrée Figure 16, les forces associées aux nombreux impacts s'annulent mutuellement, et aucun mouvement ni contribution aux contraintes est donc à prévoir.

D'autres impacts, perpendiculaires à la fibre, sont également à considérer (un impact en général peut être décomposé dans une composante dans la direction de la fibre et une autre perpendiculaire à celle-ci). Cette nouvelle situation est illustrée Figure 17.

<figure><img src="image_19.png" /><figcaption>Figure 17. Fibre soumise à des impacts dirigés perpendiculairement à son axe. Les flèches représentent les forces hydrodynamiques associées.</figcaption></figure>

On pourrait penser que cette situation va nous conduire également à une compensation de tous les effets, mais cela est loin d'être le cas comme nous allons le voir. Effectivement, quand une fibre subit un impact, vers le haut par exemple, elle tournera d'un petit angle $\theta$ dans le sens des angles positifs. Une force hydrodynamique normale à la fibre l'accompagnera dans son mouvement (Fig. 18) en s'opposant à celui-ci.

<figure><img src="image_20.png" /><figcaption>Figure 18. Fibre soumise à un impact vers le haut. Les flèches représentent la force hydrodynamique qui apparaît.</figcaption></figure>

La force hydrodynamique peut être calculée en imposant l'équilibre des couples, puisque les forces sont supposées auto-équilibrées.

Dans la configuration tournée d'un angle $\theta$ positif (dans celle non tournée, les efforts seront auto-compensés), la contribution aux contraintes est donnée par :

$$ \tau_{+\theta} = \mathbf{p}_k \otimes \mathbf{F}_H = \begin{pmatrix} \beta \sin\theta \cos\theta & -\beta \cos^2\theta \\ \beta \sin^2\theta & -\beta \sin\theta \cos\theta \end{pmatrix} \quad (67) $$
où $\beta$ représente l'amplitude de la force : $\beta = \|\mathbf{F}_H\|$.

Quand l'impact va vers le bas, il en résulte que :

$$ \boldsymbol{\tau}_{-\theta} = \mathbf{p}_k \otimes \mathbf{F}_{lt} = \begin{pmatrix} \beta \sin\theta \cos\theta & \beta \cos^2\theta \\ -\beta \sin^2\theta & -\beta \sin\theta \cos\theta \end{pmatrix} \quad (68) $$

En considérant tous les impacts possibles, on aura donc pour une fibre orientée selon l'axe horizontal $\varphi = 0$ :

$$ \boldsymbol{\tau}(\varphi=0) = \begin{pmatrix} a & 0 \\ 0 & -a \end{pmatrix} \quad (69) $$

qui peut à son tour être tournée pour nous donner l'expression des extra-contraintes d'origine brownienne pour toutes les orientations $\varphi$ :

$$ \boldsymbol{\tau}(\varphi) = a \begin{pmatrix} \cos^2 \varphi & \sin \varphi \cos \varphi \\ \sin \varphi \cos \varphi & \sin^2 \varphi \end{pmatrix} - \frac{1}{2} \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \quad (70) $$

Ce tenseur est également symétrique.

## 4.4 Simulation de l'écoulement d'une suspension semi-diluée.

L'effet des forces browniennes sera donc introduit via un couple aléatoire schématisé sur la Figure 19.

<figure><img src="image_21.png" /><figcaption>Figure 19. Efforts agissant sur une fibre dans une suspension semi-diluée.</figcaption></figure>

Maintenant, l'équilibre des couples s'écrit :

$$ 2L\mathbf{p} \times \mathbf{F}_H + \mathbf{M}_R = \mathbf{0} \Rightarrow 2L\mathbf{p} \times (\mathbf{F}_H + \mathbf{F}_R) = \mathbf{0} \quad (71) $$

qui implique :

$$ \mathbf{F}_H + \mathbf{F}_R = \alpha \mathbf{p} \Rightarrow \xi L (\nabla \mathbf{V} \cdot \mathbf{p} - \dot{\mathbf{p}}) + F_R \mathbf{p}^{\perp} = \alpha \mathbf{p} \quad (72) $$

où $\mathbf{p}^{\perp}$ est le vecteur unitaire dans la direction normale à l'axe de la fibre. En multipliant l'équation (72) par $\mathbf{p}$ et en prenant en compte son inextensibilité comme précédemment, on obtient à nouveau :

$$ \xi L (\mathbf{p} \cdot \nabla \mathbf{V} \cdot \mathbf{p}) = \alpha \quad (73) $$

La force hydrodynamique qui en résulte est :

$$ \mathbf{F}_H = \xi L (\mathbf{p} \cdot \nabla \mathbf{V} \cdot \mathbf{p}) \mathbf{p} - F_R \mathbf{p}^{\perp} \quad (74) $$

La définition des extra-contraintes nécessite à nouveau de supposer une population de fibres suffisamment large et un schéma d'intégration semi-implicite (« mid-step »).

Concernant la cinématique, la vitesse de rotation des fibres s'obtient en combinant les équations (72) et (73) :

$$ \begin{aligned} \dot{\mathbf{p}} &= \nabla \mathbf{V} \cdot \mathbf{p} - (\mathbf{p} \cdot \nabla \mathbf{V} \cdot \mathbf{p} - \dot{\mathbf{p}})\mathbf{p} + \frac{F_R}{\xi L} \mathbf{p}^{\perp} \\ &= [\dot{\mathbf{p}}]_{Jeffery} + \frac{F_R}{\xi L} \mathbf{p}^{\perp} \end{aligned} \quad (75) $$

Maintenant, nous sommes prêts pour donner l'algorithme de simulation d'une suspension semi-diluée de fibres courtes, issu de la dynamique brownienne :

Soient $\mathbf{p}_1, \dots, \mathbf{p}_N$ les vecteurs définissant l'orientation d'une population de $N$ fibres à l'instant $t_n$. En supposant que le tenseur des extra-contraintes $\boldsymbol{\tau}^n$ est connu au temps $t_n$, nous pouvons :

- calculer le champ de vitesses en résolvant les équations de conservation de la quantité de mouvement, de la masse et l'équation constitutive :

$$ \begin{aligned} \nabla \cdot \boldsymbol{\sigma}^{n+1} + \mathbf{G} &= \rho \frac{d\mathbf{V}^{n+1}}{dt} \\ \nabla \cdot \mathbf{V}^{n+1} &= 0 \end{aligned} \quad (76) $$

$$ \boldsymbol{\sigma}^{n+1} = -\mathbf{p}^{n+1} \mathbf{I} + 2\eta \mathbf{D}^{n+1} + \boldsymbol{\tau}^n $$

- actualiser la position des centres de gravité de toutes les fibres, ainsi que leurs orientations :

$$ \left. \begin{aligned} \mathbf{X}_i^{n+1} &= \mathbf{X}_i^n + \mathbf{V}^{n+1} (\mathbf{X}_i^n) \cdot \Delta t \\ \mathbf{p}_i^{n+1} &= \mathbf{p}_i^n + \dot{\mathbf{p}}^{n+1} (\mathbf{p}_i^n) \cdot \Delta t \end{aligned} \right\} i=1, \dots, N \quad (77) $$

- calculer les extra-contraintes :

$$ \begin{aligned} \boldsymbol{\tau}^{n+1} &= \sum_k \mathbf{p}_k^{n+1} \otimes \mathbf{F}_H^{n+1} = \\ & \sum_k \lambda \mathbf{D}^{n+1} : (\mathbf{p}_k^{n+1} \otimes \mathbf{p}_k^{n+1} \otimes \mathbf{p}_k^{n+1} \otimes \mathbf{p}_k^{n+1}) + \\ & \sum_k a \left( (\mathbf{p}_k^{n+1} \otimes \mathbf{p}_k^{n+1}) - \frac{\mathbf{I}}{d} \right) \end{aligned} \quad (78) $$

où $d = 2$ en 2D et $d = 3$ en 3D.
## 5. Théorie cinétique

Malgré la simplicité conceptuelle des simulations de type dynamique brownienne, les temps de calcul continuent à être leur principal handicap. Les outils de la théorie cinétique viennent nous apporter un coup de main. Mais, parfois, nous verrons que le prix à payer est simplement excessif !

La première grande simplification de la théorie cinétique est d'ignorer l'individualité des entités (particules, fibres, molécules, ...) en introduisant à sa place une fonction de distribution de probabilité.

Avant d'introduire cette fonction et l'équation qui gouverne son évolution, nous avons besoin d'introduire la notion de coordonnée configurationnelle. Nous allons le faire sur quelques exemples.

### 5.1. Systèmes composés de particules

Dans des systèmes composés de particules, position, vitesse et temps suffisent pour décrire le système. Ainsi, nous pourrons définir une fonction $Ψ(x, v, t)$ dont $Ψ(x, v, t) \, dx \, dv$ représente la fraction de particules qui au temps $t$ ont une position appartenant à l'intervalle $[x - dx/2, x + dx/2]$ avec des vitesses appartenant à $[v - dv/2, v + dv/2]$.

Maintenant, nous pouvons définir les deux flux :

$$ \left\{ \begin{aligned} \varphi_x &= \frac{d\mathbf{x}}{dt} \cdot \mathbf{\Psi} \\ \varphi_v &= \frac{d\mathbf{v}}{dt} \cdot \mathbf{\Psi} \end{aligned} \right. \qquad (79) $$

et écrire la conservation de $Ψ(\mathbf{x}, \mathbf{v}, t)$ :

$$ \frac{\partial \Psi}{\partial t} + \nabla_x \cdot \varphi_x + \nabla_v \cdot \varphi_v = C \qquad (80) $$

où $\nabla_x = \partial/\partial\mathbf{x}$, $\nabla_v = \partial/\partial\mathbf{v}$ et C est un éventuel terme source, dit de collision.

L'équation (80) peut être réécrite sous la forme :

$$ \frac{\partial \Psi}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\mathbf{v} \cdot \mathbf{\Psi}) + \frac{\partial}{\partial \mathbf{v}} (\mathbf{a} \cdot \mathbf{\Psi}) = C \qquad (81) $$

qui n'est rien d'autre que l'équation de Boltzmann [15]. Il faut souligner que, pour résoudre l'équation (81), il ne faut rien connaître au préalable sur le champ de vitesses puisque celui-ci est l'une des coordonnées du modèle. Par contre, on aura besoin d'une modélisation nous permettant d'avoir accès au champ d'accélérations $\mathbf{a}(\mathbf{x}, \mathbf{v}, t)$. C'est là où la physique intervient. Voyons quelques situations :

### 5.1.1 Gaz composé de particules chargées.

Dans ce cas, l'accélération peut être calculée à partir de la force électrostatique :

$$ \mathbf{a}(\mathbf{x}) = -\frac{\mathbf{F}_{\text{elec}}(\mathbf{x})}{m} \qquad (82) $$

Les équations (81) et (82) définissent le modèle de Vlasof-Poisson-Boltzmann.

### 5.1.2 Solution diluée de particules non-inertielles.

L'équation de conservation se réduit à :

$$ \frac{\partial \Psi}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\mathbf{v} \cdot \mathbf{\Psi}) = 0 \qquad (83) $$

où la détermination des vitesses fera intervenir des argumentations physiques. Par exemple, pour une suspension de particules ou une suspension colloïdale, on pourrait dériver la vitesse des particules à partir de l'équilibre des forces agissant sur chaque particule (Fig. 20).

<figure><img src="image_22.png" /><figcaption>Figure 20. Schéma de l'équilibre des forces.</figcaption></figure>

Bien entendu, en l'absence de forces browniennes ($\mathbf{F}_R = \mathbf{0}$) et de forces externes ($\mathbf{F}_{\text{Ext}} = \mathbf{0}$), la vitesse des particules est celle du fluide: $\mathbf{v} = \mathbf{V}$.

### 5.1.3 Solution diluée de particules inertielles.

Dans ce cas, l'équation de conservation donne :

$$ \frac{\partial \Psi}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\mathbf{v} \cdot \mathbf{\Psi}) + \frac{\partial}{\partial \mathbf{v}} (\mathbf{a} \cdot \mathbf{\Psi}) = 0 \qquad (84) $$

Comme nous l'avons signalé précédemment, aucune connaissance n'est requise concernant les vitesses, par contre le champ d'accélérations doit être défini. Il peut être déduit, comme précédemment, de l'équilibre des forces qui cette fois devra intégrer l'inertie (Fig. 21).

<figure><img src="image_23.png" /><figcaption>Figure 21. Schéma de l'équilibre des forces.</figcaption></figure>
L'équilibre nous permettra donc de calculer $\mathbf{a}(\mathbf{x}, \mathbf{v}, t)$, permettant la résolution de l'équation (84).

## 5.2. Suspensions de bâtonnets inextensibles.

Cette situation est spécialement intéressante parce qu'elle englobe les suspensions de nanotubes de carbone, nanofibres, fibres courtes, molécules « linéaires » rigides (cristaux liquides), ...

Dans ce cas, la configuration du système est définie par la position du centre de gravité des particules, le temps et l'orientation des fibres (en général les effets d'inertie peuvent être ignorés). Ainsi, la fonction de distribution $\Psi(\mathbf{x}, \mathbf{p}, t)$ nous donne la fraction de fibres qui sont, au temps $t$ et au point $\mathbf{x}$ de l'espace physique, alignées dans la direction $\mathbf{p}$.

L'équation d'évolution associée, connue sous le nom d'équation de Fokker-Planck, s'écrit :

$$ \frac{\partial \Psi}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\mathbf{v} \cdot \Psi) + \frac{\partial}{\partial \mathbf{p}} (\dot{\mathbf{p}} \cdot \Psi) = 0 \quad (85) $$

où la vitesse $\mathbf{v}$ est celle du fluide, tandis que la vitesse de rotation résulte de l'équilibre des couples que nous avons formulé dans la section précédente et l'équation (64). La composante brownienne de la vitesse de rotation des fibres est souvent exprimée à partir du gradient de la fonction de distribution elle-même, ce qui donne comme résultat un terme diffusif dans l'équation de Fokker-Planck.

## 5.3. Polymères en solution.

Une molécule au sein d'un solvant peut être modélisée par une série d'haltères élastiques, comme nous l'illustrons sur la Figure 22.

<figure><img src="image_24.png" /><figcaption>Figure 22. Représentation d'une molécule en solution.</figcaption></figure>

L'action du solvant sur la molécule est supposée concentrée dans les $N + 1$ « perles ». Les ressorts peuvent être linéaires ou non linéaires (FENE). Comme d'habitude, l'équilibre des forces au niveau

de chaque « perle » permet d'accéder à sa vitesse et, par conséquent, au calcul des vitesses d'extension et de rotation de chaque connecteur élastique.

Cette vitesse résulte des effets de la traîné imposée par l'écoulement, de la contribution élastique induite par la raideur des ressorts élastiques et finalement des actions browniennes (voir le magnifique ouvrage de Bird et al. [16] pour plus de détails):

$$ \dot{\mathbf{q}}_i - \nabla \mathbf{v} \mathbf{q}_i + \Xi \mathbf{Q} - \frac{D}{\Psi} \frac{\partial \Psi}{\partial \mathbf{q}_i} \quad (86) $$

où $\Xi$ représente l'élasticité de la molécule et $\mathbf{Q}^T = (\mathbf{q}_1, \dots, \mathbf{q}_N)$.

Ainsi, une représentation adéquate de la conformation moléculaire peut se faire à l'aide de la fonction de distribution $\Psi(\mathbf{x}, \mathbf{q}_1, \dots, \mathbf{q}_N, t)$. Cette distribution peut être calculée en résolvant l'équation de Fokker-Planck associée, qui prend la forme :

$$ \frac{\partial \Psi}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\mathbf{v} \cdot \Psi) + \sum_{i=1}^{N} \frac{\partial}{\partial \mathbf{q}_i} (\dot{\mathbf{q}}_i \cdot \Psi) = 0 \quad (87) $$

avec la vitesse de chaque connecteur $\dot{\mathbf{q}}_i$ donnée par l'équation (86).

## 5.4. Autres modèles cinétiques.

Le même type de procédure peut être suivi dans de nombreux domaines dont nous allons citer quelques exemples.

### 5.4.1 Les polymères fondus

Souvent, les modèles sont basés sur la théorie de la reptation de De Gennes [17] et Doi et Edwards [18]. Des modèles plus évolués ont été proposés et définis dans le contexte de la théorie cinétique, incorporant la double reptation, la « convective contraint release », ... [19, 20].

<figure><img src="image_25.png" /><figcaption>Figure 23. Illustration du concept de tube.</figcaption></figure>

Ainsi, l'un des modèles pionniers, celui de Doi-Edwards, introduit la notion de tube défini à partir
des restrictions cinématiques latérales que les molécules environnantes exercent sur la molécule considérée. Ainsi, une molécule ne peut que bouger dans son tube, comme le fait un verre de terre ! (Fig. 23).

La tête de la molécule a « a priori » la liberté de prendre un chemin ou un autre, et cela permet à la molécule d'abandonner progressivement son tube original (la situation est schématisée Figure 24).

<figure><img src="image_26.png" /><figcaption>Figure 24. Illustration du mouvement de reptation.</figcaption></figure>

Une façon simple pour décrire la conformation consiste à donner le vecteur tangent **u** à chaque abscisse curviligne **s** (Fig. 25).

<figure><img src="image_27.png" /><figcaption>Figure 25. Coordonnées conformationnelles d'une molécule dans un système enchevêtré.</figcaption></figure>

La fonction de distribution est donnée par $\Psi(\mathbf{x},\mathbf{u},s,t)$. Dans le modèle de Doi-Edwards, nous pouvons définir les deux flux associés aux coordonnées conformationnelles :

$$ egin{cases} \varphi_u = \frac{du}{dt} \cdot \Psi \\ \varphi_s = \frac{ds}{dt} \cdot \Psi \end{cases} \qquad (88) $$

où $du/dt$ est donné par la déformation affine, tout en tenant compte de la norme unitaire du vecteur **u** qui se traduit par :

$$ \frac{du}{dt} = \nabla V \cdot \mathbf{u} - (\mathbf{u} \cdot \nabla V \cdot \mathbf{u}) \mathbf{u} \qquad (89) $$

qui coïncide avec l'équation de Jeffery pour des fibres infiniment élancées, que nous avons déduite précédemment (Eq. (64)).

Concernant le terme $ds/dt$, il est supposé de type diffusif et il prend donc la forme habituelle :

$$ \frac{ds}{dt} = - \frac{D}{\Psi} \frac{\partial \Psi}{\partial s} \qquad (90) $$

Avec ces expressions des flux, nous pouvons établir l'équation de conservation de $\Psi(\mathbf{x},\mathbf{u},s,t)$, c'est-à-dire l'équation de Fokker-Planck du modèle de Doi-Edwards :

$$ \frac{\partial \Psi}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\mathbf{V} \cdot \Psi) + \frac{\partial}{\partial \mathbf{u}} (\dot{\mathbf{u}} \cdot \Psi) - D \frac{\partial^2 \Psi}{\partial s^2} = 0 \qquad (91) $$

Des modèles plus sophistiqués existent [21].

### 5.4.2 Théorie cinétique microstructurale du mélange

Elle prend en compte l'existence d'une micro-structure représentée par une longueur caractéristique, une orientation et une morphologie [22]. Pour cela, on part d'une situation comme celle illustrée sur la Figure 26. On définit $\Psi(\mathbf{x}, \mathbf{n}, t)$ qui représente l'aire de l'interface par unité de volume pour une position $\mathbf{x}$, au temps $t$ et orientée selon la direction $\mathbf{n}$ [23].

<figure><img src="image_28.png" /><figcaption>Figure 26. Microstructure lors du mélange passif.</figcaption></figure>

Sans tenir compte des effets de coalescence ni de la tension de surface, l'interface évolue à cause du changement d'orientation et également à cause de l'étirement (« stretching ») qu'elle subit. Dans le cadre du modèle de Chella et Ottino [24], l'équation cinétique qui en résulte est :

$$ \frac{\partial \Psi}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\mathbf{V} \cdot \Psi) + \frac{\partial}{\partial \mathbf{n}} (\dot{\mathbf{n}} \cdot \Psi) = C \Psi \qquad (92) $$

où l'évolution du vecteur normal à la surface prend une fois de plus la forme :

$$ \frac{d\mathbf{n}}{dt} = -\nabla V \cdot \mathbf{n} + (\mathbf{n} \cdot \nabla V \cdot \mathbf{n}) \mathbf{n} \qquad (93) $$

et le terme $C\Psi$ est introduit pour prendre en compte le « stretching » de l'interface :

$$ C = (\mathbf{n} \cdot \nabla V \cdot \mathbf{n}) \qquad (94) $$

Cette modélisation peur être enrichie pour prendre en compte d'autres phénomènes physiques comme celui de la coalescence ou de la tension de surface. Nous avons prouvé [23] que le modèle de Doi-Ohta
[25] peut être retrouvé en ajoutant un terme de diffusion et un terme source dans l'équation (92).

### 5.4.3 Description cinétique de la rhéologie des suspensions de nanotubes de carbone

Ici, l'orientation mais également l'état d'agrégation du système doivent être prises en compte [26-32].

Dans les suspensions de nanotubes de carbone, il faut distinguer le cas où les nanotubes sont fonctionnalisés pour éviter leur agrégation et le cas où ils ne le sont pas. Cette dernière situation est susceptible d'entraîner leur agrégation avec d'importantes conséquences côté rhéologique.

Dans le cas des nanotubes fonctionnalisés, le modèle cinétique développé pour les fibres courtes (Eqs. (75) et (85)) s'est avéré pertinent, permettant de décrire la rhéologie non linéaire [32] comme nous pouvons le constater sur la Figure 27.

<figure><img src="image_29.png" /><figcaption>Figure 27. Suspension de nanotubes de carbone fonctionnalisés : rhéologie non-linéaire.</figcaption></figure>

Dans le calcul de la viscosité apparente, il nous a fallu définir la contrainte associée à une population de tubes dans le cadre de la théorie cinétique. Pour cela, il suffit de reprendre l'expression que nous avons déduite dans le contexte de la dynamique brownienne (équation (78)), que nous reprenons ici par souci de clarté :

$$ \tau = \sum_k \lambda \mathbf{D} : (\mathbf{p}_k \otimes \mathbf{p}_k \otimes \mathbf{p}_k \otimes \mathbf{p}_k) + \sum_k a \left( (\mathbf{p}_k \otimes \mathbf{p}_k) - \frac{\mathbf{I}}{d} \right) $$

Maintenant, pour passer d'une approche brownienne, qui procède avec la prise en compte de toutes les entités individuelles, à une description cinétique, où nous n'avons qu'une fonction de distribution, il n'y a qu'à substituer les sommes par des intégrales :

$$ \tau = \lambda \int_{S(0,1)} \mathbf{D} : (\mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p}) \Psi(\mathbf{p}) d\mathbf{p} + a \left( \left( \int_{S(0,1)} (\mathbf{p} \otimes \mathbf{p}) \Psi(\mathbf{p}) d\mathbf{p} \right) - \frac{\mathbf{I}}{3} \right) \quad (95) $$

avec $S(0, 1)$ la surface de la sphère unité.

En introduisant les tenseurs d'orientation de quatrième et deuxième ordre :

$$ \left\{ \begin{aligned} ^4\mathbf{a} &= \int_{S(0,1)} (\mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p}) \Psi(\mathbf{p}) \, d\mathbf{p} \\ ^2\mathbf{a} &= \int_{S(0,1)} (\mathbf{p} \otimes \mathbf{p}) \Psi(\mathbf{p}) \, d\mathbf{p} \end{aligned} \right. \quad (96) $$

et en définissant $\lambda = 2\eta N_p$ et $a = D$, on obtient :

$$ \tau = 2\eta N_p \mathbf{D} : ^4\mathbf{a} + D \left( ^2\mathbf{a} - \frac{\mathbf{I}}{3} \right) \quad (97) $$

La viscosité apparente est définie (dans un écoulement de cisaillement simple) à partir de :

$$ \eta_a = \frac{\sigma_{12}}{\dot{\gamma}} = \eta \left( 1 + 2N_p a_{1212} \right) \quad (98) $$

Le modèle a été enrichi afin de reproduire les résultats obtenus en régime linéaire [32]. Pour cela, il nous a suffit de conjecturer que la cinématique induite par des oscillations harmoniques de faible amplitude représente une perturbation de l'état d'équilibre isotrope induit par la fonctionnalisation des nanotubes qui elle-même induit un réseau d'interactions. Si la raideur de ce réseau est proportionnelle à la racine carrée de la fréquence de l'oscillation appliquée à l'échantillon, on retrouve l'évolution des modules $G'$ et $G''$ mesurés expérimentalement [32].

Quand les nanotubes ne sont pas fonctionnalisés, on observe une augmentation de la viscosité apparente (Fig. 28) qui ne peut s'expliquer que par l'existence d'une structure interne qui évolue dans le temps pour atteindre un état stationnaire dont la valeur est fonction des taux de déformation imposés.

De plus, cet état d'agrégation, ainsi que son évolution avec les taux de déformation appliqués, a pu être observé au microscope [31].

Il a donc fallu enrichir le modèle cinétique qui ne prend en compte que l'orientation des nanotubes (85) pour y introduire les effets d'une cinétique d'agrégation-désagrégation [31].
<figure><img src="image_30.png" /><figcaption>Figure 28. Comparaison entre la rhéologie des suspensions de nanotubes de carbone fonctionnalisés et non fonctionnalisés.</figcaption></figure>

La configuration est, dans cette situation, décrite par la fonction de distribution $\Psi(x, \mathbf{p}, n, t)$ qui représente la fraction de tubes au point $x$ et au temps $t$ appartenant à une population $n \in [0,1]$ ($n = 0$ pour les tubes libres de contacts avec les autres et $n = 1$ pour l'état d'enchêvêtrement maximal) et qui sont orientés selon $\mathbf{p}$. Dans [31], nous avons établi l'équation de Fokker-Planck suivante pour un écoulement homogène stationnaire :

$$ \frac{\partial \Psi(\mathbf{p}, n)}{\partial t} = \left. \frac{\partial \Psi(\mathbf{p}, n)}{\partial t} \right|_{\text{Orient}} + \left. \frac{\partial \Psi(\mathbf{p}, n)}{\partial t} \right|_{\text{Agrég}} \quad (99) $$

avec:

$$ \left. \frac{\partial \Psi(\mathbf{p}, n)}{\partial t} \right|_{\text{Orient}} = - \left. \frac{\partial}{\partial \mathbf{p}} \left( \dot{\mathbf{p}} \cdot \Psi(\mathbf{p}, n) \right) \right| \quad (100) $$

et

$$ \left. \frac{\partial \Psi(\mathbf{p}, n)}{\partial t} \right|_{\text{Agrég}} = v_c \int_{0}^{n} \Psi(\mathbf{p}, r) \, dr - v_c \Psi(\mathbf{p}, n) \int_{n}^{1} dr \\ v_d \int_{n}^{1} \Psi(\mathbf{p}, r) \, dr - v_d \Psi(\mathbf{p}, n) \int_{0}^{n} dr \quad (101) $$

Ce modèle ne comporte que deux paramètres rhéologiques additionnels : la vitesse d'agrégation $v_c$ et celle de désagrégation $v_d$. La richesse du modèle vient de la possibilité de faire dépendre les différents paramètres rhéologiques (diffusion brownienne et $N_p$) avec la population considérée (état de la microstructure) et en même temps de faire dépendre les vitesses d'agrégation et de désagrégation des taux de déformation. Une

modélisation extrêmement simple nous a permis un recalage très fin des mesures expérimentales, comme le prouve la Figure 29.

<figure><img src="image_31.png" /><figcaption>Figure 29. Suspension de nanotubes de carbone non-fonctionnalisés : rhéologie non-linéaire.</figcaption></figure>

Différentes questions restent ouvertes. En voici quelques unes :

- la compréhension de la rhéologie linéaire, même pour les nanotubes fonctionnalisés, nécessite probablement de nouveaux ingrédients (flexion des tubes, modélisation des interactions dues à la fonctionnalisation des tubes, ...);

- il faudrait retarder la cinématique de Jeffery pour des tubes appartenant à des agrégats et cela,
bien entendu, sans compromettre l'objectivité des expressions résultantes;

- la proposition d'une cinétique d'agrégation-désagrégation avec un grand contenu physique reste à faire;

- enfin, le couplage entre les nanotubes de carbone et un solvant viscoélastique est aussi une question récurrente et encore ouverte.

## 5.4. Discussion : avantages et difficultés.

Effectivement, l'utilisation d'une description cinétique simplifie énormément la représentation du système. Par contre, comme nous allons le voir, la résolution numérique de tels systèmes est loin d'être simple.

Premièrement, à la place de décrire le système à partir d'un nombre exorbitant de particules (ce qui n'empêchera pas malgré tout le bruit statistique), le système en théorie cinétique est décrit par une (parfois plusieurs) fonction scalaire : la fonction de distribution des configurations du système. La résolution de l'équation de Fokker-Planck, qui n'est qu'une équation différentielle aux dérivées partielles de type convection-diffusion, permet de trouver l'évolution de cette fonction de distribution définissant la conformation du système. Cependant, cette résolution n'est pas sans difficulté numérique.

En effet, si l'on revient par exemple à l'équation (87), nous voyons qu'il s'agit d'une équation définie dans un espace multidimensionnel. Ainsi, en supposant des chaînes composées d'une dizaine de segments, immergées dans un écoulement homogène et stationnaire, le modèle résultant (équation de Fokker-Planck) comporte 10 dimensions dans un espace physique unidimensionnel, 20 dimensions dans un espace physique 2D et 30 dimensions dans un espace physique tridimensionnel.

Cette difficulté a motivé d'une part la simplification des modèles qui, habituellement, ne comportent qu'un seul ressort (associé au vecteur bout-à-bout) et, d'autre part, la recherche d'algorithmes de résolution capables de s'affranchir, du moins en partie, des difficultés liées à la malédiction de la dimensionnalité que nous avons détaillée au début de ce travail.

Sur la première possibilité, nous n'allons pas insister : il existe une vaste bibliographie autour des modèles mono-segment avec un ressort linéaire (qui conduit au modèle macroscopique de Maxwell) ou non-linéaire : FENE, FENE-P, ... [33]. En tout cas, et même pour ce type de modèles, il a fallu chercher des techniques capables de procéder à leur résolution avec efficacité.

Une autre possibilité dans la même ligne a été de dériver des modèles macroscopiques. Cependant, on rencontre souvent dans ce passage des problèmes de fermeture qui ont besoin d'un choix parfois assez arbitraire ! Nous reviendrons sur ce point dans la section suivante.

Enfin, pour résoudre les différentes équations de Fokker-Planck que nous avons rencontrées tout au long de cette section, toutes multidimensionnelles, des stratégies avancées de calcul sont nécessaires. Les techniques de discrétisation classiques basées sur l'emploi de maillages subjacents ne peuvent être appliquées que pour une dimensionnalité suffisamment réduite [34].

Historiquement, la première solution proposée a été de revenir de là où nous venons : les simulations browniennes !

## 5.5. Résolutions stochastiques : une solution de dépannage !

On savait depuis longtemps que les équations de convection-diffusion peuvent être résolues avec une technique particulaire stochastique. Pour l'illustrer, considérons l'équation de convection-diffusion suivante :

$$ \frac{\partial C}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\mathbf{V} \cdot C) = D \frac{\partial^2 C}{\partial \mathbf{x}^2} \quad (102) $$

qui peut être écrite sous une forme purement convective :

$$ \frac{\partial C}{\partial t} + \frac{\partial}{\partial \mathbf{x}} \left( \left( \mathbf{V} - \frac{D}{C} \frac{\partial C}{\partial \mathbf{x}} \right) \cdot \mathbf{C} \right) = 0 \quad (103) $$

ou encore :

$$ \frac{\partial C}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\mathbf{V}^* \cdot \mathbf{C}) = 0 \quad (104) $$

avec la vitesse équivalente $\mathbf{V}^*$ définie par

$$ \mathbf{V}^* = \mathbf{V} - \frac{D}{C} \frac{\partial C}{\partial \mathbf{x}} \quad (105) $$

Cette équation (104) a été à la base des techniques particulaires déterministes de résolution SPH (« Smooth Particle Hydrodynamics »). Cependant, ces techniques cachent de nombreux paramètres difficiles à estimer pour avoir des bons résultats [35]. Encore sur la base des techniques particulaires, des stratégies stochastiques peuvent également être appliquées. Celles-ci n'entrainent pas de paramètres ajustables, mais par contre ont besoin de populations de particules (trajectoires du processus stochastique) suffisamment larges [36].
A titre d'illustration, nous allons procéder à la résolution de l'équation (102) par une technique stochastique. Pour cela, nous supposons la condition initiale concentrée dans un certain nombre de points $\mathbf{x}_i^0$, $i=1,2,\dots,N$ :

$$C(\mathbf{x}, t=0) \approx c \cdot \delta(\mathbf{x} - \mathbf{x}_i^0) \qquad (106)$$

où $\delta$) est la masse de Dirac.

Maintenant, on peut appliquer la convection et les effets diffusifs à toutes ces masses ponctuelles afin de suivre l'évolution du système :

$$\begin{cases} \mathbf{x}_i^* = \mathbf{x}_i^n + \mathbf{V}(\mathbf{x}_i^n) \cdot \Delta t \\ \mathbf{x}_i^{n+1} = \mathbf{x}_i^* + N(0, \sqrt{2D\Delta t})\mathbf{n} \end{cases} \qquad (107)$$

où $N(0, \sqrt{2D\Delta t})$ est une variable aléatoire normale de moyenne nulle et d'écart type $\sqrt{2D\Delta t}$. $\mathbf{n}$ est un vecteur unitaire défini sur la surface de la sphère de rayon unitaire :

$$\mathbf{n} = U(S_{B(0,t)}) \qquad (108)$$

Des choix plus performants que le schéma explicite d'intégration du premier ordre (107) existent, mais nous avons préféré utiliser celui-ci parce qu'à ce stade nous préférons la clarté par rapport à l'efficacité !

L'analyse que nous venons d'effectuer peut s'étendre aux différentes équations de Fokker-Planck que nous avons présentées dans la section précédente. Le lecteur intéressé peut consulter le magnifique ouvrage de H. C. Ottinger [37] ainsi que les références incontournables [20, 38].

Le grand avantage des solveurs stochastiques est le fait de faire augmenter la complexité du calcul de façon linéaire avec la dimension des espaces où les modèles sont définis (dans les méthodes déterministes basées sur la construction d'un maillage, la complexité augmente exponentiellement avec la dimension de l'espace de travail). Cette évolution linéaire est possible parce que ce que l'on cherche n'est pas la reconstruction de la fonction de distribution elle-même (qui aurait entrainé une complexité exponentielle) mais plutôt le calcul de certains moments de cette distribution qui définissent la contrainte (voir l'expression (95) à titre d'exemple). Quand il s'agit de calculer des intégrales, un nombre assez réduit de points suffisent, même en dimension élevée: voici le miracle des techniques stochastiques !

Cependant, et malgré sa relative efficacité, deux inconvénients majeurs demeurent :

- la nécessité de résoudre le problème transitoire, même quand on n'est intéressé que par le régime stationnaire;

- l'inévitable bruit statistique qui rend difficile le calcul de sensibilités ou des dérivées de la solution (celles-ci sont nécessaires dans le cas où l'on s'intéresse à l'optimisation ou à l'identification rhéologique inverse). La Figure 30 illustre l'évolution de quelques trajectoires du problème de convection-diffusion résolu dans un espace unidimensionnel.

<figure><img src="image_32.png" /><figcaption>Figure 30. Position (en ordonnées) fonction du temps (en abscisses) pour 7 trajectoires du solveur stochastique de l'équation de convection-diffusion unidimensionnelle.</figcaption></figure>

Certes, le bruit peut être diminué en utilisant des techniques dites de réduction de la variance, mais même si on l'attenue, on ne peut pas l'éliminer. Une alternative aux techniques stochastiques que nous venons de décrire, qui ont été présentées comme des techniques permettant de réduire la variance, est la technique BCF (« Brownian Configuration Fields ») [39]. Cette stratégie stochastique mais eulérienne a permis de la combiner avec des techniques avancées de réduction de modèles [40].

Pour « by-passer » ces difficultés, une résolution déterministe de l'équation aux dérivées partielles aurait été souhaitable, mais il a fallu attendre à cause de la malédiction de la dimensionnalité. Les algorithmes « gloutons » sont venus à notre secours !

## 5.6. Techniques avancées de résolution.

Nous avons présenté dans la Section 2.3 une technique, basée sur la séparation de variables, capable de résoudre des modèles fortement multidimensionnels. Nous allons illustrer dans cette section l'une de ces résolutions. Pour cela, nous allons considérer le modèle d'une macromolécule en solution (schématisée Fig. 23) et composée de 9 segments dans un espace physique unidimensionnel
(encore pour des rasions de clarté). Un écoulement d'élongation homogène est imposé à la suspension et l'on cherche à calculer la distribution stationnaire.

La fonction de distribution s'écrit sous la forme :

$$ \Psi(q_1, \dots, q_9) \approx \sum_{i=1}^{Q} F_1^i(q_1) \times \dots \times F_9^i(q_9) \quad (109) $$

qui permet de résoudre l'équation de Fokker-Planck associée :

$$ \sum_{i=1}^{i=9} \frac{\partial}{\partial q_i} (\dot{q}_i \cdot \Psi) = 0 \quad (110) $$

Pour les détails concernant la résolution, le lecteur peut s'adresser à la Section 2.3 ou bien à la référence [3] où différentes situations plus complexes sont traitées. L'algorithme de résolution converge après 7 enrichissements, c'est-à-dire $Q = 7$ dans l'équation (109). La Figure 31 illustre les 7 fonctions pour les segments aux extrémités de la chaîne $F_1^1(q_1), \dots, F_1^7(q_1)$, tandis que la Figure 32 trace ces 7 fonctions pour le segment central de la chaîne $F_5^1(q_5), \dots, F_5^7(q_5)$.

Une fois que la fonction de distribution a été calculée, les contraintes peuvent être calculées par intégration.

Pour effectuer une analyse transitoire, il suffit de considérer le temps comme une nouvelle coordonnée [41] :

$$ \Psi(q_1, \dots, q_N, t) \approx \sum_{i=1}^{Q} F_1^i(q_1) \times \dots \times F_9^i(q_9) \times G^i(t) \quad (111) $$

Ce type d'analyse transitoire nous permet de tracer les courbes d'évolution de la contrainte, des différences de contraintes normales, ... lors d'un calcul transitoire, au démarrage, par exemple, d'un écoulement rhéométrique.

Cette technique de séparation de variables nous a permis également de traiter les modèles des polymères fondus basés sur la théorie de la reptation [42].

Quand la dimension des modèles est modérée, la construction de la base séparée peut être effectuée dans le cadre d'une simulation standard (éléments finis, ...) par l'introduction de décomposition de Karhunen-Loève (« Proper Othogonal Decomposition ») combinée avec des enrichissements basés sur l'emploi des sous-espaces de Krylov. Cette stratégie a été appliquée avec succès dans la résolution des écoulements rhéometriques de suspensions de fibres

<figure><img src="image_33.png" /><figcaption>Figure 31. Séparation de variables : fonctions</figcaption></figure>

$F_1^1(q_1), \dots, F_1^7(q_1)$

<figure><img src="image_34.png" /><figcaption>Figure 32. Séparation de variables : fonctions</figcaption></figure>

$F_5^1(q_5), \dots, F_5^7(q_5)$

courtes [43], la résolution du modèle FENE monohaltère associé aux polymères en solution [44], aux cristaux liquides-polymères [13], ainsi que pour l'évaluation des relations de fermeture des différentes équations constitutives macroscopiques [45-47].

## 6. Modélisation micro-macro : enfin l'écoulement !

Pour illustrer la démarche de passage à l'échelle macroscopique, nous avons décidé de nous concentrer sur les suspensions de fibres courtes que nous avons développées dans les sections précédentes.

### 6.1. Modèles macroscopiques

La première possibilité consiste à dériver un modèle complètement macroscopique. Le principal avantage d'une modélisation entièrement à l'échelle macroscopique est la réduction des temps de calcul [45, 47]. Le prix à payer, comme nous allons le voir,
est la nécessité d'introduire certaines hypothèses dont la validité est parfois mise en échec [48-49].

Nous allons retracer la construction d'un tel modèle. Le point de départ est l'expression de Jeffery que nous avons obtenue précédemment. Nous allons continuer à supposer les fibres à rapport de forme infini. Egalement, et pour des raisons de simplicité, nous allons restreindre notre analyse aux suspensions diluées pour lesquelles nous pouvons négliger les interactions et donc les termes browniens.

A partir de la définition du tenseur d'orientation de second ordre, on obtient son évolution :

$$ 
\begin{aligned}
\dot{\mathbf{a}} &= \int_{S(0,1)} (\dot{\mathbf{p}} \otimes \mathbf{p}) \Psi(\mathbf{p}) \, d\mathbf{p} + \int_{S(0,1)} (\mathbf{p} \otimes \dot{\mathbf{p}}) \Psi(\mathbf{p}) \, d\mathbf{p} \\
&= \nabla \mathbf{V} \cdot {}^2\mathbf{a} + {}^2\mathbf{a} \cdot (\nabla \mathbf{V})^T - 2 {}^4\mathbf{a} : \mathbf{D}
\end{aligned}
\quad (112) $$

Côté contrainte, comme nous l'avons vu, le tenseur de Cauchy peut s'écrire :

$$ \sigma = -p\mathbf{I} + 2\eta\mathbf{D} + 2\eta N_p {}^4\mathbf{a} : \mathbf{D} \quad (113) $$

Le modèle doit être complété par les équations de conservation de la masse et de la quantité de mouvement où dans cette dernière, on néglige les effets massiques et d'inertie :

$$ \left\{ \begin{array}{l} \nabla \cdot \boldsymbol{\sigma} = 0 \\ \nabla \cdot \mathbf{V} = 0 \end{array} \right. \quad (114) $$

Le modèle (112)-(114) permet d'obtenir la cinématique de l'écoulement ainsi que l'évolution de la microstructure, celle-ci décrite à partir des tenseurs d'orientation. Cependant, la résolution n'est pas encore possible du fait de l'existence de deux tenseurs d'orientations inconnus. Pour parvenir à fermer le modèle, et donc à pouvoir procéder à sa résolution, le tenseur de quatrième ordre doit s'exprimer en fonction de celui de deuxième ordre. Il n'existe aucune procédure exacte pour y parvenir. Différents choix sont possibles, l'un des plus employés est la fermeture quadratique :

$$ {}^4\mathbf{a} = {}^2\mathbf{a} \otimes {}^2\mathbf{a} \quad (115) $$

Enfin, d'autres alternatives existent, mais aucune n'a de validité générale [50-51]. Voici encore une raison justifiant l'emploi des procédures, dites micro-macro, qui peuvent se passer de la nécessité d'introduire des fermetures à faible contenu physique et d'une validité, dans la plus part des cas, très mitigée.

## 6.2. Modèles micro-macro

Pour éviter d'être confrontés au choix d'une fermeture, nous allons coupler une modélisation macroscopique de la cinématique de l'écoulement (Eqs. (113) et (114)) avec une modélisation microscopique de l'évolution de la microstructure. Cette dernière est décrite par l'équation de Fokker-Planck correspondante (85), que nous reprenons ici par souci de clarté :

$$ \frac{\partial \Psi}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\mathbf{V} \cdot \Psi) + \frac{\partial}{\partial \mathbf{p}} (\dot{\mathbf{p}} \cdot \Psi) = 0 \quad (116) $$

La résolution de cette dernière équation nous permet d'accéder à la fonction de distribution $\Psi(\mathbf{x}, \mathbf{p}, t)$, nous permettant de remonter au calcul du tenseur de quatrième ordre :

$$ {}^4\mathbf{a}(\mathbf{x}, t) = \int_{S(0,1)} (\mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p}) \Psi(\mathbf{x}, \mathbf{p}, t) d\mathbf{p} \quad (117) $$

L'équation de Fokker-Planck peut être résolue avec une technique stochastique, ce qui nous renvoie à la simulation de dynamique brownienne décrite dans la Section 4.4. ou avec une technique déterministe comme nous venons de le voir.

<figure><img src="image_35.png" /><figcaption>Figure 33. Evolution du domaine fluide.</figcaption></figure>

Les Figures 33 à 36 présentent quelques images de la simulation micro-macro du remplissage d'un
moule 3D avec une suspension de fibres. Tous les détails concernant la discrétisation du modèle continu ainsi que la gestion du front de matière se trouvent dans la référence [52].

<figure><img src="image_36.png" /><figcaption>Figure 34. Evolution de l'orientation des fibres.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>Figure 35. Evolution de l'orientation des fibres.</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>Figure 36. Evolution de l'orientation des fibres.</figcaption></figure>

Les Figures 37 et 38 présentent des résultats similaires pour un autre moule tridimensionnel. Dans les figures concernant l'évolution de l'orientation des fibres, nous avons procédé au calcul du tenseur d'orientation du deuxième ordre à partir de la fonction de distribution $\Psi(x, p, t)$, elle-même solution de l'équation de Fokker-Planck (116). La diagonalisation de $^2a(x, t)$ nous offre trois valeurs propres ($\lambda_1, \lambda_2, \lambda_3$) et trois vecteurs propres associés. Dans les figures où nous montrons l'évolution de l'orientation, nous avons représenté le vecteur propre associé à la plus grande valeur propre. Le code de couleur traduit l'intensité de l'orientation

dans la direction du vecteur tracé, c'est-à-dire, la valeur propre la plus grande. La couleur rouge signifie que la fibre se trouve dans la direction du vecteur représenté, par contre la couleur bleue indique une orientation plutôt isotrope ($\lambda_{max} = 1/3$).

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

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

<figure><img src="image_41.png" /><figcaption>Figure 37. Evolution du domaine fluide.</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>Figure 38. Evolution de l'orientation des fibres.</figcaption></figure>
Pour simuler des écoulements à surface libre, comme ceux qui ont lieu lors de l'extrusion de polymères, nous avons proposé le couplage entre la résolution de l'équation de Fokker-Planck et des techniques numériques adaptées à la résolution des équations du mouvement (Eqs. (113) et (114)), tout en mettant l'accent sur la description très précise des trajectoires de l'écoulement. Nous allons justifier la nécessité de cette représentation et le choix des techniques pour y parvenir. La nécessité est simple : les équations constitutives (macroscopiques ou celles basées sur l'approche microscopique issue de la théorie cinétique) sont en général purement convectives dans l'espace physique. Ainsi, si l'on procède avec une technique eulérienne, la principale difficulté sera la stabilisation de la formulation discrète qui introduit, dans le meilleur des cas, une inévitable diffusion numérique. L'intégration par la méthode des caractéristiques des équations de convection semble être la plus précise, mais pour la mettre en œuvre, nous avons besoin d'une description très fine des trajectoires de l'écoulement.

Pour y parvenir nous nous sommes focalisés sur deux stratégies. La première concerne l'utilisation de la technique dite des tubes de courant, proposée il y a quelques années par J-R. Clermont et ses collaborateurs au Laboratoire de Rhéologie de Grenoble. Cette technique permet de transformer le domaine physique dans un domaine de référence où les lignes coordonnées correspondent aux trajectoires de l'écoulement, permettant une intégration très précise des équations constitutives, dans notre cas l'intégration de l'équation de Fokker-Planck par une technique de séparation de variables. Les principaux résultats de ce couplage sont décrits dans la référence [53].

La deuxième alternative que nous avons explorée est celle d'utiliser une description lagrangienne actualisée des équations du mouvement. Ce type de description est délicat si l'on procède avec la méthode des éléments finis. En effet, dans cette description, les nœuds sont amenés à suivre l'écoulement et donc le maillage associé à ces nœuds devient très rapidement trop distordu. Nous rappelons que la précision des solutions obtenues avec les méthodes basées sur l'emploi d'un maillage (éléments finis, différences finies, volumes finis, ...) est subsidiaire de la qualité géométrique du maillage subjacent. Ainsi, dès que le maillage devient trop distordu, une procédure de remaillage doit être appliquée.

Les opérations de remaillage introduisent de nombreuses difficultés :

- le remaillage est par nature délicat en trois dimensions ;

- il est cher à calculer ;

- les variables connues aux nœuds de l'ancien maillage doivent être projetées sur les nouveaux nœuds et cette opération est délicate puisqu'elle induit une diffusion numérique significative ...

Pour remédier à ces difficultés, nous avons proposé d'utiliser des méthodes dites « meshless » qui nous permettent d'aller plus loin dans les simulations malgré les fortes distorsions des maillages subjcents qui peuvent, pour certaines techniques « sans maillage », exister (comme c'est le cas de la méthode des éléments naturels). La « quid de la question » n'est pas l'existence ou non de maillage subjacent; le point clé est le fait que la qualité de ce maillage (quand il existe) joue ou non un rôle décisif dans la qualité des résultats obtenus. De mon point de vue, cet aspect fait bien la différence entre les méthodes standard et celles dites « sans maillage », mais ce n'est que ma vision personnelle des faits ! Dans le cas des éléments finis, on sait que la qualité du maillage est décisive, d'où la nécessité de remailler, et donc l'impossibilité de mener à terme une simulation avec le même nuage de nœuds et sans aucune opération de repositionnement des nœuds. Au contraire, les méthodes dites « sans maillage » sont « a priori » capables d'un tel pari ! Cependant, de nombreuses difficultés subsistent. Voici les plus importantes :

- l'imposition des conditions aux limites ;

- l'intégration numérique ;

- le choix du support des fonctions de forme des approximations résultantes ;

- les approximations mixtes et les conditions de stabilité LBB ; ...

La méthode des éléments naturels, introduite dans le domaine des simulations des écoulements complexes par E. Cueto et ses collaborateurs à l'université de Saragosse en Espagne, permet de surpasser la plupart de ces difficultés que nous venons de souligner [54].

Elle nous a permis l'analyse fine des écoulements à rhéologie complexe et à surfaces libres. Pour une description détaillée de cette méthode « sans maillage » et de son application à la simulation des écoulements de mise en forme, le lecteur peut s'adresser aux références [54-61].

Bien entendu, dans certains de nos travaux, nous avons comparé les prédictions numériques de l'orientation de fibres avec celles observées
expérimentalement [62]. Comme les bons numériciens, nous n'attribuons qu'une validité limitée à nos prédictions !

## 7. Conclusions

Nous voici à la fin de notre promenade où nous avons traversé les différentes échelles de description des fluides complexes, avec le « coarse graining » comme fil conducteur. En effet, les simulations atomistiques ne sont possibles que dans des systèmes de taille très (trop) réduite et, même dans ce cas, comme nous l'avons prouvé dans la deuxième section, nous ne pouvons pas compter sur les vrais potentiels interatomiques provenant des calculs à l'échelle quantique.

Il a fallu donc simplifier les descriptions, premièrement en éliminant certaines entités (par exemple les particules de solvant) ou bien en simplifiant les interactions des particules retenues dans la modélisation. Ainsi, nous sommes arrivés aux simulations de dynamique brownienne qui, malgré leur simplicité conceptuelle, sont encore très (trop) gourmandes en temps de calcul.

L'échelle suivante est celle de la théorie cinétique, où l'on revient au milieu continu, avec une fonction (parfois plusieurs) continue décrivant la conformation microstructurale du système et dont l'évolution est gouvernée par l'équation de Fokker-Planck, qui n'est rien d'autre que l'équation de conservation dans un espace qui regroupe les coordonnées physiques (espace et temps) et les coordonnées dites de conformation, qui décrivent la configuration microscopique du système étudié. A nouveau, l'approche est très simple côté conceptuel mais le prix à payer est la multidimensionnalité des équations résultantes de Fokker-Planck. Deux possibilités existent et ont été recensées plus haut : revenir à la dynamique brownienne en profitant de l'équivalence existant entre l'équation de Fokker-Planck et les équations stochastiques d'Ito [37] ; ou bien rester dans le domaine déterministe et utiliser des techniques novatrices et performantes pour traiter ce type de problèmes à grande dimension. Nous avons présenté l'une de ces techniques, à notre connaissance la plus puissante, que nous avons illustrée sur un exemple en chimie quantique et sur un exemple en théorie cinétique, bien que de nombreuses autres applications aient été abordées tout au long de notre recherche. Enfin, la Figure 39 caricature la démarche suivie dans ce travail qui, malgré tout, reste ouvert à toutes les échelles !

<figure><img src="image_43.png" /><figcaption>Figure 39. De l'atome à la pièce : une caricature !</figcaption></figure>

## Remerciements

Bien entendu, les travaux reportés dans cet article ne sont pas l'œuvre d'une seule personne, loin de là. Ce travail témoigne d'une étroite collaboration avec plusieurs chercheurs reconnus, qui s'est prolongée tout au long de plusieurs années. Je voudrais remercier tout d'abord mon principal co-équipier, le Professeur Amine Ammar du laboratoire de Rhéologie (CNRS-UJF-INPG) de Grenoble

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

ainsi que (par ordre alphabétique) mes principaux collaborateurs tout au long de ces années et auxquels je dois beaucoup : Elias Cueto (Saragosse, Espagne), Roland Keunings (Louvain la Neuve, Belgique), Malcolm Mackley (Cambridge, Royaume Uni), Arnaud Poitou (Centrale Nantes) et Claude Wolff (CNAM Paris, à présent à la retraite).

Je voudrais également remercier encore d'autres collaborateurs avec lesquels j'ai cosigné de nombreux travaux : J. Azaiez (Calgary, Canada), K. Chiba (Osaka, Japon), S. Cescotto (Liège, Belgique), Th. Coupez (CEMEF - Mines de Paris), N. El Kissi (Grenoble), J.A. Garcia (Valence, Espagne), M.
Grmela (Montréal, Canada), A. Hamdouni (La Rochelle), A. Huerta (Barcelone, Espagne), P. Joyot (Bayonne), M. Kroger (ETH Zurich), M. Laso (Madrid, Espagne), M. Pasquali (Houston, Etats-Unis), T. Phillips (Cardiff, Royaume Uni), G. Regnier (ENSAM Paris), D. Ryckelynck (Mines de Paris), N. Sukumar (Californie, Etats-Unis), R. Torres (Valence, Espagne), F. Trochu (Montréal, Canada), P. Villon (UTC Compiègne), K. Walters (Aberystwyth, Royaume Uni), ... et j'en passe, en m'excusant auprès des oubliés !

Finalement je tiens à remercier le Groupe Français de Rhéologie de m'avoir décerné cette distinction dont je ne suis pas sûr de l'avoir méritée, quand je me compare à ceux qui l'ont reçue avant moi, ou encore à ceux qui devront sans doute la recevoir les prochaines années.

Références

[1] Weiner J.H., *Statistical Mechanics of Elasticity*, Dover (2002).

[2] Laughlin R.B., The theory of everything, Proceedings of the USA National Academy of Sciences, 97 (2000).

[3] Ammar A., Mokdad B., Chinesta F., Keunings, R., A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids, *J. Non-Newt.Fluid Mech.*, 139, 153-176 (2006).

[4] Bungartz H.J., Griebel M., Sparse grids, *Acta Numer.*, 13, 1–123 (2004).

[5] Achdou Y., Pironneau O., Computational methods for option pricing, *SIAM Front. Appl. Math.*, (2005).

[6] Ladevèze P., *Nonlinear Computational Structural Mechanics*, Springer, New York (1999).

[7] Ammar A., Chinesta F., Circumventing curse of dimensionality in the solution of highly multidimensional models encountered in quantum mechanics using meshfree finite sums decomposition, *Lecture Notes Comput. Sci. Eng.*, Springer, 1-17 (2008).

[8] Chinesta F., Ammar A., Joyot P., The nanometric and micrometric scales of the structure and mechanics of materials revisited: An introduction to the challenges of fully deterministic numerical descriptions, *Intern. J. Multiscale Comput. Eng.*, 6/3, 191-213 (2008).

[9] Cook D.B., *Handbook of Computational Quantum Chemistry*, Dover (2005).

[10] Cancès E., Defranceschi M., Kutzelnigg W., Le Bris C., Maday Y., Computational quantum chemistry: A primer. In: *Handbook of Numerical Analysis*, Vol. X. Elsevier, New York, 3–270 (2003).

[11] Le Bris C., *Handbook of Numerical Analysis*, Vol. X: Computational Chemistry, Elsevier, New York (2003).

[12] Chantrenne P., Volz S., *Thermique à l'échelle submicronique. Introduction à la dynamique moléculaire*. Techniques de l'Ingénieur, Traité Génie Energétique, BE-290, 1-20.

[13] Laso M., Muneta L.M., Müller V., Alcazar V., Chinesta F., Ammar A., Hierarchical approach to flow calculations for polymeric liquid crystals. In: *Integrated Multiscale Modelling for Polymer Properties Design*, E.A. Perpète, M. Laso, (Eds.), Elsevier, Vol. 22, 359-402 (2006).

[14] Jeffery G.B., The motion of ellipsoidal particles immersed in viscous fluid, *Proc. R. Soc., A102*, 161 (1922).

[15] Hakim R., Introduction à la Mécanique Statistique, Masson, Paris (1996).

[16] Bird R.B., Curtiss C.F., Armstrong R.C., Hassager O., *Dynamic of Polymeric Liquids*, Vol. 2: Kinetic Theory, Wiley, New York (1987).

[17] de Gennes P.G., Reptation of a polymer chain in the presence of fixed obstacles, *J. Chem. Phys.* 55 572-579 (1971).

[18] Doi M., Edwards S.F., *The Theory of Polymer Dynamics*, Oxford Science Publications (1978).

[19] Fang J., Kroger M., Ottinger H.C., A thermodynamically admissible reptation model for fast flows of entangled polymers. II. Model predictions for shear and extensional flows, *J. Rheol.*, 40, 1293-1318 (2000).

[20] Keunings R., Micro-macro methods for the multiscale simulation viscoelastic flow using molecular models of kinetic theory, *Rheol. Rev.*, DM Binding, K Walters (Eds.), Brit. Soc. Rheol., 67-98 (2004).

[21] Dealy, J.M., Larson, R.G., Structure and Rheology of Molten Polymers: From Structure to Flow Behavior and Back Again, Hanser, Munich (2006)

[22] Wetzel E.D., Tucker Ch. L. III, Area tensors for modelling microstructure during laminar liquid-liquid mixing. *Int. J. Multiph. Flow*, 25, 35–61 (2000).

[23] Chinesta F., Mackley M., Microstructure evolution during liquid-liquid laminar mixing: A kinetic theory approach, *Intern. J. Mater. Form.*, 1, 47-55 (2008).

[24] Chella R., Ottino J.M., Stretching in some classes of fluid motions and asymptotic mixing efficiencies as a measure of flow classification, *Arch. Ration. Mech. Anal.*, 90, 15–42 (1985).

[25] Doi M., Ohta T., Dynamics and rheology of complex interfaces, *J. Chem. Phys.*, 95, 1242–1248 (1991).

[26] Ma A., Chinesta F., Mackley M., The microstructure and rheology of carbon nanotube suspensions, In: *Science and Applications of Carbon Nanotubes*, M. Chipara (Ed.), Elsevier. A paraître.
[27] Ma A., Chinesta F., Tuladhar T., Mackley M., Filament stretching of carbon nanotube suspensions, *Rheol. Acta*, 47, 447-457 (2008).

[28] Cueto E., Ma A., Chinesta F., Mackley M., Numerical simulation of spin coating processes with functionalized carbon nanotubes suspensions, *Intern. J. Mater. Form.*, 2, 89-99 (2008).

[29] Ma A., Chinesta F., Mackley M., Ammar A., The rheological modelling of carbon nanotube (CNT) suspensions in steady shear flows, *Intern. J. Mater. Form.*, 2, 83-88 (2008).

[30] Ma A., Mackley M., Chinesta F., The microstructure and rheology of carbon nanotube suspensions, *Intern. J. Mater. Form.*, 2, 75-81 (2008).

[31] Ma A., Chinesta F., Ammar A., Mackley M., Rheological modelling of carbon nanotube aggregate suspensions, *J. Rheol.*, 52/6, 1311-1330 (2008).

[32] Ma A., Chinesta F., Mackley M., The rheology and modelling of chemically treated carbon nanotube suspensions, *J. Rheol.*, à paraître.

[33] Larson R.G., *The Structure and Rheology of Complex Fluids*, Oxford University Press (1999).

[34] Chauviere C., Lozinski A., Simulation of dilute polymer solutions using a Fokker-Planck equation, *Comput. Fluids*, 33, 687-696 (2004).

[35] Ammar A., Chinesta F., Ryckelynck D., Deterministic particle approach of multi-bead-spring polymer models, *Europ. J. Comput. Mech.*, 15, 481-494 (2006).

[36] Chinesta F., Chaidron G., Poitou A., On the solution of the Fokker-Planck equations in steady recirculating flows involving short fiber suspensions, *J. Non-Newt. Fluid Mech.*, 113, 97-125 (2003).

[37] Ottinger H.C., Stochastic processes in polymeric fluids, Springer (1996).

[38] Ottinger H.C., Laso M., Smart polymers in finite element calculation, Int Congr. on Rheology, Brussel, Belguim (1992).

[39] Hulsen M.A., van Heel A.P.G., van der Brule B.H.A.A., Simulation of viscoelastic flows using Brownian Configuration Fields", *J. Non-Newt. Fluid Mech.*, 70, 79-101 (1997).

[40] Chinesta F., Ammar A., Falco A., Laso M., On the reduction of stochastic kinetic theory models of complex fluids, *Model. Simul. Mater. Sci. Eng.*, 15, 639-652 (2007).

[41] Ammar A., Mokdad B., Chinesta F., Keunings R., A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Part II: Transient simulation using space-time separated representation, *J. Non-Newt. Fluid Mech.*, 144, 98-121 (2007).

[42] Mokdad B., Pruliere E., Ammar A., Chinesta F., On the simulation of kinetic theory models of complex fluids using the Fokker-Planck approach, Appl. Rheo., 17, 26494, 1-14 (2007).

[43] Ryckelynck D., Chinesta F., Cueto E., Ammar A., On the "a priori" model reduction: Overview and recent developments, *Archives Comput. Methods Eng.*, State of the Art Reviews, 13, 91-128 (2006).

[44] Ammar A., Ryckelynck D., Chinesta F., Keunings R., On the reduction of kinetic theory models related to finitely extensible dumbbells, *J. Non-Newt. Fluid Mech.*, 134, 136-147 (2006).

[45] Pruliere E., Ammar A., Chinesta F., Empirical natural closure relation for short fiber suspension models, *Intern. J. Form. Proc.*, 10, 361-385 (2007).

[46] Kroger M., Ammar A., Chinesta F., Consistent closure schemes for statistical models of anisotropic fluids, *J. Non-Newt. Fluid Mech.*, 149, 40-55 (2008).

[47] Pruliere E., Ammar A., El Kissi N., Chinesta F., Multiscale modelling of flows involving short fiber suspensions, *Archives Comput. Methods Eng.*, State of the Art Reviews, 16, 1-30 (2009).

[48] Chiba K., Ammar A., Chinesta F., On the fiber orientation in steady recirculating flows involving short fibers suspensions, *Rheol. Acta*, 44, 406-417 (2005).

[49] Keunings R., On the Peterlin approximation for finitely extensible dumbbells, *J. Non-Newt. Fluid Mech.*, 68, 85-100 (1997).

[50] Folgar F.P., Tucker III Ch.L., Orientation behavior of fibers in concentrated suspensions, *J. Reinf. Plast. Comp.*, 3, 98-119, 1984.

[51] Dupret F., Verleye V., Modelling the flow of fiber suspensions in narrows gaps, In: *Advances in the Flow and Rheology of Non-Newtonian Fluids*. D.A. Siginer, D. De Kee, R.P. Chabra (Eds.), Rheology Series, Elsevier, 1347-1398 (1999).

[52] Azaiez J., Chiba K., Chinesta F., Poitou A., State-of-the-Art on numerical simulation of fiber-reinforced thermoplastic forming processes, *Archives Comput. Meth. Eng.*, 9, 141-198 (2002).

[53] Mokdad B., Ammar A., Normandin M., Chinesta F., Clermont J.R., A fully deterministic micro-macro simulation of complex flows involving reversible network fluid models, *Math. Comput. Simul.*, à paraître.

[54] Chinesta F., Cueto E., Lorong Ph., Cescotto S., *La Méthode des Eléments Naturels dans la Simulation des Structures et des Procédés*, Hermes, Paris, à paraître.

[55] Martinez M.A., Cueto E., Doblaré M., Chinesta F., Natural element meshless simulation of injection processes involving short fiber suspensions, *J. Non-Newt. Fluid Mech.*, 115, 51-78 (2003).

[56] Martinez M.A., Cueto E., Alfaro I., Doblaré M., Chinesta F., Updated Lagrangian free surface flow
simulations with the natural neighbour Galerkin methods,
*Intern. J. Numer. Meth. Eng.*, 60, 2105-2129 (2004).

[57] Chinesta F., Cueto E., Ryckelynck D., Ammar A., Alpha-NEM and model reduction: Two new and powerful numerical strategies for describing flows involving short fiber suspensions, *Revue Europ. Eléments Finis*, 14, 903-923 (2005).

[58] Alfaro I., Gonzalez D., Bel D., Cueto E., Doblare M., Chinesta F., Recent advances in the meshless simulation of extrusion and other related forming processes, *Archives Comput. Methods Eng.*, State of the Art Reviews, 13, 3-44 (2006).

[59] Gonzalez D., Cueto E., Chinesta F., Doblaré M., An updated Lagrangian strategy for free-surface fluid dynamics, *J. Comput. Physics*, 223, 127-150 (2007).

[60] Garcia J.A., Gascón L., Cueto E., Ordeig I., Chinesta F., Meshless methods with application to resin transfer molding simulation, *Comput. Meth. Appl. Mech. Eng.*, à paraître.

[61] Cueto E., Laso M., Chinesta F., Meshless stochastic simulation of micro-macro kinetic theory models, *Intern. J. Multiscale Comput. Eng.*, à paraître.

[62] Poitou A., Chinesta F., Bernier G., Orienting fibers by extrusion in reinforced reactive powder concrete, *J. Eng. Mech.*, 127, 593-598 (2001).