National audienceUne méthode de simulation numérique directe a été développée pour l'écoulement de fluides fortement chargés, avec plusieurs populations de particules (fibres et sphères), pouvant entrer en collision. Ce fluide complexe modélise les écoulements de BMC (Bulk Molding Compound) utilisés par Schneider Electric. Une approche par domaines fictifs, avec des conditions limites périodiques, a été utilisée

Beaume, Grégory

Coupez, Thierry

Laure, Patrice

HAL-MINES ParisTech

Simulations nume´riques d’e´coulements de fluidescomplexes a` l’e´chelle microscopique : un nouvel outil derhe´ologieGre´gory Beaume, Patrice Laure, Thierry CoupezTo cite this version:Gre´gory Beaume, Patrice Laure, Thierry Coupez. Simulations nume´riques d’e´coulements defluides complexes a` l’e´chelle microscopique : un nouvel outil de rhe´ologie. Groupe Franc¸ais deRhe´ologi. 40e`me colloque annuel du Groupe Franc¸ais de Rhe´ologie, Oct 2005, Nice, France. 4p., 2005. <hal-00578626>HAL Id: hal-00578626https://hal-mines-paristech.archives-ouvertes.fr/hal-00578626Submitted on 21 Mar 2011HAL is a multi-disciplinary open accessarchive for the deposit and dissemination of sci-entific research documents, whether they are pub-lished or not. The documents may come fromteaching and research institutions in France orabroad, or from public or private research centers.L’archive ouverte pluridisciplinaire HAL, estdestine´e au de´poˆt et a` la diffusion de documentsscientifiques de niveau recherche, publie´s ou non,e´manant des e´tablissements d’enseignement et derecherche franc¸ais ou e´trangers, des laboratoirespublics ou prive´s.Simulations numériques d’écoulements de fluides complexes à l’échelle microscopique : un nouvel outil de rhéologie  G. Beaume1,2, P. Laure3 et T. Coupez2  1  Schneider Electric-Technocentre, 38 TEC Grenoble  2 CEMEF, Ecole Nationale Supérieure des Mines de Paris, UMR CNRS 7635, 06904, Sophia Antipolis 3  Institut Non-Linéaire de Nice-UMR 6618 CNRS- Université de Nice Sophia Antipolis, 06560 Valbonne Résumé : Une méthode de simulation numérique directe a été développée pour l’écoulement de fluides fortement chargés, avec plusieurs populations de particules (fibres et sphères), pouvant entrer en collision. Ce fluide complexe modélise les écoulements de BMC (Bulk Molding Compound) utilisés par Schneider Electric. Une approche par domaines fictifs, avec des conditions limites périodiques, a été utilisée. Mots-clé : Rhéologie, BMC, Fluides chargés, Fibres, Simulation numérique directe, Domaines fictifs, Eléments finis  1. Introduction  1.1 Objectif  L’étude présentée dans cet article porte sur la modélisation et la simulation d’écoulements de fluides complexes. L’objectif est de développer un outil de simulation pour l’écoulement d’un fluide modèle, représentatif d’une catégorie de matériaux utilisés par Schneider Electric pour la fabrication de pièces moulées par le procédé d’injection. Les matériaux considérés sont les BMC (Bulk Molding Compound). Ce sont des résines thermodurcissables contenant divers types de charges, en grande concentration : de l’alumine, des charges minérales, des nodules thermoplastiques, des fibres de renforts … Les BMC sont ici modélisés par une matrice fluide newtonienne, fortement chargée de plusieurs populations de particules : des particules longitudinales (fibres) et des sphères de diverses tailles (les autres charges). 1.2 Simulation Numérique Directe : principe et intérêt Il existe dans la littérature de nombreux modèles de comportement macroscopiques pour les fluides newtoniens chargés de fibres et/ou de particules sphériques, mais uniquement pour des concentrations faibles ou modérées (un récapitulatif de ces modèles est proposé dans [1]). Ces modèles sont construits par homogénéisation : les différentes grandeurs (définies à l’échelle microscopique) sont moyennées sur un volume élémentaire représentatif (macroscopique) ou sur toutes les orientations possibles (description macroscopique de l’orientation). Aucun de ces modèles ne parvient à représenter le comportement de fluides fortement chargés. En effet, ils ne prennent pas en compte les interactions à courte portée (hydrodynamiques ou de contact) existant entre les particules. Une première alternative à cette difficulté théorique consiste à construire des modèles de comportement approximatifs (fluide à seuil, loi puissance,…) à partir d’observations expérimentales. Cette approche est aussi limitée avec les BMC, dont la rhéologie est particulièrement difficile (matériaux très pâteux, vieillissement de la résine…). Une seconde approche, retenue et présentée dans cet article, consiste à construire une loi de comportement à partir de simulations numériques directes : il s’agit ici de simuler directement l’écoulement du fluide et des particules qu’il contient à l’échelle microscopique [2]. On peut ensuite moyenner les différentes grandeurs qui nous intéressent (viscosité effective, tenseur d’orientation des fibres,…) sur un volume élémentaire et obtenir une loi de comportement macroscopique : on parle alors d’homogénéisation numérique. Dans cet article nous présentons la méthode numérique que nous avons développée pour la simulation numérique directe d’écoulements fortement chargés, avec plusieurs populations de particules, de formes différentes. Nous présentons ensuite les résultats préliminaires obtenus en 2D, avec des concentrations solides importantes (de 30 à 65 %) et une seule population de charges. 2. Positionnement du problème  On va résoudre l’écoulement considéré sur un domaine spatial Ω  de géométrie simple (généralement parallélépipédique), séparé  en deux sous-domaines, désignant respectivement les domaines fluide fΩ  et solide sΩ . Le domaine solide s’écrit comme la réunion des domaines   associées à chaque particule i présente dans le domaine total : Uniiss1= Ω=Ω      (1) Chaque sous-domaine est régi par des équations d’évolution qui lui sont propres. 2.1 Domaine fluide  Le fluide est supposé newtonien et incompressible. En injection, les forces d’inertie et de gravité sont négligées et le mouvement du fluide est correctement décrit par les équations de Stokes : 0. =∇ ur      (2) 0.r=∇σ      (3)) uIp fr&(2 εησ +−= )    (4) 2.2 Domaine solide  Chaque particule est caractérisée par sa masse Mi, la position de son barycentre iXr (et sa vitesse barycentrique Ui) et, pour les particules non-sphériques (donc les fibres), par son vecteur orientation ipr (et le vecteur rotation iωr ). En négligeant les forces d’inertie et de gravité pour les particules, nous supposons que nous sommes en régime hydrodynamique, malgré les grandes concentrations en charges considérées : on parle alors de régime lubrificationnel. Les particules ne sont donc soumises, en théorie, qu’aux seules forces hydrodynamiques. Néanmoins, dans la littérature [3], d’autres auteurs ont constaté la formation d’amas de particules au cours de leurs simulations. Ces amas n’ont aucun sens physique; ils sont d’origine numérique : pour des pas de temps trop grands, les forces lubrificationnelles, fortement répulsives à très courte portées, sont sous-estimées lors du rapprochement de deux particules. Pour s’affranchir de ce phénomène, une méthode consiste à rajouter des forces correctives entre particules. Ainsi, les équations du mouvement pour la particule i sont les suivantes : Pour une sphère : 0²² rrrr ≈+= icpihii FFdtXdM    (5) Fhi désigne la force hydrodynamique et Fcpi la force corrective subie par la particule i du fait de sa proximité des particules voisines. Pour une fibre, deux équations supplémentaires sont nécessaires pour suivre l’évolution de l’orientation : ⎪⎪⎩⎪⎪⎨⎧ ×=≈+=iiiicpihiipdtpdΓΓdtJdrrrrrrrωω0)(   (6) Γih et Γcpi désignent respectivement les couples associés aux forces hydrodynamiques et correctives. 2.3 Conditions limites  Sur la frontière Γ du domaine, les conditions limites choisies correspondent à un écoulement de Couette. On impose donc la vitesse du fluide aux parois horizontales de façon à avoir un taux de cisaillement moyen indépendant du temps dans le domaine : ⎩⎨⎧ == 0yxu yu γ&      (7) On impose par ailleurs des conditions limites pseudo-périodiques dans la direction x. En réalité, la périodicité ne porte que sur l’évolution des particules : toute particule sortant d’un côté du domaine est réintroduite de l’autre côté avec la même vitesse. Le champ de vitesse sur les parois verticales est laissé libre (d’où le terme pseudo périodicité). De manière générale, les conditions limites périodiques permettent de ramener un calcul sur un domaine spatial infini (donc sans effets de bord) à un domaine spatial plus petit. Ce type de conditions est parfaitement adapté pour l’homogénéisation numérique. Enfin, les conditions limites aux interfaces fluide/particules traduisent le mouvement de corps rigide des charges. On a donc la relation suivante : isiii XxUu Ω∂−×+= sur)( rrrrr ω   (8) 3. Méthode numérique  On simule donc l’écoulement du fluide et des particules qu’il contient, sur le domaine et sur un intervalle de temps donné, qui sont discrétisés pour la résolution  numérique. Le problème spatial à chaque pas de temps est résolu par éléments finis. 3.1 Méthode avec remaillage La méthode intuitive pour ce type de problème consiste à découper la résolution à chaque pas de temps en deux étapes : le calcul du champ de vitesse sur fΩ , puis le déplacement des particules à partir du champ de vitesse fluide et des forces correctives à courte portée. La formulation variationnelle en éléments finis mixtes s’écrit de façon classique [4] : Trouver ),( pur tel que pour tout ),( qvr dans )()( 20 ff LH Ω×Ω  : ∑ ∫∫∫∫Γ−Γ−=Ω∇−ΩΩ∂ΓΩΩifdvndvndvpdvuisff rrrrrr&r&).().(.)(:)(2σσεεη(9) avec : { }isiiiff XxUvHvH Ω∂−×+=Ω∈=Ω sur)()()( 1 rrrrrr ω  Cette approche nécessite un remaillage du domaine fluide à chaque pas de temps, puisque fΩ  évolue, et est donc lourde numériquement. C’est pourquoi on lui préfère la méthode sans remaillage dite "par domaine fictif". 3.2 Méthode sans remaillage : technique des domaines fictifs La technique des domaines fictifs consiste à étendre la formulation variationnelle à Ω tout entier, en faisant disparaître les intégrales aux interfaces fluide/particule [4]. En effet, la condition limite à ces interfaces, et les approximations de l’injection permettent d’écrire : icpiicpiihiihi FUFUvnisΓ−−=Γ+=∫Ω∂ rrrrrrrrrr ....)..( ωωσ  (10) Pour étendre le champ de vitesse au domaine solide, on peut remarquer que la condition de mouvement de corps rigide est équivalente aux deux conditions d’incompressibilité et d’indéformabilité [5]. Sur sΩ  nous avons donc : ⎩⎨⎧ ==∇ 0)( 0.uur&rε      (11) Ces 2 conditions sont des contraintes à imposer au champ de vitesse sur sΩ  pour avoir un mouvement de corps rigide. La première est déjà imposée sur fΩ et doit donc être vérifiée sur tout le domaine. Elle est imposée à l’aide du multiplicateur pression p. Nous avons choisi de pénaliser la contrainte d’indéformabilité, à l’aide du facteur de pénalisation α. Nous obtenons ainsi la formulation suivante : Trouver ),( pur tel que pour tout ),( qvr dans )()( 201 Ω×Ω LH  : ∑∑∫∫∫∫Γ++Γ−=Ω∇−Ω+ΩΓΩΩ ΩΩ ΩiicpiiicpifFUdvndvpdvuIdvuIsfrrrrrrrr&r&r&r&..).(.)(:)(2)(:)(2ωσεεαεεη(12) Aux termes correctifs près, nous avons à résoudre un problème de Stokes multi-domaines ; le facteur de pénalisation joue le rôle d’une viscosité, grande devant celle du fluide. Le problème étant formulé sur un domaine fixe, le remaillage n’est pas nécessaire. En contrepartie, l’évolution des sous-domaines est décrite à l’aide des fonctions caractéristiques sfII ΩΩ , que l’on doit mettre à jour à chaque pas de temps.  Dans notre formulation variationnelle, nous n’avons pas tenu compte des deux derniers termes du membre de droite de (12), correspondant aux forces correctives à courte portée. Nous avons  pris en compte ces interactions au cours du déplacement des particules, en ajoutant simplement des termes correctifs dans les équations d’évolution de iXr et de ipr. 3.4 Evolution du domaine solide  Le calcul des fonctions caractéristiques au pas de temps suivant passe par deux étapes importantes : le calcul des nouvelles positions et orientations des particules, puis la projection des particules dans le maillage global. Les vitesses des particules sont simplement interpolées à partir du champ de vitesse calculé dans toute la cavité, et corrigées en fonction des interactions courte portée. Nous avons  pour la position: [ ]icpiii uXuttXttX rrrrr δ+∆+=∆+ )()()(    (13) avec   ∑≠ →∆= ij ijiicp cpFMtu rrδ    (14) jicpF→rdésigne la force corrective à courte portée exercée par la particule j sur la particule i. Pour des sphères, nous avons pris la force modèle suivante : ⎪⎩⎪⎨⎧ +≥− +<−+−−= →→ dRXXsi dRXXsindXXRF ij ijijijjicp 20 2)2( rrrrrrrrr κ   (15) κ est un paramètre ajustable, R le rayon des particules et d la distance limite entre les deux points les plus proches de la particules, en dessous de laquelle l’interaction est non nulle. Pour des fibres, le centre iXrest remplacé par le point ijPrde l’axe de la fibre i le plus proche de la fibre j. De la même façon, on peut évaluer la nouvelle orientation des particules. Si 1iXr et 2iXrdésignent les extrémités de la fibre i (avec, par convention, 12 iii XXprrr −= ), nous avons : [ ]1212 )()()()()()()(iiiiiicpiiicpiiiXXXuXuttpttptpttpttprrrrrrrrrrrrrr−−∆+×∆+= ×+∆+=∆+ ωδ ωδω      (16)   avec : ∑≠ →×−∆= ij ijiijiicp cpFXPJt rrrr )(ωδ    (17) L’impression des particules dans le maillage du domaine (c’est-à-dire le calcul direct des fonctions caractéristiques) se fait par une méthode de voxelisation [6]. La fonction caractéristique sIΩ est calculée par élément, comme étant le taux de remplissage de l’élément par la phase solide (approximation P0).  4. Résultats préliminaires en 2D Les premières simulations ont été effectuées en 2D, avec une seule population de charge (uniquement des fibres ou  des sphères). A concentration modérée (jusqu’à 30 %), l’ajout de termes correctifs pour les forces de courte portée permet d’empêcher le chevauchement de particules. En l’absence de ces forces, on observe une fusion progressive des particules, qui n’a aucun sens physique.   (a)  (b) Fig. 1 : Ecoulement de sphères (a) et de fibres (b); en haut, nous avons l’évolution avec prise en compte des interactions de courte portée, en bas sans. A grande concentration (64 %), ces chevauchements de particules réapparaissent (Fig. 2). La vitesse de cisaillement est dirigée vers la droite. On s’aperçoit que les chevauchements apparaissent préférentiel-lement en entrée de zone. La formation d’amas est donc probablement liée aux conditions limites choisies ("pseudo périodiques").   Fig. 2 : Ecoulement de sphères à grande concentration; en haut, nous avons l’évolution avec prise en compte des interactions de courte portée, en bas sans. 5. Conclusion et perspectives  Les premières simulations ont permis de valider notre outil de simulation à concentration modérée en charges. Nos travaux actuels visent à améliorer notre méthode, pour pouvoir travailler à des concentrations encore plus élevées. Une majeure partie de ces travaux concerne les conditions limites : ajout de conditions limites pseudo-périodiques sur les autres bord (pour s’affranchir de l’effet de confinement dû à la proximité des parois horizontales), périodicité du champ de vitesse selon la direction de l’écoulement (imposé sous la forme d’une contrainte supplémentaire dans la formulation variationnelle). Une autre partie concerne les forces de courte portée : modèles de forces correctives plus réalistes (notamment en fonction des vitesses relatives) à tester, et prise en compte des forces dans la formulation variationnelle. 6. Références  [1] Ancey, C. Coussot, P. Rhéophysique des pâtes et des suspensions, EDP Sciences (1999). [2] Hu, H.H., Joseph, D.D., Crochet, M.J. Direct simulations of flows of fluid particle motions, Theor. Comput. Fluid Dyn., 3: 285-306 (1992). [3] Hu, H.H. Direct simulations of flows of solid-liquid mixtures, Int. J. Multiphase Flow, 22, 335-352 (1996). [4] Glowinski, G., Pan, T-W., Hesla, T.I., Joseph, D.D. A distributed Lagrange multiplier/fictitious domain method for particulate flow. Int. J. Multiphase Flow, 25, 755-794 (1999). [5] Diaz-Goano, C., Minev, P. D., Nandakumar, K. A fictitious domain/finite element method for particulate flows. J. Comput. Physics, 192, 105-123 (2003).  [6] Megally A. Etude et modélisation de l’orientation de fibres dans les thermoplastiques renforcés, Thèse de Doctorat, Ecole des Mines de Paris, 2005.  

Simulations numériques d'écoulements de fluides complexes à l'échelle microscopique : un nouvel outil de rhéologie

HAL-UNICE

Simulations nume´riques d’e´coulements de fluides
complexes a` l’e´chelle microscopique : un nouvel outil de
rhe´ologie
Gre´gory Beaume, Patrice Laure, Thierry Coupez
To cite this version:
Gre´gory Beaume, Patrice Laure, Thierry Coupez. Simulations nume´riques d’e´coulements de
fluides complexes a` l’e´chelle microscopique : un nouvel outil de rhe´ologie. Groupe Franc¸ais de
Rhe´ologi. 40e`me colloque annuel du Groupe Franc¸ais de Rhe´ologie, Oct 2005, Nice, France. 4
p., 2005. <hal-00578626>
HAL Id: hal-00578626
https://hal-mines-paristech.archives-ouvertes.fr/hal-00578626
Submitted on 21 Mar 2011
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
destine´e au de´poˆt et a` la diffusion de documents
scientifiques de niveau recherche, publie´s ou non,
e´manant des e´tablissements d’enseignement et de
recherche franc¸ais ou e´trangers, des laboratoires
publics ou prive´s.
Simulations numériques d’écoulements de fluides complexes à 
l’échelle microscopique : un nouvel outil de rhéologie 
 
G. Beaume
1,2
, P. Laure
3
 et T. Coupez
2
  
1
  Schneider Electric-Technocentre, 38 TEC Grenoble  
2 CEMEF, Ecole Nationale Supérieure des Mines de Paris, UMR CNRS 7635, 06904, Sophia Antipolis 
3  
Institut Non-Linéaire de Nice-UMR 6618 CNRS- Université de Nice Sophia Antipolis, 06560 Valbonne 
Résumé : Une méthode de simulation numérique directe a été développée pour l’écoulement de fluides fortement 
chargés, avec plusieurs populations de particules (fibres et sphères), pouvant entrer en collision. Ce fluide complexe 
modélise les écoulements de BMC (Bulk Molding Compound) utilisés par Schneider Electric. Une approche par 
domaines fictifs, avec des conditions limites périodiques, a été utilisée. 
Mots-clé : Rhéologie, BMC, Fluides chargés, Fibres, Simulation numérique directe, Domaines fictifs, Eléments finis 
 
1. Introduction  
1.1 Objectif  
L’étude présentée dans cet article porte sur la 
modélisation et la simulation d’écoulements de 
fluides complexes. L’objectif est de développer un 
outil de simulation pour l’écoulement d’un fluide 
modèle, représentatif d’une catégorie de matériaux 
utilisés par Schneider Electric pour la fabrication de 
pièces moulées par le procédé d’injection. Les 
matériaux considérés sont les BMC (Bulk Molding 
Compound). Ce sont des résines thermodurcissables 
contenant divers types de charges, en grande 
concentration : de l’alumine, des charges minérales, 
des nodules thermoplastiques, des fibres de renforts 
… Les BMC sont ici modélisés par une matrice 
fluide newtonienne, fortement chargée de plusieurs 
populations de particules : des particules 
longitudinales (fibres) et des sphères de diverses 
tailles (les autres charges). 
1.2 Simulation Numérique Directe : principe et 
intérêt 
Il existe dans la littérature de nombreux modèles de 
comportement macroscopiques pour les fluides 
newtoniens chargés de fibres et/ou de particules 
sphériques, mais uniquement pour des 
concentrations faibles ou modérées (un récapitulatif 
de ces modèles est proposé dans [1]). Ces modèles 
sont construits par homogénéisation : les différentes 
grandeurs (définies à l’échelle microscopique) sont 
moyennées sur un volume élémentaire représentatif 
(macroscopique) ou sur toutes les orientations 
possibles (description macroscopique de 
l’orientation). Aucun de ces modèles ne parvient à 
représenter le comportement de fluides fortement 
chargés. En effet, ils ne prennent pas en compte les 
interactions à courte portée (hydrodynamiques ou de 
contact) existant entre les particules. Une première 
alternative à cette difficulté théorique consiste à 
construire des modèles de comportement 
approximatifs (fluide à seuil, loi puissance,…) à 
partir d’observations expérimentales. Cette approche 
est aussi limitée avec les BMC, dont la rhéologie est 
particulièrement difficile (matériaux très pâteux, 
vieillissement de la résine…). Une seconde 
approche, retenue et présentée dans cet article, 
consiste à construire une loi de comportement à 
partir de simulations numériques directes : il s’agit 
ici de simuler directement l’écoulement du fluide et 
des particules qu’il contient à l’échelle 
microscopique [2]. On peut ensuite moyenner les 
différentes grandeurs qui nous intéressent (viscosité 
effective, tenseur d’orientation des fibres,…) sur un 
volume élémentaire et obtenir une loi de 
comportement macroscopique : on parle alors 
d’homogénéisation numérique. 
Dans cet article nous présentons la méthode 
numérique que nous avons développée pour la 
simulation numérique directe d’écoulements 
fortement chargés, avec plusieurs populations de 
particules, de formes différentes. Nous présentons 
ensuite les résultats préliminaires obtenus en 2D, 
avec des concentrations solides importantes (de 30 à 
65 %) et une seule population de charges. 
2. Positionnement du problème  
On va résoudre l’écoulement considéré sur un 
domaine spatial Ω  de géométrie simple 
(généralement parallélépipédique), séparé  en deux 
sous-domaines, désignant respectivement les 
domaines fluide fΩ  et solide sΩ . Le domaine 
solide s’écrit comme la réunion des domaines 
  
associées à chaque particule i présente dans le 
domaine total : 
U
n
i
i
ss
1= Ω=Ω      (1) 
Chaque sous-domaine est régi par des équations 
d’évolution qui lui sont propres. 
2.1 Domaine fluide  
Le fluide est supposé newtonien et incompressible. 
En injection, les forces d’inertie et de gravité sont 
négligées et le mouvement du fluide est 
correctement décrit par les équations de Stokes : 
0. =∇ ur      (2) 
0.
r=∇σ      (3)) 
uIp f
r&(2 εησ +−= )    (4) 
2.2 Domaine solide  
Chaque particule est caractérisée par sa masse Mi, la 
position de son barycentre iX
r
 (et sa vitesse 
barycentrique Ui) et, pour les particules non-
sphériques (donc les fibres), par son vecteur 
orientation ip
r
 (et le vecteur rotation iωr ). En 
