Blog

Cette simulation s’inscrit dans la continuité de la série sur l’allée de tourbillons de Von Kármán. Cette fois-ci, au lieu d’une sphère, nous plaçons un prisme à l’intérieur du cylindre.

Comme OpenFOAM ne propose pas de forme de prisme prédéfinie (de type searchable), nous devons utiliser un fichier STL. Deux approches sont envisageables pour spécifier la géométrie :

1. Créer un prisme et l’intégrer dans un searchableCylinder.
2. Créer l’ensemble de la géométrie combinée dans un logiciel de CAO.

La seconde approche permet un bien meilleur contrôle des arêtes vives.

La géométrie a été créée sous FreeCAD puis exportée au format STL. Avant de configurer snappyHexMesh, il est nécessaire de préparer le maillage surfacique.

snappyHexMesh est capable de traiter les fichiers STL multi-solides. La première étape consiste donc à séparer le STL en différentes surfaces. Cela simplifie l’étape de placage du maillage (snapping) ainsi que la génération des couches limites prismatiques, car chaque surface du STL devient soit lisse (le cylindre), soit plane (toutes les autres faces).

Pour découper le fichier STL en surfaces distinctes, nous utilisons l’utilitaire surfacePatch. Cet utilitaire est une évolution de surfaceAutoPatch et permet de spécifier les surfaces ainsi que les algorithmes de découpe via le fichier surfacePatchDict :

geometry
{
    TubePrism.stl
    {
        type triSurfaceMesh;
    }
};

surfaces
{
    TubePrism.stl
    {
        type        autoPatch;
        featureAngle    30;
    }
}

L’exécution de cet utilitaire génère un fichier TubePrism_patched.stl contenant plusieurs solides :

solid Tube
 facet normal -0.999975 -1.57288e-18 0.00708072
  outer loop
   vertex -1.9998 12 0.0282381
   vertex -2 12 2.66454e-15
   vertex -2 -4 -8.88178e-16
  endloop
 endfacet
 ...
endsolid Tube
solid Outlet
...
endsolid Outlet
...
solid Prism-3
...
endsolid Prism-3

Remarque : Par défaut, surfacePatch nomme les solides selon un motif générique du type zone*. Les noms de patch Tube, Outlet et Prism-3 ont été attribués manuellement en éditant directement le fichier texte STL.

Après découpage, le fichier STL comprend 6 solides distincts : la paroi du tube, l’entrée (inlet), la sortie (outlet) et les trois faces du prisme.

Extraction des arêtes caractéristiques

La deuxième étape de préparation du maillage consiste à extraire les arêtes caractéristiques explicites à l’aide de l’utilitaire surfaceFeatureExtract.

Ces arêtes caractéristiques permettent de guider le raffinement du maillage de fond et d’améliorer le placage de la géométrie lors de la génération du maillage.

Génération du maillage

Une fois la préparation de la surface terminée, nous pouvons générer le maillage. La première phase consiste à créer le maillage crénelé (castellated mesh). Nous raffinons le maillage de fond le long des arêtes caractéristiques extraites :

features
(
    {
        file "TubePrism_patched.extendedFeatureEdgeMesh";
        level 2;
    }
);

Ainsi que sur l’ensemble de la surface géométrique :

refinementSurfaces
{
    TubeWithPrism
    {
        level (2 2);
    }
}

Le placage (snapping) s’appuie sur des paramètres standards ciblant à la fois les caractéristiques surfaciques explicites et implicites.

Ensuite, nous ajoutons trois couches limites surfaciques le long de la paroi du tube et des surfaces du prisme :

layers
{
    "TubeWithPrism_(Tube|Prism.*)"
    {
        nSurfaceLayers 3;
    }
}

L’insertion des couches limites repose sur le nom des patchs. Étant donné que nous avons précédemment renommé nos solides (Inlet, Outlet, Tube et Prism{1,2,3}), nous pouvons les cibler directement à l’aide d’expressions régulières dans la configuration des couches. Les autres paramètres d’ajout de couches conservent leurs valeurs par défaut.

L’exécution de snappyHexMesh permet d’obtenir le maillage final. Ci-dessous sont illustrées les couches prismatiques générées autour des patchs spécifiés :

Le maillage capture l’ensemble des détails géométriques requis, ce qui nous permet de passer à la configuration du solveur.

Configuration de la simulation

