Cavité entrainée/Lid-driven Cavity Flow

Dans cette partie, nous nous intéressons au cas d'étude de la cavité entraînée: un écoulement de fluide incompressible dans un domaine carré 2D, où le mur supérieur se déplace vers la droite à la vitesse U. On cherche alors à résoudre les champs de vitesse et pression à l'aide d'OpenFoam.

Géométrie de la cavité

Figure 1: Géométrie de la cavité
Source:http://www.openfoam.org/docs/user/cavity.php#x5-40002.1

On génère initialement un maillage de 20*20 cellules sous OpenFoam grâce au générateur de maillage blockMesh (commande à entrer dans un terminal), dont les spécifications se trouvent dans le fichier constant/ployMesh/blockMeshDict. Tout les maillages sont automatiquement générés en trois dimensions, mais il suffit de spécifier une condition au limite de type "empty" pour les murs perpendiculaires à la troisième dimensions (murs de devant et de derrière, perpendiculaire à l'axe Oz). Les autres murs de la cavité seront du type "fixedWalls" ou "movingWalls".

Les conditions initiales des champs de vitesse et de pression sont rentrés dans les fichiers 0/P et 0/U. Les propriétés physiques du fluide (ici, la viscosité $ \nu $ en m²/s) sont spécifiées dans le fichier transportProperties. On peut ainsi définir le nombre de Reynolds de l'écoulement: $ Re =\frac{d \cdot |U|} {\nu} $

  • Si l'écoulement est laminaire, on utilise le solveur icoFoam (laminaire, isotherme, incompressible)
  • Si l'écoulement est turbulent, on utilise le solveur pisoFoam (incompressible, modèle RAS)

On définit dans system/controlDict les constantes de contrôle du calcul, comme la pas de temps, le pas d'écriture, le temps de début et de fin du calcul... On fera en sorte que le calcul s'arrête alors que l'écoulement est établi (on règle endTime), et que le nombre de Courant $ Co =\frac{\delta t \cdot |U|} {\delta x} $$ \delta t $ est le pas de temps, et $ \delta x $ est la taille d'une cellule dans la direction de la vitesse est inférieur à 1 (on règle deltaT).

 

  1. Cavité carrée

On commence par lancer une simulation en écoulement laminaire sur la géométrie décrite précédemment (on se place dans le répertoire : /tutorials/incompressible/icoFoam/cavity). Il est d'ailleurs conseillé de copier ce répertoire à un autre emplacement - dans le home par exemple, sauf si les calculs prennent trop de place - nous pourrons ainsi faire toute les modifications qui vont suivre sans risquer "d'endommager" le tutorial.

On choisit $ endTime = 0.5 s $  , $ |U| = 1 m/s $ et une pression initiale uniforme dans tout le domaine. De plus $ \nu = 0.1 m²/s $ , l'écoulement est donc bien laminaire ( $ Re =\frac{0.1 \times 1} {0.1}=1 $ ).  On lance alors les calculs (commande icoFoam & dans un terminal, le & servant à mettre le calcul en arrière plan), et il vient:

Figure 2: Champs de Pression à $ T_{final} $ dans la cavité carrée

On constate bien une dépression dans le coin supérieur gauche, et une surpression dans le coin supérieur droit, la plaque supérieur avançant vers la droite.

Figure 3: Champs de Vitesse à $ T_{final} $ dans la cavité carrée

Nous sommes ici en régime établi (0.05 s suffise). Il est intéressant d'observer les vecteurs vitesse (filtre Cell Centers - qui permet de prendre la valeur au centre de la cellule - puis Glyph dans ParaView) dans la cavité:

Figure 4: Vecteurs Vitesse à $ T_{final} $ dans la cavité carrée

L'écoulement ainsi visualisé correspondant bien à nos attentes théoriques. La plaque supérieure mise en mouvement vers la droite met en rotation le fluide visqueux (petit Reynolds) par un effet cisaillement (c'est le même phénomène que pour un écoulement de couette), pour créer un tourbillon dans la partie supérieure de la cavité. Le centre de celui-ci est à vitesse nulle. On aurait pu attendre une décentralisation du tourbillon vers la droite, pour une cavité plus large par exemple.

