Assessment of insulating package performance by mathematical modelling

A mathematical model has been developed in the present work to describe the temperature change in a typical insulated shipping container as a function of time. The model was created by combining steady state and transient models in a 2D geometry of a typical shipping container and was subsequently validated by an ice melt test and comparison of temperature change obtained from the model and experimental measurement. An excellent agreement was obtained between the computational model developed in this work and experimental results. In addition, a parametric study was also carried out to investigate various factors in controlling the insulation performance of the packaging. It was found that the model has capability of evaluating the effect of a wide range of packaging design parameters such as thermal conductivity, surface emissivity, packaging geometry, and sounding temperature.

, Thickness of air gaps without reflective foil [mm]; ε in , Emissivity of package outer surface; ε out , Emissivity of package inner surface; ε product , Emissivity of transported product postal service as a mean of goods transportation, which is often much cheaper than active cold chain. The performance of a passively chilled package is influenced by multiple factors including packaging shape, insulation materials, quantity, and type of PCM used.
A number of studies have been undertaken to evaluate the performance of passively chilled packages. [3][4][5][6][7][8] Burgess derived equations describing the thermal resistance (R-value) of a package based on packing material, thickness of packing walls, geometry of the container, and other design features. 9 The equations were obtained after conducting multiple ice melt tests and fitting the formulas to match the observed results through linear regression. The model was then used to predict the shipment time for only selected packages and was compared with physical measurement. The difference was found within 20%. In addition, this approach presented limited utility as multiple tests were required to manipulate packaging models in order to create a database. 9 8,11 Both incorporated detailed descriptions of the package and shipment products into their models. Despite the fact that the models neglected the heat transfer through radiation or convection, they were still able to achieve excellent correlation between experimental and simulation results. This approach showed great potential for insulated package design, although it also requires the access to expensive software and significant computational resources. This could significantly limit its accessibility for ordinary passive cold chain users. 8,11 The experimental approach has also been commonly used in the past in order to design and test the performance of the shipping con-wide range of factors affecting the outcome of the measurements make it extremely difficult to compare the results from different studies. 12  In the present work, we have developed a transient mathematical model to assess insulation performance of passive insulated packages.
In addition, the effect of various insulating materials on package insulating performance was also investigated through ice melt tests and temperature profiling. The results obtained from the model and experimental measurement were subsequently compared to validate the model. Finally, the parametric study was carried out to investigate sensitivity of the model.

| MATHEMATICAL MODELLING
A passively insulated package typically includes a container with single or multiple layers of insulating materials to for part of packaging wall.
As a result, such package may be viewed as a heat sink with external energy attempting to travel through the container wall. It is reasonable to assume that conduction will dominate overall heat transfer in a per- The mathematical model developed in this work is based on thermal energy balance, which describes basic thermal processes such as heat transfer, conversion of heat to different forms of energy as well tainers. 12 However, large variety of insulated container design and a ronment inside the package. Such system opens up the possibility of as heat accumulation. It is assumed that the temperature inside the package must be known, and it is exactly the same as the temperature of the transported goods. Starting with the balance equation, we have Assuming that the direction of net energy flow is from the external environment towards the inner space of the container as shown in Figure 1, Equation (1) will be reduced to Assuming the temperature change of the coolant and the transported product will change at the same rate with time, the rate of energy change might be described as 4 Moreover, total heat flow entering the package can be written as The relation between HPR and the thermal resistance of the package can be described as HPR may be evaluated using following equation 5: In order to obtain RW, the following equation was used: 0:217*ta 3 5:918 þ ta 3 : Moreover, heat transfer coefficients of outside and inside of the package were calculated using Substituting Equations 3 and 4 into Equation 2 leads to below ordinary linear differential equation: The analytical solution to Equation 10 is in a form of Equation 11: When phase changing materials (eg, ice) are present inside the container, the latent heat of fusion of coolant must be taken into account.
For instance, our model assumes that the energy transferred into the container is totally consumed to sustain the process of ice melting at a temperature of 0°C. The beginning of phase change is defined at the ice melting temperature, and the process is assumed to end when

