N°d’ordre NNT : 2024ISAL0034

THÈSE de DOCTORAT DE L’INSA LYON,

membre de l’UNIVERSITÉ DE LYON
délivré en partenariat international avec
UNIVERSITÉ DE SHERBROOKE

École Doctorale N° 162

MÉCANIQUE, ÉNERGETIQUE, GÉNIE CIVIL, ACOUSTIQUE

Spécialité / discipline de doctorat :

Thermique-Énergétique

Soutenue publiquement, le 24/04/2024 par :

Jihyuk Jeong

Modélisation CFD du transfert de

chaleur et de masse dans une remorque
de camion réfrigérée équipée de
plaques eutectiques

CFD Modeling of the Heat and Mass

Transfer in a Refrigerated Truck Trailer
Equipped with Eutectic Plates

Devant le jury composé de :

FOURNAISON, Laurence, Directrice de Recherche, INRAE, Présidente
SAFDARI SHADLOO, Mostafa, Maître de Conférences HDR, INSA Rouen, Rapporteur
GOSSELIN, Louis, Professeur, Université Laval (Québec), Rapporteur
FERTEL, Camille, Docteure, Tecnea Canada, Examinatrice
MICHEL, Benoît, Maître de Conférences, INSA Lyon, Examinateur
BONJOUR, Jocelyn, Professeur des Universités, INSA Lyon, Directeur de thèse
PONCET, Sébastien, Professeur, Université de Sherbrooke, Co-directeur de thèse

Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Département de la Formation par la Recherche
et des Études Doctorales (FEDORA)
Bâtiment INSA direction, 1er étage
37, av. J. Capelle
69621 Villeurbanne Cédex
fedora@[Link]

Référence : TH1093_JEONG Jihyuk

L’INSA Lyon a mis en place une procédure de contrôle systématique via un outil de
détection de similitudes (logiciel Compilatio). Après le dépôt du manuscrit de thèse,
celui-ci est analysé par l’outil. Pour tout taux de similarité supérieur à 10%, le manuscrit
est vérifié par l’équipe de FEDORA. Il s’agit notamment d’exclure les auto-citations, à
condition qu’elles soient correctement référencées avec citation expresse dans le
manuscrit.

Par ce document, il est attesté que ce manuscrit, dans la forme communiquée par la
personne doctorante à l’INSA Lyon, satisfait aux exigences de l’Etablissement concernant
le taux maximal de similitude admissible.

Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Département FEDORA – INSA Lyon - Ecoles Doctorales

SIGLE ECOLE DOCTORALE NOM ET COORDONNEES DU RESPONSABLE

CHIMIE DE LYON M. Stéphane DANIELE

C2P2-CPE LYON-UMR 5265
ED 206 [Link] Bâtiment F308, BP 2077
Sec. : Renée EL MELHEM 43 Boulevard du 11 novembre 1918
CHIMIE Bât. Blaise PASCAL, 3e étage 69616 Villeurbanne
secretariat@[Link] directeur@[Link]

ÉVOLUTION, ÉCOSYSTÈME, MICROBIOLOGIE, MODÉLISATION Mme Sandrine CHARLES

Université Claude Bernard Lyon 1
[Link]
ED 341 UFR Biosciences
Sec. : Bénédicte LANZA Bâtiment Mendel
E2M2 Bât. Atrium, UCB Lyon 1 43, boulevard du 11 Novembre 1918
Tél : [Link].62 69622 Villeurbanne CEDEX
secretariat.e2m2@[Link] [Link]@[Link]

INTERDISCIPLINAIRE SCIENCES-SANTÉ Mme Sylvie RICARD-BLUM

Laboratoire ICBMS - UMR 5246 CNRS - Université Lyon 1
[Link]
ED 205 Bâtiment Raulin - 2ème étage Nord
Sec. : Bénédicte LANZA 43 Boulevard du 11 novembre 1918
EDISS Bât. Atrium, UCB Lyon 1 69622 Villeurbanne Cedex
Tél : [Link].62 Tél : +33(0)4 72 44 82 32
[Link]@[Link] [Link]-blum@[Link]

MATÉRIAUX DE LYON M. Stéphane BENAYOUN

Ecole Centrale de Lyon
[Link]
ED 34 Laboratoire LTDS
Sec. : Yann DE ORDENANA 36 avenue Guy de Collongue
EDML Tél : [Link].44 69134 Ecully CEDEX
[Link]-ordenana@[Link] Tél : [Link].37
[Link]@[Link]

ÉLECTRONIQUE, ÉLECTROTECHNIQUE, AUTOMATIQUE M. Philippe DELACHARTRE

[Link] INSA LYON
ED 160 Sec. : Philomène TRECOURT Laboratoire CREATIS
Bâtiment Direction INSA Lyon Bâtiment Blaise Pascal, 7 avenue Jean Capelle
EEA Tél : [Link].70 69621 Villeurbanne CEDEX
Tél : [Link].63
[Link]@[Link] [Link]@[Link]

INFORMATIQUE ET MATHÉMATIQUES M. Hamamache KHEDDOUCI

Université Claude Bernard Lyon 1
ED 512 [Link] Bât. Nautibus
Sec. : Renée EL MELHEM 43, Boulevard du 11 novembre 1918
INFOMATHS Bât. Blaise PASCAL, 3e étage 69 622 Villeurbanne Cedex France
Tél : [Link].46 Tél : [Link].69
infomaths@[Link] [Link]@[Link]

MÉCANIQUE, ÉNERGÉTIQUE, GÉNIE CIVIL, ACOUSTIQUE M. Etienne PARIZET

INSA Lyon
[Link]
ED 162 Laboratoire LVA
Sec. : Philomène TRECOURT Bâtiment St. Exupéry
MEGA Tél : [Link].70 25 bis av. Jean Capelle
Bâtiment Direction INSA Lyon 69621 Villeurbanne CEDEX
mega@[Link] [Link]@[Link]

ScSo1 M. Bruno MILLY (INSA : J.Y. TOUSSAINT)

Univ. Lyon 2 Campus Berges du Rhône
ED 483 [Link] 18, quai Claude Bernard
Sec. : Mélina FAVETON 69365 LYON CEDEX 07
ScSo Tél : [Link].79 Bureau BEL 319
[Link]@[Link] [Link]@[Link]

Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
1
ScSo : Histoire, Géographie, Aménagement, Urbanisme, Archéologie, Science politique, Sociologie, Anthropologie
 6
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Résumé

Le système de refroidissement par plaques eutectiques (matériau à changement de phase,

MCP) est une possible alternative aux systèmes conventionnels de réfrigération alimen-
tés par des combustibles fossiles pour le transport de produits alimentaires surgelés ou
réfrigérés dans des remorques de camions. Des modèles numériques ont été dévelop-
pés pour évaluer sa faisabilité et ont été validés avec succès par rapport à des résultats
numériques et expérimentaux issus de la littérature.
Initialement, un modèle d’infiltration d’air dans la remorque a été développé utilisant un
modèle de turbulence k − ω SST pour prédire le comportement thermoaéraulique de l’air
pendant les périodes d’ouverture des portes. Différentes configurations des plaques eu-
tectiques et des ventilateurs ont été analysées. Sans cargaison, les plaques disposées en
série le long du plafond de la remorque ont montré un temps de renouvellement plus élevé
que celles disposées en parallèle à l’arrière, en raison des zones de recirculation. Cepen-
dant, une fois la cargaison introduite, les deux configurations offrent des performances
similaires car les zones de recirculation n’ont pas pu se former.
Par ailleurs, un modèle multiphasique granulaire eulérien-eulérien a été développé pour
prédire la formation et la croissance du givre sur les plaques eutectiques. Le modèle de
turbulence k − ω SST a été intégré pour étendre l’applicabilité du modèle de givre à une
gamme plus large de vitesses d’air. Un modèle de solidification et de fusion a également
été implémenté et couplé aux modèles précédents. Avec ce modèle combiné, les perfor-
mances d’un système eutectique ont été étudiées pour des conditions estivales typiques à
Montréal, Canada. Environ 2.3 % du MCP change effectivement de phase au cours des
120 premières secondes. Globalement, le système eutectique offre une alternative viable
au système conventionnel, bénéficiant des mécanismes de prévention des infiltrations ou
de dégivrage.

Mots-clés : Réfrigération, Matériau à changement de phase, Transfert de chaleur et de

masse, Dynamique des fluides numérique.

7
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 8
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Abstract

An eutectic refrigeration system is proposed to replace conventional refrigeration systems

powered by fossil fuel consumption for transporting frozen food products via refrigerated
truck trailers. To study the feasibility of the eutectic system, numerical models have been
first developed and favorably validated against published numerical and experimental re-
sults. The consolidated model has then been used to accurately predict the performance of
the eutectic plates, while being compared to new measurements on a lab scale system.
Initially, an infiltration model was developed using a k − ω SST turbulence model to
predict the thermoaerolic behavior of air during door-opening periods. Two fan modes and
different configurations of eutectic plates were analyzed. Without any cargo, plates placed
in series along the truck trailer’s roof had a higher renewal time than those placed in parallel
at the rear due to recirculation zones created by the pre-existing flow. However, as cargo
is introduced into the system, both configurations behaved similarly as recirculation zones
couldn’t form.
Furthermore, a granular multiphase Eulerian-Eulerian model was developed to predict
frost formation on top of a cold plate. The k − ω SST turbulence model was incorporated
into the frost model to extend its applicability to a wider range of inlet velocities. Finally,
the solidification and melting models were consolidated into the model as a user-defined
function (UDF). With the consolidated model, the performance of a eutectic system was
studied for typical summer conditions in Montreal, Canada. The analysis showed that ap-
proximately 2.3% of the phase change material undergoes a phase change within the first
120 seconds. Overall, the performance of the eutectic system can be considered a viable
replacement for the conventional system. Moreover, the system will be greatly enhanced
through infiltration prevention and frost removal mechanisms.

Keywords: Refrigeration, Phase Change Materials, Heat and Mass Transfer, Computational
Fluid Dynamics.

9
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 10
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Acknowledgement

I am immensely grateful to the many individuals and organizations whose unwavering sup-
port and guidance have paved the way for my doctoral journey.
First and foremost, I want to express my heartfelt appreciation to my supervisors, Pro-
fessor Sébastien Poncet and Professor Jocelyn Bonjour. Their expertise, encouragement,
and unwavering support have been the foundation of my research efforts. Their insightful
feedback, constructive criticism, and infinite patience have greatly influenced the direction
and quality of my dissertation.
I would also like to express my gratitude to my dissertation committee member, Dr.
Benoit Michel for his invaluable guidance, scholarly insights, and thoughtful critique through-
out the various stages of my research. Additionally, I would like to express my gratitude to
Professor Hachimi Fellouah for his guidance during my experiments.
This project has been funded by the NSERC chair on industrial energy efficiency, estab-
lished in 2019 at Université de Sherbrooke, with the support of Hydro-Québec (laboratoire
des technologies de l’énergie), Natural Resources Canada (CanmetEnergy-Varennes) and
Copeland Canada Inc. All calculations have been done using the HPC facilities of the Digital
Research Alliance of Canada. They are all here grateful acknowledged.
I extend my deepest gratitude to my family for their unwavering love, unshakeable belief
in my abilities, and constant encouragement during this challenging journey. Their sacrifices
and unwavering support have been a constant source of inspiration and motivation.

11
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 12
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Contents

1 Introduction 15
1.1 Context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.2 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.3 Novelty . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.4 Structure of the manuscript . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2 Literature Review 19
2.1 Refrigeration Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.1.1 Vapour Compression Systems . . . . . . . . . . . . . . . . . . . . . 19
2.1.2 Eutectic Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.1.3 Cryogenic Cooling System . . . . . . . . . . . . . . . . . . . . . . . 21
2.2 Phase Change Material . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.2.1 Solidification and Melting Model . . . . . . . . . . . . . . . . . . . 22
2.2.2 PU-PCM Truck Trailer Walls . . . . . . . . . . . . . . . . . . . . . 25
2.3 Frost Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.4 Infiltration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.4.1 Refrigerated Chamber . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.4.2 Refrigerated Truck Trailer . . . . . . . . . . . . . . . . . . . . . . . 32
2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3 Infiltration Model 39
3.1 Numerical modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.1.1 Geometrical modeling . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.1.2 Numerical method . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.1.3 Numerical parameters . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.2 Validation of the numerical model . . . . . . . . . . . . . . . . . . . . . . . 43
3.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.3.1 Configurations without cargo . . . . . . . . . . . . . . . . . . . . . 46
3.3.2 Influence of the cargo . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

4 Frost Model 57
4.1 Geometrical modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.2 Numerical modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.2.1 Numerical method . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.2.2 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.3 Grid independence and time sensitivity analysis . . . . . . . . . . . . . . . 61
4.4 Validation of the Frost model . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

13
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 5 Consolidated Model 69
5.1 Numerical modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.1.1 Geometrical modeling . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.1.2 Numerical method . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
5.1.3 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . 71
5.1.4 Solidification and melting . . . . . . . . . . . . . . . . . . . . . . . 71
5.2 Grid independence and time sensitivity analysis . . . . . . . . . . . . . . . 72
5.3 Validation of the UDF solidification and melting model . . . . . . . . . . . 72
5.4 Consolidated Model Validation . . . . . . . . . . . . . . . . . . . . . . . . 75
5.4.1 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.4.2 Experimental setup adjustment . . . . . . . . . . . . . . . . . . . . 79
5.5 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
5.6 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

6 Conclusions and future views 87

6.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.2 Future Works and Challenges . . . . . . . . . . . . . . . . . . . . . . . . . 90

7 Conclusions et perspectives (Version Française) 93

7.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
7.2 Perspectives de travail . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

14
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Chapter 1

Introduction

1.1 Context
Around 2 % of the total greenhouse gas emissions in developed countries are due to food
transportation, including motive power and refrigeration [1]. For road transport vehicles,
different types of refrigeration units are employed mainly depending on the size of the
vehicle and its load: vapour compression system with a vehicle or an auxiliary alternator
unit, a direct belt drive or an auxiliary diesel unit. For temperature-controlled vehicles,
legislations such as Agreement on the International Carriage of Perishable Foodstuffs
and on the Special Equipment to be used for such Carriage (ATP) [2] exist and provide
common standards and certifications. Road transport refrigeration equipment needs to
operate reliably in harsher environments and its operating conditions required to preserve
the cold chain, constraints due to available space, noise and weight makes them less effi-
cient than stationary systems. For medium to large road transport vehicles, refrigeration
units are primarily powered by a self-contained diesel engine. To guarantee −20 °C within
the trailer, the fuel consumption is typically up to 4 liters per hour for a semi-truck trailer
with an inside volume of 78.79 m3 [1].
Eutectic systems typically consist of tubes, beams or plates filled with a Phase Change
Material (PCM) and appear as an alternative to reduce or annihilate this fuel consump-
tion. An example of eutectic refrigeration system inside a refrigerated truck trailer is
shown in Figure 1.1. During the typical delivery cycle, the eutectic plates absorb heat in-
filtrated into the trailer as latent heat and produce cooling to ensure the desired regulated
temperature inside the trailer. However, heat losses through frequent door openings are a
major concern for a medium sized truck trailers with day-deliveries, especially in humid
climates. Frequent door openings account for the majority of the heat infiltration during
delivery and induce condensation and formation of ice along the plates [3; 4]. Thick con-
densed frost layers developed on the eutectic plates may interfere with the heat transfer
with the surrounding air and reduce its effects.

1.2 Objectives
The main objective of this manuscript is to analyze the feasibility of the eutectic refrig-
eration system to replace the conventional refrigeration system powered by fossil fuels
for the transport of frozen foodstuff. To analyze the feasibility of the eutectic system,
overall performance of the system has to be accurately predicted. For the prediction of
the overall system performance, the following analyses are required:

15
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Figure 1.1: Example of an eutectic refrigeration system inside a refrigerated truck trailer.
Courtesy of FrygyCube.

• Analysis of the thermoaeraulic behavior of the air during the door opening period
for a refrigerated truck trailer equipped with eutectic plates. Parameters such as
the renewal time, temperature distribution of the truck trailer, and the maximum
temperature observed by the cargo are few key parameters of the analysis.

• Prediction of the frost development on top of a cold plate under a wide range
of initial and boundary conditions. As the infiltrating hot atmospheric air comes
in contact with the cold plate, frost development is inevitable. However, most
published numerical models has only been validated for very specific ranges of initial
and boundary conditions. Therefore, a more generalized numerical model is required
for a wider range of initial and boundary conditions.

• Prediction of the solidification and melting of the PCM. The eutectic refrigeration
system is subjected to the heat infiltrated during the door opening periods of the
refrigerated truck trailer. During the door opening period, eutectic mixture inside
the eutectic refrigeration system will undergo a phase change. As the heat is not
uniformly distributed throughout the eutectic refrigeration system, there will be an
inhomogeneous phase change of the eutectic mixture. Therefore, a numerical model
is required to predict the solidification and melting of the eutectic mixture during
the door opening periods.

1.3 Novelty
One original aspect of this work lies with the choice of identifying the feasibility of the
sub-zero eutectic system with a substantially low phase change temperature of T =−29 ◦ C
for the transport of frozen foodstuff. Due to the phase change temperature of the eutectic
system, frost development occurs during the door opening period as the hot humid air
comes in contact. Frost formed on top of the eutectic system interferes with the heat
transfer between the eutectic mixture and the surrounding air. Additionally, as the frost

16
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 formation is an exothermic process, it will release heat to the eutectic system and attribute
towards the phase change of the system.
Furthermore, the analysis of the thermoaeraulic behavior of the air during the door
opening of a refrigeration truck trailer equipped with eutectic plates includes two different
fan operating modes and configurations of the eutectic plates. Presence of the cargo is
also introduced to the system to analyze the temperature distribution of the cargo during
the infiltration period to confirm it meets the regulatory temperature imposed by the
ATP [2]. Additionally, the complete analysis of the frost development, solidification and
melting of the PCM, and the thermoaeraulic behavior of the air for a refrigerated eutectic
system has never been studied experimentally nor numerically in detail.

1.4 Structure of the manuscript

To achieve the objectives listed in Section 1.2, this manuscript is organized as follows:

1. Literature Review (Chapter 2): summarization of the state of the art of the re-
searches related to the current objectives of the manuscript.

2. Infiltration model (Chapter 3): a model developed to accurately predict the infil-
tration behavior of the air during the door opening period of a refrigerated truck
trailer equipped with eutectic plates.

3. Frost model (Chapter 4): an extended model of the works of Wu et al. [5] to predict
the formation of the frost developed on top of a cold plate with its interaction with
the air flow above the frost.

4. Consolidated model (Chapter 5): a model consolidating the aforementioned infiltra-

tion and frost model with the solidification and melting model to predict the overall
performance of the eutectic plate system.

5. Conclusions and Future Works and Challenges (Chapters 6 in English and 7 in

French): summarization of the objectives achieved within this manuscript. It con-
tains also the possible future works that could be performed and the challenges
associated with the current and future works.

17
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 18
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Chapter 2

Literature Review

This chapter presents the results and methodologies employed by the researchers to in-
vestigate the key components (infiltration, frost development, solidification and melting)
required to develop the comprehensive Computational Fluid Dynamics (CFD) model for
a refrigerated truck trailer equipped with eutectic plates. Refrigeration systems currently
employed for a refrigerated truck trailer are described in Section 2.1. Then the infiltra-
tion dynamics during the door opening for a refrigerated chamber and a refrigerated truck
trailer are presented in Section 2.4. Finally, the physical phenomena involved with the
eutectic refrigeration systems, frosting and solidification and melting of the PCM, are
presented in Section 2.2.

2.1 Refrigeration Systems

2.1.1 Vapour Compression Systems
The majority of refrigeration systems currently employed for a refrigerated truck trailer
are the vapour compression systems [3]. Depending on the factors such as weight, noise,
maintenance, environmental concerns, and regulation, an appropriate vapour compression
system is selected. Nevertheless, the Coefficient of Performance (COP) of the refrigeration
systems for a truck trailer is relatively low. The value of COP typically ranges from 0.5
for −20 ◦ C inner temperature and 1.5 ∼ 1.75 for 3 ◦ C inner temperature with the ambient
temperature of 30 ◦ C [3]. Four vapour compression systems most commonly employed
for a refrigerated transport are [3; 4]: vehicle alternator unit, direct belt drive, auxiliary
alternator unit, and auxiliary diesel unit (Figure 2.1).
Vehicle alternator unit is commonly employed for a small delivery vans. The vehicle
engine crankshaft is used to drive an upgraded single alternator. Then, the alternator
charges the battery and powers the refrigeration system. Direct belt drive is mostly
used for a van-sized vehicles. The vehicle engine directly drives the compressor of the
refrigeration unit through a belt. Auxiliary alternator unit uses the main traction engine
to drive a dedicated large alternator by a belt. The power generated from the main
engine is used to drive the electric motor in the refrigeration unit. Auxiliary diesel unit
is a system with a diesel engine incorporated into the refrigeration unit. This system is
employed in the majority of the medium to large refrigerated vehicles [3; 1].
The vehicle alternator unit, direct belt drive, and auxiliary alternator unit provide
limited power to the refrigeration system as the main source of power is from the main
engine. Also the power losses through mechanical transmission by a crankshaft or a

19
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Figure 2.1: Schematic diagram of vapour compression transport refrigeration unit driven by
a diesel engine [6].

belt contributes to their limited usage to a small to medium sized transport vehicle.
Also, to meet the requirements of the ATP agreement [2] and operate over wide range of
conditions, vapour compression refrigeration systems needs to be oversized by up to 1.75
time its calculated load [3]. Therefore, auxiliary diesel engines accounts for approximately
90 % of the market for the transport refrigeration system [6].
While the auxiliary diesel engine is capable of reliably providing power for a medium
to large sized transport vehicle refrigeration system, it can also account up to 40 % of the
environmental impact of the vehicle engine [3]. To reduce the environmental impact of
the transport refrigeration systems, eutectic refrigeration systems are employed.

2.1.2 Eutectic Systems

Eutectic systems consists of tubes, beams or plates filled with a Phase Change Material
(PCM) to store energy as latent heat and provide a cooling effect to maintain the desired
temperature inside a transport vehicle (Figure 1.1). The main difference to the conven-
tional system (vapour compression) is that the cold source is provided by the phase change
of the material rather than the direct expansion of the refrigerant gas [3]. The eutectic
systems operate silently and provide cooling for a specific duration before the PCM is
completely melt. Also, static eutectic refrigeration units are capable of charging the PCM
by a connector, eliminating the need to have its own condensing unit and leading to a
decrease in the overall cost of the system [3].
However, for a medium-sized vehicle with multiple day-deliveries, heat losses through
frequent door openings is a major concern. The majority of the heat infiltration during
delivery is due to door openings and it can also induce frost and ice along the eutectic
system [3; 4]. Frost and ice formed on the eutectic system could interfere with the heat
transfer between the system and the surrounding air. Also, depending on the length of
the delivery, frequency of the door openings, and time period, PCM inside the eutectic
systems could also completely melt before the end of its delivery cycle. Therefore, more
thorough analysis of the eutectic transport refrigeration system is required.

20
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 2.1.3 Cryogenic Cooling System
Alternative to the eutectic refrigeration system, total loss systems using cryogenic fluids
can also be employed. The cryogenic fluids are stored in tanks and connected to a spraying
system through the length of the vehicle [3]. Then, it is released into the vehicle and
vaporises to reduce the temperature inside. Cryogenic cooling system is capable of rapidly
decreasing the temperature inside the vehicle with very low noise emissions. However, it
is rather expensive to operate.

2.2 Phase Change Material

Latent heat storage provided by the PCMs is one of the most efficient methods to store
thermal energy [7]. PCM takes advantage of the material’s phase transition to store
or release energy due to the enthalpy change caused by the phase change. Due to the
phase change, the thermal energy storage operate almost isothermally at the phase change
temperature of the material [7]. Initially, as the temperature increases or decreases, energy
is stored or released by the temperature change of the storage material, acting as the
sensible heat storage [7; 8]. The amount of sensible thermal energy stored is defined as:
Z T2
Qsens = m Cp ∆T (2.1)
T1

where m is the mass of the storage material, Cp its specific heat, T1 and T2 are the initial
and final temperatures, respectively.
As the temperature reaches the phase change temperature of the substance, the mate-
rial consumes significant quantities of energy as latent heat. The amount of latent thermal
energy stored is defined as:
QL = m∆hL (2.2)
with ∆hL the phase change enthalpy (Latent heat of fusion/solidification).
After the complete phase change of the substance, the energy is stored once again in
sensible fashion as shown in Figure 2.2.

Figure 2.2: A schematic of sensible and latent heat storage principle [7].

Typically, PCM classifications are divided into four categories: solid-solid, solid-liquid,
solid-gas, and liquid-gas [7; 8; 9; 10]. From these four types of PCMs, solid-liquid PCMs

21
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Figure 2.3: Overview of PCMs [10].

are the most suitable for thermal energy storage, and they are further classified as organic,
inorganic and eutectic as shown in Figure 2.3.
Paraffin is typically composed of long straight chain of alkanes and the melting point
depends on the length of the chain [7]. It has reasonably high heat of fusion and can be
used over a wide range of temperature [10]. Fatty acid typically has higher heat of fusion
compared to the paraffin wax with the ability to reproduce melting and freezing with little
to no supercooling [10]. However, fatty acids can cost up to 2.0 to 2.5 times more than
the cost of paraffin wax [10]. Salt hydrates are one of the most extensively studied PCM
due to their moderately high thermal conductivity, cost and safety [7]. However, the salt
hydrates exhibits incongruent melting when the salts are insoluble to its associated water,
leaving the salt solution supersaturated after the phase transition [7]. Metallic PCMs have
not been seriously considered due to their heaviness [10], but have gained some attraction
recently due to some of their favourable characteristics such as high thermal conductivity,
heat of fusion per unit volume and small volume change associated with phase change
[7].