Comme pour les simulations précédentes, nous définissons la vitesse d’entrée ainsi que le nombre de Reynolds ciblé dans le dictionnaire personnalisé settings situé à la racine du cas :

// Flow conditions
Uin 1;
Re 400;

Nous calculons ensuite les propriétés dynamiques du fluide dans le fichier constant/transportProperties :

nu              #eval "$Uin*$D/$Re";

Script d’exécution HPC

Pour soumettre le calcul sur un cluster HPC géré par le gestionnaire de tâches Slurm, le script d’exécution en lot (batch) se présente comme suit :

#!/bin/bash

#SBATCH -p partition
#SBATCH -N 1
#SBATCH -n 32
#SBATCH -J VK/P
#SBATCH -t 8-00:00:00

module purge
module load openfoam/v2606

. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions

runApplication -o foamDictionary \
                       -set $SLURM_NTASKS \
                       -entry numberOfSubdomains system/decomposeParDict

[ -d processor0 ] || {
runApplication blockMesh
runApplication surfaceFeatureExtract
runApplication decomposePar
srun snappyHexMesh -overwrite -parallel > log.snappyHexMesh 2>&1
runApplication reconstructParMesh -constant
rm -rf processor*

restore0Dir
runApplication decomposePar
}

srun pimpleFoam -parallel > log 2>&1

Le script configure l’environnement à l’aide des modules, puis met à jour le paramètre numberOfSubdomains dans decomposeParDict pour correspondre au nombre de cœurs Slurm alloués (32 cœurs). Il génère ensuite le maillage de base avec blockMesh, extrait les arêtes caractéristiques, génère le maillage en parallèle avec snappyHexMesh, le reconstruit, et enfin exécute pimpleFoam sur l’ensemble des cœurs.

Post-traitement & Visualisation

Visualisation par lignes de courant

Une fois la simulation terminée, nous passons à la phase de visualisation. L’extraction du champ de vitesse dans le plan médian du tube confirme l’instabilité de l’écoulement :

Bien que la coupe dans le plan médian confirme le lâcher de tourbillons, il s’agit d’une vue 2D. Pour appréhender la structure tourbillonnaire globale en 3D, nous pouvons émettre des lignes de courant.

Le maillage au niveau de l’entrée étant relativement fin, utiliser l’intégralité du patch d’entrée comme source d’injection générerait un nombre excessif de lignes de courant. À la place, nous utilisons un Programmed Source dans ParaView afin d’injecter des points de départ personnalisés.

En supposant que le centre de l’entrée se situe en (0, 0), nous distribuons les points sources en coordonnées polaires (r*cos(theta), r*sin(theta)).

Remarque : Si le centre de votre plan d’entrée est décalé, ajoutez simplement le décalage en translation correspondant.

Dans le script du Programmed Source, nous définissons le nombre total de points, le rayon de distribution ainsi que le facteur multiplicateur de l’angle d’or pour assurer une répartition uniforme :

import vtk
from math import pi, sqrt, cos, sin

num_points = 200
radius = 0.5
golden_angle = pi * (3.0 - sqrt(5.0))

Nous créons ensuite les objets de sortie VTK :

output = self.GetOutput()
points = vtk.vtkPoints()
verts = vtk.vtkCellArray()

Les coordonnées des points sont calculées selon une distribution en spirale de Fibonacci :

for i in range(num_points):
    r = radius * sqrt(i / float(num_points))
    theta = i * golden_angle
    
    x = r * cos(theta)
    y = r * sin(theta)
    
    point_id = points.InsertNextPoint(x, -3.99, y)

    verts.InsertNextCell(1)
    verts.InsertCellPoint(point_id)

Enfin, nous assignons les points et sommets générés à l’objet de sortie VTK :

output.SetPoints(points)
output.SetVerts(verts)

Cela génère une distribution équilibrée de points d’injection sur le plan d’entrée :

L’application du filtre Stream Tracer à partir de ces points d’injection produit une animation qui révèle clairement la dynamique 3D des tourbillons en aval du prisme :

Iso-surfaces du critère Q

Une autre méthode très efficace pour identifier les structures tourbillonnaires est le critère Q, défini comme la moitié de la différence entre la norme au carré du tenseur de vorticité (Ω) et celle du tenseur des taux de déformation (S) :

Q = 0.5( ||Ω||² − ||S||² )