| Ice melt test
Ice was placed in the room at a controlled temperature (20°C) until reaching melting temperature. Five hundred grams of ice was placed inside the container whose dimension are presented in Figure 2. The sealed container was left in a temperature-controlled environment at 20°C. After 10 hours, the remaining ice was taken out of the container, and the amount of water produced in the process of ice melting was weighted. Assuming that the energy transferred into the box was entirely consumed to melt the ice, HPR may be calculated by using following equation: The ice melt test can help to evaluate whether the model devel-

| Measurement of the maximum insulation time
This measurement was conducted in accordance with ASTM D3103.
The outside conditions were kept constant with a temperature of 15°C (±1°C) and a relative humidity of 54% (±2%). The thermocouple was a Lascar Electronics EasyLog EL-USB-1 with measurement accuracy of ±0.5°C. Two 400-g ice packs were sourced in order to keep the volume and shape of the ice constant throughout all experiments.
The ice packs were placed in a freezer for 48 hours before being posi-     Nevertheless, there are still a number of areas presenting some distinctive differences. Firstly, the rate of temperature increase is steeper for calculated values than experimental results during the postmelting stage, and such temperature increase also occurs earlier in the experiment than in the modelling. This is mainly due to temperature discrepancy between the measurement and the modelling. Temperature measurement is taken from the surface of the ice packs, which have a temperature gradient from the surface to the bulk in nonequilibrium state. The model, however, assumes a uniform temperature throughout the ice packs and therefore results in a delayed phase change followed by a rapid increase in temperature after completion of the phase change. One should expect that such discrepancy can be addressed if the model takes into account temperature gradient in the coolant itself. This improvement will require greater computational resources due to the increased calculation difficulty. Since the difference in the MIT is relatively small as shown in Figure 4, the above assumption is kept for the remaining investigation in this work.

| Heat penetration rate
Although only a single size of container was used in the present study, previous work had showed that the HPR could be calculated with reasonably good accuracy for a wide range of containers with external sizes varying between 23.5 and 109.2 cm. 5 As the HPR is a key parameter for the calculation of the MIT, it is reasonable to assume the model derived in this work should also be applicable to packages with different dimensions.  Table 1. In addition, Figure 5 shows packages with four different insulating liners, namely low density polyethylene, polystyrene foam, PIR board, and silica aerogel blanket.
It can be found in Table 1 that the aerogel-lined package gives a MIT nearly 2.5 times greater than that obtained from the package lined with polyethylene foam, despite the fact that the former is 5 mm thinner than the latter. This is mainly due to the extremely low thermal conductivity in the aerogel. In addition, the results for EPS-lined containers in Table 1 also clearly indicate the benefit of incorporating reflective inner surface, which gave rise to 10% increase in MIT. Furthermore, the results in Table 1 have also revealed that the model is reasonably accurate to estimate the MIT with an average difference of 10% relative to the measured values. This accuracy is   It was previously proven that both coefficients are relatively constant within the expected temperature range. 5 Similarly, the thermal conductivity of the insulation materials is expected to rise in accordance with an increasing temperature, but thermal conductivity as a function of temperature was not introduced within the mathematical model. A decision was undertaken following comparative study presenting minor differences between constant and varying thermal conductivity for the temperature range we expect the package to be exposed to.

| PARAMETRIC STUDY
A parametric study was carried out to investigate the influence of design parameters of a typical insulated package on its insulation performance, which was characterised by the MIT. These parameters included surface emissivity, thermal conductivity of the material, mass of coolant, thickness of insulation layer, and external temperature.
Their relation to the MIT was expressed by rearranging Equations 11 to 13 as below.
The package considered in the parametric study was a cardboard box that contained 1-kg ice packs with initial temperature at −20°C.
The external temperature was set to be at 20°C, and as defined earlier

