Procédé d’injection : étude de l’écoulement du flux et des paramètres critiques

 

Elèves: Lounes Bouhadoun - Romain Alis

 

Encadrant ENSEEIHT : Eric CLIMENT

Encadrant industriel : Aurélie Leonardi
 

 


Objectif: Etudier l'injection de colle entre deux plaques planes afin de connaître le temps final d'injection, la répartition des pertes de charge et l'influence des différents paramètres physiques.


 

  

 

Introduction

Pour la fabrication d'un nouvel objet, l'entreprise CLIX-Industries souhaite étudier le collage de deux plaques planes. C'est dans ce cadre que le projet de BEI suivant a été proposé.

Objectifs:

Notre but principal est de trouver un moyen numérique de calculer le temps de remplissage de l'espace entre deux plaques de dimensions données.

CLIX-Industries a réalisé au préalable des tests expérimentaux qui lui ont permis de connaître les temps nécessaires pour que la colle recouvre tout l'espace entre les deux plaques. Ces essais ont été réalisés pour des plaques de petites tailles et nous devons dans un premier temps nous assurer que notre méthode de calcul permet de prévoir les mêmes temps de remplissage.

Lorsque la méthode est validée, nous devons l'utiliser pour prévoir les temps nécessaires correspondants à des configurations de plaques dont il n'y a pas eu de test.

Enfin, nous devons regarder l'effet de l'injection de la colle sur les plaques. En effet, celles-ci ne doivent pas être déformées ni abîmées par l'injection du fluide.

 

Phénomènes indésirables possibles :

Une recherche succincte permet de dégager deux phénomènes qui seraient néfastes à l'injection de colle.

- L'instabilité de Saffman-Taylor (ou phénomène de digitation). L'interface entre les deux phases se déstabilise pour former des doigts qui avancent. Ce phénomène se produit quand la viscosité du fluide qui pousse est inférieure à celle du fluide poussé. Dans notre cas, un liquide (colle) pousse un gaz (air), le phénomène est donc inexistant. Par contre si les rôles étaient inversés, il faudrait étudier cette instabilité pour éviter qu'elle ne se développe.

- La formation d'un film de gaz entre les parois et la colle qui ensuite entraîne l'emprisonnement de bulles de gaz dans le fluide. Ce phénomène dépend fortement de la mouillabilité du fluide, de sa tension de surface et de la vitesse à laquelle il avance. Nous ferons donc attention dans notre étude à ne pas déclencher cet effet. L'idéal serait de trouver les conditions critiques pour lesquels il apparaît.

 

Solutions possibles

Pour répondre aux besoins de CLIX-Industries, nous avons envisagé trois méthodes différentes:

-une simulation de l'injection de colle avec un maillage en trois dimensions de la structure entière.

-une simulation simplifiée de l'injection où l'on simule en trois dimensions les sites importants (buses d'entrée, de sortie, ...) et  une simulation à deux dimensions pour le reste de l'écoulement.

-une méthode où l'on ne calcule pas l'écoulement, on cherche juste à connaître le temps final de remplissage en connaissant le volume total à remplir et le débit de la colle.

 

Simulation complète

C'est la solution qui permet d'avoir les résultats les plus précis car elle permet de simuler l'écoulement dans son ensemble. Toutefois, cette méthode est la plus coûteuse en calculs car elle nécessite de mailler le domaine en entier. Nous n'avons pas pu la mettre en oeuvre car le nombre de mailles requises pour effectuer un calcul précis  dépasse les ordinateurs que nous utilisons à l'ENSEEIHT.

Voici à titre d'exemple, un calcul du nombre de maille d'une configuration dont les dimensions sont: 21 cm en largeur, 30 cm en longueur et 1 mm d'espacement entre les plaques.

Premièrement nous allons mettre 10 mailles sur l'épaisseur (même si cela n'est pas suffisant pour avoir un calcul précis), ce qui donne des mailles d'une hauteur de 0,1 mm.

Donc, pour conserver des mailles de formes carrées (et donc un calcul propre), il faut mettre 3000 mailles sur la longueur et 2100 mailles sur la largeur.

