Kinetic parameter estimation in prepackaged foods subjected to dynamic thermal treatments


Material Information

Kinetic parameter estimation in prepackaged foods subjected to dynamic thermal treatments
Physical Description:
Welt, Bruce, 1967-
Publication Date:

Record Information

Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 23827072
oclc - 35790433
System ID:

Table of Contents
    Title Page
        Page i
        Page i-a
        Page ii
        Page iii
    Table of Contents
        Page iv
        Page v
        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 code--finite 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
        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







UNIVERSiTY OF F.,' 'A L. :,.-


(1c) i /( i/

-. SCI. NC.


The author would like to express sincere gratitude to his major advisor, Dr.

Arthur A. Teixeira and unofficial co-chairman 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


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.


ACKN OW LED GM ENTS ........................................................................................... ii

ABSTRACT ........................................................................................... ........ ........... vi

1. INTRODU CTION ................................................................................................. 1

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

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 on-Capacitance Surface N odes (N CSN )......................................................... 49
Results and Discussion......................................................................................... 51
Lethality Calculation from Discrete Time-Temperature Data................................ 54

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

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


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



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


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

non-uniformly under non-isothermal 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 x-rays at 0, 1.0 and 3.0 kGy

at an electron linear accelerator facility prior to thermal treatment in a vertical still-cook

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.


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


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


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 in-situ. 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.



The Equivalent Point Method (EPM) was developed by Swartzel (1982) in

order to facilitate design of continuous flow ultra-high 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

non-isothermal methods for kinetic parameter estimation to the food industry.

Two approaches are commonly used to obtain thermal kinetic parameters that

involve either isothermal or non-isothermal 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 first-order kinetics with Arrhenius

temperature dependence (Labuza, 1979). Since thermal process design often involves

the assumption of first-order kinetics with Arrhenius temperature dependence, there is

room for expanded use of the non-isothermal approach in kinetic parameter estimation.

Non-isothermal 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 non-isothermal methods is limited, however,

due to increased complexity of data analysis. Several non-isothermal 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 easy-to-use 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, in-and-of-itself, 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,


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 non-lethal (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



The functional form of a first-order rate process is

d- _kC (2-1)

where C is the concentration of a reactant at time, t, and k is the reaction rate constant.

Solving Eq. 2-1 by integration gives

C -=te (2-2a)

or equivalently,

ln(- )= -k.(t-to) (2-2b)

The temperature dependence of the rate constant, k, may be described by the Arrhenius


k= k0 -exp E.} (2-3)
LR-T j

where kinetic parameters ko and Ea are the pre-exponential 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. 2-3 into Eq. 2-2:

ln(C ) =-ko exp{- E. (t to). (2-4)
Co c Rex{ .(T t

In order to determine the extent of a reaction under non-isothermal conditions Eq. 2-4


ln(-) exp -Ea dt. (2-5)
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

right-hand side of Eq. 2-5 is equated to the right-hand side of Eq. 2-4. The isothermal

temperature, T, in Eq. 2-4 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. 2-4. 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. 2-5 to Eq. 2-4, but introduces the

extra step of normalizing the resulting equation by -ko. The resulting equation takes the


In( C
-E- G } exp E dt = t,-exp -Ea (2-6)
-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 straight-line equations of the form

ln(t) = ln(G)+ Ea (2-7)

resulted when these G values were equated to the isothermal part of Eq. 2-6 (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 (2-8)

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, (2-9a)

-ln(G2)=-ln(te)+- E.2

T (Ea2 Eal) (2-9b)
R.ln(Gj G2)

Equivalent time, te, is found by back-substitution 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. 2-6, Eq. 2-7 and Eq. 2-9. 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. 2-7

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. 2-7 by solving them

simultaneously. Since the two equations in Eq. 2-9a are independent and unique, their

solution, or point of intersection, must also be unique. It is not possible for Eq. 2-9 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. 2-6, using the actual thermal history.

3. Equate G values to the isothermal form of the equation (right most side of Eq.

2-6), using respective Ea values.

4. Generate straight line curves on a ln(te) versus l/To plot, for each Ea G


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 non-isothermally heated from 0C to 140C over a total process time of 200

minutes, according to the following linear equation (see Fig. 2-1, Process a):

T(Kelvin) = 273 + 0.7t (2-10)

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 = expR---Ea } dt (2-11)
0 R -R(273 + 0.7t)J

The integrals were numerically evaluated to a tolerance of 1 x 10-6, 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 2-1.

By plotting the three respective forms of Eq. 2-7, three points of intersection,

representing three different isothermal processes, were obtained for each pair of Ea

values (i.e. for E, value pairs 1-2, 2-3 and 1-3). 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. 2-2. 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 reaction-apparent 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 reaction-apparent

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


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 extent-of-reaction 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


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 extent-of-reaction data are exactly known. Implications of various

Ea2 specification schemes for the typical case of uncertain extent-of-reaction

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. 2-2b using the EIE

specifications (from step 3) and respective extent-of-reaction data (from step



5. For each pair of kT values determined in step 4, simultaneously solve Eq. 2-3, 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 1-2, 1-3, and 2-3).

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 extent-of-reaction

data from dynamic thermal experiments. Consider experiments involving thiamin in pea

puree using the three dynamic thermal exposures depicted in Fig. 2-1. Expected

extents of thiamin degradation from each exposure may be calculated by substituting

representative kinetic parameters into Eq. 2-5 (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 2-2. The PEIE method was applied to these

data as follows:

Step 1: See Table 2-2.

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 2-3).

Step 3: Equivalent isothermal exposures for each process pair (a-b, a-c and b-c)

were determined as previously described (see columns labeled "Te" and "t,"

in Table 2-3).

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 2-3:

k3a86435 =l-n(0.421)/ 101.194 8.548 x 10-3 min-'
377.662 3 -
k77662 in(0.600) /135.230 = 3.773 x 10-3 min-1 (2-12)
k368882 = ln(0.773)/160.974 = 1.600 x 10-3 min1

Step 5. Activation energy, Ea, estimates were calculated for each pair of kT values

by substitution into respective Arrhenius equations (Eq. 2-3). Solving

these simultaneously yields Eq. 2-13

R.ln(k k2)
Ea= -- ('2-13)
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 2-4). Details from the third and final

iteration are shown in Table 2-5.

Step 7: Regression analysis was performed on ln(kT) versus 1/Te data from the final

iteration (Table 2-5). Resulting kinetic parameters were Ea 113,000

J/mol and k, = 1.00 x 1013 min1.

Table 2-6 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. 2-6 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 2-6).

PEIE Example #2

The same dynamic thermal processes (Fig. 2-1) and kinetic parameters (Ea =

113 kJ/mol, ko = 1.OE+13 min-1) 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 2-2), plus or minus 5

percent. Equations of the thermal histories and obsen'ed extent-of-reaction data were

provided in Table 2-2 (column labeled 'Observed,' Table 2-2).

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 2-10. The best parameters given the observed data were Ea = 120,289 J/mol and

ko = 9.270E+13min"1.


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. 2-3. The magnitude of these differences

may be assessed by comparing extent-of-reaction predictions from these sets of

parameters to the exact experimental data in Table 2-6. The SSQ errors obtained using

parameters from the first and last iterations were 6.79E-02 and 1.93E-14, respectively,

indicating that the PEIE method provided closer agreement with exact experimental

data. Figure 2-4 and Table 2-10 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. 2-4), 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 2-7). 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.

