Motivation
Un jour, j’ai simulé l’allée de tourbillons de Karman : le maillage avait été
généré avec Gmsh, j’avais joué avec différents schémas de discrétisation et
contrôles de résidus, puis tracé les résultats de post-traitement VTK avec
Matplotlib. Cette fois, j’ai décidé de créer un cas test d’évaluation
(benchmark) qui exploite toutes les nouvelles fonctionnalités pratiques
introduites dans OpenFOAM depuis lors.
Géométrie & Maillage
La géométrie est un canal rectangulaire 2D simple de 16D × 8D, où D
représente le diamètre de notre cylindre (voir figure ci-dessous). Le cylindre
est situé à 4D en aval de l’injecteur (inlet). Cela correspond à la géométrie
de l’article (https://onlinelibrary.wiley.com/doi/abs/10.1002/fld.1683), où l’on
trouve également les valeurs des coefficients de force. Bien que le but de ces
simulations soit de montrer les capacités des dictionnaires OpenFOAM, disposer
de valeurs de référence est un plus appréciable.
Mailler un cylindre dans un canal
Pour mailler une géométrie avec l’utilitaire blockMesh, nous devons d’abord la
découper en blocs hexaédriques. Nous plaçons ensuite les sommets des blocs dans
la liste vertices, les blocs dans la liste blocks, la description des arêtes
courbes dans la liste edges, et enfin les patchs de frontière du maillage dans
la liste boundary.
Entrées de fonction dans les dictionnaires OpenFOAM
Les dictionnaires dans OpenFOAM utilisent une syntaxe de type C ; par
conséquent, il existe des entrées de fonction qui commencent par # pour imiter
le préprocesseur C. La première entrée de fonction que nous utiliserons est
include, qui inclut simplement le contenu d’un autre dictionnaire à l’endroit
de son appel. Pour simplifier la gestion du cas, nous déplaçons toutes les
variables modifiables par l’utilisateur vers un dictionnaire settings situé
dans le dossier racine du cas, qui est ensuite inclus via l’entrée
#include "<case>/settings".
En plus du stockage des variables, nous devons effectuer certains calculs
(densité du maillage, nombre de mailles dans les zones de progression, etc.).
Auparavant, OpenFOAM fournissait une entrée calc à cet effet. Celle-ci créait
un fichier C++ avec l’expression, le compilait, et insérait le résultat de
l’exécution du code dans le dictionnaire. Dans les versions modernes, l’entrée
eval évalue directement les expressions sans compiler de bibliothèque
dynamique externe.
Enfin, il existe des entrées conditionnelles : if et ifeq. La première
vérifie si l’expression est vraie (au sens d’OpenFOAM, l’expression peut donc
être yes ou on). La seconde compare deux expressions. Dans les deux cas, si
l’évaluation est vraie, la partie après if est insérée dans le dictionnaire ;
sinon, la partie optionnelle else est utilisée.
Dictionnaire de configuration
Dans le dictionnaire settings, nous conservons le diamètre du cylindre, la
taille des mailles et les conditions d’écoulement. Le premier définit la
géométrie du canal, la deuxième définit la densité du maillage, et les dernières
définissent les propriétés du fluide. Comme le détachement tourbillonnaire
dépend du nombre de Reynolds, nous calculons la viscosité du fluide à partir du
diamètre du cylindre, de la vitesse d’entrée et du nombre de Reynolds.
Nous utilisons un maillage 2D dans la simulation : nous avons donc une taille
verticale de maillage (dz) et deux tailles moyennes de mailles : dans la zone
d’écoulement libre (dx) et près du cylindre (dx_fine). Pour réduire le
nombre de mailles, nous pouvons utiliser le rafinement progressif (grading).
Dictionnaire blockMeshDict
Avant de rédiger le dictionnaire, nous devons découper notre géométrie en blocs
hexaédriques et définir les sommets d’ossature qui serviront à construire le
maillage.
Le schéma des blocs du maillage est présenté dans la figure ci-dessus. Afin
d’obtenir une taille de maille normale uniforme au niveau du patch du cylindre,
nous introduisons des sommets supplémentaires autour de celui-ci. Notre maillage
compte ainsi 20 blocs : 4 près de l’entrée, 4 près de la sortie et 12 dans la
zone du cylindre.
Le fichier blockMeshDict commence par l’inclusion du dictionnaire de
paramètres, puis nous calculons des variables pratiques en utilisant l’entrée de
fonction eval.
#include "<case>/settings"
R #eval "0.5*$D";
a #eval "$R/sqrt(2)";
b #eval "$D/sqrt(2)";
c #eval "4*$D";
d #eval "12*$D";
z0 -$dz;
z1 $dz;
Ensuite, la densité du maillage et les facteurs de progression sont calculés.
Ici, nous utilisons l’entrée conditionnelle #if et une solution analytique
pour le nombre de mailles dans la zone progressive. Comme nous avons des
directions de progression différentes selon les zones, nous calculons deux
facteurs de progression : un pour le raffinement et un autre pour le
déraffinement.
#if $grading
q #eval "$dx_fine/$dx";
p #eval "$dx/$dx_fine";
Ng #eval "ceil(log($q)/log(($b -$dx)/($b -$dx_fine)))";
#else
q 1;
p 1;
Ng #eval "ceil($b/$dx)";
#endif
Puis nous définissons le nombre de mailles dans les différentes parties du
maillage ; comme nous avons 4 subdivisions dans les directions X et Y, 4 valeurs
sont définies.
Nx1 #eval "ceil(($c - $b)/$dx)";
Nx2 $Ng;
Nx3 $Ng;
Nx4 #eval "ceil(($d - $b)/$dx)";
Ny1 #eval "ceil(($c - $b)/$dx)";
Ny2 $Ng;
Ny3 $Ng;
Ny4 #eval "ceil(($d - $b)/$dx)";
Nous avons besoin d’arêtes courbes dans notre maillage : les arêtes du
cylindre et les arêtes des blocs autour de celui-ci. Dans les anciennes versions
d’OpenFOAM, la définition des arêtes nécessitait deux points d’extrémité et un
troisième point milieu le long de l’arc (définissant l’arc par trois points).
Désormais, nous pouvons projeter une arête droite sur la géométrie courbe. Pour
cela, nous devons d’abord définir la géométrie dans le dictionnaire geometry.
Deux cylindres y sont définis : le premier avec la vraie géométrie du cylindre
et le second, deux fois plus grand, pour les blocs d’« uniformisation ».
geometry
{
cylinder1
{
type cylinder;
point1 (0 0 -100);
point2 (0 0 100);
radius $R;
}
cylinder2
{
type cylinder;
point1 (0 0 -100);
point2 (0 0 100);
radius $D;
}
}
Pour la définition des sommets, nous utilisons des sommets nommés. Dans les
anciennes versions d’OpenFOAM, il était nécessaire d’utiliser l’index du sommet
dans la définition des blocs (en commençant par 0). Cette approche d’indexation
manuelle est fastidieuse et source d’erreurs. Ainsi, pour les sommets de cette
simulation, nous utilisons des noms. Et pour simplifier la description des blocs
de maillage, nous utilisons des mnémoniques pour le nommage :
v<décalage en X><décalage en Y><décalage en Z>
Les numéros de décalage correspondent aux valeurs sur le schéma ci-dessus.
L’utilisation de sommets nommés nous évite de devoir suivre l’index des sommets,
et la convention de décalage simplifie la construction des blocs.
vertices
(
name v000 (-$c -$c$z0)
name v100 (-$b -$c$z0)
name v300 ( 0 -$c$z0)
...
);
La description des blocs utilise les sommets nommés et les variables
précalculées (densités de blocs, facteurs de progression).
blocks
(
hex (v000 v100 v110 v010 v001 v101 v111 v011) ($Nx1$Ny1 1) simpleGrading ( 1 1 1)
hex (v100 v300 v310 v110 v101 v301 v311 v111) ($Nx2 $Ny1 1) simpleGrading ($q 1 1)
hex (v300 v500 v510 v310 v301 v501 v511 v311) ($Nx3 $Ny1 1) simpleGrading ($p 1 1)
hex (v500 v600 v610 v510 v501 v601 v611 v511) ($Nx4$Ny1 1) simpleGrading ( 1 1 1)
...
);
Les arêtes courbes du maillage sont décrites via la projection des arêtes
droites sur la géométrie définie.
edges
(
project v110 v310 (cylinder2)
project v310 v510 (cylinder2)
project v510 v530 (cylinder2)
project v530 v550 (cylinder2)
project v550 v350 (cylinder2)
project v350 v150 (cylinder2)
project v130 v110 (cylinder2)
project v130 v150 (cylinder2)
...
);
Après avoir défini les blocs et les arêtes, nous pouvons passer aux patchs de
frontière. blockMesh possède un patch spécial nommé defaultFaces, où toutes les
faces non attribuées sont ajoutées. Nous pouvons configurer le comportement de
ce patch via le dictionnaire defaultFaces.
defaultPatch
{
name cylinder;
type wall;
}
Le reste des patchs comprend inlet, outlet, walls, ainsi que le patch
frontAndBack de type empty, classique pour une simulation 2D.
inlet
{
type patch;
faces
(
(v000 v010 v011 v001)
(v010 v030 v031 v011)
(v030 v050 v051 v031)
(v050 v060 v061 v051)
);
}
Pour générer le maillage, il suffit donc de définir le diamètre du cylindre
ainsi que la taille de maille souhaitée, puis de lancer blockMesh (voir la
figure ci-dessous pour un exemple de maillage généré).
Le maillage résultant offre une résolution fine dans la zone d’intérêt et
conserve une non-orthogonalité raisonnable (autour de 40 degrés).
Checking geometry...
Overall domain bounding box (-4 -4 -0.01) (12 4 0.01)
Mesh has 2 geometric (non-empty/wedge) directions (1 1 0)
Mesh has 2 solution (non-empty) directions (1 1 0)
All edges aligned with or perpendicular to non-empty directions.
Boundary openness (-6.7961e-19 3.04482e-18 -1.66941e-13) OK.
Max cell openness = 4.34617e-16 OK.
Max aspect ratio = 9.60439 OK.
Minimum face area = 2.52457e-06. Maximum face area = 0.000632222. Face area magnitudes OK.
Min volume = 5.04914e-08. Max volume = 1.26444e-05. Total volume = 2.54429. Cell volumes OK.
Mesh non-orthogonality Max: 44.1274 average: 7.02609
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 0.433471 OK.
Coupled point location match (average 0) OK.
Schémas, solveurs linéaires, convergence et contrôle du pas de temps
Les schémas de discrétisation sont configurés dans le fichier fvSchemes situé
dans le sous-dossier system du cas. En raison de la non-orthogonalité du
maillage, nous utilisons le schéma de reconstruction de gradient leastSquares.
Le schéma upwind de second ordre est employé pour le terme de convection, et les
schémas linear corrected sont utilisés pour le laplacien et l’interpolation.
gradSchemes
{
default leastSquares;
}
divSchemes
{
default none;
div(phi,U) Gauss linearUpwind grad(U);
div((nuEff*dev2(T(grad(U))))) Gauss linear;
}
Pour le solveur pimpleFoam, nous avons besoin de deux jeux de solveurs
linéaires : un pour les itérations externes avant la convergence, et un
autre pour l’itération finale, effectuée une fois la convergence atteinte.
Pour la pression, nous utilisons le solveur multigrille algébrique géométrique
(GAMG). La convergence du solveur linéaire est contrôlée par deux tolérances :
relative et absolue. La première mesure la réduction du résidu ; dans notre
cas, si le résidu initial de l’équation de pression est réduit d’un facteur 10,
nous considérons que le solveur a convergé. Pour l’itération finale, nous
basculons sur le solveur de gradient conjugué préconditionné (PCG) avec GAMG
comme préconditionneur. Cette fois, nous réglons la tolérance relative à 0 afin
d’utiliser la tolérance absolue comme critère de convergence du solveur
linéaire. GAMG utilise le réagglomération du maillage pour résoudre un système
linéaire plus petit avant de propager la solution sur le maillage dense. Notre
maillage étant statique, nous pouvons mettre en cache l’agglomération grâce au
drapeau cacheAgglomeration.
p
{
solver GAMG;
smoother GaussSeidel;
tolerance 1E-8;
relTol 0.1;
cacheAgglomeration true;
}
pFinal
{
solver PCG;
preconditioner
{
preconditioner GAMG;
smoother GaussSeidel;
tolerance 1E-8;
relTol 1E-2;
cacheAgglomeration true;
}
tolerance 1E-8;
relTol 0;
}
L’utilisation de la tolérance relative lors de la convergence initiale permet
d’accélérer la résolution. Réalisant un couplage pression-vitesse, nous
n’utilisons pas de tolérances trop strictes au départ pour les champs de
pression ou de vitesse, ces deux champs n’étant alors qu’une solution
intermédiaire. Une fois arrivés à la convergence globale, nous appliquons la
tolérance absolue pour garantir une solution finale bien convergée.
Pour contrôler la convergence de la solution, on utilise le sous-dictionnaire
residualControl du dictionnaire PIMPLE. On y indique le nom des champs, la
tolérance relative et la tolérance absolue.
residualControl
{
"(U|p)"
{
relTol 0;
tolerance 1e-6;
}
}
À chaque itération externe, le solveur vérifie le résidu initial du système
linéaire pour chaque champ configuré (dans notre exemple, nous avons utilisé une
expression régulière pour cibler les champs de vitesse et de pression, mais ils
peuvent être définis dans des sous-dictionnaires distincts). Si le résidu
initial passe sous la tolérance absolue ou si le rapport entre résidu actuel et
résidu initial passe sous la tolérance relative, la solution est considérée
comme convergée.
La stabilité de la simulation dépend de la valeur du pas de temps. Une façon de
la garantir consiste à fixer un pas de temps très court, mais cela rend le
calcul extrêmement lent. Il est donc bien plus pratique d’adapter le pas de
temps au champ d’écoulement. À cette fin, pimpleFoam propose les paramètres
suivants dans le fichier controlDict :
adjustTimeStep yes;
maxCo 0.5;
Une valeur assez conservatrice de 0,5 est choisie pour le nombre de Courant
maximal en raison de l’utilisation de schémas de discrétisation de second
ordre ; un schéma upwind classique permettrait d’utiliser une valeur
plus élevée.
Post-traitement en cours de calcul (Run-time)
Nous devons non seulement calculer le champ d’écoulement, mais aussi extraire
certaines données de la simulation. Cela se fait au fil du calcul grâce aux
objets fonctions (function objects). Ils sont configurés dans le
sous-dictionnaire functions du fichier controlDict. Pour en faciliter la
gestion, chaque description est isolée dans un fichier séparé, tous inclus
ensuite dans le sous-dictionnaire correspondant.
functions
{
#include "<case>/system/forceCoeffs"
#include "<case>/system/surfaces"
#include "<case>/system/timeInfo"
}
Coefficients de forces
forces
{
type forceCoeffs;
libs (forces);
writeControl adjustableRunTime;
writeInterval 0.1;
patches (cylinder);
rho rhoInf;
log true;
rhoInf 1;
liftDir (0 1 0);
dragDir (1 0 0);
CofR (0 0 0);
pitchAxis (0 0 1);
magUInf $Uin;
lRef $D;
Aref #eval "2*$dz*$D";
}
Afin de comparer les résultats de simulation avec les données expérimentales,
nous devons extraire les forces de traînée et de portance s’exerçant sur le
cylindre. Nous utilisons pour cela l’objet fonction forceCoeffs. La
configuration spécifique débute à partir du mot-clé patches. Nous y indiquons
les patchs soumis aux forces (ici, la frontière cylinder). Nous fixons la masse
volumique du fluide via rho. Notre simulation étant incompressibles, nous
renseignons une valeur constante rhoInf. Pour le calcul des forces, nous devons
définir les directions de traînée et de portance, le centre de rotation, l’axe
de tangage, la vitesse amont, la longueur de référence et la surface de
référence. Dans notre cas, la traînée est portée par l’axe X, la portance par
l’axe Y, le centre de rotation est le centre du cylindre, l’axe de tangage est
Z, et la vitesse amont est la vitesse d’entrée. La longueur de référence
correspond au diamètre du cylindre et la surface de référence est cette longueur
multipliée par l’épaisseur du maillage selon Z.
Cet objet fonction enregistre dans le fichier
postProcessing/forces/0/coefficients.dat les valeurs des coefficients :
traînée (Cd), portance (Cl), force latérale (Cs), ainsi que les moments de
roulis (CmRoll), tangage (CmPitch) et lacet (CmYaw). Il fournit également
ces coefficients décomposés selon leurs composantes avant et arrière (suffixes
(f) et (r)).
Échantillonnage (Sampling)
velocity
{
type surfaces;
libs (sampling);
writeControl writeTime;
surfaceFormat raw;
fields
(
U
);
interpolationScheme cellPoint;
surfaces
{
centerPlane
{
type cuttingPlane;
point (0 0 0);
normal (0 0 1);
interpolate true;
triangulate true;
}
}
}
Nous extrayons également les valeurs de vitesse le long du plan central du canal
en utilisant l’objet fonction cuttingPlane. Il permet d’exporter les données
échantillonnées sous différents formats. Pour le tracé des courbes, nous avons
retenu le format VTK, le cylindre créant un “trou” dans les données qu’il
convient de trianguler correctement.
timeInfo
{
type timeInfo;
libs (utilityFunctionObjects);
perTimeStep yes;
writeControl timeStep;
writeInterval 1;
}
Pour recueillir quelques métriques sur les performances du solveur, nous faisons
appel à l’objet fonction timeInfo. À chaque pas de temps, il affiche le temps
simulé, le temps CPU total, le temps réel (wall clock), ainsi que les temps CPU
et réel consommés par pas de temps. La résolution du temps réel étant de 1
seconde, il s’avère un peu moins exploitable.
Résultats
Pour donner une dimension de benchmark à nos simulations, nous avons choisi de
comparer les résultats de différentes versions d’OpenFOAM. Les versions v2306,
v2312, v2412 et v2606 étant installées sur le cluster, les calculs ont été
exécutés sur l’ensemble de ces versions avec trois finesses de maillage : 5k,
30k et 300k mailles.
L’influence de la densité du maillage sur le résultat est comparée à partir de
la version v2606. La vitesse d’exécution du solveur est quant à elle comparée
entre les différentes versions d’OpenFOAM.
Coefficients de traînée et de portance
Bonne nouvelle : toutes les versions d’OpenFOAM testées calculent exactement la
même valeur pour les coefficients de traînée et de portance :
L’écart constaté sur les coefficients Cd et Cl par rapport à l’article cité
provient principalement des conditions aux limites appliquées sur les parois du
canal. Si l’on remplace les parois par des plans de symétrie (ce qui revient à
simuler un réseau infini de cylindres, comme dans l’article), les valeurs de
Cd et Cl se rapprochent nettement des données publiées.
Pour évaluer les écarts de performance entre versions, le solveur a été exécuté
sur un intervalle allant de t=150 s à 160 s avec un pas de temps fixe (0,002 s)
sur le maillage de 30k mailles.
Le temps CPU de la version v2606 a été pris comme référence basique, et les
différences de temps CPU ont été tracées. Il s’avère que la version v2606 est la
plus lente. La plus rapide est la v2306, avec un gain de 14 secondes par rapport
à la v2606. L’ensemble de la simulation ayant nécessité 1157 secondes CPU, la
v2306 est donc 1 % plus rapide : un écart qui tient à la marge de tolérance des
mesures et peut être négligé.
Fichiers du cas
Les fichiers complets du cas sont disponibles sur le dépôt d’exemples :
https://code.o2m.solutions/o2m/show-cases/src/branch/trunk/von-karman-vortex-street
