Post-traitement des données

OpenFOAM

Cette rubrique a pour but d'illustrer les post-traitements que l'on a utilisé à l'aide des outils OpenFOAM. En effet, OpenFOAM dispose d'un certain nombre d'outils, majoritairement regroupé sous le nom de "Post-processing - Standard Utilities", qui permet d'analyser l'écoulement simulé.

 

foamCalc

foamCalc est un outil simple qui permet de calculer des grandeurs nécessaires au post-traitement. Pour la détermination de la vitesse en norme, on pourra utiliser la commande suivante.

rmartin@tanche : /work/rmartin/MaSimulation/system$    foamCalc mag U

Pour l'examination des composantes du tenseur des contraintes, on pourra utiliser la commande suivante.

rmartin@tanche : /work/rmartin/MaSimulation/system$    foamCalc sigmaxx

Pour connaître toutes les possibilités de l'utility foamCalc, on peut utiliser la commande suivante.

rmartin@tanche : /work/rmartin/MaSimulation/system$    foamCalc xxx

 

 

Sample

L'utility "sample" permet d'enregistrer des données de calcul le long de lignes pour tracer des graphes ou sur des surfaces pour l'affichage  d'isograndeurs en deux dimensions.
Pour bénéficier de cet outil, il suffit de spécifier les éléments nécessaires dans un fichier nommé "sampleDict" à placer dans le dossier "system".
Comme pour la plupart des fichiers OpenFOAM, il est intéressant d'aller chercher un exemplaire "sampleDict" dans les tutoriels proposés par OpenFOAM. Pour éviter de perdre du temps à chercher ce "dictionary", la commande suivante permet de recenser tous les fichiers du même nom, c'est une commande de recherche.

rmartin@tanche : /mnt/hmf/OpenFOAM/OpenFOAM-2.1.0$    find -name sampleDict

 

Une fois le fichier placé dans le répertoire /work/MaSimulation/system, ouvrir le fichier avec gedit.

rmartin@tanche : /work/rmartin/MaSimulation/system$    gedit sampleDict &

 

Le "SampeDict" ci-dessous permet de créer un répertoire /MaSimulation/sets contenant la grandeur spécifiée, ici Umag, sur les 3 lignes line_R=0.*m.

 

