# Mise en forme des composites à matrice polymère renforcée avec des fibres courtes : Vingt ans de simulation numérique

Francisco Chinesta, Amine Ammar, Kunji Chiba

► To cite this version:

Francisco Chinesta, Amine Ammar, Kunji Chiba. Mise en forme des composites à matrice polymère renforcée avec des fibres courtes : Vingt ans de simulation numérique. Revue des composites et des matériaux avancés = Journal of Composite and Advanced Materials, 2005, 15, pp.323-337. 10.3166/rcma.15.323-337. hal-01577584

HAL Id: hal-01577584

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

Submitted on 26 Aug 2017

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Mise en Forme des Composites à Matrice
Polymère Renforcée avec des Fibres Courtes :
20 Ans de Simulation Numérique.

Francisco Chinesta<fnref n="1" />, Amine Ammar<fnref n="2" />, Kunji Chiba<fnref n="3" />

RESUME. Les suspensions de fibres courtes sont abondamment impliquées dans l'industrie des matériaux composites. La connaissance de la microstructure finale des pièces est essentielle pour la prédiction des propriétés mécaniques. De même, la prédiction d'une telle microstructure pourra aussi servir à l'optimisation des procédés afin d'obtenir les propriétés finales souhaitées. Dans ce papier nous faisons un tour d'horizon sur les dernières avancées dans la simulation numérique des suspensions de fibres courtes.

ABSTRACT. Short fiber suspensions are widely used in industry for manufacturing composites parts. The structures involving short fiber composites exhibit anisotropy induced by the forming process, because of the coupling between the flow kinematics and the fiber orientation evolution. The final microstructure determines the mechanical properties of the conformed part. Thus, accurate simulation tools allow for optimizing the forming processes in order to improve the mechanical properties. In this work we summarize some of the most significant recent advancements in this field

MOTS CLES : Suspension polymère - fibres courtes, Polymère, Simulation Numérique, Orientation des fibres, Composites à fibres courtes.

KEYWORDS: Short fiber polymer suspension, Polymer, Numerical modeling, Fiber orientation, Short fiber composites.

<footnote n="1">¹ LMSP: UMR CNRS – ENSAM – ESEM
151 Boulevard de l’Hôpital, F-75013 Paris, France.
francisco.chinesta@paris.ensam.fr</footnote>

<footnote n="2">² Laboratoire de Rhéologie
1301 rue de la piscine, BP 53 Domaine universitaire
38041 Grenoble Cedex 9, France
Amine.Ammar@ujf-grenoble.fr</footnote>

<footnote n="3">³ Faculty of Education, Shiga University
2-5-1 Hiratsu, Otsu, Shiga 520-0862, Japan
kchiba@sue.shiga-u.ac.jp</footnote>
# 1. Introduction

