Experimental and numerical study of insulation packages containing dry ice pellets

Shippers of healthcare and food products, in most cases, have regulatory responsibility to ensure quality by maintaining product’s temperature during transportation. Dry ice is commonly used as cold source inside the insulated cold-storage packages which are often transported via air freights. But dry ice sublimation can lead to asphyxiation of crew members or passengers due to which, several regulatory bodies have recommended safe limits of dry ice aboard aircraft mainly based on rule-of-thumb. Prior to establishing the limits, understanding dry ice sublimation is of paramount importance. Therefore the current research, experimentally and numerically explores dry ice sublimation packed in commonly encountered containers made of polystyrene foam. A three dimensional quasi-steady model is developed to predict variation of dry ice mass inside the insulation packages. The model is validated against mass and temperature measurements performed on two types of insulation packages made of different material and geometry. The model is able to fairly predict the reduction in dry ice mass over time and the end-of-sublimation time. Both experimental and numerical results show that the sublimation rate of dry ice can be reduced by covering the inner walls of insulation package with a reflective layer like aluminized mylar foil. The good agreement between model and experimental data of dry ice mass variation makes it possible for the shippers to use the modeling approach in estimating sublimation rates of dry ice.


Introduction
Since the beginning of twentieth century, the developments in transport of temperature controlled goods has given rise to a cold chain ecosystem which contributes significantly to the global economy. Mainly the pharmaceuticals, bio-pharmaceuticals and perishable food products require controlled temperature conditions to maintain quality during their transport. The temperature control spectrum of these products range from ambient (15 • C to 25 • C), refrigerated (2 • C to 8 • C), frozen (-80 to − 20 • C), to cryogenic (< − 150 • C) temperatures. As a result of the regulatory requirements set by the governing bodies like FDA, ICH and ATP, the distributors are required to transport the temperature sensitive product without braking the cold chain. Air freighting is increasingly being utilized for perishables, for instance, Emirates SkyCargo transported over 400,000 tonnes of perishables in the year 2018 [1]. Temperature sensitive shipments of air cargo often utilizes dry ice as a cold source due to its high sublimation enthalpy (571*10 3 J kg − 1 ) which helps in maintaining the desired temperature for longer period of time. The products along with dry ice are commonly transported into insulated container consisting of expanded polystyrene (EPS) foam box and an outer cardboard covering, see Fig. 1. The main purpose of the cardboard covering is to carry safety labels and provide protection to the foam box.
The gaseous carbon dioxide (CO 2 ) formed during the sublimation of dry ice is a potential risk for airline operators as high concentration of CO 2 gas can lead to asphyxiation of passengers and crew members. Dejohn et al. [2] reported an incident of respiratory impairment of the crew members present in the Douglas DC-8 aircraft while taxiing for takeoff in the year 1998. The impairment was caused by the accumulation of CO 2 fumes in the cockpit which originated from dry ice present in the cargo compartment of DC-8 aircraft. Out of the ten cases of sudden in-flight incapacitation of Air France pilots between 1968 and 1988 reported by Martin-Saint-Laurent et al. [3], the cause of one case was CO 2 intoxication due to inadequate packaging of container refrigerated using dry ice. Although none of them resulted into an accident, dry ice sublimation can be dangerous in enclosed spaces of an aircraft where ventilation rates are low. Therefore, dry ice is classified as hazardous material under Code of Federal Regulations (CFR) -49 CFR 172.101 [4], and the concentration limit of 0.5% by volume for gaseous CO 2 is specified under 14 CFR 25 [5].
In the past, several recommendations, mainly based on experiments are proposed to establish quantity of dry ice that can be safely carried on an aircraft. Federal Aviation Administration (FAA), in its advisory circular AC 103-4 [6] of 1974 listed some precautionary measures for carrying dry ice aboard aircraft. The document suggested to use the sublimation rate of 1% per hour based on the experiments performed by Pan American airlines in 1963 using a 45.3 kg (100 lb) block of dry ice [7]. However, in 2009, the circular AC 103-4 was subsided by a new circular AC 91-76A [8] which recommended sublimation rates for commonly encountered dry ice packages. It recommended sublimation rate of 2% per hour based on the experiments performed by Caldwell et al. [7]. They placed Thermosafe insulation box (0.27 m x 0.22 m x 0.18 m), packed with 2.2 kg (5 lb) of dry ice pellets, in an altitude chamber depressurized to an altitude of 8000 ft. The sublimation rate of 2% per hour was determined by measuring the weight loss of dry ice pellets over the period of six hours. In 2013, Murphy and McSweeney [9] assessed dry ice limits for aircraft by analyzing sublimation of dry ice placed inside two types of commonly used containers, namely, unit load devices and EPS packages. They measured weight and temperature data for the EPS package (0.29 m x 0.24 m x 0.2 m) containing dry ice pellets. The measurement data was recorded manually at different intervals of time and only weight data was presented in the report. The heat transfer rate into the system was estimated from a one dimensional steady state heat transfer model using Eq. (1). Where U is the overall heat transfer coefficient, A is the total surface area of the EPS box, and ΔT (T ∞ − T s ) is the temperature difference between dry ice and ambient temperature. The overall heat transfer coefficient takes into account several convective and conductive resistances as shown in Eq. (2). The convective resistances are considered at outside wall ( 1 ho ), inside wall ( 1 hi ), dry ice surface ( 1 hs ), and at the gap between cardboard and EPS box ( 1 hg ). The conductive thermal resistances are accounted for EPS foam box ( x f k f ) and cardboard box ( xc kc ). The assumptions made for evaluating the heat transfer coefficients and thermal conductivity values of Eq. (2) can be found in their report [9].
The energy content of sublimating dry ice is given by Eq. (3) where ṁ and l are the sublimation rate and specific sublimation enthalpy of dry ice. Finally, the sublimation rate of dry ice is calculated from heat balance by equating Eq. (1)  However, there are some limitations associated with the methodology of the studies in the past which resulted into inaccurate prediction of sublimation rate. Several studies determined the sublimation rate experimentally which in most of the cases was evaluated by measuring the initial, final weight of the package and time period of the experiment. The variation of mass of dry ice inside the packages until all the dry ice has sublimated (end-of-sublimation time) has not been reported. Currently, the dry ice limits on an aircraft are established based on simple rule of thumb and simplified heat transfer models. The state-ofthe art models are often based on steady state one dimensional heat transfer as described in Eqs. (1)-(3). Murphy and McSweeney [9] reported several uncertainties with the steady state one dimensional model. Namely, the model does not consider reduction in level of dry ice. Moreover, it is assumed that the gas between the foam box and the dry ice is at sublimation temperature of dry ice. These assumptions lead to a constant temperautre at the inner walls of the EPS package equal to the dry ice sublimation temperature. As a result the heat transfer rate (Eq. (1)) and in turn the sublimation rate (Eq. (3)) are overpredicted by the state-of-the art models.
In this paper, our objective is to study the sublimation phenomenon inside an insulation package to preclude hazard due to dry ice sublimation aboard aircraft. To develop understanding of this phenomenon, we performed mass and temperature measurements until all the dry ice is sublimated. Further, we propose a modeling approach that can be useful for the aircraft carriers/shippers in predicting the mass variation of dry ice and estimating the end-of-sublimation time. We also seek to shed light on the manner how radiation affects the sublimation of dry ice by covering the inner walls with highly reflective aluminumized mylar foil.