/*--------------------------------*- C++ -*---------------------------------------------------*\
| =====                        |                                               
| \\      /  F ield               | OpenFOAM: The Open Source CFD Toolbox          
|  \\    /   O peration       | Version:  2.1.0                                
|   \\  /    A nd                | Web:      www.OpenFOAM.org                    
|    \\/     M anipulation    |                                               
\*----------------------------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      sampleDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

interpolationScheme       cellPoint;              //choix d'interpolation des grandeurs demandées
                                                                // cellPoint : interpolation linéaire basée sur la taille des cellules

setFormat                     raw;                      //choix du format d'écriture 1D : raw - format ASCII en colonne                                                                           //exploitable sur matab. surfaceFormat pour une écriture 2D
sets                                                         //nom du répertoire dans lequel sera écrit les données
(

    line_R=0.2m                                          //nom d'écriture des fichiers
    {
        type         uniform;                             //répartition uniforme des 146 points sur la ligne
        axis          z;                                      //coordonnées de référence
        start         ( 0 0.2 -0.25 );                   //coordonnées du premier et du dernier point de la ligne
        end          ( 0 0.2 0.75 );
        nPoints     146;                                  //nombre de points sur la ligne
    }

    line_R=0.3m
    {
        type         uniform;
        axis          x;
        start         ( 0 0.3 -0.25 );
        end          ( 0 0.3 0.75 );
        nPoints     134;
    }

    line_R=0.4m
    {
        type         uniform;
        axis          x;
        start         ( 0 0.4 -0.25 );
        end          ( 0 0.4 0.75 );
        nPoints     135;
    }
);

fields          (magU);                               //liste des grandeurs à enregistrer sur les lignes

// ******************************************************************************* //

 

 

L'image ci-dessus illustre bien la position de ces trois lignes dans la cuve agitée.

 

Dans le cas de l'hélice en cuve réelle, nous tracerons les données selon quatre lignes horizontales à différentes hauteurs selon l'image suivante.

 

 

 

stressComponents

Comme nous avons expliqué dans la partie Caractérisation des Mobiles, nous pouvons calculer la puissance dissipée par la formule ci-dessous :

\begin{equation}
 P=\iiint\limits_V  \mu \phi_v \, dv
\end{equation}

Dans les coordonées cartésiennes, nous pouvons impliciter le terme de $ \phi_v $ comme suivant :

\begin{equation}
 \phi_v = [2 {\tau_{xx}} ^2+ 2 {\tau_{yy}} ^2 + 2 {\tau_{zz}} ^2 + {\tau_{xz}} ^2 + {\tau_{yz}} ^2 + {\tau_{xy}} ^2] / {\mu ^2}
\end{equation}

Avec :

\begin{equation}
{\tau_{xx}}= - 2 {\mu} \frac { \partial U} {\partial x}
\end{equation}

\begin{equation}
{\tau_{yy}}= - 2 {\mu} \frac { \partial V} {\partial y}
\end{equation}

\begin{equation}
{\tau_{zz}}= - 2 {\mu} \frac { \partial W} {\partial z}
\end{equation}

\begin{equation}
{\tau_{xy}}= - {\mu} (\frac { \partial U} {\partial y}+\frac { \partial V} {\partial x})
\end{equation}

\begin{equation}
{\tau_{xz}}= - {\mu} (\frac { \partial W} {\partial x}+\frac { \partial U} {\partial z})
\end{equation}

\begin{equation}
{\tau_{yz}}= - {\mu} (\frac { \partial V} {\partial z}+\frac { \partial W} {\partial y})
\end{equation}

Pour pouvoir calculer le puissance, il faut demander à OpenFOAM à sortir les valeurs des contraintes. Il existe une fonction auxiliaire sur l'OpenFOAM : stressComponents qui permet de sortir toutes les valeurs de $ \sigma_{ij} $.

L'équation qui relie $ \sigma_{ij} $ et $ \tau_{ij} $ est comme suit :

\begin{equation}
\sigma_{ij} = -p I_{ij}+\lambda (div U_i) I_{ij} +2 \mu D_{ij}
\end{equation}

avec :

$ I_{ij} $ est la matrice d'indentité

$ div (U_i) = 0 $ car c'est incompressible

$ D_{ij} = \frac {1} {2} (grad U_i + grad^t U_i) $

Pour sortir les composantes du tenseur sigma, lancez la commande suivante :

rmartin@tanche : /work/rmartin/MaSimulation$    stressComponents

 

Afin de mieux comprendre comment OpenFoam calcule $ \sigma_{ij} $ à partir des gradients de vitesse, nous allons sur la librairie d'OpenFOAM :

/mnt/hmf/OpenFOAM/OpenFOAM-2.1.0/applications/utilities/postProcessing/stressField/stressComponents/

 

Essayons de comprendre le script de stressComponents.C :

volSymmTensorField sigma
            (
                IOobject
                (
                    "sigma",
                    runTime.timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
                ),
                laminarTransport.nu()*2*dev(symm(fvc::grad(U)))   // sigma est le produit de $ \nu $ et des gradients de U. OpenFoam
            );                                                                               

L'unité de la puissance est le Watt et sa dimension est $ [M{L^2} {T^{-3}}] $. La dimension de la viscosité dynamique est $ [M L^{-1} T^{-1}] $ et la dimension de volume est $ [L^3] $. A l'aide de l'analyse dimensionnelle, on peut déduire que la dimension de $ \phi_v $ est $ [T^{-2}] $. 

Nous pouvons donc déduire qu'il suffit de normaliser le terme de $ \sigma_{ij} $ obtenu avec OpenFOAM par $ \nu $ pour avoir le terme de $ \tau_{ij} $.

Sur paraview, afin d'obtenir le terme de $ \phi_v $, il faut normaliser le terme de $ \sigma_{ij} $ calculé par OpenFOAM comme suit :

$$ \phi_v= [2 \sigma_{xx}^2 +2 \sigma_{yy}^2 + 2 \sigma_{zz}^2 + \sigma_{xy}^2 +\sigma_{xz}^2 + \sigma_{yz}^2 ]/ {\nu^2}  $$

 

 

pyFoamPlotWatcher

pyFoamPlotWatcher.py est un script python très utile puisqu'il peut afficher les résidus et la continuité d'une simulation en cours avec gnuplot.

Pour une visualisation des résidus, il faut rentrer les deux commandes suivantes.

rmartin@tanche : /work/rmartin/MaSimulation/system$ MRFSimpleFoam > log &
rmartin@tanche : /work/rmartin/MaSimulation/system$ pyFoamPlotWatcher.py log &

 

Voici ce que l'on peut obtenir :

Matlab

L'ensemble des données extraites par les utilities OpenFOAM ont été traité sur Matlab. Nous allons détailler en quelques lignes le script Matlab utilisé.

Nous avons précisé une extension en .xy pour toutes les données que le sampleDict sortirait grâce à l'option setformat : raw. L'exploitation sous Matlbab se fait sans soucis.

L'exploitation des données se fait en deux étapes :

  • l'importation sous Matlab des données
  • le tracé des graphes correctement renseignés

 

Importation des données

L'importation des données se fait avec l'outil "importdata

F02_SIMPLE = importdata('500_Umag_0.2m_Fluent_SIMPLE.xy');
OF02 = importdata('line_R=0.2m_magU.xy')

 

 

Tracé des graphes

Afin de disposer de graphes lisibles sur la majorité des supports, optez pour les commandes set, fontsize et linewidth pour grossir et passer un gras les caractères.

figure(4);
hold on;
grid on;
plot(F02_SIMPLE(:,1),F02_SIMPLE(:,2), 'g','linewidth',3);
plot(F02_SIMPLE(:,1),OF02(:,2), 'b','linewidth',3);
xlabel('Height Z  m ', 'fontsize',22);
ylabel('Velocity Mag  m/s', 'fontsize',22);
set(gca, 'fontsize',22);
grid on;
title('Velocity Mag for R=0,2m at it=500', 'fontsize',22);
legend('Fluent 13.0 - SIMPLE', 'OpenFoam 2.1.0 - SIMPLE');
set(legend, 'fontsize',22);

 

Le graphe associé :

 

Paraview

Une partie du post-traitement des données calculées par OpenFOAM peut se faire à l'aide de paraview, notamment pour la détermination des nombres adimensionnels caractéristiques.

En effet, d'après la définition des nombres adimensionnels de débit, l'utilisation de paraview peut se révéler utile pour la création de surfaces à travers desquelles peuvent être déterminés des débits.

Cette rubrique a pour but d'illustrer la méthode employée pour la détermination grandeurs globales.

 

Méthode à suivre :

Commencer par ouvrir un terminal et aller dans le répertoire ou se trouve l'ensemble des fichiers de la simulation OpenFOAM :

rmartin@tanche :   cd /work/rmartin/MaSimulation
rmartin@tanche : /work/rmartin/MaSimulation$   ls
0   10   20   30   40   constant   system   log

Pour pouvoir lancer les commandes OpenFOAM depuis n'importe quel répertoire :

rmartin@tanche : /work/rmartin/MaSimulation$    source /mnt/hmf/OpenFoam/OF21.sh

Pour bénéficier de tous les outils de Paraview, on préférera lancer Paraview avec les deux commandes suivantes plutôt que paraFoam directement :

rmartin@tanche : /work/rmartin/MaSimulation$   touch .OpenFoam
rmartin@tanche : /work/rmartin/MaSimulation$   /usr/local/bin/paraview

Pour pouvoir exploiter les données OpenFOAM avec paraview, il faut convertir les données dans un format que peut lire Paraview. Le format .vtk est adapté. Pour convertir les données en .vtk il suffit de lancer la commande suivante.

rmartin@tanche : /work/rmartin/MaSimulation$   foamToVTK

Ouverture de paraview avec la commande suivante.

rmartin@tanche : /work/rmartin/MaSimulation$   paraview

L'ouverture des fichiers .vtk se fait manuellement dans un premier temps : => Files => Open => *.vtk

Pour la suite des commandes, on divisera les cas en fonction du type de pompage de l'agitateur : axial ou radial.

 

 

Calcul du débit de pompage

 

Cas des mobiles d'agitation à pompage axial :

Afin de déterminer le débit de pompage d'un agitateur à pompage axial, il faut créer une surface spécifique qui servira de support à l'intégration de la vitesse axiale. La détermination de ce débit de pompage se divise en 3 parties.

Voici le cas d'une cuve à agitation axiale (une hélice qui tourne selon l'axe-y) qui servira d'illustration :

 

Création d'un cylindre avec le Programmable Filter (attention, cette option n'est pas disponible dans paraFoam)

L'outil Programmable Filter est un outil puissant puisqu'il permet de créer un filtre coder en python correspondant à la fonction recherchée. Nous allons créer un Clip Cylinder, filtre qui n'existe pas dans la version utilisée, à savoir la 12.6.0 :

Pour bénéficier de l'outil Programmable Filter : => Filters => Alphabetical => Programmable Filter

Une fenêtre permet de rentrer l'ensemble des lignes de code nécessaire. Pour ce qui concerne la création d'un cylindre par un Clip, nous avons rentré l'ensemble des lignes de code suivante :

input = self.GetInputDataObject(0, 0)
inp_copy = input.NewInstance()
inp_copy.ShallowCopy(input)
inp_copy.UnRegister(None)

transf = vtk.vtkTransform()
transf.RotateY(90)

cyl = vtk.vtkCylinder()
cyl.SetCenter(0,0,0)
cyl.SetRadius(1)
cyl.SetTransform(transf)

clipper = vtk.vtkClipDataSet()
clipper.SetClipFunction(cyl)
clipper.SetInput(inp_copy)
clipper.InsideOutOn()
clipper.Update()

self.GetOutputDataObject(0).ShallowCopy(clipper.GetOutputDataObject(0))

Ce programme permet de créer un cylindre  de dimension infini selon Y grâce à une fonction implicite.

Le paramètrage des dimensions du cylindre dans les deux autres directions se fait via les deux commandes suivantes :

cyl.SetCenter(0,0,0) permet de fixer le centre du cylindre
cyl.SetRadius(1) permet de choisir le rayon du cylindre

On obtient la visualisation suivante :

C'est bien un cylindre de rayon de 0.95 m ayant un diamètre plus petit que celui de la cuve.

 

Création du slice disque :

Le Clip_Cylindre précédent a permis de délimiter la zone dans les plans normaux à Y. La création de la surface en forme de disque se fait à l'aide du filtre Slice :

=> Filters => Alphabetical => Slice

On choisit de créer un Slice_plane qui, délimité par le Clip_cylindre, revient à la création d'un disque.

Il faut bien choisir la bonne normale qui est selon l'axe Y pour créer un disque.

On obtient la visualisation suivante dans le cas d'un slice fait trop prêt de l'hélice :

 

Calcul de l'intégrale de la vitesse verticale sur le disque avec Integrate Variables

On peut utiliser la fonction IntegrateVariables sur paraview pour calculer le débit volumique de fluide qui passe à travers le disque défini dans la partie précédent.

=> Filters => Alphabetical => Integrate Variables

La fonction Integrate Variables est capable de calculer la grandeur souhaitée sur une surface ou un volume prédéfini. Attention à bien appliqué le filtre IntegrateVariables sur le Slice : ceci se fait en cochant le Slice défini dans la colonne de droite avant d'utiliser le filtre.

Le calcul par IntegrateVariables montre que le débit de fluide qui passe par le disque selon l'axe Y est d'environ  0,21 $ m^3 /s $.

 

 

Cas des mobiles d'agitation à pompage radial :

 

Afin de déterminer le débit de pompage d'un agitateur à pompage radial qui tourne selon l'axe-Z, il faut créer une surface d'un cylindre de la même hauteur de l'agitateur qui servira de support à l'intégration de la vitesse radiale.

La géométrie du cas test 3D sera utilisée pour l'illustration de la procédure à suivre dans le cas des agitateurs à pompage radial.

 

Calcul de vitesse radiale avec calculator

Le logiciel de l'OpenFOAM travaille dans un système de coordonnées cartésiennes. Afin de calculer le débit de pompage à travers la surface d'un cylindre, il nous faut tout d'abord la vitesse radiale à l'aide de l'outil calculator sur paraview.

=> Filters => Alphabetical => Calculator

On a des équations comme ci_dessous qui nous permettent de passer du système de coordonées cartesiennes au système de coordonées cylindriques :

\begin{equation}
U_r=U_x cos (\theta)+U_y sin(\theta)
\end{equation}

\begin{equation}
U_\theta=- U_x sin (\theta) +U_y cos(\theta)
\end{equation}

avec :

$$ \theta= arctg (\frac {coordY} {coordX} ) $$

Attention à la fonction non définie pour coord-x $ =0 $ car la valeur de $ \theta $ tend vers l'infini. Pour éviter ce problème, on peut définir les équations de $ U_r $ et $ U_{\theta} $ autrement :

\begin{equation}
U_r=U_x \frac {coordX} {\sqrt {{coordX}^2+{coordY}^2}} +U_y \frac {coordY} {\sqrt {{coordX}^2+{coordY}^2}}
\end{equation}

\begin{equation}
U_\theta=- U_x \frac {coordY} {\sqrt {{coordX}^2+{coordY}^2}} +U_y \frac {coordX} {\sqrt {{coordX}^2+{coordY}^2}}
\end{equation}

On peut modifier la nouvelle variable qu'on définit par la fonction calculatrice comme U_r (la vitesse radiale).

 

Création d'un cylindre avec le Programmable Filter (attention, cette option n'est pas disponible dans paraFoam)

input = self.GetInputDataObject(0,0)

inp_copy = input.NewInstance()

inp_copy.ShallowCopy(input)

inp_copy.UnRegister(None)

cutter = vtk.vtkCutter()

transf = vtk.vtkTransform()

transf.RotateX(90)

cyl = vtk.vtkCylinder()

cyl.SetCenter(0,0,0)

cyl.SetRadius(0.26)

cyl.SetTransform(transf)

cutter.SetCutFunction(cyl)

cutter.SetInput(inp_copy)

cutter.Update()

self.GetOutputDataObject(0).ShallowCopy(cutter.GetOutputDataObject(0))

Ce script de python ne marche que sur Paraview 6.14.

On obtient la visualistion suivante :

 

Création d'une boîte d'une hauteur 0.5m contenant l'hélice

On se rappelle qu'on veut faire l'intégral de la vitesse radiale sur la surface transversale d'un cylindre ayant la même hauteur que les pâles d'agitateur. Il faut donc enlever la partie haute et la partie basse de la cuve où il n'y pas d'agitateur à l'aide de la fonction clip sur paraview:

=> Filters => Alphabetical => Clip

On choisit le clip type : box et il faut absolument cocher Inside Out pour garder la partie contenant l'hélice.

On obtient donc la visualisation suivante :

 

Calcul de l'intégrale de la vitesse radiale sur la surface transversale du cylindre avec Integrate Variables

On peut utiliser la fonction Integrate Variables sur paraview pour calculer le débit volumique de fluide qui passe la surface transversale du cylindre.

=> Filters => Alphabetical => Integrate Variables

 

 

Calcul de la puissance pour des mobiles à pompage axial ou radial

 

Comme nous avons expliqué dans la partie OpenFOAM, il faut intégrer le champ de sigma (les contraintes de cisaillement) ou $ \phi_v $ avec la viscosité dynamique sur tout le volume de la cuve.

 

 

Calcul de $ \mu \phi_v $ à l'aide de calculator

Afin de pouvoir utiliser la fonction Integrate Variables pour calculer l'intégral sur tout le volume, il faut calculer abord le produit du terme de $ \phi_v $ normalisé avec la viscosité dynamique à l'aide de la calculator. Il faut rentrer la formule ci-dessous sur la fonction de calculator :

$$ \mu [2 \sigma_{xx}^2 +2 \sigma_{yy}^2 + 2 \sigma_{zz}^2 + \sigma_{xy}^2 +\sigma_{xz}^2 + \sigma_{yz}^2 ]/ {\nu^2}  $$

Avec : $ \mu=0,04 SI $ et $ \nu=2,67.10^{-5} SI $

L'image illustre la formule de $ \mu \phi_v $ entrée sur la fonction calculator sur Paraview.

 

Calcul de l'intégrale de $ \mu \ phi_v$ sur tout le volume avec Integrate Variables

Une fois le produit de $ \mu \phi_v $ calculé, nous utilisons la fonction Integrate Variables pour calculer l'intégrale de $ \mu \phi_v $ sur tout le volume de la cuve.

=> Filters => Alphabetical => Integrate Variables