# A Hybrid Numerical Approach for Characterising Airflow and Temperature Distribution in a Ventilated Pallet of Heat-Generating Products: Application to Cheese

Dihia Aguenihanai, Denis Flick, Steven Duret, J. Moureh

► To cite this version:

Dihia Aguenihanai, Denis Flick, Steven Duret, J. Moureh. A Hybrid Numerical Approach for Characterising Airflow and Temperature Distribution in a Ventilated Pallet of Heat-Generating Products: Application to Cheese. Journal of Food Engineering, 2025, 387, pp.112323. 10.1016/j.jfoodeng.2024.112323. hal-04703396

HAL Id: hal-04703396

https://hal.inrae.fr/hal-04703396v1

Submitted on 20 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Hybrid Numerical Approach for Characterising Airflow and Temperature Distribution in a Ventilated Pallet of Heat-Generating Products: Application to Cheese

Dihia AGUENIHANAI*(a,b), Denis FLICK(c), Steven DURET(a), Jean MOUREH(a)

(a) Université Paris-Saclay, INRAE, FRISE, 92761 Antony, France

(b) CNIEL, 75009 Paris, France

(c) Université Paris-Saclay, INRAE, AgroParisTech, UMR SayFood, 91120 Palaiseau, France

*Corresponding author: dihia.aguenihanai@inrae.fr

## Abstract

Temperature control throughout the cold chain is of crucial importance in the preservation of the quality of cheese. As a result of cheese heat generation, both natural and forced convection need to be considered. This numerical study aimed to characterise the airflow and temperature fields within a ventilated pallet of heat-generating cheeses. An original computational fluid dynamics (CFD) hybrid approach was developed. This approach is based on a combination of a porous media approach for the contents of the boxes and a direct CFD approach for the outer cardboard walls, including vent size and position. The computational domain is limited to one pallet level. The simulations were conducted on a steady state for two upwind air velocities 0.31 m/s and 0.73 m/s and three generated heat fluxes 0.05 W, 0.15 W, and 0.3 W per product item (250 g). The model was validated by comparison with experimental results related to velocity and product temperature profiles obtained on a full-scale experimental set-up. The hybrid approach shows good accuracy while reducing the mesh size and the computational time in comparison with the direct CFD approach.

**Keywords:** Heat-generating products, mixed convection, porous medium, direct CFD, hybrid approach

## Nomenclature

<table><tr><td>C<sub>p</sub></td><td>Heat capacity at constant pressure, J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>D</td><td>Cheese diameter, m</td></tr><tr><td>F</td><td>Forchheimer coefficient, [-]</td></tr><tr><td>K</td><td>Permeability, m<sup>2</sup></td></tr><tr><td>k</td><td>Turbulent kinetic energy, m<sup>2</sup>/s<sup>2</sup></td></tr><tr><td>Ø<sub>z,in</sub></td><td>Vertical airflow rate, m<sup>3</sup>/s</td></tr><tr><td>p</td><td>Pressure, Pa</td></tr><tr><td>Q</td><td>Heat generation flux per product item (250 g of camembert), W</td></tr><tr><td>Q<sub>heat.tot</sub></td><td>Heat generation for one carton (30 camemberts), W</td></tr></table>
<table><tr><td>q<sub>heat.tot</sub></td><td>Volumetric heat generation for one carton, W/m<sup>3</sup></td></tr><tr><td>S</td><td>Momentum source term, kg/(m s)<sup>2</sup></td></tr><tr><td>S<sub>interface</sub></td><td>Interface surface between the porous zone and the air gap, m<sup>2</sup></td></tr><tr><td>T</td><td>Temperature, °C</td></tr><tr><td>u</td><td>Air velocity, m/s</td></tr><tr><td>u<sub>z</sub></td><td>Vertical component of the airflow velocity, m/s</td></tr><tr><td>V</td><td>Volume, m<sup>3</sup></td></tr><tr><td>β<sub>T</sub></td><td>Thermal expansion coefficient, K<sup>-1</sup></td></tr><tr><td>ν</td><td>Kinematic viscosity, m<sup>2</sup>/s</td></tr></table>

Abbreviations

<table><tr><td>CFD</td><td>Computational Fluid Dynamics</td></tr><tr><td>FV</td><td>Front vent</td></tr><tr><td>LDV</td><td>Laser Doppler Velocimetry</td></tr><tr><td>REV</td><td>Representative Elementary Volume</td></tr><tr><td>RMSE</td><td>Root Mean Square Error, °C</td></tr></table>

Dimensionless numbers

<table><tr><td>Gr</td><td>Grashof number, [-]</td></tr><tr><td>Pr</td><td>Prandtl number, [-]</td></tr><tr><td>Ra</td><td>Rayleigh number, [-]</td></tr><tr><td>Re</td><td>Reynolds number, [-]</td></tr><tr><td>Ri</td><td>Richardson number, [-]</td></tr></table>

Greek symbols

<table><tr><td>ε</td><td>Turbulence dissipation rate, m<sup>2</sup>/s<sup>3</sup></td></tr><tr><td>ε<sub>f</sub></td><td>Fluid porosity of stacked cheese: ε<sub>f</sub> = V<sub>f</sub>/V<sub>tot</sub>, [-]</td></tr><tr><td>λ</td><td>Thermal conductivity, W/(m K)</td></tr><tr><td>μ</td><td>Dynamic viscosity, kg m<sup>-1</sup>/s</td></tr><tr><td>ρ</td><td>Density, kg/m<sup>3</sup></td></tr><tr><td>ω</td><td>Specific turbulent dissipation rate, s<sup>-1</sup></td></tr></table>

Subscripts

<table><tr><td>air.in</td><td>Upwind airflow</td></tr><tr><td>d</td><td>driving</td></tr><tr><td>exp</td><td>Experimental</td></tr><tr><td>f</td><td>Fluid</td></tr><tr><td>max</td><td>Maximum</td></tr><tr><td>num</td><td>Numerical</td></tr><tr><td>out</td><td>Outlet</td></tr></table>
<table><tr><td>p</td><td>Porous media</td></tr><tr><td>t</td><td>Turbulent</td></tr><tr><td>tot</td><td>Total</td></tr></table>

# 1. Introduction

