Calcul numérique le long du domaine (transition laminaire/turbulent)

Pour ne pas imposer un quelconque profil expérimental en entrée qui pourrait modifier l'écoulement, nous choisissons de réaliser un calcul dans le domaine en entier, c'est à dire de calculer la couche limite depuis le bord d'attaque, en passant par la zone de transition du régime laminaire vers le régime turbulent, jusqu'au bord de fuite (fin de la plaque).

Pour ce, nous procédons de deux façons:

a) Un seul calcul dans tout le domaine en considérant cette fois, un entrée libre.

b) Connaissant le point de transition laminaire/turbulent en x=0.8 m, nous effectuons deux calculs : Le premier va concerner la zone où la couche limite est laminaire, et ce en annulant la viscosité turbulente dans tout le domaine, ensuite nous passerons à un second calcul qui va concerner uniquement le domaine où la couche limite est turbulente, avec comme conditions en entrée les résultats en sortie du premier calcul. De cette façon, nous pourrons forcer la transition.

Mise en place des calculs

Calcul en un seul bloc

Le calcul se fera comme précédemment avec le solveur buoyantBoussinesqSimpleFoam, avec cette fois-ci une modification de la géométrie ( l'entrée ne commence plus à x=1.44m mais plutôt à 0). Le solveur potentialFoam n'est plus utilisé puisqu'il n'y a plus de problème de discontinuité des conditions de vitesses, et les conditions aux limites sont résumés dans le tableau ci-dessous :

  Entrée Plaque Sortie Bord
U (m/s) nonuniform List<vector>

fixedValue uniform (0 0 0)

zeroGradient zeroGradient
p (Pa) zeroGradient zeroGradient

fixedValue uniform

fixedValue uniform
T (K) nonuniform List<scalar>

fixedValue uniform 333.15K

zeroGradient zeroGradient

 

Calculs en deux blocs

Dans le premier domaine (de x=0 à x=0.8 m), nous lançons un premier calcul sous buoyantBoussinesqSimpleFoam avec cette fois-ci une viscosité turbulente nulle pour forcer le régime laminaire, et les conditions aux limites seront exactement identiques à celles utilisées précedemment.

Une fois le calcul convergé, nous récupérons les résultats au dernier pas de temps en x= 0.8m, pour ensuite les imposer comme conditions d'entrée pour le second domaine/calcul.

La récupération de ces profils en entrée peut se faire en passant d'abord par le logiciel "pararaview", pour ensuite enregistrer les datas sous format texte (data.dat), et enfin faire une des opérations de tri via quelques commandes linux :

  • Afin d'extraire les lignes contenant uniquement les données relatifs à la position x=0.8 et z=0.001

> grep -w 0.8 data.dat | grep -w 0.001  > Y.dat

  • Pour extraire les composantes des vitesses (colonnes 3, 4 et 5)

> cut -d',' -s -f3,4,5 Y.dat > U.dat

  • Pour extraire le profil de température (colonne 2)

> cut -d',' -s -f2 Y.dat > T.dat

  • Pour extraire le profil des pressions (colonne 1)

> cut -d',' -s -f1 Y.dat > p.dat

  • Pour respecter la forme avec laquelle s'écrit un profil de vitesses de type nonuniform List<vector>

> sed 's/^/(/' U.dat > U.dat (ouvrir une parenthèse au début de chaque ligne)

> sed 's/$/)/' U.dat > U.dat (fermer les parenthèses à la fin de chaque ligne)

> sed -e 's/,/ /g' U.dat > U.dat (remplacer le virgules par des espaces)

Après avoir récupérer les profils de vitesse, pression et température en sortie du premier bloc, on passe au calcul dans le deuxième bloc où les conditions en entrée seront de type nonuniform. Ainsi, il va falloir copier les dernier profils de vitesse, température et pression, lancer un premier calcul sous potentialFoam pour résoudre les problèmes de discontinuité, et ensuite réaliser un dernier calcul via buoyantBoussinesSimpleFoam.

 

 

 

 

Étude de la couche limite laminaire

Étude théorique

Dans cette partie, nous allons exploiter les simulations qui ont été faites sur le bloc en entier, en particulier la zone de couche laminaire, pour deux raisons. La première est de s'assurer que notre simulation est bien représentative de la réalité en comparant les résultats numériques à  la théorie, la deuxième raison est de s'assurer que les profils expérimentaux en entrée fournis sur le site d'ercoftac collent bien avec la simulation.

Tout d'abord, nous allons présenter les résultats théoriques concernant la partie laminaire.

Pour déterminer les profils de température et de vitesse théoriquement dans cette couche, nous allons nous baser bien évidemment sur les équations de Navier Stokes en couche limite, avec le modèle de Prandtl.

Ce modèle se base sur les mêmes hypothèses que nous avons mentionnées précédemment (voir http://hmf.enseeiht.fr/travaux/bei/beiep/content/g18/analyse-theorique-et-etude-preliminaire).

Nous obtenons donc les mêmes équations de conservation sauf le terme turbulent qui disparaît dans ce cas.

Conditions aux limites :

Vitesse :

$ \overrightarrow{U} (x,0)= \overrightarrow{0}$

$ u(x, \infty )=0 $

$ u(0,y)=0 $

 

Température :

T(x,0)=Tw ,

$ T(x,\infty)=Ta $

T(0,y)=Ta 

Donc finalement, notre problème est bien fermé.

Pour le résoudre, l'astuce est d'utiliser des variables de similitude qui permettent de passer des trois équations à dérivées partielles, à deux équations différentielles ordinaires.

 La variable de similitude adéquate pour ce cas est  : $ \eta (x,y)= (\frac {Gr_x} {4})^{1/4} \frac {y}{x} $

 Avec $ Gr_x=\frac{\beta g x^3 (T_w-T_a)}{\alpha^2} $.

Pour simplifier la résolution des équations, au lieu de considérer T, u et v , nous allons travailler avec les grandeurs suivantes :

$ \theta(x,y)=\theta(\eta)=\frac{T-T_\infty}{T_s-T_\infty} $

la fonction de courant $ \phi  $ telle que : $ u=\frac{\partial \phi}{\partial y}  $

 et $  v=\frac{\partial \phi}{\partial x}$

Nous pouvons écrire cette fonction de courant de la forme suivante : $ \phi=4 \nu \frac {Gr_x}{4}^{1/4}  \xi(\eta) $

Donc au lieu de travailler avec des EDP dont les inconnues sont u, v et T, nous allons considérer des EDO dont les inconnues sont $\theta$ et $\phi$ qui ne dépendent que de $\eta$.

Les équations ordinaires à résoudre sont :

 $\frac {d^3 \xi}{d \eta^3}+3 \xi \frac {d^2 \xi}{d \eta^2}-2 (\frac {d \xi}{d \eta})^2+ \theta=0 $

$ \frac {d^2 \theta}{d \eta^2}+3 Pr \xi (\frac {d \theta}{d \eta})=0 $

Pour les conditions aux limites, de la même façon, il faut les écrire en fonctions des nouvelles variables $ \xi $, $ \theta $ et $ \eta $.

Nous avons donc affaire à deux équations différentielles ordinaires (non-linéaires), avec les conditions aux limites qu'il faut, nous pouvons les résoudre numériquement facilement (avec matlab par exemple)..

A partir de cette résolution, nous pouvons accéder à des paramètres importants comme le nombre de Nusselt, l'épaisseur de la couche limite, la vitesse maximale en fonction de x, ..etc.

Le nombre de Nusselt :

Par définition du nombre de Nusselt : $ Nu= - \frac {\frac {\partial T}{\partial y} x }{Tw-Ta}$

Ce qui donne après calcul :  $ Nu= (\frac {Gr_x}{4})^{0.25} g(Pr).$

Avec $g(Pr)$=$ \frac {0.75 Pr^{0.5}}{(0.609+1.221 Pr^{0.5}+1.238 Pr)^{1/4}}$

L'épaisseur de la couche limite :

A partir de la définition de cette grandeur, on peut déduire facilement que cette grandeur varie en $ x^{1/4} $ : 

$ \delta (x) = x (\frac{Gr_x}{4})^{0.25}$

La vitesse maximale :

Une autre grandeur intéressante est la vitesse maximale Umax.

Théoriquement, on trouve que $ Umax \approx 0.5 \sqrt{g \beta (Tw-Ta) x} $

Exploitation des résultats & Comparaison

A présent nous avons les résultats théoriques ainsi que ceux de simulation. Un dernier point à préciser est que pour les différentes grandeurs que nous allons comparer, nous n'avons pas accès directement par Parafoam. Pour chaque grandeur, il y en a des manipulations à a faire.

Le nombre de Nusselt :

Numériquement, nous avons les profils de température en tout point, pour calculer le nusselt, il suffit de calculer le gradient de température au niveau de la paroi. Pour cela, nous avons calculé la différence entre la température de la paroi et celle de la maille qui suit suivant la direction y. Pour que le calcul soit correct, il faut avoir des mailles suffisamment fines.

Remarque Pour tout calcul sur Parafoam, il ne faut surtout pas exporter les résultats pour des lignes données utilisant plotoverline, car ceci ne permet d'avoir que des valeurs pour un pas d'espace de 0.005 quelque soit le maillage, l'une des faiblesses de cet outil. Pour surmonter cette difficulté, il faut exporter les résultats de tout le domaine, et après faire un tri sur les zones qui intéressent l'utilisateur.

La comparaison des deux approches donne :

                                 

Nous remarquons que les deux profils ne collent pas du tout, quelque soit le maillage. L'origine de cet écart ne provient pas de la simulation ni du maillage, car nous avons bien vérifié que la simulation est correcte, les profils expérimentaux à l'entrée de la couche limite turbulente collent bien avec les profils à la sortie du premier bloc. Le problème vient plutôt de la définition du gradient de température que nous avons considéré pour calculer le Nusselt numériquement. En effet, théoriquement, on montre que le Nusselt augmente avec x, ce qui est bien physique. Cependant, Si pour notre définition du gradient, comme la température suivant x augmente avec x, donc le gradient diminue forcément, par suite Le Nusselt diminue avec x. Nous avons beau cherché pour résoudre ce problème, mais malheureusement nous nous sommes pas arrivés à le résoudre. Les idées que nous avons eu pour surmonter cette difficulté étaient de chercher si OpenFOAM permet de calculer automatiquement ces grandeurs, en les codant par exemple (comment ce qu'on peut faire pour les coefficients de portance et de traînée). L'autre idée était de se servir du Python pour les coder, car ce langage s'avère très puissant pour l'automatisation des calculs et l'obtention des résultats. Malheureusement, nous n'avons pas eu le temps pour faire ces manipulations.

Vitesse maximale :

Pour cette grandeur, le calcul est moins délicat. La comparaison des deux approches nous donne :

 

 

                                      

 

Nous constatons que les deux profils ont la même forme avec un "léger" écart.

Epaisseur de la couche limite dynamique :

Concernant cette grandeur, c'était la galère pour la déterminer numériquement. En effet, si nous faisons un plotoverline pour tous les points du maillage, le calcul de cette épaisseur est impossible avec Parafoam, vu que cet outil considère des pas d'espace de 0.005m et non pas les points du maillage, ce qui est quasiment pas faisable, puisque cette épaisseur est de l'ordre du millimètre. Nous avons essayé d'exporter les résultats de tout le domaine, cependant, vu le nombre de points du maillage, et les vitesses qui étaient très sensibles, le calcul de cette grandeur était imprécis.

Exploitation des résultats

Un bloc

Un premier calcul basé sur des grandeurs turbulentes mises en place par défaut (cas test hotroom) nous permet d'obtenir les courbes ci-dessous :

Nous constatons que les profils numériques sont plus ou moins satisfaisants puisqu'ils s'approchent des profils expérimentaux, surtout en x = 2.535 m, où non seulement la simulation numérique prévoit la pente de décroissance mais également  le maximum de vitesse.


Un second calcul avec cette fois-ci des valeurs pour les grandeurs turbulentes calculées à partir des données expérimentales, donne les résultats ci-dessous :

 

Le tracé des courbes de vitesse et température pour cette deuxième simulation permet d'abord de confirmer la sensibilité du calcul aux grandeurs turbulentes, ainsi que de remarquer la non-amélioration des résultats, contrairement à ce que nous pourrions nous attendre.

Cette remarque pourrait être justifiée par la manière dont le calcul des grandeurs turbulentes s'est effectué ou imposé dans le domaine de calcul. En effet, le fait d'imposer une valeur moyenne dans tout le domaine ne peut être physiquement correcte surtout au niveau de l'entrée, d'ailleurs  cette hypothèse peut être confirmée par l'amélioration observée des profils quand on s'éloigne de l'entrée.

Deux blocs

 

La dernière configuration de calcul, à savoir la simulation numérique en deux blocs permet d'obtenir les profils ci-dessous :

Nous constatons que les profils numériques sont d'autant plus satisfaisants que ceux obtenus pour un calcul en un seul bloc, ce qui est parfaitement cohérent avec la théorie/expériences de la transition de la couche limite laminaire vers la couche limite turbulente en x = 0.8 m. 

Le remarque faite précédemment sur l'influence des valeurs des grandeurs turbulentes est revérifiée.

Raffinement du maillage

Afin d'étudier l'influence du maillage sur la précision des résultats du calcul en un seul bloc, nous avons considéré deux maillages différents, un maillage grossier (500*50*1) avec un coefficient d'expansion suivant y de 10, et un autre raffiné (1000*50*1) avec un coefficient de 20 suivant y. Les résultats obtenus, toujours en les comparant avec la référence expérimentale, sont affichés ci-dessous pour les différentes hauteurs 1.438m, 1.918m, 2.535m et 3.244m :

Nous remarquons que les courbes obtenus avec le maillage grossier sont légèrement plus proches des valeurs expérimentales , ce qui est surprenant.  Comme pour le calcul ERCOFTAC, l'hypothèse pour expliquer cette remarque  pourrait se baser sur le coefficient d'expansion. En effet, le raffinement en proche paroi a pour effet d'étirer les cellules de gauches, ce qui pourrait induire des erreurs de calculs.

 

Comparaison

 

Enfin nous comparons les différentes simulations mises en place pour ce calcul le long de la plaque en traçant comme précédemment les profils de vitesse et de température aux niveaux des quatre positions de mesures expérimentales: x=1.438m, x=1.918m, x=2.535m et x=3.244m.

Remarque : La courbe en rouge correspond aux résultats d'un calcul en deux blocs avec pour manière d'établir l'écoulement laminaire en bas du domaine, la désactivation de la turbulence (turbulence OFF) dans le fichier transportpropeties.

La comparaison des différentes simulations à partir des profils de vitesse/température permet de valider le calcul en deux bloc. Cette méthode permet donc de d'obtenir les résultats les plus proches de la réalité et s'avère la plus logique puisqu'elle permet de prendre en compte la zone de transition, ce qui n'est pas le cas du calcul en "un bloc". Le seul inconvénient de cette méthode s'avère être le positionnement de la zone de transition!