Nous pouvons cependant noter que cette représentation laisse penser qu'il existe une vitesse au niveau des parois fixe. Ce n'est cependant pas le cas si on s'y intéresse de près.

 

NB: Si l'écoulement avait été turbulent (vitesse plus grande, où viscosité plus petite), comme c'est la cas dans la partie suivante, il aurait fallu utiliser la tutorial présent dans /tutorials/incompressible/pisoFoam/ras/cavity, et toujours en prenant soin de copier le dossier ailleurs, pour ne pas endommager le tutorial préparé.

 

  1. La cavité tronquée

Nous souhaitons à présent modifier notre cavité, en raffinant son maillage et en tronquant sa géométrie. Pour ce faire, présentons d'abord le fichier à l'origine de la génération de celui-ci: blockMeshDict (dans constant/polyMesh).

Figure 5: Fichier BlockMeshDict pour la cavité carrée en maillage $ 20 \times 20 $

Il est nécessaire ici de souligner certains points qui vont nous être utile par la suite:

  • Ligne 17: La fonction "convertToMeters" indique qu'il faut multiplier toutes les distances par un certain nombre, ici 0.1, pour l'avoir en mètre.
  • Ligne 19 à 29: On définit ici les "vertices", ie les différents points de la géométrie. Nous avons bien 8 points, correspondant au 8 sommets du cube formant notre cavité. En effet, la génération du maillage est automatiquement en 3D sous OpenFoam, même si notre étude se limite à 2 dimensions. La ligne 17 nous indique bien que 1 = 0.1 m = d. La numérotation des points commence à 0.
  • Ligne 31 à 34: On définit maintenant les "blocks", à partir des numéros des points précédemment définis. La géométrie que nous avons étudiée dans la partie précédente était faite d'un seul bloc: le cube définissant notre cavité.

C'est également à dans ce groupement que nous définissons la résolution de notre maillage (le nombre de cellule le composant), grâce à la parenthèse entourée sur l'image: (20 20 1) correspondant au nombre de cellule selon x, y et z. Nous avons ici $ 20 \times 20 $ cellules. Nous n'avons qu'une seule cellule selon z, notre étude étant 2D. Il suffit alors, pour raffiner le maillage, de remplacer (20 20 1) par (50 50 1) par exemple, nous obtenons un maillage de $ 50 \times 50 $ , et les résultats seront plus précis.

A présent, nous allons complexifier la géométrie, en la définissant grâce à quatre blocs au lieu d'un seul, comme le montre la figure ci-dessous.

Figure 6: Définition de la nouvelle géométrie en 4 blocs, avec numérotation de la cavité tronquée
Source:http://www.openfoam.org/docs/user/cavity.php#x5-40002.1

Pour cela, il faut tout d'abord définir 10 nouveaux points dans "vertices", ceux sur les arêtes et sur les faces du cube, on obtient alors 18 points, numérotés de 0 à 17, correspondant aux sommets de nos quatre blocs. Puis, nous définissons les quatres nouveaux blocs dans "blocks".

Cette manière de définir une géométrie un bloc permet dans un premier temps de définir une finesse de maillage propre à chaque bloc, mais surtout, elle va nous permettre de tronquer notre cavité au niveau du coin inférieur droit (par exemple), en enlevant le bloc 1. Pour ce faire, nous enlevons les points A et B du groupement "vertices" (respectivement (1 0 0.1) et (1 0 0)). Nous ne définissions donc plus que 16 points. Il ne nous reste qu'à définir les trois nouveaux blocs (attentions, la numérotation des points change ! ), comme suit:

      hex (0 1 3 2 8 9 11 10) (24 16 1) simpleGrading (1 1 1) bloc 0
      hex (2 3 6 5 10 11 14 13) (24 24 1) simpleGrading (1 1 1)