2.2.1 Solidification and Melting Model

For the latent heat storage systems such as PCMs, melting and solidification are the two
processes generally involved. Heat is transferred to a solid PCM and the melting occurs
as the latent heat is stored. Reversely, when the latent heat is released, solidification
occurs and the heat is transferred from the PCM. Therefore, the PCM transition involves
a solid, a mushy, and then a liquid state. The heat transfer associated can be either by
conduction or convection (mostly by natural convection) [11].
The solidification and melting model based on the enthalpy porosity method of Brent
et al. [12] is already readily available and well established in ANSYS FLUENT [13]. Zara-
jabad and Ahmadi [14] numerically investigated the effective PCM volume to be employed
for a cold storage by the solidification and melting model. Zarajabad and Ahmadi [14]
employed an eutectic solution of N aCl − H2 O as the PCM with −21◦ C melting tem-
perature. Gou et al. [15] numerically simulated the cold energy release process of a cold
storage plate in a container for temperature control by employing the Volume of Fluid
(VOF) multiphase flow model and the solidification and melting model in ANSYS FLU-
ENT. The impacts of the inlet velocity of 1.8 m s−1 to 4.2 m s−1 and the spacing of the
cold storage plates of 10, 20 and 30 mm on the release of the cold energy were investi-
gated. Calati et al. [16] employed the solidification melting model in ANSYS FLUENT
to numerically investigate the effects of the PCM system for food transportation with the

22
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 melting range of PCM being between 1 to 3 ◦ C.
The enthalpy-porosity approach proposed by Voller et al. [17] and Brent et al. [12] has
been shown to adequately describe the natural convection effect in the melting region.
It is capable of accurately predicting the transient position and shape of the melt front.
In this model, the fraction of the cell volume that is in liquid form is computed at each
iteration based on an enthalpy balance. Liquid fraction of 0 indicating a completely solid
state and 1 indicating a completely liquid state. The mushy zone is modelled as a porous
medium where the porosity decrease from 1 to 0 as the solidification occurs. As the
solidification process complete, the porosity becomes zero and so does the velocities.
The enthalpy-porosity model based on the energy equation is described below [17; 12]:

∂
(ρh) + ∇ · (ρvh) = ∇ · (λ∇T ) + S (2.3)
∂t
with h the enthalpy, ρ the density, v the fluid velocity, λ the thermal conductivity and S
the source term. Z T
h = href + Cp dT + ϕ∆hl (2.4)
Tref

Enthalpy considers both the sum of sensible enthalpy and the latent heat of fusion (hl ),
with ϕ the liquid fraction. The liquid fraction is identified by:

0, T < Ts
T −Ts
ϕ = Tl −Ts , Ts < T < Tl (2.5)
0, T > Tl


with Ts and Tl the solidus and liquidus temperatures.

The momentum sink is also applied for the reduced porosity in the mushy zone as
following:
(1 − ϕ)2
Su = 3 Amush v (2.6)
ϕ +χ
where χ is a small number (default value of 0.001) to prevent division by zero, Amush the
mushy zone constant, v the fluid velocity and subscript u refers to the momentum.
The solidification and melting model based on the enthalpy-porosity approach has
been employed by the researchers for the PCM phase changes in various applications.
The described models in Table 2.1 focuses on the rectangular geometry, to be applied for
the eutectic plates inside a refrigerated truck trailer.
All of the aforementioned models were able to accurately predict the solidification and
melting of the PCM using the existing "Solidification and Melting" model in FLUENT.
However, the obtained results depended heavily on the value of the mush zone parameter,
Amush . The mushy zone parameter (Amush ) is set to the default value of 105 with the rec-
ommended operating range between 103 and 108 . The mushy zone parameter determines
the amplitude of the damping, the higher the value, the faster the velocity of the material
becomes zero as it solidifies. Furthermore, the mushy zone significantly influences the
heat transfer and flow characteristics during the solidification and melting process. Fadl
and Eames [11] investigated the effect of the value of the mushy zone parameters and
compared the results against the experimental results of Shokouhmand and Kamkari [22].
In Figure 2.4, predicted melt front locations with different values of Amush obtained
by Fadl and Eames [11] is presented. From the obtained results, resulting shape of the
solid PCM for Amush = 5 × 105 presented with the most agreeable results against the

23
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Figure 2.4: Experimental and predicted melt front locations with time for different values of
Amush indicated on the top [11].

24
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Author(s) PCM Mesh Time Material
step Property
Vogel et al. [18] Eutectic Non-uniform 0.0125 s Boussinesq
mixture wall refinement (Constant)
Kim et al. [19] Gallium Not Not Temperature-
Available Available Varying
Fadl and Eames [11] Lauric Uniform 0.2 s Not
acid 1 mm × 1 mm Specified
Oliveski et al. [20] Lauric Non-uniform 0.1 s Temperature-
acid Acell ≈ 1.45 mm2 Varying
Abdulmunem et Paraffin Uniform 0.5 s Temperature-
al. [21] wax 0.0625 mm × Varying
0.25 mm
Zarajabad and Ahamdi N aCl − Quad element 0.1 s Boussinesq
[14] H2 O (Constant)
Gou et al. [15] Water Rectangular mesh Not Boussinesq
8.3 × 105 elements Available

Table 2.1: Summary of recent investigations on rectangular PCM solidification and melting
problems.

experimental results. However, as the value of Amush increased, the numerical predictions
of the melt front lagged behind the experimental results and had a sharper curvature
along the melt interface. On the contrary, as the value of Amush decreased, the melt front
exceed that of the experimental results and overpredicted the melting process. However,
Fadl and Eames [11] noted that the extent of correlation between the values of Amush and
the material properties are unknown and an approach to identify the optimal value of the
mush parameter does not exists as of yet.

2.2.2 PU-PCM Truck Trailer Walls

In the recent years, application of PCMs for the walls of refrigerated vehicles have been
proposed and investigated by the researchers [3; 1; 23; 24; 25]. Tinti et al. [23] investigated
the characteristics of rigid polyurethane (PU) foams typically used as refrigerated truck
trailer walls, containing MicroPCM. They found that with increasing PCM content, the
thermo-regulating and thermal energy capabilities of PU-PCM foams increased. Coper-
taro et al. [24] and Michel et al. [25] investigated the thermal performance of a sandwich
structured PU-PCM insulation walls for the refrigerated vehicles. The effect of the design
of the PU-PCM walls (PU/PCM layer thickness, PCM mixture, ...) on their thermal
performances over 24 hours was analyzed. However, none of the aforementioned studies
analyzed the heat transfer inside the refrigerated vehicle. Most recently, Calati et al. [26]
numerically investigated the performance of the PU-PCM walls for a refrigerated truck
trailer, considering the solar radiation from 6 AM to 4 PM of a typical summer day in
Vicenza (Italy). They simulated the solidification and melting of the PCM by using the
existing package in FLUENT. While the heat transfer inside a refrigerated vehicle was
considered, the analysis was carried out for a closed configuration.

25
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Figure 2.5: Schematic diagram of the physical processed involved in frost development [29].

2.3 Frost Model

To the best of the author’s knowledge, there has yet to be an investigation on the frost
development on the PCM. Berdja et al. [27] experimentally and analytically investigated
the characteristics and thickness effect of PCM and frost for a conventional refrigerator.
However, their analytical model did not consider the formation of frost on top of the PCM.
But, rather incorporated the thermal resistance of the frost to the overall heat transfer
coefficient and investigated the influence of the frost on the heat transfer depending on
the varying thickness. Therefore, comprehensive models capable of simulating the devel-
opment of frost on the PCMs (eutectic plates) as well as the solidification and melting of
the PCMs (described in Section 5.1.4) are required.
Frost is a thin layer of ice formed on top of a solid surface with a below freezing tem-
perature when the water vapor in an above-freezing condition comes in contact. Figure 2.5
presents generalized schematic of the frost formation on a cold plate. Hayashi et al. [28]
divided the frost formation into three steps: the crystal growth period, the frost layer
growth period and the frost layer full growth period (densification). The frost formation
steps are detailed below [28]:
1. Crystal growth period: initially, a thin layer of frost covers the surface. Then the
frost crystals located relatively far apart from each generate on the thin frost surface
and grow in vertical direction at an approximately same rate. This crystal growth
is approximately characterized by a growth in a linear direction.
2. Frost layer growth period: Rough frost (cluster of rod type crystals) changes in
shape by generating branches around the top of a crystal or by the interaction of
each crystal. Then it grows gradually into a meshed and progressively more uniform
structure until the frost surface becomes nearly flat. During this period, density of
the frost increases with the frost layer growth.
3. Frost layer full growth period: the shape of the frost layer does not change until
the surface temperature becomes 0 ◦ C due to its increased thermal resistance. Then
the melted water from the frost surface soaks into the frost layer and freezes. This
phenomenon increases the frost layer density and the frost deposition occurs again.
This cycle continues periodically until the equilibrium condition of the heat transfer
is reached and the frost formation stops.
In the last decades, there has been extensive research on the development of frost,
both experimentally and theoretically. Hayashi et al. [28] and Tao et al. [30] modeled the

26
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 fully developed growth period as a homogeneous porous medium with diffusion leading
to the frost growth and densification. Jones and Parker [31] developed a model based on
the mass and energy balances and diffusion of the water vapour molecules at the frost
interface. O’Neal and Tree [32] calculated the mass transfer by modeling the diffusion
of molecules in porous media and proposed frost thickness correlation based on their
experimental results. Then, Na and Webb [33], Lee et al. [34], and Yao et al. [35] further
improved upon the model of O’Neal and Tree. The interested readers are referred to the
work of Léoni et al. [36] for further details on the empirical and theoretical models for
the frost formation developed over the years.
Frost formation models can be classified into three different categories.

1. Frost formation without the consideration of the influence of the air flow [33; 31; 37;
38; 39]: the heat and mass transfers were calculated using the empirical correlation
for the diffusion on the air side. Na and Webb [33; 37] assumed the humid air to be
supersaturated, while Barron et al. [39], Jones and Parker [31], and Lee et al. [38]
assumed that the humid air was saturated near the frost surface.

2. Frost formation with consideration of both the frost region and the humid air region.
The interface condition was employed to connect the two regions [34; 40; 41; 42].
Lee et al. [34], Yang and Lee [41], and Yang et al. [42] calculated the frost growth
on a cold plate with laminar and turbulent flow with this model.

3. Frost formation modeled by considering the mass transfer between the water vapour
within the humid air to the ice. Multi-phase flow methods are employed to solve the
governing equations for the humid air and the ice phases. Frosting processes were
simulated by adding the mass transfer source terms to the governing equations. Cui
et al. [43; 44] proposed a frosting mass transfer source term based on the nucleation
theory to model the frost formation on the flat plate and the fin-and-tube heat
exchanger surfaces. Zhaung et al. [45] implemented the mass transfer source term
model to simulate the condensing droplet formation. Wu et al. [5; 46] proposed
a frosting mass transfer source term by assuming the phase change driving force
from the water vapour to the ice as the difference between the water vapour partial
pressure in the humid air and the water vapour saturation pressure based on the frost
surface temperature. Wu et al. [5; 46] used the model to simulate the frost formation
on a flat plate with a local cooling and fin-and-tube heat exchanger surfaces and
obtained good agreement with the experimental results. Afrasiabian et al. [47] later
extended the model of Wu et al. [46] by separating the velocity and the water vapour
concentration parameters of the frosting condition to simulate the frost growth on
one channel of a plate-fin evaporator.

However, the models belonging to the first and the second categories determine the
thickness and the weight of the frost without considering the interaction between the
humid air flow field change and the frost growth [46]. The frosting model based on the mass
transfer source terms can simulate the interaction between the humid air flow field change
and the frost growth, predicting the localized frost and humid air flow characteristics
(Figure 2.6). Therefore, among the three categories of frost formation models, the ones
with the mass transfer during the frost formation contributing to the increase in both
the frost thickness and density showed best agreement with the experimental frost data
[5; 46].

27
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Figure 2.6: Schematic diagram for the computing domain in the frosting process [46].

Cui et al. [43; 44], Wu et al. [5; 46], and Afrasiabian et al. [47] employed a Eulerian-
Eulerian multiphase model with the mass transfer source terms between the primary
phase of humid air and the secondary phase of ice. Based on the mass transfer, the
momentum and energy transfers from water vapor to ice can occur and are reflected as
the momentum and energy source terms of the governing equations. Also a source term
is included in the species conservation equation for the water vapour diffusion driven by
the concentration gradient from the humid air to the frost layer. Therefore, two sets of
conservation equations for mass (Equation 2.7, 2.8), momentum (Equation 2.9, 2.10),
energy (Equation 2.11, 2.12), and a species conservation equation (Equation 2.13) were
solved:
∂
(αc ρc ) + ∇ · (αc ρc uc ) = Smc (2.7)
∂t
∂
(αa ρa ) + ∇ · (αa ρa ua ) = Sma (2.8)
∂t
∂
(αc ρc uc ) + ∇ · (αc ρc uc uc ) = Kac (ua − uc ) − αc ∇P + ∇ · τ¯c + αc ρc g + Suc (2.9)
∂t
∂
(αa ρa ua ) + ∇ · (αa ρa ua ua ) = Kca (uc − ua ) − αa ∇P + ∇ · τ¯a + αa ρa g + Sua (2.10)
∂t
∂ ∂pc
(αc ρc hc ) + ∇ · (αc ρc uc hc ) = Q̄ac + αc + τ¯c : ∇uc − ∇ · q⃗c + Shc (2.11)
∂t ∂t
∂ ∂pa
(αa ρa ha ) + ∇ · (αa ρa ua ha ) = Q̄ca + αa + τ¯a : ∇ua − ∇ · q⃗a + Sha (2.12)
∂t ∂t
∂
(αa ρa wv ) + ∇ · (αa ρa wv ua ) = ∇ · (ρa DH2 O ∇wv ) + Sva (2.13)
∂t
with α the volume fraction, ρ the density, u the velocity, K the momentum transfer
coefficient, g the gravity, τ̄ the stress-strain tensor, Q̄ the interphase heat transfer, ⃗q the
heat flux, w the mass fraction, and DH2 O the diffusivity of the water vapour.
The mass (Equations 2.14, 2.15), momentum (Equations 2.16, 2.17), energy (Equa-
tions 2.19, 2.18) and species(Equations 2.20) sources employed for the defined mass trans-
fer rate (ṁac ) are listed below:

Smc = ṁac − ṁca = ṁac (2.14)

Sma = −ṁac (2.15)

Suc = ṁac uc (2.16)
Sua = −ṁac uc (2.17)

28
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Sha = ṁac hva (2.18)

Shc = ṁac (hva + ∆hL,va ) (2.19)

Sva = −ṁac (2.20)

Author(s) Mesh Time Method Simulation

step parameters
Cui et al. [43; 44] Non-uniform, 1s Nucleation 0.6 ≤ uin ≤ 2.5 m s−1
wall refinement Tin = 19.8°C
nx × ny = 300 × 30 Tplate = −20.5°C
RH = 58%
Wu et al. [5] Uniform, 1.0 s Phase uin = 1.5 m s−1
1 mm × 0.05 mm change Tin = 2 ◦ C
driving Tplate = −9 ◦ C
force RH = 83 %
Wu et al. [46] Uniform, 1s Phase 0.31 ≤ uin ≤
0.1 mm × 0.05 mm change 0.92 m s−1
driving −5 ≤ Tin ≤ 21.4 ◦ C
force −20 ≤ Tplate ≤ −5 ◦ C
38 ≤ RH ≤ 85 %
Afrasiabian et Not Available 1s Phase uin = 0.7 m s−1
al. [47] change 1.7 ≤ Tin ≤ 22 ◦ C
driving −16 ≤ Tplate ≤ −4 ◦ C
force 80 ≤ RH ≤ 85 %

Table 2.2: Summary of the frost model based on the mass transfer rate.

Cui et al. [43; 44] defined the mass transfer rate, ṁac , based on the classical homoge-
neous nucleation theory as:

4 3 δr̄
ṁac = πρc Ir∗ + 4πρc ηr̄2 (2.21)
3 δt
with I the rate of nucleation per unit volume, ρi the density of ice, r∗ the Kelvin-Helmholtz
critical radius, η number of droplets per unit volume, r̄ droplet radius, δr̄δt
the nucleation
growth rate.
Wu et al. [5] defined the mass transfer rate based on the phase change driving force
as: 
τ α ρ (w − wvs ), wva − wvs > Bua
ṁac = v a a va (2.22)
0, otherwise
Where, τv is the time relaxation coefficient, αa the volume fraction of humid air, ρa the
density of humid air, wva mass fraction of water vapour, and wvs saturated water vapour
mass fraction.
Then Wu al. [46] modified the correlation by replacing the phase change driving force
(wva − wvs ) with wva wvaw−w
vs
vs
as following:

τv αa ρa wva wvaw−w , wvaw−w ≥ B uuina

 vs vs
ṁac = vs vs (2.23)
0, otherwise

29
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Finally, Afrasiabian et al. [47] modified the frosting criteria by separating the velocity
component from the water vapour concentration:

τv αa ρa wva wvaw−w , wva − wvs > 0 and uuina ≤ B

 vs
ṁac = vs (2.24)
0, otherwise

Afrasiabian et al. [47] reported a reduced sensitivity of the mass transfer to the value of
the non-dimensional frosting coefficient, B.
Compared to the experimental results of Lenic et al. [40], numerical results of Cui et
al. [43] showed a satisfactory agreement with a maximum error of approximately 13 % for
the frost layer thickness. The results of Afrasiabian et al. [47] demonstrated a maximum
difference of 7 % against the experimental data of Sommmer et al. [48] after 30minutes for
the average frost thickness and density. However, the frosting criteria parameter, B, is not
defined by Afrasiabian et al. [47]. While the non-dimensional coefficient, B, from Table 2.2
is developed only for the simulation conditions of the cooling surface temperature between
−20 ◦ C and −5 ◦ C, air velocity between 0.31 and 0.92 m s−1 , and the water vapour mass
fraction between 0.002 and 0.006 by Wu et al. [46]. Cui et al. [43; 44], Wu et al. [5; 46], and
Afrasiabian et al. [47] performed their simulations using FLUENT with the User Defined
Function (UDF) DEFINE_MASS_TRANSFER to define the mass transfer source term.
Table 2.2 presents the simulation parameters employed for the frost models based on the
mass transfer rate.

2.4 Infiltration
2.4.1 Refrigerated Chamber
The heat and mass transfers inside a refrigerated truck trailers exhibit some similarities
with those in a refrigerated chamber. However, the heat and mass transfers inside a refrig-
erated chamber have received significantly more attention by the researchers. Therefore,
infiltration dynamic for a refrigerated chamber is firstly introduced in this section. Dur-
ing the door opening period for a refrigerated chamber, numerically, infiltration does not
present with constant infiltration with time. Rather, infiltration can be split into three
separate stages as shown in Figure 2.7: lag, steady state, and tail off stages [49]. Lag
stage corresponds to the initial few seconds after the door opening when the flow takes
some time to fully develop. It is then proceeded by the steady-state stage where there is a
constant flow rate through the entrance. Then finally, tail off stage where the temperature
difference (driving force) between the cold store and the surroundings is reduced.

Figure 2.7: Predicted infiltration for the 2.3 m wide entrance for different door opening
times [49].

30
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Figure 2.8: Deformation of temperature (left) and velocity (right) profiles with respect to the
exterior temperature [50].

For the velocity and temperature profiles during the door opening period at the door-
way vertical centerline, there exists a "neutral point" where the internal and external
potential pressure are equal and the mean velocity between the infiltration and exfiltra-
tion is null [50]. This point is typically situated approximately at the location y ≈ H/3
(Figure 2.8). Where H corresponds to the height of the doorway. In the region near the
neutral point, there is a continuous mixing and swirling of the cold and the hot air present.
But as the velocity magnitude in this region is minute, it is challenging to identify the
directions of the movement of the air in this region.
In the region above the neutral point, hot exterior air enters into the chamber with a
velocity that increases when y increases from the neutral point. Similarly, cold air from
the chamber exits through the zone below the neutral point with a velocity that increases
when y decreases from the neutral point. The velocity of the air strongly correlates
to the temperature difference between the hot atmospheric air and the cold air inside
the chamber. Another main flow characteristic is that the flow variables vary mainly in
the vertical direction, but exhibits small to negligible differences horizontally [50]. This
behaviour is used for the justification of a 2D model to accurately simulate the infiltration
during the door opening period by Croquer et al. [4], and Bonaventure et al. [51].
Azzouz et al. [50] experimentally determined the cold losses caused by the door open-
ings for a refrigerated chamber with an inner volume of approximately 13 000 m3 with
four identical doors of 3 m in height and 2.8 m in length. The speed of air during the door
opening period varied between 0 to 1.5 m s−1 with the temperature ranging between −25
and 25 ◦ C. The velocity field reached steady state 3.0 s after the door opening, while the
temperature field reached steady state more than 13.0 s after the door opening.
Foster et al. [49] experimentally and numerically investigated the heat and mass trans-
fers between a refrigerated store with an inner volume of approximately 106 m3 with a
large single opening of 7.36 m2 . Temperature inside the refrigerated store was set at
−20 ◦ C with the external temperature of 20 ◦ C. Before the door opening period, the
evaporator fans were switched off to allow for the air movement to settle for 30 s. Then
the infiltration was measured for the door opening times of 10, 20, 30, and 40 s. The infil-
tration rate increased linearly with the opening time until reaching an asymptotic value
corresponding approximately to 90 % of the inside store volume at t ≈ 18 s. Foster et
al. [49] calculated infiltration rate based on the concentration of CO2 before the door
opening and after closing the door.
Furthermore, many expressions for the steady-state heat infiltration rate based on the

31
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 ideal flow theory for natural convection of fluids at different densities through openings
have been developed over the years [52; 53; 54; 55; 56; 57].

Author(s) Infiltration Expression [m3 ]

−ρo 0.5
Brown and Solvason [53] I = 0.343A(gH)0.5 [ ρρiavg ] [(1 − 0.498) Hb ]
Tamm [54] I = 0.333A(gH)0.5 [ ρi −ρ
ρi
] [ 1+(ρo /ρ2 i )0.333 ]1.5
o 0.5

Fritzsche and Lilienblum [55] I = IT amm Kf,L , Kf,L = 0.48 + 0.004(To − Ti )

Gosney and Olama [56] I = 0.221A(gH)0.5 [ ρi −ρ
ρi
] [ 1+(ρi /ρ2 o )0.333 ]1.5
o 0.5

Pham and Oliver [57] I = 0.226A(gH)0.5 [ ρi −ρ

ρi
] [ 1+(ρo /ρ2 i )0.333 ]1.5
o 0.5

Table 2.3: Comparison of the expressions for the infiltration rate used in various analytical
models.

In Table 2.3, various expressions for the infiltration through a doorway by natural
convection are summarized. Brown and Solvason [53] assumed that the height of the neu-
tral level is half the height of the doorway, attributing to the substantial over-prediction
(52.1 to 122.7 %) of the infiltration rate compared to the experimental results of Foster
et al. [49]. Whereas Tamm [54] improved upon the model by calculating the height of the
neutral level and replacing ρavg with ρi and resulted in better prediction of the infiltration
rate [49]. Fritsche and Lilienblum [55] added a correction factor (Kf,L ) accounting for the
contraction of the flow, friction and thermal effects based on the experiments performed
using vane anemometers. It assumed that the infiltration volume flow rate is equal to the
exfiltration volume flow rate. However, this is only if the air entering into the refrigerated
chamber does not cool. Therefore, Gosney and Olama [56] instead provided an equation
for the constant mass flow rate (as the mass of air inside refrigerated chamber remains con-
stant) by replacing the term (ρo /ρi ) with (ρi /ρo ). Finally, Pham and Oliver [57] modified
Tamm’s equation based on their experiments.
However, it was identified that these analytical models based on the ideal flow theory
generally over-predicts the experimental results [58; 59]. Foster et al. [49] observed that
the Gosney model was the best analytical model from Table 2.3, with the model accurately
predicting the infiltration within the experimental error for two of their experiments.
Interestingly, their CFD predictions showed significant improvement in accuracy over the
fundamental analytical equations of Brown and Tamm, but less accurate to those of
Gosney and Fritzsche. The authors suggested that a more detailed CFD model would
likely yield a better results, but it could not be done due to limited computing resources.

2.4.2 Refrigerated Truck Trailer

While the heat and mass transfer inside a refrigerated truck trailer exhibit some similari-
ties with those in a refrigerated chamber, it is not possible to simply apply the analytical
models for the chambers mentioned in Section 2.4.1 such as the Gosney model directly
to the refrigerated truck trailer. The main cause is due to the significant difference of
the opening area to enclosure volume ratio between the refrigerated truck trailer and the
refrigerated chamber. Therefore, heat and mass transfer inside a refrigerated truck trailer
requires separate analysis from the refrigerated chamber. However, as there has not been
as much attention to the refrigerated truck trailers than for the refrigerated chambers, it
is still a topic of interest for the researchers.

32
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Figure 2.9: Evolution of the infiltration flow rate for the baseline case [60].

For the refrigerated chamber, infiltration is solely driven by the temperature difference
between the cold storage and the surrounding air. However, Lafaye de Micheaux et al. [60]
observed that during the door opening period for a truck trailer, two different mechanisms
were observed for the heat and mass exchanges (Figure 2.9):

1. The volume flow rate increases progressively while the flow pattern develops over the
doorway. Then, an unsteady mass exchange ("buoyancy-driven flow" [60]) occurs
driven by the density differences between the inner and external air volumes. The
volume flow rate reaches its peak value then quickly decreases as the temperature
of the inner volume increases, effectively reducing the driving force for the unsteady
mass exchange.

2. Heat exchange between the outside air and the inner walls ("boundary layer flow" [60]).
As the inner volume is filled with the atmospheric air, natural convection over the
cold walls of the container drives the flow inside the trailer. This phenomena was ob-
served roughly when the inside air has reached a ceiling temperature of a few kelvins
below the external temperature. Afterwards, quasi-steady state heat transfer was
observed until the doors were closed.