The temperature of cheese within pallets must be controlled and maintained below a recommended value throughout the cold chain to preserve quality and reduce losses (Pham et al., 2019). In the forced convection regime at the beginning of the cold chain (precooling), temperature heterogeneity can be encountered within stacks such as pallets, depending on several parameters such as the airflow rate (Wang et al., 2020), stacking pattern (Moureh et al., 2009b), ventilated box configuration including vent hole design and the total area of the openings (Berry et al., 2021). Heat generation induced by the respiration activity of the products can enhance natural convection and affect the level and heterogeneity of the product temperature. It is important to take this aspect into account, particularly in the case of periods during which the air velocity is low, such as during storage, where forced convection can be dominated by natural convection (Pham et al., 2021), of polylined products (O'Sullivan et al., 2017), and during long precooling processes (Redding et al., 2016).

Soft cheese products such as “Camemberts” contain lactic ferments and continue their metabolism once the products have been conditioned. As a result of the metabolic activity (respiration) of the micro-organisms, they generate heat (Hélias et al., 2008; Pham, 2021). This heat generation due to product respiration depends on the product temperature (Delele et al., 2013; Han et al., 2015) and can vary from one product to another. Recently, Pham et al. (2021) studied the effect of the natural convection induced by the respiration heat of cheese on the temperature of products within one level of a pallet. The study showed that products under mixed convection regime are better cooled than those under a forced regime. Chourasia and Goswami (2007) obtained similar results and showed that recirculation due to natural convection helps remove heat from the products, leading to a lower product temperature. Thus, natural convection exerts a positive influence on the heat transfer within a product stack.

Several experimental studies have been conducted to better understand and evaluate the cooling behaviour of stacks of products (Duret et al., 2014; Praeger et al., 2020). These studies represent real conditions and consider the random arrangement of products within a cardboard box and the spaces between ventilated boxes which can affect the connection between successive holes under real ventilation conditions. Nevertheless, there are few such studies in the literature because of their cost and the required manual handling time (Zou et al., 2006a). Computational Fluid Dynamics (CFD) modelling is therefore used to avoid the complexity of implementing experiments or to reduce the number of experiments.

Two main CFD methods are the most widely used: the direct CFD and porous media approaches. Both numerical approaches have their advantages and drawbacks. The direct CFD approach explicitly
reproduces both ventilated box geometry and product shapes and their arrangements. As it explicitly solves the local heat transfer and fluid transport, it accurately predicts the characteristics of local airflow through vent holes and around products as well as the product temperature level and heterogeneity within the stack. However, it involves a large number of meshes and a high computational cost. Hence, it is limited to a small computational domain, such as one box (Wang et al., 2020), one pallet level (Han et al., 2018), and rarely one pallet (Wu et al., 2019).

To manage computational resource issues, a porous medium approach has been used to model large domains such as cold rooms (Hoang et al., 2015) and transport vehicles (Moureh et al., 2009a). This method simplifies the geometry and models the flow roughly with a reduced computational time and an acceptable accuracy. According to Hoang et al. (2015), results obtained from modelling four pallets within a cold room using porous media, the relative error between numerical and experimental values of Nusselt number and velocity are 15% and 44%, respectively. However, this approach lacks sensitivity in modelling local airflow behaviour, issuing vent holes and temperature characteristics within boxes or stacks (Zou et al., 2006b). It involves integrating the conservation equations in the porous region. The airflow in the porous medium is characterized by its superficial velocity determined through local volume averaging. In addition, it requires the determination of additional parameters related to the pressure drop coefficients (often by using the Darcy-Forchheimer equation) and the heat transfer coefficient in the case of Local Thermal Non-Equilibrium (LTNE) modelling (Zhao et al., 2016; Alvarez et al., 2003; Hoang et al., 2015). Furthermore, (Pham et al., 2021) and (Ambaw et al., 2017), showed that pressure losses are located at the vents of boxes for a stack of products. This highlights the importance of modelling the packaging details (vents position, size and shape) to better characterise local airflow behaviour, which in turn has an impact on temperature distribution within the stack.

To overcome the limitations of both the direct CFD (mainly related to calculation resources limitations for a large-scale modelling) and porous medium approaches (lack of sensitivity in modelling local airflow behaviour), this study aimed to develop an original hybrid approach applied to a single pallet level of heat-generating cheese. This new method combines the porous media method for the contents of the boxes (cheeses) and the direct CFD approach to explicitly model the ventilated cartons, including vent design, size, and positions. A previous numerical study was carried out by Pham et al. (2021) on the same computational domain (one level pallet of cheese generating heat) using direct CFD approach. It involved a computational time of 48 hours and a mesh size of $10^7$ cells. Both studies take into account forced convection and natural convection promoted by cheese heat generation.

The objectives of this study were to evaluate the ability of the hybrid approach in reducing calculation time while ensuring accuracy of the airflow and temperature prediction. The hybrid approach will help to gain a better insight into local airflow behaviour, notably for the interaction of small jets through vent holes with the porous stack of product, while reducing mesh size and computing time. The results obtained using this new approach were evaluated and compared to the experimental and
numerical results obtained by Pham et al. (2021). Different conditions of air velocity and heat generation were tested to evaluate the approach performances, discuss its limitations and identify future research.

# 2. Materials and methods

## 2.1 Experimental device and measurements

The experimental device consists of one cheese pallet with industrial dimensions (800 mm × 1200 mm × 1455 mm) placed within a cold room with 5100 mm × 3200 mm × 2860 mm dimensions subjected to ambient parameter (temperature and velocity) control. The temperature of the cold room was set at $T_{air,in} = 4$°C. The pallet was composed of nine levels with six cardboard boxes of product per level. A space of about 10 mm separated the boxes placed at the same level. Each cardboard box contained 30 cheeses superposed in three layers and arranged in two rows (see Fig. 2). Each product (D = 110 mm, H = 40 mm) weighed 250 g and generated heat flux $Q$. In order to ensure both experimental stability conditions and repeatability, cylindrical plaster products equipped with controllable heating resistances were used to represent cheese products (camembert) and their heat generation; more details can be found in (Pham et al., 2019).

The experimental measurements used in this study were obtained by Pham et al. (2021). The transversal air velocity was measured at the mid-height of the air gap using Laser Doppler Velocimetry (LDV). The temperature of the products was measured on the axis of the cheese at 45 mm from the top using T-Type thermocouples with an accuracy of ± 0.1 K (precision: ± 0.05 K) after individual calibration between 273 and 293 K (Pham et al., 2021) (see Fig. 1).

## 2.2 Numerical modelling

### 2.2.1 Computational domain

In this study, a three-dimensional (3D) CFD model was developed and applied at steady state. Based on the same computational domain as Pham et al. (2021), the interactions between the different levels of the pallet were not taken into account. The model is limited to one level (see Fig. 2). Therefore, the bottom vents of the cardboard box were not considered.

The space between the cardboard boxes was included in the 3D geometry. A previous steady-state study by Pham et al. (2021) was performed on the same computational domain (one pallet level) using a direct CFD approach. This study considered a cardboard box with the same shape and size, including its vents, as that used by Pham et al. (2021).

### 2.2.2 Hybrid approach versus direct CFD approach

The simulations applied in this study were conducted using a new method called the hybrid approach. This method aims to combine in the same model the two most commonly used modelling methods: direct CFD and porous media approaches. As shown in Fig. 3, this original approach involves modelling the internal part of the cheese package (i.e., cheese and surrounding air) as a porous zone,
while the cardboard box, including its vent sizes and positions, is explicitly modelled using the direct CFD approach. The modelling of the box vent details (design and positions) allows better predictions for local jet trajectories and preferential air paths. According to Pham et al. (2021) and Ambaw et al. (2017), pressure drop due to product stacking is negligible, whereas the pressure drop across vents is more significant. Also, the validation of this approach allows the modelling of an entire pallet or a load of pallets, including their mutual interactions (which is hardly possible with the direct CFD method).

In addition, it is important to underline that to accurately capture the behaviour of a system, the representative elementary volume (REV) must be significantly larger than the pore scale or the fruit diameter (Defraeye et al., 2022) but still much smaller than the overall dimensions of the macroscopic flow domain (Nield and Bejan, 2017). According to Zhao et al. (2016), the continuous-medium assumption becomes questionable, especially when the ratio of package-to-produce diameter is less than 10 which often occurs in the case of individual packages of cheese. In these cases, the heterogeneity of the local airflow pattern within packages exerts a major impact on the transport phenomena (Ferrua and Singh, 2008). More the product diameter is high, more the deviation with real configuration become potentially high.

The calculated fluid porosity of one box of products is 34.5% (porosity calculation details were added to the supplementary material section). In order to simplify the geometry and reduce computational time, the thickness of the walls of the cardboard box was not considered in the modelling domain. Nevertheless, its thermal resistance effect on the convective heat transfer between the product and the ambience was taken into account.

The direct CFD method used by Pham et al. (2021) is based on detailed modelling of box and product geometry of one pallet level. It, therefore, provides results with a good accuracy in comparison with experimental results but requires considerable computing time (48 hours for the modelling of a single pallet level).

Product heat generation (Q per cheese item) and induced natural convection are taken into account in both direct CFD and hybrid approaches. For the direct CFD method, the heat generation is considered uniform in each camembert. For the hybrid method, heat generation is considered as uniform within the porous zone domain.

To validate the model and compare the two approaches, hybrid and direct CFD, the results obtained using this new approach related to the airflow velocity and product temperature were compared to the experimental and numerical results obtained by Pham et al. (2021). Fig. 1 summarises the positions at which temperature and velocity measurements were performed for this comparison:

* The transversal air velocity profiles along the lines $L_i \in [1, 3]$ for the three cartons (or boxes).

* The temperature of the products along the line $L_T$.
### 2.2.3 Modelling conditions and thermophysical properties of the materials

The modelling conditions of this numerical study were identical to those used by Pham et al. (2021). The simulations were performed for two air upwind velocities $u_{air,in}$: 0.31 m/s and 0.73 m/s and three heat generation per product item Q: 0.05 W, 0.15 W and 0.3 W. According to an internal calorimetry study of the heat generation of 250 g of camembert products, the measured heat generation of this type of cheese is estimated between 0.1W and 0.15W (unpublished report (Delahaye et al., 2019)). This result is in agreement with the calculations using the data of Hélias et al. (2007) and Gaucel et al. (2012), where the calculated heat-generated flux of 250 g of cheese is estimated to be 0.1 W.

Depending on the pallet position and the cold chain facilities, such as transport and storage stages, the airflow velocity can vary from 0.1 m/s (or even less) to 1 m/s (Hoang et al., 2015; Moureh et al., 2009b). Thus, the Richardson number during the storage can vary between 0.14 and 13.61 ($Q = 0.1$ W and $0.1$ m/s $\le u_{air,in} \le 1$ m/s). However, the upwind air velocities in the experimental cold room cannot be lower than 0.31 m/s. In order to: 1) obtain a Richardson number closer to the range of industrial values and 2) take into consideration the thermal runway conditions of cheese products during the cold chain (Leclercq-Perlat et al., 2012), the choice was made to adjust the heat flux. Three heat fluxes were chosen: 0.05W, 0.15W and 0.3W ($0.13 \le Ri \le 4.25$).

