NUMERICAL MODELING OF NATURAL CONVECTION
AND CONDUCTION HEAT TRANSFER IN CANNED FOODS
WITH APPLICATION TO ONLINE PROCESS CONTROL
By
ASHIM KUMAR DATTA
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1985
Dedicated to the memory of my
grandfather and father
ACKNOWLEDGEMENTS
The author is indebted to Dr. Arthur A. Teixeira, his major
professor, for his encouragement, guidance and patience throughout
the course of the work and especially for his assistance in
preparing the manuscript.
The author is also highly grateful to Dr. Tom IP. Shih, Dr.
Chung K. Hsieh and Dr. Ruey J. Shyu for their time and extremely
valuable suggestions in the study of the natural convection heat
transfer. Further thanks to Dr. Khe V. Chau and Dr. Spyros Svoronos
for their interest in the study and helpful suggestions.
The importance of a good computing facility (number crunching
to be particular) in developing the natural convection model cannot
be overemphasized. The availability of the Cyber 170 computer
through Dr. John Gerber and the IFAS administration was vital for
that purpose. The author also greatly appreciates help from Dr.
Ellen Chen, who introduced him to the Cyber facility, and Victor,
Jimmy and Margaret for endless support on problems with the Cyber
system.
The mountain of data available from the transient natural
convection heating model would have been of little use without the
plotting capability provided by the excellent graphics software
(PLOTPAK) made available from Dr. William E. Dunn from the
Mechanical and Industrial Engineering Department of the University
of Illinois at UrbanaChampaign. Help by Subhasis Laha in physically
transporting the plotting software to Gainesville is also
appreciated. Availability of a hard copy plotter was the next
necessity and the author appreciates the equipment provided by Dr.
James W. Jones in the Agricultural Engineering Department.
Finally, the author deeply acknowledges the assistance from his
wife Anasua in typing this manuscript, and her encouragement and
perseverence during the study.
TABLE OF CONTENTS
ACKNOWLEDGMENTS .................... ...........................i
LIST OF SYMBOLS................................................ vi
LIST OF TABLES....................................................xi
LIST OF FIGURES............ ................................ xii
ABSTRACT ........................................ .......... xvi
INTRODUCTION...................... .............................
REVIEW OF LITERATURE.............................................5
Computer Control of Food Thermal Processing..................5
Overview of Food Thermal Processing........................5
Online Control of Thermal Processing....................10
Conduction Heating of Canned Foods.........................12
Studies Based on Analytical Solution.....................12
Numerical Studies .......................................13
Natural Convection Heating in Canned Foods..................14
Previous Works on Canned Liquid Foods..................15
Other Similiar Works on Natural Convection
Heating in Cylindrical Enclosures........................19
PROBLEM FORMULATION...........................................24
Natural Convection Heat Transfer in Canned Foods............24
Governing Equations and Boundary Conditions..............24
Boussinesq Approximation.................................26
Transformation of the PDEs.............................27
Nondimensionalization of the Variables.................28
Conduction Heat Transfer in Canned Foods...................30
Governing Equation..................................31
Boundary Conditions.............................. .....31
Computer Control of Conduction Heating Process..............32
METHODOLOGY .................................................... 35
Modeling of Natural Convection Heating in Cans..............35
Grid Generation and Grid Stretching......................36
Discretization of the Parabolic
Temperature and Vorticity Equation......................39
Time derivatives .............................. .....42
Linearization. ........................................44
Convection terms......................................45
Diffusion terms ...................................... 50
Discretization of the Elliptic
Stream Function Equation..................................61
Computational Boundary Conditions for the FDEs...........69
Convergence Criteria for Iteration .......................74
Algorithm for Iterative Solution.........................76
Coding of the Computer Program and Hardware..............80
Modeling of Conduction Heating in Cans......................80
Discretization and Solution of the Equation..............81
Coding of the Computer Program..........................82
Experimental Studies ....... ........ ....................83
Use of the Conduction Heating Model for Online Control.....86
RESULTS AND DISCUSSION.... ......................... ...... ......93
Natural Convection Heating in Canned Foods..................93
Transient Flow Patterns and Temperature Profiles.........94
Start of flow and conduction layer....................95
Radial temperature and velocity profiles.............106
Axial temperature profiles........................... 113
Slowest Heating Zone................ ..... ............. 115
Contrast with Conduction................................120
Assessment of the Numerical Method......................124
Convergence..........................................124
Selection of the grid size...........................125
Selection of the time step...........................130
Convergence of boundary vorticity...................134
A note on the half grid points.......................134
Conduction Heating in Cans.................................138
Analysis of Boundary Conditions........................139
Comparison of the Transient Temperature Values.......... 141
Effect of Sudden Pressure Drops in the Can.............. 141
Performance of the Online Control Logic...................141
CONCLUSIONS AND RECOMMENDATIONS..............................151
Conclusions................................................151
Recommendations............................................ 153
APPENDIX
A ALTERNATIVE FINITE DIFFERENCING OF
STREAM FUNCTION EQUATION................................154
B ALTERNATIVE FINITE DIFFERENCING OF BOUNDARY VORTICITY...155
C INPUT DATA FILE FOR NATURAL CONVECTION MODEL............156
D ADE FORMULATION FOR THE CONDUCTION EQUATION.............157
LIST OF REFERENCES.......... ..... .............. ................ 161
ADDITIONAL REFERENCES............................. ....... ....... 167
BIOGRAPHICAL SKETCH................................... ....... 172
LIST OF SYMBOLS
A = constant
a,b,c = coefficients in tridiagonal matrices
f = inverse slope of timetemperature graph
g = gravity
Fo = sterilization Fo value
H3
Gr = Grashof number = gB (Ti To) 7
h = outside heat transfer coefficient
H = height of cylinder
j = j value (a constant)
JO = Bessel function of order zero
k = thermal conductivity
n = normal direction
nz = number of grid points in z direction
nr = number of grid points in r direction
p = pressure
Pr = Prandtl number = pC /k = v/a
r = non dimensional distance in radial direction
r = distance in radial direction
R = radius of cylinder
= Rayleigh number = Gr*Pr
= time
= temperature at any point inside at time t > 0
= initial uniform temperature inside (at t = 0)
= boundary temperature at time t > 0
= retort temperature at time t
H 
= u = nondimensional velocity in vertical direction
= velocity in vertical direction
H
= v = nondimensional velocity in radial direction
= velocity in radial direction
 nondimensional distance in vertical direction
