Blast-cooling of beef-in-sauce catering meals: numerical results based on a dynamic zero-order model

Beef-in-sauce catering meals under blast-cooling have been investigated in a research project which aims at quantitative HACCP (hazard analysis critical control point). In view of its prospective coupling to a predictive microbiology model proposed in the project, zero-order spatial dependence has proved to suitably predict meal temperatures in response to temperature variations in the cooling air. This approach has modelled heat transfer rates via the a priori unknown convective coefficient h c which is allowed to vary due to uncertainty and variability in the actual modus operandi of the chosen case study hospital kitchen. Implemented in MS Excel®, the numerical procedure has successfully combined the 4 th order Runge-Kutta method, to solve the governing equation, with non-linear optimization, via the built-in Solver, to determine the coefficient h c . In this work, the coefficient h c was assessed for 119 distinct recently-cooked meal samples whose temperature-time profiles were recorded in situ after 17 technical visits to the hospital kitchen over a year. The average value and standard deviation results were h c = 12.0 ± 4.1 W m - 2 K - 1 , whilst the lowest values (associated with the worst cooling scenarios) were about h c » 6.0 W m - 2 K - 1 .


Abbreviations and nomenclature
HACCP hazard analysis critical control point TM-0 thermal model of zero-order dependence on spatial coordinates TM-1 thermal model of first-order dependence on spatial coordinates Latin symbols A area for heat transfer from meal to air (m 2 ) Bi Biot number (dimensionless) c specific heat capacity of meals (J kg −1 K −1 ) D h hydraulic diameter (m) h c convective heat transfer coefficient (W m −2 K −1 ) L shelf level in the supporting trolley (dimensionless) I 1 /2 half-length separating meal core from horizontal (i.e., top or bottom) edges (m) m mass of meals in food containers (kg) N laid number of standardized containers laid in the supporting trolley (dimensionless) Nu D Nusselt number based on hydraulic diameter (dimensionless) Pr Prandtl number (dimensionless) Q heat transfer rate (W) Re D Reynolds number based on hydraulic diameter (dimensionless) T temperature (K or • C) t time (s) v f luid reference velocity of cooling air flow (m s −1 ) Greek symbols β lumped parameter (s −1 ) λ thermal conductivity (W m −1 K −1 ) υ kinematic viscosity (m 2 s −1 ) ρ density (kg m −3 ) Subscripts and superscripts air cooling air flowing around a given food container final shut-down of blast-cooling operation (final instant) food beef-in-sauce catering meal in a food container at a given position in the trolley inf below the food container (bottom) init start-up of blast-cooling operation (initial condition) meat meat portion of the catering meal num numerical meal temperature rec recorded meal temperature sauce sauce portion of the catering meal sup above the food container (top)