bloc 2
      hex (3 4 7 6 11 12 15 14) (16 24 1) simpleGrading (1 1 1) bloc 3

A noter que les points nouvellement définis sur les segments et faces du cube ne le sont pas aux milieux de ceux-ci, ce qui explique la répartition des cellules du maillage à 40/60% - (24 16 1) au lieu de (20 20 1):

Numéros des points Coordonnées
Numérotation des nouveaux points en fonction de leurs coordonnées
0 (0 0 0)
1 (0.6 0 0)
2 (0 0.4 0)
3 (0.6 0.4 0)
4 (1 0.4 0)
5 (0 1 0)
6 (0.6 1 0)
7 (1 1 0)
8 (0 0 0.1)
9 (0.6 0 0.1)
10 (0 0.4 0.1)
11 (0.6 0.4 0.1)
12 (1 0.4 0.1)
13 (0 1 0.1)
14 (0.6 1 0.1)
15 (1 1 0.1)

Nous obtenons alors le maillage suivant (celui de droite, celui de gauche étant le maillage de la cavité étudiée dans la première partie):

Figure 7: Comparaison entre le maillage de la cavité carrée, et celui de la cavité tronquée, plus fin

 


 

Une fois le maillage ainsi obtenu, nous définissons les conditions initiales et aux limites de notre simulation.

En conséquence, il convient de définir les types de murs de notre nouvelle géométrie, et nous le faisons en modifiant la suite du fichier blockMeshDict, (plus précisemment, le groupement "boudary"):

Figure 8: Fin du fichier BlockMeshDict correspondant à la cavité tronquée

En plus de simplement tronquer la cavité, nous avons décidé d'imposer une vitesse sur le mur vertical de la troncature.

  • Ligne 53 à 70: Nous définissons donc deux "movingWalls". Un mur est nécessairement défini par 4 points au minimum. Le mur supérieur, aux vues de la disposition des blocs, est défini par les deux surfaces (5 13 14 6) et (6 14 15 7). Ce type de définition ("movingWalls") va permettre d'imposer une vitesse initiale non nulle dans 0/U sur ces deux murs, comme nous le verrons par la suite.
  • Ligne 72 à 82: Ce sont les murs inférieurs et ceux de côté, qui sont définis comme des parois fixes, "fixedWalls".
  • Ligne 84 à 94: Ce sont les murs de devant et de derrière, de type "empty" pour avoir un problème 2D.

Nous avons pris le parti de séparer en deux groupements les deux "movingWalls". Cette décision va nous permettre de définir des conditions initiales différentes sur ces deux murs, ce qui ouvre beaucoup plus de perspective en pratique, surtout concernant la vitesse:

Figure 9: Fichier O/U correspondant à la cavité tronquée

On applique ici une vitesse initiale $ U_1 $ de 5 m/s vers la droite (axe Ox) sur le mur supérieur, et une vitesse initiale $ U_2 $ de 5 m/s vers le bas (axe -Oy) sur le mur de la troncature. On visualise la vitesse de l'écoulement de la cavité tronquée à t=0 s:

Figure 10: Contour et Vecteur Vitesse de la cavité tronquée (zoom sur mur de la troncature) à t = 0 s

De plus, on applique un gradient de pression uniforme sur la totalité de la cavité.

NB : En faisant cette séparation des "movingWalls", nous avons changé la structure du fichier blockMeshDict. Il faut donc faire attention à conserver cette même structure dans les autres fichiers.

 


 

Nous allons maintenant lancer nos calculs. On souhaite se placer en écoulement turbulent, et donc utiliser le solveur pisoFoam (cf introduction de ce paragraphe).  Nous avons modifier le tutorial déjà crée dans /tutorials/incompressible/pisoFoam/ras/cavity, et dans lequel $ \nu = 10^{-5} m²/s $ .