Lorsque Q est positif, le mouvement de rotation l’emporte sur la déformation par cisaillement, faisant des iso-surfaces de Q un excellent indicateur du cœur des tourbillons. Dans ParaView, le critère Q peut être calculé directement à l’aide du filtre Gradient :

Le tracé des iso-surfaces de Q positif donne l’animation d’écoulement suivante :

Bien que ces deux méthodes de rendu apportent des éclairages très intéressants, la représentation par lignes de courant offre une vision plus lisible des trajectoires de l’écoulement.

Comme d’habitude, l’ensemble des fichiers de configuration du cas ainsi que les géométries CAO sont disponibles dans le dépôt showcases d’O2M.

Cet article fait suite à la série de simulations sur l’allée de tourbillons de Von Kármán. Cette fois-ci, nous étendons le domaine en 3D et générons le maillage de calcul à l’aide de snappyHexMesh. La géométrie étant relativement simple, nous utilisons les surfaces recherchables (searchable surfaces) intégrées à OpenFOAM au lieu de fichiers STL externes.

Géométrie et maillage

Cette simulation est une extension en 3D de la classique allée de tourbillons de Von Kármán. Au lieu d’une géométrie plane en 2D, nous utilisons un tube cylindrique d’un diamètre de 4D et y plaçons une sphère de diamètre D. Nous appliquons une vitesse uniforme en entrée et observons le développement tridimensionnel de l’écoulement dans le sillage.

Le maillage est généré avec snappyHexMesh. Puisque notre domaine est constitué de formes géométriques de base, nous utilisons des objets géométriques intégrés (searchableCylinder et searchableSphere) plutôt que d’importer des fichiers STL externes :

geometry
{
    tube
    {
        type searchableCylinder;
        point1 (0 0 -100);
        point2 (0 0 100);
        radius #eval "2*$D";
    }

    sphere
    {
        type searchableSphere;
        origin (0 0 0);
        radius #eval "0.5*$D";
    }
}

Pour l’étape de maillage créneaux (castellated mesh), nous ne raffinons que les limites des surfaces. Les formes géométriques étant lisses et continues, l’extraction des arêtes n’est pas nécessaire (même l’entrée et la sortie sont définies par la boîte englobante du maillage de fond) :

    refinementSurfaces
    {
        sphere
        {
            level (1 1);
        }

        tube
        {
            level (1 1);
        }
    }

Pour résoudre les gradients de la couche limite le long des surfaces lisses, nous ajoutons des couches d’inflation prismatiques :

    layers
    {
        sphere
        {
            nSurfaceLayers 3;
        }

        tube
        {
            nSurfaceLayers 3;
        }
    }

Le reste du fichier snappyHexMeshDict utilise des paramètres standards.

L’exécution de snappyHexMesh produit un maillage propre avec une non-orthogonalité maximale d’environ 50° :

Note : Pour éviter que ParaView n’affiche de manière erronée les cellules limites coupées, l’option de découpe Crinkle clip est activée.

Checking geometry...
    Overall domain bounding box (-2 -2 -4) (2 2 12)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (-3.15769e-16 2.9073e-16 7.81195e-17) OK.
    Max cell openness = 4.07388e-16 OK.
    Max aspect ratio = 4.04672 OK.
    Minimum face area = 0.000805246. Maximum face area = 0.013633.  Face area magnitudes OK.
    Min volume = 4.0276e-05. Max volume = 0.00125246.  Total volume = 200.51.  Cell volumes OK.
    Mesh non-orthogonality Max: 50.0132 average: 7.10234
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.377543 OK.
    Coupled point location match (average 0) OK.

Conditions d’écoulement et propriétés du fluide

Comme pour le cas 2D, nous spécifions une vitesse constante en entrée (Uin) et faisons varier la viscosité cinématique (nu) pour atteindre le nombre de Reynolds ciblé (Re). Ceci est configuré de manière dynamique via l’entrée #eval dans constant/transportProperties :

#include "<case>/settings"

transportModel Newtonian;

nu              #eval "$Uin*$D/$Re";

Les variables Uin, D et Re sont définies dans le dictionnaire settings à l’intérieur du dossier du cas.

Visualisation de l’écoulement

Dans les simulations 2D, tracer les coefficients de force ou les cartes de contours de vitesse sur le plan médian fournit une image claire du champ d’écoulement. Cependant, en 3D, les coupes dans le plan médian ne parviennent souvent pas à capturer le détachement tourbillonnaire complexe autour de la sphère.