Materials and methods
The examination performed in this study is based on an extensive comparison with experiments. Tests were performed using two packages, a cuboid box made of expanded polystyrene foam (EPS) and a cylinder package made of graphite filled expanded polystyrene foam (GEPS), see Fig. 2. The primary interest of studying dry ice sublimation inside two foam packages of different material and geometry is to get a robust validation of the modeling approach with the experimental results. This would allow us to conclude if the modeling approach presented in Section 3 is valid irrespective of insulation package geometry and material. In engineering applications like insulaion for walls of the building, GEPS is increasingly becoming popular due to its low thermal conductivity. Lakatos et al. [10] measured thermal conductvity of five different samples of GEPS. The thermal conductivity value measured at 293 K for a sample with density of 14.35 kg m − 3 is approximately equal Koru [11] reported the thermal conductivity of EPS approximately equal to 0.0381 W m − 1 K − 1 . Since there is only a marginal difference (10 %) in the thermal conductivity values of EPS and GEPS, the vlaues considered in this work are similar for both materials to predict dry ice mass vairiation and end-of-sublimation time up to engineering accuracy. The value and the rationale for the thremal conductivity chosen for EPS and GEPS in this work are stated in Section 4.2.
The cuboid box shown in Fig. 2 (a) has an external dimension of 0.265 x 0.215 x 0.190 m and a wall thickness of 0.04 m. The outer diameter, height and thickness of the cylinder package shown in Fig. 2 (b) are 0.148 m, 0.2 m and 0.0265 m, respectively. The cuboid box is from a commercial supplier and comes with a fitting cardboard box of 3 mm thick walls as shown on the right side of Fig. 2 (a). The cylinder does not have a fitting cardboard cover as it is made in-house from a slab of graphite filled polystyrene foam.
The experiments were designed with an objective to validate the numerical solutions, namely, the variation of dry ice mass and temperature inside the packages with time. For this purpose, the packages were completely filled with dry ice pellets and were placed on a weighing scale to measure the mass of remaining dry ice pellets as the CO 2 gas escaped the box after sublimation. A dry ice pellet is in the shape of a small cylinder having approximately 9 mm diameter. The packages were filled with many such pellets until the packages were full. The temperatures inside the packages were measured using E-type thermocouples whose exact location is depicted in the schematic shown in Fig. 3. The experimental procedure is similar for both packages and consists of the following steps: (1) Thermocouple tip is attached at their locations using a 1 x 1 cm aluminium tape; (2) Empty package is placed on the weighing scale and dry ice pellets are filled in the package; (3) The package is covered with a lid and the weighing scale is tared to zero; (4) A LabView program measures the mass and temperature every 3 s and 10 s, respectively; (5) Measurement is stopped when the temperature readings of the thermocouples are close to the ambient temperature. For the EPS box, the CO 2 gas generated during the experiment escapes the EPS package through the thin gap between the lid and its body. Whereas a small hole is drilled on the lid of the GEPS package, where thermocouple 3 is placed, for the gas to escape.
Additional set of experiments are carried out for the cuboid package with its inner walls covered with a single layer of aluminized mylar foil. The purpose of these experiments was to investigate the effect of foam surface reflectivity on the sublimation of dry ice pellets. Fig. 4 (a) shows an experimental arrangement for two cuboid packages, where one of the two packages contains mylar foil on its inner walls as shown in Fig. 4 (b). The methodology similar to the one discussed in the previous paragraph is followed and the measurements are started simultaneously for the two cuboid packages.