Introduction
In food service operations, cooling frequently raises concerns about Clostridium perfringens outbreaks (Steele & Wright, 2001;Juneja & Marks, 2002;Kalinowski, Tompkin, Bodnaruk, & Pruett, 2003;de Jong, Beumer, & Zwi-etering, 2005). In order to prevent such risks, public health authorities have set specific time-temperature combinations which are achievable via blast-cooling equipment (Doyle, 2002;Crouch & Golden, 2005). Yet, other combinations can be approved if operators are able to demonstrate that appropriate levels of IJFS October 2014 Volume 3 pages 213-227 Blast-cooling of catering meals: results from zero-order model 215 food safety are ensured. Accordingly, modelling and simulation are helpful for risk analysis as well as useful for process improvement and re-engineering (Bellara, McFarlane, Thomas, & Fryer, 2000;Huang, 2003;Amezquita, Weller, Wang, Thippareddi, & Burson, 2005;Almonacid, Simpson, & Teixeira, 2007). Following such a rationale, a research project on quantitative risk-based HACCP (hazard analysis critical control point), named Quant'HACCP Project, has proposed to model food chain processes on a modular basis, combining predictive microbiology and physical phenomena. Coldchain catering performed at a hospital kitchen was chosen as the case study and the present work focused on beef-in-sauce meals undergoing blast cooling. As part of the aforesaid project, two dynamic thermal models were tested to predict meal temperatures in response to temperature variations in the cooling air (Rabi, Trezzani-Harbelot, Morelli, & Guilpart, 2012). While both models relied on a convective coefficient h c to assess heat transfer between hot meals in containers and the cooling air flowing around them, they differed from each other with respect to the spatial dependence of temperatures. In the zero-order model, referred to as TM-0, temperatures could only depend on time t whereas one-dimensional variation was additionally allowed in the firstorder model (TM-1). In view of its prospective coupling to a predictive microbiology model proposed in another part of the Quant'HACCP project (Jaloustre, Cornu, Morelli, Noel, & Delignette-Muller, 2011), TM-0 was a suitable choice that could be implemented in MS Excel ® (Rabi et al., 2012). As claimed in Corradini, Amezquita, Normand, and Peleg (2006), MS Excel ® is a general-purpose software, finding widespread use among industrial food microbiologists, when compared to relatively advanced mathematical software like Mathematica ® , MathCAD ® or Maple ® . In TM-0 (and also in TM-1), coefficient h c is a priori unknown and allowed to vary due to uncertainty and variability in the usual (nonacademic) modus operandi of the kitchen. Coefficient h c is obtained by best-fitting the numerical meal temperature against the profile recorded in loco along with the temperature-time profile of the cooling air flowing nearby (Rabi et al., 2012). After 17 technical visits to the kitchen over a year, it was possible to apply the numerical procedure to 119 meal samples and this work presents and discusses the average and lowest values of those best-fitting coefficients h c .

Blast-Cooling Operation as Case Study
This section briefly describes the blast-cooling operation habitually carried out at the hospital kitchen taken as the case study. It was arbitrarily chosen from a list of Parisian hospitals performing cold-chain catering at their central kitchens. Filled up with recently-cooked beef-in-sauce meals, food containers are top-sealed with plastic film and laid in a trolley with 20 shelves, as shown in Fig. 1(a). For identification purposes, shelves were numbered upwards from L = 1 (bottom) to L = 20 (top). As Figure 1(b) shows, the trolley is placed in a blast-cooling cell, so that a given container is further identified as "front" or "back" whether it is near or opposite the cell door, respectively. Air currents are induced by four fans vertically lined up in one side of the cell, as shown in Fig.  1(c), so that trolley sides are referred to as "fan" or "wall". As described in (Rabi et al., 2012), warm air is aspired to an evaporator behind the fans while chilled air is blown back to the cell interior through two narrow vertical openings, one at each side of the fans from floor to ceiling. With standard dimensions (CEN, 1993), two sizes of equal-height rectangular containers were used, namely "1/1" (dimensions in mm: 530 × 325 × 65) and "1/2" (dimensions in mm: 265 × 325 × 65). On each trolley shelf, one may place up to 2 containers "1/1" or 4 containers "1/2". Each container "1/2" counts as 0.5 in the number N laid of containers laid in the trolley and N laid was noted to vary from one blast-cooling operation to another. As the idea was to let food operators pursue their customary (and, by the way, speedy) procedure, no suggestions were given about containers' placement. In other words, no pre-defined pattern was followed so that containers were laid at random in the trolley. A pictorial sketch was pro-  Fig. 2 for an invented load of N laid = 3×1 + 3×0.5 = 4.5 containers. In this illustrative example, containers placed at 5 th and 18 th shelf levels are identified as L-05 front and L-18 back, respectively. Over the investigation period, three types of beef-in-sauce meals were prepared, namely: "boeuf bourguignon", "boeuf mironton" and "goulash". Temperatures were measured through Proges-Plus small disc-shape recorders (diameter = 20 mm, height = 6 mm) and temperature-time profiles were recovered from disc memory via Thermotrack hardware and software. Recorders were placed in meat pieces at simmer start, Fig. 3(a), and those pieces were wrapped by a soft net for recognition, Fig. 3(b). For sauce temperature measurements, recorders were supported by a small plaque and wrapped by a soft net as well, Fig. 3(c). After containers were filled up with recently-cooked meals, those identifiable meat and sauce temperature recorders were placed at the container geometric centre, which tends to be the hottest point throughout the process. It is worth noting that not all containers had recorders therein. For air temperature measurements, recorders were placed during the very short time gap between containers' top-sealing with plastic film and their placement in the trolley.