2-4). Since there was uncertainty in the extent-of-reaction 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. 2-4).

Ideally, the value of Ea2 should not influence the result. Figure 2-5 is an

Arrhenius plot for example #1, where extent-of-reaction 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. 2-5 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 2-6 is an

Arrhenius plot for example #2 (uncertainty in extent-of-reaction 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. 2-6).

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. 2-6). Both extremes may artificially influence regression calculations,

resulting in poorer quality parameter estimates. As expected, all points in Fig. 2-6

tended to represent one particular Arrhenius curve, however, results from separate

regressions show some variation. Table 2-11 provides separate regression results and

resulting SSQ error values for all cases considered in Fig. 2-6. 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. 2-6 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 2-1. Ea and corresponding G values calculated from Eq. 2-8.

Index Ea Value G Value
1 30,000 8.914E-03
2 100,000 4.277E-12
3 400,000 1.281E-50

Table 2-2. Hypothetical dynamic thermal experiments used to demonstrate
calculations of the Paired Equivalent Isothermal Exposures (PEIE)
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 2-3. First iteration of the PEIE method for example #1. Progress from each
iteration is shown in Table 2-6.

Process Ea G Te te k Ea
Pair Guesses Estimate
(J/mol) (K) (min) (1/mmin) (J/mol)
a 30,000 8.914E-03
a 60,000 7.853E-07 386.435 101.194 8.548E-03
b 30,000 9.590E-03
b 60,000 5.137E-07 377.662 135.230 3.773E-03 113,101
a 30,000 8.914E-03
a 60,000 7.853E-07 386.435 101.194 8.548E-03
c 30,000 9.094E-03
c 60,000 5.137E-07 368.882 160.974 1.600E-03 113,162
b 30,000 9.590E-03
b 60,000 6.801E-07 377.662 135.230 3.773E-03
c 30,000 9.094E-03
c 60,000 5.137E-07 368.882 160.974 1.600E-03 113,220