Model
The theoretical analysis is based on quasi-steady approximation. This approximation is valid for small values of Stefan number. Dandforth [12] showed in his doctoral dissertation that for Stefan numbers less than 0.2, errors less than 10 % can be expected in predicting the ice depth of freezing lakes. Moreover, the author also indicated possiblity to apply the quasi-steady approximation to other moving boundary problems involving high latent heat of phase change as it will lead to smaller values of Stefan number.
For a slab of dry ice sublimating in one dimension, the expression for  Stefan number is shown in Eq. (4). Where c p is the specific heat capacity of carbon dioxide gas; T ∞ and T s are respectively the ambient temperature and sublimation temperature; and l is the latent heat energy per unit mass added during sublimation. The numerator of Eq. (4) represents the sensible heat added to a unit mass of gas to raise its temperature from sublimation temperature to the ambient temperature. Using the values: .65 K and l = 571*10 3 J kg − 1 , we get Ste = 0.13. Since the latent heat of sublimation of dry ice is very large, the small Stefan number implies that dry ice interface moves very slowly and the temperature distribution at each moment can be obtained by solving a steady state heat transfer model.
The quasi-steady assumption for the sublimation of dry ice placed inside insulation package can be made based on an estimate of its sublimation rate. This can be done using the heat balance shown in Eq. (5) is calculated to be 18.24 W. Consequently, the sublimation rate ṁ = 0.03 gs − 1 is calculated. Since the sublimation rate is small, it is assumed that quasi-steady assumption is applicable in this case.

