domain: sdv spi

text:
Tanathep Leungtongkum
email: tanathep.leungtongkum@inrae.fr

Onrawee Laguerre

Steven Duret

Denis Flick

How to choose a model to address practical issues encountered during food transport in an insulated box equipped with phase change material

Keywords: Insulated box, Phase change material, Modelling, Heat transfer, Airflow, Temperature prediction

Highlights

 Lumped, zonal and CFD models were developed for transport in insulated boxes with PCM.  The demonstration of how to use these models to deal with practical issues are shown.  The lumped model is suitable where the temperature heterogeneity is not a concern.  The zonal model, more complex, provides temperature evolution in different zones.  The CFD model, the most complex, provides temperature and air velocity fields.

Introduction

Insulated boxes equipped with a Phase Change Material (PCM) has an important role in food cold chain [START_REF] Nie | Review on phase change materials for cold thermal energy storage applications[END_REF][START_REF] Zhao | Application and research progress of cold storage technology in cold chain transportation and distribution[END_REF], particularly for the last mile delivery to consumers [START_REF] Robertson | Innovations in cold chain equipment for immunization supply chains[END_REF]. The advantages of the insulated box are its simplicity of use, the flexibility related to several box designs available, and the low cost. However, product waste caused by temperature abuse were often observed because of spatial and temporal temperature variations during food transport [START_REF] Laguerre | Experimental investigation and modelling in the food cold chain: Thermal and quality evolution[END_REF][START_REF] Liu | Fibroin treatment inhibits chilling injury of banana fruit via energy regulation[END_REF][START_REF] Mercier | Time-Temperature Management Along the Food Cold Chain: A Review of Recent Developments: Food preservation along the cold chain[END_REF].

The temperature evolution inside an insulated box equipped with PCM was investigated experimentally and numerically by several authors and summarized in Leungtongkum et al. [START_REF] Leungtongkum | Insulated box and refrigerated equipment with PCM for food preservation: State of the art[END_REF]. However, no studies provided a general solution on how to efficiently transport food products in an insulated box with PCM.

Modelling is an important tool to help the users to choose the suitable box design (volume and insulation) and PCM (type and amount) for transporting food products under certain ambient conditions and transport durations. Navaranjan et al. [START_REF] Navaranjan | Thermal insulation requirements and new cardboard packaging for chilled seafood exports[END_REF] proposed a model to predict fish spoilage based on thermal resistance of the box but it was a data-based model which could be applied only under specific ranges of operating parameters. This work focuses on thermal based models by taking into account the physical phenomena in an insulated box with PCM. Thus, they are applicable to a wider range of operating conditions.

There are different thermal modelling approaches described in the literature: 1D conduction models inside the product [START_REF] Laguerre | Modelling of fish refrigeration using flake ice[END_REF][START_REF] Kozak | Experimental and comprehensive theoretical study of cold storage packages containing PCM[END_REF][START_REF] Anand | Performance improvement of a PCM cold box by two bilayers configuration[END_REF][START_REF] Laguerre | How to predict product temperature changes during transport in an insulated box equipped with an ice pack: Experimental versus 1-D and 3-D modelling approaches[END_REF], lumped models considering the average product temperature and overall heat transfer resistances [START_REF] Davey | Calculating quantities of ice for cooling and maintenance of freshly harvested fish at sea[END_REF][START_REF] Zeng | Study on the effect of aluminium foil on packaging thermal insulation performance in cold chain logistics[END_REF], zonal models distinguishing warmer and colder regions [START_REF] East | Combining a hybrid genetic algorithm and a heat transfer model to optimise an insulated box for use in the transport of perishables[END_REF], and Computational Fluid Dynamic (CFD) models describing detailed temperature and velocity fields [START_REF] Laguerre | Modelling of fish refrigeration using flake ice[END_REF][START_REF] Laguerre | How to predict product temperature changes during transport in an insulated box equipped with an ice pack: Experimental versus 1-D and 3-D modelling approaches[END_REF][START_REF] Burgess | Optimisation of a portable phase-change material (PCM) storage system for emerging cold-chain delivery applications[END_REF][START_REF] Margeirsson | Numerical modelling of temperature fluctuations in superchilled fish loins packaged in expanded polystyrene and stored at dynamic temperature conditions[END_REF][START_REF] Du | Cooling performance of a thermal energy storage-based portable box for cold chain applications[END_REF][START_REF] Xiaofeng | Simulation and experimental investigation of a multitemperature insulation box with phase change materials for cold storage[END_REF][START_REF] Tiari | Numerical study of a latent heat thermal energy storage system enhanced by varying fin configurations[END_REF][START_REF] Ray | Comparative study of cooling performance for portable cold storage box using phase change medium[END_REF][START_REF] Purandare | Experimental and numerical study of insulation packages containing dry ice pellets[END_REF][START_REF] Fu | Study on the storage time of a cold box based on conduction-convection-radiation coupling[END_REF][START_REF] Calati | CFD analyses for the development of an innovative latent thermal energy storage for food transportation[END_REF][START_REF] Rahimi-Khoigani | Application of an improved latent heat storage system in the food packaging[END_REF].

The 1D models predict the product temperature at different depths from the surface based on only heat conduction equation. According to our experimental results of an insulated box with PCM, natural convection should be taken into account in air gaps [START_REF] Leungtongkum | Experimental investigation of airflow and heat transfer by natural convection in an insulated box with a Phase Change Material using a Particle Image Velocimetry technique[END_REF], thus these models are not adapted to our work.

The lumped models predict the effect of operating parameters on average temperature evolution under wide ranges of box insulation, PCM (type and amount) and product. In literature, the lumped models considered only external heat convection and heat conduction within box's wall [START_REF] Davey | Calculating quantities of ice for cooling and maintenance of freshly harvested fish at sea[END_REF][START_REF] Zeng | Study on the effect of aluminium foil on packaging thermal insulation performance in cold chain logistics[END_REF], which is not adapted to our work. Thus, we developed a lumped model including internal heat convection as explained above. To our knowledge, no studies have proposed such type of model.

The zonal models provide information of spatial temperature variation inside an insulated box with PCM by dividing the domain into several zones and applying heat balance equations between them. This approach is unique for an insulated box with PCM compared to other modelling approaches. Our previous work has developed the zonal model for an insulated box with PCM on a side wall [START_REF] Leungtongkum | Simplified heat transfer model for real-time temperature prediction in insulated boxes equipped with a phase change material[END_REF]. Hence, we present here further applications of this approach for real practice and a new model version adapted for PCM at top.

The CFD models are widely used to predict product temperature thanks to their flexibility although they require high computing resource [START_REF] Paquette | Modeling the thermal performance of a multilayer box for the transportation of perishable food, Food Bioprod[END_REF]. We consider CFD model as an approach for comparison with lumped and zonal models in terms of applicability and limitation.

Each model needs different types of input data, and thus has a different degree of complexity and provides different outputs. The users of the model (scientists, manufacturers, and stakeholders) have various levels of expertise and different practical questions: the amount of PCM needed to maintain the recommended temperature throughout a supply chain, the product temperature evolution, the warmest/coldest temperatures and their positions, etc.

This article aims to present three different modelling approaches from basic to advanced: lumped, zonal and CFD models to answer these questions. The capabilities and limitations of each model related to the insulated boxes with PCM for food transport are discussed. This study is original since it is the first time that three modelling approaches are compared, and their ability to answer certain questions is demonstrated. It is to be emphasized that the zonal and CFD model development was presented in our previous studies [START_REF] Leungtongkum | Simplified heat transfer model for real-time temperature prediction in insulated boxes equipped with a phase change material[END_REF][START_REF] Leungtongkum | CFD modelling of heat transfer and airflow in an insulated box equipped with Phase Change Material[END_REF].

Modelling approaches 2.1 Lumped model 2.1.1 Model descriptions and assumptions

An insulated box equipped with PCM loaded with product can be schematized by using electrical analogy (Fig. 1). The following assumptions were applied to develop the lumped model:

-The thermal inertia of air is neglected.

-The product is assumed to be a lumped object i.e., uniform temperature. is the heat transfer resistance between the product and the PCM [K‧W -1 ] 𝑅 𝑝/𝑝𝑐𝑚

Product temperature evolution

The energy balance of the product can be written as indicated in Eq. ( 5).

C p.p is specific heat of product [J‧kg -1 ‧K -1 ]
If T e and T pcm are constant and the characteristic time [s] is defined in Eq. ( 6), 𝜏 At steady state, if PCM is still melting (T pcm = T m ), the product temperature (T p ) reaches an equilibrium value (T p.eq ) related to the unknown thermal resistance R pcm/p , which can be represented by Eq. ( 9). ( 9)