Lafaye de Micheaux et al. [60] and Croquer et al. [4] have calculated infiltration rate
as: ZZ
I˙ = vn dAdoor,inf iltration (2.25)

where vn is the velocity component normal to the doorway, ρda the density of dry air, hda
the enthalpy of dry air, ρva the density of water vapour, and hva the enthalpy of water
vapour.
At t = 10 s after the door opening (Figure 2.10), flat temperature profiles are observed
for the upper part (y/H ≈ 0.4 to y/H = 1) and the lower part (y/H = 0 to y/H ≈ 0.4) of
the doorway [4; 60]. At approximately y = H3 , transition between the profiles is sudden,
representing the "neutral level" observed for the refrigerated chamber [50]. Also for the
velocity profile, the two peaks of velocity magnitudes are observed at the ceiling and the
floor of the trailer with minimum obtained at the neutral level [60; 4].

33
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Figure 2.10: The temperature and velocity profiles at the truck trailer door vertical centerline
[60].

At t = 20 s, temperature profiles remain the same above the neutral level. Below the
neutral level, temperature has increased but maintains its flat profile while the region
close to the floor keeps its initial temperature. This temperature change occurs as the
outside air mixes with the inside air that stays in the truck. For the velocity profile, upper
region shows a flat profile while the lower region retain its profile. The magnitude of the
velocity decreased throughout the doorway and the neutral point remains at the same
height. The change in the profile of the upper region may be due to the outside air being
sucked into the truck trailer after the renewal of the inside air [60]
At t = 40 s, temperature profile no longer exhibits distinguishable neutral level and
an increasing profile along the door opening from the floor to the ceiling can be observed
(Figure 2.9). Generally decreased velocity magnitudes and increased temperatures are
observed along the doorway. Finally, at t = 60 s , the temperature and velocity profiles
are almost identical to the ones at 40 s, indicating quasi-state has been reached.
Importantly, both Lafaye de Micheaux et al. [60] and Croquer et al. [4] numerically
demonstrated that the 3D features of the flow does not significantly influence the tempera-

34
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Figure 2.11: Total and sensible infiltration heat loads compared to energy received by the
internal air volume [60].

ture and velocity distributions at the midplane of the doorway as presented in Figure 2.10
and Figure 2.9. The greatest difference for the 3D model is attained at the t = 11 s to 22 s
period, coinciding with the passing of a large eddy across the door. This is not properly
observed from the experimental results of Lafaye de Micheaux et al. [60] as the values are
averaged over the door centerline, while the 3D results comprise a full integration over
the door surface. However, while the 3D model shows in general a better agreement with
the experimental data, it does not bring a significant improvement over the results of the
2D models [4].
During the door opening period, the sensible heat load related to the mass interchange
is mainly responsible for the temperature increase of the internal air volume [60]. The
sensible infiltration load (Qsens
˙ ) and the total heat load (Qtotal
˙ ) are obtained by:
ZZ
˙ =
Qsens vn ρda hda dAdoor (2.26)

vn + |vn |
ZZ
Q̇total = vn ρda hda + ρva hva dAdoor (2.27)
2
where vn is the velocity normal to the doorway, ρda the density of dry air, hda the enthalpy
of dry air, ρva the density of water vapour, and hva the enthalpy of water vapour.
From Figure 2.11, the discrepancies between Q̇total and Qsens
˙ were deemed insignificant
as the uncertainty ranges overlap [60]. Therefore, latent heat load analysis and predictions
were not considered by Lafaye de Micheaux et al. [60].
In recent years, various researchers experimentally and numerically investigated the
heat and mass transfers inside a refrigerated truck trailer. Table 2.4 presents the summary
of the recent studies investigating the performance of refrigerated vehicles, focusing on
the experimental and CFD approaches.
For the experimental approaches, heat and mass transfers inside a truck trailer with
deployment of air curtains and varying aperture door openings were investigated by the
researchers [61; 62; 60]. Tso et al. [61] experimentally investigated the influence of air,
fan and plastic strip air curtains for a 7.2 m3 refrigerated truck body with an opening of
0.9 m2 , subjected to the typical conditions encountered in Singapore. The deployment of
an air curtain saved up to 40 % and 11 % of energy compared to cases without an air

35
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Author(s) Cargo Turbulence Door Mesh Wall
Model & Volume function
Tso et al. [61] Experimental * V = 7.2 m3
Clavier et al. [62] Experimental *
Moureh et al. [63] * Experimental V = 81.8 m3 High
& RSM 749, 632 cells y+,
Log-law
Moureh et al. [64] * RSM V = 81.8 m3 High
749, 632 cells y+,
Log-law
Lafaye de Experimental, * Unstructured High
Micheaux et Realizable polyhedral y+,
al. [60] k−ϵ V = 32.4 m3 Log-law
21, 000 cells
Rai et al. [65; 66] * SST k − ϵ * Structured Standard
hexagonal
V = 44.2 m3
6.1 million
cells
4.6 million
cells
Getahun et * SST k − ω V = 86.3 m3 Low y+,
al. [67] 16 million Low-Re
cells
Jara et al. [68] SST k − ω Triangular Low y+,
V = 8.2 m3 Low-Re
222, 463 cells
Croquer et al. [4] SST k − ω * Unstructured Low y+,
triangular Low-Re
A = 19.2 m2
903, 084 cells
Bonaventure * SST k − ω * Unstructured Low y+,
et al. [51] triangular Low-Re
A = 19.2 m2
903, 084 cells

Table 2.4: Summary of recent studies focusing on the performance of refrigerated vehicles.
Symbol: * Cargo present/Door openings.

36
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 curtain and with a plastic strip curtain, respectively. Clavier et al. [62] experimentally
investigated the infiltration through the doorway for a medium-size refrigerated truck
trailer placed in a wind tunnel and compared the heat load with and without two types
of air curtain. They observed an average infiltration heat load of approximately 5.7 kW
for an unprotected door with an opening time of 18 min. Lafaye de Micheaux et al. [60]
numerically and experimentally investigated the velocity and temperature distributions
in a truck trailer with an inner volume of 32 m3 for full, 2/3 lateral, 1/3 lateral and 1/3
central aperture door openings. Their findings are presented in Section 2.4.2.
For the numerical approaches, many researchers investigated the heat and mass trans-
fer inside a refrigerated truck trailer with CFD to overcome many restrictions in experi-
ments, such as manpower, material and time requirements [63; 64; 65; 66; 60; 67; 68; 4; 51].
Moureh et al. [63; 64], Getahun et al. [67] and Jara et al. [68] numerically investigated
the airflow patterns and temperature distributions inside a refrigerated vehicle with and
without cargo for the closed configuration. Moureh et al. [63; 64] developed the model
by using FLUENT and second-moment closure with the Reynolds Stress Model (RSM)
for a truck trailer with an inner volume of 81.795 m3 containing 32 polystyrene slotted
boxes (1.728 m3 ) treated as a porous medium. By deploying air ducts on airflow above the
load, better diffusion of injected air along the container was observed. As the recirculating
structure located in the inlet area was split into three distinct structures with less intensity
by the air ducts, better aerodynamic deployment for each inlet jet, thus more homogeneity
of ventilation within the entire trailer was observed. Getahun et al. [67] modeled a reefer
with an inner volume of 86.3 m3 containing 20 pallets each holding 768 kg of apple fruit
by using FLUENT and calculating the turbulence by using k − ω Shear Stress Transport
(SST) model. They observed monotonous decrease in the airflow from the inlet side to
the door side, dropping by 50 % between the first pallet to the middle of the trailer and
by 42 % between the middle of the trailer and the last pallet located near the doorway.
Jara et al. [68] investigated the internal temperature distribution and the airflow patterns
without any load with an inner volume of 8.16 m3 by using CFX with k − ω SST model
to solve for the turbulence. The maximum temperature deviation of 0.43 ◦ C at the air
outlet of the system was observed compared to the experimental results.
For the open configuration, Lafaye de Micheaux et al. [60], Rai et al. [65; 66], Croquer
et al. [4] and Bonaventure et al. [51] investigated the heat and mass transfer across
the doorway. Lafaye de Micheaux et al. [60] numerically investigated the velocity and
temperature distributions in a truck trailer with an inner volume of 32 m3 with full,
2/3 lateral, 1/3 lateral and 1/3 central aperture door opening areas. They used STAR
CCM+ with realizable k − ϵ SST turbulence model. Rai et al. [65; 66] investigated the
performance of an air curtain for a 44.2 m3 refrigerated truck trailer using FLUENT with
a 3D standard k − ϵ SST model. They considered door opening period of 15 min and
trailer configurations with and without cargo. Similar to Lafaye de Micheaux et al. [60],
natural infiltration was found to be the main cause of cold air flowing out from the lower
part of the door. An air curtain at an optimum averaged velocity of 3.1 m s−1 placed
inside the trailer resulted in the reduction of energy consumption by 48 %. Croquer et
al. [4] developed a preliminary CFD model using CFX with k − ω SST model for a truck
trailer equipped with three eutectic plates placed in parallel at the back of the trailer.
Their k − ω turbulence model demonstrated slight improvement to the predictions of
the realizable k − ϵ model of Lafaye de Micheaux et al. [60]. They also optimized the
interplate distance of the eutectic plates and the geometry of the separating wall. For the
loading region of the trailer, an interplate distance of 0.06 m resulted in the lowest average

37
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 temperature over the cargo loading region compared to cases with an interplate distance
of 0.036 m and 0.10 m. Their work was extended later by Bonaventure et al. [51] for the
cases with the plates placed in series at the roof of the trailer and for a configuration with
cargo.

2.5 Conclusion
The literature review encompasses the latest knowledge, methods, and techniques related
to the various aspects of the research. It explores the potential of PCM systems to
replace conventional refrigeration units on a refrigerated truck trailer that rely on fossil
fuels. Moreover, it discusses the physical phenomena associated with sub-zero eutectic
systems, such as the solidification and melting of PCM and the formation of frost on
the eutectic system. Advanced numerical models in this field are also presented. Finally,
considering the manuscript’s focus on the application of eutectic systems in refrigerated
truck trailers, the thermoaeraulic behavior of the air during the door opening period is
examined and discussed.

38
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Chapter 3

Infiltration Model

To investigate the feasibility of the eutectic system for a refrigerated food transport sys-
tem, it is foremost important to accurately predict the thermoaeraulic behavior of the
air during the infiltration (door-opening) period. This chapter presents the obtained
results submitted in Applied Thermal Engineering and presented during the 18th In-
ternational Refrigeration and Air Conditioning Conference held at Purdue University,
Lafayette, United States.

3.1 Numerical modeling

In the following section, the geometrical configuration is first described in detail. The
numerical methodology (algorithms, numerical schemes) as well as the numerical param-
eters, including the mesh grid and the boundary conditions are then presented.

3.1.1 Geometrical modeling

The computational domain is a truck trailer of length 7.13 m and height 2.62 m. The
refrigeration system includes three eutectic plates of thickness 0.054 m and height/length
1.76 m and fans, which maintain the air circulation through the plates and the container
when the doors are closed. Figure 3.1a displays the two configurations assessed here:
plates placed vertically in parallel at the back of the trailer (Configuration A) and plates
in series located at the roof (Configuration B). For each configuration, two air circulation
modes are employed, namely the blowing and suction modes. For the blowing mode, air
is blown out from the plate region to the trailer and for the suction mode, air is sucked
into the plate region. The truck trailer is placed in a wind tunnel of length 31 m and
height 5.43 m maintained at 293 K with a zero velocity field (Fig.3.1b). The wind tunnel
is large enough that its impact on the infiltration behavior is negligible. However, some
interactions between the doorway of the truck trailer and the floor of the wind tunnel are
expected.
The plates are maintained at 244 K, which corresponds to the melting point of a typical
sodium-chloride mixture as phase change material (PCM). The plates have been designed
such that the system can operate during a whole working day of 8 hours. During this
time, the PCM is not fully melted. The assumption of a constant temperature for the
plates is assumed to be valid as it is a mixture of solid and liquid in thermal equilibrium.
The wall that separates the eutectic plates from the loading area of the trailer is 1.81 m
in length and 0.01 m in thickness for Configuration A and 5.45 m in length and 0.01 m in

39
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 (a)

(b)

Figure 3.1: (a) Trailer diagram with the position of the eutectic plates for Configuration A
(top) and Configuration B (bottom); (b) General schematics of the computational domain
with the relevant boundary conditions.

thickness for Configuration B. For both configurations, a distance of 0.06 m is maintained

between the wall and the plates based on the former results of Croquer et al. [4].
Figure 3.2 shows the details of the optimized plate system design obtained by Croquer
et al. [4]. The dimensions, which are based on commercially available plates and truck
container models, are summarized in Table 3.1. The inter plate spacing is set to 0.06 m
as for the distance between the plates and the walls.

3.1.2 Numerical method

The conservation equations for mass, momentum and energy are solved using the finite-
volume method. Second-order upwind and centered schemes are preferred respectively
for the spatial discretization of the advection and diffusion terms, while a second-order

40
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Figure 3.2: Details and dimensions of the eutectic plates.

Feature Distance [m]

Container length (a) 7.13
Container height (b) 2.62
Plate thickness (c) 0.054
Plate height (d) 1.76
(e), (f ), (g), (h) 0.2, 0.345, 0.2, 0.206
(i), (j), (k) 0.06, 0.06, 0.06

Table 3.1: Container and eutectic plate dimensions.

Backward-Euler scheme with implicit time-stepping is used for the temporal discretization.
The pressure-velocity coupling is overcome using a Rhie-Chow fourth-order algorithm.
The flow field is modeled by solving the (un)steady-state Reynolds Averaged Navier-
Stokes equations (URANS) including the thermal energy equation, as the kinetic energy
effects are negligible. Buoyancy air-driven flows with a maximum temperature difference
of 29 K are expected, therefore air may be assumed to behave as an ideal gas with vari-
able density and constant transport properties and the Boussinesq model accounts for the
density change [69]. Turbulence effects are modeled using the Shear Stress Transport k-ω
model developed by Menter [70] with production terms representing buoyancy contribu-
tions to the turbulent field. The k-ω SST model applies the robust k-ϵ formulation in the
free turbulence regions and the more near-wall accurate k-ω model within the boundary
layers, depending on a wall-distance weighted function. All calculations are performed
using the commercial software CFX ANSYS.v19.

41
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 3.1.3 Numerical parameters
Unstructured meshes of triangular (resp. tetrahedral) elements are generated for the
spatial discretization of the 2D (resp. 3D) domain, with a minimum of 10 prismatic layers
close to all solid boundaries with a growth rate set to 1.1 as well as local refinement in
the door interface region. To perform 2D simulations on CFX ANSYS.v19, the 2D mesh
was extruded in the normal direction by two elements for a total thickness of 0.02 m.
Translational periodic boundary conditions are imposed on both faces of the extrusion.
A mesh-sensitivity analysis was carried out in the 2D configuration to determine the
minimum element size. A mesh grid composed of about ∼ 9.0 × 105 elements has proven
to provide grid-independent results while guaranteeing a maximum dimensionless wall
coordinate y + always lower than 1, a prerequisite for the use of a low-Reynolds number
approach. Details of these mesh grids are displayed in Figure 3.3. A symmetry condition
is imposed at the mid-plane of the 3D case to reduce the computational cost. For 3D
cases, the mesh grid is composed of ∼ 1.3 × 107 elements.

Figure 3.3: Overview of the mesh grids for Configuration A (top) and Configuration B (bot-
tom) with a focus on the fan and the plate regions on the right.

Steady-state calculations are performed on a 2D domain as shown in Figure 3.3. An

overall external heat transfer coefficient E = 0.4 W m−2 K−1 corresponding to a heavily
insulated equipment is imposed on the external walls [51], with an outside atmospheric
temperature of 293 K to model thermal resistance of the walls and the air outside trailer
computational domain. The separating wall is made of polypropylene with a thermal
conductivity of 1.4 × 10−4 W m−1 K−1 and a specific heat of 1800 J kg−1 K−1 . Regarding
the ventilation system, the need for accurately modeling the fans is avoided by inserting
a source term in the momentum equation over a surface representing the fan diameter of
0.2 m , such that the air passing through the fans is accelerated to an averaged horizontal
velocity of 8.4 m s−1 in a prescribed direction. Hoang et al. [71] modeled also the fans
through a source term in the momentum equation but based on a quadratic law giving
the pressure losses overcome by the fan as a function of the superficial velocity. The
steady-state calculations are initialized with a constant pressure of 101.325 kPa, an initial
temperature of 253 K and a zero-velocity field (V⃗0 = ⃗0 m s−1 ).
For the transient simulations during the door opening period, the truck trailer is
placed inside the wind tunnel displayed in Figure 3.1. All wind tunnel walls are assumed
to be adiabatic with an opening boundary condition in CFX imposed on the rightmost

42
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 surface, with Topening = 293 K and Popening = 101.325 kP a. The direction of the flow is
specified to be normal to the boundary condition and the magnitude of the velocity is
solved as a part of the solution. The trailer is located 10.0 m away from the leftmost wind
tunnel wall and 0.335 m above the floor of the wind tunnel. A wall equivalent thermal
resistance of 1.84 m2 K W−1 is used to model the heat flux through the trailer walls. The
initial conditions of the wind tunnel are fixed to be 101.325 kPa, 293 K and 0 ms−1 at
t = 0 s. For the transient simulations, results from the steady-state simulations with the
corresponding operating mode and configuration are used to initialize the flow field inside
the trailer. Instantaneous door opening is assumed at t = 0 s and the fans are turned off.
An adaptive time-step approach is adopted with 3 to 5 target inner-loop iterations and
initial and maximum time steps equal to of 10−7 and 10−2 s, respectively. The convergence
criteria for mass, momentum and energy RMS residual levelss RMS are fixed under 10−4 .
This results in a typical time step of about 5 × 10−3 s. Simulations are performed using
the HPC facilities of Calcul Québec using 20 to 40 CPU nodes with 12 AMD Opteron
6172 cores 32 GB RAM per node.

3.2 Validation of the numerical model

In order to determine the optimal numerical setup offering good accuracy at the lowest
computational cost, a preliminary benchmark is carried out involving two aspects: (i) 3D
domain vs. 2D domain and (ii) heat conducting vs. adiabatic container walls. The present
numerical model is validated against the experimental and numerical data of Lafaye de
Micheaux et al. [60], which provided temperature and velocity profiles of the air flow
inside an empty truck trailer without the eutectic plates. The configuration is shown in
Figure 3.1. Its inner dimensions are 2.35 × 2.5 × 5.51 m. It is placed midway in a tunnel
of 5.43 × 5.5 × 31 m at a distance of 0.35 m from the floor. In the baseline experiment,
the container and wind tunnel initial temperatures were respectively 253 K and 293 K. At
t = 0 s, the container door was fully opened, and velocity and temperature measurements
inside the container and at the door were logged for about 60 s and compared to the
results obtained by the model of Lafaye de Micheaux et al. [60].
Unstructured meshes of triangular ([Link]) elements are generated as de-
scribed in Section 3.1.3. It results in mesh grids composed ∼ 1.46 × 105 and ∼ 1.34 × 107
elements for the 2D and 3D configurations, respectively. The boundary conditions are
similar to the case with the eutectic plates (Fig.3.1b). The container walls are assumed
to be either adiabatic or allowing for heat conduction with a wall equivalent thermal
resistance of 1.84 m2 K W−1 , value based on typical refrigerated container insulation lay-
ers [4]. Thus one could compare the results against the regulatory requirements [2]. A
general-connection interface is set at the container door to connect the meshes across the
doorway and allow for the air to enter and exit through the interface and a symmetry
condition is imposed at the mid-plane of the 3D case to reduce the computational cost.
For the benchmark study, three numerical configurations are considered: (i) 3D domain
without heat transfer, (ii) 2D domain with heat conduction through the container walls
and (iii) 2D domain with adiabatic walls. A comparison of the door velocity profiles at
t = 10 s, 20 s, 40 s and 60 s after the door opening is shown in Figure 3.4. A good
general agreement is observed in terms of both the instantaneous shape profiles and their
temporal evolution. The 3D numerical setup seems to better fit the experimental data, in
particular at t = 10 s and t = 40 s, nonetheless both 2D approaches show good accuracy
and even a better capture of the velocity profile at t = 60 s.

43
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 1.0 1.0

0.8 0.8

y / hdoor [-]

y / hdoor [-]
0.6 0.6

0.4 0.4

0.2 0.2

0.0 0.0
0.0 0.5 1.0 1.5 0.0 0.5 1.0 1.5
-1
Velocity [m s ] Velocity [m s-1]

(a) (b)
1.0 1.0

0.8 0.8
y / hdoor [-]

y / hdoor [-]
0.6 0.6

0.4 0.4

0.2 0.2

0.0 0.0
0.0 0.5 1.0 1.5 0.0 0.5 1.0 1.5
-1 -1
Velocity [m s ] Velocity [m s ]

(c) (d)
Figure 3.4: Comparison of the velocity magnitude profiles at the container door vertical
centerline between the present simulations and the published data of Lafaye de Micheaux
et al. [60] at (a) t = 10 s, (b) 20 s, (c) 40 s and (d) 60 s. Legend: • Experiment - Lafaye de
Micheaux et al. [60], CFD - Lafaye de Micheaux et al. [60], CFD 3D, CFD 2D with
no heat transfer, CFD 2D with heat transfer.

Figure 3.5 compares the door temperature profiles for the same conditions described
above. At t = 10 s, the air flow is yet in an early stage not far from the initial step
profile, therefore no conclusive differences are seen. Nonetheless as the time passes, the
profiles evolve into a curved shape with a maximum temperature of about 289 K, which all
numerical approaches seem to over predict by about 1 K. At t = 20 s and t = 40 s, both 2D
setups fail to fully replicate the experimental profile at the lower section (y/Hdoor ∼ 0.2),
this is particularly noticeable in the adiabatic walls setup at t = 20 s. Nonetheless the
differences are of the same range as the dispersion of the experimental data. Deviations are
also observed in the near-wall regions, particularly for y/Hdoor → 1, where the numerical
models predict a constant wall temperature of 293 K whereas the experimental data reveal
an important temperature gradient. The temperature gradient at the bottom is well
captured, specially by the 2D models.
The numerical results showed a recirculation located at the door region which leads to
turbulent fluctuations. The instantaneous profiles of Figures 3.4 and 3.5 show differences
between the numerical and experimental results. This scattering might be attributed to
the natural turbulent fluctuations, thus implying that the three numerical models are
quite capable of capturing the main flow features across the door.
Further validation is shown in Figure 3.6 in terms of the infiltration rate. It is defined
as the volume flow rate into the trailer across the doorway as follows:

44
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 1.0 1.0

0.8 0.8

y / hdoor [-]

y / hdoor [-]
0.6 0.6

0.4 0.4

0.2 0.2

0.0 0.0
250 260 270 280 290 250 260 270 280 290
Temperature [K] Temperature [K]

(a) (b)
1.0 1.0

0.8 0.8
y / hdoor [-]

y / hdoor [-]
0.6 0.6

0.4 0.4

0.2 0.2

0.0 0.0
250 260 270 280 290 250 260 270 280 290
Temperature [K] Temperature [K]

(c) (d)
Figure 3.5: Comparison of the temperature profiles at the container door vertical centerline
between the present simulations and the published data of Lafaye de Micheaux et al. [60]
at (a) t = 10 s, (b) 20 s, (c) 40 s and (d) 60 s. Legend: • Experiment 1 - Lafaye de Micheaux
et al. [60], Experiment 2 - Lafaye de Micheaux et al. [60], CFD - Lafaye de Micheaux et
al. [60], CFD 3D, CFD 2D with no heat transfer, CFD 2D with heat transfer.

ZZ
I˙ = vn dAdoor,inf iltration (3.1)

with vn the velocity normal to the doorway with the positive direction towards the trailer,
and Adoor,inf iltration the area of the doorway occupied by the air infiltrating into the trailer
of the doorway. .
In Figure 3.6, the three numerical setups are able to properly capture the main features
of the experimental results, in particular the initial flow acceleration and the time needed
for reaching the steady state. The greatest difference is attained by the 3D model in the
t = 11 s to 22 s period, which coincides with the passing of a large eddy across the door,
implying recirculation and thus negative infiltration values for a short period. It is unclear
if the experimental procedure properly reports these structures as the reported values are
averaged over the door centerline, whereas the 3D results comprise a full integration over
the door surface. This would also explain the better fit of the 2D models, which essentially
represent the flow across the door vertical centerline.
In conclusion, although the 3D approach shows in general a better agreement with
the experimental data, both in terms of velocity and temperature profiles across the
door, it does not bring a significant improvement over the results obtained with the 2D
models. Inclusion of the heat conduction across the container walls introduces additional
complications to the model. However, it has non-negligible effects during an entire delivery

45
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 3.0

Infiltration Rate [m3s-1]

2.5
2.0
1.5
1.0
0.5
0.0
0 10 20 30 40 50 60
Time [s]

Figure 3.6: Comparison in terms of the infiltration rate at the container door between the
present simulations and the experimental and numerical data of Lafaye de Micheaux et al.
[60]. Legend: • Experiment 1 - Lafaye de Micheaux et al. [60], CFD - Lafaye de Micheaux
et al. [60], CFD 3D, CFD 2D with no heat transfer, CFD 2D with heat transfer.

cycle and is specified by regulatory bodies [2]. Therefore the 2D model with the heat
conduction across the container walls was chosen for modeling the refrigerated container
with the eutectic plates.

3.3 Results and discussion

Results are analyzed in this section for the truck container with the eutectic plates de-
scribed in Section 4.2. Steady-state calculations are initially performed for the two config-
urations with the plates located at the back (Configuration A) and the roof of the trailer
(Configuration B) with two operating modes (Blowing and Suction). The aforementioned
configurations are compared without and with a given cargo load. Afterwards, results
obtained from the steady-state simulations are imposed as the initial values inside the
trailer before the door opening for the transient calculations to investigate the interaction
between the developed fluid flow inside the trailer and the outside-infiltrating atmospheric
air.

