Text-only Preview

"Food for Thought - Cool" AIRAH Conference, Darling Harbour, Sydney, Australia, 24 May, 2002
Q. Tuan Pham
School of Chemical Engineering and Industrial Chemistry
University of New South Wales, Sydney 2052, Australia
email: [email protected]
Refrigeration can be costly in terms of equipment and energy, and if not done correctly will fail to
achieve its objectives and lower the quality and safety of the product. To ensure that refrigeration is
effective, we need to be able to calculate processing times, product temperatures, heat loads and water
diffusion into and out of the product.
The calculation of food cooling and freezing may be complicated by phase change. The product
undergoing cooling often has very complex shape and composition. Heat transfer may be coupled with
moisture transfer and the equations governing the two processes should be solved simultaneously. The
heat transfer coefficient is often difficult to determine for the infinite variety of real-life situations such
as packaged products, cryogenic cooling, highly turbulent flow, swirling and non-parallel flow, etc. and
within the same chiller the heat transfer coefficient will vary greatly from place to place. The heat
transfer coefficient may also vary along the surface of a product due to the complex turbulent flow
pattern and the development of the flow boundary layer, a problem still largely unsolved even by the
most advanced computational fluid dynamics packages. Food products often have inconsistent
compositions, shapes and sizes, which lead to variable thermal behaviour and quality after processing,
so that it is difficult to quantify accurately the effects of different processing practices.
This paper will review methods for modelling foods with complex shapes, from highly simplified
models to detailed numerical models including CFD. Typical results on the accuracy of various
methods will be reported.
Cooling time, freezing time, product heat load, modelling
From the refrigeration engineer's point of view, the most important parameters determining the design
of cooling and freezing processes and equipment are the processing time and the heat load. The former
determines the size of the processing equipment (chiller, freezer) since holding capacity is the product
of throughput and residence time, while the latter determines the capacity of the refrigerating
equipment (compressors, condensers) and the energy consumption.
Twenty years ago, the refrigeration engineer would estimate processing time and product heat load
from experience or, if he was really "nerdy", by using some simple empirical equation or chart. There
is a lot of uncertainty in using such methods but as products and processing methods tended to change
only slowly, design engineers had time to experiment and learn from their experience.
In the last twenty years, things have changed greatly. New markets and products are being developed
more and more quickly, while computers have also progressed at an increasing rate, doubling in speed
and capacity every three years. Commercial calculation packages such as finite element methods and
computational fluid dynamics have appeared in the market, giving the designer new tools that are very

powerful but also at times dangerous to use. This paper will attempt to give a review of the range of
computational methods that are now available.
During the chilling and freezing of solids, heat must be conducted through layers of material before
reaching the surface, then cross over to the surrounding fluid either directly or through layers of
packaging. Processing time and heat load depend on the properties of the product (thermal
conductivity, density and specific heat), the surface heat transfer coefficient and the thermal resistance
of any packaging, and also on the size and shape of the product. The latter is often a major obstacle to
the accurate prediction of processing parameters, due to the complex shape of many products of
The processing time is usually defined as the time for the thermal centre of the product to reach a
certain temperature. In a complex geometry, it depends on what happens in the thickest part of the
object, for example the rump of a carcass. The heat load varies greatly during processing, being highest
at the start, and this peak heat load, which may determine refrigeration requirements, often depends on
what happens in the thinnest parts (biggest surface area per unit volume). Only in objects of simple,
regular shapes (spheres, long cylinders, large slabs) will both be determined by the same dimension.
When the internal resistance to heat transfer is negligible compared to the surface resistance, the whole
object will be at uniform temperature and the calculation of processing time and heat load becomes
trivial. In the more general case the Biot number, Bi, measures the ratio of internal to external
Bi = hR / k
where h is the surface heat transfer coefficient (including the effect of any packaging, radiation and
evaporative cooling), R the smallest half-dimension and k the thermal conductivity of the product. The
Biot number is of fundamental importance to any study of heat transfer to solid objects. Internal
resistance can be neglected if Bi << 1 (uniform internal temperature). Another limiting case is when Bi
tends to infinity, for which analytical solutions have been found for many cases. Many simplified
methods use interpolation to predict cooling times and heat loads for intermediate cases.
It is vital for practicing engineer to be aware of the Biot number of their thermal processes. If it is much
larger than 1 (i.e. if internal resistance is controlling heat transfer), then there is little use trying to
improve heat transfer say by reducing packaging or increasing air velocity, and more advantage will be
gained by trying to reduce the dimensions of the product (if possible); the inverse applies when Bi is
much smaller than 1.
Cooling time
When a solid object is cooled from a uniform initial temperature Ti in an medium at constant
temperature Ta , the temperature T at any location in the solid as well as the mean temperature obey the
following equation, which is a sum of exponential decay terms (Carslaw and Jaeger [1]):
(T - Ta )/(Ti - Ta ) = c1 exp(-b1t) + c2 exp(-b2t) + ...
where the coefficients bi and ci depend on the Biot number and the location in the product. The series
solution is difficult to compute but graphs of the solutions are provided in many heat transfer textbooks
for some simple shapes (e.g. Holman [2] p.148). However, the situation is greatly simplified if the
temperature fall is more than about 30% of the maximum that can be achieved (i.e. (T - Ta )/(Ti - Ta ) ?
0.7). In that case, all the terms in the series solution become negligible except for the first, which has