Une approche de visualisation bien supérieure consiste à tracer les lignes de courant en 3D.

Génération des lignes de courant

Approche sous ParaView

ParaView propose deux méthodes principales pour l’ensemencement (seeding) des lignes de courant :

1. Stream Tracer : Place les points de départ en utilisant des formes géométriques de base (ligne ou sphère).
2. Stream Tracer With Custom Source : Utilise une géométrie source personnalisée pour définir l’emplacement des points de départ.

Pour étudier les structures d’écoulement prenant naissance juste devant l’obstacle, il est utile d’ensemencer les lignes de courant à l’intérieur d’une zone circulaire sur le plan d’entrée correspondant à la projection de la sphère. Nous pouvons créer cette source d’ensemencement personnalisée dans ParaView grâce au pipeline suivant :

1. Extract Block: Isole le patch de la limite d’entrée.
2. Delaunay 2D: Triangule la géométrie du patch pour la préparer au raffinement.
3. Loop Subdivision: Raffine le maillage de surface triangulé.
4. Clip (Cylinder): Découpe le patch raffiné pour créer un disque d’ensemencement circulaire du rayon souhaité.

En réduisant le rayon de cette source personnalisée, nous pouvons visualiser clairement les motifs tourbillonnants complexes qui se forment dans le sillage de la sphère :

Enfin, l’application du filtre Stream Tracer permet de produire des animations fluides de l’écoulement :

Objet de fonction OpenFOAM streamLine

Les lignes de courant peuvent également être calculées directement pendant l’exécution ou en post-traitement à l’aide de l’objet de fonction streamLine d’OpenFOAM. Cet objet de fonction génère des particules lagrangiennes, suit leurs trajectoires à travers le champ de vitesse et enregistre les parcours résultants au format VTK.

Cependant, les options d’ensemencement intégrées sont assez limitées. L’option native la plus proche est patchSeed, qui distribue un nombre fixe de points sur l’ensemble d’un patch limite. Si on l’applique directement à notre patch d’entrée, patchSeed génère de nombreuses lignes de courant rectilignes sans grand intérêt près des parois extérieures du tube.

Pour restreindre l’ensemencement à une région spécifique, il nous faut une méthode du type faceZoneSeed. En combinant cela avec topoSet, nous pouvons définir une petite sous-zone sur l’entrée.

Implémentation de faceZoneSeedSet

L’examen de la classe patchSeedSet d’OpenFOAM montre qu’elle s’appuie sur triangulatedPatch pour distribuer les points d’ensemencement sur les faces d’un patch. Puisqu’une zone de faces (faceZone) peut être évaluée de manière similaire, nous pouvons facilement implémenter une classe personnalisée faceZoneSeedSet.

Nous adaptons patchSeedSet en remplaçant patchSet_ par faceZoneName_ et en supprimant selectedLocations_ :

class faceZoneSeedSet
:
    public sampledSet
{
    // Private Data
        //- Face zone name to sample
        const word faceZoneName_;
        
        //- Maximum number of zone faces to seed
        const label maxPoints_;
        
        //- Random number generator
        autoPtr<Random> rndGenPtr_;
        
    // Private Member Functions
        //- Calculate distribution points on the face zone
        void calcZoneSamples
        (
            const label nAvailable,
            const label nZonePoints,
            DynamicList<point>& samplingPts,
            DynamicList<label>& samplingFaces,
            DynamicList<label>& samplingSegments,
            DynamicList<scalar>& samplingCurveDist
        );

public:
    //- Runtime type information
    TypeName("faceZoneSeed");
    
    // Constructors
    faceZoneSeedSet
    (
        const word& name,
        const polyMesh& mesh,
        const meshSearch& searchEngine,
        const dictionary& dict
    );
};

Le positionnement des points d’ensemencement s’effectue dans calcZoneSamples, qui simule calcPatchSamples de la classe patchSeedSet :