3.3.1 Configurations without cargo

Configurations A and B - suction and blowing modes are first compared for an empty
trailer.

Closed trailer without cargo

From the results obtained from the steady-state simulations of Configurations A and B
without the cargo, it appears that the temperature distribution inside the trailer is homo-
geneous with a near-zero velocity zone near the center of the trailer. This is observed in
Figure 3.7 for both configurations and fan directions. Slight inhomogeneity of the temper-
ature distribution is observed due to the recirculation zones at each corner of the trailer
with the maximum temperature inside the trailer. Table 3.2 reveals that the difference in
the plate location and the operating mode has negligible impact on the temperature inside
the trailer for the steady-state simulations, as the maximum difference is 0.5 K for the
area-averaged temperature over the cargo loading area (excluding the area occupied by

46
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 (a) (b)
Figure 3.7: (a) Temperature [K] and (b) velocity magnitude [ms−1 ] contours for Configura-
tions A and B. The blowing mode is always on the top and the suction mode at the bottom.

Maximum Area-averaged Area-averaged

Configuration temperature temperature velocity
[K] [K] [m s−1 ]
A - blowing 254.5 245.1 3.64
A - suction 252.1 245.6 4.21
B - blowing 253.5 245.4 4.30
B - suction 254.0 245.2 3.33

Table 3.2: Temperature and velocity magnitude obtained from the steady-state simulations.

the eutectic plates). The counterclockwise flow inside the trailer (Configuration A - suc-
tion and Configuration B - blowing) results in an overall lower area-averaged temperature
and in a higher velocity inside the trailer. For all configurations and operating modes, the
maximum temperatures inside the trailer are below the regulated cargo temperature of
frozen foodstuffs at 261.15 K and quick (deep) frozen foodstuffs at 255.15 K [2]. Also, the
homogeneity of the temperature field inside the cargo loading area respects the regulated
maximum temperature difference of 2 K [2].
Also, there are few distinctions that can be observed between each configuration:
1. For Configuration A with the blowing mode (Fig. 3.7), there is no flow between the
first plate and the separating wall, effectively reducing the influence of the plate.
For Configuration B, all the plates are fully utilized.

2. For Configuration B with the blowing mode, the flow velocity from the plates to

47
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 the trailer is reinforced by gravity. On the contrary for the suction mode, the
flow being dragged from the bottom of the trailer to the plates is weakened by
gravity. Therefore, the highest area-averaged velocity is obtained for Configuration
B - blowing mode and the lowest is obtained for Configuration B - suction mode. For
Configuration A, gravity similarly reinforces and weakens the strength of the flow
for the suction and blowing modes, respectively. However, the influence is weaker
as the fans are located closer to the bottom of the trailer.

Door opening period without cargo

Config A - blowing Config A - blowing

Config A - suction Config A - suction
Config A - suction - Cargo Config A - suction - Cargo
Config B - blowing Config B - blowing
Config B - suction Config B - suction
Infiltration rate [m3 s -1 ]

Config B - blowing - Cargo Config B - blowing - Cargo

Heat load [kW]

Time[s] Time[s]

(a) (b)
Figure 3.8: Temporal evolutions of the (a) infiltration rate and (b) heat load through the door
for the two configurations and modes.

The infiltration rate (I)

˙ is defined as the volume flow rate of the air infiltrating into
the trailer and the heat load (Q̇) is defined as the amount of heat exchange occurring
through the entire doorway. The calculation of the infiltration rate is defined in Eq. 3.1.
The heat load across the doorway is calculated as:
ZZ
P
Q̇ = vn cp,da (T − 273) dAdoor (3.2)
RT
where vn is the velocity normal to the doorway with the positive direction towards the
trailer, P the pressure, R = 287.085 J kg−1 K−1 the specific gas constant of dry air and
cp,da = 1006 J kg−1 K−1 its specific heat. In Eq. 3.2, the reference temperature is chosen
as 273 K. Note that the infiltration rate is defined as the volume flow rate into the trailer.
However, for the heat load induced through the door, both exfiltration and infiltration
through the door are accounted. To obtain the infiltration rate and the heat load across
the entire doorway, the 2D values were scaled to that of 2.5 m in width as in Lafaye de
Micheaux et al. [60] and Croquer et al. [4].
The renewal time is defined as the time required by the internal air to be completely
renewed [60]. However, it should be noted that the simulations are initialized by a pre-
existing flow from the steady-state results before the door opening. Therefore, the infil-
tration rate does not directly correspond to the amount of pure atmospheric air entering
into the trailer. Especially in the first few seconds (t = 0 to ≈ 5.0 s) mixing of cold high-
velocity air and the hot atmospheric air is prominent near the doorway and can be the

48
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 main source of infiltration. For Configuration B (Fig.3.8b) with the blowing and suction
modes at t = 5.3 s and t = 2.6 s respectively, the heat load induced through the door is
negative as the exfiltration of high-temperature air occurs due to this interaction. How-
ever, it appears that for the majority of the simulation, the heat load induced through
the door directly corresponds to the infiltration rate. Therefore, the infiltration rate and
the renewal times from Figure 3.8a and Table 3.3 are used as the parameters to determine
the performance of the configurations and the operating modes.
Normalized infiltration

Config A - blowing
Config A - suction
Config A - suction - Cargo
Config B - blowing
Config B - suction
Config B - blowing - Cargo
CFD - Foster et al.
CFD - extended - Foster et al.

Time[s]

Figure 3.9: Time evolution of the normalized infiltration compared against the CFD results
obtained from Foster et al. (Extended: extended boundary CFD model) [49].

Foster et al. [49] split the infiltration profile into three regions: lag, steady-state, and
tail-off regions. The time evolution of the normalized infiltration, normalized by their
given values at t = 40 s, exhibits comparable profiles (Fig. 3.9). Extended boundary
CFD model of Foster et al. [49] refers to the model with the extended domain from 3 m
to 6 m outside the walls of the modelled cold room [49]. For the lag region, from t = 0
to t = 5 s, the current model demonstrates higher normalized infiltration gradients at
the door vertical centerline as the trailer is initialized with the steady-state results with
the pre-existing velocity. In the model of Foster et al. [49], the fans were turned off 30 s
before the door opening. Otherwise, the current model similarly demonstrates a steady-
state region from t = 5 to t = 30 s with a constant normalized infiltration gradient as the
infiltration flow is fully developed. The tail-off region from t = 30 to t = 40 s demonstrates
a decreasing normalized infiltration gradient due to the reduced temperature difference
between the air inside the trailer and the outside.
Figure 3.10 shows the temperature at 5 instants and the velocity at t = 10 s. At
t = 2.5 s, the pre-existing flow inside the trailer dominates the flow. Depending on the
direction of the rotation of the pre-existing flow inside the trailer, some atmospheric air
is dragged into the trailer through the top or bottom regions of the trailer. As infiltration
occurs at a relatively high velocity (Fig. 3.8), peak infiltration occurs within this time
interval. However, further infiltration and the penetration of the atmospheric air is limited
by the exiting flow acting similar to an air curtain investigated by Foster el al. [72], Hayes
and Stoecker [73], or Tso et al. [61].

49
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Configuration Renewal Peak infiltration Peak infiltration
time rate time
[s] [m3 s−1 ] [s]
A - blowing 14.5 6.7 3.3
A - suction 13.0 10.5 0.7
B - blowing 27.7 3.1 0.3
B - suction 24.5 7.9 0.2

Table 3.3: Infiltration data from the transient simulations without cargo.

Figure 3.10: (top) Time evolution of the temperature contours [K] for both configurations
and fan modes at five instants; (bottom) Corresponding velocity magnitude contours [ms−1 ]
at t = 10 s.

At t = 5 s, atmospheric air begins to penetrate through the top region of the doorway
into the trailer. The flow regime inside the trailer is still largely dominated by the pre-
existing flow and its interaction with the incoming atmospheric air. Mixing of hot and
cold air is still very prominent near the doorway due to this interaction. For Configuration
B with blowing mode, infiltration of the atmospheric air through the top region of the

50
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 doorway cannot be observed until t = 10 s due to the exfiltration of the cold high-velocity
air through the top of the doorway. From t = 10 s, further advancement of the hot
atmospheric air into the trailer can be observed. However, due to the presence of the
recirculation zones created by the pre-existing flow and the atmospheric air dragged into
the trailer before t = 2.5 s, further advancement of the atmospheric air into the trailer is
limited. This effect can be observed until t = 15 s and visualized by the velocity contours
in Figure 3.10b.
At t = 20 s, infiltration of the atmospheric air has reached the back of the trailer
as the recirculation zones have diminished. During this time, the flow is mostly driven
by buoyancy as the influence of the pre-existing flow progressively diminishes over time.
The time evolution of the infiltration rate from Figure 3.8 supports this statement as
it is relatively constant until t = 35 s. This infiltration dynamics was also observed by
Lafaye de Micheaux et al. [60] and named "buoyancy driven flow". Moreover, at t = 35 s,
the infiltration rate decreases slightly as the trailer is now completely filled with the
atmospheric air and the infiltration rate is caused solely by the temperature difference in
the plate region. This was also observed by Lafaye de Micheaux et al. [60] and named as
"bounday layer flow."
Table 3.3 and Figure 3.10 show that for either configuration, the operating mode has
a minor effect on the renewal time. The blowing mode only takes 11.5% and 13.1% more
time to reach the renewal time compared to the suction mode for Configurations A and
B, respectively.
However, significant differences are observed between Configurations A and B both
from the temperature contours and the time evolution of the infiltration rate. Table 3.10
shows that Configuration B takes 91.0% and 88.5% longer to reach the renewal time
compared to Configuration A with the blowing and suction modes, respectively. The
main cause of these differences can be seen from the temperature contours in Figure 3.10
from t = 0 s to t = 15 s. For both operating modes, for Configuration B, infiltration of
the atmospheric air is mainly due to the highest velocity of the cold air near the door,
where an air curtain behavior is observed. This limit infiltration and significantly increase
the renewal times for Configuration B.

3.3.2 Influence of the cargo

Similar simulations are performed for the closed and open situations with the presence of
cargo.

Closed configurations with cargo

To identify the influence of the cargo loads on the flow development and infiltration, 5
boxes full of frozen meats are placed in the trailer. The first cargo box is located 0.2 m
away from the separating wall and 0.15 m above the floor of the trailer to allow for the
flow to pass through between the cargo and the floor of the trailer. Proceeding boxes
are separated by 0.1 m from each other. The thermophysical properties of the cargo
consisting of frozen meat and cardboard packaging are extracted from Paquette et al. [74]
and displayed in Table 3.4. Each box is 1 m long and 1.2 m high. Conduction through the
boxes is accounted for in the simulations and the initial temperature of the cargo is set
to 258 K, as the maximum temperature of the frozen foodstuffs is regulated to be below
260.15 K [2].

51
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Material Density Specific heat Thermal conductivity
[kg m−3 ] [J kg−1 K−1 ] [W m−1 K−1 ]
Frozen meat 1050 3000 0.25
Cardboard packaging 276 1700 0.058

Table 3.4: Thermophysical properties of the cargo.

For the transient simulations with cargo, Configuration B - blowing mode was chosen
for its better performance in terms of the renewal time and Configuration A - suction
mode as a comparison to keep the direction of the pre-existing flow consistent.

(a) (b)
Figure 3.11: (a) Temperature [K] and (b) velocity magnitude [ms−1 ] contours of Configura-
tion A with cargo (top) and Configuration B with cargo (bottom).

Similar to Section 3.3.1, the temperature distribution inside the trailer remains rela-
tively homogeneous throughout the entire trailer (Fig. 3.11a). However, the main differ-
ence is that for Configuration B, due to the presence of the cargo, the main flow is blocked
by the narrow passage between the cargo and the bottom of the trailer. Due to this block-
age, a deflection of the flow and a strengthened recirculation zone at the back of the trailer
can be observed on Figure 3.11b. Also, the area-averaged velocity for Configuration B is
31.7 % lower than Configuration A. Consequently, the area-averaged temperature inside
the trailer for Configuration B is higher than Configuration A by 0.2 K. On the contrary,
for Configuration A, the flow from the fan is directly guided into the passage between the
floor of the trailer and the cargo boxes, resulting in neither flow blockage nor strengthened
recirculation zones. The homogeneity of the temperature distributions and the maximum
temperatures observed inside the trailer similarly meets the regulation criteria discussed
in Section 3.3.1.
As shown in Table 3.5, temperature inside the cargo is nearly homogeneous with a
difference of less than 0.1 K. Similarly, the maximum temperature observed inside the
cargo for Configuration A is 0.3 K lower than Configuration B.

Door opening period with cargo

The doors are opened at t = 0 s for the delivery of the frozen meat and Configurations A
and B are compared in terms of infiltration rate.

52
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Configuration Area- Area Maximum Minimum
averaged averaged cargo cargo
temperature velocity temperature temperature
[K] [m s−1 ] [K] [K]
A - suction 244.9 3.24 244.9 244.8
B - blowing 245.1 2.46 245.2 245.1

Table 3.5: Temperature and velocity magnitude data from the steady-state simulations with
cargo.

Figure 3.12: Temperature contours [K] for Configurations A - suction (top) and B - blowing
(bottom) at t = 2.5, 5, 10 and 20 s (from left to right).

Configuration Renewal Peak infiltration Peak infiltration

time rate time
[s] [m3 s−1 ] [s]
A - suction - no cargo 13.0 10.5 0.7
A- suction - cargo 25.5 10.7 0.6
B - blowing -no cargo 27.7 3.1 0.3
B - blowing - cargo 25.4 5.3 0.7

Table 3.6: Infiltration data for Configurations A and B with and without cargo.

The presence of the cargo results in an increase of the renewal time for Configuration
A by 96.2% (Table 3.6) compared to Configuration A without cargo. However, for Con-
figuration B, the renewal time decreases by 8.3% with an increase of 71.0% in the peak
infiltration rate. As shown in Figure 3.8, the time evolutions of the infiltration rate and
the heat load induced through the door exhibit similar profiles for Configurations A and
B with cargo. Peak infiltration occurs at t = 0.6 ∼ 0.7 s with a near constant infiltration
rate from t ∼ 7.5 s.
Due to the presence of the cargo, recirculation zones seen without the cargo in Fig-
ure 3.10b are no longer observed. As shown in Figure 3.12, During the first few seconds
after the door opening, the cold high-velocity air near the doorway exits the trailer. Then
the infiltration of the atmospheric air occurs along the top region of the trailer with the

53
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 exfiltration of the cold air above the cargo. The area of the infiltration progressively in-
creases before the trailer is completely filled with the atmospheric air. For Configuration
A, the renewal time increases due to the infiltration of the atmospheric air being limited
by the exiting cold air on top of the cargo. However, for Configuration B, the renewal time
decreases as the cargo prevents the formation of the recirculation zones. From the tem-
perature contours on Figure 3.12, the cargo limits the variation of flow structures that
can form inside the trailer and therefore both configurations demonstrate similar flow
pattern throughout the simulation. Due to minor temperature changes within the cargo
compared to the air inside the trailer, temperature contours of the cargo are excluded
from Figure 3.12.
However, there are minor differences that can be observed between Configurations A
and B with cargo. At t = 2.5 s, the exiting cold high-velocity air causes stronger mixing
to occur for Configuration A. This mixing limits the infiltration of pure atmospheric air
into the trailer until t ∼ 10 s. This is due to the velocity of the pre-existing flow being
higher than that of Configuration B as discussed in Section 3.3.2 (Table 3.5).
Temperature [K]

Config A - blowing
Config A - suction
Config A - suction - Cargo
Config B - blowing
Config B - suction
Config B - blowing - Cargo

Time[s]

Figure 3.13: Time evolution of the area averaged temperature inside the trailer.

From Figure 3.13, the temperature inside the trailer for Configuration B with cargo
is higher than Configuration A with cargo for the time interval t ∼ 15 s to t ∼ 25 s. As
the exiting cold high-velocity air mixing with the atmospheric air limits the infiltration
of pure atmospheric air for Configuration A. Regardless, the temperature change of the
air inside the trailer for all configurations and operating modes exhibits similar profile to
the temperature change observed by Rai et al. [65].

Configuration Cargo 1 Cargo 2 Cargo 3 Cargo 4 Cargo 5

A - suction 246.7 247.3 247.4 247.6 247.9
B - blowing 247.5 249.1 249.9 249.8 249.8

Table 3.7: Maximum temperature [K] of the cargo at t = 40 s.

54
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 As the temperature change is higher for Configuration B after the cold high-velocity
air exits the trailer, from Table 3.7, the maximum temperature observed at each cargo
box after 40 seconds is significantly higher for Configuration B (locations of each cargo
box described in Section 3.3.2). The temperature difference is greater than the initial
maximum temperature difference of 0.3 K (Table 3.5). A maximum temperature differ-
ence of 2.5 K for Cargo 3 and a minimum temperature difference of 0.8 K is observed
for the inner-most cargo (Cargo 1) between Configurations A and B, respectively. For
Configuration B, the highest cargo temperature is observed on Cargo 3 as the exiting cold
air on top of the cargo and the recirculation zones at the back of the trailer mixing with
the infiltrated atmospheric air limits other cargo boxes from making contact with the hot
air. For Configuration A, the highest cargo temperature is observed at the outer-most
cargo (Cargo 5), as the exiting cold air mixed with the atmospheric air makes the earliest
contact and increases its temperature. The maximum cargo temperatures observed are
still below the regulated temperature for quick (deep) frozen foodstuffs at 255.15 K [2].

3.4 Conclusion
To analyze the feasibility of the eutectic system, a numerical model has been developed
by employing the k − ω SST turbulence model to predict the thermoaeraulic behavior of
the air during the infiltration period for a refrigerated truck trailer equipped with eutectic
plates. The numerical model has been validated against the published experimental results
of Lafaye de Micheaux et al. [60]. During the validation, it was identified that while the
3D approach showed a better agreement with the experimental data, it does not bring
a significant improvement over the results obtained with the 2D models. Therefore, 2D
model with the heat conduction across the container walls was chosen for modelling the
refrigerated truck trailer with the eutectic plates.
With the numerical model, two different configurations were tested with the eutectic
plates placed in parallel at the back of the trailer (Configuration A) and in series lo-
cated along the roof of the trailer (Configuration B). Without the presence of the cargo,
Configuration B took 91.0 % and 88.5 % longer to reach the renewal time compared to
Configuration A with the blowing and suction modes respectively. For both fan oper-
ating modes, Configuration B exhibits an air curtain behavior and significantly limits
the infiltration of the hot atmospheric air. For both configurations, recirculation zones
created by the pre-existing flows within the truck trailer also limits the infiltration rate
and its advancement. However, when the cargo is introduced, recirculation zones are no
longer observed and both configurations demonstrate similar flow patterns throughout
the simulation with similar renewal times. During the simulation, maximum cargo tem-
peratures observed for either configurations are still below the regulated temperature for
quick (deep) frozen food stuff at 255.15 K [2]. Therefore, infiltration analysis validates
the viability of the eutectic system for the transport of frozen foodstuff. However, as the
eutectic plates were modelled as constant temperature plates without the phase change of
the eutectic mixture and the formation of the frost on top of the plates, further analysis
with the incorporation of the aforementioned phenomena is required.

55
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 56
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Chapter 4

Frost Model

For the transport of the frozen foodstuff, temperature is regulated to be below 255.15 K
for the duration of the transport [2]. To ensure the regulatory requirements are satisfied,
cold plates with the constant temperature of 244 K was employed for the infiltration model
in Chapter 3. Given the duration of the door opening period of less than two minutes
with the initial cold plate temperature of 244 K, desublimation process was assumed to be
the sole phase change phenomena occurring on top of the plate. Therefore, a numerical
model capable of accurately predicting the development of frost on top of a cold plate
was developed. This chapter presents the obtained results presented during the 7th IIR
International Conference on Sustainability and the Cold Chain held in United Kingdom.

4.1 Geometrical modeling

In the following section, the geometrical configuration is first described in detail. The nu-
merical methodology (algorithms, numerical schemes) as well as the numerical parameters,
including the mesh grid, the boundary conditions and the time step are then presented.
Figure 4.1 depicts the domain of the numerical simulation with the boundary con-
ditions developed for the frost formation on a cold flat plate. Since the validation was
performed against a number of published experimental results, readers are referred to the
published dimensions of each experimental setup. For each validation case, an adiabatic
wall of length L1 is added to the leading edge of the cold aluminum plate for the flow
to better resolve before interacting with the plate. For the validation against the results
of Nascimento et al. [75], a cold plate is also located on the top surface of the channel
at the same horizontal location as the bottom plate. Symmetry condition is imposed on
the top surface of the channel for the validation cases against the experimental results of
Cheng and Wu [76], Hermes et al. [77], and Wang et al. [78]. Heights of the experimental
setup for these authors are significantly greater than the frost thickness, ranging from
150 mm to 300 mm. Therefore, to reduce the computational cost, appropriate height with
symmetry condition is adopted for the simulation. Frost simulated at the opposite side
of the symmetry condition is not of importance, symmetry is solely utilized to simulate
the accurate air flow far away from the frost layer. As a matter of a fact, a channel with
50 mm in height with symmetric boundary condition applied on the top surface was tested
and showed negligible difference to the results obtained from the complete 150 mm height
channel.

57
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Figure 4.1: Numerical domain for the frost formation on a cold plate.

4.2 Numerical modeling

For the validation of the frost model, the air channel is initially assumed to be full of
humid air at rest with the inlet temperature and water vapour mass fraction, as shown
in Equations 4.1 and 4.2.

ua0 = 0, va0 = 0, Ta0 = 0, wva = win (4.1)

uc0 = 0, vc0 = 0, Tc0 = 0, wc0 = 0 (4.2)
where u, v are the horizontal and vertical velocitiy components of the humid air, T is
the temperature, w is the mass fraction, α the volume fraction, va the water vapour, a
the humid air, c the ice, and subscript 0 refers to the initial condition.
For the boundary conditions, all walls of the air channels are modelled as adiabatic
walls. Velocity inlet is applied to the inlet of the channel. Outflow condition is applied
at the outlet of the channel as shown in Equation 4.3 and Equation 4.4. Pressure outlet
boundary condition provided similar results as the outflow condition in terms of obtained
frost thickness and density. However, outflow boundary condition provided better numer-
ical stability of the solver.
∂ua ∂va ∂Ta ∂wva
= 0, = 0, = 0, =0 (4.3)
∂x ∂x ∂x ∂x
∂ua ∂va ∂Ta ∂αc
= 0, = 0, = 0, =0 (4.4)
∂x ∂x ∂x ∂x

4.2.1 Numerical method

The Eulerian-Eulerian granular multiphase model in ANSYS FLUENT 2020R2 is adopted.
The ice is formed directly from the water vapour desublimation and all other phase change
processes such as condensation are neglected. The primary phase is set to be the humid
air, aluminum and the secondary phase is set to be the ice droplets with the diameter of
10−5 m.
All the physical properties of the humid air are taken as the weighted averages of
the dry air and water vapour properties, except for the thermal conductivity. As the
concentration of the water vapour is relatively low, the thermal conductivity of the humid
air is assumed to be the value of the dry air. For the ice and the aluminium, the physical
properties within the FLUENT database is employed.
The conservation equations for mass, momentum and energy are solved using the finite-
volume method. A second-order Backward-Euler scheme with implicit time-stepping
is used for the temporal discretization. The pressure-velocity coupling is overcome us-
ing SIMPLE phase-coupled scheme. PRESTO! scheme was used for the pressure while

58
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 QUICK scheme was used for the momentum, volume fraction, energy and species conser-
vation. All the under-relaxation factors were set to be 0.2.
Depending on the Reynolds number, the k − ω Shear Stress Transport (SST) model
is adopted with second-order upwind discretization schemes for the turbulence kinetic
energy and specific dissipation rate. Otherwise, the flow is considered to be laminar. The
convergence criteria for all variables were chosen to be a relative error of 10−5 Simulations
are performed using the HPC facilities of Compute Canada or Calcul Québec’s networks
using 4 CPU nodes with 32 2 × Intel E5-2683 v4 Broadwell @ 2.1GHz cores per node.
Due to the nature of the Eulerian-Eulerian multiphase model, it is difficult to define
the boundary of the frost by an arbitrary limit of the ice volume fraction for each cell.
Therefore, thenumerical values of the frost thickness are determined by the MATLAB
Image Processing Tool. The Sobel gradient operator function (‘Imgradient’) is applied
to the normalized greyscale image (‘rgb2gray’) of the frost to define the outline of the
frost layer. The frost thickness is calculated as an average over the plate length. For
the validation against the results of Cheng and Wu [76], the frost thickness is measured
30 mm from the leading edge of the plate, as specified.

4.2.2 Governing equations

The governing equations include the mass, momentum, and energy conservation equations
with the species conservation equation for the water vapour. As only the water vapour
desublimation is being simulated, the mass transfer from the humid air to the ice governs
the frosting process. Based on the mass transfer, the momentum and energy transfers from
the water vapour to ice can occur and are reflected as the momentum and energy sources
of the governing equation. To account for the water diffusion driven by the concentration
gradient from the humid air to the frost layer, a source term is included in the species
conservation equation. Interested readers are referred to the works of Wu et al. [5; 46] for
the detailed description of the frost model.

Source terms

The source terms in the governing equations correspond to the mass transfer rate from
the humid air to the ice, ṁac , As only the desublimation is taken into consideration, the
mass transfer rate from the ice to the humid air, ṁca , is assumed to be zero. Therefore,
the source terms in the mass, momentum, energy and species conservation governing
equations are shown in Equations 4.5 to 4.8, respectively.

Smc = ṁac , Sma = −ṁac (4.5)

Suc = ṁac uc , Sua = −ṁac ua (4.6)