Model description
In practice, the heat transfer rate into the insulation package will reduce with time as opposed to the constant value calculated in the above estimate. Therefore, to predict the reduction in heat transfer rate and dry ice mass, the approach used in this paper includes solving a quasi-steady three dimensional numerical heat conduction and radiation model using the commercial software package COMSOL 5.4.
The modeling approach is similar for both geometries and it is described for the cuboid box, made of EPS material. Fig. 5 shows a two   dimensional representation (front view cut section) of the three dimensional geometry used in the numerical model. The geometry constitutes of four domains, namely, EPS, dry ice, CO 2 and cardboard. The dry ice used in the experiments is in the pellet form but a continuous solid domain of dry ice is assumed in the model shown as shaded area in Fig. 6. The cold CO 2 gas formed after sublimation occupies the head space above dry ice until slight over pressure is created in the head space for it to escape through the gap between the lid and the body of the EPS box. The dynamic process of gas accumulating and leaving the head space is not simulated in the model. It is assumed that head space between dry ice and lid is occupied by freshly formed CO 2 gas and is stagnant at all instance of time. It is also important to note that thermal resistances of thin air layers between cardboard box and EPS, and other different contact resistances are considered negligible as compared to the large thermal resistance of the thick EPS foam box.
The energy transfer due to radiation occurs between the boundaries that bound CO 2 gas domain which is schematically shown with the help of arrows in Fig. 5. The CO 2 gas does not participate in radiation and is transparent to the radiative energy transfer between boundaries which are assumed as diffuse-gray surfaces.
In radiative heat transfer between two diffuse-gray surfaces (i,j), the heat transfer that leaves the surface i is given by Eq. (6), where ∊ i and A i are the emissivity and the area of surface i, respectively. In Eq. (6), J i is the surface radiosity which is the sum of reflected portion of the incoming radiation (ρG) and the emitted radiation (∊ i E b,i ). According to the Stefan-Boltzmann law, E b,i depends on the fourth power of the temperature as shown in Eq. (7), where σ = 5.67*10 − 8 W m − 2 K − 4 is the Stefan Boltzmann constant.
Similarly, the equation for heat transfer that leaves the surface j can be obtained by replacing the subscript i with j. The net heat transfer rate between two surfaces i and j can be calculated using Eq. (8), the detailed derivation for which can be found in the book of Heat transfer by Bejan [13]. Where F ij is called the view factor which represents the fraction of radiation from surface i being intercepted by surface j. In COMSOL, the view factors are automatically calculated from the mesh for each subsurfaces of the mesh as noted in [14]15.
The radiative heat transfer presented above for two surfaces can be generalized to an enclosure formed by n diffuse gray surfaces. The net heat transfer rate by any surface i is given by Eq. (9) where J i is the radiosity of surface i depending on the geometry and radiative properties of all the surfaces of the enclosure as shown in Eq. (10) where i = 1, 2, …, n.
In COMSOL, the heat transfer rate by radiation is coupled to heat conduction through a source term added to the heat equation [14]. The heat equation is derived from first law of thermodynamics (Eq. (11)) which states that the variations in internal energy and kinetic energy are caused by mechanical power or by the exchanged heat rate. Where the exchanged heat rates can account for thermal conduction (q), radiation (q r ) and other heat sources like Joule's heating and exothermic chemical reactions (Q) as shown in Eq. (12). The vector n is normal to the boundary ∂Ω of the domain Ω.
In the present work, we perform steady state simulations where a simplified heat Eq. (13) is solved in COMSOL. The conductive heat flux, q, is obtained using Fourier's law and the radiative heat flux, q r , is calculated based on radiosity method discussed above.