Model validation by experiment

Materials

The insulated box for model validation was a 45-L commercialized multilayer insulated box (Manutan SA, Gonesse, France) with 500 mm x 310 mm x 300 mm internal dimensions (Fig. 2). The walls, containing four layers, consisted of three materials: expanded polystyrene (25 mm thickness), polypropylene (inner and outer layers with a thickness of 3.5 mm) and air gap between the expanded polystyrene and the inner layer (estimated thickness: 5 mm). The measured heat transmission coefficient (K) of this box, by internal heating method [START_REF] Atp | International Carriage of Perishable Foodstuffs and on the Special Equipment to be Used for Such Carriage: (ATP) as amended 6 July 2020[END_REF], was 0.90 W•m -2 •K -1 .

To allow uniform heat exchange with the ambient, the box was put on a 50-mm height wooden support. This box was placed on a 0.7-m height table located in the center of the temperaturecontrolled test room (Width 3.4 m x Length 3.4 m x Height 2.5 m). The PCM plate (external dimensions 460 mm x 280 mm x 47 mm) with an enclosure made of polypropylene (2.5-mm thickness) and filled with 3.5 kg of tap water (melting point ~ 0°C). 16 kg of Tylose slabs (dimensions 400 mm x 200 mm x 200 mm) was used as a test product. It contained 23% methyl hydroxyethylcellulose, 76.4% water and 0.5% NaCl (Refrigeration Development and Testing Ltd., North Somerset, UK).

From the box geometry and the material conductivities, two heat transfer resistances can be determined: R e/pcm = 3.68 K‧W -1 and R e/p = 1.67 K‧W -1 . Fig. 2: Insulated box used in the experimental validation (Source: Leungtongkum et al. [START_REF] Leungtongkum | Experimental investigation of airflow and heat transfer by natural convection in an insulated box with a Phase Change Material using a Particle Image Velocimetry technique[END_REF]).

Experimental protocol

The PCM was frozen in a freezer (-2°C) for at least 48 h prior to the experiment. Although it was horizontally placed, the thickness of the PCM slab was not uniform (35 mm to 50 mm). The test products were placed in a polystyrene box and put into a domestic refrigerator with setting temperature of 4°C for at least 24 h. Temperatures inside the PCM, product and air at different positions were measured using calibrated thermocouples with a precision of ± 0.1°C [START_REF] Leungtongkum | Experimental investigation of airflow and heat transfer by natural convection in an insulated box with a Phase Change Material using a Particle Image Velocimetry technique[END_REF]. The average product temperature was determined from 16 measurement positions.

Product equilibrium temperature determination

The equilibrium temperature was determined experimentally by replacing the PCM with a completely frozen one every 12 h over a period of 72 h. For the box with PCM on a side wall and 20°C ambient temperature, the average product temperature at equilibrium (T p.eq ) was 8.6°C. Thus, the thermal resistance between product and PCM (R p/pcm ) was 1.26 K‧W -1 .

Product temperature evolution validation

The temperature evolution inside the box was measured with a completely frozen PCM slab for the product initial temperature of 4.4°C under 20°C ambient. Fig. 3a compares the measured and calculated temperature evolution. It can be seen that the lumped model can reliably predict the product average temperature evolution. The possible deviation from the measured values can be caused by the neglect of external heat convection coefficient in thermal resistance calculation which leads to overestimation of the heat flux coming to product.

Due to the assumption of the lumped object, the model does not provide spatial temperature distribution data.

Validation of the evolution of mass of melted PCM

The amount of melted PCM from 0 h to 24 h was measured immediately after taking it out of the box. Fig. 3b shows good agreement between the measured and calculated amounts of melted PCM by the lumped model (Eq. 13). In our experiment, the characteristic time (τ) was about 8.8 h; after this time, the asymptote lumped model (Eq. 14) gives a good approximation. Since the PCM melted mass measurement provided only one value at each time point, slight differences of product initial temperature between the experiments could lead to deviation from the predicted values.

Model description and assumptions

This model applies for an insulated box with PCM on a side wall and a gap below the product. As suggested by our previous experimental study [START_REF] Leungtongkum | Influence of use conditions on heat transfer in an insulated box equipped with a phase change material[END_REF], the following assumptions were applied:

-2D-airflow path is considered (Fig. 4).

-There is a temperature difference between the top and the bottom of the box because of thermal stratification, and between the cold side (near the PCM) and the warm one.

-By taking conductive and convective heat exchange into account, there is also a temperature difference between product surface and product core.

-The thermal inertia (mass multiplied by the heat capacity) of air is considered as negligible compared with that of the product and box walls.

van der Sman [START_REF] Van Der Sman | Simple model for estimating heat and mass transfer in regular-shaped foods[END_REF] suggested that, when heat conduction and convection occur simultaneously, distinguishing the shell and core temperatures gave more accurate results than assuming a uniform temperature. In our zonal model, the product (parallelepiped shape with a length L) was divided into four blocks. Each block had a height and width of C with a shell (C/4) and core (3C/4). Following state variables describe 12 solid zones at a given time:

-T sh.n when n  [START_REF] Nie | Review on phase change materials for cold thermal energy storage applications[END_REF][START_REF] Laguerre | Experimental investigation and modelling in the food cold chain: Thermal and quality evolution[END_REF] represents average temperature in the product shell.

-T c.n when n  [START_REF] Nie | Review on phase change materials for cold thermal energy storage applications[END_REF][START_REF] Laguerre | Experimental investigation and modelling in the food cold chain: Thermal and quality evolution[END_REF] represents average temperature in the product core.

-T w.n when n  [START_REF] Nie | Review on phase change materials for cold thermal energy storage applications[END_REF][START_REF] Laguerre | Experimental investigation and modelling in the food cold chain: Thermal and quality evolution[END_REF] represents wall temperature, where T w.1 is the surface temperature of PCM and T w.2 to T w.4 are the temperatures of the internal walls.

T pcm and  define PCM temperature and its ice fraction, respectively. Air temperature consists of eight values representing the evolution along the airflow path. By exchanging with the product shells, air temperature varies from T a.n to T' a.n when n  [START_REF] Nie | Review on phase change materials for cold thermal energy storage applications[END_REF][START_REF] Laguerre | Experimental investigation and modelling in the food cold chain: Thermal and quality evolution[END_REF]. While its temperature changes from T' a.n to T a.n+1 when n  [START_REF] Nie | Review on phase change materials for cold thermal energy storage applications[END_REF][START_REF] Robertson | Innovations in cold chain equipment for immunization supply chains[END_REF] and from T' a.4 to T a.1 by exchanging with the internal walls.

More details regarding model development are presented in Leungtongkum et al. [START_REF] Leungtongkum | Simplified heat transfer model for real-time temperature prediction in insulated boxes equipped with a phase change material[END_REF]. Appendix 1 presents some equations representing air or product heat balances.

Model validation

Fig. 5 compares the measured and calculated values (by zonal model) under steady state (Fig. 5a) and transient state (Fig. 5b). The model provides a good prediction of the temperature distribution and temperature evolution. It is to be highlighted that steady state was the condition when PCM temperature was always at its melting temperature (T pcm = T m ), thus temperatures inside an insulated box would reach a constant value after a certain time. While transient state includes 3 periods: at the beginning when PCM temperature increases from its initial temperature to its melting point, then its temperature remains constant and at the end after PCM is completely melted, temperature increases again. Under transient state, temperatures inside the insulated box always change.

The discrepancies of the predicted values from the measured ones can be due to i) the difference between point value measurement and zone-averaged value prediction, ii) underestimation of thermal stratification at the top of the box and iii) loading procedure during which the product surface exposed to the ambient temperature. More detail can be found in our previous article [START_REF] Leungtongkum | Simplified heat transfer model for real-time temperature prediction in insulated boxes equipped with a phase change material[END_REF].

As stated in the previous section, this model applies for an insulated box with PCM on a side wall and a gap below the product. Thus, changing the PCM position or the product arrangement inside the box significantly impacts the model assumptions, hence heat balance equations. The demonstration of applying zonal model for a different PCM position is shown in Section 3.7.1.

Steady state

Model description and assumptions

In our CFD model, the following assumptions were applied:

-Laminar airflow as its Rayleigh number is lower than 10 9 (calculation not shown).

-Boussinesq approximation (density is assumed constant except in the gravity term).

-Viscous dissipation into heat is neglected.

-PCM density is assumed to be constant as that of its liquid state.

Governing equations of the CFD model

For air, laminar flow caused by natural convection were applied (the parameters without index in this section are for air and ρ is the air density at the reference temperature):