Nous avons donc un total $3000*2100*10=6,3.10^7$ mailles. Le calcul est tout simplement impossible car les licences sont limitées à 512 000 mailles.

Il y a encore la possibilité de dégrader la qualité du maillage (qui est déjà très faible avec 10 mailles sur l'épaisseur) en utilisant des mailles rectangulaires pour diminuer le nombre de mailles sur la longueur et la largeur. Mais cela est inutile car il faudrait diviser chaque coté par environ 10 pour passer sous la limite des licences. Les mailles serait donc trop déformées pour que les calculs soient  valables.

 

 

Simulation partielle

Cette solution a été employée lors de projets numériques réalisés l'ENSEEIHT en 2012/2013. Voici les sites sur lesquels sont présentés les résultats:

http://hmf.enseeiht.fr/travaux/projnum/node/912

http://hmf.enseeiht.fr/travaux/projnum/node/913

La méthode utilisée ici n'a pas permis d'obtenir des résultats concluants et nous allons donc utiliser la dernière méthode pour essayer de prévoir les temps de remplissages des différentes configurations.

Prévision du temps

Le but de cette méthode est de prévoir le temps de remplissage en regardant simplement quel est le temps nécessaire pour remplir l'espace entre les plaques avec un débit $Q(t)$.

Il est donc nécessaire de connaître le débit arrivant par la buse d'entrée pour en déduire le temps final. On connaît le $\Delta P$ imposé par l'industriel, on choisit donc d'utiliser une loi de Poiseuille pour connaître  le débit qui entre dans le domaine.

L'écoulement de Poiseuille est valable pour des cas monophasiques. Il est donc nécessaire de réaliser des simulations pour connaître la relation entre le débit et le gradient de pression en fonction des paramètres de l'écoulement.

Prévision du temps final

Pour connaître le temps final de remplissage nous avons écrit un script Matlab qui calcule le temps final en fonction des paramètres de l'écoulement. L'idée principale est de remplir petit à petit le volume que doit occuper la colle et l'on regarde à quel temps ce volume est plein.  Les étapes du code sont décrites ci -dessous.

On considère le domaine rectangulaire suivant qui est rempli initialement d'un petit volume de colle, l'air est en bleu et la colle en rouge. De plus, on admet qu'il n'y a pas de variation selon l'axe z, c'est à dire la répartition air/colle est la même sur toute l'épaisseur.

On suppose qu'il n'y a aucune perte de charge dans la phase gazeuse (air). Cela permet donc d'estimer le gradient de pression associé à ce petit volume: $\nabla P_0=\frac{P_s-P_e}{L_0}$

Grâce à une formule de type loi de Poiseuille, on calcule le débit induit par ce gradient de pression: $Q_0 = - K*\nabla P_0$

On suppose que ce débit va rester constant pendant un temps $dt$ très petit devant le temps final. On peut donc calculer le nouveau volume de colle dans le domaine: $V_1=V_0 + Q_0*dt$

On obtient ainsi ue nouvelle estimation du gradient de pression qui donne un nouveau débit $Q_1 = - K*\nabla P_1$ qui permet de trouver le volume suivant $V_2=V_1 + Q_1*dt$ et ainsi de suite ...

De manière générale on applique donc l'algorithme suivant à l'instant n:

--> $t^n = t^{n-1} + dt$

--> Estimation de $\nabla P_n$

--> $Q_n = - K*\nabla P_n$

--> $V_{n+1}=V_n + Q_n*dt$

Lorsque $V_n=V_{max}$, on arrête le calcul et on regarde le temps obtenu.

 

Dans cette méthode, il y a deux difficultés à surmonter pour bien estimer le temps final:

- la détermination du coefficient K

- l'estimation de $\nabla P_n$

Les parties suivantes se penchent sur la résolution de ces problèmes.

Détermination du coefficient K

Pour déterminer le coefficient K, nous avons comparé un écoulement plan diphasique dans un à un écoulement monophasique dans le même plan.

En effet, la relation entre le débit et le gradient de pression pour un écoulement monophasique dans un plan est connue théoriquement, c'est un écoulement de Poiseuille. Donc pour calculer de manière théorique le $K_{diphasique}$, il faudrait trouver une relation entre le coefficient diphasique et le coefficient monophasique de la forme: $K_{diphasique} = f( K_{monophasique})$

Quatre paramètres critiques ont été retenus pour l'étude:

- l'angle de contact $\theta$ (en radians ou degrés)

- la tension de surface $\sigma$ (en $N.m^{-1}$)

- la viscosité dynamique $\mu$ (en Pa.s)

- l'épaisseur de la conduite $e$ (en m)

La relation entre les coefficients est donc choisie de manière à avoir la forme suivante :$\frac{K_{diphasique}}{K_{monophasique}} = f(\theta ,\mu ,e ,\sigma)$

Pour identifier la fonction f, on fait varier ces différents paramètres autour d'un cas de référence dont les propriétés sont les suivantes:

- l'angle de contact $\theta_{ref}=30°$

- la tension de surface $\sigma_{ref}=0.072 N.m^{-1}$

- la viscosité $\mu_{ref}=10^{-3}  Pa.s$

- l'épaisseur du plan $e_{ref}=1mm$

 

Détermination du rapport des coefficients:

On lance sous Fluent une simulation d'un écoulement plan (les paramètres sont décrits dans les parties suivantes). Ensuite, on relève au cours de la simulation:

- la pression le long de l'axe central

- le débit massique à l'entrée

- la position de l'isosurface pour laquelle la fraction volumique vaut 0,5 (elle est représentative de l'interface entre le liquide et l'air).
 

 

Lorsque la simulation est terminée, on calcule le gradient de pression dans la phase liquide à l'aide des données de pression. Enfin, on trace le débit en fonction du gradient de pression et on relève le coefficient directeur de la droite obtenue.

 

Il faut toutefois faire attention aux points que l'on utilise pour trouver le coefficient directeur. En effet, il faut exclure les points qui correspondent à l'établissement d'un régime permanent car la relation linéaire est valable seulement quand l'écoulement est établi. En principe, ce sont les points correspondant aux forts gradients de pression car c'est pendant l'établissement que les gradients sont les plus forts.

On calcule théoriquement le coefficient de Poiseuille: $K_{monophasique} = K_{Poiseuille} = \frac{e^3}{12 \nu}$

Enfin, on fait le rapport  de ces coefficients $r = \frac{K_{diphasique}}{K_{monophasique}} $

On répète l'opération en faisant varier un paramètre $p$ puis on trace $r$ en fonction du paramètre que l'on a fait changer. On obtient ainsi une fonction qui nous permet de prédire le rapport des coefficient en fonction de la valeur du paramètre:

$r = f_p (p)$

 

Simulation plan

Pour faire le moins de calculs possible et donc gagner du temps, les simulations sont réalisées avec l'hypothèse de symétrique au niveau de l'axe. C'est-à-dire que nous allons simuler l'écoulement seulement une moitié du plan. Ensuite par symétrie, le plan est reconstitué.

Hypothèse des simulations:

Le $\Delta P$ entre l'entrée et la sortie de la conduite est de 373 Pa. Il a été choisi de manière à avoir des gradients de pression similaires à ceux obtenus avec la colle.

De manière à ce que l'écoulement s'établisse correctement, un rapport de 100 est toujours conservé entre la longueur et l'épaisseur: $\frac{L}{e} = 100$

Les simulations ont été réalisées sous Fluent avec les paramètres suivants:

- Écoulement 2D transitoire

- Écoulement laminaire

- Formulation VOF implicite pour les phases: un liquide et un gaz

- Option "Wall adhésion" activée pour pouvoir modifier l'angle de contact

- Entrée = "Pressure Inlet"

- Sortie =  "Pressure Outlet"

- Mur = "Wall"

- Axe = "Symetry"

 

Dans les parties suivantes se trouvent les résultats des simulations organisés de la manières suivante:

- Le paramètre qui a changé et ses différentes valeurs.

- Le déplacement de l'interface dans le plan en fonction du temps pour les différentes valeurs. Cela permet de se faire une idée tangible de l'influence du paramètre. La position de l'interface est adimensionnalisée par l'épaisseur et le temps par $t_c = \frac{e \mu}{\sigma}$ (analyse dimensionnelle).

- La fonction $r=f(p)$ qui ressort des simulations.

 

Attention: lorsqu'un paramètre varie, les autres prennent les valeurs de référence notées dans la partie Détermination du coefficient K.

Influence de l'angle de contact

Les simulations ont été réalisées avec les angles de contact suivants:

  • 150° (cas de référence)
  • 120 °
  • 90°
  • 60°
  • 30°

 

 

 On remarque que l'écoulement est plus rapide pour des angles de contact élevées. Cela est très utile car un fluide mouillant (c'est-à-dire avec un fort angle de contact) aura un double avantage: il évite de piéger des bulles et il est celui qui avance le plus rapidement.