On a donc $ Re =\frac{0.1 \times 5} {\nu}=5000 $, le régime est bien turbulent. On se limite à une durée de simulation de 2s, suffisante pour que l'écoulement s'établisse.

Il faut prendre soin de régler deltaT de sorte que le nombre de Courant $ Co \le 1 $. On a vu que:

$$ Co =\frac{\delta t \cdot |U|} {\delta x} $$

$ \delta t $ est le pas de temps, et $ \delta x $ est la taille d'une cellule dans la direction de la vitesse (ie $ \delta x = \frac{d}{n} $ , n étant le nombre de cellule le long de d, ici 40 (24+16)).

On doit donc avoir $$ \delta t \le \frac{d}{n \cdot |U|} $$ 

$$ \Rightarrow \delta t \le \frac{0.1}{40 \cdot 5} $$ 

$$ \Rightarrow \delta t \le 0.0005 $$ 

On prendra cette valeur de $ \delta t $ .

On obtient la simulation suffisante:

Animation 1: La cavitée tronquée, contour de vitesse pour $ U_1 $ = 5 m/s selon Ox, $ U_2 $ = -5 m/s selon Oy

Mettons ce résultat en comparaison avec la cavité tronquée à $ U_2 $ = 0 m/s (on retrouve ici un écoulement similaire à celui d'une conduite avec une marche, soit l'écoulement sur une discontinuité):

Animation 2: La cavitée tronquée, contour de vitesse pour $ U_1 $ = 5 m/s selon Ox, $ U_2 $ = 0 m/s

Certains phénomènes sont alors flagrants. Le tourbillon supérieur, celui dans la partie pleine de la cavité, est du au cisaillement entrainé par la plaque supérieure. C'est le même tourbillon que celui observé dans la cavité carré, mais "perturbé" par la présence de la marche, que nous retrouvons dans les deux simulations.