Dans la famille des fluides complexes, dont font partie les suspensions polymères - fibres courtes, les matériaux ont des longueurs caractéristiques nettement plus grandes que celles des fluides simples. Le couplage entre les microstructures (concentration des fibres et leur état d'orientation) et les écoulements conduit souvent à des comportements rhéologiques particuliers.

Les suspensions de fibres courtes sont abondamment impliquées dans l'industrie des matériaux composites. La connaissance de la microstructure finale des pièces est essentielle pour la prédiction des propriétés mécaniques. De même, la prédiction d'une telle microstructure pourra aussi servir à l'optimisation des procédés afin d'obtenir les propriétés finales souhaitées.

La qualité d'une telle prédiction repose sur la maîtrise de différents aspects :

1. Modélisation mécanique de la suspension,

2. Identification des paramètres rhéologiques,

3. Utilisation des techniques numériques appropriées.

Bien que de grands progrès aient été faits depuis quelques années, on n'est pas encore aujourd'hui en mesure de donner une réponse définitive à ces trois aspects.

En terme de modélisation mécanique, on préfère encore travailler dans une échelle macro plutôt que dans une échelle micro même si différentes approches micro-macro ont été proposées ces dernières années (Keunings, 2004). Les approches purement macro reposent sur une homogénéisation spatiale des équations de comportement où la présence des fibres est prise en compte par une homogénéisation à caractère statistique (au travers de l'introduction des tenseurs d'orientation). En général ce type de modèles nécessite la considération d'un fluide newtonien. Du côté micro, différents travaux ont porté sur la simulation directe de l'orientation d'une fibre dans un écoulement de fluide newtonien à partir de la résolution de l'équation de Jeffery. Le passage micro-macro est alors fait en moyennant, dans chaque point de l'espace, l'orientation des différentes fibres (méthodes de type statistique et de simulation directe). La prise en compte du caractère non-newtonien du solvant nécessite une simulation directe. En effet il n'existe pas d'équations offrant l'évolution des fibres dans l'écoulement d'un tel fluide. De la même façon, la prise en compte des interactions entre fibres, fibres-parois, ... est encore mal maîtrisée du point de vue macro.

En terme d'identification, la problématique que nous venons d'évoquer reste toujours d'actualité. Aux difficultés relatives à l'obtention d'un état de contraintes homogène dans un fluide, il faut ajouter la présence des fibres. La validité des essais rhéologiques les plus habituels est mise en cause comme nous l'avons prouvé dans nos derniers travaux. D'autres essais proposés dans le domaine des polymères
pourront être étendus aux composites à fibres courtes à l'état fondu, mais sans doute des difficultés de mise en place seront au rendez-vous.

Pour finir, et concernant les aspects purement numériques, de grands progrès ont été faits ces dernières années. La puissance du calcul informatique (pour un coût de plus en plus faible) permet la simulation de problèmes de grande taille. Ainsi, différents calculs, à différentes échelles, sont aujourd'hui possibles. Les nouvelles techniques numériques (méthodes sans maillage, méthodes issues de la dynamique moléculaire, techniques de réduction de modèles ...) étendront encore plus le domaine de simulation. Aujourd'hui, sans aucun doute, les limitations de la modélisation résident plus dans les aspects physiques que dans les aspects numériques.

## 2. Modélisation mécanique

Les équations de conservation de la masse et de la quantité de mouvement pour un matériau incompressible s'écrivent :

$$ \mathrm{Div} \underline{v} = 0 \qquad [1] $$

$$ \rho \frac{dv}{dt} = -\mathrm{Grad} p + \mathrm{Div} \underline{\tau} \qquad [2] $$

où $\underline{v}$ est le vecteur vitesse, $\rho$ la densité, $p$ la pression, $\underline{\tau}$ le tenseur des extra-contraintes et $d/dt$ la dérivée matérielle. La théorie des milieux continus développée par Lipscomb et al. (1988) dans le cadre des suspensions diluées de particules rigides avec un rapport de forme élevé est souvent utilisée pour prendre en compte l'existence des fibres en suspension

$$ \underline{\tau} = 2\eta \underline{D} + \phi \underline{\mu} \underline{D} : \underline{a} = 2\eta \underline{D} + 2\eta N_p \underline{D} : \underline{a} \qquad [3] $$

où $\eta$ est la viscosité de la matrice fluide supposée Newtonienne, $\underline{D}$ est le tenseur des taux de déformation, $N_p$ un paramètre matériau qui dépend de la concentration de fibres ainsi que de leur rapport de forme, $\phi$ est la fraction (en volume) de fibres, $\underline{\mu}$ est un paramètre du matériau et $\underline{a}$ est le tenseur d'orientation d'ordre 4

$$ \underline{a} = \oint \underline{p} \otimes \underline{p} \otimes \underline{p} \otimes \underline{p} \psi(\underline{p}) d\underline{p} \qquad [4] $$
où ρ̄ est le vecteur unitaire qui défini l'orientation d'une fibre et ψ(ρ̄) est la fonction de distribution de l'orientation des fibres dans un élément de fluide. Pour des fibres à rapport de forme élevé (rapport entre la longueur et le diamètre) λ ≥ 10, μ peut être donnée par (Lipscomb et al., 1988)

$$ \mu = \frac{\eta \lambda^2}{\ln(\lambda)} \qquad [5] $$

La fonction de distribution de l'orientation ψ(ρ̄) doit vérifier la condition de normalité :

$$ \int \psi(\underline{\rho}) d\underline{\rho} = 1 \qquad [6] $$

Si ψ(ρ̄) = δ(ρ̄ - ρ̂), avec δ() la distribution de Dirac, toute la probabilité d'orientation se trouve concentrée dans la direction définie par ρ̂, et dans ce cas le tenseur d'orientation qui en résulte est : $\hat{a} = \rhô \otimes \rhô \otimes \rhô \otimes \rhô$.

Nous pouvons également définir le tenseur d'orientation du second ordre :

$$ \underline{a} = \int \rho \otimes \rho \psi(\underline{\rho}) d\underline{\rho} \qquad [7] $$

et dans le cas où ψ(ρ̄) = δ(ρ̄ - ρ̂), nous pouvons vérifier facilement l'égalité suivante

$$ \underline{a} = \underline{a} \otimes \underline{a} \qquad [8] $$

dont les composantes sont données par $a_{ijkl} = a_{ij}a_{kl}$.

Dans le cas général la relation [8] n'est pas exacte, elle définie alors une fermeture dite quadratique, que nous notons par $a_{ijkl}^{quad} = a_{ij}a_{kl}$. Il y a cependant d'autres relations de fermeture (Advani et Tucker, 1990 ; Dupret et al. 1998 ; Cintra et Tucker 1995), parmi lesquelles nous pouvons citer la fermeture linéaire

$$ a_{ijkl}^{\text{lin}} = -\frac{1}{24}(\delta_{ij}\delta_{kl} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) + \frac{1}{6}(a_{ij}\delta_{kl} + a_{ik}\delta_{jl} + a_{il}\delta_{jk} + a_{kl}\delta_{ij} + a_{jl}\delta_{ik} + a_{jk}\delta_{il}) \qquad [9] $$

la fermeture hydride

$$ a_{ijkl}^{\text{hyb}} = f a_{ijkl}^{\text{quad}} + (1-f) a_{ijkl}^{\text{lin}} \qquad [10] $$
ou encore la fermeture naturelle (Dupret et al., 1998) qui prend en 2D la forme

$$a_{ijkl}^{nat} = \frac{1}{6}\det(\underline{a})(\delta_{ij}\delta_{kl} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) + \frac{1}{3}(a_{ij}a_{kl} + a_{ik}a_{jl} + a_{il}a_{jk}) \quad [11]$$

Pour des distributions isotropes, $\Psi(\underline{\rho}) = \frac{1}{4\pi}$ (en 3D), le tenseur d'orientation du second ordre qui en résulte est $\underline{a} = \frac{I}{3}$ et nous pouvons vérifier facilement que dans ce cas la fermeture linéaire est exacte.

Du point de vue physique, nous pouvons imaginer que les valeurs propres du tenseur d'orientation du second ordre nous donnent l'intensité de l'orientation sur la direction des vecteurs propres associés.

Jusqu'à présent nous avons discuté l'effet de l'orientation des fibres sur les extra-contraintes, et par conséquent sur la cinématique de l'écoulement. Dans les paragraphes qui suivent, nous allons nous focaliser sur l'effet réciproque, c'est à dire, sur l'effet de la cinématique de l'écoulement soit sur l'orientation des fibres individuellement, soit sur l'évolution des différents tenseurs d'orientation.

Si on considère que les fibres admettent une forme ellipsoïdale, l'évolution de leurs orientations est alors donnée par l'équation de Jeffery (Jeffery, 1922)

$$\frac{d\underline{\rho}}{dt} = \underline{\Omega}\underline{\rho} + k(\underline{D}\underline{\rho} - (\underline{D}: (\underline{\rho}\otimes\underline{\rho})))\underline{\rho} \quad [12]$$

où $\underline{\Omega}$ est le tenseur de vorticité et $k$ une constante qui dépend du facteur de forme des fibres $\lambda$ :

$$k = (\lambda^2 - 1)/(\lambda^2 + 1) \quad [13]$$

Par ailleurs, l'évolution de la fonction de distribution de l'orientation des fibres est donnée par l'équation de Fokker-Planck

$$\frac{d\psi(\underline{\rho})}{dt} + \frac{\partial}{\partial \underline{\rho}} \left\{ \psi(\underline{\rho}) \frac{d\underline{\rho}}{dt} \right\} = 0 \quad [14]$$

où la dérivée matérielle est donnée par :

$$\frac{d\psi}{dt} = \frac{\partial\psi}{\partial t} + \underline{v} \text{Grad}\psi \quad [15]$$
Maintenant en combinant les équations [7], [12] et [14], nous pouvons en déduire :

$$ \frac{da}{dt} = \underline{\underline{\Omega}} a - a \underline{\underline{\Omega}} + k \left( \underline{\underline{D}} a + a \underline{\underline{D}} - 2 \left( a : \underline{\underline{D}} \right) \right) \quad [16] $$

Nous pouvons également en déduire une équation similaire gouvernant l'évolution du tenseur d'orientation du quatrième ordre. On fera alors intervenir le tenseur d'orientation du sixième ordre. Il faut donc avoir recours à une relation de fermeture.

Pour prendre en compte les effets d'interaction entre fibres, Folgar et Tucker (1984) ont proposé l'introduction d'un terme de diffusion dans l'équation de Fokker-Planck :

$$ \frac{d\psi(\underline{\underline{\rho}})}{dt} + \frac{\partial}{\partial \underline{\underline{\rho}}} \left\{ \psi(\underline{\underline{\rho}}) \frac{d\underline{\underline{\rho}}}{dt} \right\} = \frac{\partial}{\partial \underline{\underline{\rho}}} \left\{ D_r \frac{\partial \psi(\underline{\underline{\rho}})}{\partial \underline{\underline{\rho}}} \right\} \quad [17] $$

L'équation gouvernant l'évolution du tenseur d'orientation prend alors dans ces conditions la forme :

$$ \frac{da}{dt} = \underline{\underline{\Omega}} a - a \underline{\underline{\Omega}} + k \left( \underline{\underline{D}} a + a \underline{\underline{D}} - 2 \left( a : \underline{\underline{D}} \right) \right) - 4 D_r \left( a - \frac{I}{N} \right) \quad [18] $$

avec

$$ N = \begin{cases} 2 & \text{en } 2D \\ 3 & \text{en } 3D \end{cases} \quad [19] $$

L'utilisation du formalisme de Fokker-Planck évite la nécessité d'introduire des relations de fermeture, par contre, l'équation devient à caractère multidimensionnel (la fonction de distribution est définie tant dans l'espace physique que dans celui des configurations). Pour cette raison, le nombre de travaux portant sur la résolution de l'équation de Fokker-Planck reste assez réduit. Nous y reviendrons dans la section 4.