| Effect of surface emissivity
As shown in Figure 7, surface emissivity can significantly affect the insulation performance of a shipping container. By merely incorporating a low surface emissivity surface on the outside of the container, the MIT was increased by 12%. More increases were observed at surface emissivity below 0.8 when the inner surface was equipped with reflective foil. The maximum improvement in the MIT was found to be almost 40%. It can be also seen that in both cases the MIT increases exponentially as surface emissivity reduces and graphs converge at the surface emissivity value of 0.83.
Differences in surface emissivity are mainly attributed to the radiation part of the energy transfer, and, as described by Stefan-Boltzmann equation, decreasing surface emissivity will lead to a decrease in the energy flux through the container walls. 18 However, the impact of such a layer will differ at different positions within the container. In the case of the inner wall, the surface with low emissivity creates conditions in which elevated surface temperature causes only a small value of radiant thermal energy being radiated towards the shipping product. On the other hand, the application of a low emissivity surface at the outside of the shipping box will also increase the reflectivity of the container. As a result, a significant part of energy radiated towards the package from the outside environment will be reflected and will not contribute towards the increase in temperature within package. Furthermore, it can be found that much better effects upon MIT are noticed when low surface emissivity material is applied to the inside surface of the container. This happens, due to the thickness of the wall which reduces the inner surface compared with the outer one. As a result, even though the outside surface area and property controls the amount of energy being absorbed by the container, it is the inside surface which limits the heat passage towards the product inside the container and enhances the effects the surface emissivity has on MIT.  Figure 9. There appears to be a drastic increase in the MIT after thermal conductivity falls below a critical point as seen in Figure 9.

| Effect of insulation material properties
In order to identify this critical point, the linear regression method was applied to the linear regions before and after nonlinear region for each insulation thickness. A thermal conductivity was then obtained at the intersection point of both lines. For all analysed thicknesses, the critical value falls in the vicinity of 14 mW/mK, although it is susceptible to changes if significant alterations of design parameters are introduced. Nevertheless, it is clear from Figure 9 that an attempt at increasing packaging MIT will be much more effectively achieved when the thermal conductivity of a lining material is below this critical point. In practice, such low thermal conductivity can be obtained by using materials such as aerogels. Furthermore, the effect achieved by decreasing thermal conductivity also scales with the thickness of the insulating liner. On the other hand, increasing thermal conductivity leads to convergence of the MIT approaching 20 hours as seen in Figure 9. This is expected since for a given wall thickness its insulation effect will diminish as the wall becomes more and more heat

| Effect of mass of coolant
A common practice in transportation of perishable goods is the addition of a phase-change material (eg, ice) inside the package in order to elongate possible shipment time by sustaining low temperature inside a container. The model derived in this work was used to investigate the effect of the quantity of ice pack on the MIT. Figure 11 shows that the addition of PCM can elongate the MIT linearly as indicated by Equation 13. This is caused by the addition of heat sinks inside the package, which effectively absorb the entering energy to transform coolant phase at the phase-change temperature before changing the temperature of the product. It should be also noted that the MIT can be adjusted by the amount of coolant as well as type of coolant. The model can help to select the optimum amount of cooling substances to keep the products refrigerated whilst minimising package weight and associated transportation costs. Figure 12 shows the modelling results for the MIT as a function of external temperatures. As expected, the increase of outside temperature gives a monotonic decrease in the MIT since increasing the temperature difference across the packaging walls results in higher heat flux entering the package. It is worth mentioning that the nonlinear correlation between the external temperature and the MIT is implied by Equation 13. The external temperature affects all modes of heat transfer, although in previous study, it was shown that radiation and convection were not significantly affected by a temperature change as selected in this investigation. 6 The results in Figure 12 also shows that the model can help to evaluate the insulation capacity of a package experiencing a wide range of temperature change during transportation.

| CONCLUSION
The present work presents a mathematical model capable of evaluating the insulation performance of passive insulated packages during their transportation of perishable goods. The model showed a good agreement with experimental measurement giving an average of 10% relative difference for the HPR and the MIT. In addition, a parametric study of the model was also carried out to investigate the influence of a variety of factors on the MIT of a typical delivery package.
These include surface emissivity, thermal conductivity of insulation, insulation thickness, mass of coolant, and external temperature. The modelling results for all these parameters followed closely to the  Modelling results for the effect of external temperature on the maximum insulation time expectation, and some of these results were also validated by comparing with those obtained from the experimental measurement. More importantly, it was clearly demonstrated that the model possessed the ability to provide a much wider spectrum for illustrating the effect of most key parameters for insulated packaging design. Furthermore, this analytical model can be implemented without engaging any numerical approach and resourceful computational tools. Thus, this mathematical model proved to be an excellent tool for design of insulated packaging as well as gaining a good understanding of temperature-time behaviour for a given insulated package.