Shc = ṁac (hva + hl,va ), Sha = −ṁac hva (4.7)

Sva = −ṁac (4.8)

where h is enthalpy, hl the released latent heat per unit mass of water vapour during
the desublimation process. Subscripts m, u and h refers to mass, u to momentum, and h
to energy, respectively.

59
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Interphase Transfer
For the interphase momentum exchange, Wen-Yu’s model for the fluid-solid momentum
exchange coefficient is adopted [79; 46]. The momentum exchange coefficient K between
two phases K is defined as:
18ρa αa va αc f
K= (4.9)
d2c
where f is the drag function from the Wen-Yu’s experimental correlation:
f = (1 + 0.15Re0.687
c )αa−2.65 (4.10)

For the interphase heat exchange, the Ranz-Marshall correlation is preferred to identify
the Nusselt number N uca [80; 46]. The intensity of heat exchange is defined by the
temperature difference between the two phases:
⃗qca = −⃗qac = Eca (Tc − Ta ) (4.11)

where Fca is the heat transfer coefficients between the two phases defined as:
6λa αa αc N uca
Eca = (4.12)
d2c
where λa is the thermal conductivity of the humdid air and the Nusselt number N uca is
defined by the Ranz-Marshall correlation as:
N uca = 2.0 + 0.6Re1/2 1/3
c P ra (4.13)

Mass transfer rate

The mass transfer rate from the humid air to the ice ṁac , governs the frosting process
as the governing equations and the source terms are defined by it. Wu et al. [5; 46]
defined the mass transfer rate by a non-dimensional phase change driving force based on
Gibbs free energy as shown in Equation 2.23. Then Afrasiabian et al. [47] separated the
frosting condition in Equation 2.23 by two independent physical parameters as shown in
Equation 2.24. As the authors noted the modification of the frosting condition leading to
the results being less sensitive to the non-dimensional coefficient B, Equation 2.24 with
the modified frosting condition is employed.
While the numerical results obtained by Wu et al. [46] and Afrasiabian et al. [47]
demonstrated a good agreement with the experimental data, their models and the defi-
nition of the non-dimensional coefficients were limited to the inlet air velocities between
0.31 and 0.92 m s−1 , inlet air temperatures between −5 and 22 ◦ C, relative humidities be-
tween 38 and 85 % and cooling surface temperatures between −20 and −4 ◦ C. Therefore,
while adopting the mass transfer rate defined by Afrasiabian et al. [47] with the sepa-
rated frosting condition for the velocity and the water vapour concentration, equation for
the non-dimensional coefficient B was identified to minimize the difference between the
numerical results and the published experimental data [78; 77; 75; 76], as to be valid in
their extended range of operating conditions.
The modified non-dimensional coefficient for the frosting condition is defined as shown
in Equation 4.14.

B = (0.0000005328Tw4 − 0.0005446910Tw3 + 0.2085719116Tw2 − 35.4547217993Tw

(4.14)
+2257.4626692424) × (−0.0282087u2in + 0.868663uin + 0.4363793)

60
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 This definition of the modified non-dimensional coefficient does not apply for the vali-
dation cases against Nascimento et al. [75], as it involves a parallel cold plate configuration
and does not follow the configuration of the other validation cases. While the validation
against Nascimento et al. [75] demonstrates the applicability of the current model to
the parallel cold plate configuration, it would require further validations against simi-
lar configurations for the adjustment of the non-dimensional coefficient B. Furthermore,
definition has changed from the previous work of Jeong et al. [81] after larger sample
size of the non-dimensional coefficient values were tested against the published results via
trial-and-error.

4.3 Grid independence and time sensitivity analysis

For the air channel, rectangular meshes with a cell size of ∆x =0.1 mm and ∆y =0.5 mm
is used for all validation cases. Wu et al. [46] observed grid independence with a mesh of
same cell sizes with the inlet air velocities ranging from 0.31 to 0.92 m s−1 , with a time
step of ∆t =1.0 s. For the extreme conditions of Cheng and Wu [76] with the air velocity
of 10 m s−1 , grid independence analysis showed that a finer mesh of ∆x × ∆y = 0.05 mm
× 0.025 mm resulted in a negligible difference to the mesh used in the validation cases.
Through the time sensitivity analysis, time-step independent results were observed for
∆t = 0.05 s, 0.01 s, 0.005 s for the cases with uin = 1, 5 and 10 m s−1 , respectively.

4.4 Validation of the Frost model

The numerical model is validated for a wide range of operating conditions against the
published experimental results of Lee et al. [38], Cheng and Wu [76], Hermes et al. [77],
Wang et al. [78] and Nascimento et al. [75]. The operating conditions of each validation
case are summarized in Table 4.1.
Figure 4.2 illustrates an example of the numerical results obtained from the frost
model. Figure 4.2(a) and Figure 4.2(b) presents the obtained numerical volume fraction
of ice and temperature distributions for case No.12, thick condensed layer of frost with
the same temperature as the cold plate at 262.6 K has been formed near the plate with
the layers of the frost on top with increasing temperature. Figure 4.2(c) shows the humid
air velocity distributions at various times. As the velocity of the humid air near the frost
surface was small enough at the bottom of the velocity boundary layer of the humid air
flow, frosting condition of Equation 2.23 was satisfied for the mass transfer from the water
vapour in the humid air to the ice. Therefore, frost thickness increased over time as the
ice formed by the phase change deposited on the frost layer [46]. Also, the local increase
of the air velocity caused by the deflection of the air flow due to the frost can be observed,
demonstrating the capability of the current model to predict the interaction between the
humid air flow field changes and the frost growth.
Figure 4.3 displays the time evolution of the frost thickness for the validation cases
(as presented in Table 1) both obtained numerically and experimentally. It also presents
the experimental and numerical frost densities att = 3600 s. From t = 1800 s to t = 3600 s
(t = 900 s to t = 1800 s for cases No.15 and No.16), numerically obtained frost thicknesses
exhibit a maximum deviation of ±18.3 % against the experimental values. As the frost
thickness is relatively small for lower time intervals, percentage differences are not shown.
Numerically obtained frost densities are compared against the experimental results of

61
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Authors No. Cooling Surface Inlet Air Water Content Inlet Air
Temperature Temperature Level Velocity
Tw Tin w uin
 
kg
(◦ C) (◦ C) kgair
(m s−1 )
Hermes 1 −4 16 0.009 0.7
et al. [77] 2 −8 16 0.09 0.7
3 −12 16 0.09 0.7
4 −16 16 0.09 0.7
Nascimento 5 −10.3 11.8 0.006 0.9
et al. [75] 6 −24.5 11.7 0.0058 1
7 −24.8 6 0.0038 1
8 −10.3 5.7 0.0037 1
Lee 9 −15 25 0.0157 1
et al. [38] 10 −15 25 0.0098 1
Wang 11 −10.5 16.5 0.093 5
et al.[78] 12 −10.5 19 0.0109 5
13 −10.5 10.5 0.0063 5
14 −16 10.5 0.0063 5
Cheng 15 −8.2 28.6 0.0178 10
and Wu [76] 16 −8.2 28.6 0.0098 10

Table 4.1: Frosting conditions used to validate the frost model.

Hermes et al. [77], Wang et al. [78] and Nascimento et al. [75]. Experimental densities for
the other cases were not available. At t = 3600 s, the model predicts the frost densities
within 25.0 % of their experimental values except for case No.13 with 29.4 % difference.
A density difference of 29.4 % observed for case No.13 is suspected to be due to the
low temperature difference between the plate and the air at 22 K. Conventionally, frost
models are seen to have difficulty for the cases with low temperature difference between
the air and the plate [78]. For such a temperature difference, condensation could affect
the frosting process, which is not considered by the current model. However, such density
difference between the numerical results and the experimental data is not observed for
the low velocity cases with low temperature difference between the air and the plate. Due
to the high velocity of the air for case 13, the effect of condensation affecting the frosting
process might have been exacerbated. As the results of Wu et al. [46] reported absolute
deviation of 30
p % and 25 % for the averaged frost thickness and the weight, maximum
deviation of (30%)2 + (25%)2 ≈ 39% is approximated for the frost density. All the
obtained results fall within the deviations observed by Wu et al. [46]. The current model
with the modified frosting criterion coefficient provides improved results with an absolute
deviation of 18.3 % and 29.4 % for the frost thickness and density, respectively.

62
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 (a) (b)

(c)

Figure 4.2: Numerical results for Case No.12: (a) Volume fraction of ice at t = 3600 s, (b)
Temperature distribution (K ) at t = 3600 s, and (c) Air velocity distribution(ms−1 ).

Figure 4.4 displays the comparison between the experimental data and the numerical
results for both the frost thickness and the frost density. Generally, numerically obtained
results for both the frost thickness and the frost density exhibit a good agreement with
the experimental data. While the numerical results for case 13 provide inadequate frost
density with a 29.4 % difference compared to the experimental data, Figure 4.4 suggests
that the difference is still within an acceptable range.
By extending the models of Wu et al. [5; 46], and Afrasiabian et al. [47] a new nu-
merical frosting condition was developed for a wider range of operating conditions. The
modification was performed to the non-dimensional frosting criterion coefficient, B, and
k − ω SST turbulence model was applied. The present model was favorably validated in
terms of frost thickness and density against the experimental results of e et al. [38], Cheng
and Wu [76], Hermes et al. [77], Wang et al. [78] and Nascimento et al. [75]. However,
the model showed difficulty when the temperature difference between the air and the cold
plate was below 25 K for the tested operating conditions. The model can be further im-
proved by optimizing the equation of the non-dimensional frosting criterion coefficient, B.
However, significant efforts are required to repeat the simulations with varying coefficient
values. It would also require many validation cases with various operating conditions and
configurations to be applied as a generalized model for the frosting process. Furthermore,
identifying the exact range of Reynolds number for the application of the k − ω SST
turbulence model will be a difficult task.
While there are still improvements to be made to the frost model, it was deemed
that the validated frost model above is applicable to be used for the combined frost and

63
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Figure 4.3: Results for the frost thickness and density at t = 3600 s.

64
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Figure 4.3: Results for the frost thickness and density at t = 3600 s (cont.).

65
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 (a)

(b)

Figure 4.4: Comparison between experimental data and numerical results for the (a) Frost
thickness and (b) frost density.

solidification and melting model for the system. However, it should be noted that the
numerical errors from the frost model will be applied to the combined model as the frost
thickness and density both contribute towards the formation of heat transferred to the
surrounding air inside the air channel and to the PCM below the plate.

4.5 Conclusion
Eulerian-Eulerian granular multiphase model was developed based on the mass transfer
source term from the works of Wu et al. [5] and Afrasiabian et al.[47]. The numerical
model has been validated against the published numerical and experimental results. While
the numerical results obtained by Wu et al. [5] and Afrasiabian et al. [47] demonstrated
a good agreement with the experimental data, their models were limited to the specific

66
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 ranges of initial and boundary conditions. Especially, their models were limited to the inlet
air velocities below 1 m s−1 , far below the air velocities observed around the eutectic plates
from the infiltration model in Chapter 3. To extend the ranges of applicable initial and
boundary conditions, k − ω SST turbulence model introduced in the infiltration model in
Chapter 3 was applied along with the refinement of the non-dimensional frosting condition
coefficient B via trial-and-error against the published results. The validated ranges of the
operating parameters presented in this chapter are well within the close proximity to the
operating conditions of the refrigerated truck trailer equipped with eutectic plates with
the inlet air temperature of 21.6 ◦ C, inlet air velocity of 10 m s−1 , relative humidity of 69
%, and the initial plate temperature of −30 ◦ C.

67
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 68
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Chapter 5

Consolidated Model

The main objective of this chapter is to consolidate the infiltration and frost models
presented in Chapter 3 and Chapter 4, respectively, with the solidification and melting
model for the eutectic mixture inside the eutectic plates. As the Eulerian-Eulerian gran-
ular frost model presented in Chapter 4 incorporates the k − ω SST turbulence model
introduced in Chapter 3, it is capable of simulating the infiltration behavior of the air
and the frost formation along the eutectic plates for a truck trailer equipped with eutectic
plates. However, as the heat transferred to the eutectic plates from the frost formation
and the infiltrating hot atmospheric air is not uniform throughout the plate, melting of
the eutectic mixture is also assumed to be inhomogeneous. Therefore, solidification and
melting model to predict the phase change of the eutectic mixture inside the eutectic
plates must be incorporated. In this chapter, a consolidated numerical model will be
presented incorporating the infiltration model presented in Chapter 3, frost model pre-
sented in Chapter 4, and the solidification and the melting model. As the solidification
and melting model might behave unpredictably due to its interaction with the other nu-
merical models, the incorporated solidification and melting model is validated against the
published experimental and numerical results. Then the consolidated numerical model is
applied to the experimental test-bench and validated against the obtained experimental
results presented in this Chapter.

5.1 Numerical modeling

5.1.1 Geometrical modeling
The computational domain is an aluminum encasing containing the PCM with the length
of 16 in.(406.4 mm)and height of 2.0625 in.(52.3875 mm) with the inner dimension of 1.5
in. by 14 in.(38.1 mm by 355.6 mm). The inner volume is located 0.1875 in.(4.7625 mm)
away from the air channel with the 0.375 in.(9.525 mm) thickness at the top. The air
channel is located at the bottom of the aluminum encasing containing the PCM with the
length of 19.5 in.(495.3 mm) and width of 0.4 in.(10.16 mm). It should be noted that the
Figure shown in Figure 5.1 is not to scale to extenuate the PCM and the encasing.
The properties of the aluminum are taken from the ANSYS FLUENT database. The
properties of the PCM (Rubitherm SP-28, Rubitherm GmbH) are summarized in Ta-
ble 5.1.
For validation of the consolidated model, the air channel is initially assumed to be full
of humid air at rest with the inlet temperature and water vapour mass fraction, as shown

69
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 (a)

(b)

Figure 5.1: (a) Numerical/experimental schematic of the aluminum plate containing the
PCM and the air channel below; (b) Dimensions of the experimental/numerical domain.

Density, ρ 1250 kg m−3 (liquid/solid)

Heat capacity, Cp 2000 J kg−1 K−1
Thermal expansion coefficient, CT E 0.002 K−1
Solidus temperature, Ts −30 ◦ C
Liquidus temperature, Tl −28 ◦ C
Thermal conductivity, λ 0.6 W m−1 K−1
Dynamic viscosity, µ 0.12 kg m−1 s−1

Table 5.1: Properties of Rubiterm SP-28 (PCM).

in Equation 4.1 and Equation 4.2.

For the boundary conditions, all walls of the air channels are modelled as adiabatic
walls. Velocity inlet is applied to the inlet of the channel. Outflow condition is applied
at the outlet of the channel as shown in Equation 4.3 and Equation 4.4. Pressure outlet
boundary condition provided similar results as the outflow condition in terms of obtained
frost thickness and density. However, outflow boundary condition provided better numer-
ical stability of the solver.

70
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 5.1.2 Numerical method
The Eulerian-Eulerian granular multiphase model in ANSYS FLUENT 2020R2 from
Chapter 4 is adopted. For the phase change of the phase change material, as the Eulerian-
Eulerian multiphase model in ANSYS FLUENT 2020R2 does not allow for the existing
solidification and melting model to be enabled, it was applied via User-Defined Function
(UDF). The primary phase is set to be the humid air, aluminum, and the PCM and the
secondary phase is set to be the ice droplets with the diameter of 10−5 m. The solidifica-
tion and melting UDF solely functions within the primary phase of the model between
the solidus and liquidus temperatures of the PCM. It should be noted that this was only
possible due to the significant temperature difference between the humid air and the ice
compared to the liquidus and solidus temperature of the PCM.

5.1.3 Governing equations

The governing equations are identical to the frost model described in Section 4.2.2, the
source terms are listed in Equation 4.5, Equation 4.6, Equation 4.7, and Equation 4.8.
An additional source term is applied in the momentum equation to describe the flow of
the PCM due to the solidification and melting process as shown in Equation 5.1.

Ssm = −A(ϕ)v (5.1)

where A(ϕ) is the porosity function [12]. The source term Ssm is used to describe
the flow of the PCM in the porous medium in the momentum equation. It is zero in the
liquid phase to allow for the free motion while it is arbitrarily large value in the solid
phase to force the velocity values to be approximately zero. To achieve this requirement,
Carman-Kozeny equation derived from the Darcy law for fluid flow in porous media in a
modified form is used as shown in Equation 5.2 [11; 12].

Amush (1 − ϕ)2
A(ϕ) = (5.2)
ϕ3 + χ

0,
 if T < Ts
ϕ= T −Ts
, if Ts < T < Tl (5.3)
 Tl −Ts
1, if T > Tl


where Ts and Tl are solidus and liquidus temperatures, χ is a small computational

constant to prevent division by zero (χ = 0.001 for this study), Amush is the mushy zone
constant to measure how fast the fluid velocity approaches zero as it solidifies (Amush = 106
for this study [11]), and ϕ is the liquid fraction as shown in Equation 5.3.

5.1.4 Solidification and melting

For the primary phase of the numerical model, solidification and melting of the PCM
is applied via a UDF. A source term, Ssm , is introduced in the momentum equation as
shown in Equation 5.1, and the energy equation is defined as shown in Equation 5.4. The
specific enthalpy (H) is the sum of the sensible enthalpy (h) and the latent heat (∆H).
Z T
H = h + ∆h = (href + Cp dT ) + (ϕL) (5.4)
Tref

71
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 where href is the reference enthalpy at the reference temperature Tref , Cp is the specific
heat, and L the latent heat.

5.2 Grid independence and time sensitivity analysis

For the aluminium encasing and the PCM, rectangular meshes with a cell size of ∆x =
1 mm and ∆y = 1.5 mm is used for all validation cases [11]. Higher mesh density resulted
in a negligible difference in the melt front location and the average melt fraction. For the
time step, while the time independent results were obtained for ∆t = 0.1 s, time steps for
the air channel was applied as the overall system requires the smaller of the time steps to
observe time-step independent results.

5.3 Validation of the UDF solidification and melting model

Since ANSYS FLUENT v2020 R2 is incapable of applying the existing solidification and
melting model with the Eulerian-Eulerian multiphase model used for the frost model, the
solidification and melting model is applied via a UDF as described in Section 5.1.4.
For the validation of the UDF solidification and melting model, the obtained numerical
results are compared against the experimental and numerical results of Fadl and Eames
[11] and Kamkari and Shokouhmand [82]. For the details regarding the experimental
setup, readers are referred to the works of Fadl and Eames [11]. Conventionally within
the solidification and melting model present in the ANSYS FLUENT v2020 R2, the
thermophysical properties are based on the Boussinesq model for the natural convection
modeling with mean property values between liquidus and solidus temperature applied.
While the application of the mean property values is reasonable, it completely neglects
the possible unique profiles of the thermophysical properties between the liquidus and
solidus temperature of a given fluid and could lead to the significant differences between
the numerical and experimental results.
Figure 5.2 displays the predicted melt front locations for the ANSYS FLUENT version
of the solidification and melting model denoted as Model 1 in the legend of this figure with
the thermophysical properties based on the Boussinesq model for the natural convection
modeling. For the varying values of Amush ranging from 106 to 107 , the number of peaks
and troughs of the melt front varies. Locations of the peaks and troughs are inaccurately
predicted both from the ANSYS FLUENT model applied and the simulated results of
Fadl and Eams [11]. However, the value of Amush = 106 presents the most accurately
melt front location and hence the most accurate melt fraction of the PCM.
For the UDF version of the solidification and melting model, the Boussinesq model
for the natural convection modeling previous applied to the ANSYS FLUENT version of
the solidification and melting model cannot be applied as the ANSYS FLUENT does not
allow for the existing solidification and melting model to be enabled with the Eulerian-
Eulerian multiphase model. Therefore, a more simple approach is applied where the
density is linearly interpolated between the liquidus and solidus temperature and the
viscosity is assumed to be constant except a significantly large number when the PCM is
solid. Furthermore the specific enthalpy of the PCM varies according to section 5.1.4 to
incorporate the solidification and melting of the PCM as a UDF by varying specific heat
of the PCM.

72
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Figure 5.3 compares the melt front locations between the UDF and the ANSYS FLU-
ENT versions of the solidification and melting model. The UDF version of the solidifi-
cation and melting model is denoted as Model 2. For t = 600 s to t = 2400 s, negligible
difference is observed between the ANSYS FLUENT model and the UDF solidification
and melting model. While the locations of the peaks and troughs of the melt front are
different, average y location of the melt front and the number of peaks and troughs of the
melt front are nearly identical. At t = 3600 s, the accuracy of the UDF model slightly
decreases as it predicts the average y location of the melt front 13.1 % and 7.6 % higher
than that of the experimental value and the Boussinesq model, respectively. However, the
average y location of the melt front and the melt fraction is still within the acceptable
range with 4.2 % to 13.1 % difference compared to the experimental results, respectively.

(a)

(b)

Figure 5.2: Melt front of the Solidification and Melting model for Amush = 106 , 5 × 106 and
107 at (a) t = 600 s, (b) t = 1200 s, (c) t = 2400 s, (d) t = 3600 s.

73
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 (c)

(d)

Figure 5.2: Melt front of the Solidification and Melting model for Amush = 106 , 5 × 106 and
107 at (a) t = 600 s, (b) t = 1200 s, (c) t = 2400 s, (d) t = 3600 s (cont.).

It should be noted that greater differences arises for the average melt front y location at
t = 600 s for all simulated results compared to the experimental result.
The peaks and troughs of the melt front are formed due to the formation of the
Rayleigh-Bénard convection cells as shown in Figure 5.5. As the bottom of the experimen-
tal setup is exposed to the high temperature while the rest of the sides are adiabatic (cov-
ered via plexiglass and insulators), Rayleigh-Bénard convection cells are created. From
both the Boussinesq model and the UDF model, accurate prediction of the locations of
the peaks and troughs of the melt front could not be obtained. This indicates that for
an accurate prediction of the locations of the peaks and troughs of the melt front, more
accurate approach in modelling the thermophysical properties and evaluation of the ap-
propriate Amush value is required. However, as shown in Figure 5.4, the average y location
of the melt front and the melt fraction values are within the satisfactory range, the UDF

74
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 (a)

(b)

Figure 5.3: Melt front of the UDF Solidification and Melting model for Amush = 106 at (a) t =
600 s, (b) t = 1200 s, (c) t = 2400 s, (d) t = 3600 s.

solidification and melting model is applied for the complete model.

5.4 Consolidated Model Validation

5.4.1 Experimental setup
For the experimental validation of the consolidated numerical model, the experimental
setup of Rahal [83] was used. The 2D schematic of the experimental setup can bee seen
in Figure 5.1. For more detailed information regarding the experimental setup, readers
are referred to the works of Rahal [83].
The experimental investigation has been performed by a set up primarily composed

75
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 (c)

(d))

Figure 5.3: Melt front of the UDF Solidification and Melting model for Amush = 106 at (a) t =
600 s, (b) t = 1200 s, (c) t = 2400 s, (d) t = 3600 s. (cont.)

of two components:
1. An air channel with the converging air inlet and diverging outlet made out of plastic
resin by additive manufacturing.
2. A container for the PCM made out of aluminum that is connected to the test section
of the air channel. Aluminum panels of varying sizes covers the upper portion of
the PCM container as shown in Figure 5.6(d).
The aluminum encasing containing the PCM is connected to the air channel test section,
where the upper portion of the aluminum encasing acts as the floor of the test section of
the air channel. Then a plexiglass cover acts as the upper and the side walls of the test
section as show in Figure 5.6 for the observation and measurement of the frost.

76
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 (a)

(b)

Figure 5.4: (a) Melt fraction for different values of Amush using the UDF model; (b)Average
melt front y location for two applied solidification and melting schemes and Amush = 106 .

The aluminum PCM encasing is 355 mmin length, 304 mm in width and 34 mm in
depth. It is divided into four identical compartments as shown in Figure 5.6(b). As
displayed in Figure 5.6(c) for one of the divided compartments, total of 16 type T ther-
mocouples (Omega Engineering TT-T-36-25) separated into four groups of four thermo-
couples placed equidistant to each other are placed within the compartment to measure
the temperature of the PCM. The coil pipe system within the encasing is connected to the
Thermo Bath (Thermo Haake ARCTIC AC200 A40 Immersion Bath) with 900 W cooling
capacity at 20 ◦ C. For the thermofluid, silicone oil SIL180 with the thermal conductivity
of 0.12 W m−1 K−1 and dynamic viscosity of 0.041 kg m−1 s−1 from Thermofischer is used to
transfer the heat away from the PCM during the solidification phase. Two surface thermo-
couples (SA 1XL-T-SRTC, Omega Engineering) are also attached to the inlet and outlet

77
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Figure 5.5: Velocity magnitude of the melted PCM for Amush = 106 at four different times.

(a)

(b) (c) (d)

Figure 5.6: (a) CAD view of the experimental setup, and zooms on the (b) aluminum en-
casing, (c) thermocouples locations inside the aluminium encasing, (d) aluminum panels of
varying sizes for measuring the frost weight.

of the coil pipes. The walls of the aluminium encasing are covered with the polyurethane
foam (McMaster-Carr, ref. 9385K21) with the thermal conductivity of 0.022 W m−1 K−1
for insulation.
The inlet of the experimental setup is connected to the climatic chamber HSH27C at
the PIMUS laboratory of Université de Sherbrooke to control the relative humidity and

78
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Measurement Provider Model Number Range Uncertainty
Frost thickness HSI 1280 × 1080 pixels 0.1mm
PCM Temperature Omega TT-T-36-25 −40 ◦ C to 30 ◦ C ±0.2 ◦ C
Surface Temperature Omega SA 1XL-T-SRTC −29 ◦ C to 260 ◦ C ±0.05 ◦ C

Table 5.2: Uncertainty of the equipment.