Table 1 summarises the dimensionless numbers $Re$, $Re_p$ and $Ri$ for the different conditions.

<table><caption>Table 1: Dimensionless numbers of the model.</caption><thead><tr><th>u<sub>air,in</sub> (m/s)</th><th>Re (-)</th><th>Re<sub>p</sub> (-)</th><th>Q (W)</th><th>Ri (-)</th></tr></thead><tbody><tr><td rowspan="3">0.31</td><td rowspan="3">2172</td><td rowspan="3">6296</td><td>0.05</td><td>0.71</td></tr><tr><td>0.15</td><td>2.12</td></tr><tr><td>0.3</td><td>4.25</td></tr><tr><td rowspan="3">0.73</td><td rowspan="3">5115</td><td rowspan="3">14825</td><td>0.05</td><td>0.13</td></tr><tr><td>0.15</td><td>0.38</td></tr><tr><td>0.3</td><td>0.77</td></tr></tbody></table>

$$Re_p = \frac{u_{air,in} D}{\nu \epsilon_f} \qquad (1)$$

$$Re = \frac{u_{air,in} D}{\nu} \qquad (2)$$

$$Gr = \frac{g\beta_T QD^2}{\lambda\nu^2} \qquad (3)$$
$$ Ri = \frac{Gr}{Re^2} = \frac{g\beta_T Q}{\lambda u_{air.in}^2} \quad (4) $$

where: $v = 15.7 \times 10^{-6} m^2/s$; $\epsilon_f = 0.345$; $\lambda = 0.026 W/(m K)$; $\beta_T = 1/T_{air.in}=0.00336 K^{-1}$; the diameter D is considered equal to the diameter of the camemberts: $D = 0.11 m$.

As shown in Table 1, the Reynolds number related to porous medium ($Re_p$) is higher than 300. Therefore, the airflow regime in the camembert stack for both air upwind velocities is turbulent (Delele et al., 2013; Eisfeld and Schnitzlein, 2001; Nield and Bejan, 2017). Furthermore, according to Bejan (2013), the critical Reynolds number for round jet flows is about 30. In this study, the Reynolds number through round vents of the box with a minimum diameter of 15 mm (Pham et al., 2019) is $Re_{vents} = 296$ ($D_{min.vents} = 15 mm$; $u_{air.in} = 0.31 m/s$).

Mixed convection is considered when $Ri \approx 1$. The airflow is dominated by forced convection at a low Richardson number $Ri < 0.1$ and by natural convection at a high Richardson number $Ri > 10$.

The thermal conductivity of the plaster and of the corrugated cardboard are 0.35 W/(m K) (data from the manufacturer) and 0.064 W/(m K) (data from Ho et al. (2010)), respectively. In addition, it's important to mention that the thermal conductivities of plaster and cheese are close: between 0.32 and 0.38 W.m⁻¹K⁻¹ for cheese depending on the type of cheese (Iezzi et al., 2011). This leads to similar temperature levels at steady state between cheese and plaster.

### 2.2.4 Assumptions

Some simplifications and assumptions were made to reduce the computational cost while providing a good system representation.

* The airflow was considered incompressible and turbulent ($Re_p > 300$) (Delele et al., 2013; Nield and Bejan, 2017).

* Natural convection was taken into account using the Boussinesq approximation.

* The respiration heat flux per product item (Q) is considered independent of the temperature. It is used to identify the volumetric source term $q_{heat.tot} = Q_{heat.tot}/V_p$, where $Q_{heat.tot} = 30 Q$.

* The heat source term related to cheese respiration is assumed uniform on the porous zone.

* With the exception of the density variation in the buoyancy term assumed by the Boussinesq approximation, the thermophysical properties of both phases of the porous media are assumed to be constant and independent of the temperature.

* As the flow rate was high ($Re \gg 1$), the airflow resistance characteristics of the porous zone were established using Darcy-Forchheimer law (Verboven et al., 2006).

* Local Thermal Equilibrium (LTE) between the two phases of the porous media is assumed. It should be borne in mind that LTE is not strictly verified in the case of products generating heat (Verboven et al., 2006). Nevertheless, Delele et al. (2013) applied this assumption to grapes with respiration activity estimated to 49 W/m³ and 4.35 W/m³ for a product temperature of 21°C
218

and -0.5°C, respectively. Numerical results for the product temperature showed good accuracy

219

compared to experimental data, with a relative error of 17.5%.

220

## 2.2.5 Governing equations

221

