II. Ecoulement cisaillé : Approche Surfacique

Dorénavant, l’écoulement sera considéré comme cisaillé. On cherchera alors à prendre en compte les effets hydrodynamiques sur le coefficient de transfert de masse. L'objectif est de mettre en avant le fait que le coefficient de transfert de masse subit une influence de l’écoulement cisaillé et donc de l’hydrodynamique.

Dans cette approche, la réaction est considérée surfacique, c'est à dire qu'elle ne se déroule qu'à l'interface laitier/acier. C'est la diffusion de l'aluminium vers l'interface qui va alimenter la réaction.

Les concentrations des espèces dépendent maintenant de y et t mais également de x. Cette partie sera plus lourde en termes de calculs dans la mesure où l’étude se fera en 2D, sur Fluent.

 

Mis en place de l'écoulement cisaillé

Les écoulements d'aciers dans la lingotière sont soumis à des cisaillements importants. Pour notre étude il est donc nécessaire d'arriver à reproduire ces cisaillements, en les simplifiant, pour réussir à modéliser les phénomênes qui se déroulent dans le procédé. 

L'objectif est ici de créer un écoulement d'acier cisaillé dans tout le domaine étudié. Sous Fluent cela se traduit par la mise en place d'une "user-defined function" (UDF) basée sur les vitesses à l'interface acier/laitier et à la limite inférieure. Cet outil va nous permettre ici d'implémenter en entrée du domaine, une fonction linéaire de la vitesse dépendante uniquement de la distance à l'interface. La linéarité de cette fonction est la conséquence des hypothèses de stationnarité et de fluide newtonien, pour notre système.

La pente de l'UDF est calculé grâce au rapport de viscosité entre les deux phases et grâce aux données d'ArcelorMittal.

Formule pour l'injection du profil de vitesse linéaire en entrée du domaine:

$$v(y)=\frac{v_{interface}-v_{bottom}}{h}\cdot y+v_{bottom}$$

 

Après avoir mis en place l'UDF, il est nécessaire d'imposer des conditions limites pour les cotés supérieur et inférieur de la géométrie. En effet, pour que les conditions limites soient en accord avec la fonction de l'UDF on considèrera les limites supérieur et inférieur du système comme des "moving walls" ou l'on imposera successivement la vitesse à l'interface (côté supérieur) et la vitesse de la limite inférieure (côté inférieur). Ces valeurs aux limites doivent être recalculées à chaque fois que la géométrie du système change.

 

Figure 12 : Profil de la vitesse imposée par l'UDF en entrée d'un domaine d'une hauteur de 4cm

 

 

Création du mélange acier sous Fluent

Pour entamer la modélisation, il est nécessaire de bien définir le matériau "acier".

L'acier est un alliage métallique principalement composé de fer, mais aussi d'un certain pourcentages d'autres composants tel que l'aluminium, le magnésium, le manganèse etc... Dans un soucis de simplification, nous considèrerons d'abord que l'acier est composé uniquement de fer et d'aluminium. Ces deux composants, sous forme liquides, n'étant pas répertoriés dans la base de donnée Fluent, il est nécessaire de les créer et de les paramétrer en fonction de leurs propriétés physiques.

Dans un deuxième temps, il faut paramétrer les propriétés du mélange à partir des deux composants fer et aluminium. On choisit les lois "weighted-mixing-law"pour la densité et la viscosité, qui réalise une moyenne pondérée des propriétés des deux composant pour donner celle du mélange. Pour la diffusivité, on choisit une loi "dilute-approx" car l'aluminium a une fraction massique inférieure à 1% dans l'acier.

D'après les données d'ArcelorMittal, la fraction molaire de l'aluminium dans l'acier est de 32e-5, qu'il faudra renseigner à l'inlet pour obtenir le mélange souahité.

 

Tableau 1 : Tableau récapitulatif des paramètres des espèces de l'acier simplifié

Simulations Laminaires

L'écoulement en situation réel est turbulent. Cependant, en première approche, nous allons réaliser une étude du système en considérant l'écoulement laminaire. Cette étude préalable va nous permettre de valider les premières étapes des simulation Fluent pour un cas simple, en nous appuyant sur des résultats fournis par l'entreprise.

A la suite de la résolution numérique, nous pourrons analyser les profils de concentration des composants de l'acier dans l'écoulement. Ces profils de concentration devront mettre en évidence une couche de diffusion correspondant aux transfert de matière entre l'interface et l'acier.

 

 

Mise en place des simulations

Paramétrage du modèle de résolution et choix du solveur

Le modèle de résolution pour notre système doit prendre en compte l'écoulement laminaire et le transport de matière.

L'écoulement laminaire doit être résolu ici par le "Viscous Model" : "Laminar".

Pour la résolution des équations du transfert de matière, on choisi dans "Species Model", le modèle "Species Transport" avec le fluide "acier" que l'on crée. Dans le cas d'une réaction surfacique, il faut décocher l'onglet "Volumetric" dans la case "Reaction". De plus, on choisit l'option "Diffusion energy source". 

Le solveur utilisé en première approche est le solveur "SIMPLE" avec des discrétisation du premier ordre. Ces dispositions sont nécessaires pour minimiser les temps de calculs et n'influent pas sur les résultats.

 

Choix de la géométrie et du maillage

Le choix de la géométrie est basé sur les résultats du rapport ArcelorMittal. En effet, le précédent rapport se base sur une géométrie rectangulaire 2D de 4mm de hauteur. Après les premiers tests nous avons décidé de réduire cette hauteur afin d'avoir une meilleure résolution pour les phénomènes de diffusion dont la couche limite ne dépassait pas 2mm. Une longueur de 0,4m a été basé la longueur de d'établissement du profil de vitesse pour une résolution laminaire.

Le maillage choisi a pour but de bien comprendre les phénomènes à l'interface laitier/acier. Ainsi on choisit un maillage de type "geomteric2" avec un ratio de 1.05, pour l'axe vertical. Le système étant en régime stationnaire, la résolution selon l'axe x est peu élevée et "uniform".

On a finalement un domaine étudié pour bien analyser les phénomènes de diffusion à l'interface.

Figure 13 : Schéma de la géométrie

 

Figure 14 : Maillage d'étude pour la résolution laminaire

 

Définition des conditions limites

Le système est composé de quatre segments limites qu'il est nécessaire de paramètrer avant la simulation.

 

Inlet : Velocity inlet

- vitesse : définie par l'udf crée précedement

- species : d'après les données de l'entreprise, l'acier est composé d'une fraction massique d'aluminium de Yal = 32.10-5

 

Outlet : Pressure outlet

On impose la sortie du domaine en "pressure outlet" en paramètrant de manière à ne pas avoir de gradient de pression entre l'entrée et la sortie, ce qui correspond à un écoulement cisaillé.

 

Interface : Moving wall

- vitesse : vinterface= 0,1125 m/s, calculé à partir du rapport de viscosité entre le laitier et l'acier

- species : Yal = 5.10-5, qui correspond à la valeur de l'équilibre chimique de la réaction entre l'aluminium et le SiO2 du laitier

 

Bottom : Moving wall

- vitesse : vbottom = 0,206 m/s, calculé à partir de la fonction implémentée dans l'udf et vinterface

-species : Yal = 32.10-5, qui correspond à l'état de l'acier loin de l'interface

Résultats et Validation

Résultats des simulations laminaires

Les résultats des simulations sous Fluent nous donnent les profils de concentrationd'aluminuim selon la distance à l'interface. Grâce à ces profils, on peut calculer le coefficient de transfert de matière à l'interface.

 

Figure 15 : Courbe de l'évolution de la concentration en fonction de la distance à l'interface pour 4 valeurs de x

Les résultats sont confirmé par la faible monté en pression dans le domaine, ce qui correspond à un écoulement de Couette. De plus, les profils de vitesse en entrée et en sortit sont les même, on a donc réussi à imposer le cisaillement dans tout le domaine.