Methodology
The model uses an iterative approach in which the dry ice level is changed in every subsequent iterative step. Initially, at time = t 0 , the dry ice level is close to the lid and it corresponds to the situation of a fully filled package, see Fig. 6. A small gap of 1 mm in height is assumed for the CO 2 gas in the first time step. With this geometry, a steady state heat diffusion equation is solved with radiative heat transfer occurring only in the CO 2 gas domain.
As a part of post-processing, the total heat transfer rate (Q 0 ) from ambient to the package is calculated by surface integration of normal heat flux on all the six faces of the cardboard box. In the first iterative/ "time" step, t 1 − t 0 , the previously calculated heat transfer rate is used to evaluate the mass of dry ice sublimated (m d ) using Eq. (14). Consequently, the height by which the dry ice interface moves down is calculated using Eq. (15). It is assumed that the dry ice interface moves down uniformly, or in other words, the interface is flat and has a constant surface area. In Eq. (15), L i and W i represent the inner dimensions of the EPS package and, ρ d is the effective density of the dry ice (850-900 kg m − 3 ) which is experimentally obtained by dividing the initial mass of dry ice pellets by the inner volume of the EPS package. After calculating the new height of dry ice, the model geometry is changed and simulated to obtain a new value of heat transfer rate (Q 1 ). The procedure is repeated in subsequent steps using MATLAB script for every 1 mm reduction in dry ice interface height until it becomes equal to zero as shown in Fig. 6. The time at which this condition is satisfied is considered as the end-of-sublimation time (t f ).

Boundary conditions
As dry ice sublimates and its level reduces, the inner walls of EPS box are only partially in contact with dry ice. Therefore, a temperature boundary condition equal to the sublimation temperature (T s = 194.65 K) is specified on the part of boundaries which are in contact with dry ice at a given instance of time corresponding to the dry ice level, See Fig. 7. This temperature boundary condition is also specified at the dry ice interface. The heat transfer from the ambient to the package is taken into account by specifying a heat transfer coefficient between cardboard walls and the ambient. In the work of Kozak et al. [16] on foam packages containing phase changing materials, an effective heat transfer coefficient of 10 W m − 2 K − 1 , taking natural convection and radiation into account, is calculated.
Alternatively, the heat flux boundary condition on the wall can be replaced by a temperature boundary condition (T ∞ ) without having a significant effect on the end-of-sublimation time of dry ice. Consider a one dimensional heat transfer case shown in Fig. 8 (left) where the thin cardboard and thick insulation are considered as one slab with effective thermal conductivity k f and total thickness of L. The heat flux from the ambient, q ′′ sublimates the dry ice which is at sublimation temperature, T s . The energy balance at the external surface of the slab is shown in Eq. (16) and Eq. (17). The equation can be rearranged to introduce a nondimensional Biot number ( hL k f ) which signifies the ratio of heat transfer resistances inside the insulation and at the surface of the insulation. The values of wall temperature, calculated using Eq. (18), for different Biot number is shown in Fig. 8 (right). It can be seen that the wall temperature approaches the ambient temperature as the Biot number increases. The Biot number for the cuboid box considered in this work is 16 and the corresponding wall temperature is evaluated to be 287.35 K as marked in Fig. 8 (right). The measured values of the ambient temperature (T ∞ ) and the wall temperature (TC 1) for both insulation packages are shown in Fig. 12. It can be seen that they differ by few degrees from each other and consequently replacing heat flux boundary condition by a temperature boundary condition equal to T ∞ will not have a drastic effect on the heat transfer rate or the end-of-sublimation time.

Verification of numerical results
The numerical model is first verified by presenting the convergence behavior in Fig. 9. For each geometry, the numerical solution of total heat transfer rate (Q), at three different dry ice levels, is shown to vary with number of mesh elements in Fig. 9. Surface area integration of normal heat flux on the exterior walls of geometry is performed to obtain Q. It can be seen that the value of Q reduces as the number of mesh elements are increased. For instance, Q calculated for the cuboid geometry when the dry ice level is at H = 100 mm, varies with mesh elements as: 9.30 W, 8