Cependant, la vitesse de la marche introduit une perturbation supplémentaire antagoniste au mouvement de la plaque supérieure. On va donc observer, pendant la phase transitoire, des vitesses légèrement plus faible (elles s'annulent) accompagnées d'un léger mouvement oscillatoire.

Regardons à présent l'évolutions des lignes de courant pour ces deux simulations:

Animation 3: La cavitée tronquée, ligne de courant pour $ U_1 $ = 5 m/s selon Ox, $ U_2 $ = -5 m/s selon Oy

Le plus gros tourbillon est ici le phénomène principal, mais on constate la formation d'un tourbillon "parasite" sur la parti supérieure de la marche.

Animation 4: La cavitée tronquée, ligne de courant pour $ U_1 $ = 5 m/s selon Ox, $ U_2 $ = 0 m/s

En comparant ces deux animations, on voit que l'écoulement s'établi plus rapidement dans le cas de $ U_2 $ = 0 m/s. De plus, contrairement au cas précédent, il y a une recirculation dans la totalité de la cavité inférieur. Cette recirculation est empéchée par le mouvement de la marche, ce qui explique la formation du petit tourbillon dont nous avons parlé après l'animation 3.

 

  1. Ajout de la température

Il s'agit à présent de modifier le solveur, pour rajouter le paramètre température. On veut alors résoudre les équations de Navier-Stokes ci-dessous :

$ \nabla U =0 $

$ \frac{\partial U}{\partial t} + \nabla(U.U) = \nabla (\nu \nabla U) - \nabla p $

$ \frac{\partial T}{\partial t} + \nabla(U.T) = \nabla (Dt \nabla U) $

On veut donc rajouter l'équation de la chaleur dans le nouveau solveur. On va alors créer un répertoire application en exécutant dans le terminal : mkdir -p applications/solvers/incompressible

On copie ensuite dans ce répertoire le solveur icoFoam déjà existant en le renommant pour le modifier ensuite. 

On modifie dans un premier temps, le fichier createField. En effet, on y déclare le scalaire DT en rajoutant :

 

dimensionedScalar DT
(
     transportProperties.lookup("DT")
);
 
De plus, on déclare aussi le scalaire T, en ajoutant :
 
Info<< "Reading field T\n" <<endl;
volScalarField T
(
    IOobject
    (
         "T",
         runTime.timeName(),
         mesh,
         IOobject::MUST_READ,
         IOobject::AUTO_WRITE
     ),
     mesh
);

 

 

On rajoute alors dans le fichier myIcoFoam.C l'équation de transport de température qui s'écrit sous cette forme : 

 

fvScalarMatrix TEqn
        (
            fvm::ddt(T)
            + fvm::div(phi, T)
            - fvm::laplacian(DT, T)
        );
 
        TEqn.solve();

 

Il faut alors compiler le nouveau solveur en exécutant wmake.

 

Le solveur maintenant créé, il faut modifier le répertoire cavity pour adapter les conditions aux limites et les conditions initiales. Dans un premier temps, on fait une copie du répertoire IcoFoam/cavity.

On crée un fichier T dans le répertoire 0. Pour cela, on peut copier le fichier p et faire les modifications nécessaires pour obtenir le fichier T ci-dessous :

                                             figure 11: fichier T/0 correspondant à la cavité après ajout de la température.

Il faut aussi ajouter dans le fichier transportProperties la ligne : 

DT            DT [0 2 -1 0 0 0 0] 0.002;

 

On doit ensuite modifier les fichiers fvSolution et fvSchemes. Dans le fichier fvSolution, on modifie de sorte à obtenir le même code que ci-dessous: 

 

solvers
{
    p PCG
    {
        //information à propos du solveur de pression
    };
//à ajouter
    T 
    {
       solver            BICCG;
        preconditioner   DILU;
        tolerance        1e-7;
        relTol           0;
    };
 
et dans le fichier fvSchemes, il faut ajouter :
 
divSchemes
{
    default         none;
    div(phi,U)      Gauss linear; 
    div(phi,T)      Gauss upwind; 
}
 
laplacianSchemes
{
    default         none;
    laplacian(nu,U) Gauss linear corrected;
    laplacian((1|A(U)),p) Gauss linear corrected;
    laplacian(DT,T) Gauss linear corrected;
}
 
 

On peut désormais utiliser le solveur myIcoFoam sur la cavité entraînée. 

On se place à une vitesse $ U_x =1 m/s $ sur le mur supérieur. En ce qui concerne la température, au temps initial, on a 350K au niveau du mur supérieur et 300K au niveau des autres murs fixes. Après exécution de myIcoFoam, on obtient le résultat suivant :

 

Animation 5 : champs de température dans la cavité

 

On observe bien les deux phénomènes de convection et de diffusion ce qui est cohérent avec la théorie. En effet, on peut voir que le flux de température est modifié par la vitesse au niveau du mur supérieur : la température de la cavité est plus élevée en haut à droite. Par ailleurs, on observe bien le phénomène de diffusion car la température tend à s'homogénéiser dans la cavité.

 


 

Nous avons pu, à partir d'un tutorial déjà mené, étudier un problème tout à fait différent. La connaissance de la structure et du fonctionnement, même génral d'OpenFoam, nous ouvre déjà un large champs de possibilités. Uniquement grâce à la connaissance des techniques précédentes, nous aurions pu établir une tout autre géométrie, tronquer un autre coin, faire un escalier, imposer d'autres conditions aux limites, chauffer, etc. afin de résoudre un problème complètement différent de l'original.

Malgré un aspect d'abord peu attrayant d'openFoam (pas d'interface graphique en particulier), nous avons constaté qu'en se penchant sur le sujet, nous arivons à avoir rapidement une connaissance suffisante pour pouvoir travailler efficacement.

 

 

 

Suite: Rayleigh-Benard

Précédent: Barrage/Breaking of a Dam