# 3. Modélisation numérique

## 3.1. Modèle découplé

Dans les premières simulations concernant les suspensions de fibres courtes l'évolution de l'orientation des fibres était obtenue en résolvant l'équation de Jeffery tout en supposant que les fibres n'ont aucun effet sur la cinématique de
l'écoulement. Cette dernière était obtenue en résolvant les équations de Stokes. Un des travaux pionniers dans cette approche est celui de Givler et al. (1982). Connaissant l'orientation des fibres dans la frontière communicante au domaine, l'équation d'évolution de l'orientation des fibres est intégrée tout au long des trajectoires de l'écoulement. Si l'orientation initiale est donnée par une distribution, on doit évaluer l'orientation de chacune des fibres initialement placées dans chaque point du bord d'entrée définissant la distribution connue. Nous pouvons résumer la procédure numérique dans le cas où on considère une distribution isotrope à l'entrée (dans le cas 2D) :

* A chaque pas de temps

□ Pour chaque particule fluide $k$ contenant $N$ fibres, qui décrivent à l'instant intial la distribution d'orientation supposée connue, on effectue les étapes suivantes :

➤ Calculer la vitesse et le gradient de vitesse au point où la particule fluide se trouve à l'instant actuel : $v_k^n$ et $\text{Grad}v|_{\underline{x}_k^n}$

