UFDC Home  UF Institutional Repository  UF Theses & Dissertations  Internet Archive   Help 
Material Information
Record Information

Table of Contents 
Title Page
Page i Page ia Acknowledgement Page ii Page iii Table of Contents Page iv Page v Abstract Page vi Page vii Chapter 1. Introduction Page 1 Page 2 Page 3 Chapter 2. Kinetic parameter estimation in foods that heat uniformly Page 4 Page 5 Page 6 Page 7 Page 8 Page 9 Page 10 Page 11 Page 12 Page 13 Page 14 Page 15 Page 16 Page 17 Page 18 Page 19 Page 20 Page 21 Page 22 Page 23 Page 24 Page 25 Page 26 Page 27 Page 28 Page 29 Page 30 Page 31 Page 32 Page 33 Page 34 Page 35 Page 36 Page 37 Page 38 Page 39 Chapter 3. Heat transfer simulation for thermal process design Page 40 Page 41 Page 42 Page 43 Page 44 Page 45 Page 46 Page 47 Page 48 Page 49 Page 50 Page 51 Page 52 Page 53 Page 54 Page 55 Page 56 Page 57 Page 58 Page 59 Page 60 Page 61 Page 62 Page 63 Page 64 Page 65 Page 66 Page 67 Page 68 Page 69 Page 70 Chapter 4. Kinetic parameter estimation in foods that do not heat uniformly Page 71 Page 72 Page 73 Page 74 Page 75 Page 76 Page 77 Page 78 Page 79 Page 80 Page 81 Page 82 Page 83 Page 84 Page 85 Page 86 Page 87 Page 88 Page 89 Page 90 Page 91 Page 92 Page 93 Page 94 Page 95 Page 96 Page 97 Page 98 Page 99 Page 100 Page 101 Page 102 Page 103 Page 104 Chapter 5. Exploiting treatments that sensitize bacterial spores to heat in thermal process design Page 105 Page 106 Page 107 Page 108 Page 109 Page 110 Page 111 Page 112 Page 113 Page 114 Page 115 Page 116 Page 117 Page 118 Page 119 Page 120 Page 121 Page 122 Page 123 Page 124 Page 125 Page 126 Page 127 Page 128 Page 129 Page 130 Chapter 6. Summary and conclusions Page 131 Page 132 Page 133 Page 134 Appendix A. Invalid assumption of EPM Page 135 Page 136 Page 137 Page 138 Appendix B. Heat transfer codefinite cylinder Page 139 Page 140 Page 141 Page 142 Page 143 Page 144 Page 145 Page 146 Page 147 Page 148 Page 149 Page 150 Page 151 Page 152 Page 153 Page 154 Page 155 Page 156 Page 157 Page 158 Page 159 Page 160 Page 161 Page 162 Page 163 Page 164 Page 165 Page 166 Page 167 Appendix C. PEIE code for conduction heating canned food Page 168 Page 169 Page 170 Page 171 Page 172 Page 173 Page 174 Page 175 Page 176 Page 177 Page 178 Page 179 Page 180 Page 181 Page 182 Page 183 Page 184 Page 185 Page 186 Page 187 Page 188 Page 189 Appendix D. Microbial enumeration  0.0 kGy Page 190 Page 191 Appendix E. Microbial enumeration  1.0 kGy Page 192 Page 193 Appendix F. Microbial enumeration  3.0 kGy Page 194 Page 195 Page 196 Appendix G. Process pair data  0.0 kGy Page 197 Page 198 Appendix H. Process pair data  1.0 kGy Page 199 Page 200 Page 201 Page 202 Appendix I. Process pair data  3.0 kGy Page 203 Page 204 Page 205 References Page 206 Page 207 Page 208 Page 209 Page 210 Page 211 Page 212 Page 213 Biographical sketch Page 214 Page 215 Page 216 
Full Text 
KINETIC PARAMETER ESTIMATION IN PREPACKAGED FOODS SUBJECTED TO DYNAMIC THERMAL TREATMENTS By BRUCE ARI WELT A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1996 UNIVERSiTY OF F.,' 'A L. :,. LD 1780 (1c) i /( i/ . SCI. NC. 'LI AY ACKNOWLEDGMENTS The author would like to express sincere gratitude to his major advisor, Dr. Arthur A. Teixeira and unofficial cochairman Dr. Murat 0. Balaban for their confidence and enthusiasm throughout this research program. Their friendship, guidance and support were critical ingredients to the successful completion of this work. The author also wishes to thank his supervisory committee Dr. Glen H. Smerage, Dr. David E. Hintenlang and Dr. Burrell J. Smittle for their ample guidance and suggestions related to the research and resulting technical papers which have been submitted for publication. Likewise, the author wishes to thank graduate coordinator, Dr. K.V. Chau for his friendship and guidance throughout the graduate program and for his assistance with the work covered in Chapter 3. A special thanks goes to the faculty and staff of the Agricultural and Biological Engineering Department and Department of Food Science and Human Nutrition for assistance with everything from operation of photocopy machines to pilot plant equipment. Special thanks also go to Dr. Daniel S. Sage for his assistance in designing the formal mathematical proof to support the work covered in Chapter 2. The author wishes to thank his parents, Martin and Ruth Welt, his brother and sister, Andrew and Jodi Welt, and his grandmothers Syd Welt and Pauline Schact, for their love and support through all of the best and worst of times. Finally, the author wishes to thank his fiancee, Janet E. Evans, for her support and devotion throughout what will be remembered as the beginning of our journey through life together. TABLE OF CONTENTS ACKN OW LED GM ENTS ........................................................................................... ii ABSTRACT ........................................................................................... ........ ........... vi 1. INTRODU CTION ................................................................................................. 1 2. KINETIC PARAMETER ESTIMATION IN FOODS THAT HEAT UNIFORM LY ...................................................................................... 4 Introduction .. ....................................................................................................... T h e o ry ................................................................................................................... 7 EPM Example ...................................................................... ........................ 12 Development of the Method of Paired Equivalent Isothermal Exposures (PEIE).. 14 PEIE Exam ple #1 ........................................................................................... .... 16 PEIE Example #2 .............................................................................................. 18 Discussion ...................................................................................................... 19 3. HEAT TRANSFER SIMULATION FOR THERMAL PROCESS DESIGN ........40 Introduction .......................................................................................................40 M eth o d s .................................................................................................. .......... 4 3 O v e rv iew ..................4................................................................................... 4 3 N odes and Volume Elements ......................................................................... 44 Energy Balances ............................................................................................47 Capacitance Surface N odes (CSN )....................................................................... 48 N onCapacitance Surface N odes (N CSN )......................................................... 49 Results and Discussion......................................................................................... 51 Lethality Calculation from Discrete TimeTemperature Data................................ 54 4. KINETIC PARAMETER ESTIMATION IN FOODS THAT DO NOT HEAT UNIFORM LY .................................................................................... 71 Introduction ......................................................................................................... 71 T h e o ry ............................................. ........................................................... .... .... 7 6 M e th o d s ............ ................................................................................................. 8 4 N umerical H eat Transfer............................................................................. 84 Experim ental .................................................................................................. 84 Results and D iscussion................................................................................... 88 5. EXPLOITING TREATMENTS THAT SENSITIZE BACTERIAL SPORES TO HEAT IN THERMAL PROCESS DESIGN.......................... 105 Introduction ................................................................... ............................... 105 M eth o d s .......................................................................................................... 1 10 R results and D iscussion....................................................................................... 114 6. SUMMARY AND CONCLUSIONS............................................................... 131 APPENDIX A: INVALID ASSUMPTION OF EPM ............................................. 135 APPENDIX B: HEAT TRANSFER CODEFINITE CYLINDER........................ 139 APPENDIX C: PEIE CODE FOR CONDUCTION HEATING C A N N E D F O O D ....................................................................................... 168 APPENDIX D: MICROBIAL ENUMERATION 0.0 kGy................................. 190 APPENDIX E: MICROBIAL ENUMERATION 1.0 kGy.................................... 192 APPENDIX F: MICROBIAL ENUMERATION 3.0 kGy.................................... 194 APPENDIX G: PROCESS PAIR DATA 0.0 kGy................................................ 197 APPENDIX H: PROCESS PAIR DATA 1.0 kGy................................................ 199 APPENDIX I: PROCESS PAIR DATA 3.0 kGy ................................................ 203 R E F E R E N C E S .......................................................................................................2 06 B IO G R APH ICAL SK ETCH ................................................................................... 214 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy KINETIC PARAMETER ESTIMATION IN PREPACKAGED FOODS SUBJECTED TO DYNAMIC THERMAL TREATMENTS By Bruce Ari Welt August, 1996 Chairman: A. A. Teixeira Major Department: Agricultural and Biological Engineering Consumer demand for higher quality processed foods provides an impetus for designing processes that are less detrimental to product quality, while maintaining adequate levels of security from food borne illness. Thermal process design requires knowledge kinetic parameters of food constituents that affect product safety and quality. The purpose of this study was to develop a method for estimating product and process specific thermal kinetic parameters. The method, referred to as the Paired Equivalent Isothermal Exposures (PEIE) method, may be applied to products that heat nonuniformly under nonisothermal heating conditions. Such conditions are typical with thermally processed prepackaged foods. The research was performed in four phases: (1) development of the PEIE method for the case of foods and/or samples that heat uniformly, (2) assessment of approaches to heat transfer simulation for foods that do not heat uniformly, (3) development of the PEIE method for the case of foods that do not heat uniformly and (4) demonstration of one potential application of the PEIE method involving hybrid irradiation/thermal process design. The PEIE method was used to quantify the impact of irradiation pretreatment on the thermal sensitivity of Bacillus stearotherniophilhis (NCA 1518) spores in canned pea puree. Cans of inoculated puree were irradiated with xrays at 0, 1.0 and 3.0 kGy at an electron linear accelerator facility prior to thermal treatment in a vertical stillcook retort. Enhanced thermal sensitivity was apparent from changes in Arrhenius kinetic parameters that described the rate of inactivation of the spores as a function of temperature. Mean rate constants, k121.1, at 121. PC and activation energy values, Ea, for spores in pea puree pretreated with 0, 1.0 and 3.0 kGy were 0.26, 0.28 and 0.38 min' and 250, 200 and 200 kJ/mol, respectively. Thermal processes designed to provide a 5 log cycle reduction in spore count were shown to require Fo values of 44.5, 41.0 and 30.0 minutes, respectively. Retort process simulations indicated that irradiation offered the potential for significant thermal process reductions. Steam shut off times of 72, 65 and 54 minutes, corresponding to products treated with 0.0, 1.0 and 3.0 kGy, respectively, were estimated for subsequent thermal treatment at 121.1 C. CHAPTER 1 INTRODUCTION Thermal processing is one of the most widely used methods of food preservation today. The principal objective of thermal processing is microbial inactivation by heat. Thus, thermal process design requires consideration of microbial inactivation kinetics and heat transfer principles. Since foods are very sensitive to environmental factors, especially temperature, relatively poor product quality is often the consequence of thermal processing. Consumer demand for higher quality processed foods provides an impetus for designing processes that are less detrimental to product quality, while maintaining adequate levels of security from food borne illness. Improving thermal processing equipment, operational techniques, and/or process design methods can contribute to higher quality thermally processed foods. Significant improvements in equipment and operational techniques have been made throughout this century. Relatively recent developments in computer automation technology have also contributed. However, methods of designing thermal processes remain virtually unchanged. Although it is now generally known that microbial inactivation kinetic parameters are functions of various external factors, the thermal processing industry continues to treat such parameters as immutable (Perkins et al., 1975). Examples of 2 external factors that influence inactivation kinetics include (1) conditions during sporulation, (2) conditions of spore recovery, (3) physical treatment of spores, and (4) the particular chemistry of the suspending medium. Since microbial inactivation parameters are product and process specific, efforts should be made to understand how such conditions affect thermal processing requirements. Additionally, methods such as irradiation are available that enhance sensitivity of microorganisms to heat. Accurately accounting for sensitized organisms in thermal process design could allow for significant reductions in thermal process severity. Therefore, methods are required that allow determination of product and process specific thermal kinetic parameters. The purpose of this study was to develop a method for estimating thermal kinetic parameters on a product and process specific basis. This approach requires a departure from traditional procedures, because it relies on survivor experiments performed insitu. Since foods may or may not heat uniformly and processes are seldom isothermal, such methods need to account for such situations. The research was performed in four phases: (1) development of the method for the case of uniformly heating foods, (2) assessment of approaches for heat transfer simulation in foods that do not heat uniformly, (3) development of the method for the case of foods that do not heat uniformly, and (4) application of the method for designing thermal treatments when a nonuniformly heating product was pretreated with relatively mild doses of ionizing radiation. Each phase of the work is reported separately in chapters 2, 3, 4 and 5. Overall concepts are progressively built from chapter to chapter, however, each chapter was prepared as a complete conceptual unit. The purpose of this format was to facilitate submission of the chapters, as independent entities, to scientific journals for hopeful publication as refereed articles. The benefit of this format is that each chapter may be read independently. This comes at the cost of some redundancy when covering the work in its entirety, since each chapter contains independent background material. Chapter 6 provides conclusions and implications relevant to the entire work and all references were compiled into one section at the end. CHAPTER 2 KINETIC PARAMETER ESTIMATION IN FOODS THAT HEAT UNIFORMLY Introduction The Equivalent Point Method (EPM) was developed by Swartzel (1982) in order to facilitate design of continuous flow ultrahigh temperature (UHT) processes. Subsequently, it was suggested that the EPM offered a reliable means for determining kinetic parameters from dynamic thermal processing conditions (Swartzel, 1984). The EPM has been used in kinetic analysis and thermal assessment of continuous flow systems (Wescott et al., 1995; Berry, et al., 1990; Sadeghi and Swartzel, 1990; Swartzel, 1986; Swartzel and Jones, 1985). The EPM was recognized as, "a significant advance in the application of food technology to food production" (Giese, 1994, p. 94), which contributed to earning the Industrial Achievement Award from the Institute of Food Technologists in 1994. This recognition is a testament to the value of practical nonisothermal methods for kinetic parameter estimation to the food industry. Two approaches are commonly used to obtain thermal kinetic parameters that involve either isothermal or nonisothermal experiments (Lenz and Lund, 1980). The isothermal approach is used more often, due to conceptual simplicity. However, the approach is often criticized for unavoidable thermal lags, inconveniently small sample sizes and a limited temperature range from which to select sufficiently different isothermal conditions (Welt et al., 1993). In fact these problems confound each other since thermal lags become more prevalent as either sample size or selected isothermal temperature increases. The isothermal approach may be preferable when the goal is reaction order determination (Constantinou, 1994). Although reactions of primary interest in food processing often involve complex mechanisms and multiple pathways, it is well known that most of these reactions follow apparent firstorder kinetics with Arrhenius temperature dependence (Labuza, 1979). Since thermal process design often involves the assumption of firstorder kinetics with Arrhenius temperature dependence, there is room for expanded use of the nonisothermal approach in kinetic parameter estimation. Nonisothermal approaches for kinetic parameter estimation are attractive because they are specifically designed to deal with thermal transients. This not only eliminates much of the experimental tedium associated with the isothermal approach, but also allows for determination of kinetic parameters from realistic processing conditions. Since food chemistry is very sensitive to thermal exposure and food chemistry influences the kinetic behavior of food constituents (mircoorganisms), kinetic parameters obtained from realistic processing conditions may be more appropriate for use in food process design. The use of nonisothermal methods is limited, however, due to increased complexity of data analysis. Several nonisothermal methods for kinetic parameter estimation have been described (Nasri et al., 1993; Lenz and Lund, 1980; Hayakawa et al., 1969). As proposed by Swartzel (1984) the EPM represented an easytouse non isothermal method for kinetic parameter estimation. Consistent with the name, the method suggests that an equivalent isothermal process may be obtained for any dynamic thermal process. This concept, inandofitself, is not new to food engineers, as processes designed to provide a specified sterilizing value (number of log cycles that the concentration of the target organism is reduced) are designed to achieve a specific Fo value (equivalent time at 121.1 C for a given target microorganism) (Teixeira, 1992). The Fo value provides a means to equate a dynamic process to an equivalent isothermal process for a particular target organism. The isothermal process may be described as an idealized rectangular pulse, where the product is instantaneously and uniformly heated to the isothermal temperature, held at this temperature for a certain amount of time, and then instantaneously and uniformly cooled to a nonlethal (non reacting) temperature. If we say the Fo of a dynamic process is 5 minutes, we are saying that the lethality (extent of reaction) of the target organism (reactant) is equivalent to what would have been observed under isothermal conditions at 121.1 C (250F) for 5 minutes. Notice that two pieces of information are required in order to specify an Fo value, (1) a particular reactant represented by its kinetic parameters and (2) an arbitrarily specified time or temperature. Selection of the temperature 121.1 C (250F) is a matter of adopted convention, since this is the temperature at which many industrial sterilizers are operated. Mathematically, the only difference between the EPM and Fo value method is that the EPM does not involve an arbitrarily specified time or temperature. Thus, the method must provide a means to satisfy this additional degree of freedom. The EPM appears to be much more powerful, however, since it is claimed to be able to specify an isothermal process that is universally applicable to all reactants that follow Arrhenius kinetics. This universality was implied by the assumption that the equivalent point "is independent of the Ea [activation energy] value" (Swartzel, 1984, p. 804). Careful review of the development of the EPM in this work, however, revealed that this assumption is not valid, suggesting that the EPM should not be used for kinetic parameter estimation without further consideration. The goal of this work was, (1) to describe and mathematically demonstrate the basis for questioning the assumption related to the universality of an equivalent point, (2) to develop an alternative approach for kinetic parameter estimation from non isothermal conditions, and (3) to demonstrate the new method with sample calculations. Theory The functional form of a firstorder rate process is dC d _kC (21) dt where C is the concentration of a reactant at time, t, and k is the reaction rate constant. Solving Eq. 21 by integration gives C =te (22a) Co or equivalently, ln( )= k.(tto) (22b) Co The temperature dependence of the rate constant, k, may be described by the Arrhenius equation, k= k0 exp E.} (23) LRT j where kinetic parameters ko and Ea are the preexponential factor and activation energy, respectively, R is the ideal gas law constant and T is absolute temperature. Under isothermal conditions, the extent of a reaction may be determined directly by substituting Eq. 23 into Eq. 22: ln(C ) =ko exp{ E. (t to). (24) Co c Rex{ .(T t In order to determine the extent of a reaction under nonisothermal conditions Eq. 24 becomes ln() exp Ea dt. (25) In( Co =  f, lR. T(t)J At this point the EPM deviates from the Fo value method. In order to determine the Fo value for a dynamic process with respect to a target organism, the righthand side of Eq. 25 is equated to the righthand side of Eq. 24. The isothermal temperature, T, in Eq. 24 is arbitrarily set to 394.25 Kelvin (121. 1C), the known activation energy, Ea, of the target organism is substituted on both sides, and the resulting equation is solved for time, t, from Eq. 24. The equivalent time at 121.1 C is referred to as the Fo value of the dynamic process, with respect to the particular target organism (Ea value). In a similar manner, the EPM equates Eq. 25 to Eq. 24, but introduces the extra step of normalizing the resulting equation by ko. The resulting equation takes the form In( C E G } exp E dt = t,exp Ea (26) ko to p R. T(t)J I RT J where G, as defined by Swartzel (1982), is the value obtained by evaluating the integral from the actual thermal history. As with the Fo value method, the time, t., at isothermal temperature, Te, comprises the equivalent isothermal exposure which would provide the same extent of reaction (lethality), as the dynamic thermal process. Swartzel (1982) observed that straightline equations of the form ln(t) = ln(G)+ Ea (27) resulted when these G values were equated to the isothermal part of Eq. 26 (right hand side). Using several of these equations from arbitrarily selected Ea values, Swartzel (1982) plotted the lines on an ln(te) versus 1 /Te plot and noticed that the lines had a tendency to intersect in a particular region of the plot. The effect was sufficiently convincing that Swartzel (1982) claimed that all such lines intersect at one point. Swartzel (1984, p. 804) noted situations where three or more lines did not intersect at the same point, but attributed this to experimental uncertainty and loss of arithmetic precision during calculations: Any accumulated error in the process thermocouplee calibration, data collection, determination of exponential integral, etc.) would be indicated with less than perfect intersections with three or more log(t) versus T lines. The intersection points of these lines may be evaluated by regression analysis to determine the most error free statistical intersection point for all lines. Although the EPM specifies the need for more than one activation energy to define the equivalent process, algebra dictates that not more than two are actually necessary. Since an isothermal process is functionally represented as G te exp{ E (28) two sets ofEa G values provide two equations in two unknowns, Te and t,, which may be solved as follows: ln(G,)= ln(t) Eal R.T, (29a) ln(G2)=ln(te)+ E.2 R.T, T (Ea2 Eal) (29b) R.ln(Gj G2) Equivalent time, te, is found by backsubstitution of the equivalent temperature, Te, into either equation of Eq 9a. In order for the equivalent point to be independent of Ea, all additional equations would have to provide the same solution. To see why this is not possible, consider Eq. 26, Eq. 27 and Eq. 29. Equation 6 shows that for a given actual thermal history, T(t), G is a monotonic function of Ea; each different value of Ea yields a different and distinct value of G. With Ea and G defined, substitution into Eq. 27 yields one unique equation with two unknowns, te and Te. Equation 9 shows that te and Te may be specified from two different Ea G versions of Eq. 27 by solving them simultaneously. Since the two equations in Eq. 29a are independent and unique, their solution, or point of intersection, must also be unique. It is not possible for Eq. 29 to yield the same solution for any other pairs of independent equations derived from different Ea values. A rigorous proof of this statement was provided in Appendix A. Thus, it is not possible to specify a unique equivalent isothermal process from a dynamic thermal exposure that is universally applicable to all reactants. An isothermal process derived from two Ea values may be interpreted as the equivalent exposure that would result in the same extents of reaction (lethalities) for reactants (microorganisms) characterized by the two respective Ea values. An isothermal process specified in this manner may be referred to as an equivalent isothermal exposure (EIE). The definition of an EIE differs from the Fo value. The Fo value provides an equivalent time at a given temperature, such that the same extent of reaction (lethality or sterilizing value) would be observed for reactants (microorganisms) characterized by the one Ea value used to determine the Fo value. An EIE specifies an isothermal exposure (time and temperature) that would provide equivalent extents of reaction (lethality) for reactants characterized by the two E, used to define the EIE. EPM Example The following example illustrates the dependence of the equivalent point on activation energy. As proposed, the EPM involves the following steps: 1. Record the dynamic temperature history. 2. Arbitrarily select several Ea values and determine corresponding G values by evaluating the integral in Eq. 26, using the actual thermal history. 3. Equate G values to the isothermal form of the equation (right most side of Eq. 26), using respective Ea values. 4. Generate straight line curves on a ln(te) versus l/To plot, for each Ea G equation. 5. Obtain the equivalent isothermal process from a unique intersection point. Step 5 is the point of contention, as it shall be demonstrated that the lines generated in step 4 do not intersect at one point. For the purpose of this example, we assumed that a food item was uniformly and nonisothermally heated from 0C to 140C over a total process time of 200 minutes, according to the following linear equation (see Fig. 21, Process a): T(Kelvin) = 273 + 0.7t (210) where T is in Kelvin and t is in minutes. Activation energy, Ea, values of 30, 100 and 400 kJ/mol were arbitrarily selected for equivalent point calculations. The following integral was solved in order to calculate G values for each E, value: 200 E ( G = expREa } dt (211) 0 R R(273 + 0.7t)J The integrals were numerically evaluated to a tolerance of 1 x 106, using the programmable spreadsheet program, Excel (Microsoft Corp., Redmond, WA). An integration algorithm of Simpson's rule (Press et al., 1992) was used to program the spreadsheet using Visual Basic (Microsoft Corp., Redmond, WA). Calculated G values are shown in Table 21. By plotting the three respective forms of Eq. 27, three points of intersection, representing three different isothermal processes, were obtained for each pair of Ea values (i.e. for E, value pairs 12, 23 and 13). The equivalent isothermal processes obtained were 392.4 Kelvin for 87.8 minutes, 404.0 Kelvin for 67.5 minutes, and 406.8 Kelvin for 29.7 minutes, relating to activation energy pairs of 30 and 100, 30 and 400, and 100 and 400 kJ/mol, respectively. The spatial relationship of these intersection points is illustrated in Fig. 22. This example shows that one unique intersection point should not be expected to exist. Thermocouple errors were not a factor here, as this was a mathematical exercise. Additionally, the same points of intersection were obtained when integration was performed to a tolerance of 0.01, indicating that precision of numerical integration and rounding errors did not alter the result significantly. Therefore, this example supports the conclusion that a unique equivalent isothermal process, independent of activation energy, should not be expected to exist. This suggests that kinetic parameter estimates obtained with the EPM may not be accurate. Use of inaccurate kinetic parameter estimates could result in insufficient thermal process specifications. Development of the Method of Paired Equivalent Isothermal Exposures (PEIE) A method, referred to as Paired Equivalent Isothermal Exposures (PEIE), was developed for estimating kinetic parameters from dynamic thermal processes using the true nature of an equivalent isothermal exposure (EIE). For a given dynamic thermal exposure, the point of intersection of any two E, G lines on a ln(te) versus 1 /Te plot, specifies an EIE (time, to, at temperature, Te) with respect to reactants characterized by the two particular Ea values used to specify the EIE. The method of PEIE exploits this definition to allow estimation of thermal kinetic parameters of a reactant that has been subjected to at least two dynamic thermal exposures. Since two Ea, values are required to define an EIE from one dynamic thermal exposure, when one of the two Ea values is the apparent Ea value of the reaction, the corresponding isothermal rate constant, kT, when plotted on an Arrhenius plot, must fall somewhere along the desired, but unknown Arrhenius curve. The other Ea, value simply defines the location of this point along the Arrhenius curve. Thus, if the reactionapparent Ea value is used to define EIE's from any pair of dynamic thermal exposures, the desired Arrhenius curve may be specified. The crux of the method of PEIE, therefore, is to find the reactionapparent Ea value for all pairs of data related to dynamic thermal exposures. The method of PEIE's is iterative in nature, whereby an initial Ea guess is increasingly improved such that corresponding EIE's result in points that converge to the appropriate Arrhenius curve. Two examples of the PEIE method are provided to demonstrate (1) that the PEIE method is capable of finding the exact parameters used to generate simulated experiments and (2) that the method is capable of finding the best possible parameters using uncertain, but realistic extentofreaction data. The following steps outline the method of PEIE: 1. Record temperature histories and initial and final concentrations of a reactant from at least two dynamic thermal processes that result in different extents of reaction. 2. Arbitrarily select one Ea value, Eal, and establish a scheme for specifying the other Ea value, Ea2. For each iteration, the value of Ea2 may be a fixed value or some multiple ofEal. Theoretically, the choice of Ea2 does not matter, as is the case when extentofreaction data are exactly known. Implications of various Ea2 specification schemes for the typical case of uncertain extentofreaction data are discussed below. 3. Determine respective EIE's for each pair of dynamic thermal experiments using the Ea values specified in step 2. 4. Calculate isothermal rate constants, kT, with Eq. 22b using the EIE specifications (from step 3) and respective extentofreaction data (from step 1). 16 5. For each pair of kT values determined in step 4, simultaneously solve Eq. 23, to find Ea values of curves that pass through respective pairs of points. Graphically, each kT represents a point on an Arrhenius plot. Values of Ea are determined from the slopes of curves that pass through each pair of points (e.g. if three dynamic thermal experiments were performed, calculate Ea estimates for each of the pairs 12, 13, and 23). 6. Return to step 3 and replace the old E al values with respective Ea estimates obtained in step 5. Repeat the procedure until Ea 1 stops changing. Create an Arrhenius plot in the typical fashion with kT values from all final EIE specifications. Obtain kinetic parameters in the typical fashion by performing a regression analysis of all kT points on the Arrhenius plot. PEIE Example #1 The following example was designed to demonstrate the ability of the PEIE method to determine kinetic parameters that were used to generate extentofreaction data from dynamic thermal experiments. Consider experiments involving thiamin in pea puree using the three dynamic thermal exposures depicted in Fig. 21. Expected extents of thiamin degradation from each exposure may be calculated by substituting representative kinetic parameters into Eq. 25 (Ea = 113 kJ/mol and k, = L.OE+13 min ', Welt et al., 1993). Extents of thiamin degradation were calculated and shown in the column labeled 'Exact (C/Co)' in Table 22. The PEIE method was applied to these data as follows: Step 1: See Table 22. Step 2: For the first iteration, Eal was chosen to be 30 kJ/mol. A scheme which specified Ea2 as 2 x Eal was selected (see Table 23). Step 3: Equivalent isothermal exposures for each process pair (ab, ac and bc) were determined as previously described (see columns labeled "Te" and "t," in Table 23). Step 4: Isothermal rate constants, kT, were calculated using EIE and final concentration data. These calculations are illustrated below for the first iteration, and summarized in Table 23: k3a86435 =ln(0.421)/ 101.194 8.548 x 103 min' 377.662 3  k77662 in(0.600) /135.230 = 3.773 x 103 min1 (212) k368882 = ln(0.773)/160.974 = 1.600 x 103 min1 Step 5. Activation energy, Ea, estimates were calculated for each pair of kT values by substitution into respective Arrhenius equations (Eq. 23). Solving these simultaneously yields Eq. 213 R.ln(k k2) Ea=  ('213) T1 T2 VI li T .T2 ) where the subscripts, 1 and 2, refer to different dynamic thermal processes. Step 6: The E, estimates were used to replace respective Eal values. Ea2 values were appropriately adjusted and the second iteration followed by repeating the procedure from step 3 (Table 24). Details from the third and final iteration are shown in Table 25. Step 7: Regression analysis was performed on ln(kT) versus 1/Te data from the final iteration (Table 25). Resulting kinetic parameters were Ea 113,000 J/mol and k, = 1.00 x 1013 min1. Table 26 shows the progress of each iteration (see 'Iteration #1 3,' Table 2 6). Although all calculations were performed with full precision (15 digit) arithmetic, not all digits are shown in the tables. Simpson's rule was used to calculate the G value integral in Eq. 26 with a tolerance of 1 x 10.6 (Press et al., 1992). Improvements to the parameters were indicated by reductions to the sum of squared error (SSQ error) of predicted data (Table 26). PEIE Example #2 The same dynamic thermal processes (Fig. 21) and kinetic parameters (Ea = 113 kJ/mol, ko = 1.OE+13 min1) were used as a basis for this example. However, experimental uncertainty was introduced in the calculated final concentration data. For each data point, random values were selected from an evenly distributed range specified by the exact final concentration (column labeled 'Exact,' Table 22), plus or minus 5 percent. Equations of the thermal histories and obsen'ed extentofreaction data were provided in Table 22 (column labeled 'Observed,' Table 22). For the first iteration, a value of 300 kJ/mol was arbitrarily chosen for Eal. As with the previous example, Ea2 was set to 2 x E,1 for each iteration. Details of all iterations are provided in Tables 7, 8 and 9. Progress of iterations was provided in Table 210. The best parameters given the observed data were Ea = 120,289 J/mol and ko = 9.270E+13min"1. Discussion The fact that the method of Paired Equivalent Isothermal Exposures (PEIE) was capable of precisely determining kinetic parameters used to generate the data in example #1 provides confidence that the method works as proposed. An absolute performance comparison to the Equivalent Point Method (EPM) (Swartzel; 1982, 1984) is not possible, because contrary to the fundamental assumption of the EPM, the result depends on the Ea values used to determine an equivalent point. Although this work presented reasons to question the fundamental assumption regarding the universality of equivalent points, the mechanics of applying the EPM to kinetic parameter estimation, while not identical, are similar to those of the first iteration of the PEIE method. With the EPM, if a researcher has a preconception of the true Ea value and/or is fortunate enough to select a value that is close to the true value, then the EPM may provide parameters that are close to what might be considered correct parameters. This provides an explanation for the apparent success of the EPM in the literature. A notion of the differences in results obtained from the EPM and PEIE method can be seen by comparing results of the first and last iterations of the PEIE method. For example #1, Arrhenius curves that provided kinetic parameters from the first and last PEIE iterations are shown in Fig. 23. The magnitude of these differences may be assessed by comparing extentofreaction predictions from these sets of parameters to the exact experimental data in Table 26. The SSQ errors obtained using parameters from the first and last iterations were 6.79E02 and 1.93E14, respectively, indicating that the PEIE method provided closer agreement with exact experimental data. Figure 24 and Table 210 show that the same was true for example #2. The purpose of example #2 was to demonstrate, (1) how experimental uncertainty introduces additional points on the equivalent Arrhenius plot (Fig. 24), and (2) how these points introduce an additional consideration regarding selection of a scheme for specifying Ea2 for each iteration. For the first iteration, the same values of Eal and Ea2 were used for all pairs of experiments (Table 27). The resulting EIE for process a when paired with process b was the same as the EIE for a when paired with c. Thus, only three points appeared on the Arrhenius plot from the first iteration (Fig. 24). Since there was uncertainty in the extentofreaction data, Eal and E,2 values differed for each pair of experiments in subsequent iterations. Therefore, the EIE for process a when paired with process b, was different from the EIE for a when paired with c. Since three dynamic thermal processes were considered in the example, six points appeared on the final Arrhenius plot ('Iteration #3,' Fig. 24). Ideally, the value of Ea2 should not influence the result. Figure 25 is an Arrhenius plot for example #1, where extentofreaction data were known exactly. The figure shows points from final iterations when Ea2 was specified according to the following schemes: (1) 10 kJ/mol, (2) E,1 / 2, (3) El x 2 (4) EI x 4 and (5) 1,000 kJ/mol. As expected, Fig. 25 shows that the same Arrhenius curve was defined regardless of the scheme used to specify Ea2. However, in the face of experimental uncertainty, additional points appear on the Arrhenius plot and their spatial relationships may artificially bias the final regression analysis. Figure 26 is an Arrhenius plot for example #2 (uncertainty in extentofreaction data), showing points from final iterations when Ea2 was specified according to the schemes just described. When Ea2 was relatively small, EIE's involved low temperatures and long times, resulting in relatively low isothermal rate constants (lower right, Fig. 26). Additionally, small Ea2 values resulted in points that tended to have greater horizontal scatter (EIE's with larger temperature differences) than those from high Ea2 values. Conversely, high Ea2 values resulted in points that tended to be more vertically oriented (upper left, Fig. 26). Both extremes may artificially influence regression calculations, resulting in poorer quality parameter estimates. As expected, all points in Fig. 26 tended to represent one particular Arrhenius curve, however, results from separate regressions show some variation. Table 211 provides separate regression results and resulting SSQ error values for all cases considered in Fig. 26. The table shows that for the present example the best parameters, as indicated by the smallest SSQ error, were obtained when Ea2 was specified as Ea I x 4. Inspection of Fig. 26 shows that this scheme resulted in points that were fortuitously oriented with respect to the underlying kinetic behavior. Although the Eal x 4 scheme provided the lowest SSQ error value of the Ea2 schemes considered here, the Ea I x 2 scheme was chosen for demonstration in order to emphasize the fact that in most cases one reasonable Ea2 scheme is sufficient for the PEIE method to provide adequate kinetic parameter estimates. Comparison of SSQ values indicates marginal improvement in predictive ability of the kinetic parameters. Therefore, the Ea I x 2 scheme was considered to be a useful and easily recallable guideline for kinetic parameter estimation of food constituents with the PEIE method. The method of PEIE requires data from at least two dynamic thermal processes to generate kinetic parameter estimates. However, when only two experiments are used, no indication of uncertainty in the parameter estimates is evident. As with the traditional isothermal kinetics approach, a minimum of three data points (three non isothermal experiments) is recommended, as linear regression of the Arrhenius plot helps to mitigate the effect of random uncertainty. Additionally, statistical analysis of the regression provides insight into the reliability of the kinetic parameter estimates. As with any method of kinetic parameter estimation, the ultimate success of the PEIE method in generating reliable kinetic parameters depends on the quality of the original experimental data. Table 21. Ea and corresponding G values calculated from Eq. 28. Index Ea Value G Value (J/mol) 1 30,000 8.914E03 2 100,000 4.277E12 3 400,000 1.281E50 Table 22. Hypothetical dynamic thermal experiments used to demonstrate calculations of the Paired Equivalent Isothermal Exposures (PEIE) method. Process Temperature Process Exact Observed Function Time (simulated) (t in min) (min) (C/Co) (C/Co) a T(K) = 273 + 0.7t 200 0.421 0.412 b T(K) = 273 + 0.5t 260 0.600 0.588 c T(K) = 273 + 0.4t 300 0.773 0.793 Table 23. First iteration of the PEIE method for example #1. Progress from each iteration is shown in Table 26. Process Ea G Te te k Ea Pair Guesses Estimate (J/mol) (K) (min) (1/mmin) (J/mol) a 30,000 8.914E03 a 60,000 7.853E07 386.435 101.194 8.548E03 b 30,000 9.590E03 b 60,000 5.137E07 377.662 135.230 3.773E03 113,101 a 30,000 8.914E03 a 60,000 7.853E07 386.435 101.194 8.548E03 c 30,000 9.094E03 c 60,000 5.137E07 368.882 160.974 1.600E03 113,162 b 30,000 9.590E03 b 60,000 6.801E07 377.662 135.230 3.773E03 c 30,000 9.094E03 c 60,000 5.137E07 368.882 160.974 1.600E03 113,220 Table 24. Second iteration of PEIE for example #1. Progress from each iteration is shown in Table 26. Process Ea G Te te k Ea Pair Guesses Estimate (J/mol) (K) (min) (1/min) (J/mol) a 113,101 8.393E14 a 226,201 2.141E28 404.824 32.897 2.629E02 b 113,101 4.947E14 b 226,201 1.145E29 395.204 43.938 1.161E02 113,000 a 113,162 8.240E14 a 226,324 2.065E28 404.828 32.880 2.631E02 c 113,162 2.447E14 c 226,324 1.145E29 385.579 52.306 4.923E03 113,000 b 113,220 4.769E14 b 226,441 5.181E29 395.212 43.895 1.162E02 c 113,220 2.403E14 c 226,441 1.104E29 385.583 52.281 4.925E03 113,000 Table 25. Third iteration of PEIE for example #1. Progress from each iteration is shown in Table 26. Process Ea G Te te k Ea Pair Guesses Estimate (J/mol) (K) (mnin) (1/mmin) (J/mol) a 113,000 8.650E14 a 226,000 2.272E28 404.817 32.924 2.627E02 b 113,000 5.103E14 b 226,000 1.266E29 395.197 43.974 1.160E02 113,000 a 113,000 8.650E14 a 226,000 2.273E28 404.817 32.924 2.627E02 c 113,000 2.575E14 c 226,000 1.266E29 385.569 52.375 4.916E03 113,000 b 113,000 5.103E14 b 226,000 5.921E29 395.197 43.974 1.160E02 c 113,000 2.575E14 c 226,000 1.266E29 385.569 52.375 4.916E03 113,000 Table 26. Progress of PEIE iterations for example #1. Predicted from iteration # 1 2 3 ko(1/min) 1.690E+13 9.973E+12 1.000E+13 Ea (J/mol) 113,162 112,990 113,000 Process Exact (C/Co) (C/Co) (C/Co) (C/Co) a 0.421 0.24846 0.42094 0.42107 b 0.600 0.44023 0.60022 0.60035 c 0.773 0.66131 0.77289 0.77299 SSQ Error 6.790E02 4.327E08 1.932E14 Table 27. First iteration of PEIE for example #2 (Ea2 = 2 x Eal 1). Progress from each iteration is shown in Table 210. Process Ea G Te te k Ea Pair Guesses Estimate (J/mol) (K) (min) (1/min) (J/mol) a 300,000 7.541E38 a 600,000 4.354E76 409.801 13.062 6.797E02 b 300,000 1.151E38 b 600,000 9.489E80 399.952 17.426 3.044E02 111,153 a 300,000 7.541E38 a 600,000 4.354E76 409.801 13.062 6.797E02 c 300,000 1.403E39 c 600,000 9.489E80 390.100 20.731 1.121E02 121,607 b 300,000 1.151E38 b 600,000 7.601E78 399.952 17.426 3.044E02 c 300,000 1.403E39 c 600,000 9.489E80 390.100 20.731 1.121E02 131,555 Table 28. Second iteration of PEIE for example #2 (Ea2 = 2 x Eal). Progress from each iteration is shown in Table 210. Process Ea G Te te k Ea Pair Guesses Estimate (J/mol) (K) (min) (1/min) (J/mol) a 111,153 1.504E13 a 222,306 6.771E28 404.689 33.426 2.656E02 b 111,153 8.994E14 b 222,306 6.073E32 395.075 44.646 1.188E02 111,242 a 121,607 6.580E15 a 243,214 1.407E30 405.365 30.768 2.886E02 c 121,607 1.724E15 c 243,214 6.073E32 386.068 48.934 4.748E03 121,681 b 131,555 1.738E16 b 263,110 7.910E34 396.244 38.195 1.389E02 c 131,555 7.617E17 c 263,110 1.276E34 386.567 45.481 5.109E03 131,617 Table 29. Third iteration of PEIE for example #2 (E,2 = 2 x Eal). Progress from each iteration is shown in Table 210. Process Ea G Te te k Ea Pair Guesses Estimate (J/mol) (K) (min) (1/min) (J/mol) a 111,242 1.465E13 a 222,484 6.423E28 404.695 33.402 2.658E02 b 111,242 8.751E14 b 222,484 5.800E32 395.081 44.613 1.189E02 111,242 a 121,681 6.436E15 a 243,362 1.347E30 405.369 30.751 2.887E02 c 121,681 1.684E15 c 243,362 5.800E32 386.072 48.906 4.751E03 121,681 b 131,617 1.705E16 b 263,235 7.617E34 396.247 38.178 1.389E02 c 131,617 7.470E17 c 263,235 1.227E34 386.570 45.461 5.111E03 131,617 Table 210. Progress of PEIE iterations for example #2. Predicted from iteration # 1 2 3 ko (1/mmin) 2.262E+14 9.267E+13 9.270E+13 Ea (J/mol) 121,691 120,289 120,289 Process Observed (C/Co) (C/Co) (C/Co) (C/Co) a 0.412 0.23422 0.40476 0.40469 b 0.588 0.44752 0.60282 0.60277 c 0.793 0.68400 0.78533 0.78530 SSQ Error 6.308E02 3.094E04 3.092E04 Table 211. Results of example #2 with various Ea2 specification schemes. Results per Ea2 scheme Ea2 = 10.000 Eal/2 Eal x4 1,000,000 ko(1/min) 9.189E+14 2.452E+14 6.406E+13 4.696E+13 Ea(J/mol) 127,578 123,447 119,063 118,025 Process Observed (C/Co) (C/Co) (C/Co) (C/Co) (C/Co) a 0.412 0.36235 0.39421 0.40565 0.40567 b 0.588 0.58340 0.60101 0.60085 0.59856 c 0.793 0.78395 0.78874 0.78232 0.77936 SSQ Error 2.520E03 4.763E04 2.982E04 3.161E04 0 50 100 150 200 250 300 Tine (nin) Dynamic thermal processes used to demonstrate calculations for the Equivalent Point Method (EPM) and the method of Paired Equivalent Isothermal Exposures (PEIE). Figure 21. 385 390 395 400 405 410 415 Temperature (K) Spatial relationship of equivalent isothermal processes, as calculated by the Equivalent Point Method. The figure demonstrates that one unique intersection point, independent of activation energy, should not be expected to exist. Figure 22. Iteration #3: ko = 1.OOOE+13 (1/min Ea = 113,000 (J/mol) Iteration #1: ko = 1.690E+13 (1/min) Ea = 113,162 (J/mol) Arrhenius plots of first and last iterations of the PEIE method for example #1. Resulting kinetic parameters matched those used to generate the experiment. 0.100 0.001 2 Ij .45 2.50 2.55 2.60 2.65 2.70 Inverse Absolute Temperature (x 1000) Figure 23. 0.100 S0.010 0.001 2.42 2.44 2.46 2.48 2.50 2.52 2.54 2.56 2.58 2.60 Inverse Absolute Temperature (x 1000) Arrhenius plots of first and last iterations of the PEIE method for example #2. Experimental uncertainty was introduced to extentof reaction (lethality) data. Iteration #1: ko = 2.262E+13 (1/min) ' Ea = 121,691 (J/mol) Iteration #3: ko = 9.270E+13 (1/min) Ea = 120,289 (J/mol) Figure 24. A Ea2 Scheme o 10,000 ,AEal/2 o Eal x 2 xEal x4 o 1,000,000 ko = 1.00OE+13 (1/min) Ea = 113,000 (J/mol)  C 2.45 2.50 2.55 2.60 2.65 2 Inverse Absolute Temperature (x 1000) Arrhenius plot showing resulting points and regression lines from various Ea2 specification schemes for example #1. 2.50  3.00 3.50 4.00 4.50  5.00  5.50  6.00  6.50  7.00  2.40 Figure 25. ^'n 2.50  3.00 3.50 4.00 4.50 5.00 5.50 6.00 6.50 7.00  2.40 Figure 26. 2.45 2.50 2.55 2.60 2.65 2.70 Inverse Absolute Temperature (x 1000) Arrhenius plot showing resulting points and regression lines from various Ea2 specification schemes for example #2. CHAPTER 3 HEAT TRANSFER SIMULATION FOR THERMAL PROCESS DESIGN Introduction Use of numerical methods for simulating heat transfer for thermal process calculations continues to increase in sympathy with the prevalence and power of modem computers. Numerical methods are favored because, in many cases, analytical solutions do not exist, Numerical methods sacrifice the freedom of ascertaining information at any location and time, by providing data at only specific points in space and time. Efforts to mitigate this limitation spawned the development of numerous algorithms, each offering particular benefits and drawbacks. Thermal process analysis embodies both heat transfer simulation and reaction kinetics. Development of a reliable thermal process analysis tool requires proper selection of a numerical method, and a sound approach for incorporating discrete timetemperature data into lethality (extent ofreaction) predictions. Desirable aspects of any numerical method include accuracy, minimal computation time and ease of programming. Programming ease is typically sacrificed for gains in accuracy and/or overall computation time. Methods of sufficient accuracy for thermal process calculations include both explicit and implicit algorithms. Briefly, explicit algorithms are so named because all information required to generate temperature estimates at each time step is explicitly specified by temperature estimates of the previous step. On the other hand, implicit algorithms require information that must be determined concurrently with each temperature estimate. As a result, implicit algorithms are more difficult to program and require greater computational time on a per time step basis (Thibault, 1985). However, implicit methods are inherently stable, allowing greater leaps in time, resulting in fewer steps and overall time savings. Explicit methods are often criticized because the time step is limited by the laws of thermodynamics. The time step must be small enough to ensure that heat is always considered to flow from hotter to colder regions. However, many texts describe explicit methods in a manner that unnecessarily exaggerates the effect of this limitation (e.g. Lienhard, 1987; Chapman, 1984). This is caused by a seemingly innocuous habit of associating thermal mass with the twodimensional surface of the object. Both explicit and implicit heat transfer approaches involve estimating energy flows between previously defined points, which are often referred to as nodes. Since nodes are one dimensional points, volume elements constructed around nodes, provide the corresponding thermal capacity required to construct energy balance equations. The common practice of defining equally spaced nodes results in volume elements whose nodes are on the surface. In the terminology of Clausing (1981), nodes on the surface may be considered, capacitance surface nodes (CSN). As demonstrated below, this approach suffers from the fact that the description of the problem is physically inconsistent, because instantaneous temperature changes of the twodimensional surface are not consistent with actual temperature changes of associated three dimensional volume elements. This results in a situation where the maximum allowed time step for a stable solution decreases as the surface heat transfer coefficient increases, unnecessarily increasing computational time. Analyzing such problems using noncapacitance surface nodes (NCSN) results in a physically consistent depiction of the problem and a practical time step restriction. In this case, surface nodes are not considered to have thermal mass. Thus, in a convective environment, heat transfer through the surface is analyzed as a combination of resistances to heat transfer involving both convection and conduction. When treated in this manner, the time step is bounded as the convective heat transfer coefficient increases. Although very little evidence of the NCSN approach exists in the literature, it is ideally suited for many thermal processing situations. Researchers who have used the NCSN approach include Bellagha and Chau (1985), Chau and Gaffney (1990) and Silva et al. (1992). Estimation of the impact of a process on a thermally labile constituent requires integration of the temperature sensitive function (usually the Arrhenius equation) over the temperature history, including space and time. Errors are expected with all numerical methods, because data are available at only certain points in space and time. Since the process of spatial discretization is similar for all finite difference methods, no particular numerical method offers significant advantages over others. Although implicit methods do not restrict time step size, errors related to estimation of extent of reactions during each time interval, grow as the time step increases. Thus, when the goal is extentofreaction prediction, an explicit algorithm with a practical time step would be preferable to an implicit method with a similar time step, because the explicit algorithm is both easier to program and executes faster. The objectives of this work were (1) to compare methods of formulation of explicit methods using both the capacitance and noncapacitance surface node approaches, (2) to demonstrate that both methods offer comparable accuracy, but with greatly different time step restrictions, and (3) to present guidelines for using discrete timetemperature data for approximating extents of reactions (lethality) in systems whose thermal histories are estimated with numerical heat transfer methods. Methods Overview Two approaches are commonly used to develop finite difference representations, namely, (1) truncated Taylor expansions and (2) energy balance. The first approach is often described in mathematics and heat transfer texts (e.g. Stoer and Bulirsch, 1993; Ozisik, 1985; Haberman, 1983), and was employed in the earliest applications for thermal processing (Teixeira et al., 1969). The text by Ozisik (1985) provides an excellent presentation on the development of finite difference equations using both difference quotients and energy balances; however, both reflect the CSN case. The second approach offers greater flexibility, because energy balance equations are developed on a volume element basis, allowing location specific variations to be considered. Clausing (1981, p. 170) strongly recommended the second approach by stating, . if one blindly employs Taylor's expansion or various difference quotients to obtain required difference representations, a physically consistent set [of equations] will probably not result. . Although it may be possible to obtain a physically consistent set employing Taylor's expansions, it is not always clear how this is done nor how one would be assured of this consistency without making a comparison with the physically derived relationships. Application of finite difference methods are often described by one or two dimensional problems in Cartesian coordinates. Chau and Gaffhey (1990) demonstrated the spherical case, and Silva et al. (1992) considered infinite geometries. Many practical problems in thermal processing involve finite cylinders. For these reasons application of explicit finite difference methods using the energy balance approach on a finite cylinder was considered. Specifically, CSN and NCSN approaches were developed for an arbitrary conduction heating finite cylinder, subjected to convective boundary conditions. For both cases, the following steps were performed: (1) define node locations and volume element boundaries, (2) perform energy balances on each type of node, and (3) develop equations for determining the maximum allowed time step for stability. Nodes and Volume Elements Spatial descritization is a critical step in finite difference heattransfer modeling. A node is considered to reside in the central region of a volume element. Each volume element is defined by imaginary boundaries that represent separate control volumes, on which to perform energy balances. Although nodes may be placed arbitrarily, it is convenient to use a standard regimen. The most straightforward method is to separate all nodes by equal distances along the dimensions of the object. This is the typical basis of the CSN approach. Thus for Nr nodes along the radius, and Ny nodes along the halfheight, distances, dr and dy, separate nodes along the radius and halfheight, respectively. These distances are calculated as follows (see Fig. 31): dr = Radius / (Nr 1) (31a) dy = HalfHeight / (Ny 1) (31 b) The situation is slightly different for the NCSN case (Fig. 32). For Nr nodes along the radius, all interior nodes (node 1 through node Nr 1) are separated by the same distance, dr, while the distance between the last interior node, Nr 1, and the surface node, Nr, is dr/2. This provides 2Nr 3 subintervals of length dr/2 along the radius. Thus, distances, dr and dy, between interior nodes along the radius and half height are given by (see Fig. 32 and Fig. 34) dr 2 Radius / (2 Nr 3) (32a) dy = 2 HalfHleight / (2 Ny 3) (32b) Visual Basic code for discretizing a finite cylinder using the NCSN is shown in Fig. 33. Note that several lines of code (Fig. 33) have additional code in square brackets (e.g. [dr = Radius / (Nr 1)]). Substitution of the bracketed code provides the CSN case. Knowledge of surface temperatures are not required to estimate internal temperatures with the NCSN approach, nor are surface temperatures required for estimation of extent of reactions (lethality). In practice, however, it is worthwhile to develop equations for surface nodes so that they are accounted for in computer code, in case the information is subsequently desired. It is important to remember that for the NCSN approach, surface nodes are not necessary and are treated differently from interior nodes. Specification of volume element boundary locations requires consideration of the geometry of the object. Heat flowing vertically in an upright finite cylinder passes through crosssections of equal area, from the bottom to the top of the cylinder. Heat flowing along the radial direction, however, passes through nonconstant areas. Thus, tops and bottoms of volume elements are placed at the mean distance between nodes. Sides of volume elements are placed at the logmean radius between nodes (Fig. 33): rim = (ro ri) / ln(r o/ri) (33) where ro, ri, and rim are outer radius, inner radius and log mean radius, respectively. Since volume elements containing the vertical centerline have a zero inner radius, ri, the boundary may be specified at the mean distance between the first and second radial nodes. It should be noted that for small dr, very little error is introduced by simply using the mean distance between nodes (Lienhard, 1987). The area of the side of a volume element is equal to the height of the element times the circumference of the circle defined by the radial location of the boundary. The areas of the top and bottom of each node are equivalent. The area of the tops/bottoms are equal to the difference of circular areas defined by the radial locations of the outer and inner side boundaries (Fig. 33). Volumes of elements for both approaches are calculated in a similar manner. Since each node for the CSN approach resides within a volume element, Nr x Ny volume elements are created. In the NCSN case, surface nodes do not reside in volume elements, thus Nr 1 x Ny 1 volume elements are created (Fig. 33). The volumes of each element may be calculated by multiplying the area of the top boundary by the height of the volume element. Energy Balances Once nodes and volume elements are specified, energy balances may be performed. The first step is to group volume elements according to common geometries and heattransfer situations. Nine types of volume elements were identified for both the CSN and NCSN nodal arrangements for a finite cylinder (labeled A I, Fig. 31 and Fig. 32 and Table 31 and Table 32). In writing energy balance equations, it is customary to consider heat flowing toward each node through the perpendicular area of the volume element boundary. For both the CSN and NCSN approaches, at least three possibilities exist for heat flow into a volume element. For the CSN case: (1) no heat flow, (2) heat flow via conduction, and (3) convection at the surface. For the NCSN case: (1) no heat flow, (2) heat flow via conduction, and (3) convection at the surface followed by conduction to the node. Other modes of heat transfer are possible, such as radiation, but the three mentioned are of primary importance in thermal processing. Volume element type C for the CSN approach (Fig. 31 and Table 32) and the NCSN approach (Fig. 32 and Table 32), exhibits all three of the respective heat transfer possibilities. Derivation of finite difference equations using energy balances for both approaches are demonstrated below (for complete code see Appendix B): Capacitance Surface Nodes (CSN) The node is located on the surface of the finite cylinder. As viewed in Fig. 31 and in the diagram for node type C in Table 31, convection occurs through the right boundary of the volume element. Conduction occurs through the top and left boundaries, and no heat is transferred from the bottom boundary. The energy balance for this node is (i = Nr, j = 1) p Cp dVi. jT i T n_ k *rAi 1.i j + k yAi(Tn+ TP) dt ijdr i dy I ~j1 (34) +h. rAi. j(TTr.) where p, cp, and k are density, heat capacity, and thermal conductivity of the contents of the finite cylinder, rA_ 1j and rAi.j are the inner and outer side boundary areas of the volume element containing node (i, j), yA1 is the area of the top and bottom of the volume element, dVi.j is the volume of the element containing node (i, j), TILj and Tnj+i are previous temperatures of nodes immediately to the right and above node (i, j), Tf, is the ambient temperature, and T{,j and TifW are previous and current temperatures of node (i, j). For appropriate nodes, subscripts (i + 1, j) and (i, j 1) refer to nodes immediately to the left and immediately below node (i, j), respectively. Solving Eq. 34 for Ti'.1 yields Tj+l' = Tj (1 cL cT cR) + cL TI.j + cT. Tj + cR T (35) where cL, cT and cR represent coefficients related to heat flow directed toward the node from the left, top and right, respectively (definitions of cL, cR, cT were provided in Table 31). The coefficient, cB, was used for volume elements with heat flow from the bottom. The ratio k / (pcp) is thermal diffusivity, denoted by the Greek letter, ta. In order for Eq. 35 not to violate the laws of thermodynamics (i.e. specifying heat flows from cold to hot regions), the sum of terms within parentheses (1 cL cT  cR) must not be negative. Solving this for the time step, dt, provides dt~oe c < (36) 1 krA i. + k yA...j' + h. rAi.( p Cp dVi. j dr dy ) Eq. 36 shows that as the heat transfer coefficient, h, increases, the time step, dt, must decrease. NonCapacitance Surface Nodes (NCSN) In this case, the node within volume element type C is at a distance dr/2 from the surface (Fig. 32 and Table 32). The temperature at the node is affected by conduction through the left and top boundaries, and by convection and conduction from the right boundary. The energy balance is written as follows (i = Nr 1, j = 1): p dVTi.j (Tnl Tij) kr Tii.j Ti.j) + Tij+i Tinj T Tqj (37) + 1 dr/2 h .rAi+i,j k rAi.j Solving Eq. 37 for T.t1 yields Tj = Ti.j" (1 cL cT cR) + cL. TPi.j + cT. Ti.j+ + cR T( (38) Eq. 38 is similar in form to Eq. 35, but the coefficient, cR, includes terms related to both convection and conduction. Similarly, the sum of the terms in parentheses must be nonnegative. Thus, the maximum time step for the NCSN case is given by d NCNS (39) tNode C  1 k rAii.j +hk r yAu.j +k 1+I p cp dVi. j dr dy 1 dr/2 + hrAi+ij krAi.j / Aside from describing the physical situation more appropriately, the superiority of the NCSN approach is clearly demonstrated by Eq. 39. As the heat transfer coefficient, h, increases, the term representing the resistance to convective heat transfer becomes less significant; this limits the impact of the heat transfer coefficient, h, on the time step, dt. Ifh is considered to be infinite, as is often done (Shin and Bhowmik, 1990; Teixeira et al., 1969), the time step restriction reduces to that of a specified temperature on the boundary (Type I or Dirichlet boundary condition). Once interior node temperatures are known, temperatures of surface nodes may be calculated by considering the following heat flux balances: Side surface (i = Nr 1, 1 j < Ny 1): k r T]i") = h* rAi +j, (Ti+1 T+l). (310) Top surface (1 < i < Nr 1, j = Ny 1): k* yAi, (T +l T1)n+l = h ,A +l (TPt T ()+I dy i.j 2.+ lhya, i.J+ (0) Side and top surface temperatures are determined by solving Eq. 310 and Eq. 311 for, T4tji and T., respectively. Notice that the time step, dt, does not appear in Eq. 310 or Eq. 311. This demonstrates that surface temperatures are not required to specify internal temperatures, however, surface temperatures may be computed as desired. Energy balances for all nine types of nodes were performed in a similar manner. Resulting coefficients and equations for the CSN and NCSN approaches were provided in Tables 31 and 32, respectively. The maximum allowed time step depends on (1) the number of nodes, (2) can dimensions, (3) thermophysical properties of the product, and (4) convective conditions. Thus, a time step no larger than the smallest maximum time step from all nodes should be used. Once the smallest maximum time step is specified, heat transfer simulation may be performed. Before entering the simulation loop, node temperatures are set according to a known initial temperature profile. When the simulation loop is entered, the ambient temperature is retrieved, time is advanced by the time step, and the energy balance equations provide temperatures of respective nodes. Figure 34 shows an example of programming code for performing these steps. Results and Discussion It may be shown that more complicated algorithms offer smaller degrees of error in nodal temperature estimates due to higher order differential approximations (Press et al., 1992; Clausing, 1981). However, when properly employed, explicit methods have been shown to be sufficiently accurate for the purpose of thermal process calculations (Thibult, 1985). In order for the NCSN approach to offer value in thermal process design, it must provide considerable time savings as well as comparable accuracy to the CSN approach. Alleviation of the time step restriction afforded by the NCSN approach is demonstrated in Table 33. Several standard industrial can sizes were selected for comparison (Miltz, 1992), as well as a plastic can, similar to that studied by Shin and Bhowmik (1990). High and low heat transfer coefficients, h, were chosen to represent external heating and cooling by condensing steam and cooling water, respectively (Walas, 1988). Table 33 shows that the maximum time step allowed for the CSN approach was consistently much lower than that for the NCSN approach. Additionally, high h values further restricted the time step for the CSN approach. The fact that the time step was unaltered for the NCSN case demonstrated that the smallest maximum time step was specified by an interior node (for both values of h). Indeed, it is often observed that the volume element containing node type A (Fig. 32) provides the smallest maximum time step due to its smaller dimensions. Most traditional cans have thin metal walls that are often ignored in heat transfer simulations. This is due to their relatively insignificant resistance to heat transfer, as compared to the contents of the can. However, retortable plastic cans are becoming increasingly popular. Shin and Bhowmik (1990) suggested that ignoring the resistance of the plastic can may be inappropriate. They proposed a finite difference method to model this situation using two sets of nodes, one for the can wall and one for the food. The researchers used this approach because they felt that practical node spacing for the food was too course for the thin can wall. Both sets of nodes were assembled according to the CSN approach. Consequently, dimensions of volume elements in the wall became so small, the authors concluded that explicit finite difference methods were impractical, because of the time step restriction. The authors implemented a more complicated implicit method (which is unconditionally stable with regard to time step), and justified its use by stating that, "the time increment in the numerical equations can be large (10 s)" (Shin and Bhowmik, 1990, p. 170). Although the Shin and Bhowmik (1990) model generated reasonable solutions, the explicit NCSN approach would have offered a more practical and convenient method of analysis. Instead of going to great lengths to obtain unnecessary temperature distributions in the walls of the can, it is more appropriate to account for the effect of this additional resistance to heat transfer on the thermal response of the food inside the can. In this case, the last term of Eq. 37 would take the form T'pn Vpn T, T F (312) 1 dr,, dr / 2 (312) + + h rAow kw rAimw k rAi. j where rAo,,w is the outside area of the wall, dr,, is the thickness of the wall, kw. is the thermal conductivity of the can wall, and rAm,,,w is the area at the log mean radius of the wall. In most cases, dr, is relatively small, thus the term containing dr, would not be expected to significantly alter the maximum allowed time step. As seen in Table 33, the maximum time step for the plastic can using the NCSN approach was 103.1 seconds. An absolute comparison of time steps was not possible, because the number of nodes used by Shin and Bhowmik (1990) was not reported. However, an explicit finite difference method with NCSN clearly offered a practical time step. Additionally, with the same number of nodes and time step, the explicit algorithm with NCSN would execute faster than the implicit method used by Shin and Bhowmik (1990). Accuracy of the CSN and NCSN approaches were compared relative to an analytical solution for heat transfer in a finite cylinder (Fig. 35). The contents of the can were assumed to be initially uniform at 0C. It was assumed that the can was instantaneously exposed to an ambient temperature of 121 'C, with a heat transfer coefficient of 5500 W/(m2 .C). The can had a diameter and height of 60.2 and 33 mm, respectively. Figure 35 demonstrates that both methods offered comparable accuracy. Additionally, Fig. 35 shows that explicit methods offer sufficient accuracy for thermal process calculations. Incidentally, the time required to obtain solutions at 3, 10 and 20 minutes required 4.6, 15.1, and 30.1 seconds, respectively for the CSN approach, and 0.3, 1.2, and 2.3 seconds with the NCSN approach. The computer used was an IBM compatible (Zenon Corp., Industry, CA), equipped with a 90 MHz Pentium processor (Intel Corp., Santa Cruz, CA) and 16 megabytes of RAM. The program was written in Visual Basic (Microsoft Corp. Seattle, WA) running under Windows for Workgroups 3.11 (Microsoft Corp., Seattle, WA). Lethality Calculation from Discrete TimeTemperature Data Prediction of thermal profiles in foods is only half of the task in thermal process design. Thermally labile factors are usually dispersed continuously throughout the volume of the can. Ideally, prediction of the extent of any reaction requires integration of the temperature sensitive kinetic function over the thermal history for the entire volume of the can. An analytical heat transfer solution provides all of the information necessary to accurately perform this task. In the absence of an analytical solution, an appropriate approximation is required. The nature of spatially discrete data often requires the assumption that the temperature throughout each volume element is uniform and equal to the temperature at the node. However, nodal temperatures are known at only specific points in time. Thermally labile factors are affected continuously through time and therefore, between the instants in times when temperatures are known. Therefore, the question arises, for each volume element, how should discrete timetemperature data be used to predict the extent of reaction occurring during each time interval? For any given time interval, two temperatures are known for each volume element: (1) temperature at the beginning of the interval, T,.j, and (2) the temperature at the end of the interval, T!..+,1 The true path taken between these points is unknown, but is reflected in observed behavior of thermally labile factors. Fig. 36 shows several options available for predicting the extent of a particular reaction. Option a refers to the possibility of using temperatures at the end of each interval to predict lethality that occurred throughout the interval. Option b refers to the possibility of using temperatures at the beginning of each interval. Option c refers to the possibility of using the mean of beginning and ending temperatures. Option d refers to the possibility of integrating an interpolating function, such as a linear function. Option e refers to the possibility of selecting some temperature other than the mean temperature. Options a and b are the most straightforward, and are probably used most often. However, most researchers do not report this detail in their methodology. Option a was used by Teixeira et al. (1969). It turns out that options a and b are essentially equivalent, because the same temperature data series is used for lethality calculations, with the exceptions of the first and last data points. Although option d would provide the most accurate result, it is probably never used, because it would be difficult to justify the time required to not only approximate appropriate continuous functions, but also to integrate such functions for each volume element and every time step. Regardless of the actual temperature profile, a linear approximation becomes increasingly accurate as the time step is reduced. Incidentally, errors from all of the options decrease as the time step is reduced; thus, large time steps offered by implicit methods are not beneficial when the goal is extentofreaction prediction. Assuming that a time step is chosen such that a linear profile adequately describes the actual profile during each time interval, it is possible to define some temperature along that profile that results in an equivalent extent of reaction as the actual linear profile (see Chapter 2). This is the objective of option e. Since temperature affects the rates of many food related reactions in an exponential manner, reaction rates at temperatures above the mean temperature are exponentially faster than those below the mean (Labuza, 1979). Therefore, for typical thermally labile food constituents, it is expected that an equivalent temperature exists in the temperature range above the mean temperature for each time interval. Many food related reactions are apparently first order with Arrhenius temperature dependence. A first order rate process is defined by C exp{k.t} (313) Co where Co is the initial concentration of a reactant, C is the concentration at time, t, and k is the reaction rate constant (note that the rate constant, k, in Eq. 313 is different from thermal conductivity, k, in Eq. 34, Eq. 36, Eq. 37, Eq. 39 and Table 31 and Table 32). The Arrhenius equation describing the temperature dependence of the rate constant is defined by k = ko exp E.114Ea) kko.exp R Tabsolute where ko is a preexponential factor, Ea, is activation energy, R is the ideal gas law constant, and Tabsolute is the absolute temperature. An equivalent isothermal process may be defined for any dynamic process (see Chapter 2): c k. Ea dt ko exp e(tn+i tn) (315) C 0R lR T(t)J  where T, is the equivalent isothermal temperature that is represented by option e in Fig. 36. For a small time step, (tn tn), the linear temperature change approximation is T.+l T!j+ [slope] t[;n (316) i~j  i~jtil The slope in Eq. 316 is, of course, the temperature difference divided by the time step. The slope of the linear temperature change represents how rapidly temperatures change, or how dynamic the process is within a particular volume element for each time interval. The location of the equivalent temperature depends upon (1) the activation energy, Ea, of the reactant (2) the slope of the linear temperature change of a given volume element over a given time interval, (3) temperature level, T, and (4) the duration of the time step. Figures 37 and 38 show how the location of the equivalent isothermal temperature varies according to these parameters. Figure 37 is for a high activation energy constituent (500 kJ/mol), such as a bacterial spore population. Figure 38 is for a low activation energy constituent (50 kJ/mol), such as a vitamin (Lund, 1980). The horizontal axis represents the slope of the linear temperature change. High values represent highly dynamic thermal situations. Three time steps (1, 10 and 100 seconds) and three beginning temperatures (70, 120 and 150 C) were considered. As expected, Fig. 37 and Fig. 38 show that as temperature changes become less rapid (more constant or isothermal), the equivalent isothermal temperature approaches the mean value. This is intuitive since the average of two equal values is the same value; for T = Ti.j = Tlj T+T 2.T 2 2 T. (317) 2 2 The smallest time step (1 second) showed the least variation from the mean temperature. The largest time step (100 seconds) showed the largest variation. The higher activation energy value showed more variation from the mean than the lower value. In all cases, initial temperature was not too important. Simulations can be used to show that for practical food processes, time steps smaller than 10 seconds are required for a linear temperature change approximation to be valid for a majority of the process time. Although simulations of typical retort processes showed that slopes of temperature changes may be as high as 5 C/second, slopes for most volume elements, for a majority of the process time, are typically less than 1 C/second. Therefore, Fig. 37 and Fig. 38 show that for extentofreaction predictions it is appropriate to select temperatures at or within 10% above the mean, for each timeinterval. Table 31. CSN finite difference formulas (Fig. 31). Type Diagram Node(s) Coefficients A A i =1 cR = acAt / dV(i, j) rA(i,j) / dr j = 1 cT = aAt / dV(i, j) yA(i) / d, T ' = Tj" (1 cR cT) + cR Ti.j + cT Tl^,_ B i i=2..Nr1 cL = aAt / dV(i, j) rA(i1, j) / dr j = 1 cR = caAt / dV(i, j) rA(i, j) / dr cT = caAt / dV(i, j) yA(i) / dy T. = T.j" (1 cL cR cT) + cL T Ij + cR TI.j + cT C  i=Nr cL = aAt / dV(i, j) rA(i1,j)/dr j = 1 cT = aAt / dV(i, j) vA(i) / dv /4 cR = a.At / dV(i, j) hrA(i. j)/k ._= T (1 cL cT cR) + cL C L jT_ + cT. T',j + cR.T _ D i = 1 cR = aAt / dV(i, j). rA(i, j) / dr j2.. Nv cT = aAt / dV(i,j) yA(i) / dy cB = caAt / dV(i, j) yA(i) / dy T'1 TI ." T (1 cR cT cB) + cR TrI.j + cT TO,+ + cB Tnj_1 E i=2.. Nr1 cL = caAt / dV(i,j) rA(il j) / dr j = 2 ..Nyi cR = At / dV(i,j) rA(i, j) / dr cT = aAt / dV(i, j). yA(i) / dv cB = aAt / dV(i, j) yA(i) / dy Tn1 = Tnj (1 cL cR cT cB) + cL T I.j + cR JIj% + cT Tn.j+, + cB. T ..i F  i=Nr cL = aAt / dV(i,j) rA(ilI j) /dr j = 22.. NvI T = cT=aAt / dV(i,j) yA(i)/dv cB = aAt / dV(i, j) yA(i) / dy .c5\ R = aAt / dV(i, j) h. rA(i. j) / k T = T'r (1 cL cR cT cB) + cL. T,1.j + cT* T,,. + cB T.ji + cR* T G i = l cR = caAt / dV(i, j) rA(i, j) / dr j = Ny cB = aAt / dV(i, j) yA(i) / dv cT = aAt / dV(i, j) h yA(i) / k Tf;1 = T (1 cR cB cT) +cR T + cB "T,. + +cT T+c H i=2.. Nr1 cL = a.At / dV(i,j) rA(i1, j) /dr j J= Ny cR = a.At /dV(i, j) rA(i, j) / dr cB = caAt / dV(i, j) yA(i) / dv cT = aAt /dV(i, j) h yA(i) /k T ,(1 & cLcR cB +cL+cR T..j + cB.l + cT. CB I i = Nr cL = aAt / dV(i, j) rA(i1, j) / dr j1 1 = Ny cB = caAt / dV(i, j) yA(i) / d" R cR =aAt / dV(i, j). h rA(i, j) / k cT = a.At / dV(i, j) h yA(i) / k Tl' = T (I cLcB cR cT) + cL'V, +cBT +cR I +cT'T Table 32. NCSN finite difference formulas (Fig. 32). Type Diagram Node(s) Coefficients A i= cR= a.At / dV(i,j) rA(i,j) / dr j = cT = aAt / dV(i, j) yA(i) / d V' j * (1 cR cT) + cR Tlj + cT j B i=2..Nr1 cL = aAt / dV(i, j) rA(i1, j) / dr j = I J lcR=a At / dV(i, j).rA(i,j) / dr cT = aAt / dV(i, j) yA(i) / dy T.; = Tj (1 cL cR cT) + cL TV ,j + cR Tri j + cT Tji C h] i=Nr cL = aAt / dV(i, j) rA(i1, j) / dr j = 1 cT = aAt / dV(i, j) yA(i) / dv ScR = aAt / dV(i j) hrA(i. j) /k T;= j (1 cL cT cR) + cL T, + cT. T,i + cR Tr D / i = 1 cR = caAt / dV(i, j) rA(i, j) / dr j = 2 ,. Ny1 cT = aAt / dV(i, j) yA(i) / dv cB = a.At / dV(i, j) yA(i) / dy Ti1 = Tij (1 cR cT cB) + cR I T'L.j + cT* T. + cB T" E 7 i=2..Nr1 cL = ca.At / dV(i, j) rA(i1, j) / dr j = 22.. Nyl cR = aAt / dV(i, j). rA(i,j) / dr cT = aAt / dV(i, j). yA(i) / dy cB = a.At / dV(i, j) yA(i) / dy Tn;1 = T.j (1 cL cR cT cB) + cL TIj + cR T 1 + cT TjI + cB. ji F  i=Nr cL = aAt / dV(i,j). rA(iL j) / dr Sj=2.. Nyi cT = aAt / dV(i,j). yA(i) / dv cB = aAt / dV(i, j). yA(i) / dyv cR = caAt / dV(i, j). h rA(i, j) /k V."1 = Tn, (1 cL cR cT cB) + cL ril,j + cT T.j+i + cB T, + cR T7 G i = I cR = aAt / dV(i, j). rA(i, j) / dr j = Ny cB = aAt / dV(i,j) yA(i) / dy V7 cT = aAt / dV(i, j) h v A(i) / k Ti' = T'.j (1 cR cB cT) + cR. T+. + cB T.j] + cT *T H i= 2..Nr1 cL = caAt / dV(i, j) rA(il, j) / dr j j= Ny cR = at.At / dV(i, j) rA(i, j) / dr ScB = aAt / dV(i, j) yA(i) / dv cT = aAt / dV(i, j) h vA(i) / k T". = T ," (1 cL cR cB cT) + cL TTL,+ + cR Ti,. + cB T,.i+ cT T2 i i=Nr cL = aAt / dV(i, j) rA(i1, j) / dr  j = Ny cB = aAt / dV(i, j) yA(i) / dv cR = aAt / dV(i, j) h rA(i, j) / k cT = a.At / dV(i. j) h vA(i) / k T '= T (1 cL cB cR cT) + cL T,1, + cB T,1 + cR. T5 + cT T7 Table 33. Comparison of maximum time steps allowed for heat transfer simulation via explicit finite difference algorithms. Application ofNCSN offers a significant time advantage over CSN. Can Sizea`b h # Nodes Maximum Time Step for Stability CNS NCSN W / (m2.oC) (seconds) (seconds) 54 x 73 mm 5700 8 1.70 121.20 (202 x 214) 570 8 15.67 121.20 5700 12 1.08 46.45 570 12 9.50 46.45 87 x 116 mm 5700 10 3.45 233.65 (307 x 409) 570 10 32.20 233.65 157 x 178 mm 5700 12 2.93 360.01 (603 x 700) 570 12 27.88 360.01 70 x 73 mm 5700 10 1.52 103.10 (Plastic) 570 10 13.88 103.10 a Values refer to diameter and height, respectively. Values in parentheses follow industrial standards; the first digit designates whole inches, and the remaining digits represent the number of sixteenths of an inch (e.g. 202 = 2 + 2/16 inches). b Dimensions of the plastic can are similar to those reported by Shin & Bhowmik (1990). H H H H I E E a Ezi i . .^, F D E i> j . PD 0 " E AL L E *0 D F E E E /  / .> /i 7D UI r / 7 '^ 9 117 lB 7/ B B B C yA(i) rA(i1,j) rA(i,j) <.R a d i u s Spatial descritization of a finite cylinder for CSN. A Figure 31. G 'D I n0 n ii yA(i) 7 r rA(i1,j) E l rA(i,j) 0) / D E"I, # F D F SE D .dr.  / E E E F " , r  I A B B B C < Rad i u s > Spatial descritization of a finite cylinder for NCSN. Figure 32. ReDim r(Nr), y(Ny), rA(Nr, Ny), yA(Nr). dV(Nr, Ny) ReDim rALocation(Nr), yALocation(Ny) Pi = 4 *Atnm(1) dr = 2 *Radius / (2 Nr 3) CSN[dr = Radius / (Nr 1)] dv = 2 HalfHeight /(2 Ny 3) CSN[dy = HalfHeight / (Ny 1)] NODES IN RADIUS r(1) =0: r(Nr) = Radius For i = 2 To Nr 1 r(i) = r(i 1) + dr Nexti NODES IN HALFHEIGHT y(1 ) = 0: y(Ny) = HalfHeight For i = 2 To Ny 1 y(i) = y(i 1) + dy Nexti LOCATIONS FOR SIDES OF VOLUME ELEMENTS rALocation(l) = r(2) / 2: rALocation(Nr) = Radius For i = 2 To Nr 1 rALocation(i) = (r(i + 1) r(i)) / (Log(r(i + 1)/ r(i))) Next LOCATIONS FOR TOPS AND BOTTOMS OF VOLUME ELEMENTS yALocation(1l) = y(2) / 2: yALocation(Nv) = HalfHeight Forj=2ToNy I yALocation(j) = yALocation(i 1) + dy Nextj VERTICAL AREAS 1 TO RADIAL HEAT FLOW For j = 1 To Ny For i = 1 To Nr rA(i, j) = 2 Pi rALocation(i) (yALocation(j) yALocation(i 1)) Next Nextj HORIZONTAL AREAS IL VERTICAL HEAT FLOW yA(1) = Pi rALocation(l) rALocation(l) yAsum = vA( 1) For i = 2 To Nr I CSN[For i = 2 to Nr] yA(i) = Pi rALocation(i) rALocation(i) yAsum vAsum = yAsum + yA(i) Next i VOLUMES Forj = 1 To Ny 1 CSN[Forj = 1 to Ny] For i = 1 ToNr 1 CSN[For i = 1 to Nr] dVol(i, j) = yA(i) (yALocation(j) yALocation(i 1)) Next Next j Figure 33. Visual Basic code for discretizing a finite cylinder using NCSN. Code in square brackets may be substituted to yield CSN. Do TimeNow = TimeNow + TimeStep AmbientTemperature = GetAmbientTemperature() HeatTransferCoefficient = GetHeatTransferCoefficientO Forj = 1 to Ny I CSN[Forj = 1 to Ny] For i = 1 to Nr 1 CSN[For i = 1 to Nr] Select Casej Case 1 Select Case i Case 1 (insert equation for node type A, Table 2 [Table 1]) Case 2 to Nr 2 CSN[Case 2 to Nr 1] (insert equation for node type B, Table 2 [Table 1]) Case Nr I CSN[Case Nr] (insert equation for node type C, Table 2 [Table 1]) End Select Case 2 to Ny 2 CSN[Case 2 to Ny I] Select Case i Case I (insert equation for node type D, Table 2 [Table 1 ]) Case 2 to Nr 2 CSN[Case 2 to Nr I] (insert equation for node type E, Table 2 [Table I]) Case Nr 1 CSN[CaseNr] (insert equation for node type F, Table 2 [Table 1 ]) End Select Case Nv 1 CSN[Case Ny] Select Case i Case I (insert equation for node type G, Table 2 [Table 1]) Case 2 to Nr 2 CSN[Case 2 to Nr 1] (insert equation for node type H, Table 2 [Table 1 ]) Case Nr 1 CSN[CaseNr] (insert equation for node type I, Table 2 [Table 1 ]) End Select End Select Next i Nextj Forj = I to Nv 1 CSN[For j = 1 to Ny] For i = 1 to Nr 1 CSN[For i = 1 to Nr] OldTemperature(i, j) = NevTemperature(i, j) Nexti Next j Figure 34. Visual Basic code for main loop of explicit finite difference heat transfer simulation using NCSN. Code in square brackets is for CSN. 130 120 110 110 20 minutes 100 90 70 60 60 minutes 50 0 Capacitance Surface Nodes 40 13 Noncapacitance Surface Nodes 30 Analytical Solution 20 10 0~ 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Dimensionless radial location on midplane (r/Radius) Figure 35. Temperature distributions along the radius on the midplane of a finite cylinder at 3, 10 and 20 minutes. NCSN offers comparable accuracy to CSN, but with significant time savings. Time Possible temperature selections for calculating extent of reactions (lethality) occurring during each time interval, in each volume element. T  avg *^ es &* u & I I. I Figure 36 1.0 0.9 0 0 0.8 E 0.7 0.6 0.5 0.4 0. Figure 37. 0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 Slope of Linear Temperature Change (C/sec) Location of isothermal temperature that provides an equivalent extent of reaction (reactant Ea = 500 kJ/mol) as a linearly approximated non isothermal exposure, between specified initial and final temperatures. 1.0 50 KJ/mol 0.9 W. 0.8 120 C 0.7 150 C S0. Time Step= 100 seconds  S ". . ......  70oC 0.6 ow.n 5 10 seconds d .0 M  0.5 1 second 0.4 0.0 0.5 1.0 1.5 2.0 25 3.0 3.5 4.0 4.5 5.0 Slope of Linear Temperature Change (C/sec) Figure 38. Location of isothermal temperature that provides an equivalent extent of reaction (reactant Ea = 50 kJ/mol) as a linearly approximated non isothermal exposure, between specified initial and final temperatures. CHAPTER 4 KINETIC PARAMETER ESTIMATION IN FOODS THAT DO NOT HEAT UNIFORMLY Introduction Selection of a target microorganism is a critical step in thermal process design. To the process designer, the name of the organism is secondary in importance to the kinetic parameters that reflect its rate of inactivation as a function of temperature. It is well known that the Arrhenius equation is often appropriate for describing the temperature dependence of the reaction rate constant, k, for food constituents and/or microbial inactivation (Cohen and Saguy, 1985). The Arrhenius equation is written as f Ea (41) k = k. exp Tabsolute I where the kinetic parameters consist of the preexponential factor, k. [min], and the activation energy, Ea [J/mol]. These parameters allow determination of the reaction rate constant, k [min'], at any temperature, T [Kelvin]. The ideal gas law constant, R, has the value 8.314 J/(molKelvin). Thermal process design often involves the use of generally accepted kinetic parameters for predicting microbial lethality or constituent losses. Perkins et al. (1975) pointed out that for the case of low acid canned foods, it is customary to use kinetic parameters for Clostridiumni bolulimim spores that were determined in the early part of the current century (Esty and Meyer, 1922; Townshend, 1938). These data were derived from traditional isothermal experiments on bacterial spores in aqueous suspensions. The practice of using generally accepted kinetic parameters in thermal process design raises an important question, are the parameters relevant to the product and process at hand? The process engineer must rely on the notion that the parameters represent expected behavior, at least under the conditions studied. Thus, the question addresses the well known fact that bacterial spore inactivation behavior depends upon at least four factors: (1) conditions under which sporulation occurred (Sedlak et al., 1995); (2) conditions of spore recovery (Fernandez et al., 1995; Gonzalez et al., 1995; Aldaton et al., 1974); (3) physical treatment of spores prior to recovery (Mallidis and Drizou, 1991; Gould et al., 1985); and (4) the chemistry of the medium in which the spores are suspended (Ahmed et al., 1995; Juneja et al., 1995; Fernandez et al., 1994; Chai and Liang, 1992; Lynch and Potter, 1988, 1989; Gould, 1985; Perkins et al., 1975). Consequently, the practice of using generally accepted kinetic parameters in thermal process design should be reassessed. Experimental conditions of sporulation and spore recovery may be controlled through adoption of standard methodologies. Exposure of spores to nonthermal, physical injury is not only controllable, but represents an exciting opportunity for product improvement through hybrid process development, geared toward thermal process minimization. Food chemistry, however, is both product and process specific. Obviously, the chemistry of different foods is different. Also, since food chemistry is very sensitive to environmental factors, different processing conditions result in different degrees of chemical change for a given product. Accounting for nonthermal, physical injury of spores and/or variations in food chemistry may permit production of higher quality thermally processed products due to reduced thermal processing requirements. In either case, sufficient impetus exists to seek a better understanding of product and process specific thermal inactivation behavior. Increasing the complexity of microbial inactivation models is one approach to accounting for nonthermal injury and/or variations in food chemistry (Whiting and Buchanan, 1994). However, practical difficulties arise, since complex models involve a greater number of parameters, which become increasingly difficult to estimate (Sapru et al., 1992). Although complex models may provide important qualitative insight under certain experimental conditions, Welt et al. (1995) found only marginal improvement in lethality predictions in retort processed conduction heating foods, when compared to predictions from the simple first order model with apparent parameters. Thus, continued use of the simple first order model for thermal process design is justified. Another approach for dealing with variations in observed inactivation kinetics is to use simple models (i.e. apparent first order kinetics with Arrhenius temperature dependence), but with product and process specific kinetic parameters. In order to be practical, simple and reliable methods are required for determining kinetic parameters from actual food products subjected to realistic processing conditions. Two approaches are commonly used to determine thermal kinetic parameters, isothermal and nonisothermal experiments (Lenz and Lund, 1980). The isothermal approach is often criticized for several reasons: (1) limited temperature range from which to select sufficiently different isothermal conditions results in too few conditions to provide good statistical confidence in the parameters (Cohen and Sugey, 1985), (2) unavoidable thermal lags during heatup and cooldown, (3) many small samples are tedious to prepare, handle and assay, (4) since samples must be small in order to minimize thermal lags, an ideal liquid medium such as a buffer solution, is often used rather than an actual food medium, and (5) even if an actual food product is used, the question of relevance exists, since experimental conditions often differ from processing conditions. The nonisothermal approach is generally more complicated than the isothermal approach. However, the nonisothermal approach can be applied to uniform and non uniform heating situations. Uniform heating may be considered to exist in small samples (as in the isothermal approach), well mixed reaction vessels (Welt et al., 1994), or continuous flow systems with turbulent flow (plug flow) conditions. Nonuniform heating occurs with large samples, where time constants of various modes of heat transfer are significant, as with thermally processed canned foods. For the case of uniform heating, nonisothermal parameter estimation methods overcome the problem of thermal lags, because the entire exposure is considered to affect the result. Nonisothermal methods that account for nonuniform heating alleviate the need for many small samples, since one sample is typically an entire product. Additionally, the observed mass average result is a function of the complete thermal profile history of the product/package system. Nonlinear optimization techniques are often employed for nonisothermal parameter estimation. Chemical engineering literature abounds with approaches and techniques (e.g. Watts, 1994; Lobo and Lobo, 1991; Kim et al., 1990; Espie and Macchietto, 1988; Biegler et al., 1986; Agarwal and Brisk, 1984; Nash, 1977; Marquardt, 1959, 1963). However, most of those reported were applied to well mixed, uniformly heated systems. Application of nonlinear parameter estimation techniques to nonuniformly heated foods has been demonstrated by Nasri et al. (1993) and Welt et al. (1995). Overall complexity, however, probably limits widespread use of nonlinear optimization methods. A nonisothermal approach, applicable to the case of uniform heating, was proposed by Swartzel (1984). However, Maesmus et al. (1995) and Welt et al. (1996) found that a fundamental assertion critical to Swartzel's method was invalid. Welt et al. (1996) proposed a new method, referred to as the method of Paired Equivalent Isothermal Exposures (PEIE). The PEIE method was shown to define kinetic parameters accurately from nonisothermal conditions for the case of uniform heating (Chapter 2). The purpose of this work was to extend the concept of the PEIE method to nonuniform heating. The goals of this work were (1) to establish the mathematical basis for applying the PEIE method to nonuniform heating, (2) to demonstrate the efficacy of the method by estimating known kinetic parameters from simulated data; and (3) to apply the method to a canned food product inoculated with thermophilic bacterial spores. Theory The following applies to the case ofreactants that follow apparent first order kinetics with Arrhenius temperature dependence. A first order rate process is defined by dC SC kC (42) dt where, C is the concentration of a reactant, t is time, and k is the rate constant as defined by Eq. 41. Integration of Eq. 42 yields ln(C/Co)=k.(t to) (43) Thus, extent of a given reaction under isothermal conditions may be determined by substituting Eq. 41 into Eq. 43: lnC/0)~k.eP Ea }(t t.) (44) ln(C / Cj = k, exp R. o. (t t) (44) Under nonisothermal conditions, Eq. 42 must be integrated over the entire thermal history, T(t): ln(C / C.) = ko. J exp,Ea dt (45) For brevity, the integral in Eq. 45 is referred to as G (Chapter 2): G = Jexp dt (46) 0 R. T(t) Clearly, when the thermal history, T(t), is known, a distinct value of G applies to a given value of Ea. Additionally, for a given T(t) and Ea, an infinite number of equivalent isothermal processes (t,, Te) yield the same value of G via Eq. 46: G = te .exp{E} (47) Here te and Te represent the equivalent time and isothermal temperature, respectively, that provide the same extent of reaction as the nonisothermal process, for reactants characterized by the particular value of Ea. Since the value of G is calculated from Eq. 46 using a known thermal history and a known or assumed E, value, Eq. 47, with unknowns, te and Te, represents an infinite number of equivalent isothermal processes. This infinite set of equivalent isothermal processes may be observed graphically by plotting a rearranged form of Eq. 47: ln(t,) = ln(G) + E (48) R T, A plot of In(te) versus 1/Te yields a line with ln(G) as the intercept and Ea/R as the slope. A low value of Ea results in a line with a shallow slope, while a large Ea results in a steep slope. Specification of two reactants, characterized by Eai and Ea2, yields two curves from Eq. 48, for a given nonisothermal process. The point where the curves intersect may be found by solving the two equations simultaneously (Chapter 2). This point of intersection is the isothermal process, (te, Te), that would yield the same extents of reaction as the nonisothermal process for reactants characterized by the two activation energies, Eai and Ea2. This point is refereed to as an Equivalent Isothermal Exposure (EIE). For the case of uniform heating, it was shown in Chapter 2 that this relationship could be exploited to find unknown kinetic parameters from at least one pair of non isothermal experiments. The method, referred to as the method of Paired Equivalent Isothermal Exposures (PEIE), is depicted in Fig. 41. The figure depicts two iterations of the method for two nonisothermal processes (indicated by "Process 1" and "Process 2") for which extents of reaction for a given reactant are known. For the first iteration, two Ea values are arbitrarily selected (Eal a and Ealb). For each dynamic thermal exposure considered, the lower Ea value, Ealb, yields curves with shallower slopes, while the higher value, Eal a, yields curves with steeper slopes. Points of intersection of respective curves, represent EIE's for these two arbitrarily selected Ea values. An estimation of the apparent or true Ea value is obtained by combining the EIE specifications (t,, T.) to previously measured extentofreaction data. New Ea guesses, Ea2a and Ea2b, are obtained by replacing Eala with the calculated Ea estimate and setting Ea2b to either a fixed value or some multiple of Ea2a. Schemes for specifying E,2b were discussed in Chapter 2. The procedure is repeated until El stops changing. Kinetic parameters, Ea and ko, are obtained in a typical manner by constructing an Arrhenius plot from the final equivalent isothermal process specifications and extentofreaction data. The PEIE method generates a set of kinetic parameters for each pair of non isothermal experiments performed. In the context of the isothermal approach, this is equivalent to drawing a line through each pair of points plotted on an Arrhenius plot (ln(k) versus 1/Tabsolute). Thus, for n nonisothermal experiments, n(n1)/2, potential sets of kinetic parameters may be defined. For example, for 4 nonisothermal experiments (a, b, c, d), a maximum of 6 sets of kinetic parameters may be expected (one from each pair: ab, ac, ad, bc, bd, cd). Application of the PEIE method for situations ofnonuniform heating was developed for the case of conduction heating canned foods. Complications arise because nonuniform thermal profiles result in nonuniform concentration profiles of reactants within such products. Also, in actual experiments it is often only possible to measure final mass average concentration of reactants (e.g. bacterial spores). Thus, a means to predict transient thermal profiles in a product, and to estimate resulting concentration profiles of a reactant are essential components of the PEIE method for the case of nonuniform heating. Numerical heat transfer algorithms are often used to obtain transient thermal profiles throughout a product. An indepth discussion of such methods for the case of a conduction heating finite cylinder was provided in Chapter 3. Briefly, numerical heat transfer algorithms are derived by considering that the product is comprised of discrete volume elements (Fig. 42). Accurate temperature estimates are obtained for a specific point, or node, in each volume element at specific points in time. Since volume elements are relatively small, it is assumed that temperatures at the nodes are applicable to the entire volume of respective volume elements. Thus, volume elements are considered to heat uniformly according to nonisothermal processes specified by temperature histories of respective nodes. In order to apply the PEIE method to nonuniform heating when only final mass average concentration of a reactant is known, additional considerations are required. Initial concentrations of reactants in volume elements are typically equal to the initial concentration of the entire product. Since numerical heat transfer analysis yields a non isothermal process for each volume element, final concentrations of a reactant in all volume elements may be calculated. Final volume element concentrations may be converted to overall mass average concentration by multiplying each final concentration by the volume of the element (giving the absolute amount of reactant in each volume element), summing these values over all elements, and dividing by the total volume of the object. Specifically, each volume element (i, j) receives a particular thermal history, T(t)ij. Thus, values of Gij may be calculated for each volume element, for assumed Ea values. Rearranging Eq. 44 and Eq. 45 yields _._j = exp{k.0 G i} (49) Co0 where Ci.j is the concentration in volume element (i, j) at any given time. Multiplying each side by the volume of the respective element AVij, and the known initial concentration, Co, yields Ni, =C, AVij exp{ko Gi,j} (410) where N.j is the absolute amount of reactant in volume element (i, j). Summing Eq. 4 10 over all volume elements and dividing by the sum of volume element volumes, yields Nj C" .(AV.'j exp{k .G}) (4) Z AV1 I AVj where the left hand side is the discrete definition of mass average concentration. This provides the basis for comparison with actual experimental results involving final mass average concentrations. For a particular thermal exposure, and selected value of Ea, the preexponential factor, k0, is the only unknown in Eq. 411. Thus, k1 may be found with any convenient root solving algorithm. Since the value of k, is typically large, it is convenient to work with the equivalent reference rate constant form of the Arrhenius equation (Eq. 41). This form of the Arrhenius equation is derived by substituting a reference temperature, Tr, and reaction rate constant at the reference temperature, kr, into Eq. 41, and then subtracting the resulting equation from Eq. 41. The unwieldy ko value cancels out and the resulting equation takes the form k = kr exp Ea FF, T (412) where kr is the reaction rate constant at temperature, Tr. Thus, k1 may be replaced with kr in Eq. 49 through Eq. 411. Since food engineers are familiar with the orders of magnitude of rate constants of food constituents at 121.1 C, this temperature serves as a convenient reference temperature. Since the value ofk121.1Ic must be positive, and its order of magnitude is generally known, bisection may be used to determine its value to a specified level of precision (Press et al., 1992). This additional step is the key to applying the PEIE method to the case of nonuniform heating. Utilization of the PEIE method for dynamic thermal treatment of samples that heat nonuniformly, involves the following steps, which are illustrated in Fig. 43: 1. Record initial concentration, ambient thermal history, and final mass average concentration from at least two dynamic thermal processes that result in statistically different extents of reaction (Processes A and B, Step 1, Fig. 43). 2. Arbitrarily select two Ea values. Simulate thermal profiles from actual ambient process conditions. This provides unique thermal histories for all volume elements (Step 2, Fig. 42). During simulation, apply the Ea values to accumulate G values according to Eq. 47. Thus, for two processes, A & B: Eai and Ea2 give G1 and G2 for each volume element, i, j, for process A and G, and G2 for each volume element, i, j, for process B. 3. Considering each process separately, use respective experimentally determined mass average concentrations to find reference rate constants (kri & kr2) that correspond to the two Ea values. This is done by applying a root solving algorithm, such as bisection (Press et al., 1986) on Eq. 411. 4. For each experiment, select a representative volume element, (i, j), and estimate the resulting concentration ratio, Cij / Co, for this volume element, for both sets of kinetic parameters (Eaikri and Ea2kr2). 5. Using the two Ea values (from Step 2) and corresponding Gij values (from Step 3), determine Equivalent Isothermal Exposures (EIE: te, T,) for the selected volume element (from Step 5), by solving the two equations (form of Eq. 47) simultaneously. Using Eq. 43, determine respective equivalent isothermal rate constants, kT. Since kinetic parameters, Eaikri, provide (Ci,j / Co)1 and parameters, Ea2kr2, provide (Ci.j / C0)2, exposure for time, t., at the equivalent isothermal temperature, Te, would result in the same concentration ratios for the respective Ea values. 6. Using the isothermal rate constants from Step 5 from both processes, calculate new respective Ea values. Thus, Eai applied to processes A and B yields a new Ea. value. Use these new Ea values for the next iteration (return to Step 2). However, before repeating, test stopping criteria (e.g. stop when both new Ea values are not more then 1 J/mol different than the old E, values). Once the stopping criteria are satisfied, take the average of the last two sets of kinetic parameters as apparent values for that experimental pair. When apparent kinetic parameters are defined for a given pair of experiments, select another pair of experimental results and return to Step 2. The PEIE method is performed repeatedly, on all pairs of experimental data, where data from one experiment consists of a recorded ambient thermal history and final mass average concentration of a reactant. Methods Numerical Heat Transfer The explicit finitedifference method described in Chapter 3 was used to simulate heat transfer in a finite cylinder. For each time interval during simulation, a temperature value 5% above the mean was used to calculate extent of reaction in each volume element during each time interval. Can dimensions, thermophysical properties of the food, applicable heat transfer coefficients, and parameters related to numerical heat transfer simulation were provided in Table 41. Sources and methods for obtaining these data are given below in the section related to actual experiments. Experimental Simulated experiments Simulated experiments were performed to demonstrate the use of the PEIE method for the case of nonuniform heating. A food constituent was considered to be characterized by the following Arrhenius kinetic parameters: k121.1ic= 0.2 min1, and Ea = 300 kJ/mol. The constituent was considered to exist uniformly in a can of conduction heating food at an initial concentration of 1 x 106 units/ml. Six simulated retort experiments were considered (Table 42). During retort comeup and heating phases, heat transfer coefficient, hHeat (Table 41), was used. During retort cooldown, hCool (Table 41), was used. When the retort temperature reached 30C, the can was considered to have been placed into an ice bath (hIce, Table 41). Simulations were halted when the temperature at the center of the can reached S10C. Final mass average concentrations were recorded in association with respective thermal treatments (Table 42). The ability of the PEIE method to predict the Arrhenius parameters that were used to generate these observed data was tested. For each pair of processes, iterations were stopped when new E, values (see Step 6 above) changed by less than 1 J/mol. The average of the last E, and kr estimates were recorded. Since 6 experiments were considered, a maximum of 15 sets of parameters were expected. Complete code of the PEIE method for the case of conduction heating canned foods is provided in Appendix C. Actual experiments Actual experiments were performed with canned pea puree, inoculated with spores of Bacillus stearothermophilus (ATCC 7953). Frozen peas were purchased from a local supermarket. The peas were freeze dried and crushed into powder. The powder was placed in home canning jars7 sterilized with ca. 50 kGy in a 10 MeV electron beam irradiator (Florida Linear Accelerator, Gainesville, FL), and stored at approximately 20C prior to use. Pea puree was prepared by reconstituting the freezedried pea powder with distilled water (80% distilled water and 20% powder by weight). The spore suspension was included in the reconstitution formula, and was added after the powder and water were well mixed. The puree had an initial pH of approximately 6.8. Measured initial concentration of spores in the puree was 9 ( 1) x 105 CFU/ml. Cans were cold filled with 95.0 ( 0.5) grams of inoculated puree. Air pockets were expelled from the puree by lightly tapping each can prior to sealing. Absence of noise when shaking indicated minimal headspace remained in the cans. Each can was uniquely identified with permanent ink, and stored in an ice bath for at least 48 hours prior to use. Thermocouples located at the center of several cans, indicated that the center reached 0C within 24 hours. Thus, initial temperature profiles of all samples were considered to be uniform at 0C. During thermal treatments, each sample experienced four ambient conditions. Each condition was reflected in the heat transfer coefficients applied during simulation (Table 1). The four conditions were: (1) after removal from the ice bath and before application of steam in the retort (hAir), (2) during steam heating in the retort (hHeat), (3) during water cooling in the retort (hCool), and (4) placement in an ice bath (hlce), after removal from the retort. A reasonable value of hHeat was assumed since a precise value is not required when it is known to be very high. Values of hAir, hCool, and hlce were inferred by measuring the thermal response of a solid aluminum cylinder to conditions similar to those experienced during treatment. Dimensions of the cylinder were similar to those of the actual cans (Table 42). Since the thermal conductivity of aluminum is very high and the cans were relatively small, the temperature of the object was considered to be uniform. Thus, a lumped parameter heat transfer solution provided the means to obtain values of the heat transfer coefficient. The lumped parameter solution has the form (Lienhard, 1987) T h.A exp t (413) TiT. P.cpV J where A and V are the area and volume of the aluminum cylinder, p and cp are the density and specific heat of aluminum (Perry and Green, 1984), T,, and Ti are the ambient and initial temperature of the aluminum cylinder and T is the temperature of the aluminum cylinder at any time, t. Equation 413 shows that the heat transfer coefficient, h, may be inferred from the slope of a semilog plot of accomplished temperature versus time. Cans were exposed in duplicate to different retort processes. In general, three approximate setpoint temperatures were used (104.4, 112.8, and 120.6C) in conjunction with various heating times. Ambient temperatures in the retort were recorded with at least three thermocouples in proximity to the cans. Recorded ambient temperature histories and resulting survivor ratios are shown in Fig. 44 and Fig. 45. Details of microbial enumerations were provided in Appendix D. Fifteen different processes were considered; the processes are numbered in the figures. Processes 4 and 6 (Fig. 41) were not plotted because they were similar to processes 3 and 5, respectively. Survival ratios for processes 4 and 6 were 0.272 and 0.068, respectively. Thermocouples were calibrated to a mercuryinglass (MIG) thermometer installed on the retort. The pilot scale vertical stillcook retort was supplied with 0.17 MPa (25 psig) steam in the pilot plant of the Department of Food Science and Human Nutrition at the University of Florida (Gainesville, FL). All samples were kept in an ice bath after treatment, and were assayed for surviving spores within 48 hours. The contents of each can were placed into a sterile 250 ml beaker and mixed well for several minutes by hand. The puree was sampled in triplicate by transferring ca. 0.5 grams of puree into respective test tubes with sterile, disposable 10 ml pipettes. Samples were weighed on a digital scale, and immediately diluted by a factor of 10 (by weight) with sterile 0.15 M potassium phosphate buffer (pH 7.0). Serial dilutions were subsequently performed. Four 50 pl aliquots of respective dilutions were placed on petri dishes containing nonselective media. The drops were dried on the surface of the media in a 55C incubator. The plates were arranged in groups of four and covered lightly with aluminum foil, inverted, and incubated at 55C for 24 hours. Colonies were manually enumerated. The recovery medium consisted of 8 grams of BBL nutrient broth (Fisher Scientific, Pittsburgh, PA) and 20 grams ofDifco Bacto Agar (Fisher Scientific, Pittsburgh, PA) per liter of distilled water. The medium was autoclave sterilized, and poured into sterile petri dishes (Fisher Scientific, Pittsburgh, PA). After cooling and hardening, the plates were stored in a refrigerator. All plates were used within two weeks of preparation. Results and Discussion The objective of the simulated experiments was to determine whether the PEIE method could identify Arrhenius parameters that were used to generate the observed data (Table 42). As mentioned above, the parameters used were k121.1oc = 0.2 min' and Ea = 300 kJ/mol. Since 6 experiments were considered, 15 experimental pairs, and thus, a maximum of 15 sets of parameters were expected. Table 43 summarizes results from all process pairs. Although some scatter was evident, the PEIE method was capable of defining the original parameters to a level of accuracy consistent with the observed data. As shown in Table 42, extentofreaction data were rounded from 15 (double precision) to 6 decimal places. Statistical analysis of Table 3 yielded kr = 0.199997 ( 0.000003) min1 and Ea = 299,998 ( 2) J/mol (intervals were determined at the 95% confidence level). Clearly, the PEIE method was successful in identifying the original kinetic parameters to an acceptable level of precision. Results from experiments on canned pea puree inoculated with spores of B. stearothermniophilus are shown in Fig. 46 and Fig. 47. Since 14 experiments were considered, a maximum of 91 experimental pairs, and thus, 91 sets of parameters were possible. However, Fig. 46 and Fig. 47 show 51 sets. Reasons for the apparent loss of data were (1) experimental uncertainty, (2) inactivation model simplification, and (3) conditions that contributed to loss in numerical precision. These factors are discussed in detail below. Appendix G provides tabulated values and the origin of all data in Fig. 46 and Fig. 47. Experimental uncertainty and model simplification have similar effects on kinetics analysis. Ideally (for irreversible, first order kinetics), as a thermal process becomes more severe, the number of surviving spores decreases. In general, this is observed, however, when comparing survivors from comparable processes, the more severe process may yield a greater mean survivor level. This may be the result of a more complex inactivation process such as activation of dormant spores (Sapru et al., 1992), or it may simply be the result of overlapping envelopes of experimental uncertainty. In either case, if this situation arises when a simple inactivation model is assumed (i.e. Eq. 44), a negative apparent activation energy results, which is not possible. In this experiment, spores were not thermally activated prior to use, and neither are they in practical situations. Thus, thermal activation of dormant spores was a probable factor. Additionally, microbial enumeration techniques are notoriously uncertain. Therefore, it is probable that both factors influenced the results. The PEIE method computer code in Appendix C was designed to simply bypass process pairs that generated negative Ea values during execution. This accounted for 6 of 91 possible data points. Each iteration of the PEIE method involves calculation of activation energy estimates from equivalent isothermal data. Therefore, the PEIE method is restricted in a manner similar to the traditional isothermal kinetics approach. In order to construct an Arrhenius plot from traditional isothermal experiments, reaction rates from at least two temperatures must be known. Thus, the PEIE method suffers when two dynamic thermal processes yield similar equivalent temperatures, T.. This problem was mitigated by two approaches. The first was to select the volume element that provided the greatest difference in equivalent temperatures during each iteration (see Step 5, above). The second was to consider only pairs of experimental processes that had significantly different thermal treatments. Thus, processes associated with one retort setpoint temperature were analyzed only with processes associated with different retort setpoint temperatures. For example, process 1 (Fig. 44) was not paired with processes 2 through 8, but was paired with processes 9 through 14 (Fig. 45). This accounted for 34 of 91 possible data points. As previously mentioned, the initial pH was approximately 6.8. Final mass average pH tended to decrease with increasing process severity, and ranged to as low as approximately 5.8. Histograms of observed ln(k12i.1oc) and Ea values showed normal distribution behavior. The observed logmean k121.1oc value was 0.26 min', with a 95% confidence interval for the population of 0.23 to 0.30 min"1. The observed mean Ea value was 250 kJ/mol, with a 95% confidence interval for the population of+ 15 kJ/mol. Fundamental differences between the PEIE method and traditional kinetic parameter estimation methods may preclude direct comparisons of estimated parameters. The following discussion will focus on the case of uniformly heating samples, since that is the typical premise of traditional isothermal kinetics experiments. Although microbial survivor curves often possess shoulders during initial phases of inactivation (Fig. 48), it is common practice to discard the shoulder, and obtain kinetic parameters from the straight portion of the curve. This is done by thermally activating spores prior to lethal exposure and/or by performing a linear regression on data only in the straight portion of the curve. As shown in Fig. 48, use of parameters derived from the straight portion of survivor curves will constantly over predict lethality to an unknown extent when shoulders are present during lethality. As implemented, the PEIE method only considers irreversible firstorder inactivation. Additionally, the PEIE method operates by connecting two points on a survivor curve (Fig. 48). Thus, in the absence of a shoulder on the actual survivor curve, the PEIE method would provide results similar to the traditional kinetic parameter estimation approach. When a shoulder does exist, the parameters obtained with the PEIE method will differ from those of the traditional approach. Reed et al. (1951) provided an excellent example of this difference. Amongst experiments involving several bacterial spores and suspending media, they performed experiments with pea puree inoculated with the same strain of spores used in this study. The authors reported the presence of shoulders in survivor curves, and reported results in the following two ways, (1) total time to achieve 99.99% reduction in spore count (four logcycle reduction) at 121.1 C, and (2) time required for the straight line portion to traverse four logcycles at 121. 1C. They reported 36 and 25.4 minutes, respectively, for the two methods; each with z values of 11.1 C. The times on a per logcycle basis are the definition of Do values (time required to achieve a one logcycle reduction in concentration. The Do value is related to the reaction rate constant, k121.1oc, via ln(10) Do = 0 (414) k121.1 The method used by Reed et al. (1951) was analogous to the result provided by the PEIE method. On a Do and z value basis, the values reported here (Do = 8.9 min; z = 11.4C) agreed very well with data from Reed et al. (1951) (Do = 9.0 min; z = 11.1 C). 