void Foam::faceZoneSeedSet::calcZoneSamples
(
    const label nAvailable,
    const label nZonePoints,
    DynamicList<point>& samplingPts,
    DynamicList<label>& samplingFaces,
    DynamicList<label>& samplingSegments,
    DynamicList<scalar>& samplingCurveDist
)
{
    if (nAvailable < 1) return;

    Random& rndGen = *rndGenPtr_;

    globalIndex globalSampleNumbers(nAvailable);
    label nGlobalPatchPoints = returnReduce(nZonePoints, sumOp<label>());

    point pt;
    label facei;

    const faceZone& fz = mesh().faceZones()[faceZoneName_];
    triangulatedFaceZone tfz(fz);
    
    const label np = nAvailable * fz.size() / scalar(nGlobalPatchPoints);
    for (label i = 0; i < np; ++i)
    {
        tfz.randomLocalPoint(rndGen, pt, facei);
    
        samplingPts.append(pt);
        samplingFaces.append(fz[facei]);
        samplingSegments.append(0);
        samplingCurveDist.append(globalSampleNumbers.toGlobal(i));
    }
}

L’algorithme appelle randomLocalPoint via triangulatedFaceZone. Puisque les particules de lignes de courant résident dans le volume des cellules, nous associons les faces aux cellules intérieures à l’aide de faceOwner():

labelList samplingCells(samplingFaces.size(), -1);

forAll(samplingFaces, i)
{
    samplingCells[i] = mesh().faceOwner()[samplingFaces[i]];
}

Le code source complet des classes est disponible dans le dépôt d’exemples.

Configuration et Résultats

Nous configurons topoSetDict pour définir la faceZone cible (smallInlet) :

actions
(
  {
    name f0;
    type faceSet;
    action new;
    source patchToFace;
    sourceInfo
    {
      patch inlet;
    }
  }

  {
    name f0;
    type faceSet;
    action subset;
    source cylinderToFace;
    sourceInfo
    {
      point1 (0 0 -100);
      point2 (0 0 100);
      radius 0.2;
    }
  }

  {
    name smallInlet;
    type faceZoneSet;
    action new;
    source setToFaceZone;
    sourceInfo
    {
      faceSet f0;
    }
  }
);

Et la configuration suivante pour l’objet de fonction streamLine :

type        streamLine;
libs        (fieldFunctionObjects libfaceZoneSeed);
setFormat   vtk;

direction   forward;
fields      (p U);
lifeTime    10000;
nSubCycle   5;
cloud       particleTracks;

seedSampleSet
{
    type        faceZoneSeed;
    axis        xyz;
    zone        smallInlet;
    maxPoints   200;
}

L’exécution de l’objet de fonction génère les trajectoires de lignes de courant suivantes :

Note : Le maillage de fond au niveau de l’entrée étant relativement grossier, la limite de la zone d’ensemencement (faceZone) apparaît crantée/carrée. Dans ParaView, le fait de raffiner et de trianguler le patch extrait avant la découpe produit un profil d’ensemencement circulaire nettement plus lisse.

Les fichiers du cas de simulation sont disponibles dans le dépôt d’exemples.

Motivation

Il était une fois, des amis boulangers m’ont posé cette question : « Est-ce que la température d’un réfrigérateur vide met plus de temps à devenir uniforme que celle d’un réfrigérateur plein ? » Même si la question me semblait un peu contre-intuitive au premier abord (dans un réfrigérateur plein, il faut évacuer beaucoup plus d’énergie thermique emmagasinée), je me suis dit que ce serait une excellente simulation illustrative pour démontrer les capacités des solveurs de transfert de chaleur couplé (CHT) dans OpenFOAM.

Pour y répondre, nous allons réaliser deux simulations : un réfrigérateur vide et un réfrigérateur équipé d’une échelle de rangement.

Réfrigérateur vide

Le cas du réfrigérateur vide est le plus simple ; nous utilisons le solveur buoyantBoussinesqPimpleFoam. Le réfrigérateur étant un parallélépipède rectangle de dimensions 0,85×0,66×1,85 m, nous utilisons un maillage orthogonal avec une taille de maille de 1 cm. Une rapide revue de la littérature montre que les réfrigérateurs de boulangerie industrielle ont généralement un coefficient de transfert thermique surfacique constant aux parois de 6 W/(m²K). Nous appliquons donc les conditions aux limites thermiques externalWallHeatFluxTemperature partout. La température extérieure ambiante est fixée à 5 °C (278,15 K) avec une température interne initiale de 20 °C. Il serait évidemment possible d’être plus précis et d’attribuer des coefficients de transfert thermique différents selon les parois, mais une valeur uniforme convient parfaitement pour une première itération.

walls
{
    type            externalWallHeatFluxTemperature;
    mode            coefficient;
    Ta              constant 278.15;
    h               constant 6.0;
    kappaMethod     function;
    kappaValue      constant 25.87e-3;
    alphaValue      constant 2.135855960582023e-05;
    value           $internalField;
}