H
= distance in vertical direction
= Zvalue for organism
= thermal diffusivity = k/pCp
= coefficient of thermal expansion
= relaxation parameter for boundary vorticity
= designates a difference when used as a prefix
= convergence criterion for boundary vorticity
= convergence criterion for SOR
= vertical distance in computational plane = n(z)
TTi
= nondimensonal temperature
ToT
= first eigenvalue in z direction
= eigenvalues in z direction
= first eigenvalue in r direction
= eigenvalues in r direction
viii
= viscosity of liquid being heated
= kinematic viscosity of liquid being heated = p/p
= radial distance in computational plane = (r)
= density of liquid being heated
= deformation parameter in grid stretching
= t = nondimensional time
H2
= radial coefficients for kinematic consistency
= radial coefficients for kinematic consistency
= nondimensional stream function = 
Ha
= stream function defined by v = 
r ar r az
= nondimensional vorticity
av au H2 v U H2 
az ar a a
= vorticity
= optimum relaxation parameter for SOR
Subscripts
a
center
cool
heat
i
J
m
n
W
= ambient
= geometrical center of can
= cooling cycle
= heating cycle
= grid points in vertical direction
= grid points in radial direction
= grid point on vertical boundary (z = H)
= grid point on radial boundary (F = R)
= Wall
Superscripts
d = design value
f = normal to grid faces
k = iteration counts in SOR
m = iteration count within a time step
n = time step
t = tangential to grid faces
= "half" step in ADI method (same as n+l/2)
LIST OF TABLES
TABLE PAGE
1. Sequence of line types for identification of stream
function contours in figures 10ao........................96
2. Heating and cooling rates for various process
conditions........................... ..................... 140
3. Adjusted heating times and resulting lethality (Fo)
in response to process deviations using proposed on
line control logic and method of GiannoniSuccar and
Hayakawa ................................................. 148
LIST OF FIGURES
FIGURE PAGE
1. Isotherms and velocity profiles in a glass bottle at
various heating times (after Engelman and Sani, 1983).....18
2. Isotherms and velocity profiles in a can at various
heating times (after Sani, 1985).........................20
3. The grid system and the boundary conditions for the
cylindrical can...........................................40
4. Notations of various quantities defined on a non
uniform grid system in cylindrical geometry...............62
5. Algorithm for the iterative solution of the set of
equations in natural convection...........................77
6. Experimental setup for conduction heating..................84
7. Timings of various realtime computations for online
control of conductionheated food.........................88
8. Flow diagram for computer control of retort opera
tions with online correction of process deviations.......90
9. Isotherms near bottom wall predicted by the
convective model compared with the predictions
considering conduction (equation 83) only.................97
10a. Isotherm and streamlines in a cylindrical can after 1
second of heating .........................................99
10b. Isotherm and streamlines in a cylindrical can after 2
seconds of heating........ .......... ......................99
10c. Isotherm and streamlines in a cylindrical can after 3
seconds of heating.......................................100
lOd. Isotherm and streamlines in a cylindrical can after 4
seconds of heating......................... ......... ..... 100
10e. Isotherm and streamlines in a cylindrical can after 5
seconds of heating.................. ... .. ............ ..101
10f. Isotherm and streamlines in a cylindrical can after 6
seconds of heating....................................... 101
10g. Isotherm and streamlines in a cylindrical can after 7
seconds of heating......................................102
10h. Isotherm and streamlines in a cylindrical can after 8
seconds of heating...................................... 102
10i. Isotherm and streamlines in a cylindrical can after 9
seconds of heating.................... ... ....... ......... 103
10j. Isotherm and streamlines in a cylindrical can after
10 seconds of heating................. ........... ...... 103
10k. Isotherms and streamlines in a cylindrical can after
30 seconds of heating................... ................ 104
101. Isotherms and streamlines in a cylindrical can after
120 seconds of heating.......................... ....... 104
10m. Isotherm and streamlines in a cylindrical can after
300 seconds of heating.............. ..................... 105
1On. Isotherms and streamlines in a cylindrical can after
600 seconds of heating................................... 105
10o. Isotherms and velocity vectors in a cylindrical can
after 1800 seconds of heating...........................107
11. Predicted radial velocity profile at the sidewall at
mid height in a cylindrical can after 30 seconds of
heating ..................................................108
12. Predicted radial temperature profile at the sidewall
at mid height in a cylindrical can after 30 seconds
of heating ...............................................108
13. Isothermal vertical surface, showing boundary layer
with velocity profile u(x,y) and temperature profile
t(x,y) in water (Higgins and Gebhart, 1983)..............110
14. Predicted radial velocity profiles at the sidewall at
mid height in a cylindrical can after various heating
times .......................... ........................... 111
15. Observed radial velocity profiles at the sidewall at
mid height in water in a cylindrical can after
various heating times (from Hiddink, 1975)...............112
16. Comparison of numerically predicted axial temperature
profiles with observed (Hiddink, 1975) values at
various times during heating.............................114
xiii
17. Locations of the slowest heating points in water over
a heating period of 10 minutes in a cylindrical can......116
18. Migration of the slowest heating points in water over
a heating period of 10 minutes in a cylindrical can......116
19. Transient temperature at various locations during
heating of water in a cylindrical can....................119
20. Transient temperature at various locations during
initial times of heating of water in a cylindrical
can.......... ... ..... ..................... .......** 121
21. Temperature at the slowest heating point during
natural convection heating contrasted with that
during conduction heating in a cylindrical can...........122
22. Isotherms during natural convection heating contr
asted with those during conduction heating in a
cylindrical can..........................................123
23. Diagram of 39(0.9)x39(0.8) grid on the left and
58(0.85)x58(0.76) on the right...........................126
24. Isotherms after 90 seconds of heating computed using
39x39 grid (left half of the figure) compared with
computations using 58x58 grid (right half of the
figure) ..................................................127
25. Streamlines after 90 seconds of heating computed
using 39x39 grid (left half of the figure) compared
with computations using 58x58 grid (right half of the
figure) .................................................. 128
26. Transient slowest heating temperatures computed with
39x39 and 58x58 grid.................. .................... 129
27. Transient temperature at 1/3rd height on the center
line computed with 39x39 and 58x58 grid...................131
28. Time step history during a typical run of the
numerical heat transfer model for natural convection.....132
29. Isotherms in a cylindrical can calculated with two
different convergence criteria for boundary vorticity
(6 =0.0001 on the left and 6 =0.001 on the right)........135
30. Streamlines in a cylindrical can calculated with two
different convergence criteria for boundary vorticity
(6 =0.0001 on the left and 6 =0.001 on the right)........136
31. Transient center temperature during conduction
heating of Bentonite in a cylindrical can................142
32. Transient center temperature during conduction
cooling of Bentonite in a cylindrical can with and
without drop in retort pressure...........................143
33. Computer plot of reference process (no deviations)
showing 66.8 minutes of heating time and an
accomplished Fo of 6.24............ .............. ......144
34. Computer plot of a process that experienced a step
functional drop of 8.30C (150F) in temperature bet
ween 50th and 60th minute and still maintained an Fo
of 6.02.................................................. 146
35. Computer plot of a process that experienced a linear
drop in temperature between 40th and 60th minute
(from 00C to 11.1C (200F)) and still maintained an
Fo of 6.18..............................................147
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
NUMERICAL MODELING OF NATURAL CONVECTION
AND CONDUCTION HEAT TRANSFER IN CANNED FOODS
WITH APPLICATION TO ONLINE PROCESS CONTROL
By
ASHIM KUMAR DATTA
December 1985
Chairman: Arthur A. Teixeira
Major Department: Agricultural Engineering
A numerical model was developed for predicting detailed flow
patterns and temperature profiles during natural convection heating
of liquids in a can (a cylindrical enclosure). The liquid (water)
was initially stagnant at a uniform temperature of 300C, and the
sidewall, top and bottom wall temperatures were suddenly raised to a
temperature of 1210C. A free (thermally insulated) top liquid
surface was also considered. Boussinesq approximation was used and
all other fluid properties were treated as constant.
Finite difference methods were used to solve the governing
equations in axisymmetric cylindrical coordinates. A vorticity
stream function formulation was used. Parabolic temperature and
vorticity equations were solved with the alternating direction
implicit method and upwinding the convective terms. The elliptic
stream function equation was solved using successive over relaxation
method. Grid stretching was used in both radial and axial
directions. Boundary wall vorticities were used to check for
convergence of the iteration process.
Plots of transient isotherms, streamlines and velocities were
provided. Calculated radial velocities and axial temperature
profiles at different times compared quite well with available
experimental data. Slowest heating points were located in the bottom
15% of the height of the container studied. These slowest heating
points migrated over time within this lower region of the can but no
particular pattern of migration was noted.
Transient can center temperature was also studied for the case
of conductionheated solid food. An alternating direction explicit
finite difference solution to the conduction equation compared very
well with analytical solutions and experimental results. The ability
of the finite difference model to perform under arbitrary variation
of boundary temperature allowed the model to be incorporated in an
algorithm for online control in batch thermal processing of
conductionheated food. The uniqueness of this control system in
maintaining the level of desired food sterilization for arbitrary
variation in process heating temperature was shown to offer
considerable advantages over other possible methods of online
control.
xvii
INTRODUCTION
Background and Justification
Thermal preservation of food involves application of heat
sufficient to destroy the microorganisms present in food that cause
spoilage. Commercial canning processes provide such heat treatment
to extend the shelflife of foods. In a typical canning process, the
container is filled with food and sealed. The sealed container is
heated in a closed vessel with steam or hot water long enough to
kill the microorganisms and is then cooled. In a production situa
tion, the heating medium (steam or hot water) temperature can
sometimes deviate significantly from design values during the
process. Heating the product more than the required amount would un
necessarily degrade its quality while wasting energy. However,
insufficient heating could lead to serious public health hazards and
cannot be allowed. Thus, food processors always have to meet this
minimum heat treatment (sterilization). Determination of the heating
time to meet the minimum sterilization demands that the actual
transient temperature history of the food during processing be known
predicted with a heat transfer model or actually measured).
ConductionHeated Foods
Traditionally, for conductionheated foods, the required
processing time is determined using the analytical solution to the
heat conduction equation. In practice, the constant medium temper
ature assumed by the analytical solution cannot always be achieved.
Thus, experimental correction factors have been proposed to adjust
the heating time, predicted by analytical solution, when the medium
temperature goes through certain (predefined) types of deviations.
However, deviations in practice are arbitrary, and it is not
conceivable to know the correction factors for all such situations.
Of course, the problem is trivial if the transient food temperature
is directly measured using a thermocouple; but this is both
inconvenient and impractical in a production situation. Better
models for predicting food temperature in response to arbitrary
fluctuations in medium temperatures are therefore needed.
A numerical finite difference model would be able to predict
the food temperature for truly arbitrary variations in the medium
temperature. However, the calculations required are relatively much
more complicated and time consuming. Using a microcomputer to
perform this would have two advantages. In addition to correctly
predicting food temperature for truly arbitrary heating conditions,
it can perform this automatically without any worker supervision.
Introduction of lowcost microcomputers for online process control
is becoming commonplace in the 1980s. Food industries are not an
exception to this situation. Use of microcomputers for online
control of conductionheated food is thus quite conceivable.
ConvectionHeated Foods
Unlike solid packed foods that are heated by conduction, thin
liquid foods heat through convection. In general such liquid foods
are mechanically agitated to improve the heat transfer. However,
there are situations where the food material cannot be agitated for
several reasons. The food then heats primarily by natural convec
tion. Inpackage pasteurization for liquids (e.g. beer) in bottles
or cans that cannot be agitated is an example. Broth, thin soups,
evaporated milk, most fruit and vegetable juices, fruits in syrup or
water, and pureed vegetables are also examples of products where
natural convection heating can occur when no agitation is used.
Flexible pouches that cannot be agitated to keep the package integ
rity are another important example. Fermentation is another process
where a natural convection heat transfer model (with heat generation
term included) can provide valuable insight.
Just as optimum heating times are required for conduction
heated foods, they are also needed for natural convectionheated
products. Thus heat transfer models are similarly needed. Since we
are interested in the minimum heat treatment the food receives, we
need to know the point or the region that heats the slowest. Because
convection sets up complex flow patterns, the location where the
temperature is minimum is by no means obvious. Only a study of
complete transient flow patterns and temperature profiles can give
such insight. Previous studies to achieve this objective have not
been successful. The use of more advanced numerical techniques and
better computing equipment would make it possible to study in detail
the natural convection heat transfer in canned liquid foods.
Such a natural convection heat transfer model could also have
several applications outside the food area. It may be used to study
the heating of buildings (without forced circulation), cooling of
gas turbine blades, storage of cryogenic fluids (e.g. liquid rocket
propellants), and the startup of chemical reactors. Since natural
convection heat transfer would involve a much higher degree of math
ematical complexity than a conduction heating model (and its appli
cation to online control), the primary objective of this study was
directed toward developing a heat transfer model for natural convec
tion; followed by improvements to existing conductionheat transfer
models for application to online process control as secondary
objectives.
Objectives
The objectives of this study were to
1. develop a mathematical model to predict temperature and
velocity profiles caused by natural convection heat
transfer in a closed cylinder,
2. investigate the need for possible improvements on an
existing conduction heat transfer model to improve its
suitability for use in online process control appli
cations, and
3. incorporate the conduction heating model in an online
process control algorithm for thermal processing of
canned food, and compare its performance against other
possible methods of online control.
REVIEW OF LITERATURE
This chapter begins with a brief introduction to thermal
processing as applied to the sterilization of canned foods. The need
for computer control of thermal processing is noted, along with the
importance of heat transfer models in such computer control systems.
Previously reported conduction heat transfer models used in studying
sterilization of conductionheated canned foods are described next.
The few studies attempting to model natural convection heat transfer
in the case of sterilizing canned liquid foods are noted. The few
other works on natural convection heat transfer in liquid inside a
vertical cylindrical enclosure are also reported.
Computer Control of Food Thermal Processing
An overview of thermal processing is presented in this section.
Problems with traditional approaches to processing of conduction
heated foods are then noted and the need for online control is
stressed. Some of the past approaches in thermal processing of con
ductionheated foods that could also be used online are discussed.
Overview of Thermal Processing
Thermal processing consists of heating foods, filled and
hermatically sealed in containers, in pressurized steam retorts at a
constant retort temperature for a prescribed length of time followed
by rapid cooling. These process times are calculated on the basis of
sufficient bacterial inactivation in each container to comply with
public health standards and reduce the probability of spoilage to a
specified minimum. Estimation of bacterial inactivation requires
understanding of thermal death kinetics of bacteria, along with the
temperature history the food (containing the bacteria) experiences
during the process. Following is a brief introduction to the thermal
death kinetics of bacteria (Ball & Olson, 1957; Stumbo, 1965).
Thermal death of bacteria follows first order kinetics. If C is
bacterial concentration (number of bacteria per unit volume of food
material) at time t,
dC 1 (1)
dt (0.4 34)OT
where (0.434)*D is the inverse of the rate constant. If Ci is the
initial concentration and C is concentration after time At
C
Ci At
log = (2)
C
The constant D (known as 0 value), which is the inverse of the rate
constant (multiplied by the factor 0.434), is the time required to
reduce the bacterial population by one log cycle at a specified
temperature. The temperature dependency of the D value is given by
T T
D = Do 10 Z (3)
where Do is the D value at some reference temperature To. The Z
value describes the temperature dependency of the lethal rate and is
related to the activation energy. The reference temperature To is
generally taken as 1210C (2500F).
From bacteriological considerations, a thermal process is
specified by noting the time (Fo) necessary to reduce the bacterial
population, from the initial (Ci) to some final safe value (C) at
temperature 1210C (2500F). Thus
C.
F0 = D log (4)
When the heating process is not at the constant temperature of
121C, introducing the temperature dependency of Do and equation (2)
for bacterial death,
T121
Z At
Fo = D 10 z
T121
= 10 Z At (5)
At is the time required at temperature T to achieve the designed
bacterial reduction. As the food material is heated, its temperature
T becomes a continuous function of time t and the expression for Fo
(equation 5) becomes
t T121
Z (6)
F = f 10 dt (6)
When the food being heated has a spatially varying temperature
profile, the chosen temperature T, for thermal process calculation,
is the temperature at the slowest heating point in the food
material. Thus, for conduction heating in a cylindrical can, the
chosen temperature T is at the geometrical center of the can.
The designed bacterial inactivation for a process is specified
by Fd which is the time required to achieve the desired
0
inactivation of bacteria when the same is heated at a constant
temperature (T) of 1210C. When temperature T is a function of time
(T(t)) instead, the objective of thermal processing is to select
process time t for the particular heating condition, such that
the equivalent Fo received (as given by eqn. 6) is greater than or
equal to the design value F. That is,
Fo(t) Fd (7)
Application to conductionheating. When applying thermal
processing to conductionheated food in a cylindrical container,
Ball and Olson (1957) used the solution to the heat conduction
equation in cylindrical coordinates.
TTR (X2+ 2)at
R E Ampe m P Jo(mr) Cos ( z) (8a)
TOTR m p
The solution was first approximated to
T TR (X2 + X2)cat
S= A e Jo(xr) Cos(xz) (8b)
T0 TR
when sufficient time has elapsed, so that the first term in the
analytical series solution (equation 8a) dominates. Equation (8b)
was then rearranged to give
j (TR TO)
t = f log TR T(9)
TR T
where
S 2.303 (10)
(f a + (10)
and
j (r,z) = A Jo(xr) Cos (xz) (11)
From the transient heat conduction equation (9) which gives
temperature T as a function of time t and the definition of Fo
(equation 6) which gives Fo as a function of temperature T, it is
possible to find Fo as a function of time t. However this cannot be
performed analytically. Upon solving the two equations numerically,
it is possible to know temperature Th achieved at the end of heating
for a given value of Fo. Ball and Olson (1957) tabulated ThTR
values as a function of F0 for particular values of f and j. Tradi
tionally, in thermal processing of conductionheated foods in cylin
drical containers, these tables are used to find ThTR for a given
Fo value. Using this (ThTR) value and equation (9), the process
time t is found.
Online Control of Thermal Processing
In a production situation, the heating time t specified
(assuming a constant value of TR) to achieve a given Fo value would
need to be adjusted for arbitrary variations in the heating medium
temperature TR from its specified constant value. Researchers in the
past have attempted this in one of the following three ways
Precalculated tabulations. Correction factors are tabulated to
be applied to th for certain possible types of deviations in
TR(t). Commercial retort control systems often use this technique
(LOGTEC system, 1984). Such tabulations cannot be exhaustive,
since the process temperature TR(t) could go through any random
deviation. Instead of tabulating such correction factors, Giannoni
Succar and Hayakawa (1982) developed expressions for correction
factors to be applied to th when TR(t) goes through a step drop. In
their work, other types of deviations had to be approximated to
closefitting step functions, and multiple drops could not be
considered. Getchell (1980) described a control system that would
attempt to maintain the design heating medium temperature and sound
an alarm when critical low or high limits for the temperature are
exceeded. The operator would then attempt to make necessary changes
in th.
Direct Measurement. Temperature at the can center T(t) is
measured online. This eliminates the need to assume constant
heating temperature. Measured T(t) is used with equation 1 to stop
the heating cycle at time th such that
F t Fd (12)
Fo(th) F
This approach was used by Mulvaney and Rizvy (1984). They used a
thermocouplefitted test can in every retort batch to obtain actual
temperature T(t). The system proposed by Navankasattusas and Lund
(1978) also planned to measure can center temperature online. In
commercial practice these methods are costprohibitive with regard
to production efficiency, and are considered impractical.
Numerical Models. Instead of measuring T(t) at the can center,
T(t) can be predicted for arbitrary variations in TR(t). In
general, analytical solutions for T(t) are of limited use for
arbitary variations of TR(t). Teixeira and Manson (1982) described
a numerical finite difference approximation to the heat conduction
equation for predicting T(t) from TR(t). This is conceptually quite
different from the approaches mentioned earlier. With a finite
difference model, T(t) can be predicted for truly arbitrary
variations in TR(t) and, of course, use of this model eliminates the
need for having a test can with an actual thermocouple inserted in
it. This model was first described and verified against published
procedures by Teixeira et al., 1969; and later verified against
experimental data and other published procedures for timevarying
boundary conditions by Teixeira et al. (1975 a,b).
However, in the proposed retort control algorithm of Teixeira
and Manson (1982), every time a deviation took place, a new value of
th was found through simulation assuming the deviated retort temper
ature continued for the rest of the heating process. This is unnec
essary because any further changes in the retort temperature would
make the estimated th value useless and fresh estimations of th
would require repeated timeconsuming computer simulations.
Using the general idea of Teixeira and Manson (1982), the
feasibility of several control algorithms that can work online with
available computer hardware needed to be investigated, and led to
the following studies as part of the work carried out in this
project:
1) Ways to obtain the transient temperature solution faster,
using other available methods of solving the conduction
equation.
2) The need to incorporate a convection (Robin's) boundary
condition to reflect the true situation, rather than the
temperature specified (Dirichlet) boundary condition used
by Teixeira et al. (1975a,b) in solving the conduction
equation.
Conduction Heating of Canned Foods
Since the pioneering work of Ball and Olson (1957) mentioned
earlier, conduction heating of canned foods has been studied
extensively. Consequently there is a vast amount of literature
existing on the subject.
Studies Based on Analytical Solution
The majority of the works on conductionheating of canned foods
were based on one of several different forms of the analytical solu
tion to the heat conduction equation. Also most of these solutions
were derived assuming, among other things, a step change in the
boundary (can wall) temperature from the initial to the heating
temperature and maintaining the boundary heating temperature
constant throughout the process. More general forms of boundary
conditions have also been formulated. Hayakawa (1964) developed the
general form for an analytical solution that could predict the
transient temperature distribution inside a cylindrical can of
conductionheated food when the boundary (can wall) temperature was
any Laplacetransformable function of time.
In practice, the heating medium temperature, which affects can
boundary temperature, is not under complete control and can go
through unexpected variations from its design value. For such
arbitrary variations in the boundary temperature, analytical
solutions are of limited use.
Numerical Studies
A finite difference analog of the heat conduction equation can
also be used to achieve the transient temperature of food at finite
but small increments of time. Such methods can utilize the actual
boundary temperature "read" for the incremental time and calculate
the interior temperature values based on this boundary temperature.
Teixeira et al. (1969) were the first to use a forward difference
explicit finite differencing method to study the transient
temperature distribution inside the cylindrical can. Use of the
explicit finite differencing severely limited the choice of grid
size and time increment, due to stability requirements for such a
method. In their later work (Teixeira et al., 1975a) an alternating
direction explicit (ADE) method (Saul'yev, 1957; Allada and Quon,
1966) was used that was unconditionally stable but still explicit.
This method was first applied to simulate thermal processing of
canned foods by Manson et al. (1974). Teixeira et al. (1975a,b)
compared outputs from the ADE method for using constant boundary
temperature and sinusoidal timevarying boundary temperature with
the results from corresponding analytical solutions. The analytical
results for constant boundary temperature agreed within 0.30C
(0.50F) with the results from the ADE method, except at interior
points near the surface after small initial time (less than 3
minutes). Similar agreement between the results of the analytical
and the ADE formulation was also obtained when a sinusoidally
varying boundary condition was used.
Teixeira (1971) also compared the ADE finite difference outputs
for conduction heating in a can with experimental data for constant
boundary temperature and for step increases in the boundary tempera
ture during heating. The agreement was good except at initial times.
In addition to this ADE method, other methods of solution to
the conduction equation needed to be investigated from the point of
view of their suitability of computation in realtime and ability to
accommodate arbitrary boundary temperature variations.
Natural Convection Heating of Canned Foods
Internal (in enclosures) natural convection flow problems have
received much less attention than the external ones. Also, internal
flow problems are generally more complex than the external ones.
Fewer works have been reported on internal flows in cylindrical
coordinates. Of the studies on natural convection heating in
cylindrical enclosures, only four studies could be located where
they were in reference to natural convection heating during the
thermal processing of liquid foods. Such thermal processing is
characterized by a rather large Grashof number, generally small
overall dimension of the container and water as the liquid. Perhaps
the food processing situation is unique to general heat transfer
studies, which would explain the fact that only one such study could
be located in the more common heat transfer literature.
Previous Works on Canned Liquid Foods
The rate of heat transfer in canned liquid foods was studied
using correlations, and slopes of timetemperature curves with
limited success (Blaisdell, 1963; BarreiroMendez, 1979). Possible
flow patterns were postulated by Blaisdell (1963) using dye and
aluminum powder tracers.
Stevens (1972) first tried to study detailed temperature and
flow patterns during natural convection heating of canned liquid
foods. He used an explicit finite difference formulation from the
work of Torrance and Rockett (1969). The can was treated as a ver
tical cylindrical enclosure. The fluid (ethylene glycol) inside the
can was initially at rest at a uniform temperature of 300C. The can
was then placed in hot water bath at about 1000C. The top surface of
the fluid was free (the can had headspace). The actual temperature
variation on the can boundary was noted and was used in the
numerical calculations. Temperature values were recorded at ten
different points inside the container. The recorded temperature
values agreed poorly with numerically calculated temperatures. The
numerical inaccuracy was attributed to insufficient grid points
which could not resolve the details of the flow. Several sources of
experimental error included inaccurate setting of boundary
temperature, nonzero initial velocity field, errors in thermocouple
measuring circuit and inaccurate placement of the thermocouples.
Some of these factors were thought to have caused three dimensional
movements in the experiment and could not be picked up by the two
dimensional numerical model.
Detailed temperature and velocity profiles during natural
convection heating of liquid foods in a cylindrical can were
extensively investigated by Hiddink (1975). The numerical solution
technique used was the explicit finite difference technique of
Barakat and Clark (1966) discussed later. The initial temperature
was considered uniform throughout. The top, bottom and sidewall were
suddenly raised to the heating temperature. The top surface of the
liquid was considered free (having headspace) and was thus treated
as thermally insulated. Several test fluids of low and high
viscosities (water, sucrose solution and silicone fluid respec
tively) were used. The numerical results included temperature and
streamline patterns. The overall flow patterns were visualized using
a "particle streak method." In this method small glittering
particles are suspended in the liquid. The particles are illuminated
by a flat narrow slit of light, exposed against a dark background
and photographed. Detailed study of flow in the boundary layer was
observed using a laserdopper velocimeter which is an improved
version of the "particlestreak method" mentioned above. Temperature
values were measured by inserting thermocouples at several points in
the liquid. The temperature variation on the top free surface was
observed using a Thermovision Infrared camera. Radial velocity
profiles were presented at different heights in the liquid and at
different times during heating. His experiment and other works (Chu
& Goldstein, 1973) had showed convective heat transfer from the
bottom, as expected, when the bottom was also heated. His numerical
results, however, could not pick up the convective eddies generated
by the bottom heating. Thus his numerical results had poor corre
spondence with his experimental results.
A finite element method was used by Engelman and Sani (1983) to
study inpackage pasteurization of fluids (beer) in bottles. In the
process studied, bottles entered the pasteurizer at 1.70C(350F) and
passed through several progressively hotter zones of hot water
spray, which raised the product temperature to 600C(1400F). This
temperature of 600C is maintained in the holding zone. The product
then passed through several progressively cooler zones, which low
ered the temperature to 70 800F. Since the product was fluid and
the package was not agitated, it was primarily heated through nat
ural convection. Experimental details for this study was provided by
Brandon et al. (1981). The bottle geometry was complex (figure 1)
and so were the boundary conditions. The bottom of the bottle was
20 SECONDS
42 SECONDS
Figure Isotherms and velocity profiles in a. glass
bottle at various heating times (after Engelman
and Sani, 1983)
treated as insulated and wall temperature specified as a function of
time. The numerical results were in good agreement with experimental
data.
Sanil (personal communication, 1985) also developed a similar
model for pasteurization of beer in cylindrical metal cans (aspect
ratio =4.2). Water was the test fluid (Pr=6.81). The initial
temperature was 350C and the sidewall boundary condition used was
approximately 600C. The top of the water level was free (the can had
headspace) and was treated as insulated. The temperature at the
bottom of the can was a function of radial distance. Temperature and
velocity profiles from this study are given in figure 2. Further
quantitative details on this study were not available. These results
were not compared with experiments and also not published.
Other Similar Works on Natural Convection Heating in Cylindrical
Enclosures
One of the earliest of the numerical studies on the subject of
natural convection heating in cylindrical enclosures was performed
by Barakat and Clark (1966), who were looking at thermal strati
fication and associated processes in liquid propellant tanks for
application to space flight. The vertical cylindrical enclosure,
they studied, was filled with water, and arbitrary temperature and
heat flux variation were imposed on its wall, while its bottom was
kept insulated. At the top free surface of the liquid, heat transfer
to the ambient was considered negligible and thus insulated. They
1R. L. Sani, CIRES, Campus Box 449, University of Colorado, Boulder,
Colorado 80309
asseasa as a I1.1
6 l 8 * * **  
Gsilla.
\ fIlimII s alllII1*
i11 a 6 ) as t i11
(LiI~LaaliIll)i
S\* * a a ) \\ .\
IS S S I s a
allalIilalasIla
Figure 2. Isotherms and velocity profiles in a can at
various heating times (after Sani, 1985)
o
i *:: : : ;illI;
I .. .. . ..
: .: : : : .: ;: : .:: : ;i : : : : : : :
. . . . . . ..
. . .. *... .. . . .. ... il**
. . * ** . . . . . . .
I1 i
S' . * **I( * * * . i
I I * **CII' i I .. t * ** *
* * * i .*....... ..** I
* * * * * e . * * * **
S*j * i* * , * * * * **
i * * s o il* i * * * * *I*
I S S **S * * 1*1 * * ****
*. * *** * * *I ***
I . . . . .. I* . S * ** *
Ie I IIil" ** * * *
S.* * * .. * *
Ir \ B r I I I .
F r ** * *
i* i ** CI is *
aa I I * * * l* I
.i gu ,2 con ,d.,*
Figure 2. contd.
had numerically calculated and plotted the isotherms and streamlines
inside the enclosure. In their experimental work, they recorded
temperature values at ten different points inside the cylindrical
enclosure. The observed temperature values deviated from the
numerical results by 10%. This discrepancy was attributed to the
small nonsymmetry in the imposed wall temperature which might have
produced three dimensional effects and, also, to heat losses from
the bottom and top of the container.
Torrance (1968) compared five separate finite difference
formulations for the problem of natural convection in cylindrical
enclosures. The five formulations were respectively by Barakat and
Clark (1966), Wilkes and Churchill (1966), Fromm (1964) and two of
his own. The physical problem chosen was the study of natural
convection flow induced in a vertical cylindrical enclosure by a
small hot spot centrally located on the bottom. One of his two
formulations was explicit and the other was implicit. The explicit
formulation was shown to be numerically superior in many ways to the
other four and it was also reported to have good. agreement with
experimental observations (Torrance and Rockett, 1969; Torrance et
al., 1969) for lower values of Grashof numbers. This was also the
formulation used by Stevens (1972) referred to earlier.
Kee et al. (1976) performed a numerical study of temperature
distribution and streamlines in a closed vertical cylinder
containing uniformly distributed heat generating tritium gas. The
walls were kept isothermal. An upwind method was used for the
convective terms. Numerical temperature values were compared with
experimentally determined temperatures at three axial and two radial
locations. Thermistor thermometers were used for temperature
measurement. The calculated temperature values were reported to be
within the limits of experimental error, over a range of Grashof
numbers.
Evans et al. (1968) measured the transient temperature field
during natural convection of fluid in a vertical cylinder due to
uniform heat flux. The top surface of the liquid was free and no
transfer of heat or mass was assumed at the liquidvapor interface
at the top. Dye tracers were used to qualitatively observe the flow
patterns. They reported symmetry about the vertical axis in the
temperature field. On the basis of the experimental observations,
the convective liquid flow was modeled in terms of three regions:
1. a mixing region at the top
2. a central core
3. a boundary layer rising at the heated wall
They developed three different equations to describe the three
regions. These equations were solved simultaneously to predict the
overall behavior of the system. They reported good agreement between
the model and their experimental data.
Shyu and Hsieh (1985) made a numerical study of unsteady
natural convection in thermally stratified water in a cylindrical
enclosure. They studied the effect of placement of insulation on
maintainence of the stratification. Transient flow and temperature
fields were provided for the study.
PROBLEM FORMULATION
The governing equations and boundary conditions for the natural
convection and conduction heating process together with the assump
tions are first described here. The online control problem for the
conduction heating process is then formulated.
Natural Convection Heat Transfer in Canned Foods
The equations governing natural convection heat transfer in
liquid foods together with the boundary conditions are first
presented. The Boussinesq approximation to the set of equations is
then discussed. To further facilitate the solution, new variables
are defined and the governing equations together with the boundary
conditions are transformed into a new set of equations.
Governing Equations and Boundary Conditions
In natural convection heat transfer, the driving force for the
liquid motion is the buoyancy caused by the density variations due
to change in temperature. The partial differential equations (PDEs)
governing such naturally convective motion of fluid in a cylindrical
space are the NavierStokes equations in axisymmetric cylindrical
coordinates (Bird et al., 1976) and are described below.
+ Tar
9r
+ u 3T
3z
 3