the smallest decay rate. Thus, after the initial period the residual temperature difference and heat load
decay exponentially. The centre temperature Tc and mean temperature Tm are described during the
exponential decay phase by
(Tc - Ta )/(Ti - Ta ) = jc exp(-2.303 t/f)
(Tm - Ta )/(Ti - Ta ) = jm exp(-2.303 t/f)
from which cooling time can be calculated. f is the time for the residual temperature difference or heat
load to decrease by a factor of 1/10 during the decay period, while j is called the "lag factor" (it
measures how much the actual decay curve lags behind a true exponential decay curve). While the j-
factor depends on the position in the product, f remains the same for all positions.
Equations are available for calculating j and f for basic shapes (slab, cylinder and sphere). Pflug [3]
plotted these solutions on a graph. Ramaswamy et al. [4] and Lacroix and Castaigne [5] gave
approximate solutions that can be computed with a calculator.
For complex or irregular shapes, there has been two main lines of approach to the prediction of cooling
time. The first is to approximate the object by an idealised (simple) shape with an equivalent thickness
or diameter. Thus, a lamb leg or beef rump can be approximated by an equivalent infinite cylinder [6,7]
or a beef carcass can be described as a collection of slabs and cylinders [8]. As long as the real shape does
not deviate too far from the ideal, the equivalent dimension can be calculated from, say, the volume/area
The second approach attempts to derive general expressions for the j and f factors of complex shapes
from their various dimensions. The most advanced of these methods is that of Lin et al. [9-11]. The
formulas are long and involved but simple to apply on a spreadsheet.
Heat load during cooling
For simple shapes, analytical and graphical solutions can be found in most heat transfer textbooks. The
f-j method should normally be avoided since it does not apply in the initial cooling phase, which is
when heat load is highest. However, if the average heat load Qav between times 0 and time t is needed,
it can be calculated from
Qav = mc (Ti - Tm) / t
where Tm can be found by the f-j method, as long as the exponential period has been reached. For
example, the average heat load Q0-0.7 between time 0 and time f. log10 (jm/0.7), at which the residual
mean temperature difference has fallen to 70% of its initial value, is
Q0-0.7 = 0.3 mc (Ti - Ta) / [f. log10 (jm/0.7)]
Freezing time
The freezing of a food can be divided into three distinct periods:
1. A precooling period, where cooling takes place without phase change.
2. A phase change period, which starts when the surface reaches the freezing point Tf and ice forms at
the surface. The freezing front gradually advances into the product, until it reaches the thermal
centre and the whole food can be considered frozen.
3. A postcooling period, during which the whole food continues to cool below Tf.
Plank [12,13] presented an analytical solution for the phase change period only. For many years,
Plank's equation was the only one available and it is still taught in many textbooks. However, it can
lead to large errors (typically 30%). There is no analytical solution for food freezing that includes all
three periods. Furthermore food freezing is more complicated than the freezing of pure water, because
the water in the food does not freeze instantaneously but does so over a range of temperature. To