Table 2-4. Second iteration of PEIE for example #1. Progress from each iteration
is shown in Table 2-6.

Process Ea G Te te k Ea
Pair Guesses Estimate
(J/mol) (K) (min) (1/min) (J/mol)
a 113,101 8.393E-14
a 226,201 2.141E-28 404.824 32.897 2.629E-02
b 113,101 4.947E-14
b 226,201 1.145E-29 395.204 43.938 1.161E-02 113,000
a 113,162 8.240E-14
a 226,324 2.065E-28 404.828 32.880 2.631E-02
c 113,162 2.447E-14
c 226,324 1.145E-29 385.579 52.306 4.923E-03 113,000
b 113,220 4.769E-14
b 226,441 5.181E-29 395.212 43.895 1.162E-02
c 113,220 2.403E-14
c 226,441 1.104E-29 385.583 52.281 4.925E-03 113,000

Table 2-5. Third iteration of PEIE for example #1. Progress from each iteration is
shown in Table 2-6.

Process Ea G Te te k Ea
Pair Guesses Estimate
(J/mol) (K) (mnin) (1/mmin) (J/mol)
a 113,000 8.650E-14
a 226,000 2.272E-28 404.817 32.924 2.627E-02
b 113,000 5.103E-14
b 226,000 1.266E-29 395.197 43.974 1.160E-02 113,000
a 113,000 8.650E-14
a 226,000 2.273E-28 404.817 32.924 2.627E-02
c 113,000 2.575E-14
c 226,000 1.266E-29 385.569 52.375 4.916E-03 113,000
b 113,000 5.103E-14
b 226,000 5.921E-29 395.197 43.974 1.160E-02
c 113,000 2.575E-14
c 226,000 1.266E-29 385.569 52.375 4.916E-03 113,000

Table 2-6. 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.790E-02 4.327E-08 1.932E-14

Table 2-7. First iteration of PEIE for example #2 (Ea2 = 2 x Eal 1). Progress from
each iteration is shown in Table 2-10.

Process Ea G Te te k Ea
Pair Guesses Estimate
(J/mol) (K) (min) (1/min) (J/mol)
a 300,000 7.541E-38
a 600,000 4.354E-76 409.801 13.062 6.797E-02
b 300,000 1.151E-38
b 600,000 9.489E-80 399.952 17.426 3.044E-02 111,153
a 300,000 7.541E-38
a 600,000 4.354E-76 409.801 13.062 6.797E-02
c 300,000 1.403E-39
c 600,000 9.489E-80 390.100 20.731 1.121E-02 121,607
b 300,000 1.151E-38
b 600,000 7.601E-78 399.952 17.426 3.044E-02
c 300,000 1.403E-39
c 600,000 9.489E-80 390.100 20.731 1.121E-02 131,555

Table 2-8. Second iteration of PEIE for example #2 (Ea2 = 2 x Eal). Progress from
each iteration is shown in Table 2-10.