the temperature of the air. For the velocity of the air, Vortex Powerfan S-1000 with the
duct diameter of 10 in. is used to control the velocity from 0.1 m s−1 to 0.5 m s−1 . A
RDT camera from HSI with the resolution of 1280 × 1080 pixels are used to measure the
thickness of the frost. The acquisition of the photos are carried out with MIDAS image
processing software. The uncertainties of the equipment used in the experiment are listed
in Table 5.2.
For the validation of the numerical model, eutectic mixture of potassium bicarbonate
solution with the phase change temperature of −5.4 ◦ C and the latent heat of fusion of
268.54 kJ kg−1 was used [84]. Due to the limitation of the experimental setup, PCM with
the lower phase change temperature could not be explored.

5.4.2 Experimental setup adjustment

As the air channel is located on the upper portion of the PCM container for the exper-
imental setup of Rahal [83], it resulted in no observable amount of frost formed during
the experiment. Results of the numerical simulation of the experimental setup with the
initial condition of Figure 5.1 (a) are shown in Figure 5.7 and Figure 5.8. From the melt
fraction of the PCM shown in Figure 5.7, melting of the PCM is shown to be focused
primarily in the upper portion of the PCM. As the high temperature air flow is above the
plate, heat is primarily transferred from the air through the plate to the upper portion
of the PCM. Therefore, highly localized melting of the PCM in the upper region of the
PCM can be observed. From t =600 s to t =1800 s, melting of the sides of the PCM can
also be observed due to the heat transferred through the aluminum encasing, but it is
difficult to conclude that the PCM is being fully utilized as the heat transfer is primarily
via conduction alone within the PCM.

Figure 5.7: 2D contours of the PCM melt fraction at four different times.

79
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 (a) (b)

(c)

Figure 5.8: Numerical results for the (a) maximum plate temperature at t = 1800 s, (b)
maximum PCM temperature over time and (c) velocity magnitude contours of the PCM at
t = 1800 s.

Furthermore, maximum plate temperature at t =1800 s shown in Figure 5.8 (a) indi-
cate a temperature difference of the plate being up to 22 K. This indicates the existence
of highly localized high temperature region within the plate (at x ≈ 0.23 m) caused by
the localized melting of the PCM. From Figure 5.8(b), after t = 800 s, maximum PCM
temperature of over 273.15 K is observed, effectively nullifying the conditions for the frost
to be formed on top of the plate. Additionally, velocity contour of the PCM at t = 1800 s
shown in Figure 5.8(c) demonstrate no movement of the melted PCM in the upper region
of the PCM, concluding that the heat transfer is primarily via conduction with the current
setup. There is slight buoyancy driven movement of the melted PCM on the sides caused
by the temperature-driven density difference. However, it is localized to the sides of the
PCM and does not influence the upper portion of the melted PCM, where the majority
of the heat transfer from the air to the PCM occurs.
From the experiment, it was also observed that as the PCM melts, there is a loss
of contact between the PCM and the aluminum plate on top of the PCM due to the
volumetric change. With the highly localized melting of the PCM in the upper region of
the PCM, contact between the PCM and the plate is immediately lost, effectively limiting
the influence of the PCM. Therefore, the experimental setup of Rahal [83] was flipped.
The aluminum plates of varying sizes originally covering the PCM was replaced with a
single aluminium panel with the same thickness. For the new experimental setup, the air
channel is located at the bottom of the PCM container. Due to the nature of the setup,
mass of the frost cannot be measured as it was originally planned with the aluminum
plates of varying sizes of Rahal [83]. Therefore the validation will primarily focus on the
thickness of the frost and the temperature of the PCM.

80
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 5.5 Experimental results
The initial temperature of the PCM containment system was set to be at T = −7 ◦ C to
ensure the complete solidification of the eutectic solution of potassium bicarbonate PCM
for the experiment. Then two different inlet conditions were tested with temperatures
at Tinlet = 5 ◦ C and 0 ◦ C, where the relative humidity and velocity of the inlet was fixed
to be 70 % and 0.1 m s−1 , respectively. Due to the adjusted experimental setup with the
air channel located below the eutectic plate, only temperature of the PCM and the frost
thickness were validated.
Table 5.3 and Table 5.4 displays the comparison between the numerical and exper-
imental results for the frost thickness and PCM temperature, respectively. Obtained
experimental frost thicknesses for both inlet temperatures of T = 5 ◦ C and 0 ◦ C agrees
with the numerical results accounting for the uncertainty of the camera being 0.1 mm.
However, as the maximum frost thickness obtained during the experiment is 0.15 mm
with 0.1 mm uncertainty, it is difficult to conclude on the accuracy of the experimental
results and its comparison against the numerical results.

Tin = 5 ◦ C Tin = 0 ◦ C
Time Experimental Numerical Experimental Numerical
[min] [mm] [mm] [mm] [mm]
0 0 0 0 0
5 0 0 0.08 0
10 0.06 0.04 0.08 0.07
15 0.06 0.07 0.13 0.09
20 0.15 0.08 0.13 0.09
25 0.12 0.09 0.14 0.09
30 0.12 0.09 0.13 0.09

Table 5.3: Experimental and numerical frost thickness for Tin = 5 ◦ C and Tin = 0 ◦ C.