Détermination de la relation de prédiction de comportement de fluide simulé:

 

D'après la disposition des points une régression quadratique semble appropriée. Voici la fonction que l'on obtient:

$f_\theta (\theta) = 0,015964 * \theta ^2 - 0,49758* \theta + 1,4421$

L'angle $\theta$ est exprimé en radians.

 

 

Influence de l'épaisseur d'âme

Les simulations ont été réalisées avec les épaisseurs suivantes:

  • 0.1 mm
  • 0.5 mm
  • 1 mm (cas de référence)
  • 5 mm
  • 10 mm

On peut voir que plus l'espace entre les plaques est petit, plus le liquide avance lentement. Un faible écartement entre les plaques implique donc des temps de remplissage longs. Toutefois, en augmentant le gradient de pression on peut réduire ces temps. A ce moment là, il faut faire attention à l'effet du gradient sur la planéité des plaques, chose qui n'a pas été étudié ici.

Détermination de la relation de prédiction de comportement de fluide simulé:

D'après la disposition des points une régression linéaire semble appropriée. Voici la fonction que l'on obtient:

$f_e (e) = 356,68* e + 1,0074 $

Influence de la tension de surface

Les simulations ont été réalisées avec les tensions de surfaces suivantes:

  • 20 mN/m
  • 44 mN/m
  • 72 mN/m (cas de référence)
  • 100 mN/m
  • 400 mN/m