The direct CFD approach solves locally the conservation equation of momentum (Reynolds time-Averaged Navier-Stokes) and energy for the air outside and inside the boxes. The interaction between air and product is ensured by the boundary conditions at the surface of the product items. With the porous media approach, air outside the boxes is treated in the same way, but inside the boxes the local volume-averaged conservation equation of momentum is solved. The interaction between air and product is taken into account by the Darcy-Forchheimer terms whose coefficients have to be determined specifically. The hybrid model includes heat generation of products and the promoted natural convection and assumes the local thermal equilibrium for the heat transfer resolution. The steady-state heat transfer and incompressible flow were solved using Reynolds Averaged Navier-Stokes (RANS) equations combined with a turbulence model.

222

223

224

225

226

227

228

229

230

231

* RANS Mass conservation equation

$$ \frac{\partial \bar{u}_i}{\partial x_i} = 0 \qquad (5) $$

232

* RANS momentum conservation equation

233

After the application of the Boussinesq approximation, the conservation equation of the momentum is given by:

234

$$ \rho_0 \frac{\partial(\bar{u}_i \bar{u}_j)}{\partial x_j} = -\frac{\partial \bar{p}_d}{\partial x_i} + \frac{\partial}{\partial x_j} \left[ \mu_f \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i} \right) \right] - \frac{\partial}{\partial x_j} \left( \rho_0 \overline{u_i' u_j'} \right) - \rho_0 \beta_T (T_f - T_0) g_i + \bar{S}_{p,i} \qquad (6) $$

235

where:

236

$(T_0, \rho_0)$ are the fluid temperature and the density at the reference state (air inlet).

237

$\bar{S}_{p,i}$ is a momentum source term that characterises the flow resistance in the porous medium. For a homogenous porous media, it is defined by the following Darcy-Forchheimer law:

238

$$ \bar{S}_{p,i} = -\frac{\mu_f}{K} \bar{u}_i - \rho_0 \frac{F}{\sqrt{K}} |\bar{u}| \bar{u}_i \qquad (7) $$

239

Viscous and inertial pressure drop coefficients are also used instead of K and F:

240

$C_1=1/K$: Viscous pressure drop coefficient [1/m²]

241

$C_2=F/K^{1/2}$: Inertial pressure drop coefficient [1/m]

242

243

The Reynolds tensor $\overline{u_i' u_j'}$ can be approximated as below:

$$ \overline{u_i' u_j'} = -\frac{\mu_t}{\rho_0} \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i} \right) + \frac{2}{3} k \delta_{ij} \qquad (8) $$
244

• RANS energy conservation equation