➤ Actualiser la position de la particule fluide : $\underline{x}_k^{n+1} = \underline{x}_k^n + \underline{v}_k^n \Delta t$

➤ Pour chacune des $N$ fibres associées à la particle fluide $k$ :

❖ Calculer la vitesse de rotation par l'équation de Jeffery :

$$ \dot{\theta}_{i,k}^{n} = \dot{\theta}_{i,k}^{n} \left( \text{Grad} v_{\underline{x}_k^{n-1}}^{n} , k, \theta_{i,k}^{n} \right) $$

❖ Actualiser l'orientation :

$$ \theta_{i,k}^{n+1} = \theta_{i,k}^{n} + \dot{\theta}_{i,k}^{n} \Delta t $$

➤ Calculer, si l'on désire, les différents tenseurs d'orientation associés à la particule fluide $k$, par exemple $a_{\equiv k}^{n+1} = \int \underline{\rho} \otimes \underline{\rho} \otimes \underline{\rho} \otimes \underline{\rho} \, \psi_k^{n+1} \, d\underline{\rho}$, avec $\psi_k^{n+1}(\theta) = \sum_{i=1}^{i=N} \frac{1}{N} \delta(\theta - \theta_{i,k}^{n+1})$, dont la forme discrète s'écrit :

$$ \left\{ \begin{array}{l} a_{1111,k}^{n+1} = \sum_{i=1}^{i=N} \cos^4 \theta_{i,k}^{n+1} \frac{1}{N} \\ \\ a_{2222,k}^{n+1} = \sum_{i=1}^{i=N} \sin^4 \theta_{i,k}^{n+1} \frac{1}{N} \end{array} \right. \quad \left\{ \begin{array}{l} a_{1112,k}^{n+1} = \sum_{i=1}^{i=N} \cos^3 \theta_{i,k}^{n+1} \sin \theta_{i,k}^{n+1} \frac{1}{N} \\ \\ a_{1122,k}^{n+1} = \sum_{i=1}^{i=N} \cos^2 \theta_{i,k}^{n+1} \sin^2 \theta_{i,k}^{n+1} \frac{1}{N} \\ \\ a_{1222,k}^{n+1} = \sum_{i=1}^{i=N} \cos \theta_{i,k}^{n+1} \sin^3 \theta_{i,k}^{n+1} \frac{1}{N} \end{array} \right. $$

avec les relations de symétrie : $a_{1112} = a_{1121} = a_{1211} = a_{2111}$ ;

$$ a_{1122} = a_{1221} = a_{2211} = a_{2112} = a_{1212} = a_{2121} ; \dots $$

Malgré la simplicité de l'approche, le fait de supposer que l'existence des fibres et leur état d'orientation n'affectent pas la cinématique de l'écoulement peut
s'avérer inapproprié. Pour certains écoulements, le fait d'introduire une faible quantité de fibres dans un fluide Newtonien peut changer significativement le comportement cinématique, comme on peut l'observer dans la figure 1.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Visualisation de l'écoulement d'un fluide Newtonien et d'une suspension de fibres dans une contraction 4 : 1 (Lipscomb et al., 1988).</figcaption></figure>

## 3.2. Modèle couplé avec hypothèse d'alignement local

Les premières simulations couplées ont introduit l'orientation des fibres dans l'écoulement tout en supposant que les fibres sont localement (dans chaque point de l'écoulement) alignées dans la direction du vecteur vitesse (Lipscomb, 1988). De cette manière le tenseur d'orientation s'écrit :

$$ \underline{a} = \frac{\underline{v} \underline{v}^T}{\| \underline{v} \|^2} \qquad [20] $$

Pour l'introduction des effets d'orientations dans les extra-contraintes il suffit de noter que l'alignement local des fibres implique, comme prouvé précédemment, la validité de la fermeture quadratique, d'où :

$$ \underline{\sigma} = -p\underline{I} + 2\eta\{\underline{D} + N_p\left(\underline{a} : \underline{D}\right)\} = -p\underline{I} + 2\eta\{\underline{D} + N_p\left(\underline{a} : \underline{D}\right)\underline{a}\} \qquad [21] $$

Malgré les simplifications effectuées dans le modèle, ses prédictions permettent de quantifier correctement la différence de contraintes normales en cisaillement
(observée expérimentalement) ainsi que la croissance de la zone de recirculation avec la concentration des fibres (traduite par le paramètre $N_p$) dans l'écoulement en contraction, comme on peut l'observer dans la figure 2.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Lignes de courant calculées et observées : (a) fluide Newtonien et (b) suspension de fibres caractérisée par $N_p/\eta = 4$ (Chiba et al. 1990).</figcaption></figure>

Dans le cas où l'hypothèse d'alignement local est considérée, le problème résultant correspond à celui de Stokes anisotrope, qui peut être résolu avec n'importe quelle méthode adaptée à cette équation. Dans nos simulations nous avons utilisé une formulation mixte par éléments finis, avec une interpolation fonctionnelle vérifiant la condition de stabilité LBB (liée aux discrétisations des formulations variationnelles mixtes). Pour mettre en évidence l'équivalence entre le modèle précédent et celui de Stokes avec une viscosité anisotrope nous pouvons réécrire les extra-contraintes de l'équation [21], avec une représentation vectorielle, sur la forme :

$$ \begin{align} \underline{\tau} = 2\eta \left\{ \underline{D} + N_p \left( \underline{a} : \underline{D} \right) \right\} &\Leftrightarrow \underline{\tau} = 2\eta \left\{ \underline{D} + N_p \left( \underline{a}^T \underline{D} \right) \underline{a} \right\} &\Leftrightarrow \nonumber \\ &\Leftrightarrow \underline{\tau} = 2\eta \left\{ \underline{I} + N_p \underline{a} \underline{a}^T \right\} \underline{D} = 2\eta^* \underline{D} \tag{22} \end{align} $$

Cette méthode simplifiée pour modéliser le couplage donne de bons résultats. Cependant, il a été prouvé que l'alignement local des fibres avec l'écoulement n'est qu'une solution locale pour des fibres a rapport de forme infini (Poitou et al. 2000); et de plus, celle-ci peut être instable dans des écoulements à caractère divergent par exemple. Il est bien connu qu'une fibre à facteur de forme infini tend à s'aligner avec l'écoulement en cisaillement et en elongation « positive » tandis qu'elle se met perpendiculairement à l'écoulement dans le cas d'une elongation « négative », comme nous illustrons dans la figure 3.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Tendance de l'orientation d'une fibre dans trois types d'écoulements : (gauche) cisaillement ; (centre) élongation positive et (droite) élongation négative</figcaption></figure>

Dans le cas d'un écoulement complexe, il est difficile de prévoire la validité de l'hypothèse d'alignement local. Ainsi, son utilisation peut induire de grandes différences entre la cinématique calculée et observée, d'où la nécessité de procéder avec des modèles couplés.

### 3.3. Modèles couplés

#### 3.3.1. Simulations particulaires

La simulation introduite dans la section 3.1. peut être appliquée sans trop de changement à la simulation couplée de l'écoulement. Le modèle couplé est souvent résolu en forme découplée en utilisant un schéma incrémental à chaque itération (ou pas de temps) dans lequel on calcule la cinématique (avec l'état d'orientation connue à l'itération précédente) ensuite l'orientation (avec la cinématique qui vient d'être calculée). Afin de procéder au couplage effectif entre l'orientation calculée dans la section 3.1 et la cinématique, nous ajoutons à la fin de l'algorithme le calcul du terme des extra-contraintes qui se fait dans la formulation variationnelle des équations de mouvement :

$$ \int_{\Omega} \left[ 2\eta \underline{D} + 2\eta N_p \underline{a} : \underline{D} \right] : \underline{D}^* d\Omega \Leftrightarrow \int_{\Omega} (\underline{D}^*)^T [2\eta \underline{I} + 2\eta N_p \underline{A}] \underline{D} d\Omega \quad [23] $$

où la contribution des fibres est déterminée à partir de :
$$ \underline{a} : \underline{D} = \underline{a}_{ijkl} \underline{D}_{kl} = \underline{a}_{ij11} \underline{D}_{11} + \underline{a}_{ij12} \underline{D}_{12} + \underline{a}_{ij21} \underline{D}_{21} + \underline{a}_{ij22} \underline{D}_{22} $$

$$ \underline{D}^* : (\underline{a} : \underline{D}) = D_{11}^* (a_{1111} \underline{D}_{11} + a_{1112} \underline{D}_{12} + a_{1121} \underline{D}_{21} + a_{1122} \underline{D}_{22}) + \\ + D_{12}^* (a_{1211} \underline{D}_{11} + a_{1212} \underline{D}_{12} + a_{1221} \underline{D}_{21} + a_{1222} \underline{D}_{22}) + \\ + D_{21}^* (a_{2111} \underline{D}_{11} + a_{2112} \underline{D}_{12} + a_{2121} \underline{D}_{21} + a_{2122} \underline{D}_{22}) + \\ + D_{22}^* (a_{2211} \underline{D}_{11} + a_{2212} \underline{D}_{12} + a_{2221} \underline{D}_{21} + a_{2222} \underline{D}_{22}) = \qquad [24] $$

$$ (D_{11}^* \ D_{12}^* \ D_{21}^* \ D_{22}^*) \begin{pmatrix} a_{1111} & a_{1112} & a_{1121} & a_{1122} \\ a_{1211} & a_{1212} & a_{1221} & a_{1222} \\ a_{2111} & a_{2112} & a_{2121} & a_{2122} \\ a_{2211} & a_{2212} & a_{2221} & a_{2222} \end{pmatrix} \begin{pmatrix} D_{11} \\ D_{12} \\ D_{21} \\ D_{22} \end{pmatrix} = (\underline{D}^*)^T \underline{A} \underline{D} $$

Cette technique a été utilisée avec succès dans les travaux de Chiba (Chiba et al. 1990 ; Chiba et al. 2001 ; Chiba et Chinesta, 2005) pour l'analyse de différents écoulements, avec suspensions mono-, bi- ou multi-disperse. La seule limitation est a priori liée au fait que le modèle décrit dans la section 3.1 que nous venons de reprendre ici n'offre pas la possibilité d'interactions entre les fibres. Comme nous l'avons précisé plus haut, cette interaction est à l'origine de l'introduction des termes de diffusion dans l'équation de Fokker-Planck [17] ou dans l'équation d'évolution du tenseur d'orientation du second ordre associée [18]. Pour s'affranchir de cette problématique, il y a deux possibilités simples : (i) introduire le terme de diffusion de l'équation de Fokker-Planck [17] dans le terme de convection en modifiant l'équation de Jeffery (voir Ammar et Chinesta, 2005) ; ou (ii) transformer l'équation [17] en une équation stochastique (comme l'on fait d'habitude dans le cadre des techniques CONNFFESSIT qui est à l'origine des techniques dites micro-macro (Keunings, 2004 ; Laso et Ottinger, 1993). Dans le cas des fibres, cette dernière technique est équivalente à un traitement de type « random walks » de la diffusion (voir Chinesta et al., 2003). L'application de ces stratégies est loin d'être facile : dans le premier cas on aura besoin d'évaluer le gradient de la distribution des orientations (qui est connue en forme discrète) et dans la deuxième stratégie on a besoin d'un grand nombre de fibres pour bien représenter les effets de diffusion, c'est pour cela que la plupart des travaux de recherche se sont orientés vers les modèles basés sur les équations d'évolution du tenseur d'orientation.

### 3.3.2. Simulations avec des équations d'évolution du tenseur d'orientation

Dans ce cas le modèle couplé est à nouveau résolu à partir d'une stratégie découplée (Rosenberg et al. 1990 ; Altan et al. 1992 ; Azaiez et al., 1997 ; Azaiez et al. 2002):

- On résout à chaque itération (ou pas de temps) la cinématique donnée par les équations [1], [2] et [3] (avec l'orientation connue à l'itération ou pas de temps précédent).
- En utilisant ensuite la cinématique donnée par le problème de Stokes anisotrope, l'orientation peut être actualisée en résolvant l'équation [18].

Mais le modèle défini par les équations. [1], [2], [3] et [18] n'est pas fermé. Il a fallut au préalable exprimer le tenseur du quatrième ordre en fonction de celui du deuxième ordre en utilisant une relation de fermeture (quelques exemples ont été donnés dans la section 2). La seule difficulté numérique est reliée au caractère hyperbolique non-linéaire de l'équation [18]. L'utilisation avec succès de la technique de Lesaint-Raviart a permis la simulation d'écoulements complexes en extrusion ou en injection. Un exemple de remplissage est montré dans la figure 4. Nous avons constaté que la progression du front n'est pas sensiblement affectée par la présence des fibres.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Simulation du remplissage d'un moule 3D avec une suspension de fibres courtes (Monton, 2004).</figcaption></figure>

L'orientation la plus probable des fibres est tracée dans la figure 5 (direction des flèches), où la couleur des flèches indique l'intensité d'orientation dans la direction représentée.

Bien que cette modélisation a permis d'effectuer de nombreuses simulations et d'obtenir de bonnes corrélations avec les résultats expérimentaux, certains problèmes persistent encore aujourd'hui : dans les écoulements recirculants (ou dans les zones de recirculation dans des écoulements généraux) l'orientation n'est pas définie au centre de la recirculation. Cette singularité diffuse dans le domaine de l'écoulement perpendiculairement aux lignes de courant. La problématique liée aux écoulements recirculants a été traitée en profondeur dans Chinesta (2003). L'autre limitation est associée à la nécessité d'utiliser une relation de fermeture dont l'erreur, bien qu'elle ne soit pas trop significative dans la simulation des procédés,
s'avère très perturbatrice dans les écoulements rhéométriques comme prouvé dans
(Chiba, Ammar et Chinesta, 2005).

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Prédiction de l'orientation des fibres (Monton, 2004).</figcaption></figure>

**4. Nouvelles tendances.**

Le formalisme de Fokker-Planck permet de s'affranchir de l'introduction d'une relation de fermeture mais rend, en contre partie, la simulation d'une telle équation à caractère multi-dimensionnel délicate. Une possibilité consiste à effectuer le couplage entre une technique Lagrangienne de description de la cinématique (sans maillage de préférence) avec une technique de particules à partir de l'équation de Jeffery (modifiée pour prendre en compte les effets des interactions). Un exemple de couplage entre une technique d'éléments naturels (sans maillage) pour le calcul de la cinématique avec une évaluation de l'orientation basée sur l'équation de Jeffery peut être trouvé dans Chinesta et al. (2005). Une autre possibilité consiste à procéder directement avec l'équation de Fokker-Planck en utilisation une technique de réduction de modèles « a priori » afin de pouvoir rendre la simulation précise mais aussi efficace du point de vue temps de calcul (le lecteur intéressé peut s'adresser à Chinesta et al., (2005)). Finalement, si on veut étudier en détail de l'interaction entre fibres, ou entre fibres et parois, la simulation directe semble être l'approche la plus efficace, puisqu'elle permet de valider et d'identifier les
paramètres des lois macroscopiques. Les recherches entamées par Megally et al. (2005) portent sur cette problématique.

Malgré tous ces commentaires à caractère purement numérique, d'autres questions plus physiques persistent, dont un exemple est celui de la rhéologie des suspensions concentrées où, de plus, la matrice fluide est viscoélastique. Bien que certains modèles proposent une forme du tenseur des extra-contraintes, le problème général du mouvement d'une fibre immergée dans un écoulement viscoélastique est encore ouvert. Même si la simulation directe permet de prédire le mouvement des fibres, l'obtention d'une équation macroscopique (équivalente à celle de Jeffery) est encore aujourd'hui fortement attendue.

## 5. Références

Advani, S.G., Tucker III, C.L., « Closure Approximations for Three-dimensional Structure Tensors », *J. Rheol.*, vol. 34, 1990, p. 367.

Ammar, A., Chinesta, F., « A Particle Strategy for Solving the Fokker-Planck Equation Governing the Fiber Orientation Distribution in Steady Recirculating Flows Involving Short Fiber Suspensions », *Lectures Notes on Computational Science and Engineering*, Springer, vol. 43, 2005, p. 1.

Altan, M.C., Güceri, S.I., Pipes, R.P., « Anisotropic Channel Flow of Fiber Suspensions », *J. Non-Newtonian Fluid Mech.*, vol. 42, 1992, p. 65.

Azaiez, J., Guénette, R., Ait-Kadi, A., « Investigation of the Abrupt Contraction Flow of Fiber Suspensions in Polymeric Fluids », *J. Non-Newtonian Fluid Mech.*, vol. 73, 1977, p. 289.

Azaiez, J., Chiba, K., Chinesta, F., Poitou, A., « State-of-the-Art on Numerical Simulation of Fiber-Reinforced Thermoplastic Forming Process », *Arch. Comput. Meth. Engng.*, vol. 9, 2002, p. 141.

Chiba, K., Nakamura, K., Boger, D.V., « A Numerical Solution for the Flow of Dilute Fiber Suspensions Through an Axisymmetric Contraction », *J. Non-Newtonian Fluid Mech.*, vol. 35, 1990, p. 1.

Chiba, K., Yasuda, K., Nakamura, K., « Numerical Solution of Fiber Suspension Flow Through a Parallel Plate Channel by Coupling Flow Field with Fiber Orientation Distribution », *J. Non-Newtonian Fluid Mech.*, vol. 99, 2001, p. 145.

Chiba, K., Chinesta, F., « Numerical Simulation of Flow Kinematics and Fiber Orientation for Multi-Disperse Suspension », *Rheologica Acta*, 2005, A paraître.

Chiba, K., Ammar, A., Chinesta, F., « On the Fiber Orientation in Steady Recirculating Flows Involving Short Fibers Suspensions », *Rheologica Acta*, vol. 44, 2005, p. 406.

Chinesta, F., Chaidron, G., Poitou, A., « On the solution of the Fokker-Planck Equations in Steady Recirculating Flows Involving Short Fiber Suspensions », *J. Non-Newtonian Fluid Mech.*, vol. 11, n. 2-3, 2003, p. 97.
Chinesta, F., Cueto, E., Ryckelynck, D., Ammar, A., « Alpha-NEM and Model Reduction: Two New and Powerful Numerical Strategies to Describe Flows Involving Short Fiber Suspensions », *Revue Européenne des Eléments Finis*, 2005, A paraître.

Cintra, J.S., Tucker III, C.L., « Orthotropic Closure Approximations for Flow Induced Fiber orientation », *J. Rheology*, vol. 39, 1995, p. 1095.

Dupret, F., Verleye, V., Languillier, B., « Numerical Prediction of the Moulding of Short Fiber Composite Parts », *First ESAFORM Conference on Material Forming*, Sophia-Antipolis, France (1998).

Folgar, F., Tucker III, C.L., « Orientation Behaviour of Fibers in Concentrated Suspensions », *J. Reinf. Plast. Composites*, vol. 3, 1984, p.98.

Givler, C., Crochet, M.J., Pipes, R.B., « Numerically Predicted Fiber Orientation in Dilute Suspensions », *NUMIFORM*, Pineridge Press, 1982, p. 559.

Jeffery, G.B., « The Motion of Ellipsoidal Particles Immersed in a Viscous Fluid », *Proc. Roy. Soc. London*, A102, 1922, p. 161.

Keunings, R. « Micro-Macro Methods for the Multiscale Simulation of Viscoelastic Flow using Molecular Models of Kinetic Theory », *Rheology Reviews*, British Society of Rheology, 2004, p. 67.

Laso, M., Ottinger, H.C., « Calculation of Viscoelastic Flow Using Molecular Models », *J. Non-Newtonian Fluid Mech.*, vol. 47, 1993, p. 1.

Lipscomb, G.G., Denn, M.M., Hur, D.H., Boger, D.V., « The Flow of Fiber Suspensions in Complex Geometries », *J. Non-Newtonian Fluid Mech.*, vol. 26, 1988, p. 297.

Megally, A., Laure, P., Coupez, Th., « Direct Flow Simulation of Fibers in a Molten Polymer Fluid », *2nd Annual European Rheology Conference*, Grenoble, 2005.

Monton, I. Simulation 3D de l'Injection de Composites Fibres Courtes (espagnol), Thèse de Doctorat à l'Université Polytechnique de Valence, 2004.

Poitou, A., Chinesta, F., Torres, R., « Numerical Simulation of the Steady Recirculating Flows of Fibers Suspensions », *J. Non-Newtonian Fluid Mech.*, vol. 90, 2000, p. 65.

Rosenberg, J., Denn, M.M., Keunings, R., « Simulation of Non-Recirculating Flows of Dilute Fiber Suspensions », *J. Non-Newtonian Fluid Mech.*, vol. 37, 1990, p. 317.