Continuity: (15) ∇ • ( 𝝆𝒗 ) = 𝟎 Energy: ( 16
) ∂ ∂𝒕 (𝝆𝑪 𝒑 𝑻) + ∇ • (𝝆𝒗𝑪 𝒑 𝑻) = ∇ • (𝝀∇𝑻) Momentum: ( 17
) ∂ ∂𝒕 ( 𝝆𝒗 ) +∇ • ( 𝝆𝒗𝒗 ) = -∇𝑷 +∇ • ( 𝝁∇𝒗 ) + 𝝆𝜷(𝑻 -𝑻 𝒓𝒆𝒇 )𝐠
For PCM, in our previous work [START_REF] Leungtongkum | CFD modelling of heat transfer and airflow in an insulated box equipped with Phase Change Material[END_REF], we assumed that there was no convection in melted PCM.

To better represent the physical phenomena inside PCM, the approach proposed by Voller and Prakesh [START_REF] Voller | A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems[END_REF] was used. PCM was assumed to be solid below T sol , liquid above T liq and in a mushy state between T sol and T liq . The fraction of melted PCM () is given by Eq. 18. The internal energy (u pcm ) of PCM is given by Eq. 19.

(18)

𝜽 = { 𝟎 𝑻 𝒑𝒄𝒎 < 𝑻 𝒔𝒐𝒍 𝟏 𝑻 𝒑𝒄𝒎 > 𝑻 𝒍𝒊𝒒 𝑻 𝒑𝒄𝒎 -𝑻 𝒔𝒐𝒍 𝑻 𝒍𝒊𝒒 -𝑻 𝒔𝒐𝒍 𝑻 𝒔𝒐𝒍 ≤ 𝑻 𝒑𝒄𝒎 ≤ 𝑻 𝒍𝒊𝒒 (19) 𝒖 𝒑𝒄𝒎 = ∫ 𝑻 𝑻 𝒓𝒆𝒇 𝑪 𝒑.𝒑𝒄𝒎 𝒅𝑻 𝒑𝒄𝒎 + 𝜽𝑳 𝒇
When PCM is completely melted ( = 1), standard free convection laminar flow equations apply.

In the mushy state (0 <  < 1), flow is limited by a momentum source term (the opposite to PCM velocity). When PCM is completely solid ( = 0), this source term becomes very high so that the PCM velocity disappears.

Continuity:

(20)

∇ • (𝝆 𝒑𝒄𝒎 𝒗 𝒑𝒄𝒎 ) = 𝟎 Energy: (21) ∂ ∂𝒕 (𝝆 𝒑𝒄𝒎 𝒖 𝒑𝒄𝒎 ) + ∇ • (𝝆 𝒑𝒄𝒎 𝒗 𝒑𝒄𝒎 𝒖 𝒑𝒄𝒎 ) = ∇ • (𝝀 𝒑𝒄𝒎 ∇𝑻 𝒑𝒄𝒎 ) Momentum: ∂ ∂𝒕 (𝝆 𝒑𝒄𝒎 𝒗 𝒑𝒄𝒎 ) +∇ • (𝝆 𝒑𝒄𝒎 𝒗 𝒑𝒄𝒎 𝒗 𝒑𝒄𝒎 ) = -∇𝑷 +∇ • (𝝁 𝒑𝒄𝒎 ∇𝒗 𝒑𝒄𝒎 ) + 𝝆 𝒑𝒄𝒎 𝜷 𝒑𝒄𝒎 (𝑻 -𝑻 𝒓𝒆𝒇 )𝐠 + 𝑺 (22)
where ( 23)

𝑺 = - (𝟏 -𝜽) 𝟐 (𝜽 𝟑 + 𝝐) 𝑩 𝒎𝒖𝒔𝒉 𝒗 𝒑𝒄𝒎
is a small number (0.001) to prevent division by zero. 𝝐 B mush is the mushy zone constant [kg‧m -3 ‧s -1 ] which is 10 5 in our calculation as in previous studies [START_REF] Ray | Comparative study of cooling performance for portable cold storage box using phase change medium[END_REF][START_REF] Rahimi-Khoigani | Application of an improved latent heat storage system in the food packaging[END_REF]. Since it was reported that this coefficient could influence the predictions [START_REF] Ebrahimi | Sensitivity of Numerical Predictions to the Permeability Coefficient in Simulations of Melting and Solidification Using the Enthalpy-Porosity Method[END_REF], some tests were carried out with other values (from 10 3 to 10 9 ) but in our case, the impact on PCM liquid fraction and product temperature evolution was not significant.

For the test product (Tylose), conduction alone was applied: Energy:

(24)

∂ ∂𝒕 (𝝆 𝒑 𝑪 𝒑.𝒑 𝑻) +∇ • (𝝆 𝒑 𝑪 𝒑.𝒑 𝑻) = ∇ • (𝛌 𝒑 ∇𝑻)
The following boundary conditions were applied:

At the internal box walls, the following momentum and thermal conditions were applied:

No slip boundary condition: Since air was considered as transparent (optical thickness = 0), so surface-to-surface radiation was activated, and the radiative flux entering the surface k coming from all the other surfaces j was calculated from (28) 𝒒 𝒊𝒏.𝒌 = ∑ 𝒋 𝑭 𝒋𝒌 𝒒 𝒓𝒂𝒅.𝒋

The external area of the box was higher than the internal one; thus, the geometric means of an external and internal area must be used to obtain the overall heat balance. Since the boundary condition applied on the internal wall, K in governing equations was corrected by a factor of for the CFD approach. 𝑨 𝒆 /𝑨 𝒊

Numerical simulation

The geometry was drawn by using SpaceClaim and meshed with Ansys Fluent meshing. A mesh independence study was first conducted by comparing the results obtained for the mesh number varying from 4.3 x 10 3 to 1.4 x 10 6 polyhedral cells (Fig. 6a). When the mesh passed from 2.6 x 10 5 to 1.4 x 10 6 cells, the results were very similar, for example the variation of average air temperature (adiabatic boundary condition at product surface, t = 10 min) was less than 0.1°C (Fig. 6a). Therefore, the mesh with 2.6 x 10 5 polyhedral cells was chosen for the reference box (cf. Section 3).

The numerical study was performed with Ansys Fluent 2021 R1 in a transient state. The results were very similar when applying a different time step from 0.1 s to 10 s (Fig. 6b) for the average air temperature. As there was no significant variation among chosen time steps, the time step of 1 s was used as the compromise between the calculation time and the number of iterations necessary to obtain convergence for one time step.

The model was solved using the SIMPLE (Semi IMplicit Pressure Linked Equation) method for pressure-velocity coupling. The spatial discretization was Second Order Upwind. The transient formulation was First Order Implicit. The convergence criteria for residual of continuity, velocity and radiation were set to 1 x 10 -3 , while it was 1 x 10 -6 for energy residual.

Summary of input and output parameters

The input parameters needed, and the outputs provided for each model are summarized in Table 1. To develop the lumped and zonal models, some measurements and parameter estimations are necessary, but none are required for CFD model development.

Model applicability

Since each model is based on different assumptions, Table 2 shows the applicability of each model to answer the technical questions.

Results and Discussions

To demonstrate the model applicability, the conditions/parameters shown in Table 3 were used as input parameters unless otherwise indicated. The demonstrated box called the "reference box" was equipped with polyurethane insulation with a thickness of 40 mm; it was shown that this box allowed the transport of temperature-sensitive products for periods of up to 96 hours [START_REF] Kacimi | Vacuum Insulated panels (VIP) in insulated packaging[END_REF].

Temperature and PCM evolution under varying ambient conditions

The ambience usually varies according to geographical location, the time of day and the season. Thus, taking the variation of ambience into account is necessary for transport design [START_REF] Navaranjan | Thermal insulation requirements and new cardboard packaging for chilled seafood exports[END_REF][START_REF] East | A method for quantitative risk assessment of temperature control in insulated boxes[END_REF][START_REF] Fioretti | A refrigerated container envelope with a PCM (Phase Change Material) layer: Experimental and theoretical investigation in a representative town in Central Italy[END_REF].

The numerical models taking into account the variable ambient temperature can be useful for designing the box and for estimating the PCM amount required to preserve food [START_REF] Kacimi | Combination of vacuum insulation panels and phase change materials in temperature-controlled containers[END_REF].

Fig. 9 shows the results obtained with the lumped model used to predict the average product temperature (T p ), PCM temperature (T pcm ), and amount of melted PCM (m pcm.melted ) under varying ambient profiles (T e ) adapted from Fioretti et al. [START_REF] Fioretti | A refrigerated container envelope with a PCM (Phase Change Material) layer: Experimental and theoretical investigation in a representative town in Central Italy[END_REF]. The model highlights the various temperature changes when the ambient temperature alters. During the first 9.5 h (T e = 14°C), the product temperature increased slowly until it reached 4.3°C (product equilibrium temperature 4.6°C, initial temperature 4.0°C) and the melting rate of the PCM was 0.06 kg/h. Then, the product temperature increased to 7.3°C during the following 12 h (T e = 28°C) with a product equilibrium temperature of 9.2°C, and the PCM melting rate was 0.11 kg/h. This rate was almost twice that observed during the first period as the temperature difference between the ambience and the PCM doubled from 14°C to 28°C. The product temperature decreased slightly during the following 1.5 h (T e = 20°C), while the product equilibrium temperature was 6.6°C. However, when the PCM was completely melted (23 h and thereafter), the PCM and the product temperatures increased. Despite the decreasing ambient temperature which dropped to 14°C during the final 8 h, the product temperature constantly increased, unlike during the first period, because the PCM temperature rose continuously after it was completely melted.

PCM melting time and maximum storage period

Section 3.1 demonstrates the determination of the PCM melting time (t melt ) and maximum storage period (t max ) so that the product temperature remains below the maximum storage temperature (T max ) for given PCM mass, product initial temperature and external temperature.

As shown in Fig. 9, the PCM melting time for m pcm = 2 kg is 23 h, and the maximum storage period for T max = 8°C is 30 h. However, in practice, stakeholders usually have information regarding the targeted transport duration (t max ), the maximum storage temperature not to be exceeded, and the external temperature. Thus, they need to determine the PCM amount required for each transport operation.

Here, we considered a constant ambient temperature. First, we examined the case where T p.eq < T max . The melting time (t melt ) can be estimated from Eq. 13 by adding the sensible heat of the PCM from T pcm.0 to T m as shown in Eq. 29.

( One could consider that the transport duration should not exceed the melting time. Therefore, it can be assumed that t melt as t max gives the first estimation of the necessary mass of PCM; but this leads to an overestimation of the PCM mass because the thermal inertia of the product and the PCM is underestimated. Thus, a better estimation can be obtained by including product thermal inertia from T p.0 to T max and the sensible heat of melted PCM from T m to T max as described in Eq. 30.

(

) 30
𝑚 𝑝𝑐𝑚 (𝐿 𝑓 + 𝐶 𝑝𝑐𝑚.𝑠𝑜𝑙 (𝑇 𝑚 -𝑇 𝑝𝑐𝑚.0 ) + 𝐶 𝑝𝑐𝑚.𝑙𝑖𝑞 (𝑇 𝑚𝑎𝑥 -𝑇 𝑚 )) + 𝑚 𝑝 𝐶 𝑝.𝑝 (𝑇 𝑚𝑎𝑥 -𝑇 𝑝.0 ) ≈ ( 𝑇 𝑝.𝑒𝑞 -𝑇 𝑚 𝑅 𝑝𝑐𝑚/𝑝 + 𝑇 𝑒 -𝑇 𝑚 𝑅 𝑒/𝑝𝑐𝑚 ) 𝑡 𝑚𝑎𝑥 Fig. 10a shows the melting time and maximum storage period for T max = 8°C in the reference box (cited in Table 3) under 20°C ambient conditions (T p.eq = 6.6°C). It can be seen that the maximum storage period was longer than the melting period thanks to product and PCM thermal inertia. The estimations of melting time and maximum storage period by Eq. 29 and 30 are very close to that obtained by direct transient simulation as shown in Section 3.1. The advantage of Eq. 30 is that it gives an analytical expression of the required PCM amount as a function of T p.0 , T pcm.0 , T e , T max , t max and box insulation (through R e/pcm ). This equation can be applied only when T max is higher than the product equilibrium temperature (T p.eq ).

In some cases, for example, transport conducted using a poorly insulated box or under high ambient temperatures, T p.eq can be higher than T max . In this case, the product temperature reaches the maximum value (T max ) at a given time (t' max ) before it reaches equilibrium, and t' max can be determined from Eq. 31 as follows:

(31)

𝑇 𝑚𝑎𝑥 = 𝑇 𝑝.𝑒𝑞 +(𝑇 𝑝.0 -𝑇 𝑝.𝑒𝑞 )𝑒 - 𝑡 ′𝑚𝑎𝑥

𝜏

If the targeted maximum storage period (t max ) is longer than t' max , the product temperature will exceed T max whatever the PCM amount. On the other hand, if t max is lower than t' max , the required PCM amount can still be determined by using Eq. 30.

Fig. 10b shows the melting time and maximum storage period for T max = 8°C in the reference box under 30°C ambient temperature (T p.eq = 9.9°C). Since t' max is 14.5 h, it is impossible to maintain T p below T max after 14.5 h, whatever the PCM mass used. If it is necessary to achieve a longer transport period under this ambient temperature, a better-insulated box and/or PCM with lower melting point are needed.

Equations 29 and 30 can be used only under constant external temperature. For varying ambient conditions, the direct approach should be used by solving t melt or t max for various PCM masses, then interpolating to determine the required amount at the targeted period (Fig. 9). The product and PCM initial temperatures were 4°C and -2°C, respectively. The blue curves also indicate the melting period (t melt ).

Effect of the box insulation

Box insulation is the main criterion for transport performance in an insulated box with PCM as emphasized by several studies [START_REF] Du | Cooling performance of a thermal energy storage-based portable box for cold chain applications[END_REF][START_REF] Leungtongkum | Simplified heat transfer model for real-time temperature prediction in insulated boxes equipped with a phase change material[END_REF][START_REF] Paquette | Modeling the thermal performance of a multilayer box for the transportation of perishable food, Food Bioprod[END_REF][START_REF] Margeirsson | Experimental and numerical modelling comparison of thermal performance of expanded polystyrene and corrugated plastic packaging for fresh fish[END_REF]. Equations 5, 9 and 11 describe the relationship between the box insulation (via thermal resistances) and product temperature evolution, product temperature at equilibrium and PCM melting rate, respectively. From Fig. 11, box insulation exerted a significant impact on temperature change in an insulated box. The box with less effective insulation (higher K value) led to a higher rate of temperature rise since the characteristic time (τ) was lower (12.6 h and 10.6 h for the reference box and the experimental box, respectively) and the equilibrium temperature (T p.eq ) was higher (6.6°C and 8.6°C for the reference box and experimental box, respectively) (Fig. 11a). The product temperature was reduced in the box with better insulation (VIP) as the equilibrium temperature (T p.eq = 2.5°C) was lower than the initial product temperature (4°C). A higher PCM melting rate was also observed in a box with a higher K value (Fig. 11b).

Effect of melting temperature and latent heat of PCM

Cold energy storage by the PCM is another key factor enabling a low temperature to be maintained in an insulated box [START_REF] Yang | Recent development on heat transfer and various applications of phase-change materials[END_REF]. There are numerous PCM materials to choose from, e.g., water, salt solution, salt hydrate, paraffin, commercially available PCM [START_REF] Oró | Review on phase change materials (PCMs) for cold thermal energy storage applications[END_REF]. Each one has different thermophysical and chemical properties. The main concerns for transport are the PCM melting temperature and its latent heat. The melting temperature determines the equilibrium temperature inside the box (Eq. 9) and the latent heat determines the time interval during which the PCM can provide cold (Eq. 11).

Fig. 12 demonstrates, by simulation, the impact of the PCM material on the product temperature (T p.0 = 4°C) and the melted PCM mass evolution assuming that the PCM is not completely melted. Three PCM materials were chosen for comparison: ice (T m = 0°C and L f = 333700 J•kg - 1 ), an eutectic Mg 2 SO 4 solution (T m = -3.9°C and L f = 264400 J•kg -1 ) and K 2 HPO 4 salt hydrate (T m = 4°C and L f = 109000 J•kg -1 ) [START_REF] Cengel | Heat and Mass Transfer: Fundamentals & Applications[END_REF][START_REF] Li | Review of cold storage materials for subzero applications[END_REF][START_REF] Sharma | Review on thermal energy storage with phase change materials and applications[END_REF]. It can be seen that using ice and K 2 HPO 4 salt hydrate led to a product temperature increase since the equilibrium temperature for each PCM was 6.6°C and 9.3°C for ice and salt hydrate, respectively. Despite the same characteristic time in each case, the temperature rise in the box with salt hydrate was the highest. The box filled with eutectic PCM did not show a temperature change as T p.eq = 4°C, i.e. the same value as the initial product temperature. The increase in melted PCM mass is inversely related to its latent heat, e.g., ice melts at the lowest rate compared with the other two.

Effect of product mass and thermophysical properties

As explained previously, the thermal inertia of the product plays an important role in the temperature evolution and PCM melting rate. Two aspects can be considered: product mass and the thermophysical properties of the product.

Fig. 13a shows the influence of salmon mass (4, 8 and 16 kg) on the product temperature evolution, and Fig. 13b shows the influence of the nature of the product (salmon, butter and milk); thus, the specific heat varies: salmon (C p = 3360 J•kg -1 •K -1 ), butter (C p = 2080 J•kg - 1 •K -1 ), and milk (C p = 3960 J•kg -1 •K -1 ) [START_REF] Cengel | Heat and Mass Transfer: Fundamentals & Applications[END_REF]. It was found that higher thermal inertia, i.e., a higher mass and/or higher specific heat led to a lower temperature change as it allowed a higher characteristic time (τ) (defined in Eq. 6). The impact of thermal inertia on the PCM melting rate followed the same trend as that of temperature evolution (results not shown). The conditions with the lowest thermal inertia (4 kg of salmon and 16 kg of butter in Figures 13a and13b, respectively) resulted in the lowest PCM melting time, so the product temperature sharply increased at the end of the simulation, implying that the PCM was completely melted.

Changing the product mass also affects the occupied volume in the box, which in turn may impact internal convective heat transfer, hence, the thermal resistance between the product and PCM (R p/pcm ). Caution must also be taken if the product is too large since in this case it obstructs the internal airflow by natural convection [START_REF] Leungtongkum | Experimental investigation of airflow and heat transfer by natural convection in an insulated box with a Phase Change Material using a Particle Image Velocimetry technique[END_REF].

Effect of internal radiation

Since heat exchange by conduction, natural convection and radiation in a cavity with cold and warm surfaces are of the same order of magnitude [START_REF] Laguerre | Temperature prediction in domestic refrigerators: Deterministic and stochastic approaches[END_REF], the inclusion of all these heat transfer modes can improve the accuracy of the thermal model. It is to be highlighted that the lumped model could not take radiation into account as it does not apply the difference between the surface and the internal temperatures. Consequently, the zonal model is used for demonstration purposes.

Fig. 14 compares the impact of the internal wall and the PCM surface emissivity on the product core temperature at steady state. Three wall emissivity configurations were numerically studied: case 1: box and PCM walls were bright (ε = 0.03, polished aluminum); case 2: box walls were bright (ε = 0.03) while the PCM wall was not (ε = 0.97, polypropylene); case 3: none of the walls were bright (ε = 0.97). The highest core temperature was observed in case 1; indeed, when every surface was bright, the heat flux from ambient to the product was reduced, but the heat flux from the PCM to the product decreased to a greater extent. The lowest core temperature was observed in case 2; in this manner, the transfer from ambient to product to ambient was reduced but not the transfer from the PCM to the product. These predictions were confirmed by an experiment undertaken by covering different walls with aluminum foil.

The results presented in Fig. 14 confirm that internal radiation must be taken into account for food transport in an insulated box with PCM to prevent temperature abuse.

Fig. 14: Effect of internal emissivity on the product core temperature in an insulated box with PCM on a side wall predicted by the zonal model under steady state

Modelling approaches for a box with PCM at the top

The position of the PCM is another factor that affects temperatures in an insulated box [START_REF] Du | Cooling performance of a thermal energy storage-based portable box for cold chain applications[END_REF]. However, the lumped model previously developed for the box with PCM on a side wall cannot directly apply when the PCM is at the top. In fact, when the PCM position is changed, thermal resistances are changed as well, and this results in different product temperature equilibria. According to our experiment using the box loaded with the product under 20°C ambient conditions, the product equilibrium temperature decreased from 8.6°C (PCM on a side wall) to 8.2°C (PCM at the top). Hence, the zonal model and the CFD model are used to study the impact of the PCM position in this section.

Zonal model

The assumption of airflow pattern in an insulated box is essential for zonal model development; thus, the model is posteriori. Changing the PCM location and consequently the airflow pattern leads to different heat transfer and heat balance equations.

Our previous experimental data showed that airflow and temperature fields in an insulated box loaded with the product and with PCM at the top are asymmetrical in spite of the geometrical symmetry [START_REF] Leungtongkum | Influence of use conditions on heat transfer in an insulated box equipped with a phase change material[END_REF]. Thus, two zonal model approaches are demonstrated and compared in this section. First, we applied the same heat transfer and airflow as those previously developed for the insulated box with PCM on a side wall by redefining the state variables (Fig. 15a). Another approach was developed from a preliminary result of 2D CFD in a rectangular cavity with a cold surface at the top. The CFD results indicated two airflow loops in a loaded box; thus, the heat transfer and airflow diagram were modified (Fig. 15b). There was a primary airflow loop (black continuous line in Fig. 15b) with a mass flow rate of and a secondary airflow loop 𝑚 𝑎 (black dashed line in Fig. 15b) of . The air mixing was estimated to be at point 1 located at 𝛿𝑚 𝑎 the mid-height of the box, while the heat transfer by conduction, convection and radiation were identical to those obtained with the previous model.

The zonal model and the adapted zonal model gave good agreement with the average experimental values obtained at 24 positions (internal walls, internal air and product) with a root mean square error of 0.83°C and 0.79°C, respectively. Fig. 15c shows the comparison between the experimental product core temperature and that predicted by each model. It was found that the adapted zonal model gave a slightly better prediction of the product core temperature.

Zonal model

Adapted zonal model

CFD model

Since CFD model algorithm directly solves Navier-Stokes and energy partial differential equations at each position [START_REF] Söylemez | CFD analysis for predicting cooling time of a domestic refrigerator with thermoelectric cooling system[END_REF], so the model is priori, i.e., the users do not need assumptions on physical phenomena to develop this model. Thus, it can directly demonstrate the impact of PCM position in an insulated box in which the airflow is initially unknown, unlike the zonal model.

Detailed airflow pattern and effect of product compactness

As previously discussed, the CFD model can provide the detailed airflow and heat transfer phenomena without knowledge of model assumptions. This information is necessary to fully understand the phenomena, especially for complex configurations, e.g., irregular product shape, porous product and PCM on several surfaces. Fig. 17a shows the airflow pattern in the air gap near the PCM located on a side wall of an insulated box with PCM (y = 245 mm). Air flew downwards and toward the lateral walls in this gap. Fig. 17b shows a lateral view (x =15 mm): air flew from the left (near the PCM) to the right (near the lateral walls). These figures confirm that there is 3-D airflow as illustrated in Fig. 17c: air flew downward near the cold wall (PCM) and toward the side walls; then it flew upward along the opposite wall and returned near the PCM. The CFD model can also be used for more complex configurations by simply changing the input geometry.

The product to be transported in food cold chains can vary in its compactness related to the nature of the product and/or packaging geometry. Since internal convective heat transfer cannot be neglected, different levels of compactness lead to different heat transfer coefficients between the internal air and the product surface. Laguerre et al. [START_REF] Laguerre | Transient heat transfer by free convection in a packed bed of spheres: Comparison between two modelling approaches and experimental results[END_REF] conducted CFD simulation for a porous load (spherical fruits or vegetables). They compared the simulation results two approaches: a direct approach with meshing of all the product items, or a porous media approach. Both approaches gave good agreement with the experimental data.

Effect of box dimensions

The boxes used in food transport can vary in volume and aspect ratio (H/L). Fig. 18a shows the effect of different aspect ratios of the box on the average load temperature evolution. In the case of the same percentage of occupied volume inside the box, changing the aspect ratio did not affect the average temperature evolution. Thus, the same characteristic time (τ) was obtained. Fig. 18b shows the temperature evolution at the cold and warm points for different H/L. Increasing H/L caused greater thermal stratification with temperature differences between the cold and warm points at 24 h of 1.5°C and 2.3°C for the box with H/L = 1.1 and 1.7, respectively (see the warm and cold points in Fig. 18c).

Prediction of quality change and sanitary risk during transport

Many studies have combined thermal models that predict the temperature evolution in different steps throughout the cold chain with quality and microbial growth models [START_REF] García | Assessment and Prediction of Fish Freshness Using Mathematical Modelling: A Review[END_REF][START_REF] Matar | Consumer behaviour in the prediction of postharvest losses reduction for fresh strawberries packed in modified atmosphere packaging[END_REF][START_REF] Onwude | Physics-driven digital twins to quantify the impact of pre-and postharvest variability on the end quality evolution of orange fruit[END_REF][START_REF] Penchaiya | Modelling quality and maturity of 'Namdokmai Sithong' mango and their variation during storage[END_REF]. This section highlights the possibility and limitations governing the use of the thermal models we developed to predict food quality alterations in a cold chain during which the ambient temperature may vary with time. In fact, if the practical objective is to evaluate the influence of the ambient temperature variability in the cold chain on the product through a Monte Carlo process [START_REF] Duret | Combining Quantitative Risk Assessment of Human Health, Food Waste, and Energy Consumption: The Next Step in the Development of the Food Cold Chain?[END_REF], a CFD model might not be appropriate because of the high calculation costs. Although some studies combined a quality model with a CFD model, they focused only on certain operating conditions in a given cold facility [START_REF] Wu | Identifying heterogeneities in cooling and quality evolution for a pallet of packed fresh fruit by using virtual cold chains[END_REF]. These studies cannot be applied when the full variability of the operating conditions in the entire cold chain must be taken into account.

In this section, we present the coupling of a predictive Listeria monocytogenes growth model (given that such growth may occur during cold storage and transport) with the lumped model for fresh salmon transport in an insulated box with PCM on a side wall. The first-order equation was used to describe the growth of L. monocytogenes taking into account the lag time as follows [START_REF] Baranyi | A dynamic approach to predicting bacterial growth in food[END_REF][START_REF] Jia | Assessing the growth of Listeria monocytogenes in salmon with or without the competition of background microflora --A one-step kinetic analysis[END_REF].

(

) 𝑑𝑌 𝑑𝑡 = 1 1 + 𝑒 -𝐸 × 𝜂 𝑚𝑎𝑥 32
where Y is the amount of L. monocytogenes [log 10 CFU‧g -1 ] E, the physiological state of the L. monocytogenes [-], can be described by Eq. 33.

( At t = 0: E 0 = -1.05 [START_REF] Kacimi | Vacuum Insulated panels (VIP) in insulated packaging[END_REF] Fig. 19 shows the evolution of PCM and the product (fresh salmon) temperatures and the growth of L. monocytogenes on the average under varying ambient temperatures adapted from Fioretti et al. [START_REF] Fioretti | A refrigerated container envelope with a PCM (Phase Change Material) layer: Experimental and theoretical investigation in a representative town in Central Italy[END_REF]. This figure indicates that an increase of L. monocytogenes corresponded to its temperature evolution. When PCM was melting (t ≤ 23 h), the product temperature increased slowly even during a drastic change in ambient temperature (10 h < t < 23 h); thus, microbial growth was low (0.11 log 10 CFU‧g -1 increase during the first 23 h). When the PCM was completely melted (t > 23 h), the product temperature increased rapidly, and this led to a sharp increase in L. monocytogenes, and growth reached 0.7 log 10 CFU‧g -1 at 48 h. Thus, it is essential to place sufficient PCM mass in the box to prevent microbial growth.

Regarding the sanitary risk, if the objective is to conduct a quantitative microbial risk assessment (QMRA) of fish product, e.g., the risk from L. monocytogenes, then the zonal model should be used. Indeed, the results of listeriosis cases in the L. monocytogenes QMRA model are linked to the extreme right values of the exposure distribution corresponding to risky situations [START_REF] Pouillot | Quantitative risk assessment of Listeria monocytogenes in French cold-smoked salmon: II. Risk characterization, Risk Anal[END_REF]. Hence, the thermal model used to predict the risk of listeriosis must indicate the temperature heterogeneity inside the box, since several sensitivity analysis studies highlighted the impact of temperature variability on the risk associated with exposure to L. monocytogenes [START_REF] Duret | Identification of the significant factors in food safety using global sensitivity analysis and the acceptand-reject algorithm: application to the cold chain of ham[END_REF][START_REF] Ellouze | Global sensitivity analysis applied to a contamination assessment model of listeria monocytogenes in cold smoked salmon at consumption, Risk Anal[END_REF].

Fig. 19: Product (fresh salmon), PCM temperature evolution and L. monocytogenes (Y/Y 0 ) increase on the average under varying ambient temperatures (T e ); the initial product and PCM temperatures were 4°C and -2°C, respectively.

Applications and limitations of each thermal model

The demonstrations of thermal model applications to design the box, choose the product and PCM arrangement have been discussed in depth using different models. Table 4 summarizes the characteristics of the models presented in this article in terms of input parameters, computational resources, outputs and model flexibility. To develop the model, the users need to establish the airflow and heat transfer assumptions in order to develop a zonal model, whereas this is not the case for CFD model development.

All models require information on the box, PCM and product dimensions and thermophysical properties. An additional experiment (or a CFD simulation) is needed to estimate the thermal resistance between the product and PCM (R p/pcm ) to develop a lumped model, since this value depends on the PCM and product arrangement inside the box. The zonal model requires an estimation or an experiment to determine the convective heat transfer coefficient (h) and the air mass flow rate ( ).

𝑚 𝑎

The lumped model and zonal model can be coded and solved with a free-licensed coding program such as Python, and do not require high computing resources; they can be solved in less than 1 minute. On the contrary, solving the CFD model requires a licensed program and experts to operate it. The CFD model also requires a longer calculation time, so it may not be suitable for investigating in real time the effect of box design combined with the operating conditions.

The lumped model can provide the temporal evolution inside the box but cannot provide spatial variations. The zonal model and the CFD model give more details on temperature heterogeneity during transport.

To summarize, one has to bear in mind that a complex model might not provide more information [START_REF] Zwietering | Quantitative risk assessment: Is more complex always better?: Simple is not stupid and complex is not always more correct[END_REF]. The choice of the model is then dependent on the application and the users.

Conclusion

Three validated thermal models (lumped, zonal and CFD models) demonstrated their ability to answer practical questions arising during food transport using insulated boxes with PCM. The lumped model, which predicts the average product temperature, can be used to illustrate the temporal evolution, to estimate the PCM amount required to maintain the product temperature (on the average) under the recommended values for a given period. This model can also be used to determine the effect of box insulation, PCM material, product mass and nature. The lumped model does not require high computing resources and it can be used as long as the PCM and product arrangement are identical to that used in the experiment undertaken for model validation. The zonal model can give a greater insight into heat transfer and airflow. This model does not need high computing resources, but it requires estimations of input parameters relating to internal heat transfer and airflow. The model can be adjusted for different PCM and product arrangements by relying on knowledge of heat transfer and airflow phenomena. The CFD model gives the most thorough comprehension of all physical phenomena occurring during transport without any assumptions, but the computational cost is high. To predict quality change and the sanitary risk by coupling predictive models with previously presented models, the practical objectives must be clear, since different models provide different insights and have different limitations. Studies on the application of these models to predict the evolution under real transport conditions can be useful for developing better boxes and guidelines on efficient transport to avoid temperature abuse.

Yλ 4 τδn

 4 Mushy zone constant [kg‧m -3 ‧s -1 ] C Length of the cross-section of the product block for the zonal model [m] C p Specific heat capacity [J‧kg -1 ‧K -1 ] E Physiological state of microorganisms [-] F jk View factor of the surface j relating to the surface k Gravitational acceleration = 9.81 m‧s -2 g h Convective heat transfer coefficient [W‧m -2 ‧K -1 ] H Height of the box [m] K Overall heat transfer coefficient [W‧m -2 ‧K -1 ] L Length of the product block or wall [m] L f Latent heat of fusion [J‧kg -1 ] m Mass [kg] Mass flow rate [kg‧s -1 ] 𝑚 MC p Thermal inertia [J‧K -1 ] Normal unit vector [-] 𝑛 P Pressure [Pa] Incoming radiative flux [W‧m -2 ] 𝑞 𝑖𝑛 Radiative flux from the surface [W‧m -2 ] 𝑞 𝑟𝑎𝑑 q r Radiative heat exchange [W] R Heat transfer resistance [K‧W -1 ] Momentum source term [kg‧m -2 ‧s -2 ] 𝑆 t Time [s] t max Maximum storage period [s] t melt PCM melting time [s] T, T' Temperature [°C or K] T m Melting temperature of PCM [°C or K] T max Maximum storage temperature [°C] Number of microorganisms [log 10 CFU‧g -1 ] x, y, z Coordinate [m] Thermal conductivity [W‧m -1 ‧K -1 ] φ Ice fraction [-] ρ Density [kg‧m -3 ] σ Stefan-Boltzmann constant = 5.67 x 10 -8 W‧m -2 ‧K -Characteristic time [s] 𝜂 Microbial growth rate [h -1 ] Ratio of air mass flow rate between the secondary and the primary airflow loops for the zonal model (see Fig. 15b) [-] Zone number of the zonal model (see Figures 4 and 15)

Fig. 1 :

 1 Fig. 1: Lumped model structure

  where m p is product mass[kg]

Fig. 3 : 2 . 2

 322 Fig. 3: Comparison between the experimental and numerical values (by lumped model) in a box with PCM on a side wall under 20°C ambient: (a) product temperature evolution with T p.0 = 4.4°C, and (b) melted PCM mass when T p.0 = 4°C 2.2 Zonal model

Fig. 4

 4 Fig. 4 illustrates the control volumes, airflow path and heat fluxes considered. More detail of model description and assumptions can be found in Leungtongkum et al. [27].

Fig. 4 :

 4 Fig. 4: Side view of a simplified heat transfer and airflow diagram in the zonal model of an insulated box with PCM on a side wall. Source: Leungtongkum et al. [27]

Fig. 5 :

 5 Fig. 5: Comparison between the experimental and numerical (by zonal model) temperatures in an insulated box equipped with PCM on one side wall with an initial product temperature of 4°C and ambient temperature of 20°C: (a) steady state, and (b) transient state. Source: Leungtongkum et al. [27] 2.3 CFD model

  ) 𝒗 = 𝟎 Cauchy type thermal boundary condition: (26) 𝑲(𝑻 𝒆 -𝑻) = 𝝀∇𝑻 • 𝒏 + 𝒒 𝒓𝒂𝒅 where is the overall heat transfer coefficient [W‧m -2 ‧K -1 ] 𝑲 and (27) 𝒒 𝒓𝒂𝒅 = -𝜺𝛔𝑻 𝟒 + (𝟏 -𝜺)𝒒 𝒊𝒏 with ε is the wall emissivity [-]

Fig. 6 :

 6 Fig. 6: Independence study of (a) mesh number and (b) time step on the average air temperature.

Figures 7 Fig. 7 :Fig. 8 :

 778 Figures 7 and 8 compare experimental and CFD results of temperature field and temperature evolution in a box with PCM on a side wall. It can be seen that the CFD model gave accurate predictions with a maximum difference of 2.0°C with the experimental value. The discrepancy, especially for air temperature, could be caused by thermal bridges in the structure of the insulating box walls. The detail of the wall structure is not taken into account in the CFD model. Since product is not entirely in contact with the walls, this has less impact on product temperature.

Fig. 9 :

 9 Fig. 9: Product (T p ) and PCM (T pcm ) temperatures, melted PCM mass evolution (m pcm.melted ) under varying ambient conditions (T e ) predicted by the lumped model; the product and PCM initial temperatures were 4°C and -2°C, respectively, PCM mass 2 kg.

Fig. 10 :

 10 Fig. 10: Required PCM amount (orange curves) as a function of the maximum storage period (t max ) for a maximum product temperature of 8°C by direct calculation and approximation under (a) constant 20°C ambient temperature, and (b) constant 30°C ambient temperature. The product and PCM initial temperatures were 4°C and -2°C, respectively. The blue curves also indicate the melting period (t melt ).

Fig. 11

 11 Fig.11shows the impact of different types of box insulation on the product temperature evolution and the amount of melted PCM. The demonstrated boxes are: reference box (polyurethane with a thickness of 40 mm, K = 0.58 W•m -2 •K -1 ), box with a vacuum-insulated panel (VIP box, vacuum-insulated panel with a thickness of 20 mm, and polyurethane with a thickness of 20 mm, K = 0.17 W•m -2 •K -1 )[START_REF] Kacimi | Vacuum Insulated panels (VIP) in insulated packaging[END_REF], and our experimental box shown in Fig.2(Exp box, K = 0.90 W•m -2 •K -1 ).

Fig. 11 :

 11 Fig. 11: Effect of box insulation on: (a) average product temperature evolution; and (b) the amount of melted PCM estimated by the lumped model. The initial product and PCM temperatures were 4°C and -2°C, respectively.

Fig. 12 :

 12 Fig. 12: Effect of the PCM type on: (a) product temperature evolution; and (b) PCM melted mass evolution calculated by the lumped model; the product initial temperature was 4°C.

Fig. 13 :

 13 Fig. 13: Effect of: (a) a mass of salmon, and (b) the nature of the product (16 kg of salmon, butter or milk) on the product temperature evolution calculated by lumped model for a PCM mass of 2 kg. The product and PCM initial temperatures were 4°C and -2°C, respectively.

Fig. 15 :

 15 Fig. 15: Heat transfer and airflow diagram in a 2-D insulated box with PCM at the top for (a) the zonal model; (b) the adapted zonal model; and (c) comparison between the average core temperature at 4 measured positions and numerical values under steady state using an experimental box

Fig. 16 Fig. 16 :

 1616 Fig.16illustrates the CFD results of air velocity and temperature fields on the middle plane of a loaded insulated box with PCM at the top. The air flew from the right to the left above the load. The anti-clockwise circulation loop explained the lower temperature in the left gap because of the downward flow coming from the PCM. It shows that the downward flow took place in one or another preferential pathway (either the left or right gap). This result confirms that despite symmetrical geometry, the airflow is asymmetrical; hence, it leads to an asymmetrical temperature field. This is due to the non-linear term in Navier-Stokes equations resulting in a symmetry rupture (even before turbulence appears).

Fig. 17 :

 17 Fig. 17: CFD results of air velocity field at t = 24 h (a) near an PCM surface, y = 245 mm (b) near lateral wall, x = 15 mm, and (c) 3-D airflow illustration (in a loaded box with PCM on a side wall)

Fig. 18 :

 18 Fig. 18: The effect of the aspect ratio (H/L) of the box on: (a) average load temperature evolution; and (b) load temperature evolution at the coldest and warmest points; (c) position of the coldest (■) and warmest (■) points.

2 For

 2 L. monocytogenes: 𝜂ref = 0.183 h -1 , T min = -2°C and T ref = 25°C[START_REF] Baranyi | A dynamic approach to predicting bacterial growth in food[END_REF][START_REF] Jia | Assessing the growth of Listeria monocytogenes in salmon with or without the competition of background microflora --A one-step kinetic analysis[END_REF]

  𝑇 𝑝 = 𝑇 𝑝.𝑒𝑞 +(𝑇 𝑝.0 -𝑇 𝑝.𝑒𝑞 )𝑒

	𝜏 =	𝑚 𝑝 𝐶 𝑝.𝑝 𝑅 𝑒/𝑝 1 1 + 𝑅 𝑝/𝑝𝑐𝑚
	Eq. (5) becomes	𝑑𝑇 𝑝 𝑑𝑡 +	𝑇 𝑝 𝜏 =	𝑇 𝑝.𝑒𝑞 𝜏	(7)
	Thus,				-	𝑡 𝜏	(8)

  PCM evolves in 3 states from solid state until it is completely melted. The energy balances of PCM can be written for each state as follows.State 1: when the PCM is completely solid (T pcm < T m ) and m pcm.melted = 0;

				𝑇 𝑒 -𝑇 𝑝.𝑒𝑞 𝑅 𝑒/𝑝 +	𝑇 𝑚 -𝑇 𝑝.𝑒𝑞 𝑅 𝑝𝑐𝑚/𝑝 = 0
	2.1.3 PCM evolution
		𝑚 𝑝𝑐𝑚 𝐶 𝑝.𝑝𝑐𝑚.𝑠𝑜𝑙	𝑑𝑇 𝑝𝑐𝑚 𝑑𝑡 =		𝑇 𝑒 -𝑇 𝑝𝑐𝑚 𝑅𝑒 𝑝𝑐𝑚 +	𝑇 𝑝 -𝑇 𝑝𝑐𝑚 𝑅𝑝𝑐𝑚 𝑝	(10)
	State 2: when the PCM is melting (T pcm = T m );
		𝐿 𝑓 𝑑𝑚 𝑝𝑐𝑚.𝑚𝑒𝑙𝑡𝑒𝑑 𝑑𝑡	=	𝑇 𝑒 -𝑇 𝑝𝑐𝑚 𝑅 𝑒/𝑝𝑐𝑚 +	𝑇 𝑝 -𝑇 𝑝𝑐𝑚 𝑅 𝑝𝑐𝑚/𝑝	(11)
	State 3: when the PCM is completely liquid (T pcm > T m ) and m pcm.melted = m pcm ;
		𝑚 𝑝𝑐𝑚 𝐶 𝑝.𝑝𝑐𝑚.𝑙𝑖𝑞 𝑑𝑇 𝑝𝑐𝑚 𝑑𝑡 =	𝑇 𝑒 -𝑇 𝑝𝑐𝑚 𝑅𝑒 𝑝𝑐𝑚 +	𝑇 𝑝 -𝑇 𝑝𝑐𝑚 𝑅𝑝𝑐𝑚 𝑝	(12)
	From Equations (8) and (11), if T e and T pcm are constant:
		𝑚 𝑝𝑐𝑚.𝑚𝑒𝑙𝑡𝑒𝑑 =	(𝑇 𝑝.𝑒𝑞 -𝑇 𝑝𝑐𝑚 )𝑡 + 𝜏(𝑇 𝑝.0 -𝑇 𝑝.𝑒𝑞 )(1 -𝑒 𝑅 𝑝𝑐𝑚/𝑝 𝐿 𝑓	-	𝑡 𝜏 )	+	𝑇 𝑒 -𝑇 𝑝𝑐𝑚 𝑅 𝑒/𝑝𝑐𝑚 𝐿 𝑓 𝑡	(13)
	Assuming	𝑡 ≫ 𝜏	(asymptotic approximation):
		𝑚 𝑝𝑐𝑚.𝑚𝑒𝑙𝑡𝑒𝑑 ≈	𝜏(𝑇 𝑝.0 -𝑇 𝑝.𝑒𝑞 ) 𝑅 𝑝𝑐𝑚/𝑝 𝐿 𝑓	+ ( 𝑇 𝑝.𝑒𝑞 -𝑇 𝑝𝑐𝑚 𝑅 𝑝𝑐𝑚/𝑝 𝐿 𝑓 +	𝑇 𝑒 -𝑇 𝑝𝑐𝑚 𝑅 𝑒/𝑝𝑐𝑚 𝐿 𝑓 ) 𝑡	(14)

Table 1 :

 1 Input and output parameters of each model

	Model	Input parameters from	Input parameters from	Outputs
		box/PCM/product/air	measurements or
		characteristics	correlations
	Lumped	Box:	R p/pcm = heat transfer	Mean product
	model	 Dimension  Insulation (by	resistance between PCM and product.	temperature evolution: T p.mean (t)
		conductivity of the	(determined from product	Melted PCM
		material)	temperature at equilibrium	evolution: m melted.pcm
			by constantly changing	(t)
		PCM:	PCM as mentioned in
		 Latent heat of melting  Melting temperature  Specific heat  Mass  Dimensions	section 2.1.4)
		Product:
		 Dimensions  Mass  Specific heat
	Zonal	Same parameters required	 h w and h p = internal convective heat transfer	Temperature
	model	for the lumped model cited	coefficient	evolution for air,
		above.		product core and
			(determined from local	shell in 4 zones
		Additional parameters:	temperature	(top/bottom;
		 Box wall emissivity  Thermal conductivity of the product  Product emissivity  Specific heat capacity of the air	measurements as described in [26] or free convection correlations)  = mass flow rate of 𝑚 𝑎	left/right):  T p.i (t)  T a.i (t) and T' a.i (t) Melted PCM evolution: m melted.pcm
			air	(t)
			(determined from
			developed relation with
			heat transfer coefficient as
			shown in [27])

Table 2 :

 2 Model applicability to answer technical questions

	Technical questions	Lumped	Zonal	CFD
		model	model	model
	How does the product temperature change with time under varying ambient temperatures as in a real supply	 *		
	chain?
	What is the required PCM mass to maintain the average product temperature under the recommended values	 *		
	during a given duration?
	How does the box insulation affect the product temperature?	 *		
	How do PCM melting temperature and its latent heat impact the product temperature?	 *		
	How do the mass and the thermophysical properties of the product impact its temperature?	 *		

Table 3 :

 3 Input parameters for numerical studies

	Parameter	Value	Unit
	Reference Box
	Internal dimensions of the box (L x W x H)	0.5 x 0.3 x 0.31	m 3
	Wall thickness	0.04	m
	Heat transmission coefficient of box insulation	0.58	W‧m -2 ‧K -
	(K)		1
	PCM (ice) on a side wall
	PCM dimensions (L x W x H)	0.5 x 0.05 x 0.3	m 3
	PCM mass	2	kg
	Specific heat (solid state)	2070	J‧kg -1 ‧K -1
	Specific heat (liquid state)	4217	J‧kg -1 ‧K -1
	Latent heat of fusion	333,700	J‧kg -1
	Melting temperature range	-0.2 to 0.2	°C
		(0°C for the lumped
		model)
	Product (Tylose)
	Product dimensions (L x W x H)	0.4 x 0.2 x 0.2	m 3
	Product mass	16	kg
	Specific heat	3372	J‧kg -1 ‧K -1

  ) 𝑚 𝑝𝑐𝑚 (𝐿 𝑓 + 𝐶 𝑝𝑐𝑚.𝑠𝑜𝑙 (𝑇 𝑚 -𝑇 𝑝𝑐𝑚.0 )) + 𝑚 𝑝 𝐶 𝑝.𝑝 (𝑇 𝑝.𝑒𝑞 -𝑇 𝑝.0 )

	(1 -𝑒 1 + 𝑅 𝑝𝑐𝑚/𝑝 𝑡 𝑚𝑒𝑙𝑡 -𝜏 ) 𝑅 𝑒/𝑝	≈ ( 𝑇 𝑝.𝑒𝑞 -𝑇 𝑚 𝑅 𝑝𝑐𝑚/𝑝 +	𝑇 𝑒 -𝑇 𝑚 𝑅 𝑒/𝑝𝑐𝑚 ) 𝑡 𝑚𝑒𝑙𝑡

Table 4 :

 4 Advantages and limitations of each model

	Criteria	Lumped model	Zonal model	CFD model
	Airflow and heat	Not applicable	+	-
	transfer assumptions
			(Posteriori)	(Priori)
	Number of input	11	17	16
	parameters
	Preliminary	+	±	-
	experiments

Acknowledgement

The authors would like to thank the Royal Thai Government Scholarship and Chulalongkorn University, Bangkok, Thailand, for T. Leungtongkum's PhD scholarship. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Appendix 1: Governing equations of the zonal model

Estimation of air temperatures from product shell and wall temperatures

The air at position I in Fig. 4 exchanges heat with PCM surface (wall 1) and its temperature shifts from T' a.1 to T a. [START_REF] Zhao | Application and research progress of cold storage technology in cold chain transportation and distribution[END_REF] . The heat balance between the adjacent air and wall 1 can be shown as follows.

𝑚 𝑎 𝐶 𝑝.𝑎 𝑑𝑇 𝑎 = ℎ 𝑤 (𝑇 𝑤.1 -𝑇 𝑎 )𝑑𝐴 By integration: (A.1)

where is the mass flow rate of air [kg‧s -1 ] 𝑚 𝑎 C p.a is the specific heat capacity of air [J‧kg -1 ‧K -1 ] h w is the heat transfer coefficient between the internal air and the internal wall

is the dimensionless convective heat transfer coefficient between the internal air 𝛼 𝑤.1 and the internal wall 1 [-] When air exchanges heat with the product shell, the same approach is also applied, for instance from T a.2 to T' a.2 . Thus, eight linear equations involving the eight air temperatures are obtained allowing air temperatures estimation from the product shell and box wall temperatures. Full details regarding the estimation of the convective heat transfer coefficient and air mass flow rate can be found in Leungtongkum et.al. [START_REF] Tiari | Numerical study of a latent heat thermal energy storage system enhanced by varying fin configurations[END_REF].

Transient evolution of product, walls, and PCM

Eq. A.2 is the unsteady heat balance equation for the shell of the product block 1 (T sh.1 ). where MC p.sh is the thermal inertia of the product shell [J‧K -1 ] R sh.c is the heat transfer resistance from the shell to the core of the product [K‧W -1 ] is the radiative heat exchange between the surface 𝑞 𝑟.𝑠1.𝑤1 = 𝜀 𝑤.1 σ(𝑇 𝑠.1

4 -𝑇 𝑤.1

)𝐶𝐿 of product block 1 and wall 1 [W]

is the radiative heat exchange between the surface of 𝑞 𝑟.𝑠1.𝑤4 = 𝜀 𝑤.1 σ(𝑇 𝑠.1 4 -𝑇 𝑤.4

)𝐶𝐿 product block 1 and wall [W]

The same approach was applied to the 12 solid zones: the shells of the product blocks, the cores of the product blocks and box internal walls.

PCM temperature and ice fraction evolution can be written as for the lumped model (Equations 10 to 12) except that the PCM exchanges heat with wall 1 which is in fact a PCM-plate wall, not directly with the product.

Appendix 2: Material's properties

Highlight

 Lumped, zonal and CFD models were developed for transport in insulated boxes with PCM.  The demonstration of how to use these models to deal with practical issues are shown.  The lumped model is suitable where the temperature heterogeneity is not a concern.  The zonal model, more complex, provides temperature evolution in different zones.  The CFD model, the most complex, provides temperature and air velocity fields.

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: