Abstract
A mathematical model was developed to correlate the four heat penetration parameters of 57 Stumbo’s tables (18,513 datasets) in canned food: (the difference between the retort and the coldest point temperatures in the canned food at the end of the heating process), (the ratio of the heating rate index to the sterilizing value), (the temperature change required for the thermal destruction curve to traverse one log cycle), and , (the cooling lag factor). The quantities , , and , are input variables for predicting , while , and are input variables for predicting the value of , which is necessary to calculate the heating process time , at constant retort temperature, using Ball’s formula. The process time calculated using the value obtained from the mathematical model closely followed the time calculated from the tabulated values (root mean square of absolute errors RMS = 0.567 min, average absolute error = 0.421 min with a standard deviation SD = 0.380 min). Because the mathematical model can be used to predict the intermediate values of any combination of inputs, avoiding the storage requirements and the interpolation of 57 Stumbo’s tables, it allows a quick and easy automation of thermal process calculations and to perform these calculations using a spreadsheet.
1. Introduction
Featherstone [1] indicated that the nutritional value of properly processed canned food is as good as that of fresh or frozen food. To assure a safe canned product with minimal damage to organoleptic quality and nutritional value, it needs to optimize the thermal processing of the product through rigorous calculations [2].
The general method was the first method developed for thermal process calculations [3]. The fundamental concepts on which it was based served as foundation for the development of the more sophisticated procedures [4].
Because the general method lacks the predictive power needed for design purpose, the difficulties associated with this procedure inspired interest in the formula method first proposed by Ball [5]. Over the years, Ball’s formula method passed through rigorous evaluations, simplifications, and improvements [6, 7].
Additional formula methods were later developed by Hayakawa [8] and Steele and Board [9], but Smith and Tung [10], in a study on accuracy assessments of the methods, reported that Stumbo’s tables [7] gave the best estimations of process lethality under various conditions.
Although accurate, Stumbo’s method is difficult to computerize because it involves 57 tables, as compared to Ball’s method, where only one table is used. In fact, in computer calculations, during the process control of retorts, a large space is required to store 18,513 data of 57 tables, and an interpolation procedure is required.
For this reason, in the last three decades, mathematical models based on heat transfer in canned food have been developed and applied in thermal processes and their optimization [11–14]. The use of mathematical modeling by computational fluid dynamics (CFD) of the heat transfer for evaluating thermal process has been showing to be a powerful tool to assure food safety and nutritional quality [15–19].
However, these models require multiple input data related to the food product and system, such as the heat transfer coefficient of the heating and cooling medium, thermal diffusivity of the food product, can shape and dimensions, and processing conditions. Only by determining these parameters precisely can the time-temperature history at any specific location in the product be obtained through the solution of a system of nonlinear algebraic equations requiring a time-consuming iterative process [20].
As an alternative to perform the thermal process calculations for canned food, in recent years, artificial neural network (ANN), which requires a sequential solution of a number of algebraic equations, has been proposed to predict the process lethality [20] or to predict the value for a given value, or vice versa, with very good accuracy [21, 22]. In the last case, the ANN technique is an information processing system that learns from Stumbo’s datasets.
In the present work, it was proposed to transform the 57 Stumbo’s tables, which form the basis for the application of the formula method, into a mathematical model easy to handle and to be computerized.
Therefore, the overall objective of this paper was to automate, in a simple and fast way, a new mathematical procedure involved in thermal process calculations. The computerization was obtained by developing an original new mathematical model, based on various modification and improvements of Ball’s model, converging to Stumbo’s datasets, to predict the value for a given value, or vice versa.
2. Review of Canned Food Thermal Process Calculations
Even though the lethal effects of heat processing, similar to those of sterilization, have been empirically known since the beginning of the 19th century, thermal processing has been optimized through several steps over the past two centuries.
The first step involved investigating the kinetics of microbial destruction at a constant temperature. In this regard, Bigelow [23] found that the reduction in microorganisms over time is usually proportional to the number of microorganisms with a constant of proportionality , known as the rate constant (min−1), which depends on the microorganism, the temperature (°C), and the chemical and physical characteristics of the product as follows: . In this equation, is the number of microorganisms at time (min), and corresponds to the initial time . The integral of the previous equation is , which can be conveniently represented by the decimal logarithm as follows: . In this equation, is the decimal reduction time (min) required to destroy 90% of the initial microbial population.
The evaluation of the effect of temperature on the kinetics of microbial destruction was the second step. An increase in the processing temperature produces a decrease in the decimal reduction time (min) in accordance with the Arrhenius law, which, if it is properly applied, will provide . If a decrease of 90% in the decimal reduction time is assumed by increasing the temperature from to , the previous equation yields . The constant is then the increase in temperature that allows the decimal reduction time to decrease 10-fold.
For spore cells, the average of is 10°C, while, for vegetative cells, is lower because of lower temperatures. The experimental values of the decimal reduction time for the various microorganisms are often obtained at a reference temperature of 121.1°C (250°F). Therefore, the values of are known. Moreover, for various microorganisms, the values of the number of decimal reductions necessary to achieve desirable sterilization are also known: . Finally, the total heating time (min) at constant temperature (equal to 121.1°C), which is known as the process lethality, is defined by the symbol and is simply calculated by . When , the process lethality is indicated by the symbol [6].
The third step was to recognize the kinetics of alteration of the constituents (e.g., enzymes, heat-sensitive proteins, and vitamins) at a constant temperature during the thermal process. These alterations usually follow first-order kinetics: , such as those of the thermal death process, where is the concentration of the constituent in its original state. Thus, the same mathematical developments could be made, yielding the same results achieved in the destruction of microorganisms, that is, the definition of a decimal reduction time of the original constituent and the increase in temperature that allows the value to be reduced 10-fold. The parameter now becomes higher, reaching values of approximately 30°C for the B vitamins, of approximately 50°C for vitamin C and of almost 100°C for chlorophyll.
The fourth step was the study of thermal processes at variable temperatures versus time of the canned food. For example, a sterilization by a retort filled with steam condensing at the constant temperature of the retort, , on the surfaces of cans, but with the geometrical center of the can (coldest point) at variable temperature during the heating time.
The change in temperature versus time in the coldest point of the canned food produces a change in the decimal reduction time , according to the following equation:
In an infinitesimal time interval , during which it is possible to assume that the product has a lethal temperature and a corresponding value of at its coldest point, there is an infinitesimal reduction in the number of microorganisms , which corresponds to an infinitesimal increase of the decimal reduction: . Consequently, we can write: , which can be integrated as follows: Rearranging, the process lethality is obtained as follows: To solve the integral and to obtain the process lethality during the heating and cooling phases of canned food, it is necessary to know the relationship between the coldest point temperature and time , also called the heat penetration curve or temperature-time history.
The first procedure used to calculate was developed by Bigelow and is usually known as the general method [3]. The general method makes direct use of the graphical temperature-time history at the coldest point. Consequently, Bigelow’s procedure involves the choice of a small time interval that can be plotted on the graph of the heat penetration curve. At each time interval, the average temperature is read, and the decimal reduction time can be calculated.
The procedure ends with the graphical construction of the function versus time and measurement of the area under this function, which yields the decimal reduction number and the process lethality .
If higher or lower values of lethality are required, the procedure can be repeated and the estimate of the cooling curve advanced or retarded on a trial-and-error basis until the desired lethality is achieved.
Bigelow’s procedure was termed the general method because it depends only on measurement of the temperature-time history of the coldest point. Because it does not consider the mode of heat transfer, the food properties, or the can size and shape, it is very accurate. However, the processing time calculated using this method is specific for a given set of processing conditions, and the method is not suitable for computer applications because it is very long, tedious, and impractical.
Over time, several improvements to the original general method have been introduced by Ball [24], Schultz and Olson [25], Patashnick [26], Hayakawa [27], and, more recently, Simpson et al. [28].
The difficulties associated with Bigelow’s original procedure prompted interest in the formula method first proposed by Ball [5]. The formula method highlights that the difference of temperature between the retort and coldest point of the can decays exponentially over the process time after an eventual initial lag period as follows: , where is the heating rate lag factor at the can center (coldest point), (°C) is the retort temperature, and (°C) is the initial food temperature and (min) is the heating rate index (the heating time required for the log temperature versus time plot to traverse one log cycle). Thus, the heating process time (min) at a constant retort temperature is given by Ball and Olson [6] and Stumbo [7] as follows: Equation (4) is known as Ball’s formula, where (°C) is the difference between the retort temperature and coldest point temperature (°C) at the end of the heating process. Because (4) requires knowledge of the factor and the index through experimental evaluation, the formula method is considered semianalytical. About the lag factors, they may assume different values depending on the rheological and thermal properties of the product.
It may have fluid foods, at low viscosity, for which the convective heat transfer inside the can is prevalent. The product has a changeable temperature versus time, but the temperature tends to be uniform across all parts of the can. In this case, the temperature differences during the heating process and the cooling process , where is the cold-water temperature, follow the exponential equation from the beginning, and, thus, the heating and cooling rate lag factors at the coldest point, and , respectively, are approximately 1. In practice, they may be less than 1 because of changes in the internal convective coefficient [29].
The product can be solid, in which case the internal heat transfer takes place by conduction. The high thermal resistance, which is associated with conduction, impairs heat penetration, and, therefore, the temperature changes not only during the time, but also between different parts of the canned food with the generation of a coldest point. During heating, the coldest point receives the thermal wave last. In this case, the temperature of the coldest point follows the exponential equation only after some time in the beginning of the heating at a constant retort temperature. Thus, the heating rate lag factor, , and cooling rate lag factor, , respectively, are equal to approximately 2.
Finally, the product may consist of solid pieces immersed in a liquid. In this case, the heat transfer within the product becomes more complex with contributions from both convection and conduction; the values of the heating and cooling rate lag factors, at the coldest point, are between 1 and 2.
The sterilizing value can be defined as the time required at retort temperature to accomplish a heat process of some given value: . Ball [5] discovered that, for a single value of , any given ratio has a value of that corresponds to it. In computing values of , by the mathematical resolution of the integral (3), Ball [5] considered the cooling lag factor to be 1.41 and, thus, produced one table of versus for different values of (6–26°F). Over the years, Ball’s method, consisting of the determination of for a given value and the resolution of formula (4), has undergone rigorous evaluations, simplifications, and improvements [6].
It was later found that, though the relationship was virtually independent of process parameters such as the retort temperature, initial food temperature, cooling water temperature, can dimensions, and heating lag factor , it was greatly dependent upon the lag factor of the food cooling curve.
In recognition of this dependency, Stumbo [7] developed 57 tables of that covered a wide range of cooling lag factor values (0.4–2) and -values (4.4°–111.1°C; 8°–200°F), but under the condition that the cooling rate index . Because Stumbo used the general method of integration to obtain the values (including the lethal value of heat during cooling) in the relationship, it follows that for any given value of , the value of in the ratio accounts for all lethal heat during both heating and cooling. Then, as used in formula (4) of the heating curve also accounts for all lethal heat of both heating and cooling.
3. Mathematical Model
As mentioned in the previous paragraph, the semianalytical methods of the formula are based on the integration of (3). Thus, the relationship between temperature and time must be defined.
During the heating process, starting from the time when the retort temperature is considered constant, Ball [6] proposed the following rearrangement of (4): Combining (5) with (3) and with integration, the result of Ball and Olson [6] is obtained as follows: where is the lethality during the heating process and the function Ei is called the exponential integral. The values of this function can be easily calculated with the following series: The lower limit of the integration, that is, the initial temperature difference of the heating process, was imposed by Ball and Olson [6] and is equal to 44.4°C (80°F).
Considering that the retort temperature is usually not greater than approximately 140°C (284°F), the initial temperature of the coldest point , considered for the calculation of , has a value lower than 100°C. For any value, a temperature lower than 100°C produces a negligible contribution to lethality when is lower than approximately 15°C (26°F) [4]. In fact, in this case, , and then, with no significant error, the equation of can be simplified as follows: Instead, if is higher than 15°C (26°F) and with higher values, this simplification produces an error that can be eliminated by the introduction of a correction factor , as shown in (27).
As was noted in the introduction, during cooling, Ball and Olson [6] considered only the case of and, therefore, the presence of a cooling lag. This limitation means that it is possible to represent the cooling of the coldest point with a function similar to (5), only with a lag behind the introduction of cold water into the retort, at the temperature .
During this cooling lag, Ball [5] proposed to represent the temperature versus time using a hyperbolic function. The characterization of this hyperbola was proposed by Ball on the basis of various empirical observations; however, it is valid when is between 3.33° and 14.4°C (6°–26°F) and .
As was discussed in the introduction, Stumbo [7] obtained the relationship for values ranging from 0.4 to 2 and for values ranging from 4.44° to 111.1°C (8°–200°F) using the general method with computational methods, replacing the graphical integration of Bigelow.
Therefore, it is difficult to extend Ball’s hyperbolic function of temperature versus time in the cooling lag over such large Stumbo’s conditions.
Here, a parabolic relationship is proposed: or where is the cooling time ().
Equation (10) was preliminarily tested, with good results, using data presented by Simpson et al. [28].
To obtain the parameters , , and , some conditions are placed as follows.
(1) When , the following derivative of the parabolic cooling lag equation (10): must be equal to the derivative of the heating equation (5): Thus, also considering (4), is obtained as follows: (2) The temperature of the coldest point at the end of the heating process must be the same temperature as that at the start of the parabolic equation (when ). This allows for calculation of the parameter : (3) When the cooling lag factor is equal to 2 and, therefore, the heating of the product occurs by pure thermal conduction, the solution of Fourier differential equation is an infinite series of exponential and Bessel functions, which asymptotically tend to a simple exponential equation.
The use of this last equation, instead of an infinite series, is effective when the dimensionless time, called the Fourier number , is greater than [25], where is the thermal diffusivity of the product, is a characteristic dimension (shape factor) defined here as that can be used for any can geometry, is the volume, and is the surface area of the can.
By defining the corresponding cooling time: and the cooling exponential equation in a form similar to (5) as well as the typical form used in the heat transfer literature, the following identity is found: where, for and cylindrical cans, is given by Mafart [29]: For cans typically used in the food industry [28], (17) yields: .
For a better approximation, in the present mathematical model, is chosen.
Now, it is possible to establish that the last point of the cooling parabolic equation coincides with the first point of the cooling exponential equation at time :
Then, is obtained as follows: By combining (10) with (3) and with integration, the cooling process lethality is achieved as follows: As can be seen later, for convergence of results of the mathematical model to data of Stumbo’s tables showing the relationship, it is necessary to introduce a correction factor.
Because the mathematical model that is being built should be extended to all values of rather than to only value of 2, this correction factor will surely be a function of but will also depend on and on .
Through the elimination of the error function in (20), a simplification of the same equation (20) and a less problematic influence of on the correction factor are achieved: With regards to the possible contribution to the lethality of the subsequent cooling described by the exponential equation, the calculation of the temperature of the coldest point at time through (10) showed values that were always lower than 100°C for each of the and values of Stumbo’s tables.
As seen for (8), when is lower than 15°C (26°F), these low temperatures produce a negligible contribution to lethality, and, therefore, it is possible to exclude the calculation of lethality due to the cooling exponential equation. Thus, in this mathematical model, the total lethality is only due to the sum of and , given by (8) and (21), respectively.
When is higher than 15°C, the error caused by these exclusions (cooling exponential equation and factor in (8) will be eliminated by the introduction of the correction factor , defined through the next equation (27), which is in any case necessary for the mathematical model to fit to Stumbo’s datasets.
As mentioned before, the cooling parabolic curve was introduced to describe the temperature-time history during the lag under the condition of pure thermal conduction and, therefore, with . Now, for example, if , the time of lag becomes negative in the cooling time scale . Virtually, it is as if the end point of the heating would come back. This setback must necessarily be related to the function (Figure 1).
Now, to extend the mathematical model to all the values of (from 0.4 to 2), it is assumed that the argument of the logarithm is not simply but a function : where is also reasonably correlated with , as we can infer by examining Stumbo’s datasets.
Thus, the values virtually become as follows: Considering, as it was proposed by Stumbo [7], that and by combining (5) with the last equation (23), the relationship is the following: Thus, multiplying by the quantity , included in (8) of and (21) of , directly or indirectly through and and by imposing equal to of Stumbo’s tables for the various and values, the values of the correction factor were obtained (Figure 2).
Figure 2 suggests a power relationship as follows: Using a multiple nonlinear regression, and functions were obtained:
Figure 2 shows that when , is independent of and is approximately 0.89 instead of 1. This was due to the choice of neglecting the error function, transforming (20) into (21).
In any case, that is, for each value of and , was almost constant for medium and high values of but tended to increase when tended to 0. Thus, the introduction into (25) of an exponential factor, found by means of a multiple regression, was necessary.
Therefore, a more complete correction factor was defined as follows: must be used as a multiplier of the quantity present, directly or indirectly through and , in (8) of and (21) of , which are useful for calculating :
4. Results and Discussion
4.1. Calculation
The mathematical model, made of the following ten equations (26) (27), (8), (7) (with and ), (13), (14), (16) (with and ), (19), (21), and (28) to be sequentially solved, was implemented in a spreadsheet to calculate values by varying values obtained from Stumbo’s tables. The comparison between the calculated value of and desired Stumbo’s values of is shown in Figure 3. The following were found: a determination coefficient , a mean relative error , and a mean absolute error .
4.2. Calculation
The same mathematical model made of the ten equations (26), (27), (8), (7), (13), (14), (16), (19), (21), and (28) was easily implemented in the spreadsheet, to calculate values for various values obtained from Stumbo’s tables. The system of equations is implicit in the value; therefore, can be rapidly calculated using the spreadsheet by putting Stumbo’s values into (28) and searching for the corresponding value that ensures the validity of the ten equations. A comparison between the calculated value of and desired Stumbo’s values of is shown in Figure 4. The following values were found: determination coefficient , mean relative error , and mean absolute error .
4.3. Mathematical Model Validation
In order to validate the mathematical model of this work, thermal process time was calculated for various heating conditions. For comparison purposes, these conditions were the same as those imposed by Sablani and Shayya [21]: (10° and 44.4°C; 18° and 80°F), (65.5°C; 150°F), (30 and 90 min), (1 and 2), (0.4 and 2), (111.1°, 121.1°, and 140°C; 232°, 250°, and 284°F), and (5, 15, and 25 min).
To calculate the thermal process time reference, by values from Stumbo’s tables [7], the Ball’s formula (4) was used. To improve the validation of the mathematical model of this work, it was instead used a modified Ball’s formula by a corrective coefficient, which was found dependent on and as follows: . Therefore the modified Ball’s formula was Figure 5 and Table 1 provide a comparison of the process time calculated using values from the mathematical model described in this work and formula (29), and values from Stumbo’s tables [7] and formula (4).
Table 1 shows also that the process time calculated using this mathematical model is closer to the values calculated using the values from Stumbo’s tables [7] than the values calculated using both the neural networks models of Sablani and Shayya [21] and Mittal and Zhang [22]. Compared to the latter, the average absolute error was reduced from 0.466 to 0.421 min. The root mean square of absolute errors value RMS was also reduced from 0.612 to 0.567, with a SD reduced from 0.400 to 0.380 min. All of these comparisons were significant at a 95% confidence level.
5. Conclusions
A new original mathematical model was developed to calculate and starting from Ball’s formula method of thermal process calculations. The new model was based on ten equations to be sequentially solved, in an easy and fast way, even on a spreadsheet to obtain , and on the modified Ball’s formula (29) to calculate the thermal process time .
The ten equations are based on various modification and improvements of Ball’s mathematical model, such that the new model results converge on those values obtained from Stumbo’s tables [7].
The mathematical model of this work allows the computerization of 18,513 Stumbo’s data, avoiding the use of 57 lookup tables and interpolation among the three data parameters associated with the use of these tables.
The thermal process time, calculated using the from the mathematical model presented in this work and using the modified Ball’s formula (29) here proposed, closely followed the time calculated from the tabulated values ( min, average absolute min with an of 0.380 min). This high prediction accuracy, better than ANN models, allows practical applications, such as an easy evaluation of the parameters for thermal processes controlled by computers.