Sublimation rate
The numerical model is validated against the experiments performed on two insulation foam boxes having different geometry and material. The two boxes are in the shape of a cuboid and cylinder made of EPS and GEPS material, respectively. Fig. 10 shows a comparison between the measured and predicted dry ice mass over time for the two boxes. The starting mass for cuboid and cylinder boxes are different and correspond to their capacity of holding dry ice pellets. The predicted dry ice mass variation is shown for two cases of foam thermal conductivity, namely, k f = 0.025 W m − 1 K − 1 and k f = 0.023 W m − 1 K − 1 . It can be seen that foam thermal conductivity affects the accuracy of predicting sublimation rate and end-of-sublimation time of dry ice. For cuboid box made of EPS material, the prediction of end-of sublimation time is better for k f = 0.023 W m − 1 K − 1 whereas the value, k f = 0.025 W m − 1 K − 1 , gives an accurate prediction for cylindrical box made of GEPS material. However, in both the cases, the numerical model has a significant improvement over the state-of-the art models in predicting the dry ice mass variation and end-of-sublimation time.

Thermal conductivity of foam
In practice, the thermal conductivity of foam is dependent on several parameters, primary among them are the temperature, density and moisture content. At the microscopic level, it is dependent on cell size, fibers/particle arrangement, gas type and its pressure, and the type of cell (closed or open). The apparent thermal conductivity of different insulation materials, in the temperature range from 100 K to 300 K, was determined in the work of Sparks et al. [17]. It was reported that thermal conductivity of polyurethane foam (ρ = 37.6 kg m − 3 ) varied in the range from 0.0131 W m − 1 K − 1 at 106 K to 0.0245 W m − 1 K − 1 at 302 K. It's value at 194 K, which is close to the dry ice temperature, is 0.0233 W m − 1 K − 1 . For the purpose of modeling dry ice sublimation inside foam package up to engineering accuracy, the thermal conductivity value of foam is assumed to be constant whose value is evaluated as follows.
Gnip et al. [18] proposed an empirical relationship to estimate thermal conductivity of EPS foam in the interval 0 • C to 50 • C from the measured thermal conductivity of EPS at the mean temperature of 10 • C. The empirical relation is shown in Eq. (19) where k f10 • C is the measured  value of foam thermal conductivity at mean temperature of 10 • C and F 10 • C is the conversion coefficient varying linearly with temperature as shown in Eq. (20). T is the temperature at which the thermal conductivity is to be determined and the constants b 0 , b 1 are calculated according to experimental data.
The variation of F 10 • C with temperature according to Eq. (20) is shown in Fig. 11. The dashed line shows the linear extrapolation of the data up to dry ice temperature. The extrapolated value of F 10 • C at the mean temperature of dry ice and ambient i.e. at − 30 • C is also highlighted in Fig. 11. Using Eq. (19) and the value of F 10 • C at − 30 • C, we calculated thermal conductivity of EPS foam at − 30 • C to be equal to 0.025 W m − 1 K − 1 . In Eq. (19), the experimental value of thermal conductivity at 10 • C (0.0305 W m − 1 K − 1 ) was obtained from the work of Koru [11] for EPS foam having density 32 kg m − 3 which is close to the densities of the rigid EPS and GEPS foam used in our experiments. It is pointed out in Section 2 that the thermal conductivity of the GEPS and EPS foam differ marginally. Therefore, the values of thermal conductivities are assumed to be same for both the material in the model. The analysis presented in the following section considers a constant value of foam thermal conductivity equal to 0.025 W m − 1 K − 1 Fig. 8. Effect of Biot number on wall temperature.