Process Ea G Te te k Ea
Pair Guesses Estimate
(J/mol) (K) (min) (1/min) (J/mol)
a 111,153 1.504E-13
a 222,306 6.771E-28 404.689 33.426 2.656E-02
b 111,153 8.994E-14
b 222,306 6.073E-32 395.075 44.646 1.188E-02 111,242
a 121,607 6.580E-15
a 243,214 1.407E-30 405.365 30.768 2.886E-02
c 121,607 1.724E-15
c 243,214 6.073E-32 386.068 48.934 4.748E-03 121,681
b 131,555 1.738E-16
b 263,110 7.910E-34 396.244 38.195 1.389E-02
c 131,555 7.617E-17
c 263,110 1.276E-34 386.567 45.481 5.109E-03 131,617

Table 2-9. Third iteration of PEIE for example #2 (E,2 = 2 x Eal). Progress from
each iteration is shown in Table 2-10.

Process Ea G Te te k Ea
Pair Guesses Estimate
(J/mol) (K) (min) (1/min) (J/mol)
a 111,242 1.465E-13
a 222,484 6.423E-28 404.695 33.402 2.658E-02
b 111,242 8.751E-14
b 222,484 5.800E-32 395.081 44.613 1.189E-02 111,242
a 121,681 6.436E-15
a 243,3-62 1.347E-30 405.369 30.751 2.887E-02
c 121,681 1.684E-15
c 243,362 5.800E-32 386.072 48.906 4.751E-03 121,681
b 131,617 1.705E-16
b 263,235 7.617E-34 396.247 38.178 1.389E-02
c 131,617 7.470E-17
c 263,235 1.227E-34 386.570 45.461 5.111E-03 131,617

Table 2-10. 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.308E-02 3.094E-04 3.092E-04

Table 2-11. 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.520E-03 4.763E-04 2.982E-04 3.161E-04

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 2-1.

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 2-2.

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.





2.50 2.55 2.60 2.65 2.70
Inverse Absolute Temperature (x 1000)

Figure 2-3.



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 extent-of-
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 2-4.


Ea2 Scheme
o 10,000
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 -




-4.50 -

-5.00 -

-5.50 -

-6.00 -

-6.50 -

-7.00 -

Figure 2-5.


-2.50 -









-7.00 -

Figure 2-6.

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.



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 time-temperature data into lethality (extent-

of-reaction) 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 two-dimensional 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 two-dimensional

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 non-capacitance 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 extent-of-reaction 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 non-capacitance 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

time-temperature data for approximating extents of reactions (lethality) in systems

whose thermal histories are estimated with numerical heat transfer methods.



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


. 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 heat-transfer 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

half-height, distances, dr and dy, separate nodes along the radius and half-height,

respectively. These distances are calculated as follows (see Fig. 3-1):

dr = Radius / (Nr- 1) (3-1a)

dy = HalfHeight / (Ny 1) (3-1 b)

The situation is slightly different for the NCSN case (Fig. 3-2). 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 2-Nr 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. 3-2 and Fig. 3-4)

dr 2 Radius / (2- Nr 3) (3-2a)

dy = 2 HalfHleight / (2- Ny 3) (3-2b)

Visual Basic code for discretizing a finite cylinder using the NCSN is shown in Fig. 3-3.

Note that several lines of code (Fig. 3-3) 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 cross-sections of equal area, from the bottom to the top of the cylinder. Heat

flowing along the radial direction, however, passes through non-constant areas. Thus,

tops and bottoms of volume elements are placed at the mean distance between nodes.

Sides of volume elements are placed at the log-mean radius between nodes (Fig. 3-3):

rim = (ro ri) / ln(r o/ri) (3-3)

where ro, ri, and rim are outer radius, inner radius and log mean radius, respectively.

Since volume elements containing the vertical center-line 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. 3-3).

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. 3-3). 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 heat-transfer situations. Nine types of volume elements were identified

for both the CSN and NCSN nodal arrangements for a finite cylinder (labeled A I,

Fig. 3-1 and Fig. 3-2 and Table 3-1 and Table 3-2). 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. 3-1 and Table 3-2) and the

NCSN approach (Fig. 3-2 and Table 3-2), 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. 3-1

and in the diagram for node type C in Table 3-1, 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 (3-4)
+h. rAi. j-(T-Tr.)

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), TI-Lj 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. 3-4 for Ti'.1 yields