Pour mesurer la température interne, nous plaçons des sondes numériques à différentes hauteurs le long de l’axe vertical central.

probes
{
    type    probes;
    libs    (sampling);

    name    probes;

    writeControl   adjustableRunTime;
    writeInterval  0.1;

    fields (T);

    probeLocations
    (
        ( 42.5E-2 33E-2   5e-2)
        ( 42.5E-2 33E-2  15e-2)
        ...
        ( 42.5E-2 33E-2 155e-2)
        ( 42.5E-2 33E-2 165e-2)
        ( 42.5E-2 33E-2 175e-2)
    );
}

Dans le réfrigérateur vide, l’air froid à proximité des parois réfrigérées descend sous l’effet de la poussée d’Archimède. Cela crée une grande boucle de recirculation qui rend la température interne quasi uniforme (passant de 20 °C à 5 °C) en environ 10 minutes.

Cette durée est en très bon accord avec les observations réelles.

Réfrigérateur avec une échelle

Création de la géométrie de l’échelle

L’échelle est modélisée comme un assemblage de barres creuses en acier inoxydable : quatre montants verticaux, huit barres horizontales périphériques en haut et en bas, et 20 glissières latérales (10 de chaque côté) pour soutenir les plaques.

Compte tenu de la structure répétitive, OpenSCAD était l’outil idéal pour générer cette géométrie. Nous définissons d’abord les paramètres :

dx = 0.8;
dy = 0.624;
dz = 1.8;

dh = 18e-2;

thickness = 2e-3;
width = 2.54E-2;

Puis nous créons l’enveloppe externe de l’échelle et nous en soustrayons un volume intérieur légèrement plus petit pour représenter l’air coincé à l’intérieur des tubes.

    cube([dx, width, width]);
    translate([0, dy, 0])
        cube([dx, width, width]);
    translate([0, dy, dz - width])
        cube([dx, width, width]);
    translate([0, 0, dz - width])
        cube([dx, width, width]);
        
    cube([width, dy, width]);
    translate([dx, 0, 0]) cube([width, dy, width]);
    translate([dx, 0, dz - width]) cube([width, dy, width]);
    translate([0, 0, dz - width]) cube([width, dy, width]);
    
    for (i = [1:9]) {
        translate([0, 0, dh*i]) cube([dx, width, width]);
        translate([0, dy, dh*i]) cube([dx, width, width]);
    }
    
    for (i = [0:5]) {
        for (j = [1:9]) {
            translate([dw*i, 0, dh*j])
                cube([width, dy, width]);
        }
    }
}

union() {
    translate([thickness, thickness, thickness])
        cube([width_in, width_in, dz - 2*thickness]);

    translate([thickness, dy + thickness, thickness])
        cube([width_in, width_in, dz - 2*thickness]);

    translate([dx + thickness, dy + thickness, thickness])
        cube([width_in, width_in, dz - 2*thickness]);

    translate([dx + thickness, thickness, thickness])
        cube([width_in, width_in, dz - 2*thickness]);

    translate([thickness, thickness, thickness])
        cube([dx, width_in, width_in]);
    translate([thickness, dy + thickness, thickness])
        cube([dx, width_in, width_in]);
    translate([thickness, dy + thickness, dz - width_in - thickness])
        cube([dx, width_in, width_in]);
    translate([thickness, thickness, dz - width_in - thickness])
        cube([dx, width_in, width_in]);
       
    translate([thickness, thickness, thickness])
        cube([width_in, dy, width_in]);
    translate([dx + thickness, thickness, thickness])
        cube([width_in, dy, width_in]);
    translate([dx + thickness, thickness, dz - width_in - thickness])
        cube([width_in, dy, width_in]);
    translate([thickness, thickness, dz - width_in - thickness])
        cube([width_in, dy, width_in]);
    
    for (i = [1:9]) {
        translate([thickness, thickness, dh*i + thickness]) cube([dx, width_in, width_in]);
        translate([thickness, dy + thickness, dh*i + thickness]) cube([dx, width_in, width_in]);
    }
    
    for (i = [0:5]) {
        for (j = [1:9]) {
            translate([dw*i + thickness, thickness, dh*j + thickness])
                cube([width_in, dy, width_in]);
        }
    }
}
}

