# A New Mathematical Model for Food  Thermal Process Prediction

Authors: Dario Friso
Journal: Modelling and Simulation in Engineering
 (2013)

---

## Abstract

A mathematical model was developed to correlate the four heat penetration parameters of 57 Stumbo’s tables (18,513 datasets) in canned food:g (the difference between the retort and the coldest point temperatures in the canned food at the end of the heating process), fh/U (the ratio of the heating rate index to the sterilizing value), z (the temperature change required for the thermal destruction curve to traverse one log cycle), and Jcc, (the cooling lag factor). The quantities g, z, and Jcc, are input variables for predicting fh/U, while z, Jcc and  fh/U are input variables for predicting the value of g, which is necessary to calculate the heating process time B, at constant retort temperature, using Ball’s formula. The process time calculated using the g value obtained from the mathematical model closely followed the time calculated from the tabulated g 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.


---

## Body

## 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 lethalityF [20] or to predict the g value for a given fh/U 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 theg value for a given fh/U 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 kT, known as the rate constant (min−1), which depends on the microorganism, the temperature T (°C), and the chemical and physical characteristics of the product as follows: (dN/dt)=-kT·N. In this equation, N is the number of microorganisms at time t (min), and N0 corresponds to the initial time t0=0. The integral of the previous equation is ln(N0/N)=kT·t, which can be conveniently represented by the decimal logarithm as follows: log(N0/N)≅(kT/2.3)·t=t/DT. In this equation, DT=2.3/kT 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 timeDT (min) in accordance with the Arrhenius law, which, if it is properly applied, will provide log(DT1/DT2)=(T2-T1)/z. If a decrease of 90% in the decimal reduction time is assumed by increasing the temperature from T1 to T2, the previous equation yields z=T2-T1. The constant z is then the increase in temperature that allows the decimal reduction time DT to decrease 10-fold.For spore cells, the average ofz is 10°C, while, for vegetative cells, z is lower because of lower temperatures. The experimental values of the decimal reduction time DT for the various microorganisms are often obtained at a reference temperature of 121.1°C (250°F). Therefore, the values of D121.1 are known. Moreover, for various microorganisms, the values of the number of decimal reductions n necessary to achieve desirable sterilization are also known: n=log(N0/N). 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 F and is simply calculated by t=F=n·D121,1. When z=10°C, the process lethality is indicated by the F0 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:dC/dt=-kT′·C, such as those of the thermal death process, where C 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 DT of the original constituent and the increase z in temperature that allows the DT value to be reduced 10-fold. The parameter z 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,TR, 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 temperatureT versus time t in the coldest point of the canned food produces a change in the decimal reduction time DT, according to the following equation:
(1)log(DTD121.1)=121.1-Tz⟺DT=D121.1·10(121.1-T)/z.In an infinitesimal time intervaldt, during which it is possible to assume that the product has a lethal temperature T and a corresponding value of DT at its coldest point, there is an infinitesimal reduction in the number of microorganisms dN, which corresponds to an infinitesimal increase of the decimal reduction: dn=log(N/(N-dN)). Consequently, we can write: dt=dn·DT, which can be integrated as follows:
(2)n=∫0ndn=∫0tdtDT=∫0tdtD121.1·10(121.1-T)/z.
Rearranging, the process lethality F is obtained as follows:
(3)F=n·D121.1=∫0t10(T-121.1)/zdt.
To solve the integral and to obtain the process lethality F during the heating and cooling phases of canned food, it is necessary to know the relationship between the coldest point temperature T and time t, also called the heat penetration curve or temperature-time history.The first procedure used to calculateF 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 Δt 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 DT can be calculated.The procedure ends with the graphical construction of the function1/DT versus time and measurement of the area under this function, which yields the decimal reduction number n and the process lethality F=n·D121.1.If higher or lower values of lethalityF 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 (TR-T) decays exponentially over the process time t after an eventual initial lag period as follows: (TR-T)=(TR-T0)·Jch·e-(2.3·t)/fh, where Jch is the heating rate lag factor at the can center (coldest point), TR (°C) is the retort temperature, and T0 (°C) is the initial food temperature and fh (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 B (min) at a constant retort temperature TR is given by Ball and Olson [6] and Stumbo [7] as follows:
(4)B=fh·log[Jch(TR-T0)g].
Equation (4) is known as Ball’s formula, where g=(TR-Tg) (°C) is the difference between the retort temperature and coldest point temperature Tg (°C) at the end of the heating process. Because (4) requires knowledge of the factor Jch and the index fh 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(TR-T) and the cooling process (T-TW), where TW 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, Jch and Jcc, 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,Jch, and cooling rate lag factor, Jcc, 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 valueU can be defined as the time required at retort temperature TR to accomplish a heat process of some given F value: U=F·10(121.1-TR)/z. Ball [5] discovered that, for a single value of z, any given ratio fh/U has a value of g that corresponds to it. In computing values of fh/U:g, by the mathematical resolution of the integral (3), Ball [5] considered the cooling lag factor Jcc to be 1.41 and, thus, produced one table of fh/U versus g for different values of z (6–26°F). Over the years, Ball’s method, consisting of the determination of g for a given fh/U value and the resolution of formula (4), has undergone rigorous evaluations, simplifications, and improvements [6].It was later found that, though the relationshipfh/U:g was virtually independent of process parameters such as the retort temperature, initial food temperature, cooling water temperature, can dimensions, and heating lag factor Jch, it was greatly dependent upon the Jcc lag factor of the food cooling curve.In recognition of this dependency, Stumbo [7] developed 57 tables of fh/U:g that covered a wide range of cooling lag factor Jcc values (0.4–2) and z-values (4.4°–111.1°C; 8°–200°F), but under the condition that the cooling rate index fc=fh. Because Stumbo used the general method of integration to obtain the U values (including the lethal value of heat during cooling) in the fh/U:g relationship, it follows that for any given value of g, the value of U in the fh/U:g ratio accounts for all lethal heat during both heating and cooling. Then, g 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):
(5)g(TR-T0)=(TR-T)(TR-T0)=Jch·exp(-2.3·tfh),T=TR-(TR-T0)·Jch·exp(-2.3·tfh).
Combining (5) with (3) and with integration, the result of Ball and Olson [6] is obtained as follows:
(6)Fh=-fh2.3·10(TR-121.1)/z×[Ei(-2.3·gz)-Ei(-2.3·44.4z)],
where Fh 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:
(7)-Ei(-x)=0.5772-ln(x)+x-x22·2!+x33·3!-⋯+xpp·p!.
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 temperatureTR is usually not greater than approximately 140°C (284°F), the initial temperature of the coldest point T, considered for the calculation of Fh, has a value lower than 100°C. For any g value, a temperature T lower than 100°C produces a negligible contribution to lethality when z is lower than approximately 15°C (26°F) [4]. In fact, in this case, Ei(-(2.3·44.4)/z)→0, and then, with no significant error, the equation of Fh can be simplified as follows:
(8)Fh=-fh2.3·10(TR-121.1)/z[Ei(-2.3·gz)].
Instead, if z is higher than 15°C (26°F) and with higher g values, this simplification produces an error that can be eliminated by the introduction of a correction factor M, as shown in (27).As was noted in the introduction, during cooling, Ball and Olson [6] considered only the case of Jcc=1.41 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 Tw.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 z is between 3.33° and 14.4°C (6°–26°F) and Jcc=1.41.As was discussed in the introduction, Stumbo [7] obtained the fh/U:g relationship for Jcc values ranging from 0.4 to 2 and for z 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:(9)Tw-TTw-Tg=a·tw2+b·tw+c
or
(10)T=Tw+(Tg-Tw)(a·tw2+b·tw+c),
where tw is the cooling time (tw=t-B).Equation (10) was preliminarily tested, with good results, using data presented by Simpson et al. [28].To obtain the parametersa, b, and c, some conditions are placed as follows.(1) Whentw=0, the following derivative of the parabolic cooling lag equation (10):
(11)dTdtw=(2a·tw+b)(Tg-Tw)
must be equal to the derivative of the heating equation (5):
(12)dTdtw=(TR-T0)2.3fh·Jch·exp(-2.3fh(tw+B)).
Thus, also considering (4), b is obtained as follows:
(13)b=2.3fh·Jch·exp(-2.3·Bfh)·(TR-T0)(Tg-Tw)=2.3fh·g(Tg-Tw).
(2) The temperature of the coldest point at the end of the heating process Tg=TR-g must be the same temperature as that at the start of the parabolic equation (when tw=0). This allows for calculation of the parameter c:
(14)Tw-TgTw-Tg=a·tw2+b·tw+c⟹c=1.
(3) When the cooling lag factor Jcc 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 numberFo=(α·tw)/l2, is greater than 0.2=Fo* [25], where α is the thermal diffusivity of the product, l is a characteristic dimension (shape factor) defined here as l=2V/A that can be used for any can geometry, V is the volume, and A is the surface area of the can.By defining the corresponding cooling time:tw*=(Fo*l2)/α and the cooling exponential equation in a form similar to (5) as well as the typical form used in the heat transfer literature,
(15)Tw-TTw-Tg=Jcc·exp(-2.3·twfc)=Jcc·exp(-s·Fo),
the following identity is found:
(16)tw*=fc·s·Fo*2.3,
where, for Jcc=2 and cylindrical cans, s is given by Mafart [29]:
(17)s=(5.783R2+π2H2)l2.
For cans typically used in the food industry [28], (17) yields: s=3.8±0.15.For a better approximation, in the present mathematical model,Fo*=0.25 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 timetw*:
(18)a·tw*2+b·tw*+c=Jcc·exp(-2.3·tw*fc).Then,a is obtained as follows:
(19)a=1tw*2·[Jcc·exp(-2.3·tw*fc)-b·tw*-1].
By combining (10) with (3) and with integration, the cooling process lethality Fc is achieved as follows:
(20)Fc=0.584za·(Tw-Tg)·10(b2(Tw-Tg))/(4·a·z)·10(Tg-121.1)/z·erf[0.759·(b+2·a·tw)(Tw-Tg)a·z]0tw*.
As can be seen later, for convergence of results of the mathematical model to data of Stumbo’s tables showing the fh/U:g relationship, it is necessary to introduce a correction factor.Because the mathematical model that is being built should be extended to all values ofJcc rather than to only value of 2, this correction factor will surely be a function of Jcc but will also depend on z and on g.Through the elimination of the error function in (20), a simplification of the same equation (20) and a less problematic influence of g on the correction factor are achieved:
(21)Fc=0.584za·(Tw-Tg)·10(b2(Tw-Tg))/(4·a·z)·10(Tg-121.1)/z.
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 T* at time tw* through (10) showed values that were always lower than 100°C for each of the Jcc and g values of Stumbo’s tables.As seen for (8), when z is lower than 15°C (26°F), these low temperatures T* 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 F is only due to the sum of Fh and Fc, given by (8) and (21), respectively.Whenz is higher than 15°C, the error caused by these exclusions (cooling exponential equation and Ei(-2.3·44.4/z) factor in (8) will be eliminated by the introduction of the correction factor M, 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, withJcc=2. Now, for example, if Jcc<1, the time of lag becomes negative in the cooling time scale tw. Virtually, it is as if the end point of the heating would come back. This setback must necessarily be related to the fc·log(Jcc) function (Figure 1).Figure 1
WhenJcc is less than 1, the time of lag becomes virtually negative in the cooling time scale tw. The dashed line shows the intersection of the cooling exponential curve with the time axis at the negative virtual time: fc·log(Jcc).Now, to extend the mathematical model to all the values ofJcc (from 0.4 to 2), it is assumed that the argument of the logarithm is not simply Jcc but a function M′:
(22)M′=M′(Jcc,z),
where M′ is also reasonably correlated with z, as we can infer by examining Stumbo’s datasets.Thus, theg values virtually become g′ as follows:
(23)g′=(TR-T0)·Jch·exp(-2.3B-fc·logM′fh).
Considering, as it was proposed by Stumbo [7], that fc=fh and by combining (5) with the last equation (23), the g′:g relationship is the following:
(24)g′=g·M′(Jcc,z).
Thus, multiplying M′ by the quantity g, included in (8) of Fh and (21) of Fc, directly or indirectly through Tg and b and by imposing fh/((Fh+Fc)·10(121.1-TR)/z) equal to fh/U of Stumbo’s tables for the various g and z values, the values of the correction factor M′ were obtained (Figure 2).Figure 2
Correction functionM′ versus cooling rate lag factor at the coldest point Jcc and versus z (°C) value. The solid lines are power functions obtained by nonlinear regression.Figure2 suggests a M′(Jcc,z) power relationship as follows:
(25)M′=K(z)·Jccm(z).
Using a multiple nonlinear regression, K(z) and m(z) functions were obtained:
(26)K(z)=exp(z0.7590),m(z)=-z0.5414.5.Figure2 shows that when Jcc=2, M′ is independent of z 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 ofz and Jcc, M′ was almost constant for medium and high values of g but tended to increase when g 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 factorM was defined as follows:
(27)M=K(z)·Jccm(z)·[1+0.1·Jccz/100·exp(-3·gJcc3)].M must be used as a multiplier of the g quantity present, directly or indirectly through Tg and b, in (8) of Fh and (21) of Fc, which are useful for calculating fh/U:
(28)fhU=fh(Fh+Fc)·10(121.1-TR)/z.


## 4. Results and Discussion

### 4.1.fh/U Calculation

The mathematical model, made of the following ten equations (26) (27), (8), (7) (with x=2.3·g/z and p=24), (13), (14), (16) (with Fo*=0.25 and s=3.8), (19), (21), and (28) to be sequentially solved, was implemented in a spreadsheet to calculate fh/U values by varying g values obtained from Stumbo’s tables. The comparison between the calculated value of fh/U and desired Stumbo’s values of fh/U is shown in Figure 3. The following were found: a determination coefficient R2=0.9823, a mean relative error MRE=2.46±3.38%, and a mean absolute error MAE=3.38±20.49.Figure 3
Predictedfh/U values using the mathematical model versus desired Stumbo’s values of fh/U.


### 4.2.g 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 g values for various fh/U values obtained from Stumbo’s tables. The system of equations is implicit in the g value; therefore, g can be rapidly calculated using the spreadsheet by putting Stumbo’s fh/U values into (28) and searching for the corresponding g value that ensures the validity of the ten equations. A comparison between the calculated value of g and desired Stumbo’s values of g is shown in Figure 4. The following values were found: determination coefficient R2=0.9990, mean relative error MRE=2.67±2.69%, and mean absolute error MAE=0.16°±0.20°C  (0.29°±0.36°F).Figure 4
Values ofg predicted using the mathematical model versus desired Stumbo’s values of g.


### 4.3. Mathematical Model Validation

In order to validate the mathematical model of this work, thermal process timeB was calculated for various heating conditions. For comparison purposes, these conditions were the same as those imposed by Sablani and Shayya [21]: z (10° and 44.4°C; 18° and 80°F), T0 (65.5°C; 150°F), fh (30 and 90 min), Jch (1 and 2), Jcc (0.4 and 2), TR (111.1°, 121.1°, and 140°C; 232°, 250°, and 284°F), and F (5, 15, and 25 min).To calculate the thermal process timeB reference, by g 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 TR and Jcc as follows: (123/TR)0.1·Jcc. Therefore the modified Ball’s formula was
(29)B=(123TR)0.1·Jccfh·log[Jch·(TR-T0)g].
Figure 5 and Table 1 provide a comparison of the process time B calculated using g values from the mathematical model described in this work and formula (29), and g values from Stumbo’s tables [7] and formula (4).Table 1
Comparison of process time calculated usingg values from the Stumbo [7] tables, ANN model of Sablani and Shayya [21], ANN model of Mittal and Zhang [22], and mathematical model of this work.
RMS (min)
Deviations or absolute error (min)
SD
Slope of line in Figure5
Max
Min
Ave
Sablani and Shayya [21]
1.164
2.57
0.05
0.936
0.697
1.011
Mittal and Zhang [22]
0.612
1.63
0.01
0.466
0.400
0.999
This work
0.567
1.61
0.003
0.421
0.380
1.001
RMS: root mean square of deviations; SD: standard deviations.Figure 5
Comparison of the process time, calculated usingg values from Stumbo’s tables [7], and the mathematical model of this work.Table1 shows also that the process time B calculated using this mathematical model is closer to the B values calculated using the g values from Stumbo’s tables [7] than the B 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.


## 4.1.fh/U Calculation

The mathematical model, made of the following ten equations (26) (27), (8), (7) (with x=2.3·g/z and p=24), (13), (14), (16) (with Fo*=0.25 and s=3.8), (19), (21), and (28) to be sequentially solved, was implemented in a spreadsheet to calculate fh/U values by varying g values obtained from Stumbo’s tables. The comparison between the calculated value of fh/U and desired Stumbo’s values of fh/U is shown in Figure 3. The following were found: a determination coefficient R2=0.9823, a mean relative error MRE=2.46±3.38%, and a mean absolute error MAE=3.38±20.49.Figure 3
Predictedfh/U values using the mathematical model versus desired Stumbo’s values of fh/U.


## 4.2.g 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 g values for various fh/U values obtained from Stumbo’s tables. The system of equations is implicit in the g value; therefore, g can be rapidly calculated using the spreadsheet by putting Stumbo’s fh/U values into (28) and searching for the corresponding g value that ensures the validity of the ten equations. A comparison between the calculated value of g and desired Stumbo’s values of g is shown in Figure 4. The following values were found: determination coefficient R2=0.9990, mean relative error MRE=2.67±2.69%, and mean absolute error MAE=0.16°±0.20°C  (0.29°±0.36°F).Figure 4
Values ofg predicted using the mathematical model versus desired Stumbo’s values of g.


## 4.3. Mathematical Model Validation

In order to validate the mathematical model of this work, thermal process timeB was calculated for various heating conditions. For comparison purposes, these conditions were the same as those imposed by Sablani and Shayya [21]: z (10° and 44.4°C; 18° and 80°F), T0 (65.5°C; 150°F), fh (30 and 90 min), Jch (1 and 2), Jcc (0.4 and 2), TR (111.1°, 121.1°, and 140°C; 232°, 250°, and 284°F), and F (5, 15, and 25 min).To calculate the thermal process timeB reference, by g 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 TR and Jcc as follows: (123/TR)0.1·Jcc. Therefore the modified Ball’s formula was
(29)B=(123TR)0.1·Jccfh·log[Jch·(TR-T0)g].
Figure 5 and Table 1 provide a comparison of the process time B calculated using g values from the mathematical model described in this work and formula (29), and g values from Stumbo’s tables [7] and formula (4).Table 1
Comparison of process time calculated usingg values from the Stumbo [7] tables, ANN model of Sablani and Shayya [21], ANN model of Mittal and Zhang [22], and mathematical model of this work.
RMS (min)
Deviations or absolute error (min)
SD
Slope of line in Figure5
Max
Min
Ave
Sablani and Shayya [21]
1.164
2.57
0.05
0.936
0.697
1.011
Mittal and Zhang [22]
0.612
1.63
0.01
0.466
0.400
0.999
This work
0.567
1.61
0.003
0.421
0.380
1.001
RMS: root mean square of deviations; SD: standard deviations.Figure 5
Comparison of the process time, calculated usingg values from Stumbo’s tables [7], and the mathematical model of this work.Table1 shows also that the process time B calculated using this mathematical model is closer to the B values calculated using the g values from Stumbo’s tables [7] than the B 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 calculatefh/U and g 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 g, and on the modified Ball’s formula (29) to calculate the thermal process time B.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 theg 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 g values (RMS=0.567 min, average absolute error=0.421 min with an SD 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.



---

*Source: 569473-2013-12-22.xml*