Tin = 5 ◦ C Tin = 0 ◦ C
Time Experimental Numerical Experimental Numerical
[min] [◦ C] [◦ C [◦ C] [◦ C]
0 -7.3 -7.3 -7.3 -7.3
15 -6.0 -6.8 -6.1 -7.0
30 -5.8 -6.9 -5.9 -6.9
45 -5.4 -6.9 -5.2 -6.9
60 -5.3 -6.8 -5.2 -6.9

Table 5.4: Experimental and numerical PCM temperature from the thermocouple placed
closest to the plate inside the PCM for Tin = 5 ◦ C and Tin = 0 ◦ C.

81
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 For the PCM temperature, temperature difference of up to 1.7 K is observed at
t =3600 s with the uncertainty of the thermocouples for the PCM temperature being
0.2 K. The temperature difference is assumed to be attributed by the condensation and
melting of the frost observed during the experimentation. As the numerical model was de-
signed for a significantly lower initial PCM temperature of −30 ◦ C, it accounted solely for
the desublimation as the phase change phenomena within the air channel. Condensation
is an exothermic process while melting of the frost is an endothermic process. There-
fore, higher experimental PCM temperature indicates that heat was released from the air
channel into the PCM, attributing condensation as the main cause for the temperature
difference between the experimental and numerical results. This also indicate that melting
of the frost observed during the experiment is most likely to be condensation and seeping
of the water vapours into the frost layers than the melting of the frost itself.

5.6 Results and discussion

As the continuation of the works of Jeong et al. [85] with the numerical simulation of the
thermoaeraulic behaviour of the air during the opening of the door of a refrigerated truck
trailer equipped with cold plates, frost formation on top of the eutectic plates and the
solidification and melting of the eutectic plates with the same initial PCM temperature
of −30 ◦ C has been simulated. The inlet conditions and the dimensions of the numerical
domain are displayed in Figure 5.1. Additionally, inlet temperature of Tinlet = 26.6 ◦ C
with identical relative humidity and velocity was simulated to compare the results. For
both cases, air channel is located to be below the PCM to promote the formation of the
Rayleigh-Bénard convection cells during the melting of the PCM, as observed from the
experiment in Section 5.4.
Figure 5.9(a) displays the average frost thickness of the numerical simulations. For the
inlet temperature of Tinlet = 21.6 ◦ C, sharp increase of the frost thickness between t = 0 s
to 60 s can be observed. Then the frost thickness steadily increases until t = 240 s and
demonstrate minimal increase afterwards. In contrast for the inlet temperature of Tinlet =
26.6 ◦ C, rather a gradual increase in the frost thickness is observed until t = 180 s. Then
a sharp increase in the frost thickness between t = 180 s to 240 s is observed with a linear
profile increase afterwards. Until t = 480 s, frost thickness for the inlet temperature of
Tinlet = 26.6 ◦ C is lower than that of Tinlet = 21.6 ◦ C. From Figure 4.3, same behavior was
observed between case 11 and case 12. While case 12 had 2.5 K higher inlet temperature,
until t = 600 s, numerically it showed lower frost thickness than that of case 11. This
behavior might be due to the higher initial temperature of the humid air increasing the
amount of time and energy it requires to cool the air to reach the supersaturated state for
the frost to form, as the relative humidity is the same. While this difference is minimal,
as the air velocity is high at uin = 10 m s−1 , the influence might be exaggerated, especially
during the early period of the frost formation. Additionally, as displayed in Figure 5.9(b),
the frost thickness increase rate between the two cases demonstrate similar profiles.
Due to the nature of the setup, melting of the PCM is focused on the lower region of the
PCM containment as shown in Figure 5.10. At t = 540 s, formation of the Rayleigh-Bénard
convection cells in the melted region of the PCM is clearly depicted. Furthermore, sides
of the PCM has been melted by the heat conducted through the surrounding aluminum
encasing. The rightmost side of the PCM has higher region of melted PCM than the
leftmost side as the airflow is from the right to the left. Figure 5.9(c) displays the linear
relationship between the melt fraction and time for both cases.

82
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 (a)

(b)

Figure 5.9: Numerical results for the (a) frost thickness, (b) rate of increase in frost thick-
ness, (c) PCM temperature, (d) melt fraction, and (e) outlet temperature.

For Tinlet = 21.6 ◦ C, it exhibited overall higher melt fraction and maximum PCM
temperature observed. For a typical door opening time of 900 s [65; 86], 16.34% of
the PCM would have melted from the linear approximation of the numerical results.
Even for delivery cycles with short door openings of between t =120 s to 180 s, 2.3 %
to 3.4 % of the PCM melt can be observed. From the maximum PCM temperature
shown in Figure 5.9(d), maximum PCM temperature of T =−17.99 ◦ C at t =180 s is
observed. As the maximum regulated temperature for quick deep frozen foodstuff is
also at T =−18 ◦ C, ineffectiveness of certain regions of the PCM to meet the regulated
temperature requirements due to the inhomogeneous melting of the PCM can be expected
even for a short door opening period [2].

83
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 (c)

(d)

Figure 5.9: Numerical results for the (a) frost thickness, (b) rate of increase in frost thick-
ness, (c) PCM temperature, (d) melt fraction, and (e) outlet temperature (cont.).

For Tinlet = 26.6 ◦ C, lower melt fraction and maximum PCM temperature were ob-
served. At t = 540 s, only approximately 4 % of the PCM melted. Also from t = 0 s
to 540 s, mere 3 K increase in the maximum PCM temperature was observed. Fig-
ure 5.9(e) displays the outlet temperature for Tinlet = 26.6 ◦ C being lower than that
of Tinlet = 21.6 ◦ C after t = 120 s. Additionally, at t = 540 s, density of the frost was
measured to be 574.53 kg m−3 and 324.33 kg m−3 for the cases with inlet temperature of
Tinlet = 21.6 ◦ C and Tinlet = 26.6 ◦ C, respectively. Therefore from the average frost thick-
ness, Tinlet = 26.6 ◦ C reported lower frost mass of 0.04 kg than for Tinlet = 21.6 ◦ C with
0.06 kg at t = 540 s for unit meter of width. This demonstrate that for Tinlet = 21.6 ◦ C,
majority of the PCM melt was due to the frost formation, observed by higher frost

84
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 (e)

Figure 5.9: Numerical results for the (a) frost thickness, (b) rate of increase in frost thick-
ness, (c) PCM temperature, (d) melt fraction, and (e) outlet temperature (cont.).

Figure 5.10: 2D contours of the melt fraction (top), and velocity magnitude (bot) at t = 540 s.

mass, higher outlet temperature and maximum PCM temperature. In contrast, for
Tinlet = 26.6 ◦ C as the frost formed on top of the PCM is lower until t = 480 s, most of the
PCM melt was due to the cooling of the air. As the outlet temperature for Tinlet = 26.6 ◦ C
is lower than for Tinlet = 21.6 ◦ C with 5 K higher inlet temperature and significantly less
PCM melt, the adverse outcome of the frost formation is exemplified.

5.7 Conclusion
In this chapter, the infiltration model and the frost model presented in Chapter 3 and
Chapter 4 were consolidated with the solidification and melting model applied as a UDF.
As the existing solidification and melting model within ANSYS FLUENT could not be
enabled with the Eulerian-Eulerian multiphase model used for the frost model, the so-
lidification and melting model was applied via a UDF. Since the UDF version of the
solidification and melting model incorporated with the infiltration and the frost models
might behave inaccurately, it has been favorably validated against the published experi-
mental and numerical results in terms of melt front location and melt fraction.
Additionally, experiment performed at the Université de Sherbrooke was presented

85
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 in this chapter. The experimental setup of Rahal [83] was adopted. However, with the
given initial and boundary conditions of the experimental setup, frost formation could
not be observed neither experimentally nor numerically. As the air channel is located on
the upper portion of the PCM container, it was numerically concluded that it results in
a highly localized high temperature region near the upper region of the PCM. Further-
more, there is little to no heat diffusion that occurs via convection apart from the minor
convection that occurs along the side walls of the PCM due to the temperature-driven
density difference caused by the heat transferred by the aluminum encasing. Experimen-
tally, a loss of contact between the PCM and the plate was observed, as the upper region
of the PCM melts and lead to the volumetric change associated with the phase change
of the PCM. Therefore, the experimental setup of Rahal [83] was adjusted so that the
air channel is located at the bottom of the PCM container. The obtained experimental
frost thickness was in agreement with the numerical results. However, as the maximum
frost thickness obtained during the experiment is 0.15 mm with 0.1 mm uncertainty, it is
difficult to conclude on the accuracy of the experimental results and its validation against
the numerical results. For the PCM temperature, experimental results showed maximum
temperature difference of up to 1.7 K after 3600 s. This temperature difference is assumed
to be attributed by the condensation observed during the experimentation. Numerically,
as the given initial temperature of the given eutectic system is −29 ◦ C with the time of the
exposure to the infiltrating hot atmospheric typically being less than t =120 s, condensa-
tion was neglected. Due to the limitation of the experimental setup, initial temperature of
the experimental setup could only be lowered to t =−10 ◦ C and resulted in the occurrence
of condensation.
Finally, the consolidated numerical model was applied to the typical initial and bound-
ary conditions experienced by a refrigerated truck trailer during humid summer. Two
different inlet temperature were simulated with Tinlet = 26.6 ◦ C and Tinlet = 21.6 ◦ C with
identical inlet relative humidity and velocity (69 % and 10 m s−1 , respectively). The re-
sults showed that until t = 480 s, frost thickness observed on top of the eutectic system
for Tinlet = 26.6 ◦ C being lower than that of Tinlet = 21.6 ◦ C. This is assumed to be
due to the higher inlet temperature requiring longer time and higher energy to reach
the supersaturated state, exaggerated by the high inlet velocity during the early period
of the frost formation. Due to the lower amount of frost formed on top of the plate,
Tinlet = 26.6 ◦ C exhibited significantly less PCM melt and lower maximum PCM tem-
perature observed. Furthermore, the outlet temperature for Tinlet = 26.6 ◦ C was lower
than that of Tinlet = 21.6 ◦ C between t = 120 s to 480 s. Therefore, the numerical results
exemplified the adverse outcome of the frost formation for the eutectic system.

86
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Chapter 6

Conclusions and future views

6.1 Conclusions
The main objective of this manuscript was to analyze the feasibility of the eutectic refrig-
eration system to replace the conventional refrigeration units operated by fossil fuels for
a refrigerated truck trailer for the transportation of the frozen foodstuff. Therefore, to
analyze the feasibility, a numerical model capable of accurately predicting the flow field,
temperature distributions, infiltration heat load into the system, formation of frost on
top of the plates, and the phase change of the eutectic mixture during the door opening
period for a truck trailer refrigerated with eutectic plate system was required. Firstly, an
infiltration model was developed to predict the thermoaeraulic behavior of the air during
the infiltration period for a refrigerated truck trailer. Then Eulerian-Eulerian granular
multiphase model was used to develop a frost model capable of predicting the growth of
the frost on top of a cold plate in terms of the frost thickness and density. The Eulerian-
Eulerian granular multiphase frost model incorporated the k − ω SST turbulence model
presented in the infiltration model to widen its applicable range of inlet air velocities.
Finally, the solidification and melting model was consolidated with the Eulerian-Eulerian
multiphase model as a UDF to predict the phase change of the eutectic mixture due to
the heat transferred to the PCM by the energy released during the frost formation and
the energy transferred from the infiltrating hot atmospheric air.
First, a numerical model has been developed and validated to predict the turbulent
flow and heat transfer inside a refrigerated truck trailer equipped with eutectic plates.
The model has been validated against the experimental data of Lafaye de Micheaux et
al. [60]. It was shown that the 3D effects did not influence significantly the velocity and
temperature distributions at the midplane of the doorway. The k − ω SST model slightly
improved the predictions of the realizable k − ϵ model of Lafaye de Micheaux et al. [60].
Based on the developed model, a 2D model was adapted to optimize the configuration of
the refrigeration system with the eutectic plates in terms of trailer and cargo temperatures
by analyzing the air flow and the infiltration heat load into the refrigerated truck trailer
during the door opening period. Without cargo, the configuration with the plates placed
in series on the roof of the trailer noticeably improved the performance in terms of trailer
temperature with respect to the configuration with the plates in series placed at the back.
However, the presence of a cargo eliminates the recirculation zones that prevented the
infiltration of the atmospheric air. Also, for the configuration with the plates placed on
the roof, a lower initial velocity of the pre-existing flow inside the trailer was observed
caused by the blockage of the flow due to the cargo. Therefore, the atmospheric air was

87
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 able to infiltrate and increase the temperature of the trailer quicker and resulted in a
higher maximum temperature observed for the cargo.
As the number of cargo boxes will progressively decrease during the delivery cycle, it
is plausible that the configuration with the plates placed on the roof could improve the
performance as the influence of the pre-existing flow seen without the cargo will progres-
sively emerge. The influence of the cargo parameters such as cargo locations, dimensions,
arrangements, etc, should be investigated into more details as it has a significant influ-
ence over the development of the pre-existing flow inside the trailer. Further development
should include a multiphase model to account for both the humidity of air and the melting
of the phase change material inside the plates. Another major concern is that the forma-
tion of frost along the plates could interfere with the heat transfer between the plates and
the surrounding air.
Then as a second step, a numerical model has been developed and validated to predict
the frost formation and the melting and solidification of the sub-zero eutectic plate. The
Eulerian-Eulerian multiphase frost model has been validated against published experi-
mental results [77; 75; 38; 78; 76]. As the solidification and melting model within the
ANSYS FLUENT cannot be activated with the Eulerian-Eulerian multiphase model, the
solidification and melting model has been incorporated into the specific heat equation by
a UDF into the multiphase model and validated against published numerical and exper-
imental results [11; 82]. Then the combined model was experimentally validated for the
frost thickness and the PCM temperature.
From both the numerical and experimental results, the key findings are as follows:

• If the airflow is above the PCM, heat is conducted to the PCM solely by conduction,
apart from the minor buoyancy driven convection at the sides of the PCM caused by
the temperature driven density difference by the heat transferred to the sides of the
PCM by the surrounding aluminum encasing. This results in a highly localized high
temperature zone at the upper region of the PCM, effectively nullifying its effect.
Furthermore, due to the volumetric change during the phase change of the PCM,
loss of contact between the PCM and the plate is observed. This would require
an additional mechanical system to address the issue and add onto the cost of the
system. To effectively diffuse the heat and utilize the PCM system, air flow should
be at the sides or beneath the PCM in order to promote the convection caused by the
temperature-driven density difference or by the formation of the Rayleigh-Bénard
convection cells.

• For the sub-zero eutectic plates at initial temperature of T = −30 ◦ C, for typical
summer delivery temperature of T = 21.6 ◦ C, relative humidity of 69 %, and veloc-
ity of 10 m s−1 caused by the buoyancy effect of the infiltration and the remaining
influence of the air fans, linear relationship between the PCM liquid fraction and
time can be observed. For typical door opening period of t = 900 s [65; 86], liquid
fraction of up to 0.1634 can be expected. Even for a short door opening period of
t = 120 s, 2.3% of liquid fraction can be expected. Additionally, at t = 180 s, max-
imum PCM temperature of T = −18 ◦ C is observed. As the maximum regulated
temperature for quick (deep) frozen foodstuffs is T = −18 ◦ C [2], ineffectiveness of
certain regions of the plate to meet the regulated temperature requirements due to
the inhomogeneous melting of the PCM can be expected. This addresses the impor-
tance of short door opening times for the delivery of the refrigerated frozen goods
and employment of heat-infiltration deterrence mechanisms such as air curtains.

88
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 • From the comparison of the results between the inlet temperatures of Tinlet = 21.6 ◦ C
and Tinlet = 26.6 ◦ C, it was demonstrated that the majority of the PCM melting
was due to frost formation. For Tinlet = 21.6 ◦ C higher frost mass and outlet tem-
perature were observed until t = 540 s, resulting in higher melting of the PCM than
that for Tinlet = 26.6 ◦ C. For Tinlet = 26.6 ◦ C, due to the higher time and energy
required to cool the air to reach the supersaturated state caused by the higher initial
temperature of the air, less frost was formed on top of the eutectic system than for
Tinlet = 21.6 ◦ C until t = 480 s. Therefore, most of the PCM melting was due to the
cooling of air, resulting in a lower outlet temperature than that of Tinlet = 21.6 ◦ C
even with the 5 K higher inlet temperature. Additionally, only 4 % of PCM was
melt for Tinlet = 26.6 ◦ C. Therefore, the adverse outcomes of the frost formation on
the eutectic system were exemplified.

Additionally, the Eulerian-Eulerian frost model showed difficulty when the tempera-
ture difference between the air and the cold plate was below 25K for the tested operating
conditions. The model can be further improved by optimizing the equation of the non-
dimensional frosting criterion coefficient, B. However, significant efforts are required to
repeat the simulations with varying coefficient values. It would also require many val-
idation cases with various operating conditions and configurations to be applied as a
generalized model for the frosting process. Furthermore, identifying the exact range of
Reynolds number for the application of the k − ω SST turbulence model will be a difficult
task.

In conclusion, eutectic system is capable of replacing the conventional refrigeration

systems operated by consuming fossil-fuels for the refrigerated truck trailers for the trans-
portation of frozen foodstuff. However, as the results indicated, there are number of
constraints. Primarily, the door opening period must be strictly controlled and be less
than t =180 s for the simulated conditions to meet the regulated temperature. This can
be further increased by the application of the infiltration deterrence mechanisms such as
air curtains. Secondly, eutectic plates places in parallel at the back of the trailer (Con-
figuration A in Chapter 3) is more suitable than the plates places in series along the roof
of the trailer (Configuration B in Chapter 3). While Configuration B performed slightly
better than Configuration A in terms of the renewal time without the cargo, the infiltra-
tion model in Chapter 3 assumed the plate to be a constant temperature plate. From
Chapter 5 it was shown that without convection, heat transfer between the PCM and
the surrounding system is primarily via conduction. For Configuration B, convection will
only occur along the bottom and sides of the eutectic plates. For the upper portion of
the plate, heat transfer will be solely due to conduction. Additionally, if a mechanism
accounting for the volumetric change of the PCM due to the phase change is not incor-
porated, there will be a loss of contact between the PCM and the top wall of the eutectic
plate. Therefore, upper portion of the eutectic plates in any configuration can be seen as
ineffective. As for the Configuration A, eutectic plates are placed in parallel at the back
of the trailer with the length of the upper portion of the plate being minimal. Therefore,
heat transfer between the PCM and the surrounding system will be far more effective for
the Configuration A.

89
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 6.2 Future Works and Challenges
The numerical model for the frost formation developed in this work currently relies on
the non-dimensional frosting criterion coefficient, B. This coefficient has been determined
through an iterative process of trial-and-error against the published experimental results.
Therefore, the accuracy of the numerical model heavily relies on the validity of the non-
dimensional frosting criterion coefficient. To improve the model, larger dataset containing
the frost thickness and density for a wide range of variation of input parameters, such
as inlet velocity, relative humidity, inlet air temperature, cold plate temperature, etc.,
are required. Given such dataset, the non-dimensional frosting criterion coefficient can
be further refined for a model with better accuracy and wider applicability. However,
feasibility of such task is subjected by the quantity of the available published experimental
results and the computational time and cost associated with the trial-and-error process.
Furthermore, the numerical model for the frost formation only considers desublima-
tion as the sole phase change phenomena occurring during the simulation. As the scope
of this research was for the transport of frozen foodstuff with the temperature of the re-
frigerated truck trailer set to be approximately −30 ◦ C according to the ATP requirement
[2], condensation was neglected. However, given the temperature range shown during the
experimentation, condensation inevitably occurs. This is further hindered by the lowered
accuracy of the frost model with the lower difference between the plate temperature and
the air. Therefore, condensation could also be incorporated into the model. However, as
the interaction between the condensed water droplets and the frost has to be accounted
for, the complexity and the computational cost of the numerical model will significantly
increase.
Finally, for the numerical model, the infiltration model developed to predict the ther-
moaeraulic behavior of air during the door opening period of a refrigerated truck trailer
and the frost formation and the solidification and melting model can be amalgamated to
create a complete numerical model of a refrigerated truck trailer equipped with eutectic
plates. This was the ultimate goal of the given research. However, there are several dif-
ficulties that must be overcome. Majorly, the grid size required for the frost formation
model and the infiltration is significantly different. While for the infiltration model, the
size of the grid cell must be sufficiently small enough to satisfy y+ < 1 for the k − ω
SST turbulence model, the grid size is required to be significantly smaller for the frost
formation model to be accurate in terms of the frost thickness and density. For the
2D infiltration model, 9.0 × 105 elements were required for the gird-independent results.
Therefore, the complete model will be significantly limited by the computational time
and resources available. As a solution, dynamic mesh could be incorporated to satisfy
the cell size conditions and reduce the computational cost. Alternatively, the numerical
model can be separated by the eutectic system and the refrigerated truck trailer, with a
coupled iterations and solutions from each time step. However, the computational time
and the complexities associated with either of the solutions still remains significantly
high. Furthermore, the inner walls of the trailer and the cargo will also accumulate frost.
However, the existing frost model only considers the average temperature of the eutec-
tic plate. Thus, it is necessary to separately incorporate the frost model for the walls
and cargo. This implementation will incur additional computational cost, as the cell size
requirements also apply to the areas adjacent to the walls and cargo.
Experimentally, equipment acquisition and availability will the be major concern. Ini-
tially, investigation of the humidity diffusion during the door opening period of a refriger-

90
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 ated truck trailer was planned for the research at INSA Lyon. However, the time-response
of the humidity sensor was a major concern. As the infiltration phenomena occurs within
the first five to ten seconds of the door opening, humidity sensor with the response time
of less than one second was required. While such humidity sensors were acquired, major
challenge was developing a data acquisition system with the given humidity sensors. The
acquisition system needed to communicate with each sensor and relay the acquired data
in a time frame much less than the response time of the humidity sensor in order to have
no-delay in the acquired data. Developing aforementioned data acquisition system was
outside the realm of this research, and required speciality of the professionals in electrical
engineering. But given the availability of the personnel, cost of the equipment, etc., the
experiment was postponed for the future researcher.
For the frost formation experiment performed at the Université de Sherbrooke, the
major challenge was the thermal capability of the Thermo Bath required to completely
freeze the PCM. Thermo Bath (Thermo Haake ARCTIC AC200 A40 Immersion Bath)
with 900 W cooling capacity at 20 ◦ C used the during the experiment was incapable of
freezing any PCM below the phase change temperature of T =−10 ◦ C. Initially, Ru-
bitherm SP-28 with the phase change temperature of T =−28 ◦ C was planned to be used
for the experiment. However, eutectic mixture of potassium bicarbonate was adopted due
to the limitation of the Thermo Bath. As the initial plate temperature is at minimum
T =−10 ◦ C, temperature difference between the plate and the air was limited. Due to
this, condensation was observed during the experimentation and it effected the valida-
tion process between the numerical model and the experiment as the numerical model
did incorporate the condensation effect. Ideally, Rubitherm SP-28 with the phase change
temperature of T =−28 ◦ C would have been used for the experiment, where no condensa-
tion would assumed to have occurred. As a solution, a smaller version of the experimental
setup with the Peltier elements as a cooling mechanism was proposed. If the Peltier el-
ements are used as the cooling mechanism, it does not require a coil-pipe system within
the PCM container unlike Thermo Bath. Therefore, the size of the PCM container will
not be limited. With the lowered volume of the PCM, required amount of heat to be
removed from the system will also be reduced. However, the acquirement of the Peltier
elements, development of the control systems, etc., would take significant amount of time.

91
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 92
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Chapter 7

Conclusions et perspectives (Version

Française)

7.1 Conclusions
L’objectif principal de ce manuscrit était d’analyser la faisabilité d’un système de réfrigéra-
tion à base de plaques eutectiques pour remplacer les unités de réfrigération convention-
nelles fonctionnant aux combustibles fossiles dans une remorque de camion réfrigérée
pour le transport de produits alimentaires surgelés. Pour cela, un modèle numérique ca-
pable de prédire avec précision le champ d’écoulement, la distribution de température, la
charge thermique due aux infiltrations, la formation de givre sur les plaques et le change-
ment de phase du mélange eutectique pendant la période d’ouverture des portes a été
développé. Un nouveau banc expérimental de plaque eutectique a également été con-
struit à l’Université de Sherbrooke afin de valider les modèles numériques de formation
de givre et de fusion du matériau à changement de phase.
Dans un premier temps, un modèle d’infiltration a été développé pour prédire le com-
portement thermoaéraulique de l’air pendant les périodes d’ouverture de la remorque. Il a
été soigneusement validé par rapport aux données expérimentales de Lafaye de Micheaux
et al. [60]. Il a été démontré que les effets 3D n’influencent pas de manière significative les
distributions de vitesse et de température au niveau du plan médian de la porte. Le mod-
èle k − ω SST a légèrement amélioré les prédictions du modèle réalisable k − ϵ de Lafaye
de Micheaux et al. [60]. Sur la base du modèle développé, un modèle 2D a été adapté
pour optimiser la configuration du système de réfrigération avec les plaques eutectiques en
termes de températures à l’intérieur de la remorque et de la cargaison en analysant le flux
d’air et la charge thermique d’infiltration dans la remorque frigorifique du camion pen-
dant la période d’ouverture des portes. Sans chargement, la configuration avec les plaques
placées horizontalement en série au plafond de la remorque a sensiblement amélioré les
performances en termes de température de la remorque par rapport à la configuration avec
les plaques en série placées verticalement à l’arrière. Cependant, la présence d’une cargai-
son élimine les zones de recirculation qui empêchent l’infiltration de l’air atmosphérique
lors des périodes d’ouverture. De plus, pour la configuration avec les plaques placées au
plafond, une vitesse initiale plus faible du flux préexistant à l’intérieur de la remorque
a été observée, causée par la présence la cargaison. Par conséquent, l’air atmosphérique
a pu s’infiltrer et augmenter la température de la remorque plus rapidement, ce qui a
entraîné une température maximale plus élevée pour la cargaison. Comme le nombre de
caisses diminuera progressivement au cours du cycle de livraison, il est plausible que la

93
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 configuration avec les plaques placées au plafond puisse améliorer les performances au fur
et à mesure que l’influence du flux préexistant sans la cargaison émerge progressivement.

Dans un deuxième temps, un second modèle numérique a été développé et validé pour
prédire la formation de givre au-dessus des plaques eutectiques ainsi que la fusion et la
solidification du matériau à changement de phase au sein de la plaque. Le modèle de
givre de type eulérien-eulérien a été validé par rapport à des résultats expérimentaux
disponibles dans la littérature [77; 75; 38; 78; 76] pour une large gamme de conditions
d’opération. Comme le modèle de solidification et de fusion d’ANSYS FLUENT n’est pas
compatible avec le modèle eulérien-eulérien, le modèle de solidification et de fusion a été
incorporé dans l’équation de la chaleur spécifique via une fonction utilisateur (UDF). Il
a ensuite été validé par rapport à d’autres données numériques et expérimentales de la
littérature [11; 82]. Enfin, le modèle combiné a été validé en termes d’épaisseur du givre et
de température du matériau à changement de phase par des données expérimentales issues
d’un nouveau banc expérimental construit spécifiquement à l’Université de Sherbrooke.
Les principales conclusions de cette seconde phase sont les suivantes:

• Si le flux d’air est au-dessus du MCP, la chaleur est conduite vers le MCP essen-
tiellement par conduction, à l’exception de la convection mineure générée sur les
côtés de la boîte. Il en résulte une zone de température élevée très localisée dans la
région supérieure du MCP. De plus, en raison du changement de volume induit par
le changement de phase du MCP, une perte de contact entre le MCP et la plaque
est observée. Pour diffuser efficacement la chaleur et utiliser pleinement le système,
le flux d’air doit être sur les côtés ou sous la plaque afin de favoriser les transferts
par convection.

• Pour des plaques eutectiques à une température initiale de T = −30 ◦ C et une tem-
pérature extérieure en été typique de T = 21.6 ◦ C, une humidité relative de 69 % et
une vitesse de l’air de 10 m s−1 , on peut obtenir une relation linéaire entre la fraction
de MCP qui a fondu et le temps. Pour une période typique d’ouverture de porte
de t = 900 s[65; 86], on peut s’attendre à une fraction fondue allant jusqu’à 0.163.
Même pour une courte période d’ouverture de porte de t = 120 s, on peut s’attendre
à 2.3% de MCP fondu. De plus, à t = 180 s, une température maximale du MCP
de T = −18 ◦ C est observée. Comme la température maximale acceptable pour les
aliments surgelés est de T = −18 ◦ C [2], les inhomogénéités de température au sein
de la plaque peuvent rendre le système inefficace pour garantir cette température
au sein de la remorque. Cela démontre l’importance d’avoir des temps d’ouverture
des portes courts lors de la livraison des produits surgelés et/ou d’employer des
systèmes additionnels pour isoler la remorque de l’extérieur, comme des rideaux en
plastique ou des rideaux d’air.

• De la comparaison des résultats entre les températures d’entrée de Tinlet = 21.6 ◦ C

et Tinlet = 26.6 ◦ C, il a été démontré que la majorité de la fonte du MCP était due à
la formation de givre. Pour Tinlet = 21.6 ◦ C, une masse de givre et une température
de sortie plus élevées ont été observées jusqu’à t = 540 s, ce qui a entraîné une fusion
du MCP plus élevée que celle de Tentre = 26.6 ◦ C. Pour Tinlet = 26.6 ◦ C, en raison
du temps et de l’énergie plus élevés nécessaires pour refroidir l’air afin d’atteindre
l’état sursaturé provoqué par la température initiale plus élevée de l’air, du givre
s’est formé au-dessus du système eutectique, mais son épaisseur était inférieure à

94
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 celle obtenue pour Tinlet = 21.6 ◦ C jusqu’à t = 480 s. Par conséquent, la majeure
partie de la fusion du MCP était due au refroidissement de l’air, ce qui faisait que
la température de sortie était inférieure à celle de Tinlet = 21.6 ◦ C même avec le 5 K
température d’entrée plus élevé. De plus, seulement 4 % du MCP a fondu pour
Tinlet = 26.6 ◦ C. Ainsi, les conséquences néfastes de la formation de givre sur le
système eutectique ont été illustrées.

De plus, le modèle de givre eulérien-eulérien a montré des difficultés à prédire de façon

précise l’épaisseur de givre lorsque la différence de température entre l’air et la plaque
froide était inférieure à 25K pour les conditions de fonctionnement testées. Le modèle
peut être encore amélioré en optimisant l’équation du coefficient du critère de givrage
adimensionnel, B. Cependant, des efforts importants sont nécessaires pour répéter les
simulations avec des valeurs variables de ce coefficient. Cela nécessiterait également de
nombreux cas de validation avec diverses conditions de fonctionnement pour être appliqués
comme modèle généralisé pour le processus de givrage.
En conclusion, le système eutectique est capable de remplacer les systèmes de réfrigéra-
tion conventionnels dans les remorques frigorifiques destinées au transport de produits
surgelés. Cependant, comme l’indiquent les résultats, il existe un certain nombre de con-
traintes. Principalement, la période d’ouverture de la porte doit être strictement contrôlée
et être inférieure à t =180 s pour respecter la température cible au sein de la remorque.
Les infiltrations peuvent être réduites par l’ajout de systèmes comme des rideaux en
plastique ou des rideaux d’air. Deuxièmement, les plaques eutectiques placées en paral-
lèle à l’arrière de la remorque (Configuration A au Chapitre 3) sont plus adaptées que les
plaques placées en série le long du plafond de la remorque (Configuration B au Chapitre 3).
Alors que la configuration B a donné des résultats légèrement meilleurs que la configura-
tion A en termes de temps de renouvellement sans la cargaison, le modèle d’infiltration
du Chapitre 3 a supposé que la plaque était une plaque à température constante. Au
chapitre 5, il a été montré que sans convection, le transfert de chaleur entre le MCP et
le système environnant se fait principalement par conduction. Pour la configuration B,
la convection ne se produira que le long du fond et des côtés des plaques eutectiques.
Pour la partie supérieure de la plaque, le transfert de chaleur se fera uniquement par
conduction. De plus, si un mécanisme prenant en compte le changement volumétrique
du MCP dû au changement de phase n’est pas incorporé, il y aura une perte de contact
entre le MCP et la paroi supérieure de la plaque eutectique. Par conséquent, la partie
supérieure des plaques eutectiques, quelle que soit leur configuration, peut être considérée
comme inefficace. Dans la configuration A, les plaques eutectiques sont placées verticale-
ment parallèlement à l’arrière de la remorque. La longueur de la partie supérieure de la
plaque étant minimale, le transfert de chaleur entre le MCP et l’air dans la remorque sera
beaucoup plus efficace que pour la configuration A.

7.2 Perspectives de travail

Concernant le modèle d’infiltration développé pour prédire le comportement thermoaéraulique
de l’air pendant la période d’ouverture des portes de la remorque, l’influence des paramètres
liés à la cargaison (emplacement, dimensions, disposition, etc.) devrait être étudiée plus
en détail car elle a une influence significative sur le développement du flux préexistant
à l’intérieur de la remorque. L’évolution temporelle de la température du chargement
reste, de plus, le paramètre principal à prédire pour vérifier qu’elle reste dans une plage

95
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 acceptable par rapport aux exigences ATP [2]. Le rayonnement entre surfaces au sein de
la remorque, bien qu’à priori faible, pourrait être également intégré au modèle final.
Le modèle numérique de formation de givre développé dans ce travail repose actuelle-
ment sur un coefficient de critère de givrage non dimensionnel, B. Ce coefficient a été
déterminé par un processus itératif d’essais et erreurs par rapport à une série de résultats
expérimentaux. La précision du modèle numérique dépend fortement de la validité de ce
coefficient B. Pour améliorer le modèle, un ensemble de données plus grand contenant
l’épaisseur et la densité du givre pour une large gamme des paramètres d’opération, comme
la vitesse d’entrée, l’humidité relative, la température de l’air d’entrée, la température de
la plaque froide, etc, est nécessaire afin d’affiner le coefficient B.
De plus, le présent modèle numérique de formation de givre ne considère que la dé-
sublimation comme seul phénomène de changement de phase. Comme le champ de cette
recherche concernait le transport de produits alimentaires surgelés avec une température
de la remorque du camion réfrigérée fixée à environ −30 ◦ C selon l’exigence ATP [2], la
condensation a été négligée. Cependant, étant donné la plage de température mesurée
lors des expériences, de la condensation se produit inévitablement. Ceci est encore en-
travé par la précision réduite du modèle de givre, avec une différence plus faible entre
la température de la plaque et celle de l’air. La condensation pourrait donc également
être intégrée au modèle. Cependant, comme les interactions entre les gouttelettes d’eau
condensée et le givre doivent être prise en compte, la complexité et le coût de calcul du
modèle numérique vont augmenter considérablement.
Enfin, un modèle complet incorporant les différents modèles déjà développés représen-
tait l’objectif principal de cette étude. Cependant, plusieurs difficultés doivent encore
être surmontées. La principale est due à la différence d’échelles entre les phénomènes
aéroliques à l’échelle de la remorque et la formation de givre à l’échelle millimétrique.
Cela nécessite des maillages extrêmement fins et des coûts de calcul [Link] guise
de solution, un maillage dynamique pourrait être incorporé pour satisfaire les conditions
de taille de cellule et réduire le coût de calcul. Sinon, le modèle numérique peut être
séparé en deux parties, une pour le système eutectique et une pour la remorque, avec des
itérations couplées. Cependant, le temps de calcul et les complexités associées à l’une ou
l’autre des solutions restent encore significativement élevés. De plus, les parois intérieures
de la remorque et de la cargaison vont également accumuler du givre. Cependant, le
modèle de givre existant ne considère que la température moyenne de la plaque eutec-
tique. Ainsi, il est nécessaire d’incorporer séparément le modèle de givre pour les parois
et la cargaison. Cette mise en œuvre entraînera un coût de calcul supplémentaire, car les
exigences de taille de cellule s’appliquent également aux zones adjacentes aux parois de
la remorque et de la cargaison.
Expérimentalement, une étude sur la diffusion de l’humidité pendant les périodes
d’ouverture des portes de la remorque était prévue initialement à l’INSA de Lyon. Cepen-
dant, la réponse temporelle du capteur d’humidité constituait une préoccupation majeure.
Étant donné que le phénomène d’infiltration se produit dans les cinq à dix secondes suiv-
ant l’ouverture des portes, un capteur d’humidité avec un temps de réponse inférieur à une
seconde était nécessaire. Le défi majeur consistait à développer le système d’acquisition
de données correspondant. Le système d’acquisition devait communiquer avec chaque
capteur et relayer les données acquises dans un laps de temps bien inférieur au temps
de réponse des capteurs d’humidité afin de n’avoir aucun retard. Le développement du
système d’acquisition de données reste à développer avant d’initier cette étude.
Pour l’expérience sur la formation de givre réalisée à l’Université de Sherbrooke, le

96
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 défi majeur était la puissance du bain thermostaté requise pour solidifier complètement
le MCP. Le bain d’immersion (Thermo Haake ARCTIC AC200 A40) ayant une capacité
de refroidissement de 900 W à 20 ◦ C était incapable de solidifier un MCP en dessous de la
température de changement de phase de T =−10 ◦ C. Initialement, il était prévu d’utiliser
le Rubitherm SP-28 avec une température de changement de phase de T =−28 ◦ C pour
l’expérience. Cependant, un mélange eutectique de bicarbonate de potassium a été adopté
en raison des limites du bain thermique. Comme la température initiale de la plaque est
au minimum de T =−10 ◦ C, la différence de température entre la plaque et l’air reste lim-
itée. Pour cette raison, de la condensation a été observée au cours de l’expérimentation.
Idéalement, le MCP Rubitherm SP-28 avec une température de changement de phase de
T =−28 ◦ C aurait été utilisé pour l’expérience. Comme solution, une version plus petite
du dispositif expérimental avec des éléments Peltier comme mécanisme de refroidissement
a été proposée. Si des éléments Peltier sont utilisés, ils ne nécessitent pas de système
de serpentin-tuyau dans le conteneur MCP contrairement à la solution avec le bain ther-
mostaté. Par conséquent, la taille du conteneur MCP ne sera pas limitée. Avec le volume
réduit du MCP, la quantité de chaleur requise à éliminer du système sera également ré-
duite. Ce nouveau banc expérimental serait intéressant à développer afin d’étendre la
campagne de mesures à d’autres conditions d’opération tout en mesurant la masse de
givre produite.

97
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 98
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Appendix A

List of Publications
Journal Article:

• J. Jeong, A.E. Benchikh Le Hocine, S. Croquer, S. Poncet, B. Michel, J. Bonjour.

Numerical analysis of the thermoaeraulic behavior of air during the opening of the
door of a refrigerated truck trailer equipped with cold plates. Applied Thermal
Engineering, 206, 118057, 2022.

Conferences:

• J. Jeong, S. Poncet, B. Michel, J. Bonjour. Numerical simulation of the frost forma-

tion on a flat plate cooled by a phase change material, 26th International Congress
of Refrigeration (ICR23), Paris, 2023.

• J. Jeong, S. Poncet, B. Michel, J. Bonjour. Eulerian-Eulerian Multiphase Frost

Model Based on Phase Change Driving Force. 7th IIR International Conference on
Sustainability and the Cold Chain, Newcastle, 2022.

• J. Jeong, A.E. Benchikh Le Hocine, S. Croquer, S. Poncet, J. Bonjour, B. Michel.

Numerical Simulation of the Heat Transfer in a Refrigerated Trailer Equipped with
Eutectic Plates for Frozen Food Delivery. 18th International Refrigeration and Air
Conditioning Conference, Lafayette, 2021.

99
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 100
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 List of figures

1.1 Example of an eutectic refrigeration system inside a refrigerated truck

trailer. Courtesy of FrygyCube. . . . . . . . . . . . . . . . . . . . . . . . . 16

2.1 Schematic diagram of vapour compression transport refrigeration unit driven

by a diesel engine [6]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.2 A schematic of sensible and latent heat storage principle [7]. . . . . . . . . 21
2.3 Overview of PCMs [10]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.4 Experimental and predicted melt front locations with time for different
values of Amush indicated on the top [11]. . . . . . . . . . . . . . . . . . . 24
2.5 Schematic diagram of the physical processed involved in frost develop-
ment [29]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.6 Schematic diagram for the computing domain in the frosting process [46]. . 28
2.7 Predicted infiltration for the 2.3 m wide entrance for different door opening
times [49]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.8 Deformation of temperature (left) and velocity (right) profiles with respect
to the exterior temperature [50]. . . . . . . . . . . . . . . . . . . . . . . . . 31
2.9 Evolution of the infiltration flow rate for the baseline case [60]. . . . . . . . 33
2.10 The temperature and velocity profiles at the truck trailer door vertical
centerline [60]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.11 Total and sensible infiltration heat loads compared to energy received by
the internal air volume [60]. . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3.1 (a) Trailer diagram with the position of the eutectic plates for Configura-
tion A (top) and Configuration B (bottom); (b) General schematics of the
computational domain with the relevant boundary conditions. . . . . . . . 40
3.2 Details and dimensions of the eutectic plates. . . . . . . . . . . . . . . . . 41
3.3 Overview of the mesh grids for Configuration A (top) and Configuration B
(bottom) with a focus on the fan and the plate regions on the right. . . . . 42
3.4 Comparison of the velocity magnitude profiles at the container door vertical
centerline between the present simulations and the published data of Lafaye
de Micheaux et al. [60] at (a) t = 10 s, (b) 20 s, (c) 40 s and (d) 60 s. . . . . 44
3.5 Comparison of the temperature profiles at the container door vertical cen-
terline between the present simulations and the published data of Lafaye
de Micheaux et al. [60] at (a) t = 10 s, (b) 20 s, (c) 40 s and (d) 60 s. . . . . 45
3.6 Comparison in terms of the infiltration rate at the container door between
the present simulations and the experimental and numerical data of Lafaye
de Micheaux et al. [60]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

101
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 3.7 (a) Temperature [K] and (b) velocity magnitude [ms−1 ] contours for Con-
figurations A and B. The blowing mode is always on the top and the suction
mode at the bottom. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.8 Temporal evolutions of the (a) infiltration rate and (b) heat load through
the door for the two configurations and modes. . . . . . . . . . . . . . . . . 48
3.9 Time evolution of the normalized infiltration compared against the CFD
results obtained from Foster et al. (Extended: extended boundary CFD
model) [49]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.10 (top) Time evolution of the temperature contours [K] for both configu-
rations and fan modes at five instants; (bottom) Corresponding velocity
magnitude contours [ms−1 ] at t = 10 s. . . . . . . . . . . . . . . . . . . . . 50
3.11 (a) Temperature [K] and (b) velocity magnitude [ms−1 ] contours of Con-
figuration A with cargo (top) and Configuration B with cargo (bottom). . . 52
3.12 Temperature contours [K] for Configurations A - suction (top) and B -
blowing (bottom) at t = 2.5, 5, 10 and 20 s (from left to right). . . . . . . . 53
3.13 Time evolution of the area averaged temperature inside the trailer. . . . . . 54

4.1 Numerical domain for the frost formation on a cold plate. . . . . . . . . . . 58

4.2 Numerical results for Case No.12: (a) Volume fraction of ice at t = 3600 s,
(b) Temperature distribution (K) at t = 3600 s, and (c) Air velocity
distribution(ms−1 ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.3 Results for the frost thickness and density at t = 3600 s. . . . . . . . . . . 64
4.3 Results for the frost thickness and density at t = 3600 s (cont.). . . . . . . 65
4.4 Comparison between experimental data and numerical results for the (a)
Frost thickness and (b) frost density. . . . . . . . . . . . . . . . . . . . . . 66

5.1 (a) Numerical/experimental schematic of the aluminum plate containing

the PCM and the air channel below; (b) Dimensions of the experimen-
tal/numerical domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.2 Melt front of the Solidification and Melting model for Amush = 106 , 5 × 106
and 107 at (a) t = 600 s, (b) t = 1200 s, (c) t = 2400 s, (d) t = 3600 s. . . . 73
5.2 Melt front of the Solidification and Melting model for Amush = 106 , 5 × 106
and 107 at (a) t = 600 s, (b) t = 1200 s, (c) t = 2400 s, (d) t = 3600 s (cont.). 74
5.3 Melt front of the UDF Solidification and Melting model for Amush = 106
at (a) t = 600 s, (b) t = 1200 s, (c) t = 2400 s, (d) t = 3600 s. . . . . . . . . 75
5.3 Melt front of the UDF Solidification and Melting model for Amush = 106
at (a) t = 600 s, (b) t = 1200 s, (c) t = 2400 s, (d) t = 3600 s. (cont.) . . . 76
5.4 (a) Melt fraction for different values of Amush using the UDF model; (b)Average
melt front y location for two applied solidification and melting schemes and
Amush = 106 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.5 Velocity magnitude of the melted PCM for Amush = 106 at four different
times. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
5.6 (a) CAD view of the experimental setup, and zooms on the (b) aluminum
encasing, (c) thermocouples locations inside the aluminium encasing, (d)
aluminum panels of varying sizes for measuring the frost weight. . . . . . . 78
5.7 2D contours of the PCM melt fraction at four different times. . . . . . . . 79
5.8 Numerical results for the (a) maximum plate temperature at t = 1800 s, (b)
maximum PCM temperature over time and (c) velocity magnitude contours
of the PCM at t = 1800 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

102
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 5.9 Numerical results for the (a) frost thickness, (b) rate of increase in frost
thickness, (c) PCM temperature, (d) melt fraction, and (e) outlet temper-
ature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
5.9 Numerical results for the (a) frost thickness, (b) rate of increase in frost
thickness, (c) PCM temperature, (d) melt fraction, and (e) outlet temper-
ature (cont.). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
5.9 Numerical results for the (a) frost thickness, (b) rate of increase in frost
thickness, (c) PCM temperature, (d) melt fraction, and (e) outlet temper-
ature (cont.). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
5.10 2D contours of the melt fraction (top), and velocity magnitude (bot) at t =
540 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

103
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 104
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 List of tables

2.1 Summary of recent investigations on rectangular PCM solidification and

melting problems. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2 Summary of the frost model based on the mass transfer rate. . . . . . . . . 29
2.3 Comparison of the expressions for the infiltration rate used in various an-
alytical models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.4 Summary of recent studies focusing on the performance of refrigerated ve-
hicles. Symbol: * Cargo present/Door openings. . . . . . . . . . . . . . . . 36

3.1 Container and eutectic plate dimensions. . . . . . . . . . . . . . . . . . . . 41

3.2 Temperature and velocity magnitude obtained from the steady-state sim-
ulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.3 Infiltration data from the transient simulations without cargo. . . . . . . . 50
3.4 Thermophysical properties of the cargo. . . . . . . . . . . . . . . . . . . . . 52
3.5 Temperature and velocity magnitude data from the steady-state simula-
tions with cargo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.6 Infiltration data for Configurations A and B with and without cargo. . . . 53
3.7 Maximum temperature [K] of the cargo at t = 40 s. . . . . . . . . . . . . . 54

4.1 Frosting conditions used to validate the frost model. . . . . . . . . . . . . . 62

5.1 Properties of Rubiterm SP-28 (PCM). . . . . . . . . . . . . . . . . . . . . 70

5.2 Uncertainty of the equipment. . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.3 Experimental and numerical frost thickness for Tin = 5 ◦ C and Tin = 0 ◦ C. . 81
5.4 Experimental and numerical PCM temperature from the thermocouple
placed closest to the plate inside the PCM for Tin = 5 ◦ C and Tin = 0 ◦ C. . 81

nomenclature[z-URANS]URANSUnsteady Reynolds-averaged Navier–Stokes

105
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 106
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Nomenclature

Roman Symbols P Pressure [Pa]

A Cross sectional area of the door [m2 ] P r Prandtl number [-]

Amush Mushy zone constant [-] Q̄ Interphase heat transfer [W m−3 ]

B Non-dimensional frosting criterion Q Heat [J]

coefficient [-]
Q̇ Heat load [W]
b Thickness of the door frame [m] ⃗q Heat flux [W m−2 ]
Cp Specific heat capacity [J kg−1 K−1 ] R Specific gas constant [287 J kg−1 K−1 ]
d Diameter [m] Re Reynolds number [-]
DH2 O Diffusivity of the water vapour r̄ Droplet radius [m]
[m2 s−1 ]
r∗ Kelvin-Helmholtz critical radius [m]
E Heat transfer coefficient [W m K ]
−2 −1
S Source term
f Drag function [-]
T Temperature [K]
g Gravitational Acceleration [9.81 m s ]
−2
t Time [s]
H Height [m]
V Volume [m3 ]
h Enthalpy [J kg ]−1
v Velocity [m s−1 ]
I Infiltration [m ]
3
W Width [m]
I˙ Infiltration rate [m s ]
3 −1
w Mass fraction [-]
K Momentum transfer coefficient [-] y+ Non-dimensional wall distance [-]
Kf,L Fritsche and Lilienblum correction Greek Symbols
factor [-]
α Volume fraction [-]
k Turbulence kinetic energy [m2 s−2 ]
δ Frost thickness [mm]
L Length [m]
δr̄
δt
Nucleation growth rate
m Mass [kg]
ϵ Dissipation rate of k [m2 s−3 ]
ṁ Mass transfer rate [kg s ]
−1
η Number of droplets per unit volume
N u Nusselt number [-] [m−3 ]

107
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 λ Thermal conductivity [W m−1 K−1 ] sens Sensible

µ Dynamic viscosity [Pa s−1 ] s Solid

ρ Density [kg m−3 ] u Momentum

τ̄ Stress-strain tensor [Pa] va Water vapour

τv Time relaxation coefficient [-] vs Saturated water vapour

ϕ Liquid fraction [-] v Specie

χ Small number, 0.001 w Wall

ω Specific turbulence dissipation rate Acronyms / Abbreviations

[s−1 ]
ATP Agreement on the International Car-
Subscripts riage of Perishable Foodstuffs and on
the Special Equipment to be used for
0 Initial such Carriage

a Air CFD Computational Fluid Dynamics

ac Air to ice COP Coefficient of Performance

avg Average CTE Coefficient of Thermal Expansion

[K−1 ]
ca Ice to air
PCM Phase Change Material
c Ice
PU Polyurethane
da Dry air
RH Relative Humidity
h Energy
RMS Root Mean Square
i Inner
RSM Reynolds Stress Model
in Inlet
SST Shear Stress Transport
L Latent
UDF User Defined Function
l Liquid
VOF Volume of Fluid
m Mass

n Normal

o Outer

s Sensible

total Total

p Plate

ref Reference

108
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 Bibliography

[1] S. Tassou, J. Lewis, Y. Ge, A. Hadawey, I. Chaer, A review of emerging technologies

for food refrigeration applications, Applied Thermal Engineering 30 (4) (2010) 263–
276.
[2] UN Economic Commission for Europe – Inland Transport Committee, ATP as
amended on 6 January 2018 – Agreement on the International Carriage of Perishable
Foodstuff and on the Special Equipment to be used for such Carriage, Tech. rep.
(2018).
[3] S. Tassou, G. De-Lille, Y. Ge, Food transport refrigeration–approaches to reduce
energy consumption and environmental impacts of road transport, Applied Thermal
Engineering 29 (8-9) (2009) 1467–1477.
[4] S. Croquer, A. Benchikh Le Hocine, S. Poncet, Numerical modelling of heat and
mass transfer in a refrigerated truck trailer, in: Proc. of the 25th IIR International
Congress of Refrigeration, Montreal, 2019, pp. 3574–3581.
[5] X. Wu, Q. Ma, F. Chu, S. Hu, Phase change mass transfer model for frost growth
and densification, International Journal of Heat and Mass Transfer 96 (2016) 11–19.
[6] A. Rai, S. Tassou, Environmental impacts of vapour compression and cryogenic trans-
port refrigeration technologies for temperature controlled food distribution, Energy
Conversion and Management 150 (2017) 914–923.
[7] T. Goswami, V. Raj, Use of phase change material (PCM) for the improvement of
thermal performance of cold storage, MOJ Current Research & Reviews 1 (2018)
49–61.
[8] R. Chaturvedi, A. Islam, K. Sharma, A review on the applications of PCM in thermal
storage of solar energy, Materials Today: Proceedings 43 (2021) 293–297.
[9] P. Singh, R. Sharma, A. Ansu, R. Goyal, A. Sari, V. Tyagi, A comprehensive review
on development of eutectic organic phase change materials and their composites for
low and medium range thermal energy storage applications, Solar Energy Materials
and Solar Cells 223 (2021) 110955.
[10] M. Mofijur, T. Mahlia, A. Silitonga, H. Ong, M. Silakhori, M. Hasan, N. Putra,
S. Rahman, Phase change materials (PCM) for solar energy usages and storage: an
overview, Energies 12 (16) (2019) 3167.
[11] M. Fadl, P. Eames, Numerical investigation of the influence of mushy zone parameter
amush on heat transfer characteristics in vertically and horizontally oriented thermal
energy storage systems, Applied Thermal Engineering 151 (2019) 90–99.

109
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 [12] A. Brent, V. Voller, K. Reid, Enthalpy-porosity technique for modeling convection-
diffusion phase change: application to the melting of a pure metal, Numerical Heat
Transfer, Part A Applications 13 (3) (1988) 297–318.

[13] ANSYS FLUENT 12.0, Theory Guide (2009) 67.

[14] O. Zarajabad, R. Ahmadi, Numerical investigation of different PCM volume on cold

thermal energy storage system, Journal of Energy Storage 17 (2018) 515–524.

[15] J. Gou, J. Lin, Y. Jiang, S. Mei, J. Xia, S. Lin, Q. Shen, X. Zhang, W. Cai, J. Liang,
et al., Numerical analysis of cold energy release process of cold storage plate in a
container for temperature control, Journal of Energy Storage 71 (2023) 108230.

[16] M. Calati, G. Righetti, C. Zilio, K. Hooman, S. Mancin, CFD analyses for the de-
velopment of an innovative latent thermal energy storage for food transportation,
International Journal of Thermofluids 17 (2023) 100301.

[17] V. Voller, C. Prakash, A fixed grid numerical modelling methodology for convection-
diffusion mushy region phase-change problems, International Journal of Heat and
Mass Transfer 30 (8) (1987) 1709–1719.

[18] J. Vogel, J. Felbinger, M. Johnson, Natural convection in high temperature flat plate
latent heat thermal energy storage systems, Applied Energy 184 (2016) 184–196.

[19] Y. Kim, A. Hossain, S. Kim, Y. Nakamura, A numerical study on time-dependent

melting and deformation processes of Phase Change Material (PCM) induced by
localized thermal input, INTECH Open Access Publisher, 2011.

[20] R. Oliveski, F. Becker, L. Rocha, C. Biserni, G. Eberhardt, Design of fin structures

for phase change material (PCM) melting process in rectangular cavities, Journal of
Energy Storage 35 (2021) 102337.

[21] A. Abdulmunem, P. Samin, H. Rahman, H. Hussien, I. Mazali, H. Ghazali, Numerical

and experimental analysis of the tilt angle’s effects on the characteristics of the
melting process of PCM-based as PV cell’s backside heat sink, Renewable Energy
173 (2021) 520–530.

[22] H. Shokouhmand, B. Kamkari, Experimental investigation on melting heat transfer

characteristics of lauric acid in a rectangular thermal storage unit, Experimental
Thermal and Fluid Science 50 (2013) 201–212.

[23] A. Tinti, A. Tarzia, A. Passaro, R. Angiuli, Thermographic analysis of polyurethane

foams integrated with phase change materials designed for dynamic thermal insula-
tion in refrigerated transport, Applied Thermal Engineering 70 (1) (2014) 201–210.

[24] B. Copertaro, P. Principi, R. Fioretti, Thermal performance analysis of PCM in

refrigerated container envelopes in the Italian context - Numerical modeling and
validation, Applied Thermal Engineering 102 (2016) 873–881.

[25] B. Michel, P. Glouannec, A. Fuentes, P. Chauvelon, Experimental and numerical

study of insulation walls containing a composite layer of PU-PCM and dedicated to
refrigerated vehicle, Applied Thermal Engineering 116 (2017) 382–391.

110
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 [26] M. Calati, C. Zilio, G. Righetti, G. Longo, K. Hooman, S. Mancin, A numerical
analysis of latent thermal energy storage for refrigerated trucks, in: International
Refrigeration and Air Conditioning Conference, Purdue, 2021.

[27] M. Berdja, A. Hamid, O. Sari, Characteristics and thickness effect of phase change
material and frost on heat transfer and thermal performance of conventional refrig-
erator: Theoretical and experimental investigation, International Journal of Refrig-
eration 97 (2019) 108–123.

[28] Y. Hayashi, A. Aoki, S. Adachi, K. Hori, Study of frost properties correlating with
frost formation types, Journal of Heat Transfer 99 (2) (1977) 239–245.

[29] D. Kim, C. Kim, K. Lee, Frosting model for predicting macroscopic and local frost
behaviors on a cold plate, International Journal of Heat and Mass Transfer 82 (2015)
135–142.

[30] Y. Tao, R. Besant, K. Rezkallah, A mathematical model for predicting the densifi-
cation and growth of frost on a flat plate, International Journal of Heat and Mass
Transfer 36 (2) (1993) 353–363.

[31] B. Jones, J. Parker, Frost formation with varying environmental parameters, Journal
of Heat Transfer 97 (1975) 255–259.

[32] D. O’Neal, D. Tree, Measurement of frost growth and density in a parallel plate
geometry, ASHRAE Transactions 90 (2) (1984) 278–290.

[33] B. Na, R. Webb, New model for frost growth rate, International Journal of Heat and
Mass Transfer 47 (5) (2004) 925–936.

[34] K. Lee, S. Jhee, D. Yang, Prediction of the frost formation on a cold flat surface,
International Journal of heat and Mass Transfer 46 (20) (2003) 3789–3796.

[35] Y. Yao, Y. Jiang, S. Deng, Z. Ma, A study on the performance of the airside heat
exchanger under frosting in an air source heat pump water heater/chiller unit, Inter-
national Journal of Heat and Mass Transfer 47 (17-18) (2004) 3745–3756.

[36] A. Léoni, M. Mondot, F. Durier, R. Revellin, P. Haberschill, State-of-the-art review

of frost deposition on flat surfaces, International Journal of Refrigeration 68 (2016)
198–217.

[37] B. Na, R. Webb, Mass transfer on and within a frost layer, International Journal of
Heat and Mass Transfer 47 (5) (2004) 899–911.

[38] K.-S. Lee, W.-S. Kim, T.-H. Lee, A one-dimensional model for frost formation on
a cold flat surface, International Journal of Heat and Mass Transfer 40 (18) (1997)
4359–4365.

[39] R. Barron, L. Han, Heat and mass transfer to a cryosurface in free convection, Journal
of Heat Transfer 87 (4) (1965) 499–506.

[40] K. Lenic, A. Trp, B. Frankovic, Transient two-dimensional model of frost formation

on a fin-and-tube heat exchanger, International Journal of Heat and Mass Transfer
52 (1-2) (2009) 22–32.

111
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 [41] D. Yang, K. Lee, Modeling of frosting behavior on a cold plate, International Journal
of Refrigeration 28 (3) (2005) 396–402.

[42] D. Yang, K. Lee, D. Cha, Frost formation on a cold surface under turbulent flow,
International Journal of Refrigeration 29 (2) (2006) 164–169.

[43] J. Cui, W. Li, Y. Liu, Z. Jiang, A new time-and space-dependent model for predicting
frost formation, Applied Thermal Engineering 31 (4) (2011) 447–457.

[44] J. Cui, W. Li, Y. Liu, Y. Zhao, A new model for predicting performance of fin-and-
tube heat exchanger under frost condition, International Journal of Heat and Fluid
Flow 32 (1) (2011) 249–260.

[45] D. Zhuang, G. Ding, H. Hu, H. Fujino, S. Inoue, Condensing droplet behaviors on fin
surface under dehumidifying condition: Part I: Numerical model, Applied Thermal
Engineering 105 (2016) 336–344.

[46] X. Wu, F. Chu, Q. Ma, Frosting model based on phase change driving force, Inter-
national Journal of Heat and Mass Transfer 110 (2017) 760–767.

[47] E. Afrasiabian, O. Iliev, S. Lazzari, C. Isetti, Numerical simulation of frost formation

on a plate-fin evaporator, in: Proceedings of the 3rd World Congress on Momentum,
Heat and Mass Transfer (MHMT’18), Budapest, 2018, pp. 12–14.

[48] A. Sommers, N. Truster, A. Napora, A. Riechman, E. Caraballo, Densification of

frost on hydrophilic and hydrophobic substrates–examining the effect of surface wet-
tability, Experimental Thermal and Fluid Science 75 (2016) 25–34.

[49] A. Foster, M. Swain, R. Barrett, S. James, Experimental verification of analytical and

CFD predictions of infiltration through cold store entrances, International Journal
of Refrigeration 26 (8) (2003) 918–925.

[50] A. Azzouz, J. Gossé, M. Duminil, Détermination expérimentale des pertes de froid

occasionnées par l’ouverture d’une porte de chamber froide industrielle, Revue Inter-
nationale du Froid 16 (1) (1993) 57––66.

[51] M. Bonaventure, A. Benchikh Le Hocine, S. Croquer, K. Huchtemann, S. Poncet,

Heat and mass transfer in a loaded truck trailer equipped with eutectic plates: a
comparative numerical study, in: Proc. of the 6th IIR Conference on Sustainability
and the Cold Chain (ICCC2020), Nantes, 2020, pp. 494–499.

[52] J. Emswiler, The neutral zone in ventilation, in: Annual Meeting of the American
Society of Heating and Ventilation Engineers, Buffalo, 1926, pp. 59–74.

[53] W. Brown, K. Solvason, Natural convection through rectangular openings in parti-

tions—1: Vertical partitions, International Journal of Heat and Mass Transfer 5 (9)
(1962) 859–868.

[54] W. Tamm, Kalterveluste durch Kuhlraumoffnungen, Kältetechnik-Klimatisierung 18

(1966) 142–144.

[55] C. Fritzsche, W. Lilienblum, Neue Messungen zur Bestimmung der Kälteverluste an

Kühlraumtüren, Kältetechnik-Klimatisierung 20 (1968) 279–286.

112
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 [56] W. Gosney, H. Olama, Heat and enthalpy gains through cold rooms doorways, in:
Proc. of Institute of Refrigeration, Environmental Science and Technology, The Poly-
technic of the South Bank, London, 1976.

[57] Q. Pham, D. Oliver, Infiltration of air into cold stores, in: Proc. of the 16th Interna-
tional Congress of Refrigeration, Vol. 4, West-Lafayette, 1983, pp. 67–72.

[58] W. Hendrix, D. Henderson, H. Jackson, Infiltration heat gains through cold storage
room doorways, ASHRAE Transactions 95 (CONF-890609) (1989).

[59] P. Chen, D. Cleland, S. Lovatt, M. Bassett, An empirical model for predicting air
infiltration into refrigerated stores through doors, International Journal of Refriger-
ation 25 (6) (2002) 799–812.

[60] T. Lafaye de Micheaux, M. Ducoulombier, J. Moureh, V. Sartre, J. Bonjour, Experi-

mental and numerical investigation of the infiltration heat load during the opening of
a refrigerated truck body, International Journal of Refrigeration 54 (2015) 170–189.

[61] C. Tso, S. Yu, H. Poh, P. Jolly, Experimental study on the heat and mass transfer
characteristics in a refrigerated truck, International Journal of Refrigeration 25 (3)
(2002) 340–350.

[62] F. Clavier, V. Sartre, J. Bonjour, Infiltration heat load through the doorway of a
refrigerated truck protected with an air curtain, in: Proceedings of the 23th Interna-
tional Congress of Refrigeration, Prague, 2011.

[63] J. Moureh, S. Tapsoba, E. Derens, D. Flick, Air velocity characteristics within vented
pallets loaded in a refrigerated vehicle with and without air ducts, International
Journal of Refrigeration 32 (2) (2009) 220–234.

[64] J. Moureh, M. Tapsoba, D. Flick, Airflow in a slot-ventilated enclosure partially

filled with porous boxes: Part II – measurements and simulations within porous
boxes, Computers & Fluids 38 (2) (2009) 206–220.

[65] A. Rai, J. Sun, S. Tassou, Numerical investigation of the protective mechanisms of

air curtain in a refrigerated truck during door openings, Energy Procedia 161 (2019)
216–223.

[66] A. Rai, J. Sun, S. Tassou, Three-dimensional investigation on the positioning of air

curtain on its effectiveness in refrigerated vehicles used for food distribution, Energy
Procedia 161 (2019) 224–231.

[67] S. Getahun, A. Ambaw, M. Delele, C. Meyer, U. Opara, Analysis of airflow and

heat transfer inside fruit packed refrigerated shipping container: Part I – model
development and validation, Journal of Food Engineering 203 (2017) 58–68.

[68] P. Torres Jara, J. Aguirre Rivera, C. Buenano Merino, V. E., G. Abad Farfán, Ther-
mal behavior of a refrigerated vehicle: Process simulation, International Journal of
Refrigeration 100 (2019) 124–130.

[69] D. Gray, A. Giorgini, The validity of the Boussinesq approximation for liquids and
gases, International Journal of Heat and Mass Transfer 19 (5) (1976) 545–551.

113
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 [70] F. Menter, Two-equation eddy-viscosity turbulence models for engineering applica-
tions, AIAA Journal 32 (8) (1994) 1598–1605.

[71] M. Hoang, P. Verboven, J. De Baerdemaeker, B. Nicolai, Analysis of the air flow

in a cold store by means of computational fluid dynamics, International Journal of
Refrigeration 23 (2000) 127–140.

[72] A. Foster, M. Swain, R. Barrett, P. D’Agaro, S. James, Effectiveness and optimum jet
velocity for a plane jet air curtain used to restrict cold room infiltration, International
Journal of Refrigeration 29 (5) (2006) 692–699.

[73] F. Hayes, W. Stoecker, Design data for air curtains, ASHRAE Transactions 75 (2)
(1969) 168–180.

[74] J. Paquette, S. Mercier, B. Marcos, S. Morasse, Modeling the thermal performance

of a multilayer box for the transportation of perishable food, Food and Bioproducts
Processing 105 (2017) 77–85.

[75] V. Nascimento Jr, F. Loyola, C. Hermes, A study of frost build-up on parallel plate
channels, Experimental Thermal and Fluid Science 60 (2015) 328–336.

[76] C.-H. Cheng, K.-H. Wu, Observations of early-stage frost formation on a cold plate
in atmospheric air flow, Journal of Heat Transfer 125 (1) (2003) 95–102.

[77] C. Hermes, R. Piucco, J. Barbosa Jr, C. Melo, A study of frost growth and den-
sification on flat surfaces, Experimental Thermal and Fluid Science 33 (2) (2009)
371–379.

[78] W. Wang, Q. Guo, W. Lu, Y. Feng, W. Na, A generalized simple model for predicting
frost growth on cold flat plate, International Journal of Refrigeration 35 (2) (2012)
475–486.

[79] C. Wen, Y. Yu, A generalized method for predicting the minimum fluidization veloc-
ity, AIChE Journal 12 (3) (1966) 610–612.

[80] W. Ranz, W. Marshall, Evaporation from droplets, Chemical Engineering Progress

48 (3) (1952) 141–146.

[81] J. Jeong, S. Poncet, B. Michel, J. Bonjour, Eulerian-eulerian multiphase frost model

based on phase change driving force, in: 7th IIR International Conference on Sus-
tainability and the Cold Chain, Newcastle, 2022.

[82] B. Kamkari, H. Shokouhmand, Experimental investigation of phase change material

melting in rectangular enclosures with horizontal partial fins, International Journal
of Heat and Mass Transfer 78 (2014) 839–851.

[83] M. Rahal, Effet de revêtements glaciophobes bioinspirés sur la formation du givre

sur une plaque froide, Master’s thesis, Université de Sherbrooke (2022).

[84] G. Li, Y. Hwang, R. Radermacher, H. Chun, Review of cold storage materials for
subzero applications, Energy 51 (2013) 1–17.

114
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 [85] J. Jeong, A. Benchikh Le Hocine, S. Croquer, S. Poncet, B. Michel, J. Bonjour, Nu-
merical analysis of the thermoaeraulic behavior of air during the opening of the door
of a refrigerated truck trailer equipped with cold plates, Applied Thermal Engineering
206 (2022) 118057.

[86] M. Ben Taher, T. Kousksou, Y. Zeraouli, M. Ahachad, M. Mahdaoui, Thermal per-

formance investigation of door opening and closing processes in a refrigerated truck
equipped with different phase change materials, Journal of Energy Storage 42 (2021)
103097.

115
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés
 FOLIO ADMINISTRATIF

THESE DE L’INSA LYON, MEMBRE DE L’UNIVERSITE DE LYON

NOM : JEONG DATE de SOUTENANCE : 24/04/2024

Prénoms : JIHYUK

TITRE : Modélisation CFD du transfert de chaleur et de masse dans une remorque de camion réfrigérée équipée
de plaques eutectiques

NATURE : Doctorat Numéro d'ordre : 2024ISAL0034

Ecole doctorale : Mécanique, Énergétique, Génie Civil, Acoustique

Spécialité :Thermique – Énergétique

RÉSUMÉ :
Le système de refroidissement par plaques eutectiques (matériau à changement de phase, MCP) est une possible alternative
aux systèmes conventionnels de réfrigération alimentés par des combustibles fossiles pour le transport de produits alimentaires
surgelés ou réfrigérés dans des remorques de camions. Des modèles numériques ont été développés pour évaluer sa
faisabilité et ont été validés avec succès par rapport à des résultats numériques et expérimentaux issus de la littérature.

Initialement, un modèle d'infiltration d’air dans la remorque a été développé utilisant un modèle de turbulence − SST pour
prédire le comportement thermoaéraulique de l'air pendant les périodes d'ouverture des portes. Différentes configurations des
plaques eutectiques et des ventilateurs ont été analysées. Sans cargaison, les plaques disposées en série le long du plafond
de la remorque ont montré un temps de renouvellement plus élevé que celles disposées en parallèle à l'arrière, en raison des
zones de recirculation. Cependant, une fois la cargaison introduite, les deux configurations offrent des performances similaires
car les zones de recirculation n'ont pas pu se former.

Par ailleurs, un modèle multiphasique granulaire eulérien-eulérien a été développé pour prédire la formation et la croissance du
givre sur les plaques eutectiques. Le modèle de turbulence − SST a été intégré pour étendre l'applicabilité du modèle de
givre à une gamme plus large de vitesses d’air. Un modèle de solidification et de fusion a également été implémenté et couplé
aux modèles précédents. Avec ce modèle combiné, les performances d'un système eutectique ont été étudiées pour des
conditions estivales typiques à Montréal, Canada. Environ 2.3 % du MCP change effectivement de phase au cours des 120
premières secondes. Globalement, le système eutectique offre une alternative viable au système conventionnel, bénéficiant
des mécanismes de prévention des infiltrations ou de dégivrage.

Mots-clés : Réfrigération, Matériaux à changement de phase, Transfert de chaleur et de masse, Dynamique des fluides
numérique.

Laboratoire (s) de recherche : CETHIL – INSA Lyon

LMFTEUS – Université de Sherbrooke (Cotutelle)

Directeur de thèse:
BONJOUR, Jocelyn, Professeur des Universités, INSA Lyon, Directeur de thèse
PONCET, Sébastien, Professeur, Université de Sherbrooke, Co-directeur de thèse

Président de jury : FOURNAISON, Laurence

Composition du jury :
SAFDARI SHADLOO, Mostafa, Maître de Conférences HDR, INSA Rouen, Rapporteur
GOSSELIN, Louis, Professeur, Université Laval (Québec), Rapporteur
FERTEL, Camille, Docteure, Tecnea Canada, Examinatrice
FOURNAISON, Laurence, Directrice de Recherche, INRAE, Examinatrice
MICHEL, Benoît, Maître de Conférences, INSA Lyon, Examinateur
BONJOUR, Jocelyn, Professeur des Universités, INSA Lyon, Directeur de thèse
PONCET, Sébastien, Professeur, Université de Sherbrooke, Co-directeur de thèse

Thèse accessible à l'adresse : [Link] © [J. Jeong], [2024], INSA Lyon, tous droits réservés