Après l’exécution du script et l’exportation en STL, nous obtenons cette belle structure :

L’échelle étant creuse, l’illustration ci-dessous montre la surface externe en beige, les cavités internes des glissières en gris et les arêtes du maillage STL en lignes blanches :

Implémentation du contrôle de convergence multi-région

Lorsque je développe des solveurs personnalisés, j’y intègre généralement des contrôles de convergence, qu’il s’agisse de configurations mono-région ou multi-régions. À ma grande surprise, il s’est avéré que la famille de solveurs CHT d’OpenFOAM ne dispose pas de contrôles de convergence intégrés. Les solveurs basés sur l’algorithme SIMPLE utilisent un nombre fixe de pas de temps spécifié via le paramètre endTime, tandis que les solveurs transitoires (famille PIMPLE) exécutent un nombre fixe d’itérations externes.

Une solution de contournement rapide consisterait à définir un nombre d’itérations externes suffisamment élevé pour garantir la convergence physique. Mais implémenter des contrôles de convergence natifs me semblait beaucoup plus amusant - et bien plus efficace, car cela permet de stopper les calculs dès que la convergence est atteinte.

OpenFOAM cherchant à réutiliser du code, il crée un ensemble de champs pour chaque région dans la famille des solveurs CHT. À chaque itération, il associe ces champs régionaux à des variables courantes (issues du fichier UEqn.H), puis réutilise ce fichier pour la résolution. Pour les solveurs basés sur l’algorithme SIMPLE, nous pouvons adopter la même approche : nous instancions des objets simpleControl dans chaque région, et à la fin de l’itération SIMPLE, nous bouclons sur ces objets pour vérifier si toutes les régions ont convergé. Si c’est le cas, nous interrompons les itérations et sauvegardons l’état final.

Pour les solveurs PIMPLE, la démarche est un peu plus subtile. On effectue des itérations externes au sein d’un même pas de temps. Lorsque les résidus passent sous un seuil configuré, la solution est considérée comme convergée et le solveur lance une itération finale supplémentaire avec un jeu distinct de solveurs linéaires. En général, on utilise une tolérance relative nulle et une tolérance absolue très faible lors de cette itération finale afin de garantir une excellente convergence de la solution à la fin du pas de temps.

Plutôt que d’évaluer directement les paramètres d’algorithmes dans chtMultiRegionFoam :

for (int oCorr=0; oCorr<nOuterCorr; ++oCorr)
{
	const bool finalIter = (oCorr == nOuterCorr-1);
	...
}

Nous reformulons la boucle pour ressembler aux solveurs mono-région comme pimpleFoam :

while (pimple.loop())
{
	...
	while (pimple.correct())
    {
		...
    }
}

L’idée derrière l’implémentation est simple : nous créons une classe de contrôle de solution par région (qui évalue la convergence locale) et une classe de contrôle multi-région qui parcourt la liste des contrôles régionaux, vérifie la convergence globale, puis stoppe la simulation ou passe au pas de temps suivant.

Contrôles de solution SIMPLE

Pour l’algorithme SIMPLE, l’implémentation est triviale : nous réutilisons la classe standard simpleControl, nous stockons les instances dans une liste et nous appelons la méthode criteriaSatisfied pour contrôler la convergence. Dès qu’elle est atteinte, nous arrêtons l’exécution.

Les modifications apportées au code de chtMultiRegionSimpleFoam sont minimes. Dans un premier temps, au lieu de :

while (runTime.loop())
{
	...
}

nous utilisons (msimple étant l’instance de notre classe de contrôle SIMPLE multi-région) :

while (msimple.loop())
{
	...
}

Dans un second temps, au lieu de créer manuellement la boucle de couplage énergétique, nous appelons la méthode energyCoupling du contrôle multi-région qui s’en charge en interne.

Contrôles de solution PIMPLE

Pour l’algorithme PIMPLE, l’implémentation demande un peu plus d’ajustements. Nous exploitons toujours la classe standard pimpleControl pour vérifier les résidus (via la méthode criteriaSatisfied), mais la méthode loop de la classe de contrôle PIMPLE multi-région doit non seulement vérifier la convergence, mais aussi déclencher l’itération finale lorsque les critères sont satisfaits dans toutes les régions.

Les modifications ultérieures dans le solveur restent similaires à celles de SIMPLE : nous utilisons l’instance du contrôle multi-région dans la boucle des correcteurs externes, puis nous nous en servons pour récupérer les instances régionales PIMPLE dans les boucles de correction internes et non-orthogonales.