r +
DF
K l1 a 3
 ~r ( (r
p r ar
= + {
8z r
u")
av
ai )
0
aT ) 2T
r) + )
ar az2
3 3u
ar ar
+ 313
3r r 3r
1 (pV) + (p) =
r ar 8z
The boundary conditions are
sidewall
T=T1
32u
+ +
Dz2
(r )) + }
aZ2
u = 0
v = 0
centerline
O z < H
3T 0
3F
bottom wall
0 < < R
top wall
0 < < R
z= 0
z=H
T=T1
T=T1
=0 = 0
(17ad)
u =0
u = 0
v = 0
v = 0
Initially the fluid is at rest at uniform temperature
0 < < R
0 < < H
T=TO
u =0
v =0 (18)
Refer to list of symbols for definition of individual symbols.
As expected, the flow field (equations 14 & 15) is coupled with
the temperature field (eq. 13) through the density term. These
a
+ +
aF
av
p ( +
at
(13)
(14)
(15)
(16)
equations have been approximated for no viscous dissipation. Jaluria
(1980) discussed various situations under which viscous dissipation
can be neglected. All other fluid properties are treated as
constants except variation of density which is discussed below.
Boussinesq Approximation
The solution to the set of equations can be considerably
simplified using the Boussinesq approximation. In the Boussinesq
approximation; the density is treated as constant except in the body
force term. The density difference in the body force term is approx
imated as
p p = pB (T T ) (19)
and the density variation in the continuity equation is ignored.
This approximation requires both
B (T Ta) << 1 (20)
and ( )T g << 1 (21)
to be true. The second inequality is generally satisfied for
liquids. So is the situation with the first inequality. For
example, for natural convection in water, as in this problem
B = 0.0002 /C
T = wall temperature of the can = 1210C
Ta= initial temperature of liquid in can = 30C
27
Thus,
B (T Ta) = (0.0002)(12130) = 0.0182 << 1 (22)
Almost all the previous works on natural convection in similar
areas have used the Boussinesq approximation.
Transformation of the PDEs
To further facilitate numerical solution to the set of
equations, the primary variables are transformed to a new set of
variables. Vorticity (J) and stream function ( ) are defined in
terms of the primary variables u (axial velocity) and v (radial
velocity)
:a u ,a(23)
3z aF
u =1a v = (24a,b)
rF Fr a
Using vorticity so defined, the pressure (p) is eliminated from
the u and v momentum equations 14 and 15. These two equations gave
the parabolic transport equation for vorticity Next, using the
definition of (equation 23) and the definition of (equations
24a and 24b), the equation for is derived. This gave the vorticity
streamfunction formulation for the natural convection problem.
28
Nondimensionalization of the Variables
The variables in the equation are then nondimensionalized as
below.
z
H
H
T H2t
H2
r
r 
H
v H
a
(25)
T T1
1
eH
T0 T1
H2 
1
ir~F :
The resulting equations are the governing equations for natural
convection in the form used for the present study.
ae a(ue) 1 (rve) _a2 1 a aeo
T + z } = { TaiZ + r ) }
(26)
aw a(u) a3(vw) D a2W a 1 a(rw)
9T + { + } =T GrPrr + Pr{ yz+ F r ar (7)
(27)
 T + ( T )
1 a
u = 1 F
(29a,b)
1 a2p
r az
(28)
The boundary
variables c.
sidewall11
r=R/H
centerline
r=O
bottom wall
0
top wall
0
y conditions
an be written
0 < z < 1
0 z < 1
z=0
z=l
in the transformed and nondimensionalized
as
= 0
=a 0
@r= 0
o 0=O
0=0
W = 0
S= 0 u=0
^ 0 .0
S=0 u = 0
ar
S= 0 u=0
S= 0 u=0
(30)
with the initial conditions
R
0
W = 0
S= 0 u=0
The relative merits of solving the transformed variables (w,p)
as compared to the (u,v,p) system of primary variables were discuss
ed by Roache (1982). To solve the (u,v,p) system, another elliptic
Poisson equation would need to be formed for the pressure (p) vari
able, using the two equations for velocity, regardless of whether
transient pressure solutions are desired. Following are the relative
evaluations of the (u,v,p) system and the (w, ) system for two
1Explicit boundary conditions for vorticities on solid walls are
unavailable.
(31)
dimensional flows when no free surface is present, iterative methods
are used to solve the Poisson equation, and transient pressure solu
tions are not required (Roache, 1982).
1. The (w,( ) system needs to solve one parabolic equation for
w and one elliptic Poisson equation for i as compared to the
(u,v,p) system which requires the solution of two parabolic
equations for u and v and one elliptic equation for p.
2. With the (0, 4) system, two additional differential equations
(29a and 29b) are needed to get the velocities u and v, but
the u and v momentum transport equations are more complicated
than the w transport equation.
3. At noslip walls, the boundary conditions for u and v are
known explicitly walll = 0, wall = 0) as compared to wall
which is not known explicitly. This would be a great advantage
for using implicit methods for the (u, v, p) system but non
linear instabilities of the pressure term prevents the (u, v,
p) system from using implicit methods.
4. The elliptic equation that needs to be solved for the pressure
takes much longer to converge than the elliptic equation for
p (eq. 28).
Conduction Heat Transfer in Canned Foods
The PDE describing conduction heat transfer in a cylindrical
body of solidpacked food material is presented here. The boundary
conditions present during the thermal processing of canned foods are
described next.
Governing Equation
For strict conduction heating (no movement of the material) of
a homogeneous isotropic cylindricalshaped body of the food materi
al, the heat transfer is described by the classical heat conduction
equation
ae a2e 1 eo 9a2
TF ar + F F + z (32)
This is the energy equation (26) approximated for zero velocities
(diffusion term only). This equation is quite well behaved and easy
to solve. The resistance to heat transfer of the metallic walls of
the can was ignored.
Boundary Conditions
Canned foods are generally heated with hot water or steam and
cooled with water. Thus, the heat transfer coefficient at the out
side of the can wall could be comparatively small (for water) or
quite large (for condensing steam). In general, the conditions at
the boundaries can be written as
sidewall
r =R/H 0 z 1 3 + hH =0
+ ar o =O
centerline
r = 0 0 < z 1 = 0 (33ad)
bottom wall
3e hH
0 < r < R/H z = 0 az e= 0
top wall
0 < r < R/H z = 1 + hH = 0
z k
with the initial conditions
0 < r < R/H 0 < z < 1 9 = 1 (34)
Computer Control of ConductionHeating Processes
In controlling thermal processes, the objective is to meet the
designed level of bacterial sterilization (Fd) for the process,
irrespective of any retort temperature variation TR(t), and with a
minimum of overprocessing. The lethal effects of thermal processing
are achieved during the heating as well as cooling times (th and
tc). Thus, the objective function in thermal processing is to
S T121 + T121
minimize F(th +t) = O 10 dt + fth 10 dt (35a)
h
subject to the constraint: F (th+tc) > Fd (35b)
where the integral from equation 6 has been broken into separate
heating and cooling times. Fo(t) is the accumulated sterilization at
any time t. The first integral above is the contribution to the
sterilization Fo(t) from heating, and the second integral is the
contribution to the sterilization during cooling. Temperature T(t)
is taken to be the temperature of the slowest heating point in the
product. This is so that when the design F (F ) is satisfied at
this point, all other points in the product have also been satis
fied. For a cylindrical can, T(t) is the temperature at the can
geometric center at time t.
This transient temperature T(t) is a function of transient
retort temperature TR(t), the can wall outside heat transfer coeffi
cient (h), dimensions of the can (R and H), and thermal diffusivity
(a) of the product for a conductionheating food. Symbolically,
T(t) = f (TR(t), R, H, a) (36)
For conduction heating in a cylindrical can, the function f is the
heat conduction equation in cylindrical coordinates.
For a given can size and product, the two variables that can
possibly be controlled are the retort temperature TR(t) during
heating (saturated steam or hot water under pressure in the retort),
and the time period of heating th. Cooling is generally done using
available water at ambient temperature. This makes the contribution
to Fo from cooling (given by the second integral in equation 35a)
immune from further control once the cooling process is under way.
Also, both the heating and cooling temperatures can only be applied
to the boundary of the product making it a boundary control problem.
In the case of a process deviation, the heating medium tempera
ture may go through unexpected variations beyond the capacity of the
temperature controller. Thus in reality, only one variable can be
controlled. It is the time period of heating (th). The problem in
thermal processing is to specify th for arbitrary variations in
TR(t) during heating while the process is under way in realtime.
Instead of searching for the proper th through simulation every
time a deviation takes place as suggested by Teixeira and Manson
(1982), Fo(t) can be monitored in realtime using T(t) calculated
from the finite difference approximation to the heat conduction
equation (32) with measured boundary (retort) temperature TR(t). The
heating process could continue until time th such that
Fo(th) > F (12)
thus satisfying the required sterilization Fd for arbitrary varia
tions in retort temperature TR(t).
However, the contribution to Fo from cooling (the second
integral in equation 35a) should not be neglected for a conduction
heating product. Depending on the can center temperature at the
start of cooling, size of the can, and other factors in a conduction
heating situation, the contribution to Fo from cooling could be as
much as 40% of the total. This is'unlike convection heated products,
where the container is sometimes agitated and the product is
rapidly cooled, making a small contribution to Fo during cooling. To
avoid gross overprocessing of conduction heated product, therefore,
the cooling Fo must be considered. Also, as mentioned earlier, there
is no control over the possible contribution to Fo from cooling once
cooling is under way.
Thus, the contribution from cooling cannot be neglected, it is
not a constant, and it cannot be controlled. Therefore, it can only
be estimated before cooling is actually started. Such estimation of
cooling lethality can be done through simulation of the cooling
cycle assuming a constant cooling water temperature in the retort.
The objective here is to develop an algorithm to perform this simu
lation in realtime when other control actions are going on.
METHODOLOGY
Detailed development of the solution techniques for natural
convection and conduction heat transfer models are described in this
chapter. For the natural convection heating problem, formation of
finite difference equations, setting up of computational boundary
conditions, and development of the algorithm for the iterative
solution of the complete problem, are described first. Along the
same lines, solutions of the finite difference analog to the
conduction equation are described next. Experimental study of
conduction heating is then reported. With the thermal models
working, use of the conduction heating model for online control of
conductionheated food is then presented.
Modeling of Natural Convection Heating in Cans
Because of the inevitable coupling of the governing equations,
except for extremely simplified cases, most natural convection flows
of practical interest are too complex to be solved analytically.
Thus numerical computations are generally employed to obtain the
desired solution. The present problem of natural convection in a
cylindrical enclosure is no exception. Numerical solutions to this
and similar problems demand high speed digital computers with large
amounts of memory that were unavailable until the 1960s. Since then
numerical studies of natural convection heating in enclosures are
being reported constantly. Among the numerical methods, finite
difference methods are probably the most popular in computational
fluid dynamics. The same was chosen for the present study.
Several choices were available for the finitedifference
analogs, computational boundary conditions (treatment of boundary
data in the numerical solution, as compared to real physical
boundary values), and the technique and order of solving individual
PDEs in the overall computational cycle. The relative merits of
these choices are discussed, and appropriate selections are
described in this section.
Grid Generation and Grid Stretching
The discrete points at which finite difference equations are to
be solved constitute the grid system for the problem. For computa
tional efficiency, the number of grid points should be the minimum
that is required to resolve all the significant spatial variations
of the flow. Since this problem involves a rather large step
increase in the wall temperatures together with a high Rayleigh
number (Ra), strong boundary layers are expected to form near the
walls. This is true particularly for the side (vertical) wall. In
order to resolve these boundary layers adequately, a large concen
tration of grid points would be necessary near the walls. An alge
braic method was used to generate the nonuniform grid system. The
particular algebraic expression chosen was evaluated for the proper
ties listed below (Anderson et al., 1984; Ka'lney de Rivas, 1971).
Here (r,z) refers to the physical domain and (,n) refers to the
computational domain with AE and An as constants.
1. The mapping from (,n) to (r,z) must be one to one.
2. Grid lines should be smooth to provide continuous
transformation derivatives. This requires r(E) and z(n)
as well as r and d should be continuous over the closed
intervals 0
functions defined on a closed interval are bounded, this in
turn implies that dr and dz should be finite over the
whole interval. If dr becomes infinite at some point, the
d
mapping r(S) would give poor resolution near that point,
which cannot be improved by increasing the number of points,
since A r r () (AE) Similarly for dz/dn.
3. Close spacing of grids where large numerical errors are ex
dr 0 at the
pected. For this problem, that means = 0 at the
r=R/H
sidewall which ensures high resolution of the boundary layer
near r = R/H. Similarly at z = 0 and z = 1.
4. Absence of excessive grid skewness which can sometimes
exaggerate truncation errors (Raithby, 1976).
An algebraic transformation equation suggested by Kublbeck et
al. (1980) for natural convection in cavities was
1 (+ tan[T/2 (2z1) al) (37)
2 tani/2 a
where a is
stretching.
ultimately
mation z(n)
below.
a deformation parameter that varies the degree of
However, since the computer coding of the problem is
done in the physical domain, the inverse transfor
is more appropriate in this context which is given
z =1/2{1 +  tan ((2n 1) tan a)}
Tr a 2
(38)
This concentrates the grids on either end z = 0 and z = 1.
In the r direction, the above transformation is modified to
concentrate grids on one side (side wall, r = R/H) only. Thus
r = 2 tan1 ( tan o) 
r 2 H
(39)
To check if the transformations indeed satisfy requirements 1
through 4, we first see that both the transformations are one to
one. Next, looking into the transformation derivatives,
dr 2 R tanw/2 a
d 77 IFT 1 + 2 tan2f/2 a
(40a,b)
dz 1 2 2 tan n /2
n w a 1 + (2n 1)2 tan2Tr/2 a
which are continuous in 0 < 5 < 1 and 0 < n < 1 respectively.
Also
dz
Lt d I
a+1 dn z=0
S2 tanira/2
S 1 + tan2Lt
a+1 we 1 + tan2nr/2
=Lt 2tan (1 2) by L'Hospital's rule
a+1
= 0
Similarly
dz 2 tanWra/2
Lt = Lt = 0
a+ z=l a lr ao 1 + tanZna'/2
and
dr 2 tani/2 a
Lt = Lt 1+tan2/2= 0
a 2l ir=R/H = 1 7r (J 1 + tan27r/2 a
Thus the transformations provide closer spacing of grids near z = 0,
1 and r = R/H. Figure 3 shows the grid system with close spacing of
grids near the top, bottom and sidewall. As can be seen, there is no
excessive grid skewness.
Ka'lney de Rivas (1971) also showed that when conditions 2 and
3 are satisfied for any smooth function, "extra truncation errors"
due to the nonuniform grid system are improved from first order to
second order in An and AE Vinokur (1983) also discussed similar
types of stretching functions.
Discretization of the Parabolic Temperature and Vorticity Equations
Of the several finite difference formulations that are
available for any set of PDEs, careful choices needed to be made
depending on many factors including boundary conditions, geometry of
the problem, type of solution (steadystate or transient) desired,
range of parameters involved (particularly the Grashof number in
e .= 1
mj r2t Z
m, m1 ,j
m0, ?Z Zi2
or amj 0
ei ,0= Oi I
i,0o O1
iO
e~0= ei
= 0
SIr i 1 ~ I I I I I II
I 111111
j = 0 O0,j=
SO,j
O,j
I t I IIIIllllll lii
tI IIII IItlll I
i ,n=
e. =
i ,n=
i ,nf
0
1
24i ,nI
rn/2 r n r
n12n ni
j = n
1
29
S1,j
r~ Z
ii1
The grid system and the boundary conditions
for the cylindrical can
ii
I I I I I
Figure 3.
'"'
SI I I H IIIIIIIIIIIIlIII
.t  It  I I ill
4++++1
ilIFill
I I I I
natural convection problems), as well as realistic constraints on
available resources (Roache, 1982). For this study, the transient
temperature and velocity solutions are needed and the range of
Grashof number is 108 1010 Since the temperature and vorticity
equations are both parabolic and are very similar, they will be
discussed simultaneously in this section.
The performance of the various finite difference methods was
noted in terms of stability, truncation error, and conservation
property. Stability was defined as the controlled growth of the
roundoff errors (O'Brien et al., 1950). An uncontrolled growth of
the roundoff errors leads to instability in the finite difference
methods. The stability requirements often put restrictions on the
time step as a function of grid spacings. Empirical stability re
striction on grid size of the type Ax<8/jui,jl and Ay< 8/Ivi,jl
were also reported by Torrance (1968). Truncation error is the error
resulting from the discretization of the PDE. Conservation of any
physical entity (e.g. energy, vorticity) within the finite differ
ence grid system exists if the difference equations for the entity
are summed over all interior grid points and no spurious sources or
sinks of these quantities are found. Torrance (1968) found that
formulations by Barakat and Clark (1966) and by Wilkes and Churchill
(1966) did not conserve energy or vorticity and he did not recommend
the use of their equations. Another important property of finite
difference formulations is the transportive property. A difference
formulation is transportive if the effect of a perturbation in a
transport property is convected only in the direction of velocity.
Time derivatives. The treatment of the time derivative in
equations 26 and 27 can be broadly classified into implicit and
explicit methods. Explicit finite difference methods are generally
easier to understand and program. Explicit methods have been used by
Barakat and Clark (1966), Torrance (1968), Stevens (1972) and
Hiddink (1975) in studying natural convection heating inside a
cylindrical container. For most of the explicit methods, however,
there are stringent restrictions on the allowed time step, from the
stability considerations. In situations where large temperature
gradients, eddies and recirculations are present like in this study
grid points need to be fine enough to resolve the physical phen
omena. When using most of the explicit methods, such a fine grid
size could severely limit the allowed time step and thus require
prohibitive amounts of computer time.
Implicit finite difference methods involve more than one
advance timelevel value in the same equation and thus require
matrix inversion to calculate function values at a new time level.
These methods generally allow larger time steps due to better
stability properties. However, even though a larger time step is
allowed in the implicit method, it requires many more calculations
to solve that step. For an implicit scheme to save computer time in
the overall computation, it should allow several times larger a time
step than that allowed by an explicit method to achieve a given
accuracy critrion.a
aT. IP. Shih, Mech. Engg. Dept., Univ. of Florida. Personal Comm.
Instead of using a fully implicit set of equations, alternating
direction implicit (ADI) methods were introduced in companion papers
by Peaceman and Rachford (1955) and Douglas (1955). The advantage of
the ADI method over a fully implicit method is that the set of
equations, although implicit, is tridiagonal. Further the ADI method
as applied to linear equations has a formal order of accuracy of
O(At2,Az2,Ar2). The stability of this method is also unconditional
as determined by VonNeumann stability analysis. For irregularly
shaped regions, programming could get complicated (Roache, 1982) but
for simple rectangular regions, as is the present situation, that is
no particular problem.
The practical advantage of the ADI method over the explicit
method is, however, nothing like that indicated by VonNeumann
analysis (Roache, 1982). It is unstable both at large time steps
and/or insufficient convergence of boundary vorticity. In general,
many researchers have indicated that ADI methods do allow larger
time steps and faster overall computation, by a factor of two or
more, and furthermore allow second order accuracy in time. The ADI
method is a very popular method in computational fluid dynamics.
While solving the set of PDEs (2629) Barakat and Clark (1966)
compared explicit and implicit methods for their study and decided
against the use of the implicit method due to unavailability of the
vorticity boundary conditions at the wall (to be explained later).
They allowed vorticity boundary values to lag one time step, neces
sary for such explicit formulations. The lagging of boundary
vorticity is considered undesirable (Roache, 1982) and iterations
are required because of this situation. The iterations can also
serve the purpose of preserving second order accuracy in time of the
ADI method. The second order accuracy may otherwise deteriorate due
to the presence of nonlinear convective terms (discussed below).
Since the iterations are needed anyway, an implicit set of equations
can be solved, instead, utilising the iterations. As noted earlier,
Stevens (1972) and Hiddink (1975) also used the explicit
formulations of Barakat and Clark (1966) to study the present
problem and could not resolve the convective eddies present due to
the bottom heating. Since the ADI method may allow finer grid
spacings with less stricter restrictions on At, it may be able to
pick up the convective eddies present, without taking prohibitive
amounts of computer time. Thus the ADI method was chosen to solve
the temperature and vorticity equation in this study. The ADI
formulation is written below for the temperature equation 26. At the
"half" step
en+1/2 en a(uo) n+1/2 a(rve) n a2 n+ /2 1 n
AT/2 8z ar a z +r ar Ir r
and at the "full" step
n+1 en+ 1/2 (u) n+ /2 n(rvI+/2 1 n+1
e o +(ue) n+ a(rve) n al^e 1 n
At/2 + az I + r az= r arar
Linearization. The velocities in nonlinear convective terms
e.g a(ue)In+1/2
e.g. (u) should ideally be evaluated at the same time level
(n+ 1/2 in this case) as the (or ). However, since u and v
(n+ 1/2 in this case) as the 9 (or w). However, since u and v
depends on i which in turn depends on u this would couple all the
equations together and make it practically incomputable. To avoid
this, the nonlinear terms are linearized. Various ways of lineari
zing these terms have been discussed by researchers (Anderson et
al., 1984; Peyret and Taylor, 1983; Roache, 1982). A very simple and
common strategy in linearizing is to evaluate the velocity coeffi
cient at n time level. This is known as "lagging" the coefficients.
However in this case, the accuracy in time drops from second order
(for ADI) to first order. Roache (1982) noted that such "lagging" of
the velocities could still be second order accurate if the velo
cities u and v were slowly varying. An improvement over this
linearization procedure is to update the coefficients u and v as
iterations proceed. The procedure is repeated until the values u and
v converge. One way to achieve this would be to update un+1/2 and
vn+1/2 by solving the stream function equation at the "half" time
step of the ADI (Roache, 1982). Since solution of the Poisson stream
function equation is the most time consuming operation, it was
decided against this updating procedure.
Anderson et al. (1984) noted that "lagging" the coefficients
has been widely and satisfactorily used by previous investigators,
and, recommended the use of the same for its simplicity. Thus the
velocity coefficients u and v were "lagged" in this finite differ
ence study.
Convection terms. The presence of the Lonvection terms makes
the temperature and the vorticity equation qualitatively different
from other simple parabolic equations. Proper numerical treatment of
the convection term is important. Several different ways of finite
differencing these terms have been presented in the previous works.
These include spacecentered differencing and upwind or upstream
differencing. The upwind methods have several desired properties
that are not shared by the other finite difference formulations, as
applied to the convection term. Both the upwind methods are
transportive, meaning any perturbation is carried along only in the
direction of velocity. Both the methods are also conservative so
that they preserve the integral conservation relations in the
continuum equation. The first of the two upwind differencing methods
is as follows (Lilly, 1965; Forsythe and Wasow, 1960; Frankel, 1956)
(u) = uO. uO.
(ue) (u > 0)
az Az
(41)
uOi+1 ueO
u= i+1 i u(u < 0)
AZ
Accuracy of this method is 0(At,Az). This was used by Barakat and
Clark (1966) to solve equations 26 through 29. The second variation
of this upwind method is to treat the convection term as follows
(Gentry et al., 1966)
(ue) URR L UL
az AZ
where
OR= i (uR > 0)
= i+ (UR < 0)
(42)
eL = i (uL > 0)
= eO (uL < 0)
where uR and uL are velocities as explained below. This latter
upwind method also has a formal accuracy of 0(Az). However, it
behaves like it is secondorder accurate when the quantity
( 0 or w ) is slowly varying (Roache, 1982). Torrance (1968) demon
strated that the second upwind method is indeed superior to the
first upwind method (Roache, 1982). The second upwind method was
used in the present formulation.
The velocities uR and uL are at the right and left face
respectively of the grid volume as shown in figure 4. As will be
explained later, a kinematically consistent velocity formulation
(Parmentier and Torrance, 1975) was used. In this formulation,
transport velocities (uj and vij in figure 4) are defined normal
to the grid volumes. Accordingly,
f f
UR= ui,j UL= Uil,j
f f
R= i,j L ,jI
Following is a development of the finite differences for the
convection terms in the temperature equation. In the z direction,
a(ue) a(ue) an
3z an az
9(ue) 1
an az/an
R i,j L i1,j 1
An z z
Si+1 i1
2 An
z +2 z( UReij uLi,j
zi+1 zi_1
UR > O, UL > 0 (a)
Similarly
a(uo) 2
az zi+1 z1 (uR i+l,j uLi1,j)
z 2 z (URi+l,j uLei,j)
i+1 *i
z i+ z i l i, uL i,j)
^ i+1 1 1
UR < 0
UL > 0 (b)
UR < 0 UL < 0 (c)
UR > 0 UL < 0 (d)
To avoid the use of FORTRAN IF statements for checking the sign of
UR and uL in the equations (a) through (d), the equations are
combined as below
a(ue) (RUR ) i+l,j+(UR+uRlUL+ u)i,j(uL+ i1,j 2
t 2 z+ z i
1+1 
(43)
In the r direction,
3(rve) ag a(rve) 1 (rve)
ar ar 9 2r/la aE
1
r 1r
j+1 j1
2A9
vRrj+ I/i j vLrj 1/i ,j1
2
rj+l rj1Rrj+ /2 ,j
Simliarly
6(rve)
6r
2
=r+1 2 (vRr1Nj+ 1,j +l rJ. 1
rj+1 rj1 Rr+ 1/2 i + vLj ij
rj+12 rj1vR N / v i Lrj. If 1,]
VR < 0 vL > 0
VR < 0 vL < 0
VR > 0 VL < 0
Therefore, combining the equations
a(rvo)
ar
(VRIVR )r l 1/201 j+1+{VR+IVRI )r l/(LIvL )rr ji j(,L+ vLI )r J/i j1
rj+lrj1
(44)
VR>0
VL>0
VL
 v Lj 1i ,j1
Diffusion terms. Central differencing was used for the
diffusion terms in both the temperature and vorticity equation. The
diffusion terms in the temperature equation are finite difference
below.
a ao an a ae
3 (7) = Tz Tn az
0ae
1 TZ i+ 12,j
Saz/an An
a i
Tz 1/2,j
9i+1,j 1i,j i1,j
1 i+1 i zi zil
z z. i1 An
2i+1
2An
2 i__ Ei+l, ,j eil
z i+ zi1 i+ zi zi zi_
i+1 i +1 1 
(45)
1 a ae 1 ac a ae
ar a(r) F r ar (r
ae
rliJ^
A&
ao
ri/?
11
r ar
5
1 1
r rj+1 rj
j 2A
S ,j+1 ij ,j ,j1
rJI2 rj+ r 12 r. rj
4/ rJ+1 r j12 j 1
2. O e. .
r j(rj+1 r ) /2 j+1 rj/2 r 1
a20
az
(46)
The finite difference formulations for the time and space deriva
tives are then put together in the ADI formulation. For each of the
temperature and vorticity equations, the equations are written for
the "half" time step and for the "full" time step of the ADI method.
Algebraic simplifications are done to convert the equations at each
of the time steps to a tridiagonal systems of equations. A glimpse
of the tridiagonal matrices for the "half" and "full" step is then
included.
The tridiagonal systems of equations are solved using the back
substitution procedure (Smith, 1978). The time step At and the
mesh sizes Az, Ar were chosen in such a way that the coefficients
(ai, bi, ci in equation 48a) are all positive and the matrices
diagonally dominant (bi > ai + ci). Richtmyer and Morton (1967)
showed that these two conditions are sufficient to keep the round
off errors down, in the back substitution procedure. Introduction of
the second upwind differencing method satisfies the criterion that
coefficients ai, bi, and ci are all positive (see Kublbeck et al.,
1980). The second requirement of diagonal dominance was not investi
gated any further in this study. However, at all times whenever the
numerical values of the coefficients ai, bi, ci were visually
checked, the diagonal dominance was found to be satisfied.
ADI formulation for the temperature equation is now developed.
For convenience, the differential equation for temperature is
repeated. The finite differences corresponding to each of the
derivatives in the temperature equation are then assembled in the
same order.
a3 a(ue) 1 a(rve)
TF az +r ar
a2e 1 a a
S= z +TF (rr) }
At the "half" or intermediate step ("n+l/' written as "*" in short)
* n
9. .e. .
AT/2
{(eluR)Oi+1,jUn+lu uI +lul),_(j +lul)i *,j
Zi+lZi_1 I1uR R i U,
(V nlv n n n ,n+1+(v+i
'r(r j+lrj
*
n 0 n zn zn
r r j+ rjl1 rj+1 rj rj r j 1
r. r r r r r.r
3 j+1 j1 +1 j j j1
(47a)
At the "full" step
on+le .
0. J e
1,j 1,3
AT/2
1 n n + n n n n n 0*
z{( l i+1,j + u+ "I u"+ RuL ) i,j ul ( L i'1,j }
(V/20 +l I )r n+1 n r +1 n n n+1
S(v vj+ R Ij RL1 L r ij J 2vl )rj i 2ij1
+ rj (rj+lrj_1
z z
S2 il~ji
Zi+lzi L Zi+l zi
* *
. . .
1,3 i1,j
i ii
n+l n+1l
+ ( 2 r.0, .
2J j++1/ j+1 ij
r (rj+l r ) rj+ r.
3 j+1 31 j+1 3
,n+1 n+1
l j j1
'J'Jl
(47b)
Vn n ) r .? n ) n
(v L L /r iJ (vL L rl j 1
After algebraic manipulation, temperature equations at the "half"
time step
ai il ,j+bi i,jci i+,j = rhsi
(48a)
where
S zi+1 Zi+l i Zii1
1 + 1/ (n+U n n+1nj) + 1 1 +
S 1/ n 1 1 1
b = R RI) 2 i
i T z +1z. R R L + iz. zi zzi
C i : ,)
zi+1 11 '+'z +1 i
rhsi = dn +ee +fe
i1 ij1 1 i,j+l
1/2v n+ vnjl)r ' rj
d = r r I + 1FTr r rr 1
(r+l j+ 1 J 1 j J+ J 1
e =
vVc nl,_ rj _li( nl n r 1/
R R L v r i 1 rj + 2+.r
rj rj+, rj_1 rj ( j +I j ) I
5 1 r+rj rjrj
1 jVR )rj rj +/2
f andR
and
f
uR = ui,j
f
uL = ui
L 1,j
f
R= Vi,j
f
L ,j1
After algebraic manipulation, temperature equations at the "full"
time step
a n+l b l c +l = rhs
a i,j1 b1,j j i,j+l rhsj
(48b)
where
1/2(vn +n I )r r
j r.r. rr.r r. + r.
n Lin n
S(v +vnR )r (vL _V )rv j/ 1 r lr r 
S. r.(r + r ) jr. r. r.r.j+1
1I Rn n ]r r
j r. r. r. r. jl r. r .1rj)
j j3+1 1 J 3+11 j +1 3
rhsj = dl,j + eOj + fi+,j
u l n+ l 1
d = +
Zi+1zi1 Zi+lzii)(zZiz)
1 12 n n i n 1 1 1
(URn 1 RI ) i
f R +
zi+lzil (zi+lZi1)(zi+1zi)
and
f
uR = u j
R = ij
f
uL = ui.lj
L 1,j
f
VR = Vi j
f
L i ,j1
The tridiagonal matrices for the temperature equations at the
two ADI time steps are given below. At the "half" time step, the
equations are written for a column (going from the bottom wall to
the top wall). Specified boundary temperature is used both at the
top and bottom walls. *
b1 c1
a2 b2
a3
am2 bm_2 cm2
am1 bm1
rhsl+a1 0,j
rhs2
rhs3
rhsm2
rhsml+aml m,j
1 /I
(49a)
At the "full" step, equations are written for the rows (from
the centerline to the sidewall). At the centerline, the insulated
boundary condition (equation 75) is used due to symmetry. At the
sidewall, the specified boundary temperature is used.
blal cl i,1 rhs1
a2 b2 c2 i,2 rhs2
a3 b3 c3 ,i3 rhs3
an2 bn2 cn2 0i,n2 rhsn2
n+
an1 bn1 in1 rhsnl+cn1 i,n
n1i,n1 i,n
(49b)
Upwind differencing of the convective terms in the vorticity equation
a(uw) a(uo)
az) = similar to (
az az
(uRluRl)"i+l,j+ (uR+IURluL+IuLI)wi,j (uL+IuL)wil,j
zi+1 zi1
a(vw)
ar
(50)
(vR IVRI)wi,j+l + (VR+ IVRI VL+ IVLI)wi,j (VL+ IVLI)&i,j1
rj+1 rj1
j+1 j1
(51)
Upwind differencing of the diffusion terms in the vorticity equation
(52)
+2 i+ mij iij
zi+1 z1 i+1 zi zi z i1
a1 a(rM)3
as a a(r))
= 1' aTr ar
1
rj+1 rj1"
2A9
1 a(rw)
r+ 1/ ar
,j+ 1/?
AE
1 a(r O)
rj._ 1ar i J1/
2 rj+1i,j+l r mi
rj+ rjl rj+ 1/2(r+1 r.)
2 j+3+1 J
r.. j r. l )
 J I/2, j rj1
(53)
Central differencing of the buoyancy term
ar r+1 rj1
(54)
a2
az2
ADI for the vorticity equation
a + a(u) + (v+) a2 2 + a 1 a(rw)
S+ I ar } = GrPrD + Pr az ar+ r 3r()
At "half" or intermediate step,
* n
13i ,jWi ,j
AT/2
Zi+lZ 1Lj
1_ n nn l vn v nn nI
Iv v *v l ,
+(r 1rj 1I ,j+R1 R R L vn+ v ,j (vn j1
j+1 :1
n+1 n+1
1 e P
SGrPr2( i,j+1 i 2 Pr
j+1 J1 i+1 Zi1
n n
2 Pr rj+1 ij+l i,
rj+r ) r j+/2rj+ r )
'j+1 j+1 + /2 j+1 3
* *
mi+ljti ,j
i+lzi
i+1
* *
i ii
n n
rj i ,rj1 "ij1
rj. /2 rrj )
(55a)
At "full" step,
n+1 *
ST/2o.
AT/2
+ 1n n un+nln (un+unI
Z + (u I u I )uR i+l,j+(UR URI ULI i,_ UL il,j
i+1 1 l ,
1 n n n+1 n+lnln + n n+1 n lin+l
(rj +lrj1 )(vRv )wij+1 vR R L L ij L L L ij 1
n+1 n+1 *
GrPr2 ij 1 i ,j1 2 Pr Oi+1,j 'ij "Oi, i1,j
Sr +1 1 zi+zi1 zi+1zi z z i1
n+1 n+1
2 Pr rj+oi j+1 r ij
(rj+lrj1 r +1/2 (rj+lr) i
n+1 n+1
rjoij rj1 "i ,j1
r i/2 r j1r
(55b)
After algebraic manipulation
*
ai j+ bimj cii+,j = rhsi (56a)
where
1/2 (u +jul ) Pr
ai Z i+lZil zi+lziilzizi)
bi E z z (un+iuiun+tunI) + rz 1 +z 1
1 AT zi+ i1R R L i+ i1 +i+1 1i
1/2 (uR R ) Pr
c. +
1 z 
ci z i+1z1 ii+1z 1)(i+lzi
(n+l n+1
rhs. = D2. +Ewm .+Fwn 0.5GrPr2 ij+l j
S ,j1 1,i 3 i j+1 rj+1 rj
1/2 (v + vjI) Pr rj1
d = (v I+ L J
j+r 1 rj1/2 rj+1 rj_1 j r1
n n n nlr
1 i (VR+IVI IvI) Pr rj 1
1 R R L L Pr 1
e T= 1/2 r r r1 r r r rj _j/
S T/ j+1 j1 j+ rj1 2 rjlr] rj
1/2 (VRVR Pr rj+1
f=1 +
rj+1r1_ rj+1lrj1 Jrj /2 rj+lr )
After algebraic manipulation
aj n+l b mn(l n+l
a ,1 + b i, + c ,j+ rhs (56b)
where
1/2 (vn+ vj ) Pr rj_1
aj= (r + r. r r )Lr.r. )
j+1 j1 1/2 j+ j j1
(v n n n n
1 R+IVR vL+Iv L + Pr rj 1 +
bj= A + 2(j+lrjr (r rj )Irj (rj+rrj rj/2 1 rj1r )
j+1 j1 j+1 1
/2 (VR R Pr rj+
c.= r +
J r rj+l rjr.jlrjl/ rj rj
(en+1 n+1
rhs. dnw +en .+fn 0.5GrPr2 J+1 ,J
1,j 1, i+1,j r j+1 r
1/2 (u+I u~n I)
Zi+ +1_Zi 1 11)
z i+1z z. zi+1zi1)(zz
i+1 ii LI) 11
1/2 (n URI Pr
f =I +
Zi+lZ iI zi+lzi Zi+li)
The tridiagonal matrices for the vorticity equations at the two
ADI time steps are given below. At the "half" time step, equations
are written for the columns (from the bottom wall to the top wall).
Both the bottom and top wall vorticity values are used from previous
iteration.
b1 c1 01,j rhsl+alw ,j
a2 b2 C2 "2,j rhs2
 a3 b3 c3 m3,j rhs3
am2 bm2 cm2 "m2,j rhsm2
am_1 ,bm ml,j rhsml+aml om,j
(57a)
At the "full" time step, equations are written for the rows
(from centerline to sidewall). The centerline boundary vorticity is
specified (=0) and the sidewall boundary vorticity is used from
previous iteration.
n+1+l
b1 cl i,1 rhsl+al wi,O
a2 b2 C2 "i,2 rhs2
a3 b3 3 ,3 rhs3
an2 bn2 cn2 "i,n2 rhsn2
n+1
an1 bn i,n1 rhSnl+cnI "in
(57b)
Discretization of the Elliptic Stream Function Equation
The stream function equation (28), obtained from the definition
of a u is also a statement of the circulation theorem
f ~z ar
which is
f w dA = J ch (58)
area perimeter
element of area
Thus the difference formulation of the stream function equation has
to satisfy the circulation theorem in addition to conserving mass.
Difference approximations that conserve mass are more difficult to
achieve in a curvilinear coordinate system as compared to the
Cartesian coordinate system.
To satisfy both the circulation theorem and the conservation of
mass, Parmentier and Torrance (1975) introduced the concept of
tangential velocities parallel to the control volume faces and
transport velocities perpendicular to control volume faces (Figure
4). In what follows, their approach has been adapted for the non
uniform grid system.
The circulation theorem (equation 58) is written for a grid
volume around node (i,j) as shown in figure 4
ij ( 1/2z 1(rj rj
v t v t t u
S(vj i1,j)(j+/2r j 1/) (u juij1)(zi+1i/2z1/
t t t t
= 1i,j i ,j i,3 i',j1 (59)
ij zi+/2z 1/2 r.+ 2 r.
i + W zi 1/21/2
j+1/2
  Zi/2
I zjI
r]. l
Notations of various quantities defined on a
nonuniform grid system in cylindrical
geometry
Zi++I
iI
0
ii,'
I 1
r i./2 r] r
Figure 4.
iJ Z+Z/2
u 
ee
rj
The tangential velocities u ,j and v,j at the faces of the
grid volume are obtained by finite differencing equation 29
ut ij+ ij
ij rj+l/2 rj+ rj j
(60a,b)
v  ip
,t i+lj i ,j
i,j r (z i+ zi)
Substituting these definitions in equation (59)
1 i+l,jj 3J i,ji,j
I i+ Zi i
1 ri j+l ij 1iji,j1
rI r )I
rj+ 1/2r 1/2 rj+ 1/2rj+1r ) r.1/ rj 1
(61)
Note that the same finite difference equation (61) can be obtained
from central differencing the stream function equation (28) without
introducing the circulation theorem. This is shown in Appendix A.
After algebraic manipulation,
 "i = Ai,j+ B i,j + CJi + Di,j1 + Eij+1
(62)
where
B (z zi j z i ii rj { (rj+lrj) jr jr jl
+ 1 1/22+ 1 '/1
C+ r.zi.1
C [zi+ /2zi_ 1/2 )Zi+i )
D (rj+ 1/2rj. 1/2 Jrj. 1/2 [rrj1
E =1
(r+ 1/2r 1/2 rj+ 1/2 rj+1r
The grid point velocities are
interpolation of the transport veloc
then achieved
t
ities u. j and
J
th
r. r ? t 1 r1/21 r t
rj +/ j .r /2u 'i rj+ 1/2 rj 1/2 u 1
. i z l/? vt
vi 1j Zi+l/2 .i1/2 j,
+ i+z 1 zi t
zi+ 1/ zi 1/,
rough linear
t
i ,j
(63a)
(63b)
ui j
1 ,3
r i (z / J zi 1)
Second, transport velocities are defined normal to the grid volumes
*
in terms of a new set of stream function values i,
f i. ,j .,j1 ,f i, 1 i1,j (64a,b)
1, r ( 1/2j 1,j j z 1/2i zi1/2
The stream function values i,j are defined at the corner of the
grid volumes. These i,j values are related to the previous stream
function values ij through a consistency criterion which says
f
Lt u u = 0
Ar+0 i,
(65a,b)
Lt v v = 0
Ar+O
ti,j was defined as
,J 2 i,j + j) + X ( ,j+1+ i+1,j+1) (66)
then, using definitions of vj and v (equations 64b and
63b), it was shown that the consistency criterion (equation 65b)
require
j + xj = 1 (67)
J J
Similarly, using definitions of uj and u it was shown that
II >J ui j
the weakest condition that satisfies equation 65a was the recursion
relation
1 + (2j 1) xj (68)
j, 1(68)
Xj 2j + 1
Using boundary conditions on the axis, it was found that
X0 = 0.25 (69)
The transport velocities defined this way conserve volume flow
exactly and are consistent with the grid point velocity field as the
mesh is refined.
The transport velocities uj and v are the velocities
used in the energy and vorticity transport equations. The grid point
velocities u. i and vi, are the velocities plotted in the
velocity plots.
To solve for I from the vorticitystream function equation
61, we note that it is an elliptic Poisson equation as compared to
the parabolic equations for temperature and vorticity. Direct (non
iterative) methods to solve the equation are available that are
considerably (10 times) faster than an iterative method like the
successive overrelaxation technique (to be discussed later).
However, direct methods usually require large amounts of storage,
and available computer codes are often limited to rectangular or
simple domains or have restrictions on the type of boundary
conditions. Shih (in press) discussed suitability of direct methods
for "large" systems of equations. Direct methods also have round off
error propagation problem. With direct methods, roundoff errors can
be incurred at each mathematical operation, and simply accumulate
until final answers are obtained. When iterative techniques are
used, the presence of roundoff errors at the end of any given
iteration simply results in those unknowns being somewhat poorer
estimates for the next iteration. For practical purposes, roundoff
error in the final converged values of an iterative scheme is only
that accumulated in the final iteration (Hornbeck, 1975). A very
efficient direct method for solving Poisson's equation was described
by Schumann and Sweet (1976) and the code was available from the
National Council of Atmospheric Research (NCAR) software library.
However, on closer look the algorithm (GENBUN) was found to be
limited to coordinates transformed in one direction only. Since in
this study the boundary layers were expected to be present both
along sidewall as well as top and bottom wall, coordinate stretching
were considered necessary in either direction. This made the above
algorithm unsuitable for the present study.
A popular approach to solve such a steady state equation is to
formulate it as an unsteady equation and solve for the asymptotic
steady state solution. ADI methods have been used for this purpose.
These methods were known to be faster than the successive over
relaxation method (to be discussed later) for a rectangular region
(Birkoff et al., 1962). However, the optimum sequence of parameters
aR. Sweet. U. S. Dept. of Commerce. National Bureau of Standards.
325 Broadway, Boulder, Colorado. Personal Communication.
that makes the ADI method advantageous is not available in general.
Programming the ADI method for a nonrectangular region may get
complicated, and they are not known for sure to be faster (Roache,
1982).
Another popular method to solve the unsteady formulation of the
stream function Poisson equation is the successive overrelaxation
(SOR) method (Frankel, 1950). It is very simple and effective. In
SOR, the set of solutions are over relaxed by a factor 1 < n < 2
to speed up convergence. The exact value is found through simple
numerical experimentation. It is easier to program compared to ADI
methods. Most of the earlier researchers have used the SOR to solve
the stream function equation. In all the five formulations of
natural convection in a cylindrical geometry, studied by Torrance
(1968), SOR was used to solve the stream function equation. However,
it is time consuming and could end up taking 90% of the total time
spent to solve the problem. Nevertheless, SOR was chosen to solve
the stream function equation in the present study for its
simplicity. Following is the development of the SOR formulation.
Rewriting equation 62,
ij = C il ,j+Ci+,j+Dijl+Eij+l+ij
To get the SOR formulation,
k+1 1 k+1 k Dk+ + k (70)
i,j B 1,j i+l,j ij1 1i,j+l+ i,j
k+l k+ (71)
i,j ,j i,j
Using equation 71, stream function values ,j at all interior
points are updated. The entire sweep is repeated until the
convergence criterion
k+l k
Max ij 1 ,j (72)
i ,J
is satisfied. The value of e chosen for this study was 105. In a
*
transient calculation like this study, the optimum Q may change
*
from one time instant to another. However, the optimum Q found
through calculations during the very first time step was used
throughout the entire transient calculations.
Computational Boundary Conditions for the Finite Difference
Equations
A PDE such as equation 27 for w can describe a wide variety
flow situation. The flows (solutions) are distinguished only by the
boundary and initial conditions and flow parameters like the Grashof
number. Thus, specification of computational boundary conditions,
besides affecting numerical stability, greatly affects the accuracy
of the solution (Roache, 1982). Using the boundary conditions
already noted for the partial differential equations in the previous
chapter, the boundary conditions for the finite difference equations
are now derived.
The boundary condition for temperature is rather simple. When
the can wall temperature is specified, the computational wall
temperature value of the fluid next to the wall is made equal to the
specified value. On the adiabatic walls and at the centerline (which
is also considered adiabatic due to symmetry), a Taylor series
expansion gives
1 W + 0 An +1/2 j An2 + O(An3) (73)
ao
where stands for the derivative of temperature in the direction
normal to the wall and An being distance normal to the wall.
for = 0, (74)
for = +20W
OW+1 = W +1/2 820' An2
an2
= e + 0(An2) (75)
Thus if eW+1 = OW is set, wall temperature is evaluated from the
adjacent temperature and it is second order accurate even though
ao is only first order accurate (Roache, 1982).
In the ADI method, boundary conditions are needed at the
computational "half" step. Since wall temperature 0W is not a
function of time, the boundary temperature value at this "half" time
step is simply equal to eW.
Unlike temperature and velocities, vorticity and stream
function are not primary variables. This makes specification of
their boundary values difficult, particularly for vorticity.
Improper specification of boundary vorticity is found to be
destabilizing in general and the same was found in the present
problem. Vorticity is produced at noslip (wall) boundaries. It is
then diffused and convected into interior points. No explicit
expression for wall vorticity is available. It is related to the
interior stream function values by the vorticity streamfunction
equation (equation 28). However, this would make the vorticity
equation implicit with the stream function equation and the number
of simultaneous equations increase drastically. To avoid this,
Barakat and Clark (1966) chose explicit methods over implicit and
used vorticity boundary values calculated at the previous time step.
Wilkes and Churchill (1966) used the ADI method but similarly used
boundary vorticity values from the previous time step. Thus, in
their formulations, wall vorticity values lagged by one time step at
all times.
Roache (1982) discussed the effects of such lagging of boundary
vorticity values. It forces the use of a smaller time step since
larger time steps could make the solution inaccurate and unstable.
Since larger time steps are a major motivation for ADI, updating of
the boundary values through iteration is preferred. From the work of
Torrance (1968), Roache also noted that the convergence of the
boundary vorticity values actually imposed a timestep restriction
of the form AT < a / Az2 where "a" was some number dependent on the
problem and convergence requirements. Thus, due to the implicitness
of the boundary vorticity values, ADI actually has a timestep
restriction like the simplest of the explicit methods. The
definition of vorticity o was used at the boundaries to obtain
expressions for boundary vorticities. At the bottom boundary,
t
SO,j v,j
O,j (z l zO)
t t
UO,j UO,j1
( rj+ 1/rj 12
 1,j
r.z z12
since
t
Vo,j
t
uj= 0
_ ( ,jOj)
rj(z1 z0)
S0 from eqn 60a
and uj=j 0 from eqn. 60a
At the top boundary,
t
SVmj Vmlj
m,j= (z z
t t
Smj um,j1
( rj/ r.j
_ ml,j
r. (z z z z 1
J m m1 m m2
since
v .= 0
m ,j z 
S m ,jm1,j
Vml,jjr rn(ZmI)
 im1,j
 r.(z z )
I m m1
utj and u t , 0 from eqn. 60a
mij mll=,
(76a)
Ip,j
r z1
^J
(76b)
At the centerline,
t t t
i,0 Vil,0 ui,0 Ui,
i,0o z. i+ z.1, r i/ 0
=0 (76c)
since
t O and v, from eqn. 60b
1,0 1a1,0
t i 4,I i,0 2_ 1
io rl O) rli
ui,0=0
At the sidewall,
t t t
S,n Vil,n i,n i,n
i,n= (z+l/2 z l/ (rn rnl/
= in1 (76d)
r ni(r r )(r r
n n n1 n n
since
Vi, and vi 0 from eqn. 60b
1i ,n 1i l,n=
t i,n i,n1 1,n1
1i,n1 r/rn r r r n rl)
n /2 n n1 n_/2 n n1
u. =0
Identical expressions can be derived using Taylor series expansions
of equation 28 at the walls. This is shown in the appendix B.
These expressions are first order accurate. But it is the only
form consistent with the kinematic velocity formulation (Paramentier
& Torrance, 1975). The second order form using Taylor series
(Jensen, 1959) has also been used. However, sometimes the second
order form destabilized the calculations (Wilkes & Churchill, 1966)
and sometimes it was less accurate than the first order formulation
(Beardsley, 1969).
The vorticity boundary value at the computational intermediate
or "half" step of the ADI method is more difficult to set than the
corresponding temperature value because the vorticity boundary value
is also a function of time. The following simple way to define this
vorticity value on the boundary is often used and gives sufficiently
accurate results (Bontoux et al., 1978)
n+1/2m+1 n n+1m+l
W = 1/2(W+W ) (77)
where
n = wall vorticity value at time step n
m+l
n+1 = most recent estimate of wall vorticity at n+1
step using equation.
m+1
nmW2 = intermediate vorticity value between n and
(n+l)th time level at the (m+l)th iteration.
Israeli (1972) also used this equation to define intermediate (half
step) boundary vorticity though he used a different equation to
nm+l
calculate boundary vorticity wW
Convergence Criteria for Iteration
As discussed earlier, the set of finite difference equations
were solved iteratively until convergence. The measure of conver
gence could have several possibilities. Ideally, at convergence, all
the variables should have converged. For practical reasons, it would
be sufficient to check on the variable that converged the slowest.
This is so that when this variable has converged, all others would
have converged. However there are no general guidelines to choose
apriori the slowest converging variable. Interior velocity
components were suggested by Roache (1982) as a measure of
convergence. Frequently, researchers in the past have decided on the
vorticity boundary values as the variables to watch for the
iteration convergence of the entire problem. Roache (1978) referred
to tests on boundary vorticities as the most sensitive test for
convergence. The same variables were used to check for iteration
convergence in this study even though they may not have been the
variables that converged the slowest.
For the chosen variable, there are no definitive criteria for
the measure of convergence. A relative error criteria was used for
this study, which is
n+m+1 ~n+lm
< 6 (78)
n+1
at all solid boundaries. The boundary vorticities n+1 were found
from the interior stream function value using equation 76.
As an example of other measures of convergence, while studying
natural convection in a rectangular cavity, Kimura and Bejan (1984)
used the following convergence criteria on interior points
E E n+lm+1 n+1i
i j i ij < 6 (79)
1I1m+11 6
n+1
i j
together with a convergence criterion on the calculated Nusselt
number.
The choice of 6 is also problem dependent and no general
guidelines exist for it either. Even for very small values of 6,
the iteration of a solution may be stopped prematurely for certain
types of solutions which are not uncommon (Roache, 1982). Values
ranging from 103 to 104 were used in their study, and the
chosen value was the largest of these values that did not cause any
perceptible changes in the final solution.
Algorithm for the Iterative Solution
The algorithm for the iterative solution of the complete system
of finite difference equations (equations 47,55,61,64) is now
described. The grid system is generated first. The vorticity, stream
function and velocity values are set to zero which is their initial
condition (t = 0) since initially the liquid was stationary with
uniform initial temperature throughout the body of the liquid. The
nondimensional temperature values are set to one (initial temper
ature). The temperature boundary condition for t>0 is then enforced.
The set of equations are then solved for finite increments of time
(At) with iterations performed on boundary vorticity at each time
increment. This is described in figure 5 and explained as below. The
selection of time step the At is discussed under the RESULTS
section.
REPEAT OVER n
*n n"
lm
2=1.0
r
ne m 1.0
ITERATE OVER m
n
I r
=.
P +
r e1
n/2m= f, n, un+1m vn+lm, /2m)
Q 1 Q I
n+1m 2 n12m
0 Q f 2( .
n+1m n+1m
u Q v
"a a
n+lm
W )
n+lm
en + 1
n /2mo
wD = 0
2n1/2' f r n n4
n +3
w = f,(w, u
n+lm I
U) =O
r
n+1" f n 1/2
^ n+lm
) = Wr
ITERATE
m
n+lm

n+lm
, Q
,n/2m
, p 2
n+1m n+l1
uS vS
n1/2m n+1m
4 w 'a
n+1m
"r
f5 (n+1m
fo 0
(,n+lm
S 7 j,
n+m+l
= n+
=yw
m
S,, n+
+ (1y)w
r
UNTIL 
u(vn+1m
(u,v) l
r: solid boundary
Q: interior points
: centerline
I ^ *
m,w,w: temporary
storage
( n+1m
f r~l!
UNTIL DESIRED TIME REACHED
Figure 5. Algorithm for the iterative solution
of the set of equations
UNTIL 
n+m+1
n+1
"r
n+lm+l
Wr
n+1m n+lm
W( 5 O
1. The temperature equations (47a and 47b) were solved using
the ADI technique. This gave new temperature values after
incremental time At.
2. The vorticity equations (55a and 55b) were then solved
(also using the ADI technique). The interior vorticity
values calculated were provisional since old values were
used for boundary vorticity.
3. The stream function equations (61) were then solved for
interior stream function values using the point SOR
technique. This required only interior vorticity values.
Previous stream function values were used as an initial
guess.
4. Boundary vorticity values ( ab) at solid walls (top,
bottom and side) are first calculatedusing equation 76
based on the newly calculated stream function values.
However, before updating, the boundary vorticity values
(b) are underrelaxed as
m+1 m+1 m
n+1 n+l1 n+ (80)
"b = b + ( y) Wb
m+l
where n+ = boundary vorticity estimated using
equation and 76, with Oy
5. Steps 2 through 4 are repeated until the estimated wb
values (using equation 76) in two consecutive iterations
meets the convergence criterion
~n+lm+1 ~n+lm
ab %b
< 6 (78)
n+l
Wb
for top, bottom and sidewall. When the convergence
criterion is satisfied, temperature, vorticity and stream
function values are considered the transient solution at
this time instant.
6. The velocities are calculated from the stream function
values using equation 63 and 64
7. Steps 1 through 6 are repeated for every time increment
until the desired time is reached.
The underrelaxation of boundary values at step 4 has often
been considered necessary in the literature (Peyret and Taylor,
1983; Roache, 1978) for convergence of the solution. Under
relaxation is usually employed to make a nonconvergent iterative
process converge (Hornbeck, 1975). This was found quite true for the
present study. Without underrelaxation, the time step allowed for a
stable solution got severely reduced to the point that solution to
the set of equations became impractical because of the amount of CPU
time required. As the computations proceeded, there were two sets of
values in storage. These were the estimates of the temperature,
vorticity and stream function values at the present time step and
converged temperature, vorticity and velocity values at the previous
time step.
Coding of the Program and the Hardware
The algorithm presented in figure 5 for solving the vorticity
stream function formulation of the natural convection problem was
coded in FORTRAN 77 with individual modules for solving temperature,
vorticity, stream function, velocity and boundary vorticity values.
Structured programming practices were used. Use of GO TO statements
were avoided except in two situations to simulate a WHILE statement.
The WHILE statement was unavailable in the particular version of
FORTRAN 77. Single precision of the variables was known to carry 14
significant digits for the particular computer and the same was
taken to be sufficient for this problem.
The input data file for a typical run is in Appendix C.
The computer used was a CDC Cyber 170/730, with two 60 bit CPUs
crunching 60 bit words at 400 MHz, 20 parallel processors, and 1091K
of 60 bit word central memory.
The plotting was done using a plotting package from the
University of Illinois1. The convenience offered by the plotting
package cannot be overemphasized. The plots were viewed on a
Tektronix 4006 terminal and hardcopies were made on Tektronix 4662
plotters.
Modeling of Conduction Heating in Cans
The partial differential equation (PDE) governing conduction
heating of a cylindricalshaped solid was discretized and finite
1W. E. Dunn, Mechanical and Industrial Engineering Department,
University of Illinois at UrbanaChampaign. Personal Communication.
difference using an alternating direction explicit (ADE)
differencing method, which is explicit but still unconditionally
stable. The transient temperature values predicted by performance of
this ADE finite difference method was compared with analytical
solutions and experimental measurements.
Discretization and Solution of the Equation
Keeping in mind that the equation was needed to be solved in
realtime for online control purposes and the ADE method was
already tested (Teixeira et al., 1975; Teixeira, 1971), the same was
selected. The ADE formulation was first reported by Saul'yev (1957).
This was later translated and described by Allada and Quon (1966),
and was first introduced for applications to thermal process
simulation by Manson et al. (1974). The method as applied to the
diffusion equation is unconditionally stable (Saul'yev, 1964) for
temperaturespecified boundary condition. It has a formal order of
accuracy of O(At2, Ar2, AZ2 ). Being explicit, it does not require
inversion of matrices as would be required by an implicit method. In
a fluid flow problem, the ADE formulation is generally not preferred
since the presence of the convection term (in any difference form)
makes the complete formulation unstable or returns it to simple
explicit stability limits ( a / for a one
dimensional problem; Roache, 1982). Another shortcoming of this ADE
method is the fact that even though it is an explicit method, it is
actually implicit in boundary conditions. To handle this implicit
ness, either explicit expressions for the temperature values need to
be derived algebraically for each of the boundaries or the computer
has to solve a simultaneous set of equations at the boundaries (see
Appendix D). For Robin's boundary condition, the application of the
ADE method required separate algebraic expressions for the top, side
and bottom wall boundaries. Separate expressions were also required
for the top and bottom corners and the centerline. This is in sharp
contrast with the situation if the ADI method were used. Adjustments
would be necessary at the boundaries for the ADI method, but for use
of the ADE method it became clumsy (see appendix D).
For temperature specified at the boundary, this is trivial and
no additional manipulations as noted above are needed. Experiments
conducted later showed that because of the low thermal conductivity
of food (which is mostly water), the Biot number (= hR/k) is quite
large and it is the internal conduction resistance that is limiting.
Thus a convection boundary (Robin's) condition can be approximated
by a temperature specified (Dirichlet) boundary condition for this
purpose. The two shortcomings of the ADE formulation are, therefore,
of little consequence in the case of conductionheated food. The
equations developed for the convection boundary condition are given
in appendix D.
Coding of the Computer Program
The ADE formulation for transient temperature calculation using
temperature specified boundary condition was coded in PASCAL and
used as a subroutine in a real time control algorithm (to be
discussed later).
Experimental Studies
Experiments were performed to record the transient temperature
values inside cylindrical cans during conduction heating. The effect
of different surface heat transfer coefficients, that were obtained
by using different media (steam and water), on the transient
temperature response was also noted.
The experimental work consisted of series of heat penetration
tests on No. 303 cans filled with a 10% bentonite suspension as a
conductionheating food model system (Niekamp et al., 1984). Eight
cans were fitted with an Ecklundtype copperconstantan thermocouple
mounted midway up the can wall, so as to have the junctions located
at the geometric can center. In order to achieve the two different
levels of headspace vacuum, four of the cans were filled at 900C
and the other four cans were filled at 670C. The cans were
identified according to fill temperature and then processed with an
extended steamcook to stabilize the bentonite suspension. An
experimental retort was used that had been designed and built at the
University of Florida with the capability of using pressurized
water, pure steam, or steam/air mixtures under controlled rates of
recirculation. The retort was fitted with a special rack designed to
hold the eight cans suspended upright in the mainstream of the
retort chamber, as shown schematically in Figure 6.
The rack consisted of a stainless steel plate with eight
cicular holes through which the cans could be inserted. Special
fixtures were welded at the edge of each hole so that the cans could
