897
Biotechnol. Rog. 1991, 7, 397-402
Mathematical Modeling of Biochemical Changes during Processing of Liquid Foods and Solutions Ashim K. Datta Department of Agricultural and Biological Engineering, Cornel1 University, Riley-Robb Hall, Ithaca, New York 14853
Thermal processing of liquid food materials and biochemical solutions always results in a distribution of biochemical changes such as bacterial and nutrient concentrations due to variation in temperature history from one liquid element to another. T o estimate such distributions, the liquid elements need to be "tagged" and followed, which is computationally difficult for most flow situations of interest. This paper describes a computational procedure for obtaining the full range of biochemical changes during processing without explicitly following liquid elements. By use of the method developed, transient distributions of bacterial death have been obtained for sterilization of liquid in a nonagitated cylindrical enclosure heated from all sides. The conventional procedure of obtaining sterilization based on temperature at the slowest heating zone is shown to be lower than the actual least sterilization of a real fluid element in the system. Nutrient distributions can be similarly calculated by following this method.
Introduction Any finite-size material would experience a spatially varying transient temperature distribution at its interior when heated or cooled at the boundary. Such unavoidable time-temperature distributions during a heat transfer process have special implications when applied to thermal processing (e.g., pasteurization, sterilization, evaporation, etc.) of food materials. The extent of chemical and biological reactions that go on in a food material during heating, e.g., destruction of microorganisms and nutrients (such as vitamin C), are temperature dependent, and therefore distributions of temperature would result in unavoidable distributions of concentrations of these reaction products. Such distributions of bacterial or nutrient concentration are relatively easier to obtain for solid food materials under conduction heating, where the food particles stay fixed in one location. To obtain a distribution of concentrations of the food particles, it is sufficient to know the transient temperature at various fixed locations within the food material. In a flowing fluid, however, the situation is quite different. Due to the movement of a fluid element, to know its temperature it needs to be tracked with time. For simple convective flows such as laminar flow in a tube where particles move in well-defined paths, temperature along the element path may not be too difficult to obtain. However, for complex heating situations as turbulent flow in a tube or recirculation flow in a container, the fluid element needs to be followed along its generally complicated path. Knowledge of a range of these biochemical changes with their associated volume fractions (a distribution) is the only way to provide complete information about the process. Single point value currently used in the food industry such as sterilization at the slowest heating zone does not provide a complete picture of the process. Distributions can more effectively compare two processes and therefore be useful for process feasibility. Such distributions also provide a lower bound on biochemical changes associated with real particles of the food system, in contrast with the current procedure of estimating the lower bound at the slowest heating zone that in reality 8756-7938/91/3007-0397802.50/0
may not correspond to any food particle. For most flow situations, however, computation of such distribution becomes quite complex. Except for the work of Paulsson and Tragardh (1984)and Ellborg (1989)for turbulent flow in a tube heat exchanger, distributions of sterilization or other quality parameters in foods have not been obtained. The objectives of this paper are to (1)develop a computationally efficient procedure for obtaining the range of biochemical changes in a nonisothermal flowing fluid and (2) apply the procedure to estimate the distribution of sterilization in unagitated heating of liquid in a container.
Problem Formulation Thermokinetics in Food Processing. Thermally induced changes in (biochemical)quality parameters such as nutrient or bacterial concentration in food material during processing are functions of the time-temperature history of the food particle. The rates of reactions in foods such as destruction of microorganisms or nutrients are generally considered first-order reactions (Karel et al., 1975) described by
To obtain the final concentration in changing temperature, eq 1 is integrated to obtain
In food literature, eq 2 is used in a different but equivalent form. Instead of referring to a final concentration, an equivalent heating time of Fo is used that gives the same final concentration when temperature T is constant at a reference value. Thus FO is a measure of extent of thermal treatment. Traditionally POis defined for concentration of bacteria (sterilization), but FOvalues can also be calculated for nutrients or other thermally vulnerable components. Using this definition of Fo at a reference temperature of TRfor which the rate constant
0 1991 American Chemical Society and American Institute of Chemical Englneers
Blotechnol. Rag., 1991, Vol. 7, No. 5
898
is k ~ yields ,
time. The most appropriate formulation for this situation is the Lagrangian formulation where flow quantities such as temperature and velocities of chosen material elements of the fluid are defined as functions of time (Batchelor, 1967). Thus the energy equation can be written as
which can be simplified for process temperature T staying close to TRas
(6) Dt = aV2T The notation D stands for the substantial derivative to remind us that it is computed for a substance. The total change in temperature D T in eq 6 can be interpreted as the change due to time variation as well as the convection of the fluid element into a zone of different temperature. The solution to eq 6 in this formulation would directly provide the transient temperature history of a fluid element, which can be combined with eq 3 to obtain ita FOvalue. This concept of following a fluid element with changing temperatures have not been discussed in great detail in the food literature, except for the work of Tragardh and Paulsson (1985) and Datta and Liu (1990). Difficulties in Using Lagrangian Formulation. Although the Lagrangian form (eq 6) is conceptually the most appropriate to use, due to the computational difficulties the flow equations are often solved in the more familiar Eulerian form, which is written as
(4)
Equation 4 is the familiar equation in food sterilization. The parameter 2 replaces the activation energy E, to describe the temperature dependence of reaction. Error involved due to linearization in eq 4 has been discussed in Datta (1991). Thus, instead of discussing the final concentration of components, equivalent heating time FO at a reference temperature TRis used for a nonisothermal process. The extent of sterilization in a process would be equivalent to its FOvalue. A reference temperature of TR equal to 121 "C is often used in the food industry and is used for this study. Distribution of Thermokinetic Parameters. Since the time-temperature history of any finite mass would vary spatially, the FOvalue as given by eq 3 would vary from particle to particle at any given time. Thus, in reality there is always a distribution of quality parameters such as Fo value due to variation in time-temperature history. This distribution is, however, different from residence time distribution (RTD). RTDs only have aduration of heating and not the temperature history and therefore cannot provide information on concentration changes in general. It is more useful for isothermal applications. For nonisothermal reactors, a quantity called thermal time has been defined by Nauman (1977) as
which includes the temperature history T(t) of a particle in the reactor. The time t , can be interpreted as the time at a constant temperature that leads to the same final concentration of reactants as the time at varying temperature T(t)for a duration oft. This definition of thermal time is quite similar to the Fo(t) value defined earlier. Both t , and F&) represent equivalent heating (or cooling) time at a reference temperature that leads to the same final concentration. Thus Fo is a measure of thermal time. Generally speaking, every particle in a nonisothermal reactor will have a different time-temperature history, producing a distribution of FO values that is a thermal time distribution. This thermal time distribution (TTD) can be studied in a manner analogous to RTDs as explained later. The area under the TTD curve between two values of thermal time FOis the fraction of liquid volume between those two FOvalues. In an enclosure, the total processing time t in the upper limit of the integral in eq 5 would be the same for all fluid elements but the value of the integral will differ depending on temperature history of the elements. Need for Lagrangian Formulation for Fluid Food Processing. Although the temperature involved in eq 1 is for a material particle, food literature almost exclusively refers to temperatures at fixed locations even in convective heating situations. For conduction heating, temperature T of a particle is identical with the temperature at its location at all times, since the food particle does not move with time. However, temperature of a fluid element during convective heating can only be known by following it with
at
+ u-dT + v-dT = aV2T ax ar
(7)
where the solution would specify a picture of spatial distribution of temperature at each instant during the motion. Information about bacterial or nutrient concentration of a particle along its flow path is not directly obtained from this most common type of flow specification. To obtain various flow quantities along a particle path, one possibility is to track locations of individual fluid elements by using calculated velocity fields as
where (ro,xo)is the initial location of the particle at t = 0 and (r,x) is the location at time t. Since the velocity field ( u p )is known (analytically or numerically), it can be substituted to find the transient locations of the particle. At these transient locations of the particle, its temperature would be found from the known temperature field a t the desired instant of time. Knowing its temperature, transient sterilization can be calculated as before by using eq 3. By repeating this procedure for a number of particles with different starting locations, a distribution of sterilization can be achieved. This is quite cumbersome and larger time steps create problems in successful integration of eq 8. Proposed Eulerian Formulation. The field of temperature obtained by using Eulerian formulation represents temperatures at the spatial locations. It is important to observe that the distribution of field values of temperature a t any given time is identical with the distribution of temperature values of the fluid elements in the system obtained by solving the energy equations in their Lagrangian formulation. Thus, distribution of a variable such as temperature can be obtained from the standard Eulerian formulation that includes information on all fluid elements in the system, but no individual fluid element (with an individual starting location) can be identified. Since the intent of this study was to obtain a distribution of sterilization, it should be possible to obtain one if the
399
Bbtechnol. Rag., 1991, Vol. 7, No. 5
sterilization variable FO can be treated the same way as temperature. One possibility is to consider the bacteria, nutrients, or any other biochemical quality variable as a different species and consider a mass balance equation for the species. For a flowing fluid, eq 1can be rewritten for concentration of bacteria as - -Dc =
Dt k f l where the contribution due to diffusion was considered very small (Tragardh and Paulsson, 1985)since many food constituents are of large molecular size or bound to structures (particles). Equation 9 can be expanded in the Eulerian form as
By the previous arguments on the interpretation of the Eulerian solution, the solution to eq 10 will provide transient concentrations at spatial locations that are identical with the distribution of concentrations obtained by explicitly following the fluid elements as in a Lagrangian formulation. Thus, concentration distributions can be obtained for bacteria or nutrients by using standard Eulerian formulation. To obtain a distribution of FOinstead of concentration c that is more common in food processing, a different procedure is now discussed. Using sterilization equation 4 as a basis for Fo,it can be rewritten in the following manner:
since it refers to a fluid element as opposed to a location. Instead of using eq 4, eq 3 may also be used to derive an expression similar to eq 11. This would not have the approximations in eq 4 and would be the direct equivalent of concentrations calculated from eq 9. In this study, eq 4 was used to be consistent with the food literature. Equation 11 can be expanded as
Although this equation was not derived from the transport equation, it is similar. In comparing with the energy equation (eq 71, the left sides of eqs 7 and 12 are identical. On the right side of eq 12, the thermal diffusion term has been replaced. For known velocities u and u and for known temperature T, eq 12 can be solved in a similar manner as the energy equation to obtain a field for Fo. By the previous arguments, this field for FO or sterilization is identical with the distribution of sterilization at that time instant obtained by following individual fluid elements explicitly during the heating (or cooling) process, by using eq 8. This procedure is fairly general since eq 12 was developed from a fairly arbitrary relationship and did not explicitly use the species mass balance equation (eq 10). The distribution of sterilization or any other time-temperature-dependent quality factor can therefore be obtained without explicitly following individual particles. Numerical Solution of the Equations. Since eq 12 is quite similar to the energy equation, it can be solved by the same technique as the energyequation. The equations were solved by using an alternating direct implicit (ADI) finite difference procedure. To write the AD1 formulation for the sterilization equation 12, it is first combined with
the continuity equation in the cylindrical coordinates
a -i -(ru) r ar
a + -(u) ax
=0
to obtain
The equations for the two steps for AD1 are written as
The boundary conditions used are
Fo(at the boundary) = ~ t l O n ( t ) - l z ' / z dt for time t
> 0 (17)
The initial condition is
Fo = 0 everywhere at t = 0
(18) Applicationto Sterilization of Canned Liquid. This procedure for calculating FOor sterilization was applied to sterilization of liquid food in a vertically placed cylindrical metallic container with a diameter of 0.0838 m and height of 0.1067 m. In absence of any agitation, the liquid would heat by natural convection heating. The energy equation, continuity equation, and sterilization equation (eq 14) were combined with the momentum equation in cylindricalcoordinates. The constant property values used were those of water. Although the Grashoff number was high (=loa) for the properties of water, laminar flow was assumed. The boundary conditions were set at constant temperature of 121 "C on top, bottom, and side walls. The thermal resistance of the walls was ignored. Fluid velocities were set to zero at the walls (no-slip condition). The initial conditions were uniform temperature (30 "C) of fluid at rest. As mentioned earlier, the equations were solved by the AD1 finite differencing method. The temperature and velocities computed by this method (Datta and Teixeira, 1987)had compared very well with experimental results and will not be discussed here. Calculations of Thermal-Time Distribution. To obtain a distribution of thermal times (Fo)from the calculated FOvalues at the nodal points, a density function for volume fraction was defined as (19) where d 9 [=p(Fo)d F o ] is the volume fraction that has a thermal time between FOand FO+ dFo. If p(F0)is plotted against Fo,the area under the curve between two values of FO would be the fraction of fluid volume that has between these values. The above equasterilization (Fo) tion (19) was computationally implemented as
where A t ? k is the summation of volume fraction of liquid having sterilization values within a particular range of
Bbtechnol. Rug..,1991, Vol. 7, No. 5
400
sterilization log Folk and log Fo(K+~. Symbolically Shaded area is the fraction 01 tdal liquidvolume having temperatures between 119.5 and 120.0%
was calculated for every interval by dividing the totalrange of sterilization. The volume (or mass) fraction At9ij associated with a grid point (ij)for a nonequidistant grid was defined for the cylindrical geometry as
-5 .-
0.55
LL
0.50
0
z
.-h l c n
A9ij = (2~rj)(O.5)(rj+l-rj-i)(O=5)(zi+1-zi-i)/V (22) where V is the total volume of the container. The logarithmic transformation of sterilization values (Fo) was needed since they covered several orders of magnitude. In another processing situation, this may not be necessary. The density function for volume fraction p(1og F&+l)was plotted against the upper limit of the sterilization range log F&+l. The temperature values were plotted similarly, except they were not transformed to their logarithms:
The volume (mass) average values of the concentrations are important, particularly when applied to nutrient retention. One can define a mass average FOvalue as
However, such an average value would lack physical interpretation. Instead, the Fo value corresponding to the average concentration of nutrients (or bacteria) can be more easily interpreted. Average concentration E is defined as
by using eq 3. Average FOcorresponding to this average concentration is
c E,(t) = -= - -In ( ce-Fo(')kTR d19) k~~ k~~
(26)
Extent of sterilization at the slowest heating point is normally used to decide the duration of processing. The slowest heating point is, by definition, the lowest temperature value in the domain. If the material at this location is sterilized to the desired extent, material at all other locations are sterilized to a value higher than or equal to the desired value, making the complete product safe from a microbiological point of view. The FO or sterilization at the slowest heating point was calculated as
where P o m t ( t ) is the lowest value of temperature at any time instant among the calculated temperature values. The lowest value of temperature did not always occur at the same location but was restricted in a zone near the bottom of the cylindrical can, which would be referred to as the slowest heating zone.
Results Distribution of Temperature. The velocities, temperatures, and associated sterilization of individual particles change constantly due to the transient nature of the
0.45 1
119.0
119 5
ti0
12b 5
12i.O
Temperature ("C)
Figure 1. Distribution of temperature values of fluid elements within the enclosureafter 10 min of heating. The area under the curve between two values of temperature provides the volume fraction with temperature between those values.
process. The distribution of temperature values within the container after 10 min of heating is shown in Figure 1. From the definition of these distributions (eq 20), the area under the curve between a range of temperature values represents the volume fraction having a temperature in that range. A transient change in the distribution is shown in Figure 2. As expected, the temperatures lie between 30 and 121 "C,which are initial and boundary temperatures, respectively. As the heating of the liquid continues, more of the fluid attains higher temperature and the distribution shifts to the right, staying below or equal to 121 "C. Also, all the temperatures become closer to the boundary temperature, making the range of the distribution progressively smaller. Distribution of Sterilization. The distribution of sterilization as shown in Figure 3 shifts to the right continuously with time since sterilization of a particle (eq 4) is an increasing positive function of temperature. The spread in sterilization values (difference between the largest and smallest values of sterilization) keeps decreasing since the spread in temperature values become progressively smaller. This is due to the presence of the exponential function of temperature in the sterilization equation (eq 4). Initially temperature values have a much wider range, which is accentuated in the sterilization calculation, making it spread over several orders of magnitude. Slowest Heating Particle and Slowest Heating Zone. The sterilization curve for 10 min from Figure 3 is enlarged in Figure 4. The volume (or mass) average sterilization calculated according to eq 26 lies within the range of the distribution, as expected. However, the sterilization calculated on the basis of the temperature at the slowest heating point lies completely outside the distribution. This is because all real fluid elements stay only some of the time in the slowest heating zone and no fluid element stays all of the time in the slowest heating zone. Thus all fluid elements in the system obtain a sterilization more than what is calculated on the basis of the slowest heating zone. For a fluid element to have the least sterilization as calculated from the slowest heating zone, the element would have to be stationary at that location, which is very unlikely. The slowest heating zone sterilization moves relatively closer to the slowest fluid element sterilization with time. However, the absolute magnitude of the differencebetween two sterilizations becomes more significant with time. In Figure 4, for example,after 10min of heating, the difference (1.4 min) between the slowest particle sterilization of 3.40 min and the slowest heating zone sterilization of 2.05 min is quite considerable. Thus, the liquid would normally be
40 1
Bbtechnol. Rug., 1991, Vol. 7, No. 5
-1
er 10 minutes
i- 0.24 50
60
70
so
90
160
li0
li0
T, teryerature ("C)
Figure 2. Distribution of temperature values in the enclosure after various times of heating, showing the expected reduction in the spread of temperature values as heating continues.
-1
Aiier 10 minutes
0.01
0.1 Sterilization. F, (minutes)
Conclusions
1
10
Figure 3. Sterilization values of fluid particles within the enclosure at various times, showing reduction in the range of sterilization values with heating due to decreasing spread in temperature values as explained in Figure 2.
31°] CL
i
d f P
Ellborg (1989) are developed for validation, the concept of slowest heated fluid element (as opposed to location) could be used by the regulatory agencies as well, resulting in possible improvement of quality. Also, these distributions (Figures 1-4) are the most appropriate way to represent the complete information on biochemical changes in processing of fluids since they vary continuously (as opposed to discretely). This procedure for obtaining a distribution of Fo values is limited to the extent that it is not possible to obtain directly the sterilization or the path of a liquid element from a specified starting location. The distribution of sterilization that was obtained includes all elements of the system, but individual fluid elements with individual starting locations cannot be identified on this distribution. Most likely this is of little consequence since one is normally interested in the overalldistribution instead of sterilization of an individual fluid element starting from a given location. Also, for turbulent flow, additional computational aspects would need to be considered.
8
(1)A mathematical procedure was developed for convective heated food materials to calculate distribution of sterilization, nutrient retention, and other quality factors (biochemical changes) by implicitly following the fluid elements during the heating process. The procedure avoids explicit Lagrangian solution technique that might otherwise be necessary to obtain such information by "tagging" fluid particles in motion. The technique developed is quite general and applies to a large class of convective flow problems. (2) By use of this technique, transient distributions of sterilization were obtained for the first time for recirculating natural convection flow of fluid in a cylindrical enclosure being heated from all sides. The results showed that the lowest sterilization achieved by any portion of the fluid in the system can be considerably more than the sterilization normally calculated by using the temperature at the slowest heating zone.
2-
Notation
B
1
10
Sterilization, Fo(minutes)
Figure 4. Sterilization values of fluid particles within the enclosure after 10 min. The calculated slowest point sterilization is shown to be lower than the actual sterilization of any portion of the fluid, including the slowest heated fluid. overprocessed, if the processing time is based on sterilization at the slowest heating point or zone. Since the nutrient degradation would behave similarly to sterilization but with a different value for 2,calculated nutrient degradation based on the slowest heating zone would be underpredicted compared to the actual nutrient degradation of the slowest heated fluid. Comments on the Procedure. The procedure developed can be used to get a distribution of biochemical changes in arbitrary nonisothermal processing situations. Whenever temperature and velocities can be calculated for a processing situation, it is possible to get a distribution of the concentrations of various components. The procedure discriminated sterilization of real particles from the sterilization at the slowest heating zone, showing the extent of overprocessing when the latter is used as a basis. Thus, the distributions obtained by this procedure would be a valuable tool for process feasibility studies. Once sophisticated experimental techniques such as those of
concentration of any component in food material undergoing changes average concentration particle or substantial derivative activation energy equivalent sterilization at 121 OC in minutes at time t average FO values as defined in eqs 24 and 26 reaction rate constant at temperature T frequency factor total number of intervals for data analysis pressure densityfunction for volume fraction of sterilization curves density function for volume fraction of temperature curves distance in radial direction temperature in degrees Celsius time of heating thermal time velocity in vertical direction total volume of the enclosure velocity in radial direction distance in vertical direction
Biotechnol. hog., 1991, Vol. 7, No. 5
402
r CY
At9 Ir
P
value for temperature dependency of reaction = 2.303R T R ~ / E , thermal diffusivity volume fraction viscosity of fluid being heated density of fluid being heated z
Subscripts 0 initial location of a particle 1 initial value 1 grid point in vertical direction grid point in radial direction i k interval grid point on the boundary r boundary R reference value Superscripts slowest slowest heating location in the enclosure n time step n+(l/2) half time step in AD1 method
research was also supported by the Hatch Project 408, U.S.Department of Agriculture.
Literature Cited Batchelor,G. K. An introduction to fluid dynamics; Cambridge University Press: Cambridge, England, 1967. Datta, A. K.Error estimatesfor approximate kinetic parameters used in food literature. J. Food Eng. 1991,in press. Datta, A. K.; Teixeira, A. A. Numerical modeling of natural convection heating in canned liquid foods. Trans. Am. SOC. Agric. Eng. 1987,30 (5),1542-1551. Datta, A. K.; Liu, J. Thermal time distributions for microwave and conventional heating of food; American Society of Agricultural Engineers: St. Joseph, MI, 1990;Paper 90-6600. Ellborg, A. Thermal time distributions in a tube heat exchanger. A study with application to sterilization of liquid foods. Doctoral Dissertation, Division of Food Engineering, Lund University, Lund, Sweden, 1989. Karel,M.; Fennema,0.R.; Lund, D. B. Principles offood science. Physicalprinciples of food preservation;Marcel Dekker,Inc.: New York, 1975. Nauman,E. B. Nonisothermal reactors. Theory and application of thermal time distributions. Chem. Eng. Sci. 1977,32,359367.
Acknowledgment This research was conducted by using the Cornel1 National Supercomputer Facility, a resource of the Center for Theory and Simulations in Science and Engineering (Cornell Theory Center), which receives major funding from the National Science Foundation and IBM Corp., with additional support from the state of New York and members of the Corporate Research Institute. This
Paulsson, B. 0.;Tragardh, C. Calculation of sterilization effect distribution in heat exchangers in turbulent flow. In Engineeringand Food; McKenna,B. M., Ed.; Elsevier Applied Science Publishers: London, 1984;Vol. 1. Tragardh, C.; Paulsson, B. Heat transfer and sterilization in continuous flow heat exchangers. In Developments in Food Preservation-$ Thorne, S., Ed.; Elsevier Applied Science Publishers: London, 1985. Accepted June 21,1991.