Vous pouvez consulter l’ensemble de l’implémentation sur mon fork GitLab : https://gitlab.com/alexey.matveichev/openfoam

Comportement du refroidissement avec l’échelle

Maintenant que le contrôle de convergence est opérationnel, nous pouvons poursuivre les simulations de notre réfrigérateur (en ayant cette fois la certitude que nos calculs sont parfaitement convergés).

Pour générer le maillage, nous utilisons snappyHexMesh dans sa forme la plus simple : nous interrompons le processus après l’étape de maillage castellé pour éviter les problèmes de non-orthogonalité. Comme nous devons capturer la géométrie des parois en acier inoxydable de l’échelle, nous raffinons le maillage autour des surfaces et des arêtes du portique. Sur un maillage de fond de 1 cm, deux niveaux de raffinement suffisent.

features
(
	{
        file "echelle.extendedFeatureEdgeMesh";
        level 2;
    }
);

refinementSurfaces
{
    echelle
    {
        level (2 2);
    }
}

Pour découper le maillage en plusieurs zones, nous utilisons la liste locationsInMesh :

locationsInMesh
(
    (( 1e-4 1e-4  1e-4) air)
    (( 0.013  0.006  0.026) inox)
    (( 0.01600001  0.01000001  0.03000001) air-inside-echelle)
);

Nous obtenons trois régions distinctes : l’air dans la cavité du réfrigérateur, l’échelle métallique et l’air emprisonné à l’intérieur de l’échelle. Après la génération du maillage, snappyHexMesh crée les zones de mailles correspondantes, qui sont ensuite exploitées par l’utilitaire splitMeshRegions pour générer les maillages régionaux définitifs. Les coordonnées un peu « particulières » sont choisies pour éviter qu’un point ne tombe exactement sur une face.

Pour éviter les difficultés numeriques , nous traitons l’air piégé à l’intérieur des tubes de l’échelle comme un solide. Dans OpenFOAM, les résidus sont normalisés par rapport aux variations de la solution : si le champ est parfaitement convergé mais uniforme, les résidus restent proches de 1. L’air interne étant quasiment immobile, résoudre la quantité de mouvement dans ces cavités alourdirait inutilement les calculs.

Au début de la simulation, la température de l’échelle et de l’air du réfrigérateur est fixée à 20 °C. Le reste de la configuration est identique à la simulation du réfrigérateur vide.

Résultats de la simulation

La présence de l’échelle modifie complètement la circulation de l’air. Au lieu d’une unique grande boucle de recirculation, on observe une multitude de petites cellules fermées, coincées entre les différents niveaux de glissières. On conserve un écoulement descendant prononcé le long des parois froides du réfrigérateur, mais il est compensé par des courants ascendants le long des barres horizontales de l’échelle.

Cette restructuration de l’écoulement affecte directement le processus de refroidissement. Non seulement il faut extraire l’énergie thermique emmagasinée dans l’acier, mais l’évacuation de la chaleur devient beaucoup plus lente en raison du cloisonnement de l’air en cellules indépendantes. La partie inférieure refroidit un peu plus vite car elle bénéficie d’un contact direct avec le fond sans obstacle horizontal. Au sommet, les légères oscillations de température s’expliquent par deux phénomènes concurrents : le refroidissement par la paroi supérieure et l’apport thermique de l’air chaud ascendant généré par le portique.

L’introduction d’une échelle perturbe donc le refroidissement global et génère des zones de températures hétérogènes selon la hauteur. Dans un contexte professionnel, cela peut se traduire par des variations de pousse de la pâte d’un étage à l’autre.

Le pré-refroidissement est-il efficace ?

Pour vérifier si un pré-refroidissement du réfrigérateur permettrait d’atténuer le problème, nous avons simulé un maintien de l’enceinte vide à 5 °C pendant 10 minutes avant d’y insérer l’échelle à 20 °C.

Le résultat est sans appel : cela n’aide absolument pas. Les cellules de convection locales formées autour de l’échelle réchauffent très rapidement l’air ambiant. En moins d’une minute, le profil de température redevient strictement identique à celui d’un réfrigérateur non pré-refroidi.

Vous trouverez les fichiers de la simulation dans le dépôt dédié : https://code.o2m.solutions/o2m/show-cases/src/branch/trunk/fridge