Temperature data
The temperature data measured and predicted at different locations inside the insulation box is shown in Fig. 12. Measured temperature at the bottom center (TC2, see Fig. 3) shows a constant behavior before it starts rising near the end of experiment. The constant behavior indicates the presence of dry ice pellets in the vicinity of TC2, and the rising up of temperature indicate the sensible heating up to the ambient temperature when all the the dry ice is sublimated. In the model, the dry ice domain is assumed to be solid and the temperature boundary condition of 194.65 K is specified on all the boundaries of dry ice domain. Therefore, the predicted value of TC2 shows a constant temperature of 194.65 K over the entire period of time in Fig. 12. Referring to the model methodology, temperature rise up as measured by TC2 is not included because the sensible heating of the insulation box, when the dry ice is sublimated, does not play a role in determining the end-of-sublimation time of dry ice.
Temperature predicted on the outer wall of the insulation box is in fair agreement with the measured temperature data by TC1. Whereas prediction of temperature data corresponding to TC3 and TC4 deviate significantly from that of measured values. The deviation is because these TCs are surrounded by the CO 2 gas for majority of the time and the simplified model does not capture dynamic process of the gas being accumulated and leaving the box. However, the heat transfer rate from the gas domain does not contribute significantly towards sublimation of dry ice as shown in Fig. 13. It shows predicted value of total heat transfer rate (Q) into the box over time. Q can be divided into three parts namely, conductive heat transfer rate through the surfaces in contact with dry ice (Q d ), conductive heat transfer rate through the surfaces in contact with the CO 2 gas (Q g ) and radiative heat transfer rate through the gas (Q r ).
Only the conductive heat transfer contribution of gas and dry ice obtained by the surface integration of conductive heat flux on their respective boundaries is shown in Fig. 13. It can be seen from Fig. 13 that Q d has major contribution (80%-85%) to the total heat transfer rate whereas the contribution of Q g is approximately 10%-18%. Hence, it will suffice to correctly evaluate the heat transfer rate through the surfaces in contact with dry ice to predict the end-of-sublimation time up to engineering accuracy.

Effect of reflective layer
A single layer of aluminaized mylar foil (reflective layer) is applied on the internal walls of the cuboid box to observe the effect of surface reflectivity on the sublimation rate of dry ice. Heat transfer due to radiation is approximately modeled to observe if experimental trends are captured in the simulation. Constant and approximate values of emissivity of radiating surfaces are considered. It is assumed that radiative surfaces behave as gray and diffused bodies which means that (1) emissivity is equal to absorptivity, and (2) emissivity is independent of radiation wavelength. The first implication is true for most opaque bodies and the second is valid when radiative power is emitted in narrow spectral and temperature range which, in this case is approximately between 194 K -293 K.
In the past scientists have evaluated emissivity of solid CO 2 present on Martial poles to study seasonal changes of Mars's polar caps. Ditteon and Kieffer [19] measured transmission spectra of solid CO 2 samples grown on a substrate placed inside a vacuum chamber. Emissivity of solid CO 2 was evaluated using the relation, ∊ = 1 − R, where R is the Fig. 11. Coefficient F 10 • C for temperature conversion of thermal conductivity of expanded polystyrene foam [18].  normal Fresnel reflectance depending on real and imaginary refractive indices of the sample. The 20 μm band emissivity is estimated to be in the range 0.62 to 0.89, and it was concluded that CO 2 deposits on Mars may be considerably lower than 1. Furthermore, Warren et al. [20] studied the effect of different parameters like grain size, dust and H 2 O contamination etc on the radiative properties of solid CO 2 . A model developed for terrestrial snow was extended to evaluate spectral albedo and emissivity of Maritian solid CO 2 . The albedo of pure CO 2 snow (grain radius of 100 μm) is calculated to be 75%-80% at visible wavelengths. It was reported that emissivity is sensitive to grain size and emission angle, and it varies significantly with wavelength. Since radiative properties for dry ice pellets are not available, we have considered an emissivity value of the order 0.1 for dry ice surface in the model.
The emissivity value of 0.6 is used for EPS as reported by Juanico [21], who designed and analyzed thermal insulation consisting of multiple air cavities inside insulation material for the roofs of buildings. In the present work, the presence of reflective layer is simulated by lowering foam emissivity value because a reflective layer has high reflectivity and low emissivity. Domen [22] measured emissivity of aluminized mylar sheet to accurately estimate the reduction in heat losses when the mylar sheet is epoxied on calorimeters that are made of graphite or A-150 plastic. The measurements were performed at 35.5 • C with an average room temperature of 23.5 • C. The emissivity value of 0.044 was determined for 6 μm thick commercially available aluminized mylar sheet. The Handbook of Optical materials by Weber [23] lists emissivity values for polished aluminum metal in the range of 0.04 to 0.06 for the temperatures varying from 50 • C to 500 • C. In our measurements, the reflective layer is exposed to the temperatures that are well below ambient temperature (20 • C). Therefore, emissivity value of 0.02 for reflective layer is assumed in the model which is slightly lower than the values reported in [22,23]. Fig. 14 shows the effect of reflective layer (RL) on the experimental and predicted values of sublimation rate and temperature inside the insulation box. Experimental data on Fig. 14 (a) shows that the presence of RL reduces the sublimation rate, and dry ice takes takes longer time to sublimate. After 50 h the amount of dry ice that sublimates in the cuboid box with RL is approximately 26% less as compared to the insulation box that does not contain RL. Although the model under predicts the sublimation time of dry ice, it shows similar trend as experiments when RL is present i.e. the sublimation rate decreases when RL is considered in the model. The temperature data is shown in Fig. 14 (b). The data of TC2 (see Fig. 3) is shown to confirm the reduced sublimation rate when the RL is present. It can be seen that the measured temperature data by TC2 remains constant for majority of the time irrespective of the RL layer because the dry ice is always in the vicinity of TC2. However, the measured rise up of temperature data by TC2, in the presence of RL, starts approximately 8 h later than the box that does not contain RL. The predicted value of TC2 does not show this behavior because the sensible heating of foam box up to the ambient temperature is not considered in the model.
The inner surface of the lid from the insulation box is facing the dry ice interface with a view factor close to 1, so a prominent effect of RL on the temperature of this surface is expected. The TC3 is located on the inner surface of the lid whose measured and predicted data are shown in Fig. 14 (b). TC3 measures significantly higher temperatures when RL is present because the RL reduces emission and enhances reflection of radiation causing a rise in its temperature. Similar trend is observed in case of results from the model. Since the heat transfer rate (Q) is proportional to the temperature difference (ΔT), a higher temperature of inner surface of lid in presence of RL reduces Q from the ambient to the box lid. The reduced heat transfer rate reaches the dry ice via the gas domain. This causes the dry ice to sublimate at a slower rate when the inner surface of insulation box is covered with a reflective layer.
To summarize, we present a comparison between existing models and the current work. The existing models are based on the heat balance shown in Eq. (5), where the left hand side i.e. the heat transfer rate into the insulation package is constant during sublimation. In practice, the heat transfer rate reduces over time becauses the total inner wall area in contact with dry ice reduces during sublimation. Consequently, the exisiting models predict higher sublimation rates and lower end-ofsublimation time as compared to the modelling approach discussed in this paper. Although the model presented in this paper is able to evaluate end-of-sublimation time accurately, it needs to be improved in estimating the temperature profile in the gas domain.