Figure 16 : Mise en évidence de la faible monté en pression dans le système

Cette légère augmentation est expliqué par le fait que Fluent se base ici sur un modèle de résolution par la pression et impose alors ce gradient. Ce faible gradient n'influence pas les résultats.En effet, le profil de vitesse reste constant selon l'axe horizontale.

Figure 17 : Analyse des profils de vitesse en entrée et en sortie du domaine

Validation des résultats

L'objectif du BEI est principalement de pouvoir donner des ordres de grandeur du coefficient de transfert de matière. Ainsi, grâce aux profils de concentration, on va pouvoir le calculer par la formule suivante :

$$k_{l} = \frac{D\frac{(\partial Y_{i}}{\partial y})_{y = 0}}{Y_{i,s}-Y_{i,\infty }}$$

Ainsi, en utilisant les valeurs fournies par Fluent, pour les concentration proche de l'interface, on peut calculer le coefficient Kl.

Kl = 1,6.10-5 m/s

Cette valeur de Kl doit être validé par la théorie. C'est pourquoi nous avons mis en place un programme Matlab qui détermine valeur de Kl à partir de la théorie d'un film tombant dans Bird et al3. En effet, à l'interface, cette théorie présente un cisaillement de type linéaire et développe une résolution analytique théorique que l'on a implémenté sous Matlab.

Matlab : Kl = 4.10-6m/s

La différence entre les deux valeurs s'explique par la présence du cisaillement dans le cas laminaire 2D.

Simulations Turbulentes

L'écoulement réel est de nature turbulent, c'est pourquoi nous souhaitons modéliser le système avec un modèle Fluent, traitant des modèles turbulents. En effet le modèle laminaire ne peut pas prendre en compte les effets de la turbulence dans l'écoulement qu'il calcul. Nous souhaitons ici savoir qu'elle sont les différences de résultats fournis par le modèle laminaire et un modèle de turbulence.

 

Mise en place des simulations

Choix du modèle de turbulence

Fluent permet de modéliser des éoculements fluides turbulents. Pour cela il utilise des modèles qu'il faut savoir paramétrer, et qui dépendent du système étudier. Chacun des modèle se base sur des simplifications bien particulières des équation de Navier-Stokes, dans le but de les résoudre. Il est important de bien choisir le modèle, en fonction du système étudié, car tous les modèle ne sont pas pertinents pour tous les type d'écoulements.

Notre système traite d'un écoulement entre deux parois fictives, à faible nombre de Reynolds. Le modèle K-Epsilon est un modèle adéquate pour ce type d'écoulement car nous n'avons pas de gradient de pression trop important, c'est pourquoi nous avons choisi de l'utiliser en premier lieu. 

Lors de la simulation pour keps Standard les résultat présente Yplus<1. Ainsi, le reste des simulations est réalisé avec l'option "Enhanced Wall Treatment" pour mieux rendre compte de ce qui se passe dans le système. Cette option permet d'avoir un meilleur traitement des calculs dans la sous couche visqueuse, proche paroi.

 

Figure 18 : Mise en évidence de la nécessité du choix du modèle "Enhanced Wall Treatment" pour les calculs à l'interface

 

Choix de la géométrie et du maillage

Après plusieurs tests avec le précédent maillage, il s'est avéré que la turbulence ne pouvait être entretenue pour une hauteur de 2mm. Il a donc fallu changer de maillage pour l'étude turbulente.

Le choix de la géométrie pour le turbulent, a été assez difficile a effectuer et résulte de nombreux tests. En effet la géométrie va influer notamment sur le Reynolds, donc sur la turbulence et donc sur le transfert de matière.

La géométrie finale est un rectangle de hauteur h = 4cm et de longueur L = 1m. La hauteur de 4 cm a été déterminée pour permettre un écoulement turbulent avec les paramètres de l'acier et le profil de vitesse imposé. La longueur a été choisie de manière à obtenir un profil de vitesse établi en fin de domaine.