négligeant les forces d’inertie et de gravité pour les 
particules, nous supposons que nous sommes en 
régime hydrodynamique, malgré les grandes 
concentrations en charges considérées : on parle 
alors de régime lubrificationnel. Les particules ne 
sont donc soumises, en théorie, qu’aux seules forces 
hydrodynamiques. Néanmoins, dans la littérature 
[3], d’autres auteurs ont constaté la formation 
d’amas de particules au cours de leurs simulations. 
Ces amas n’ont aucun sens physique; ils sont 
d’origine numérique : pour des pas de temps trop 
grands, les forces lubrificationnelles, fortement 
répulsives à très courte portées, sont sous-estimées 
lors du rapprochement de deux particules. Pour 
s’affranchir de ce phénomène, une méthode consiste 
à rajouter des forces correctives entre particules. 
Ainsi, les équations du mouvement pour la particule 
i sont les suivantes : 
Pour une sphère : 
0
²
² rrr
r ≈+= icpihii FF
dt
Xd
M    (5) 
Fh
i désigne la force hydrodynamique et Fcp
i la force 
corrective subie par la particule i du fait de sa 
proximité des particules voisines. Pour une fibre, 
deux équations supplémentaires sont nécessaires 
pour suivre l’évolution de l’orientation : 
⎪⎪⎩
⎪⎪⎨⎧ ×=
≈+=
ii
i
i
cp
i
h
ii
p
dt
pd
ΓΓ
dt
Jd
rrr
rrrr
ω
ω
0
)(
   (6) 
Γih et Γcpi désignent respectivement les couples 
associés aux forces hydrodynamiques et correctives. 
2.3 Conditions limites  
Sur la frontière Γ du domaine, les conditions limites 
choisies correspondent à un écoulement de Couette. 
On impose donc la vitesse du fluide aux parois 
horizontales de façon à avoir un taux de cisaillement 
moyen indépendant du temps dans le domaine : 
⎩⎨⎧ == 0yxu yu γ&      (7) 
On impose par ailleurs des conditions limites 
pseudo-périodiques dans la direction x. En réalité, la 
périodicité ne porte que sur l’évolution des 
particules : toute particule sortant d’un côté du 
domaine est réintroduite de l’autre côté avec la 
même vitesse. Le champ de vitesse sur les parois 
verticales est laissé libre (d’où le terme pseudo 
périodicité). De manière générale, les conditions 
limites périodiques permettent de ramener un calcul 
sur un domaine spatial infini (donc sans effets de 
bord) à un domaine spatial plus petit. Ce type de 
conditions est parfaitement adapté pour 
l’homogénéisation numérique. 
Enfin, les conditions limites aux interfaces 
fluide/particules traduisent le mouvement de corps 
rigide des charges. On a donc la relation suivante : 
i
siii XxUu Ω∂−×+= sur)( rrrrr ω   (8) 
3. Méthode numérique  
On simule donc l’écoulement du fluide et des 
particules qu’il contient, sur le domaine et sur un 
intervalle de temps donné, qui sont discrétisés pour 
la résolution  numérique. Le problème spatial à 
chaque pas de temps est résolu par éléments finis. 
3.1 Méthode avec remaillage 
La méthode intuitive pour ce type de problème 
consiste à découper la résolution à chaque pas de 
temps en deux étapes : le calcul du champ de vitesse 
sur fΩ , puis le déplacement des particules à partir 
du champ de vitesse fluide et des forces correctives à 
courte portée. La formulation variationnelle en 
éléments finis mixtes s’écrit de façon classique [4] : 
Trouver ),( pu
r
 tel que pour tout ),( qv
r
 dans 
)()( 20 ff LH Ω×Ω  : 
∑ ∫
∫∫∫
Γ
−Γ−=Ω∇−Ω
Ω∂
ΓΩΩ
i
f
dvn
dvndvpdvu
i
s
ff rr
rrrr&r&
).(
).(.)(:)(2
σ
σεεη
(9) 
avec : { }isiiiff XxUvHvH Ω∂−×+=Ω∈=Ω sur)()()( 1 rrrrrr ω  
Cette approche nécessite un remaillage du domaine 
fluide à chaque pas de temps, puisque fΩ  évolue, et 
est donc lourde numériquement. C’est pourquoi on 
lui préfère la méthode sans remaillage dite "par 
domaine fictif". 
3.2 Méthode sans remaillage : technique des 
domaines fictifs 
La technique des domaines fictifs consiste à étendre 
la formulation variationnelle à Ω tout entier, en 
faisant disparaître les intégrales aux interfaces 
fluide/particule [4]. En effet, la condition limite à 
ces interfaces, et les approximations de l’injection 
permettent d’écrire : 
i
cpi
i
cpi
i
hi
i
hi FUFUvn
i
s
Γ−−=Γ+=∫Ω∂ rrrrrrrrrr ....)..( ωωσ  (10) 
Pour étendre le champ de vitesse au domaine solide, 
on peut remarquer que la condition de mouvement 
de corps rigide est équivalente aux deux conditions 
d’incompressibilité et d’indéformabilité [5]. Sur sΩ  
nous avons donc : 
⎩⎨⎧ ==∇ 0)( 0.uur&
rε      (11) 
Ces 2 conditions sont des contraintes à imposer au 
champ de vitesse sur sΩ  pour avoir un mouvement 
de corps rigide. La première est déjà imposée sur 
fΩ et doit donc être vérifiée sur tout le domaine. 
Elle est imposée à l’aide du multiplicateur pression 
p. Nous avons choisi de pénaliser la contrainte 
d’indéformabilité, à l’aide du facteur de pénalisation α. Nous obtenons ainsi la formulation suivante : 
Trouver ),( pu
r
 tel que pour tout ),( qv
r
 dans 
)()( 20
1 Ω×Ω LH  : 
∑∑∫∫
∫∫
Γ++Γ−=Ω∇−
Ω+Ω
ΓΩ
Ω ΩΩ Ω
i
i
cpi
i
i
cpi
f
FUdvndvp
dvuIdvuI
sf
rrrrrrr
r&r&r&r&
..).(.
)(:)(2)(:)(2
ωσ
εεαεεη
(12) 
Aux termes correctifs près, nous avons à résoudre un 
problème de Stokes multi-domaines ; le facteur de 
pénalisation joue le rôle d’une viscosité, grande 
devant celle du fluide. Le problème étant formulé 
sur un domaine fixe, le remaillage n’est pas 
nécessaire. En contrepartie, l’évolution des sous-
domaines est décrite à l’aide des fonctions 
caractéristiques 
sf
II ΩΩ , que l’on doit mettre à jour 
à chaque pas de temps.  
Dans notre formulation variationnelle, nous n’avons 
pas tenu compte des deux derniers termes du 
membre de droite de (12), correspondant aux forces 
correctives à courte portée. Nous avons  pris en 
compte ces interactions au cours du déplacement des 
particules, en ajoutant simplement des termes 
correctifs dans les équations d’évolution de iX
r
 et de 
ip
r
. 
3.4 Evolution du domaine solide  
Le calcul des fonctions caractéristiques au pas de 
temps suivant passe par deux étapes importantes : le 
calcul des nouvelles positions et orientations des 
particules, puis la projection des particules dans le 
maillage global. Les vitesses des particules sont 
simplement interpolées à partir du champ de vitesse 
calculé dans toute la cavité, et corrigées en fonction 
des interactions courte portée. Nous avons  pour la 
position: [ ]icpiii uXuttXttX rrrrr δ+∆+=∆+ )()()(    (13) 
avec   ∑≠ →∆= ij ijiicp cpFMtu rrδ    (14) 
ji
cpF
→r
désigne la force corrective à courte portée 
exercée par la particule j sur la particule i. Pour des 
sphères, nous avons pris la force modèle suivante : 
⎪⎩⎪⎨
⎧ +≥− +<−+−−= →→ dRXXsi dRXXsindXXRF ij ijijijjicp 20 2)2( rrr
rrrrrr κ   (15) 
κ est un paramètre ajustable, R le rayon des 
particules et d la distance limite entre les deux points 
les plus proches de la particules, en dessous de 
laquelle l’interaction est non nulle. Pour des fibres, 
le centre iX
r
est remplacé par le point ijP
r
de l’axe de 
la fibre i le plus proche de la fibre j. De la même 
façon, on peut évaluer la nouvelle orientation des 
particules. Si 1iX
r
 et 2iX
r
désignent les extrémités de 
la fibre i (avec, par convention, 12 iii XXp
rrr −= ), nous 
avons : [ ]
12
12 )()(
)()(
)()()(
ii
ii
i
i
cpi
i
i
cpiii
XX
XuXu
t
tpttp
tpttpttp
rr
rrrr
rrr
rrrrr
−−∆+
×∆+= ×+∆+=∆+ ωδ ωδω      (16) 
  
avec : 
∑≠ →×−∆= ij ijiijiicp cpFXPJt rrrr )(ωδ    (17) 
L’impression des particules dans le maillage du 
domaine (c’est-à-dire le calcul direct des fonctions 
caractéristiques) se fait par une méthode de 
voxelisation [6]. La fonction caractéristique 
s
IΩ est 
calculée par élément, comme étant le taux de 
remplissage de l’élément par la phase solide 
(approximation P0).  
4. Résultats préliminaires en 2D 
Les premières simulations ont été effectuées en 2D, 
avec une seule population de charge (uniquement 
des fibres ou  des sphères). A concentration modérée 
(jusqu’à 30 %), l’ajout de termes correctifs pour les 
forces de courte portée permet d’empêcher le 
chevauchement de particules. En l’absence de ces 
forces, on observe une fusion progressive des 
particules, qui n’a aucun sens physique.  
 
(a) 
 
(b) 
Fig. 1 : Ecoulement de sphères (a) et de fibres (b); en 
haut, nous avons l’évolution avec prise en compte des 
interactions de courte portée, en bas sans. 
A grande concentration (64 %), ces chevauchements 
de particules réapparaissent (Fig. 2). La vitesse de 
cisaillement est dirigée vers la droite. On s’aperçoit 
que les chevauchements apparaissent préférentiel-
lement en entrée de zone. La formation d’amas est 
donc probablement liée aux conditions limites 
choisies ("pseudo périodiques"). 
 
 
Fig. 2 : Ecoulement de sphères à grande concentration; en 
haut, nous avons l’évolution avec prise en compte des 
interactions de courte portée, en bas sans. 
5. Conclusion et perspectives  
Les premières simulations ont permis de valider 
notre outil de simulation à concentration modérée en 
charges. Nos travaux actuels visent à améliorer notre 
méthode, pour pouvoir travailler à des 
concentrations encore plus élevées. Une majeure 
partie de ces travaux concerne les conditions 
limites : ajout de conditions limites pseudo-
périodiques sur les autres bord (pour s’affranchir de 
l’effet de confinement dû à la proximité des parois 
horizontales), périodicité du champ de vitesse selon 
la direction de l’écoulement (imposé sous la forme 
d’une contrainte supplémentaire dans la formulation 
variationnelle). Une autre partie concerne les forces 
de courte portée : modèles de forces correctives plus 
réalistes (notamment en fonction des vitesses 
relatives) à tester, et prise en compte des forces dans 
la formulation variationnelle. 
6. Références  
[1] Ancey, C. Coussot, P. Rhéophysique des pâtes et des 
suspensions, EDP Sciences (1999). 
[2] Hu, H.H., Joseph, D.D., Crochet, M.J. Direct 
simulations of flows of fluid particle motions, Theor. 
Comput. Fluid Dyn., 3: 285-306 (1992). 
[3] Hu, H.H. Direct simulations of flows of solid-liquid 
mixtures, Int. J. Multiphase Flow, 22, 335-352 (1996). 
[4] Glowinski, G., Pan, T-W., Hesla, T.I., Joseph, D.D. A 
distributed Lagrange multiplier/fictitious domain method 
for particulate flow. Int. J. Multiphase Flow, 25, 755-794 
(1999). 
[5] Diaz-Goano, C., Minev, P. D., Nandakumar, K. A 
fictitious domain/finite element method for particulate 
flows. J. Comput. Physics, 192, 105-123 (2003).  
[6] Megally A. Etude et modélisation de l’orientation de 
fibres dans les thermoplastiques renforcés, Thèse de 
Doctorat, Ecole des Mines de Paris, 2005. 
 


https://hal-mines-paristech.archives-ouvertes.fr/hal-00578626

Simulations numériques d'écoulements de fluides complexes à l'échelle microscopique : un nouvel outil de rhéologie

Abstract

Similar works

Full text

Available Versions

HAL-MINES ParisTech

HAL-UNICE