Conclusion
In this paper, the problem of dry ice sublimation in an insulation package, commonly encountered in shipping temperature sensitive products, was investigated experimentally and numerically. Mass and temperature measurements were carried out for insulation packages of two different geometries, namely cuboid and cylinder made of expanded polystyrene and graphite filled expanded polystyrene, respectively. Since the sublimation enthalpy of dry ice is large, the Stefan number is small compared to one. Also, based on the first estimate, it is shown that the expected sublimation rate of dry ice in the insulation package is low. Thus, the dry ice interface moves slowly and the temperature at each instant during sublimation corresponds to steady state. Based on this observation, a simplified three dimensional iterative model was developed in COMSOL to predict the dry ice mass variation and temperature inside the insulation package. The variation of dry ice mass inside the insulation package depends on the accurate prediction of heat transfer rate from the ambient to dry ice. From the model results it was shown that majority of heat transfer rate occur through the inner wall area of the insulation package in contact with dry ice. Further, it was shown experimentally and numerically that a reflective layer on the inner walls of the insulation package reduces the rate at which dry ice sublimates. Thus, creating a possibility to keep the perishables cold for a longer Fig. 14. Effect of reflective layer on the sublimation rate and temperature data inside the cuboid box containing dry ice pellets. period of time. Within the scope of this work, the model is able to predict the variation of dry ice mass and end-of-sublimation time accurately as compared to the existing models. However, the temperature prediction in the gas domain needs to be improved by determining the mechanisms of flow and heat transfer to better understand sublimation inside insulation packages. The proposed modeling approach is relatively easy to implement which can be used by aircraft carriers to predict dry ice mass variation over time and establish CO 2 limits systematically as opposed to the current limits based on rule-of-thumb.