Dynamic Thermal Model with Zero-Order Spatial Dependence (TM-0)
Seeking a compromise between application versatility and realistic description, a dynamic zeroorder model (TM-0) was proposed in Rabi et al. (2012) for the blast-cooling operation described in the previous section. It relies on lumpedparameter analysis (Özışık, 1985) so that temperatures may vary with time t only, i.e., no spatial dependence is allowed. This thermal model is briefly presented here for quick reference. Let m f ood and c f ood be the mass and specific heat capacity of a catering meal within a container at a given position in the trolley. By supposing that m f ood and c f ood remain constant, an energy balance yields: where T num = T num (t) is the meal temperaturetime profile to be numerically assessed where aṡ Q f ood→air,sup andQ f ood→air,inf are heat transfer rates to cooling air flowing above and below the container. Equation (1) assumes that (i) heat is predominantly transferred across the bottom and top horizontal surfaces of the container (of equal areas A inf = A sup = A) and (ii) temperature T num refers to the meal as a whole (i.e., zero-order model with respect to spatial dependence). In view of air flows induced in the blast-cooling cell, meals were assumed to ultimately lose heat by forced convection, whose comprehensive modelling is complex due to mutual interference between velocity and temperature in the fluidfluid (Kays & Crawford, 1993). Yet, computational fluid dynamics (Norton & Da-Wen, 2007) was considered beyond the modelling compromise and development strategy of Quant'HACCP Project so that heat transfer rates were linearized via a convective coefficient h c . No distinction was made between convective coefficients at top and bottom surfaces and, thus, the same value h c,inf = h c,sup = h c was assumed for a given container (Rabi et al., 2012). Still, coefficient h c was allowed to vary not only from one container to another in the trolley but also from one blast-cooling operation to another. It is worth recalling that lumped-parameter analysis is suitable only if the resulting Biot number is sufficiently low, namely Bi < 0.1 (Özışık, 1985). In line with an electrical network analogy (van der Sman, 2003), such a constraint is relaxed in TM-0 by accounting for thermal resistances between the meal core and each horizontal surface exposed to cooling air. By reasoning on the relative magnitude of such resistances (Rabi et al., 2012), one may sufficiently consider a conductive resistance for each half-length I 1 /2 between the meal core and each horizontal edge, being λ f ood meal thermal conductivity. If T air,sup (t) and T air,inf (t) are the temperatures of the cooling air flowing in the vicinity of the container, Eq. 1 then becomes: In the previous equation, lumped parameter β is defined as: while T air (t) is the temperature-time profile averaged from those recorded for the cooling air flowing below and above the container, as indicated in Eq. 3. One may interpret T air (t) as a representative temperature-time profile of the cooling air around the container. Solution of Eq. 3 requires an initial condition, namely: where T init is obtained from the temperatures recorded at operation's start-up.