A partir de ces résultats, on remarque que l'interface liquide/air avance très rapidement lorsque la tension interfaciale est faible. Cela peut se comprendre car plus la tension de surface est forte, plus l'interface aura tendance à "résister" face à l'avancement du liquide Cela peut éventuellement entraîner des déformations des plaques mais seulement si celles-ci sont très fragiles.

Détermination de la relation de prédiction de comportement de fluide simulé:

D'après la disposition des points une régression linéaire semble appropriée. Voici la fonction que l'on obtient:

$f_{\sigma} (\sigma) = 1,689 * \sigma + 1,4771 $

 

Influence de la viscosité

Ce paramètre n'a pas été étudié par manque de temps. Les autres paramètres étaient prioritaires car nous n'avons qu'un seul résultat expérimental avec une viscosité différente. De plus, ayant très peu d'informations sur l'effet de la tension de surface et l'angle de contact, nous avons préféré étudier ces deux derniers paramètres.

Estimation du gradient de pression

Pour trouver le gradient de pression, il est nécessaire de connaître la distance $L_n$ entre l'interface et la buse d'entrée. Cette distance est toujours prise le long de la diagonale qui relie l'entrée et la sortie.

Pour la calculer, on suppose d'abord que la répartition de la colle se fait en trois étapes. Ensuite, on se place dans le repère orthonormé dont l'origine est la buse d'entrée (coin inférieur gauche). Puis on cherche le point d'intersection entre la diagonale entrée-sortie et la diagonale qui forme l'interface en utilisant les équations de ces diagonales. La longueur $L_n$ est donnée par la distance entre le point d'intersection et l'origine.