Le maillage est ici raffiné selon y à l'interface, dans une zone de 3 mm. Pour le reste du domaine, il est choisi uniforme selon x et y.

 

Figure 19 : Zoom du maillage utilisé pour les simulations turbulentes

 

Définition des conditions limites

Pour cette étude, on reprend les mêmes conditions limites que pour le laminaire en terme de concentration et de vitesse, excepté pour le "bottom".

En effet, comme la taille du domaine change et que l'on souhaite garder le même cisaillement que dans le laminaire, on doit recalculer vbottom. Pour cela on utilise la pente précédente et la taille du nouveau domaine. Ainsi vbottom = 1,987 m/s.

 

Calcul de l'énergie cinétique turbulente du système

Dans l'étude turbulente, il est indispensable de connaître les valeur de l'énergie cinétique turbulente(kt) et de la dissipation visqueuse (epsilon) de manière à les renseigner dans les paramétrages du modèles de calcul.

En première approche, on utilise l'hypothèse de Boussinesq pour calculer ces deux paramètres à partir des vitesses moyennes du système.

kt ~ (0.05vm)2 ~ 2.10-3 m2/s2

Eps ~ 3.kt1.5/h ~ 9.10-3 m2/s3

Cependant ces valeurs ne correspondent pas vraiment au contexte de la coulée de l'acier. Pour tenter de calculer une énergie cinétique turbulente en lien avec la géométrie de la lingotière, nous avons réalisé l'étude suivante. Le but est ici de déterminer, grâce à une géométrie approchant celle de la lingotière, la valeur du kt entrant dans le système considéré, proche de l'interface qui nous intéresse.

 

Figure 20 : Géométrie des simulations de détermination de l'énergie cinétique en entrée du domaine turbulent

La simulation qui est faite pour ce maillage permet de déterminer l'ordre de grandeur réel de l'énergie turbulente dans le domaine qui nous intéresse (cadre rouge). Les résultats pour les différentes profondeurs du domaines sont fournis ici:

Figure 21 : Evolution de l'énergie cinétique turbulente en fonction de x pour différentes profondeurs y

 

Les résultats indiquent que l'ordre de grandeur de kt et de epsilon sont les suivants :

kt ~ 10-4 m2/s2