Tj+l' = Tj (1 -cL- cT cR) + cL T-I.j + cT. Tj + cR T (3-5)

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 3-1). The coefficient, cB, was used for volume elements with heat flow from

the bottom. The ratio k / (p-cp) is thermal diffusivity, denoted by the Greek letter, ta.

In order for Eq. 3-5 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 -< (3-6)
1 -k-rA -i. + k- yA...j' + h. rAi.(
p- Cp dVi. j dr dy )

Eq. 3-6 shows that as the heat transfer coefficient, h, increases, the time step, dt, must


Non-Capacitance Surface Nodes (NCSN)

In this case, the node within volume element type C is at a distance dr/2 from

the surface (Fig. 3-2 and Table 3-2). 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) k-r Ti-i.j Ti.j) + Tij+i Tinj

T Tqj (3-7)
+ 1 dr/2
h .rAi+i,j k- rAi.j

Solving Eq. 3-7 for T.t1 yields

Tj = Ti.j" (1 cL cT cR) + cL. TP-i.j + cT. Ti.j+ + cR T(


Eq. 3-8 is similar in form to Eq. 3-5, but the coefficient, cR, includes terms related to

both convection and conduction. Similarly, the sum of the terms in parentheses must

be non-negative. Thus, the maximum time step for the NCSN case is given by

d NCNS (39)
tNode C -
1 k rAi-i.j +hk r yAu.j +k 1+I
p- cp dVi. j dr dy 1 dr/2
h-rAi+ij k-rAi.j /

Aside from describing the physical situation more appropriately, the superiority of the

NCSN approach is clearly demonstrated by Eq. 3-9. 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). (3-10)

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. 3-10 and Eq. 3-11

for, T4tji and T., respectively. Notice that the time step, dt, does not appear in Eq.

3-10 or Eq. 3-11. This demonstrates that surface temperatures are not required to

specify internal temperatures, however, surface temperatures may be computed as


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 3-1 and 3-2, 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 3-4 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 3-3. 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 3-3 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. 3-2) 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. 3-7 would take the form

T'pn V-pn
T, -T F (3-12)
1 dr,, dr / 2 (3-12)
--+ +
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 3-3,

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. 3-5). 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 3-5 demonstrates that both methods offered comparable accuracy.

Additionally, Fig. 3-5 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 Time-Temperature 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 time-temperature 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. 3-6 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


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 extent-of-reaction 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} (3-13)

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. 3-13 is different

from thermal conductivity, k, in Eq. 3-4, Eq. 3-6, Eq. 3-7, Eq. 3-9 and Table 3-1 and

Table 3-2). 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 pre-exponential 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) (3-15)
C 0R lR T(t)J -

where T, is the equivalent isothermal temperature that is represented by option e in Fig.

3-6. For a small time step, (tn- tn), the linear temperature change approximation is