predict freezing time for a simple shape, many empirical equations have been presented in the last
twenty five years. A fairly simple equation which is nevertheless quite accurate (±10%) is that by the
author (Pham [14]):
? ?H
? 2 ?
f =
1 +
? ?T
? 2 ?
2 ?
where Bi = hR/kf , EFreeze is the equivalent heat transfer dimensionality (EHTD) shape factor for
freezing (1 for slab, 2 for infinite cylinders, 3 for spheres), ?H1 and ?T1 are the volumetric enthalpy
change and temperature difference respectively for the precooling period, and ?H2 and ?T2 those for
the combined freezing-postcooling period.
As in cooling, the freezing of an irregular shape can be approximated by that of an equivalent regular
shape; for example, an ellipsoid can be approximated by an equivalent sphere. The technique has also
been frequently used on specific products such as lamb legs.
Another approach for complex shapes is to first calculate the freezing time for a slab then divide it by
the shape factor EFreeze, which is available for many shapes (Hossain et al. [15, 16], Cleland et al.
[17,18], McNabb et al. 19,20]).
Freezing heat load
A simple method for calculating freezing heat loads was proposed by Lovatt et al. [21,22] (Figure1). It
is based on how the freezing front moves towards the thermal centre and involves an empirical shape
factor N, which may be different from the shape factor E for freezing time (because heat load depends
what happens in the thinnest parts of the product while freezing time depends on what happens in the
thickest parts).
Figure 1. Lovatt et al's method for freezing
heat load calculation.
Calculating Product Heat Load in an Industrial Cooler or Freezer
The methods described in this chapter enable the calculation of heat load from a single item of product.
However, in industry normally many items will be processed at one time, and it is important to take
into account the operational characteristics of the process. For a strict batch operation, when the cooler
is filled with product before refrigeration is turned on, the refrigeration must be able to cope with the
peak heat load under the starting condition. For a continuous operation, the heat load is constant and
can be calculated from the mean temperatures of the entering and leaving products and the throughput.
Many operations in the food industry operate in an intermediate mode; for example, a chiller or freezer
may be loaded/unloaded continuously during daytime and then continue to cool batchwise during the
night. At any time t, the product in the cooler will have residence times ranging from 0 to tmax and
hence varying heat loads. In such situations, the total heat load at time t is made up of the heat load
from product loaded at all times between t - tmax and t. The time from t - tmax to t can be divided into a
number of intervals, the heat load from product loaded in each interval calculated and added together.

In numerical discretization methods (finite differences, finite elements, finite volume), the product is
divided into small control volumes or elements. The heat conduction equation is written for each of
these control volumes/elements, giving a set of hundreds, thousands or sometimes millions of
equations. These equations are then solved together to calculate the change with time of the
temperature in each location.
Finite differences (FD) was the earliest such method. The product is represented by a regular grid and
the equations for heat flow between the "nodes" (grid points) are written out and solved by computers.
Nowadays, the solution of a finite difference problem on a computer is very fast, of the order of a few
second or less. The finite difference equations can be set up and solved by any mathematically
competent engineering graduate (sad to say, not all engineering graduates are mathematically
competent). However, this method is practical only for objects that have a reasonably regular shape,
such as a rectangular box or a finite cylinder.
Finite Elements (FE) is preferred to finite differences for modelling objects with complex, irregular
shapes. A grid is still used but it can be irregular, consisting of triangles, distorded rectangles or
volumes of various shapes. The grid can be in two or three dimensions. Finite elements models take
more time to solve than finite difference models due to the larger bandwith of the matrix equations that
need to be solved. On a modern PC, a two-dimensional FE model may take several minutes to solve
while a three-dimensional FE model may take an hour or more. Also, it is more difficult to set up the
model, usually requiring special graphical software. However this is not due to the nature of the model
but to the complex shape that such models are applied to. FD models are easier to set up but that is only
because the real shape is usually approximated by a more regular shape.
The theory behind finite elements is mathematically abstract, and nowadays the method of finite
volumes (FV)
, which combines the flexibility of finite elements with the conceptual simplicity of finite
differences, has become widely popular. The object is divided into small control volumes of arbitrary
shape (as in finite elements), and the equations of conservation for heat, mass and momentum are
written for each control volumes. Computational effort is similar for FV and FE. FV is widely used in
computational fluid dynamics (CFD) models.
One of the big uncertainties in calculating process time and heat load is the heat transfer coefficient at
the surface of the product. Often the air is highly turbulent, the shape is complex giving rise to
boundary layer variation and eddies, which make the heat transfer coefficient vary from place to place
and become highly unpredictable. In other cases such as cartoned product, the presence of air gaps with
irregular shapes containing natural convection cells make the effective heat transfer coefficient difficult
to predict. In these cases the technologist must often carry out a large number of experiments, or make
gross approximations based on experience.
Computational Fluid Dynamics offer the promise of eliminating these uncertainties, by calculating the
flow pattern and heat transfer at the product surface from first principles, using the basic equations of
heat, mass and momentum transfer. This has become possible owing to the huge computing power of
modern computers.
For the case of laminar flow, CFD has indeed lived up to its promise. However, most industrial
situations involve turbulent flow, for which CFD is still unable to entirely deliver the goods. This is
because turbulent flow involves stochastic variations in the flow and temperature pattern, which must
be averaged out, and during this averaging process, empirical equations and coefficients must be
introduced and the equations lose their fundamentality. They must incorporate turbulence models
which are still conceptually suspect and involve a large number of empirical coefficients. The empirical