$$ \rho_0 C_{p_f} \bar{u}_1 \frac{\partial \bar{T}}{\partial x_i} = \frac{\partial}{\partial x_i} \left[ \lambda_{\text{eff}} \left( \frac{\partial \bar{T}}{\partial x_j} \right) \right] - \rho_0 C_{p_f} \frac{\partial}{\partial x_i} (\overline{u_1' \bar{T}'}') + q_{\text{heat-tot}} \quad (9) $$

245

$$ \lambda_{\text{eff}} = (1 - \varepsilon_f)\lambda_s + \varepsilon_f\lambda_f $$

246

The turbulent heat flux (u_1'T'), can be expressed as follows:

$$ \overline{u_1'T'} = - \frac{\mu_t}{\rho_0 Pr_t} \left( \frac{\partial \bar{T}}{\partial x_j} \right) \quad (10) $$

247

248

According to the k-ε turbulence model, the turbulent viscosity μ_t is related to turbulence kinetic energy (k) and its dissipation rate (ε).

$$ \mu_t = \rho_0 C_\mu \frac{k^2}{\varepsilon} \quad (11) $$

249

The k-ω turbulence model is given as:

$$ \mu_t = \rho_0 C_\mu \frac{k}{\omega} \quad (12) $$

250

### 2.2.6 Boundary conditions

251

252

253

254

255

256

As the numerical results of this present study were to be compared to the numerical and experimental results of Pham et al. (2021), the same boundary conditions were assumed related to an upwind air velocity u\_air.in and a temperature of 4°C, at the inlet and an atmospheric pressure, at the outlet. The wall boundaries of the channel were assumed as adiabatic no-slip walls (see Fig. 2). Concerning the boxes walls, they were considered as non-slip walls with coupled heat transfer with the surrounding environment (cheese/porous zone on one side and ventilated airflow on the other side).

257

### 2.2.7 Simulation parameters

258

259

260

Anisotropic Darcy Forchheimer coefficients (Table 2) were determined numerically using a direct CFD approach. More details on the determination of these coefficients have been added in the supplementary material section.

261

<table><caption>Table 2: Isotropic Darcy-Forchheimer coefficient.</caption><thead><tr><th>Directions</th><th>C<sub>1</sub> [1/m<sup>2</sup>]</th><th>C<sub>2</sub> [1/m]</th></tr></thead><tbody><tr><td>x</td><td>1.65 &times; 10<sup>7</sup></td><td>5.35 &times; 10<sup>2</sup></td></tr><tr><td>y</td><td>5.01 &times; 10<sup>6</sup></td><td>1.68 &times; 10<sup>2</sup></td></tr><tr><td>z</td><td>1.75 &times; 10<sup>5</sup></td><td>5.26 &times; 10<sup>1</sup></td></tr></tbody></table>

262

### 2.2.8 Simulation set-up

263

264

The simulations were performed using Ansys Fluent 2021 CFD software on a 64-bit Windows 10 computer with processor Intel® Xeon® W-2133 CPU @ 3.60GHz and 256 GB installed RAM with
12 cores. The simulations were performed at steady state. The chosen solving algorithm is “Coupled”; the upwind 2nd order scheme was used for the momentum and energy discretisation, and the upwind 1st order was adopted for the turbulent dissipation and kinetic energy.

In this study, different mesh sizes were tested: 9.63× 10³, 4.90 × 10⁴, 8.29 × 10⁵, 1.60 × 10⁶, 1.92 × 10⁶, 2.49 × 10⁶ and 7.35 × 10⁶ cells. The mesh sensitivity simulations were performed for $u_{air,in}=0.73$ m/s and $Q = 0.3$ W. As a compromise between computational time and the accuracy of the results, the mesh size adopted in this study was 8.29 × 10⁵. The calculation time for the chosen mesh was about 4 hours, whereas when the direct CFD approach was used, the calculation time was about 48 hours for the same calculation domain (10⁷ cells) (Pham et al., 2021).

In addition, four turbulence models were tested (standard k-ε, RNG k-ε, realizable k-ε and SST k-ω) for two configurations (0.31 m/s & 0.3W and 0.73 m/s & 0.3W). Standard k-ε model was chosen as it offers a better compromise for both cases.

More details about mesh and turbulence model sensitivity have been added in supplementary material section.

The numerical results of this study were compared to the experimental data of Pham et al. (2021) using RMSE (Root Mean Square Error) calculated as follows:

$$ \text{RMSE} = \sqrt{\frac{1}{N} \sum_{n=1}^{N} (T_{\text{p.num,n}} - T_{\text{p.exp,n}})^2} \quad \text{for } N = 6 \quad (13) $$

$T_{\text{p.num,n}}$ is the temperature taken at a depth of 45 mm within the porous medium at the same position considered in experiments (at different positions of camemberts $n \in [1, 6]$) (see Fig. 1).

# 3. Results and discussion

## 3.1.1 Airflow distribution

Fig. 4 compares the dimensionless transversal velocities obtained using the hybrid approach ($u_x^* = u_x/u_{air,in}$) with the experimental measurements and direct CFD results obtained by Pham et al. (Pham et al., 2021). The transversal velocity is obtained along the lines $L_1$, $L_2$ and $L_3$ located at the mid-height of the air gap of the boxes $B_1$, $B_2$ and $B_3$, respectively (see Fig. 1), for the air upwind velocity, 0.31 m/s and both heat fluxes 0.05 W and 0.3 W.

According to Fig. 4, the hybrid approach, similarly to the direct CFD approach results, shows reasonable agreement with the experimental results. In addition, the differences between the two models could be explained by the fact that the porous medium induces a global homogenisation of the cheese domain, implying less peak velocities in the air headspace than the direct CFD approach due to the local jet diffusion in the porous medium.

The peaks of the velocity profiles in Fig. 4 reflect the presence of the front vents $FV_1$, $FV_2$ and $FV_3$ located in the headspace at the top of the box, which behave as air jets supplying the air layer above
297

the cheese with relatively high local velocities. For both heat fluxes, the velocity peaks in the air layer of the upstream ventilated box B₁ are of almost the same order as the inlet velocity. In the principal vent FV₂, $u_x^* \approx 0.9$ for 0.05 W and $u_x^* \approx 0.8$ for 0.3 W. Based on the analysis of the inflow and outflow rates through the front and side vents (Fig. 5b), the airflow enters through the front and side vents across the symmetry plane resulting from the 10 mm space between the boxes. When comparing velocity profiles along L₁, L₂ and L₃ (Fig. 4) and according to the distribution of velocity magnitude at mid-height of the headspace as shown in Fig. 5a, the air velocity decreases from upstream to downstream through the boxes. This decrease can be explained by the deflection of the main flow towards the lateral side and its exit through the spaces between boxes B₁, B₂ and B₃, as the total outflow from B₁ is higher than the total inflow into B₂, as well as from B₂ into B₃ (Fig. 5b). This gives rise to strong pathways between the front and side vents where an important part of the main upstream flow exits the pallet domain through side vents, which in turn reduces the peak velocities of the main flow. The strongest short circuit is observed for FV₁ due to its close proximity to the side vent located in the same box B₁. On the contrary, the central position of FV₃ allows the corresponding jet to maintain longitudinal flow without deviation and thus allows better internal ventilation of the whole pallet.

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

These strong short-circuits drive an important part of flow outside the pallet and directly affect the ventilation of the downstream part of the pallet where high temperature levels are expected.

313

314

In the upstream part of the pallet related to boxes B₁ and B₂, where forced convection predominates, quasi-similar numerical and experimental velocity profiles are obtained regardless of the intensity of the heat flux. The hybrid approach gives more accurate results than the direct CFD approach. This concerns velocity profiles, peak positions and peak levels (see Fig. 4).

315

316

317

318

In the downstream part of the pallet, in Box B₃, lower velocities were obtained due to the deflection mechanism discussed above. More accurate results were also obtained for the hybrid approach in comparison with the direct CFD approach for both heat fluxes. Whereas there was good agreement between the hybrid approach and experimental results for the low heat flux (Q = 0.05 W), there were some discrepancies for the high heat flux (Q = 0.3 W) which gave rise to reverse flow with low and negative velocities mainly induced by secondary natural airflow circulation (Pham et al., 2021). Obviously, both numerical approaches failed to predict local reverse flows associated with negative velocities. This could be explained by both the complexity of a flow highly dominated by thermal instability induced by dynamic interactions between vertical thermal plumes and the horizontal main flow occurring in relatively stagnant areas, which can be accentuated at higher heat fluxes (Q = 0.3 W) (Pham et al., 2021). However, the two-equation turbulence models (standard k-ε model) based on the eddy viscosity tend to increase the numerical diffusion to ensure their numerical stability and thus are unable to predict complex secondary flows.

319

320

321

322

323

324

325

326

327

328

329

330

331

Concerning overall differences between the experimental and numerical values, it is important to note that the roughness of the primary packaging of the cheese, the imperfections in the arrangement, and the presence of several resistance heating cables and thermocouples inside the boxes between the

332

333
cheeses influenced the experimental results. Also, the vent holes were not perfectly aligned between two successive boxes, which may have decreased the air passage area. These factors, which increase the pressure losses and limit the airflow, were not considered in the modelling, leading to an overestimation of the velocity in box B₃. Furthermore, the numerical error propagates in the model. The numerical error on L₁ will be reflected on L₂. As for the numerical velocity through L₃, besides its individual error, it will include the numerical error of L₁ and L₂.

### 3.1.2 Temperature distribution

The thermal validation of the model was performed using the product temperature profile along the line L<sub>T</sub> (see Fig. 1). The profiles of the temperature difference with the incoming air (ΔT = T<sub>p</sub> – T<sub>air.in</sub>) were investigated for two air upwind velocities, 0.31 m/s and 0.73 m/s, for three heat fluxes 0.05 W, 0.15 W, and 0.3 W. Fig. 6 shows the experimental data, the hybrid approach and the direct CFD approach results.

The apparent differences that can be observed between the direct CFD and hybrid approach plots (Fig. 6) are related to their geometrical and volumetric heat generation representation. For the direct CFD approach, the cylindrical shape of the products was modelled and the heat is generated on each cheese position. The heat transfer between the air and cheese is achieved directly at the outer surface. Thus, the product temperature is highest in the centre of each product. Meanwhile, the hybrid approach modelled the box content as a porous medium with a parallelepiped shape associated with uniform properties such as the uniform heat generation per total volume of each porous zone. It tended to uniformise the temperature plots along the line L<sub>T</sub>.

Table 3 presents the RMSE and the ΔT<sub>max.exp</sub> for each configuration. In the case of lower inlet velocity 0.31 m/s, the corresponding RMSE values are 0.30°C, 0.37°C, and 0.95°C for 0.05 W, 0.15 W and 0.3 W, respectively. The RMSE values were relatively low in comparison with the maximum temperature difference between the product and the air: RMSE/ΔT<sub>max.exp</sub> varies between 4% and 8%.

<table><caption>Table 3: Estimation of the RMSE and the maximum experimental temperature differential for the hybrid approach.</caption><thead><tr><th>u<sub>air.in</sub> (m/s)</th><th>Q (W)</th><th>RMSE (°C)</th><th>ΔT<sub>max.exp</sub> (° C)</th></tr></thead><tbody><tr><td rowspan="3">0.31</td><td>0.05</td><td>0.30</td><td>3.71</td></tr><tr><td>0.15</td><td>0.37</td><td>9.32</td></tr><tr><td>0.3</td><td>0.95</td><td>16.57</td></tr><tr><td rowspan="3">0.73</td><td>0.05</td><td>0.41</td><td>2.25</td></tr><tr><td>0.15</td><td>1.16</td><td>6.05</td></tr><tr><td>0.3</td><td>2.25</td><td>11.52</td></tr></tbody></table>
For the highest air velocity, 0.73 m/s, unlike the direct CFD approach, the temperature distribution obtained using the hybrid approach was underestimated compared to the experimental data and RMSE/ΔT<sub>max,exp</sub> was around 20%. This can be explained by the local thermal equilibrium hypothesis (LTE) assumed by the porous medium model, which supposes an infinite convective heat transfer coefficient between air and solid phases. Thus, the temperatures of these two phases are considered equal in a representative elementary volume (REV) (Nield and Bejan, 2017). However, the experimental results obtained by Pham et al. (2019) clearly indicate that for products generating heat, their temperature becomes higher than that of the surrounding air. The difference increased with heat generation and could reach about 2°C, which is comparable with the underestimation of the hybrid approach in the worst case. This suggests than an improvement could be obtained using a two-temperature porous medium approach instead of LTE. However, this implies estimating the local heat transfer coefficient between the air and solid phases, which is a complex task because it depends on local velocity magnitude, velocity orientation and turbulence (Alvarez and Flick, 1999).

As shown in Fig. 6, at a velocity of 0.73 m/s, temperature drops can be observed downstream from the inlet of each box. To explain this phenomenon, Fig. 7 depicts the pathlines of the airflow on three vertical planes, P₁, P₂ and P₃, passing by the centre of the front vents FV₁, FV₂ and FV₃, respectively. The temperature along the line L<sub>T</sub> is measured on the plane P₂ (see Fig. 7b). This figure shows a complex interaction between the flow in the air gap (z > 120 mm) and in the porous zone (air + solid). Air flows faster in the air gap and is colder than in the porous zone. This is explained by the fact that flow resistance and heat generation are located in the porous zone alone. Just after the front face of a box, there is an entrainment effect exerted by the jets generated downstream from the upper front vents (FV₁, FV₂ and FV₃). A depression appears (suction effect): as can be seen on Fig. 7c, just after the front face (x < 5 cm), the pressure is lower compared to that in the middle of the box. Warm air from the porous zone flows upwards toward the jets and is replaced by cold air coming from the air gap. This recirculation explains the temperature drops in Fig. 6. This effect (intensive mixing between the air gap and the porous zone) is perhaps over-predicted by the hybrid approach model so that overall, the temperature on line L<sub>T</sub> is underestimated.

Fig. 8 presents the inflow rates through the different front vents compared to the total inflow rate of B₁. It can be seen that the sum of inflow rates through the front vents is not 100% in box B₁ because there is also a lateral inflow (as shown in Fig. 5b). The upper front vents have the greatest effect, whereas the lower ones exert less effects because of obstacles behind them. This also explains the much lower temperature in the air gap (z > 120 mm).

It appears that the exchange between the air gap zone and the porous zone is of major importance for heat evacuation. Therefore, the airflow rate, related to the vertical velocity component, between the two zones was calculated by Eq. (14) (upward flow compensates overall downflow in a box). It decreases along the three boxes in the ventilation direction. For the case where u<sub>air,in</sub>= 0.73 m/s and Q = 0.3 W, it is 1.38 L/s, 0.99 L/s and 0.55 L/s from the first to the last box.
$$ \Phi_{z.in} = \frac{1}{2} \int |u_z| \, dS_{interface} \qquad (14) $$

It can be observed that the recirculation shown in plane P₂ has a limited spatial effect near the principal vent FV₂ (Fig. 7b) and is not observed on the P₁ and P₃ planes. Therefore, only a limited area of the first row of cheeses in the box is affected by the temperature drop, as can be seen in the temperature contours shown in Fig. 9 (horizontal plane at z = 75 mm, in the porous zone). This also means that despite this local temperature drop, the overall thermal behaviour indicates quasi-homogeneity of the temperature row by row through the three boxes. This finding agrees with the results obtained by Pham et al. (2021), which demonstrated that the temperature distribution in the same row of cheeses was almost homogeneous.

According to the results presented in Fig. 6 and the temperature contours in Fig. 9, the temperature distribution is related to the airflow rate, the product positions and the heat flow the products generate. The higher the air velocity, the more heat is extracted from the products. According to Section 3.1.1 above, air heats up, and its velocity decreases in the main airflow direction, leading to better cooling in the upstream than in the downstream box (Fig. 6). The heat flux generated by the cheeses also significantly impacts the product temperature. The higher the heat flux, the higher the product temperature. For example, at an air velocity of 0.31 m/s, the cheese in position C₆ (Fig. 1) is at 7.3°C and 21.9°C for the two heat fluxes, 0.05 W and 0.3 W respectively, for inlet air at Tₐir.in = 4°C. However, because of free convection which enhances transfer, the effect is not linear. Indeed, when the heat flux was six times higher, there was a three-fold rise in Tₚ - Tₐir.in. This aspect is of major importance, since heat generation by cheeses is influenced by the temperature of the cheeses (Leclercq-Perlat et al., 2012). Therefore, if the pallet temperature is not properly maintained at a target temperature, heat generation may increase, thereby increasing the temperature of the cheeses, leading to a thermal runaway (Zhang et al., 2018). A more detailed study is needed in order to establish a relationship between cheese heat generation, temperature distribution, and mixed convection in a pallet in order to reduce losses and ensure product quality which can be related to different factors such as cheese odour, colour and its crust (underrind) thickness and consistency (Leclercq-Perlat et al., 2015).

### 3.1.3 Hybrid approach advantages and limitations

The results of this study showed the benefits of the hybrid approach, but also highlighted some limitations that could be further studied. The major difference between hybrid and direct approach is related to dynamic interactions between local air jets issuing vent holes and the physical objects representing packed products. Small products randomly distributed behave as porous media while bigger physical objects associated to higher velocities will affect in a specific manner the air jet development leading to local obstructions, jet deviations and preferential airflow pathways within the ventilated package.
430

To better extend the use of the hybrid approach to the products, future numerical and experimental studies could investigate different extreme cases combining for example different specific ventilated boxes (centrally and edge-positioned vents), different packaging pattern (layered packed products with specific horizontal air paths and randomly packed products) and products of different size. Comparisons between experimental and predicted values obtained with hybrid approach with different velocities will help to assess and to validate the hybrid approach for other type of products while underlying its limitation.

431

432

433

434

435

436

437

Finally, hybrid approach can be applied on a larger scale (entire cold room, truck) upon validation and conduct more accurate multiscale simulation describing both airflow at the equipment scale and detailed heat transfer in the pallets.

438

439

440

## 4. Conclusion

441

442

443

444

445

446

447

An original hybrid numerical approach was developed and applied to one pallet level of heat-generating cheese under steady-state conditions. The results agreed acceptably with the direct CFD approach results and experimental data regarding air velocity and product temperature. This method reduced the computational time by a factor of 12 compared to the direct CFD approach. The hybrid approach accurately characterises the preferential flow paths and temperature heterogeneity within the pallet through the detailed modelling of the boxes (position and shape of the vent holes). It provides a clear insight into the effect of mixed convection on temperature distribution.

448

449

450

451

452

The model will be improved by incorporating the local convective heat transfer coefficient and assuming Local Thermal Non-Equilibrium. Further investigations using this approach will be conducted to study the heating and cooling kinetics of heat-generating products within an entire pallet under transient state conditions. Since the box vents are modelled, the method can be generalised for the study of the adjacent pallets effect on the distribution of airflow and temperature within the pallets.

453

## Acknowledgements

454

455

456

457

The authors acknowledge and thank the French Dairy Interbranch Organization (CNIEL) and the French National Association of Research and Technology (ANRT) for the technical and financial support that they have provided. The authors would like to thank the project coordinator Fanny Tenenhaus-Aziza for her support.

458

## References

459

460

461

https://doi.org/10.1016/S0017-9310(02)00546-X

462

463

464

Alvarez, G., Bournet, P.-E., Flick, D., 2003. Two-dimensional simulation of turbulent flow and transfer through stacked spheres. International Journal of Heat and Mass Transfer 46, 2459–2469.

Alvarez, G., Flick, D., 1999. Analysis of heterogeneous cooling of agricultural products inside bins: Part II: thermal study. Journal of Food Engineering 39, 239–245. https://doi.org/10.1016/S0260-8774(98)00166-6
465

466

467

Ambaw, A., Mukama, M., Opara, U.L., 2017. Analysis of the effects of package design on the rate and uniformity of cooling of stacked pomegranates: Numerical and experimental studies. Computers and Electronics in Agriculture 136, 13–24. https://doi.org/10.1016/j.compag.2017.02.015

468

Bejan, A., 2013. Convection Heat Transfer. John Wiley & Sons.

469

Berry, T.M., Defraeye, T., Wu, W., Sibiya, M.G., North, J., Cronje, P.J.R., 2021. Cooling of ambient-loaded citrus in refrigerated containers: What impacts do packaging and loading temperature have? *Biosystems Engineering* 201, 11–22. https://doi.org/10.1016/j.biosystemseng.2020.11.002

470

471

472

Chourasia, M.K., Goswami, T.K., 2007. CFD simulation of effects of operating parameters and product on heat transfer and moisture loss in the stack of bagged potatoes. *Journal of Food Engineering* 80, 947–960. https://doi.org/10.1016/j.jfoodeng.2006.07.015

473

474

475

Defraeye, T., Lukasse, L., Shrivastava, C., Verreydt, C., Schemminger, J., Cronjé, P., Berry, T., 2022. Is there a systematic hidden “hot spot” in refrigerated containers filled with fresh food in ventilated packaging? *Trends in Food Science & Technology* 129, 388–396. https://doi.org/10.1016/j.tifs.2022.09.005

476

477

478

479

Delahaye, A., Gahartian, J., Jouquin, C., Oignet, J., Ndoye, F.T., 2019. Dispositif de calorimétrie pour la mesure de flux de chaleur dégagés au cours de l’entreposage de fromages.

480

481

Delele, M.A., Ngcobo, M.E.K., Opara, U.L., Meyer, C.J., 2013. Investigating the Effects of Table Grape Package Components and Stacking on Airflow, Heat and Mass Transfer Using 3-D CFD Modelling. *Food Bioprocess Technol* 6, 2571–2585. https://doi.org/10.1007/s11947-012-0895-5

482

483

484

485

Duret, S., Hoang, H.-M., Flick, D., Laguerre, O., 2014. Experimental characterization of airflow, heat and mass transfer in a cold room filled with food products. *International Journal of Refrigeration* 46, 17–25. https://doi.org/10.1016/j.ijrefrig.2014.07.008

486

487

488

Eisfeld, B., Schnitzlein, K., 2001. The influence of confining walls on the pressure drop in packed beds. *Chemical Engineering Science* 56, 4321–4329. https://doi.org/10.1016/S0009-2509(00)00533-9

489

490

491

Ferrua, M.J., Singh, R.P., 2008. A nonintrusive flow measurement technique to validate the simulated laminar fluid flow in a packed container with vented walls. *International Journal of Refrigeration* 31, 242–255. https://doi.org/10.1016/j.ijrefrig.2007.05.013

492

493

494

Gaucel, S., Guillemin, H., Corrieu, G., 2012. A generalised model for cheese mass loss determination during ripening. *Journal of Food Engineering* 110, 109–116. https://doi.org/10.1016/j.jfoodeng.2011.12.002

495

496

497

Han, J.-W., Zhao, C.-J., Qian, J.-P., Luis, R.-G., Xiang, Z., 2018. Numerical modeling of forced-air cooling of palletized apple: Integral evaluation of cooling efficiency. *International Journal of Refrigeration* 11.

498

499

500

Han, J.-W., Zhao, C.-J., Yang, X.-T., Qian, J.-P., Fan, B.-L., 2015. Computational modeling of airflow and heat transfer in a vented box during cooling: Optimal package design. *Applied Thermal Engineering* 91, 883–893. https://doi.org/10.1016/j.applthermaleng.2015.08.060

501

502

503

Hélias, A., Mirade, P.-S., Corrieu, G., 2007. Modeling of Camembert-Type Cheese Mass Loss in a Ripening Chamber: Main Biological and Physical Phenomena. *Journal of Dairy Science* 90, 5324–5333. https://doi.org/10.3168/jds.2007-0272

504

505

506

Hélias, A., Trelea, I.C., Corrieu, G., 2008. Assessment of respiratory activity during surface-mould cheese ripening. *Journal of Food Engineering* 85, 632–638. https://doi.org/10.1016/j.jfoodeng.2007.09.001

507

508

509

Ho, S.H., Rahman, M.M., Sunol, A.K., 2010. Analysis of thermal response of a food self-heating system. *Applied Thermal Engineering* 30, 2109–2115. https://doi.org/10.1016/j.applthermaleng.2010.05.020

510

511

512

Hoang, H.-M., Duret, S., Flick, D., Laguerre, O., 2015. Preliminary study of airflow and heat transfer in a cold room filled with apple pallets: Comparison between two modelling approaches and experimental results. *Applied Thermal Engineering* 76, 367–381. https://doi.org/10.1016/j.applthermaleng.2014.11.012

513

514

515

516

Iezzi, R., Francolino, S., Mucchetti, G., 2011. Natural convective cooling of cheese: Predictive model and validation of heat exchange simulation. *Journal of Food Engineering* 106, 88–94. https://doi.org/10.1016/j.jfoodeng.2011.04.016

517

518

519
520

Leclercq-Perlat, M.-N., Sicard, M., Perrot, N., Trelea, I.C., Picque, D., Corrieu, G., 2015. Temperature and relative humidity influence the ripening descriptors of Camembert-type cheeses throughout ripening. *Journal of Dairy Science* 98, 1325–1335. https://doi.org/10.3168/jds.2014-8916

521

522

523

Leclercq-Perlat, M.-N., Sicard, M., Trelea, I.C., Picque, D., Corrieu, G., 2012. Temperature and relative humidity influence the microbial and physicochemical characteristics of Camembert-type cheese ripening. *Journal of Dairy Science* 95, 4666–4682. https://doi.org/10.3168/jds.2012-5368

524

525

526

527

Moureh, J., Tapsoba, M., Flick, D., 2009a. Airflow in a slot-ventilated enclosure partially filled with porous boxes: Part I – Measurements and simulations in the clear region. *Computers & Fluids* 38, 194–205. https://doi.org/10.1016/j.compfluid.2008.02.006

528

529

530

Moureh, J., Tapsoba, S., Derens, E., Flick, D., 2009b. Air velocity characteristics within vented pallets loaded in a refrigerated vehicle with and without air ducts. *International Journal of Refrigeration* 32, 220–234. https://doi.org/10.1016/j.ijrefrig.2008.06.006

531

532

533

Nield, D.A., Bejan, A., 2017. Convection in Porous Media. Springer International Publishing, Cham. https://doi.org/10.1007/978-3-319-49562-0

534

535

O'Sullivan, J.L., Ferrua, M.J., Love, R., Verboven, P., Nicolaï, B., East, A., 2017. Forced-air cooling of polylined horticultural produce: Optimal cooling conditions and package design. *Postharvest Biology and Technology* 126, 67–75. https://doi.org/10.1016/j.postharvbio.2016.11.019

536

537

538

Pham, A.T., 2021. Caractérisation aéraulique et thermique au sein d'un empilement de produits dégageant de la chaleur: application au cas des palettes de fromage. *Agriculture, alimentation, biologie, environnement et santé*.

539

540

541

Pham, A.T., Moureh, J., Belaidi, M., Flick, D., 2021. CFD modelling of a pallet of heat-generating product applied to a cheese product. *International Journal of Refrigeration* 128, 163–176. https://doi.org/10.1016/j.ijrefrig.2021.03.011

542

543

544

Pham, A.T., Moureh, J., Flick, D., 2019. Experimental characterization of heat transfer within a pallet of product generating heat. *Journal of Food Engineering* 247, 115–125. https://doi.org/10.1016/j.jfoodeng.2018.12.003

545

546

547

Praeger, U., Jedermann, R., Sellwig, M., Neuwald, D.A., Hartgenbusch, N., Borysov, M., Truppel, I., Scaar, H., Geyer, M., 2020. Airflow distribution in an apple storage room. *Journal of Food Engineering* 269, 109746. https://doi.org/10.1016/j.jfoodeng.2019.109746

548

549

550

Redding, G.P., Yang, A., Shim, Y.M., Olatunji, J., East, A., 2016. A review of the use and design of produce simulators for horticultural forced-air cooling studies. *Journal of Food Engineering* 190, 80–93. https://doi.org/10.1016/j.jfoodeng.2016.06.014

551

552

553

Verboven, P., Flick, D., Nicolaï, B.M., Alvarez, G., 2006. Modelling transport phenomena in refrigerated food bulks, packages and stacks: basics and advances. *International Journal of Refrigeration*, *Issue with Special Emphasis on Data and Models on Food Refrigeration* 29, 985–997. https://doi.org/10.1016/j.ijrefrig.2005.12.010

554

555

556

557

Wang, D., Lai, Y., Jia, B., Chen, R., Hui, X., 2020. The optimal design and energy consumption analysis of forced air pre-cooling packaging system. *Applied Thermal Engineering* 165, 114592. https://doi.org/10.1016/j.aplthermaleng.2019.114592

558

559

560

Wu, W., Beretta, C., Cronje, P., Hellweg, S., Defraeye, T., 2019. Environmental trade-offs in fresh-fruit cold chains by combining virtual cold chains with life cycle assessment. *Applied Energy* 254, 113586. https://doi.org/10.1016/j.apenergy.2019.113586

561

562

563

Zhang, K., Pu, Y.-Y., Sun, D.-W., 2018. Recent advances in quality preservation of postharvest mushrooms (Agaricus bisporus): A review. *Trends in Food Science & Technology* 78, 72–82. https://doi.org/10.1016/j.tifs.2018.05.012

564

565

566

Zhao, C.-J., Han, J.-W., Yang, X.-T., Qian, J.-P., Fan, B.-L., 2016. A review of computational fluid dynamics for forced-air cooling process. *Applied Energy* 168, 314–331. https://doi.org/10.1016/j.apenergy.2016.01.101

567

568

569

Zou, Q., Opara, L.U., McKibbin, R., 2006a. A CFD modeling system for airflow and heat transfer in ventilated packaging for fresh foods: I. Initial analysis and development of mathematical models. *Journal of Food Engineering* 77, 1037–1047. https://doi.org/10.1016/j.jfoodeng.2005.08.042

570

571

572

573

Zou, Q., Opara, L.U., McKibbin, R., 2006b. A CFD modeling system for airflow and heat transfer in ventilated packaging for fresh foods: II. Computational solution, software development, and

574
575

576

577

578

model testing. Journal of Food Engineering 77, 1048-1058.

https://doi.org/10.1016/j.jfoodeng.2005.08.043
<figure><img src="image_2.png" /><figcaption>Fig. 1. Illustration of the instrumented level of the pallet and its dimensions: (a) 3D view; (b) front view; and (c) lateral view.</figcaption></figure>

580

581

<figure><img src="image_3.png" /><figcaption>Fig. 2. Illustration of (a) the computational domain of the model; (b) the dimensions of the cardboard boxes; and (c) the dimensions of the cheese cylinders.</figcaption></figure>

582

583

584
<figure><img src="image_4.png" /><figcaption>Fig. 3. Comparison between (a) direct CFD (Pham et al., 2021), and (b) hybrid approach geometry representation.</figcaption></figure>
$Q = 0.05 \text{ W}$

<figure><img src="image_5.png" /></figure>

$Q = 0.3 \text{ W}$

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 4. Comparison of the air gap transversal velocity through the lines L₁, L₂ and L₃ for inlet velocity 0.31 m/s and two heat generation fluxes, 0.05 W and 0.3 W, with both experimental data and direct CFD approach results of Pham et al. (2021).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 5. Airflow distribution through the boxes for $u_{air,in}$= 0.31 m/s and $Q=0W$; (a) air velocity magnitude at $z = 132.5$ mm; (b) inflow and outflow rates compared to a total inflow of box B₁.</figcaption></figure>
$$u_{\text{air.in}} = 0.31 \text{ m/s}$$

$$u_{\text{air.in}} = 0.73 \text{ m/s}$$

$$Q = 0.05 \text{ W}$$

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

$$Q = 0.15 \text{ W}$$

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

$$Q = 0.3 \text{ W}$$

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 6. Comparison of the temperature difference distribution ($\Delta T = T_p - T_{\text{air.in}}$) obtained with the hybrid approach through the line $L_T$ for air inlet velocities 0.31 m/s and 0.73 m/s and three heat fluxes, 0.05 W, 0.15 W and 0.3 W, with both experimental data and direct CFD results of Pham et al. (2021). Note: Experimental data of Pham et al. (2021) (red cross); direct CFD results of Pham et al. (2021) (blue points) and hybrid approach results (black dashed line). The dashed grey, green and red vertical lines used delimit boxes $B_1$, $B_2$ and $B_3$, respectively.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 7. (a) Diagram showing the positions of the planes; (b) Airflow pathlines coloured by the temperature through three planes $P_1$, $P_2$ and $P_3$ for $u_{air,in} = 0.73$ m/s and $Q=0.3$ W; (c) pressure evolution through the line $L_T$ in the box $B_1$.</figcaption></figure>

<figure><img src="image_19.png" /></figure>

Boxes Inflow Face

<figure><img src="image_20.png" /><figcaption>Fig. 8. Airflow rates in the different front vents of the boxes compared to the total inflow rate of box $B_1$ $u_{air,in} = 0.73$ m/s $Q = 0.3$ W.</figcaption></figure>

634

635

636

637

638

639

640
<figure><img src="image_21.png" /><figcaption>Fig. 9. Temperature distribution for two airflow velocities, 0.31 m/s and 0.73 m/s, and three heat-generated fluxes, 0.05 W, 0.15 W and 0.3 W (horizontal plane at z = 75 mm).</figcaption></figure>

641

642

643