T.+l T!j+ [slope] -t[;n- (3-16)
i~j --- i~jtil

The slope in Eq. 3-16 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


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 3-7 and 3-8 show how the location of the equivalent

isothermal temperature varies according to these parameters. Figure 3-7 is for a high

activation energy constituent (500 kJ/mol), such as a bacterial spore population. Figure

3-8 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. 3-7 and Fig. 3-8 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. (3-17)
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. 3-7 and Fig. 3-8 show that for extent-of-reaction

predictions it is appropriate to select temperatures at or within 10% above the mean,

for each time-interval.

Table 3-1. CSN finite difference formulas (Fig. 3-1).

Type Diagram Node(s) Coefficients
A A i =1 cR = a-cAt / dV(i, j)- rA(i,j) / dr
j = 1 cT = a-At / dV(i, j) yA(i) / d,

T -' = Tj" (1 cR cT) + cR Ti.j + cT -Tl^,_
B ----i i=2..Nr-1 cL = a-At / dV(i, j) rA(i-1, j) / dr
j = 1 cR = ca-At / dV(i, j)- rA(i, j) / dr
cT = ca-At / dV(i, j) yA(i) / dy
T. = T.j" (1 cL cR cT) + cL T- -Ij + cR TI.j + cT
C -- i=Nr cL = a-At / dV(i, j) rA(i-1,j)/dr
j = 1 cT = a-At / dV(i, j) vA(i) / dv
/4 cR = a.At / dV(i, j) h-rA(i. j)/k
._= T (1 cL cT- cR) + cL- C L jT_ + cT. T',j + cR.T _
D i = 1 cR = a-At / dV(i, j). rA(i, j) / dr
j2.. Nv- cT = a-At / dV(i,j) -yA(i) / dy
cB = ca-At / 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.. Nr-1 cL = ca-At / dV(i,j) rA(i-l j) / dr
j = 2 ..Ny-i cR = -At / dV(i,j) rA(i, j) / dr
cT = a-At / dV(i, j). yA(i) / dv
cB = a-At / 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 = a-At / dV(i,j) rA(i-lI j) /dr
j = 22.. Nv-I T = cT=a-At / dV(i,j)- yA(i)/dv
cB = a-At / dV(i, j)- yA(i) / dy
.c5\ R = a-At / dV(i, j) h. rA(i. j) / k
T = T'r (1 cL cR cT cB) + cL. T,-1.j + cT* T,,. + cB T.j-i + cR* T
G i = l cR = ca-At / dV(i, j) rA(i, j) / dr
j = Ny cB = a-At / dV(i, j) yA(i) / dv
cT = a-At / dV(i, j) h yA(i) / k
Tf;1 = T (1 cR cB cT) +cR T + cB "T,. + +cT- T+c
H i=2.. Nr-1 cL = a-.At / dV(i,j) -rA(i-1, j) /dr
j J= Ny cR = a.At /dV(i, j) rA(i, j) / dr
cB = ca-At / dV(i, j) yA(i) / dv
cT = a-At /dV(i, j) -h- yA(i) /k
T ,-(1 & cL-cR- cB- +cL-+cR T-..j + cB.l + cT. CB
I i = Nr cL = a-At / dV(i, j) rA(i-1, j) / dr
j1 1 = Ny cB = ca-At / dV(i, j) yA(i) / d"
R cR =a-At / dV(i, j). h rA(i, j) / k
cT = a-.At / dV(i, j) h yA(i) / k
Tl-' = T (I- cL-cB- cR- cT) + cL'-V, +cB-T- +cR I +cT'T

Table 3-2. NCSN finite difference formulas (Fig. 3-2).

Type Diagram Node(s) Coefficients
A i= cR= a.-At / dV(i,j) rA(i,j) / dr
j = cT = a-At / dV(i, j)- yA(i) / d-

-V' j -* (1 cR cT) + cR Tlj + cT j-
B i=2..Nr-1 cL = a-At / dV(i, j) rA(i-1, j) / dr
j = I J lcR=a- At / dV(i, j).rA(i,j) / dr
cT = a-At / 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 = a-At / dV(i, j) rA(i-1, j) / dr
j = 1 cT = a-At / dV(i, j) yA(i) / dv
ScR = aAt / dV(i j)- h-rA(i. j) /k
T;= j (1 cL cT- cR) + cL- T-, + cT. T-,i- + cR Tr
D / i = 1 cR = ca-At / dV(i, j) rA(i, j) / dr
j = 2 ,. Ny-1 cT = a-At / 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..Nr-1 cL = ca.At / dV(i, j) rA(i-1, j) / dr
j = 22.. Ny-l cR = a-At / dV(i, j). rA(i,j) / dr
cT = a-At / 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 T-Ij + cR T 1 + cT TjI + cB. ji-
F -- i=Nr cL = a-At / dV(i,j). rA(i-L j) / dr
Sj=2.. Ny-i cT = a-At / dV(i,j). yA(i) / dv
cB = a-At / dV(i, j). yA(i) / dyv
cR = ca-At / 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 = a-At / dV(i, j). rA(i, j) / dr
j = Ny cB = a-At / dV(i,j) yA(i) / dy
V7 cT = a-At / 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..Nr-1 cL = ca-At / dV(i, j) rA(i-l, j) / dr
j j= Ny cR = at.At / dV(i, j)- rA(i, j) / dr
ScB = a-At / dV(i, j) yA(i) / dv
cT = a-At / dV(i, j) h vA(i) / k
T". = T ," (1 cL cR cB cT) + cL TT-L,+ + cR -Ti,. + cB T,.i+ cT T2
i i=Nr cL = a-At / dV(i, j) rA(i-1, j) / dr
| j = Ny cB = a-At / dV(i, j) yA(i) / dv
cR = a-At / 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 3-3. 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
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

E E a Ezi i-
. .^, F

D E i> j .
PD 0 "



- /

.> /i

-r /

7 '^







<.R a d i u s

Spatial descritization of a finite cylinder for CSN.


Figure 3-1.




n0 n ii yA(i)
7 r-- rA(i-1,j)
E l rA(i,j)
0) / D E"I, # F

SE D .-dr--.
- / E E E F

" -, r -


<---- Rad i u s >

Spatial descritization of a finite cylinder for NCSN.

Figure 3-2.

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)]