coefficients are specific to each flow situation and if they are used beyond these (i.e. in almost every
new situation practice), accuracy is no longer guaranteed.
To illustrate these points, Pham and Nguyen's [23] CFD simulation results for heat transfer coefficients
during beef chilling is presented in Figure 2, using the standard k-? model and the RNG
(ReNormalisation Group Theory) k-? model respectively. Both of these are empirical models with
several empirical coefficients.
Another obstacle against the widespread use of CFD models is the large amount of developmental time
and computation effort involved. It took about three months of an experienced user's time to build the
finite volume of a beef side model, and simulating a 20-hour run took a week on a 300-MHz personal
computer, with frequent intervention from the user to optimise relaxation factors etc. (and that's
ignoring mass transfer).
Nevertheless, computer power continues to increase quickly, and one can expect that more fundamental
turbulence models (such as the large eddy simulation model) will eventually become practicable for the
food industry and yield reliable results.
Pseudo-steady state CFD simulation
In view of the large computing time required to simulate a transient problem, a good compromise is to
use CFD to calculate the heat transfer coefficient for a steady state situation, which is much faster and
takes minutes instead of days. Since the htc varies only slowly with time, it can be assumed to remain
the same for a significant fraction of the process time. Thus, changes in temperature can be calculated
using a conduction model (FD, FE or FV) for the inside of the product only, which is again much faster
than a full transient CFD calculation. Figure 3 shows such a calculation for beef chilling (Pham and
Nguyen [23])
Temperature, deg.C
Heat transfer coefficienth (W/m2/K)
Time in chiller, h
Inlet velocity (m/s)
Figure 3. Leg centre and surface temperatures
Figure 2. Heat transfer coefficient predicted using
of beef side. ?? calculated by CFD using the
the standard k-? model (dotted curves) and the
RNG turbulence model, ?• measured.
RNG k-? model.
In view of the large range of calculation methods available it would be useful to compare them against
each other and against empirical data, where available. This will help the refrigeration engineer and
food engineer decide which to use should the need arise.