Thermophysical Parameters and Numerical Solution Method
In TM-0, meals are assumed to have uniform specific heat capacity c f ood , thermal conductivity λ f ood and density ρ f ood , which are assessed from meat and sauce values via a weighted-average procedure (Sepúlveda & Barbosa-Cánovas, 2003 They were assumed to represent the catering meals studied and to remain constant during the blast-cooling operation. For each container with meat, sauce and air temperature recordings, meal profile T rec (t) was assessed in line with Rabi et al. (2012). The numerical profile T num (t) was evaluated in response to profile T air (t) by means of Eq. 3 together with a best-fitting coefficient h c . By recalling that profiles are discrete, such an h c value is the one that minimises squared differences between T num (t) and T rec (t) according to: where t f inal is blast-cooling final instant (shutdown), which varied from one operation to another. As temperatures were recorded for every minute, time t behaves like a dimensionless counter in Eq. 6, being t f inal the final index. Consequently, ∆T 2 av has the same units as [T rec (t) − T num (t)] 2 . As stated, temperature-time profiles were IJFS October 2014 Volume 3 pages 213-227 Blast-cooling of catering meals: results from zero-order model 219 Only for the particular case where air temperature T air is constant, one is able to deduce an analytical solution of Eq. 3 subjected to Eq. 5. This is not the case of the blast-cooling operation studied so that Eq. 3 was numerically solved. Specifically, the 4 th order Runge-Kutta method (Kreyszig & Norminton, 1993) was implemented using MS Excel ® . In the same spreadsheet, the "Solver" optimization tool was used so as to determine the best-fitting coefficient h c minimising ∆T 2 av in Eq. 6. Cell assessing ∆T 2 av was defined as the target to be minimised while the cell containing h c value was defined as the sole adjustable one. The non-negative values option was checked while the linear model was left unchecked.

Results and Discussion
Temperatures were recorded for containers that food operators randomly placed at various positions in the trolley. From the whole set of measurements, 119 containers had meal temperatures recorded along with the temperature of the cooling air flowing nearby, thus enabling the determination of 119 best-fitting coefficients h c . Using the pictorial sketch of Fig. 2 along with the number N laid of containers in the trolley, Fig. 4 shows the arrangements noted in technical visits over one year. Containers with numbers are those with temperature recorders so that the best-fitting h c could be determined. Those numbers are those h c values in SI units. Although the arrangement could not be fully registered in the first 3 technical visits, it was possible to identify the position of containers whose best-fitting h c could be determined. The number N laid of containers was not the same from one operation to another because it varied in proportion to the amount of meals cooked but N laid remained below 50% of the trolley's total capacity in all visits. Apart from the unique arrangement of a full load (unlikely attainable in view of the customary quantity of meals cooked), there are distinct ways to lay containers in the trolley, which was indeed the case as shown in Fig. 4. Besides time dependence already accounted in TM-0, the number N laid of containers as well as their arrangements in the trolley are expected to affect the cooling performance (Barbin & Junior, 2011). With respect to air temperatures, Fig. 5 shows profiles T air (t) concerning the lowest and highest loads, namely (a) N laid = 12.0 containers (4 th technical visit) and (b) N laid = 18.5 containers (6 th technical visit). Profiles follow a similar trend with cooling air temperatures reaching values as low as -25 • C. Local peaks in profiles correspond to eventual openings of the cell door by food operators for inspection purposes. As each best-fitting coefficient h c can only ensure that differences between recorded and numerical temperature-time profiles are minimal, it is important to check how close T num (t) is to T rec (t). Figure 6 compares those profiles for the worst and best objective function values, namely (a) ∆T 2 av = 5.85 • C 2 (container L-16 front, 3 rd technical visit) and (b) ∆T 2 av = 0.03 • C 2 (container L-16 back, 16 th technical visit). For completeness, corresponding air temperature-time profiles T air (t) are also shown in Fig. 6. As commented, coefficient h c was expected to differ from one meal to another as inherent uncertainty and/or variability was transferred to it via model equations. Expressed as the number of counts and obtained via the MS Excel ® built-in function "Frequency", Fig. 7 presents the distribution of all 119 best-fitting coefficients h c . a) b)

Figure 6:
Recorded and numerical meal temperature-time profiles, T rec(t) and T num(t), related to (a) the worst (∆T 2 av 5.85 • C 2 ) and (b) the best (∆T 2 av 0.03 • C 2 ) objective function values, together with the corresponding cooling air profile T air (t) The average value and standard deviation results were h c = 12.0 ± 4.1 W m −2 K −1 . Table 2 shows minimum, maximum and average best-fitting h c (together with standard deviations) as grouped by meal type or by the number N laid of containers in the trolley. With respect to meal type, minimum values (linked to worst cooling scenarios) are similar, h c ≈ 6.0 W m −2 K −1 , while the standard deviation is greater for "boeuf bourguignon". With regard to the number N laid of containers, minimum h c values seem unrelated while average h c values do not seem to follow a specific trend, as shown in Fig. 8. Best-fitting coefficients h c can indeed be dissimilar for the same N laid thus ren- As far as comparison with h c values from empirical correlations is concerned, a first issue refers to whether or not a classical geometry prevails, e.g. internal flow through a rectangular duct or external flow over a flat plate. It may even resemble a combination of these two and, as suggested in Fig. 10(a), air currents are likely to follow curved split pathways. In addition, the variability in container arrangements (Fig. 4) may complicate the evaluation of necessary parameters such as cross-sectional area, wetted perimeter and reference length. For trial calculations, let the Dittus-Boelter correlation be invoked, while recalling that air is heated by containers (Özışık, 1985): where the later being comparable to the average h c value already reported. For the blast-cooling operation studied in this work, applicability of lumped-parameter analysis (i.e. zero-order modelling with respect to spatial coordinates) depends on whether the Biot number Bi results are sufficiently small (Özışık, 1985), as stated in section 3. By relying on values in Table 1 which is slightly lower than the smallest bestfitting h c (worst cooling scenario). Before jumping to the conclusion that high-order modelling with respect to spatial coordinates is strictly necessary, it is worth highlighting that the regular (and non-academic) blast-cooling operation was disturbed as little as possible. Accordingly, a balance between accuracy and process variability was attempted with reference to representative thermophysical properties and recorded profiles. Still, as TM-0 accounts for meal thermal resistance, numerical meal temperature profiles fit well to experimental counterparts, as previously suggested in Fig. 6. In order to verify how T num (t) closely follows T rec (t) for all 119 meal samples studied, R 2 values were assessed through the MS Excel ® built-in function "RSQ" by assuming a linear regression based on those two profiles. Encouraging results are shown in Fig.  11 for all 119 best-fitting coefficients h c . Though distinct convective coefficients for bottom and top surfaces (hc,inf = hc,sup) could bring an extra degree of freedom, the best-fitting procedure based on Eq. (6) would become more complex and beyond the scope of this work. Process modelling should rely on governing equations, correlations and parameters with complexity balancing between the following issues: • Variability of parameters that characterize either food products (e.g. composition, homogeneity and amount in containers) or process (e.g. equipment operation, prevailing ambient conditions and arrangement inside the equipment); • Accuracy of input data (e.g. process parameters) and desired outputs (e.g. food safety and/or performance objectives, process optimization and/or re-engineering).
As claimed in de Souza-Santos (2004), sophistication is not a guarantee for model quality and it should go only to the point where key variables can be compared with available experimental data. Moreover, complex models are prone to require information that cannot be straightforwardly measured in a non-academic environment like the hospital kitchen investigated. It is worth recalling that the Quant'HACCP Project seeks unsophisticated approaches that may render general recommendations to operators, whilst being easily adaptable to similar case studies. Bearing in mind that catering chains may follow distinct processes, care should be exercised so as to prevent over-complicating or over-simplifying a model for use. Such is a critical task as relatively simple models may satisfactorily work for predictive microbiology purposes (Buchanan, Whiting, & Damert, 1997).

IJFS October 2014 Volume 3 pages 213-227
Blast-cooling of catering meals: results from zero-order model 225 All beef-in-sauce meals 119 5.6 32.7 12.0 4.1 a The number N laid of food containers in the trolley could not be observed in one blast-cooling operation for "boeuf bourguignon" and in one blast-cooling operation for "boeuf mironton", as indicated in Fig. 4 Figure 11: R 2 values from the linear regression of numerical profiles T num (t) against the counterpart recorded T rec (t) for all best-fitting coefficients h c as calculated via MS Excel ® built-in function "RSQ".

Conclusions
It is believed that inherent variability and uncertainties were transferred to the coefficient h c so that the distribution of related best-fitting values may have been affected by non-uniformities attributable (but not restricted) to meal composition, meal quantity inside containers, position of temperature recorders with respect to meal core, number of containers placed in the supporting trolley, and their layout therein. In view of non-academic environments (as in the hospital kitchen studied), zero-order thermal model with respect to spatial dependence (TM-0) proved to be a suitable approach. Despite the corresponding Biot number being slightly above 0.1, recorded and numerical meal temperature-time profiles correlated very well for several bestfitting coefficients h c . TM-0 was implemented in MS Excel ® , enabling its built-in Solver tool to be successfully combined with the 4 th order Runge-Kutta method. As this work belongs to a research project which aims at general recommendations for operators during food thermal processing, use of widespread software like MS Excel ® is definitely of interest.