r(1) =0: r(Nr) = Radius
For i = 2 To Nr 1
r(i) = r(i- 1) + dr
y(1 ) = 0: y(Ny) = HalfHeight
For i = 2 To Ny 1
y(i) = y(i 1) + dy
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)))
yALocation(1l) = y(2) / 2: yALocation(Nv) = HalfHeight
Forj=2ToNy I
yALocation(j) = yALocation(i 1) + dy
For j = 1 To Ny
For i = 1 To Nr
rA(i, j) = 2 Pi rALocation(i) (yALocation(j) yALocation(i 1))
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
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 j

Figure 3-3. Visual Basic code for discretizing a finite cylinder using NCSN. Code
in square brackets may be substituted to yield CSN.


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

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)
Next j

Figure 3-4. Visual Basic code for main loop of explicit finite difference heat transfer
simulation using NCSN. Code in square brackets is for CSN.

110 20 minutes

60 minutes
50 0 Capacitance Surface Nodes
40 13 Non-capacitance Surface Nodes
30 Analytical Solution
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 3-5.

Temperature distributions along the radius on the mid-plane of a finite
cylinder at 3, 10 and 20 minutes. NCSN offers comparable accuracy to
CSN, but with significant time savings.


Possible temperature selections for calculating extent of reactions
(lethality) occurring during each time interval, in each volume element.

T -



-I. I

Figure 3-6



0 0.8-

E 0.7




Figure 3-7.

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.

50 KJ/mol


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 --
1 second

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 3-8. 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.



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 (4-1)
k = k. -exp Tabsolute I

where the kinetic parameters consist of the pre-exponential 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/(mol-Kelvin).

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


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 non-thermal,

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 non-thermal,

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 non-thermal 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 non-isothermal 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 heat-up and cool-down, (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


The non-isothermal approach is generally more complicated than the isothermal

approach. However, the non-isothermal 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, non-isothermal parameter estimation methods

overcome the problem of thermal lags, because the entire exposure is considered to

affect the result. Non-isothermal 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 non-isothermal

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 wide-spread use

of non-linear optimization methods.

A non-isothermal 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 non-isothermal 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.


The following applies to the case ofreactants that follow apparent first order

kinetics with Arrhenius temperature dependence. A first order rate process is defined


SC k-C (4-2)

where, C is the concentration of a reactant, t is time, and k is the rate constant as

defined by Eq. 4-1. Integration of Eq. 4-2 yields

ln(C/Co)=-k.(t- to) (4-3)

Thus, extent of a given reaction under isothermal conditions may be determined by

substituting Eq. 4-1 into Eq. 4-3:

lnC/0)~k.eP Ea }(t -t.) (4-4)
ln(C / Cj = -k, exp R. ----o. (t- t) (4-4)

Under non-isothermal conditions, Eq. 4-2 must be integrated over the entire thermal

history, T(t):

ln(C / C.) = -ko. J exp,-Ea dt (4-5)

For brevity, the integral in Eq. 4-5 is referred to as G (Chapter 2):

G = Jexp--- dt (4-6)
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. 4-6:

G = te .exp{-E} (4-7)

Here te and Te represent the equivalent time and isothermal temperature, respectively,

that provide the same extent of reaction as the non-isothermal process, for reactants

characterized by the particular value of Ea. Since the value of G is calculated from Eq.

4-6 using a known thermal history and a known or assumed E, value, Eq. 4-7, 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. 4-7:

ln(t,) = ln(G) + E (4-8)
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. 4-8, for a given non-isothermal 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 non-isothermal 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. 4-1. The figure depicts two iterations

of the method for two non-isothermal 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 (Ea-l a and Ea-lb). For each dynamic thermal

exposure considered, the lower Ea value, Ea-lb, yields curves with shallower slopes,

while the higher value, Ea-l 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 extent-of-reaction data. New Ea guesses,

Ea-2a and Ea-2b, are obtained by replacing Ea-la with the calculated Ea estimate and

setting Ea2b to either a fixed value or some multiple of Ea-2a. 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 extent-of-reaction 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 non-isothermal experiments, n-(n-1)/2, potential

sets of kinetic parameters may be defined. For example, for 4 non-isothermal

experiments (a, b, c, d), a maximum of 6 sets of kinetic parameters may be expected

(one from each pair: a-b, a-c, a-d, b-c, b-d, c-d).

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 in-depth 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. 4-2). 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 non-isothermal 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. 4-4 and Eq. 4-5 yields