Première étape: le volume de colle forme un triangle qui avance. Elle est valable jusqu'à ce que le triangle touche le coin supérieur gauche de la configuration.

Deuxième étape: le volume de colle est un trapèze qui avance. Elle est valable jusqu'à ce que la diagonale touche le coin inférieur droit de la configuration.

Troisième étape: le volume d'air est un triangle qui recule. Elle est valable jusqu'à la fin.

 

Le $\Delta P$ utilisé pour définir le gradient a d'abord été défini comme $\Delta P = P_s - P_e$. Toutefois, son calcul a été amélioré en prenant en compte le saut de pression du à l'interface:

$\Delta P = P_i - P_e$ avec $ P_i = \frac{2\sigma cos \theta}{e} + P_s$

Le gradient de pression est donc estimé à l'instant n par:

$\nabla P_n = \frac{P_i - P_e}{L_n}$

Résultats

Voici les résultats qui ont été obtenus avec la méthode précédente, ils sont comparés dans le tableau suivant avec les temps fournis par l'entreprise.

La configuration 1 correspond à des plaques de dimensions 21cmx30cm écartées de 1mm et la configuration 2 à des plaques de 30cmx42cm écartées de 1mm.

Les caractéristiques des fluides sont $\rho_1 = \rho_2 = 980 kg.m^{-3}$ et $ \mu_1 = 0.75 Pa.s, \mu_2 = 1 Pa.s $

Résultats de calcul du temps de remplissage

Configuration Fluide Diamètre de la buse(mm) Résultats expérimentaux (s) Résultats de simulation (s) Erreur(%)
1 1 0,69 3000 2500 16,66
1 1 0,84 600 -720 1083 50
1 2 1,19 1800 ? ?
2 1 0,84 3000 3057 1,9
2 1 1,19 1020 1050 2,94
2 1 1,60 720 738 2,5

 

Pour la deuxième ligne du tableau deux expériences ont été réalisées, l'écart entre les temps est d'environ 20%. Nous avons donc essayé d'obtenir une précision inférieure à ce seuil.

On peut remarquer que les résultats avec la grande configuration sont bien meilleurs que les autres. Cela peut venir du fait que les gradients mis en jeu sont plus forts pour la configuration 1. Or, de forts gradients de pression impliquent qu'une petite erreur sur K sera plus néfaste sur la précision qu'avec de faibles gradients.

Le fluide 2 n'a pas été étudié car nous n'avons pas pu effectué une étude de l'influence de la viscosité.

Annexe

Une approche théorique a été réalisé au niveau de l'écoulement plan diphasique. Elle ne permet pas de d'apporter de résultats supplémentaires mais permet de bien décrire l'avancée de l'interface dans le temps. De plus, elle valide (au coefficient près) le temps choisi pour adimensionnaliser qui avait été trouvé par analyse dimensionnelle. Les courbes ont d'ailleurs été reprises avec le temps caractéristique trouvé dans cette équation.

 

La comparaison avec les simulations est plutôt bonne au vue des hypothèses simplificatrices qui ont été réalisées.

Conclusion

Ce projet a été très intéressant pour nous car il nous a montré qu'il n'est pas toujours nécessaire d'employer de grands moyens pour arriver à des résultats concrets. En effet, une méthode simple nous a permis de remplir les objectifs principaux sans faire de lourdes simulations.

Toutefois, même si les résultats sont plutôt bons, il existe des perspectives d'améliorations. Voici les points sur lesquels il y a eu peu d'étude:

- la perte de charge singulière en entrée et en sortie du domaine, notamment quand le diamètre des buses est plus grand que l'épaisseur d'âme.

- la prise en compte des pertes de charge linéaires dans la phase gazeuse.

- l'estimation du gradient de pression peut être améliorée en affinant les étapes de répartition de la colle.

- les fonctions de correction peuvent être affinées en réalisant plus de simulations et en prenant notamment en compte l'effet de la viscosité.

- une possible utilisation de l'analyse théorique pour calculer le débit.