Epsilon ~ 7,5.10-4 m2/s3 (avec l'équation de Boussinesq)

On décide ainsi d'injecter ces valeurs dans la simulation Fluent pour paramétrer le modèle de turbulence.

 

Résultats et Validation

 

Résultat de l'étude surfacique pour un écoulement turbulent

Le but du BEI est d'analyser les profils de concentration à l'interface laitier/acier et d'en déduire les coefficients de transfert de matière entre les deux phases. Nous avons simplifié l'étude à l'analyse des concentrations d'aluminium pour en déduire une technique de simulation de la diffusion de celui ci dans le mélange acier. 

La réalité nous dit qu'un écoulement turbulent favorise les transferts de matières, ainsi les résultats attendues doivent donner des coefficient Kl turb>Kl lam.

Le profil suivant est fournis pour une simulations Fluent sous un modèle KEpsilon. Il représente l'évolution de la concentration d'aluminium en fonction de la distance à l'interface laitier/acier, à la sortie du domaine.

 

 

Figure 22 : Graphique des profils de concentration pour une étude laminaire et pour une étude turbulente

 

On souhaite comparer les valeurs des coefficients de transfert de matière entre le cas laminaire et turbulent.

Avec les données Fluent relatant l'évolution de la concentration d'aluminium, on détermine le rapport dY/dy à la paroi, et on en déduit Kl avec la formule suivante :

$$k_{l} = \frac{D\frac{(\partial Y_{i}}{\partial y})_{y = 0}}{Y_{i,s}-Y_{i,\infty }}$$

On obtient ainsi Kl= 4,3.10-5 m/s

Avec Kl lam = 1,6.10-5 m/s, on retrouve bien la condition suivante : Kl turb>Kl lam

 

De plus, les profils turbulents présentent une plus petite couche de diffusion que celle du laminaire. Ceci s'explique par l'apport infini en aluminium, provenant du "bottom". En effet, selon nos conditions limites, il n'y a pas épuisement de l'espèce, or avec les turbulences, le renouvellement de l'espèce consommée est pratiquement instantané. C'est ce qui explique la différence entre le laminaire et le turbulent, en terme de couche de diffusion.
 
Une des limites de cette simulation se trouve dans l'analyse des profils de vitesse à la sortie du domaine.

 

Figure 23 : Mise en évidence de la décroissance du cisaillement dans le domaine

 
En entrée, l'udf nous permet d'imposer un cisaillement à l'écoulement. Cependant, le profil de vitesse est modifié par la turbulence du système. Ainsi, en sortie du domaine, le profil de vitesse ne présente pas le même cisaillement qu'en entrée. Le changement d'allure des profils de vitesse, par l'injection de turbulence, reste un problème à résoudre car il influence les résultats du transfert de matière.
 

Validation des résultats par l'article de Calmet et al 1:

La validation de nos simulation se fait par la comparaison des coefficients de transfert calculés à la littérature. Pour cela nous nous basons sur les conclusions d'un article (ref) qui présente les résultats suivants.

 

Tableau 2 : Tableau récapitulatif des résultats de l'article de Calmet et al1 concernant les écoulements libre et cisaillé.

 

Les expériences de cet articles se base sur des écoulements à Sc=200 et à Re*=640. Le choix du maillage a été fait pour que les Reynolds de frottement aux parois (Re*) soient du même ordre de grandeur, en entrée, pour nos simulations, pour l'article de Calmet et al1 et pour le rapport fournit par ArcelorMittal. Ainsi les résultats pourront être comparés facilement.

 

Tableau 3 : Tableau de comparaison des nombres adimensionnels utilisés dans les différents travaux

 

Seule la deuxième colonne du tableau nous intéresse car elle traite d'une diffusion de matière dans un écoulement cisaillé, entre deux parois. L'objectif est ici de retrouver des valeur du même ordre de grandeur pour le ratio Kl/u*.

On a déjà Kl= 4,3.10-5 m/s

Il nous faut maintenant calculer le u*, vitesse de frottement à la paroi, dont la définition est la suivante :

u*=μ.(du/dy)0.5

 

Cependant la méthode pour calculer u* reste assez discutable car notre profil de vitesse en sortie n'est pas comparable à celui de l'article. Ainsi nous avons calculé cette vitesse de trois manières différentes : 

Méthode 1 : Ici on utilise la valeur de du/dy en sortie du domaine et à la paroi pour calculer le u* de notre simulation

Méthode 2 : On utilise la valeur de du/dy en entrée du domaine pour avoir une vitesse de frottement qui correspond au cisaillement initial

On obtient ainsi les résultats suivants :

Tableau 4 : Tableau de comparaison pour validation à l'article de Calmet et al1

Étant donné que l'article est basé sur des simulations ou l'écoulement conserve un cisaillement linéaire et constant tout au long du domaine, on peut penser que la méthode 2 est plus approprié. C'est en effet elle qui donne les résultats le plus proche des résultats de l'article.

La différences de valeurs peut être liées à de nombreux facteurs:

- l'espèce qui diffuse : ici aluminium liquide, gaz peu soluble pour Calmet et al1

- la différence entre le nombre de Reynolds de frottement pour nos simulation (Re*=337) et leurs simulations (Re*=640)

Cette méthode de validation est assez peu fiable mais nous permet de voir que l'ordre de grandeur de Kl trouvé est approximativement correct. Il est cependant nécessaire de trouver une meilleure manière de valider les valeurs d'éventuels simulations futures.

Discussion des résultats

kt et epsilons définissent la turbulence du courant en entrée du domaine. Dans le maillage choisi, ils n'ont pas beaucoup d'influence. En effet quelque soit leurs valeurs la turbulence dans le domaine étudié atteindra toujours des valeurs du même ordre de grandeur, autour de 0.05. Cette valeur n'est pas en adéquation avec ce qui se déroule dans la lingotière et est la conséquence de la vitesse trop élevée au bottom : vbottom>vinjection

Figure 24 : Evolution de kt dans le domaine de 4 cm de hauteur

Alternative de simulation turbulente : "fixed values"

La géométrie utilisé précédemment ne permet pas d'obtenir une énergie cinétique turbulente du même ordre de grandeur que celle de la lingotière. Elle est bien trop élevées et cela a peut avoir des conséquences sur les profils de vitesses et sur les profils de concentration.

Dans cette partie, nous étudions l'utilisation sous Fluent de l'option "Fixed Values" dans "Cell Zone Condition". Cette onglet permet de fixer des valeurs dans l'ensemble du domaine, notammment celles de kt et de epsilon. Dans notre cas, les valeurs de kt restent de l'ordre de 10-4 dans tout le domaine, d'après l'étude sur la lingotière complète. En fixant kt et la dissipation turbulente associée, on élimine le problème de surproduction d'énergie cinétique turbulente dans le cas précédent.

Figure 25 : L'option "Fixedvalues" sous Fluent

Résultats et comparaisons des deux approches

Figure 26 : Profils de vitesses pour les deux approches de simulation turbulente

 

La comparaison des profils de vitesse nous indique que dans le premier cas, c'est la forte turbulente qui a entraîné la dérive du profil de vitesse. Pour le cas "fixed values", on voit que cette dérive est moins rapide et permet de conserver un profil de vitesse similaire à l'UDF d'entrée, sur une plus longue distance. La comparaison avec l'article de Calmet & Magnaudet sera plus pertinente dans ce cas là car le cisaillement reste le même en entrée et en sortie. 

 

Figure 27 : Profils de concentrations pour les deux approches de simulation turbulente

 

Les deux profils de concentration restent sensiblement les mêmes pour les deux approches. Cependant on voit qu'il y a une différence d'épaisseur de la couhe limite de diffusion. En effet, dans le cas "Fixed value", la turbulence reste assez faible pour que cette couche soit visible sur le graphique, contrairement à la première approche.

 

Discussions

Comparaison des valeurs à l'article de Calmet et al1

Pour cette deuxième approche, on calcule les valeurs de u* et de Kl. La encore, nous comparerons les valeurs données pour les trois méthodes de calcul de u*.

Kl = 2,37.10-5 m/s

Tableau 5 : Tableaux des valeurs nécessaire à la validation por les tests en Fixed et Unfixed values

Tout d'abord, on a Kl unfixed > Kl fixed ce qui s'explique par la faible turbulence dans le deuxième cas, qui limite le transfert de matière.

Concernant les valeurs du rapport Kl/u*, on voit d'abord que les résultats sont du même ordre de grandeur dans les deux approches pour les deux méthodes. La méthode 2 reste la plus proche dans les deux approches de la turbulence, "Fixed" et "Unfixed". L'ordre de grandeur pour le ratio est un peu faible pour le cas "Fixed value", à cause de la faible valeur de turbulence dans le domaine.

Comparaison des approches "Fixed Values" et "Unfixed Values"

L'approche "Fixed value" a pour avantage de mieux retranscrire les comportement du système pour l'énergie cinétique observé dans la lingotière. En effet, dans ce cas, la valeur de kt est controlée et ne dépendant pas de la la géométrie pour la simulation. Cependant le processus de résolution se simplifie et peut  entraîner des dérives numériques. En effet, l'équation de transport de kt n'est pas résolue dans ce cas. Le véritable problème est l'absence de réalité physique que comporte le fait de bloquer la valeur de kt, ce qui peut influer sur la résolution numérique.

Il est dans l'interêt de l'industriel de connaître ces deux approches, leurs avantages et leurs inconvénients, dans l'optique de pouvoir privilégier le controle d'un paramètre plutôt qu'un autre selon les simulations souhaitées.