_._j = exp{-k.0 G i} (4-9)

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} (4-10)

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-)

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 pre-exponential

factor, k0, is the only unknown in Eq. 4-11. 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. 4-1). 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. 4-1, and then subtracting the resulting equation from Eq. 4-1. The unwieldy ko

value cancels out and the resulting equation takes the form

k = kr exp -Ea FF,- T (4-12)

where kr is the reaction rate constant at temperature, Tr. Thus, k1 may be replaced with

kr in Eq. 4-9 through Eq. 4-11. 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 non-uniformly, involves the following steps, which are illustrated in Fig. 4-3:

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. 4-3).

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. 4-2). During simulation, apply the Ea values to

accumulate G values according to Eq. 4-7. 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. 4-11.

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 (Eai-kri and Ea2-kr2).

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. 4-7)

simultaneously. Using Eq. 4-3, determine respective equivalent isothermal rate

constants, kT. Since kinetic parameters, Eai-kri, provide (Ci,j / Co)1 and

parameters, Ea2-kr2, 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.


Numerical Heat Transfer

The explicit finite-difference 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 4-1. Sources and methods for

obtaining these data are given below in the section related to actual experiments.


Simulated experiments

Simulated experiments were performed to demonstrate the use of the PEIE

method for the case of non-uniform heating. A food constituent was considered to be

characterized by the following Arrhenius kinetic parameters: k121.1ic= 0.2 min-1, 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 4-2). During retort

come-up and heating phases, heat transfer coefficient, hHeat (Table 4-1), was used.

During retort cool-down, hCool (Table 4-1), was used. When the retort temperature

reached 30C, the can was considered to have been placed into an ice bath (hIce, Table

4-1). 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 4-2).

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


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 freeze-dried 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 head-space 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 4-2). 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 (4-13)
Ti-T. 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 4-13 shows that the heat transfer

coefficient, h, may be inferred from the slope of a semi-log plot of accomplished

temperature versus time.

Cans were exposed in duplicate to different retort processes. In general, three

approximate set-point 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. 4-4 and Fig. 4-5.

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. 4-1) 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 mercury-in-glass (MIG) thermometer installed on

the retort. The pilot scale vertical still-cook 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 non-selective 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 4-2). 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 4-3 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 4-2, extent-of-reaction 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. 4-6 and Fig. 4-7. Since 14 experiments were

considered, a maximum of 91 experimental pairs, and thus, 91 sets of parameters were

possible. However, Fig. 4-6 and Fig. 4-7 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.

4-6 and Fig. 4-7.

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. 4-4), 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

set-point temperature were analyzed only with processes associated with different

retort set-point temperatures. For example, process 1 (Fig. 4-4) was not paired with

processes 2 through 8, but was paired with processes 9 through 14 (Fig. 4-5). 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 log-mean 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


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. 4-8), 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. 4-8, 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 first-order

inactivation. Additionally, the PEIE method operates by connecting two points on a

survivor curve (Fig. 4-8). 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 log-cycle reduction) at 121.1 C, and (2) time required for the straight line portion

to traverse four log-cycles 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

log-cycle basis are the definition of Do values (time required to achieve a one log-cycle

reduction in concentration. The Do value is related to the reaction rate constant,

k121.1oc, via

Do = 0 (4-14)

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).