Blog : Simuler la dynamique de refroidissement d'un réfrigérateur avec OpenFOAM
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/(m2K). 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