Errors in predicting cooling time
Prediction errors are due to the following sources:
• errors inherent the calculation method (for example, by approximating the cooling curve by the f-j
• errors in the product's thermal properties (thermal diffusivity)
• errors in the heat transfer coefficient – effect of turbulence, radiation, evaporation, wrapping etc.
• errors in predicting the effect of complex shape
The first source of errors can be eliminated by using a strict analytical or numerical method, however
the other sources of errors cannot be entirely eliminated. Surprisingly little data is available on the
extent of these errors in practice. However, we can expect that they will be of the same order as those
incurred during freezing calculations, for which data are available, as will be seen below.
Errors in predicting freezing time
Cleland [24] did some careful analysis of the most well known freezing time prediction methods (Table
1). It can be seen that the most accurate methods, such as that of Pham presented above, can predict
freezing time to about ±10%. There is also quite good agreement between simplified formulas and
numerical methods such as finite differences.
Table 1. Comparison of simplified
freezing time prediction methods
with experiments and finite
difference calculations.
It must be noted that these error values apply to experiments carried out in the laboratory under
carefully monitored conditions, with well known test materials, so that the errors due to heat transfer
coefficients, product shape and material properties are minimal. In industry, larger errors can be
expected, although with some expertise, errors of perhaps no more than ± 15% can be achieved.
Errors in predicting cooling and freezing heat load
The author and his associates at MIRINZ, Massey University (NZ) and the University of New South
Wales seem to have been the only group to have gathered reasonably accurate heat load data during
cooling and freezing. These data have been obtained using a flow calorimeter technique [25] (Figure
4). The hot product, for example a carton of meat or a beef side, is put in a wind tunnel. A thermopile
(a sensitive differential temperature sensor) is made by having an array of thermocouples distributed
over the inlet flow area and another over the outlet flow area. The thermocouple pairs are connected in
series, thus multiplying the thermoelectric signal. For example, with 20 thermocouple pairs, a voltage
signal of about 8 µV is obtained for every 0.010C temperature rise in the air, which can be easily
measured with a digital voltmeter accurate to 1µV. There is some errors introduced by heat loss
through the tunnel walls, but as long as conditions are stable these effects can be removed by a simple
baseline correction.
A typical heat load curve for beef chilling (Davey and Pham [8]) is shown in Figure 5. Davey and
Pham [26] did a comparison of predicted and measured heat load (averaged over the first two hours of

chilling) for 55 beef chilling tests (Table 2). It can be seen that the predictions for both a finite
differences model and a 2-D finite element models are quite acceptable, although there is a definite
improvement with the FE model (5.6% average error in heat load for FE vs 12.6% for FD).
FE Model
Davey and Pham (1997)
Beef Side Heat Load (W) 200
Time in Chiller (hr)
Figure 5. Comparison of heat load prediction by FE,
Figure 4. Wind tunnel cum flow calorimeter
FD and experimental results for a beef chilling
designed by Davey and Pham [25]
Table 1. Comparison of FD and FE model against experimental data for 55 tests.
FD Model
FE Model
Representation of beef side
7 simple-shaped regions
13 cross-sections
Average % error in heat removed during
-12.6 %
-5.6 %
first 2 hours
Average % error in weight loss after 20
2.32 %*
hours of chilling
Time to simulate 20 hour process on
< 1 min
4-5 hours
Pentium 166 Mhz Computer
A wide range of methods for predicting cooling/freezing times and heat loads, from simple
approximate equations to complex CFD models. Although sophisticated, theory-based methods are
becoming more practical with the fast progress in computer hardware and software, simplified methods
can still be useful in many cases. CFD is still basically a research tool in view of the major effort
involved in setting up the simulation, large computing time and continuing difficulty in dealing with
complex turbulent flow situations. For the majority of practical cases, a good practical compromise is
the use of finite differences or finite elements code for conduction inside the product, combined with
heat transfer coefficients obtained from empirical equations or CFD simulations.
[1] Carslaw, H.S., Jaeger, J.C. Conduction of Heat in Solids, 2nd edn; Oxford, Clarendon Press, 1959.
[2] Holman, J.P., Heat Transfer, 7th edn; New York, MacGraw-Hill, 1992.
[3] Pflug, I.J., Blaisdell, J.L., Kopelman, I.J. Developing temperature-time curves for objects that can
be approximated by a sphere, infinite plate or infinite cylinder. ASHRAE Transactions; 1965; 71; 238-
[4] Ramaswamy, H.S., Lo, K.V., Tung, M.A. Simplified equations for transient temperatures in
conductive foods with convective heat transfer at surface. Journal of Food Science; 1982; 47; 2042-

[5] Lacroix, C., Castaigne, F. Simple method for freezing time calculations for infinite flat slabs,
infinite cylinders and spheres. Canadian Institute of Food Science and Technology Journal; 1987;
20(4); 252-259.
[6] Earle R.L., Fleming, A.K. Cooling and freezing of lamb and mutton carcasses. Food Technology;
1967; 21; 79-84.
[7] Bailey C., Cox, R.P. The chilling of beef carcasses. Proceeding of the Institute of Refrigeration;
1976; 72; 76-90.
[8] Davey, L.M., Pham, Q.T. Predicting the dynamic product heat load and weight loss during beef
chilling using a multi-region finite difference approach. International Journal of Refrigeration; 1997;
20; 470-482.
[9] Lin, Z., Cleland, A.C., Sellarach, G.F., Cleland, D.J. A simple method for prediction of chilling
times for objects of two-dimensional irregular shape. International Journal of Refrigeration; 1996; 19;
[10] Lin, Z., Cleland, A.C., Sellarach, G.F., Cleland, D.J. A simple method for prediction of chilling
times: Extension to three-dimensional irregular shapes. International Journal of Refrigeration; 1996;
19; 107-114. Erratum, International Journal of Refrigeration; 2000; 23; 168.
[11] Lin, Z., Cleland, A.C., Sellarach, G.F., Cleland, D.J. Prediction of chilling times for objects of
regular multi-dimensional shapes using a general geometric factor. In Refrigeration Science and
Technology 1993-3. International Institute of Refrigeration, Paris; 1993; 259-267.
[12] Plank, R. Die Gefrierdauer von Eisblocken. Zeitschrift fur die gesamte Kalte Industrie; 1913;
20(6); 109-114.
[13] Plank, R. Beitrage zur Berechnung und Bewertung der Gefriergeschwindigkeit von Lebensmitteln.
Zeitschrift fur die gesamte Kalte Industrie; 1941; 3(10); 1-16.
[14] Pham, Q.T. Simplified equation for predicting the freezing time of foodstuffs. Journal of Food
Technology; 1986; 21; 209-219.
[15] Hossain, Md.M., Cleland, D.J., Cleland, A.C. Prediction of freezing and thawing times for foods
of regular multidimensional shape by using an analytically derived geometric factor. International
Journal of Refrigeration; 1992; 15(4); 227-234.
[16] Hossain, M.M., Cleland, D.J., Cleland, A.C. Prediction of freezing and thawing times for foods of
two-dimensional irregular shape by using a semi-analytical geometric factor. International Journal of
Refrigeration; 1992; 15(4); 235-240.
[17] Cleland D.J., Cleland A.C., Earle R.L. Prediction of freezing and thawing times for
multidimensional shapes by simple formulae. Part : Regular shapes. International Journal of
Refrigeration; 1987; 10; 156-164.
[18] Cleland D.J., Cleland A.C., Earle R.L. Prediction of freezing & thawing times for
multidimensional shapes by simple formulae. Part 2: Irregular shapes. International Journal of
Refrigeration; 1987; 10; 234-240.
[19] McNabb, A, Wake, G.C., Hossain, Md.M., Lambourne, R.D. Transition times
between steady states for heat conduction, Part I: General theory and some exact results. Occasional
Publications in Mathematics & Statistics No.20; Massey University; 1990.
[20] McNabb, A, Wake, G.C., Hossain, Md.M., Lambourne, R.D. Transition times
between steady states for heat conduction, Part II: Approximate solutions and examples. Massey
University, Occasional Publications in Mathematics & Statistics No.21. Massey University; 1990.
[21] Lovatt, S.J., Pham, Q.T., Loeffen, M.P.F., Cleland, A.C. A new method of predicting the time-
variability of product heat load during food cooling. Part 1: Theoretical considerations. Journal or Food
Engineering; 1993; 18; 13-36.
[22] Lovatt, S.J., Pham, Q.T., Loeffen, M.P.F., Cleland, A.C. A new method of predicting the time-
variability of product heat load during food cooling. Part 2: Experimental testing. Journal or Food
Engineering; 1993; 18; 37-62.
[23] Pham, Q.T. and Nguyen, A.V. Mean Convective Heat Transfer Coefficient of Beef Carcases
Computed by a Computational Fluid Dynamic. FOODSIM 2000: International Conference on
Simulation in Food and BioIndustries, Nantes, France, June 26-27, 2000.
[24] Cleland, A.C. Food Refrigeration Processes Analysis, Design and Simulation. London; Elsevier;

[25] Davey, L.M. and Pham, Q.T., A technique for measuring product heat load during beef chilling.
Published in Heat and Mass Transfer Australasia 1996: proc. 6th Australasian Heat and Mass Transfer
Conference, New York, Begell House, pp. 453-460, 1998.
[26] Davey, L.M. and Pham, Q.T., A comparison of three models for predicting heat and weight loss
from beef sides during chilling. 20th Int. Congress of Refrigeration, Sydney, 19-24 Sep 1999.