Title Page
 Table of Contents
 List of symbols
 List of Tables
 List of Figures
 Review of literature
 Problem formulation
 Results and discussion
 Conclusions and recommendation...
 Biographical sketch

Title: Numerical modeling of natural convection and conduction heat transfer in canned foods with application to on-line process control
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00089827/00001
 Material Information
Title: Numerical modeling of natural convection and conduction heat transfer in canned foods with application to on-line process control
Series Title: Numerical modeling of natural convection and conduction heat transfer in canned foods with application to on-line process control
Physical Description: Book
Language: English
Creator: Datta, Ashim Kumar
Publisher: Ashim Kumar Datta
Publication Date: 1985
 Record Information
Bibliographic ID: UF00089827
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: alephbibnum - 000872393
oclc - 14514671

Table of Contents
    Title Page
        Page i
        Page ii
        Page iii
        Page iv
    Table of Contents
        Page v
        Page vi
    List of symbols
        Page vii
        Page viii
        Page ix
        Page x
    List of Tables
        Page xi
    List of Figures
        Page xii
        Page xiii
        Page xiv
        Page xv
        Page xvi
        Page xvii
        Page 1
        Page 2
        Page 3
        Page 4
    Review of literature
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
    Problem formulation
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
    Results and discussion
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
        Page 124
        Page 125
        Page 126
        Page 127
        Page 128
        Page 129
        Page 130
        Page 131
        Page 132
        Page 133
        Page 134
        Page 135
        Page 136
        Page 137
        Page 138
        Page 139
        Page 140
        Page 141
        Page 142
        Page 143
        Page 144
        Page 145
        Page 146
        Page 147
        Page 148
        Page 149
        Page 150
    Conclusions and recommendations
        Page 151
        Page 152
        Page 153
        Page 154
        Page 155
        Page 156
        Page 157
        Page 158
        Page 159
        Page 160
        Page 161
        Page 162
        Page 163
        Page 164
        Page 165
        Page 166
        Page 167
        Page 168
        Page 169
        Page 170
        Page 171
    Biographical sketch
        Page 172
        Page 173
        Page 174
Full Text







Dedicated to the memory of my

grandfather and father


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 I-P. 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 over-emphasized. 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


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 Urbana-Champaign. 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.


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
On-line 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
Non-dimensionalization 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 On-line 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
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 On-line Control Logic...................141

CONCLUSIONS AND RECOMMENDATIONS..............................151
Recommendations............................................ 153

STREAM FUNCTION EQUATION................................154

LIST OF REFERENCES.......... ..... .............. ................ 161

ADDITIONAL REFERENCES............................. ....... ....... 167

BIOGRAPHICAL SKETCH................................... ....... 172


A = constant

a,b,c = coefficients in tri-diagonal matrices

f = inverse slope of time-temperature graph

g = gravity

Fo = sterilization Fo value
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 = non-dimensional velocity in vertical direction

= velocity in vertical direction
= -v = non-dimensional velocity in radial direction

= velocity in radial direction

- non-dimensional distance in vertical direction
= distance in vertical direction

= Z-value 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)
= non-dimensonal temperature
= first eigenvalue in z direction

= eigenvalues in z direction

= first eigenvalue in r direction

= eigenvalues in r direction


= 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 = non-dimensional time
= radial coefficients for kinematic consistency

= radial coefficients for kinematic consistency

= non-dimensional stream function = -
= stream function defined by v = -
r ar r az
= non-dimensional vorticity
av au H2 v U H2 -
az ar a a

= vorticity

= optimum relaxation parameter for SOR











= 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


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)



1. Sequence of line types for identification of stream
function contours in figures 10a-o........................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 Giannoni-Succar and
Hayakawa ................................................. 148



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 real-time computations for on-line
control of conduction-heated food.........................88

8. Flow diagram for computer control of retort opera-
tions with on-line 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


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




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 conduction-heated 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 on-line control in batch thermal processing of

conduction-heated 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 on-line




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

Conduction-Heated Foods

Traditionally, for conduction-heated 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 low-cost microcomputers for on-line process control

is becoming commonplace in the 1980s. Food industries are not an

exception to this situation. Use of microcomputers for on-line

control of conduction-heated food is thus quite conceivable.

Convection-Heated 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. In-package 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 convection-heated

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 on-line control), the primary objective of this study was

directed toward developing a heat transfer model for natural convec-

tion; followed by improvements to existing conduction-heat transfer

models for application to on-line process control as secondary



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 on-line process control appli-

cations, and

3. incorporate the conduction heating model in an on-line

process control algorithm for thermal processing of

canned food, and compare its performance against other

possible methods of on-line control.


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 conduction-heated 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 on-line control is

stressed. Some of the past approaches in thermal processing of con-

duction-heated foods that could also be used on-line 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

Ci At
log = (2)

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

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,

Z At
Fo = D 10 z--

= 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 T-121
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
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 conduction-heating. When applying thermal

processing to conduction-heated food in a cylindrical container,

Ball and Olson (1957) used the solution to the heat conduction

equation in cylindrical coordinates.

T-TR -(X2+ 2)at
R E Ampe m -P Jo(mr) Cos ( z) (8a)
TO-TR 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)

S 2.303 (10)
(f a + (10)

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 Th-TR

values as a function of F0 for particular values of f and j. Tradi-

tionally, in thermal processing of conduction-heated foods in cylin-
drical containers, these tables are used to find Th-TR for a given

Fo value. Using this (Th-TR) value and equation (9), the process

time t is found.

On-line 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
(LOG-TEC 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

close-fitting 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 on-line. 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

thermocouple-fitted 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 on-line. In

commercial practice these methods are cost-prohibitive 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 time-varying

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 time-consuming computer simulations.

Using the general idea of Teixeira and Manson (1982), the

feasibility of several control algorithms that can work on-line with

available computer hardware needed to be investigated, and led to

the following studies as part of the work carried out in this


1) Ways to obtain the transient temperature solution faster,

using other available methods of solving the conduction


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


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 conduction-heating 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

conduction-heated food when the boundary (can wall) temperature was

any Laplace-transformable 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 time-varying 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 real-time 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 time-temperature curves with

limited success (Blaisdell, 1963; Barreiro-Mendez, 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, non-zero 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 laser-dopper velocimeter which is an improved

version of the "particle-streak 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 in-package 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



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


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

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

\ fIlimII s alllII1*

i11 a 6 ) as t i11

S\-* |* a a ) \\ .\

IS S S- I s a

Figure 2. Isotherms and velocity profiles in a can at
various heating times (after Sani, 1985)

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 non-symmetry 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


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 liquid-vapor 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.


The governing equations and boundary conditions for the natural

convection and conduction heating process together with the assump-

tions are first described here. The on-line 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 Navier-Stokes equations in axisymmetric cylindrical

coordinates (Bird et al., 1976) and are described below.

+ Tar

+ u 3T

- 3
r +

K l1 a 3
- ~r ( (r
p r ar

= + {
8z r

ai )


aT ) 2T
r) + )
ar az2

3 -3u
ar ar

+ 313
3r r 3r

1 (pV) + --(p) =
r ar 8z

The boundary conditions are



+ +

(r )) + }

u = 0

v = 0


O z < H

3T 0

bottom wall

0 < < R

top wall

0 < < R

z= 0




=0 = 0


u =0

u = 0

v = 0

v = 0

Initially the fluid is at rest at uniform temperature

0 < < R

0 < < H


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

+ +

p (- +





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



B (T Ta) = (0.0002)(121-30) = 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


: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

stream-function formulation for the natural convection problem.


Non-dimensionalization of the Variables

The variables in the equation are then non-dimensionalized as




T H2t

r -
v H


T T1
T0- T1

H2 -

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- } = { T-aiZ + r- ) }


aw a(u) a3(vw) D- a2W a 1 a(rw)
-9T + { +- } =T -GrPrr + Pr{ -yz+ F r ar (7)


- T + ( T )

1 a
u = -1 F


1 a2p
r az


The boundary

variables c.





bottom wall

top wall


y conditions

an be written

0 < z < 1

0 z < 1



in the transformed and non-dimensionalized


= 0

=a 0
@r= 0

o 0=O


W = 0

S= 0 u=0

^ 0 .0
S=0 -u = 0

S= 0 u=0

S= 0 u=0


with the initial conditions

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


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 no-slip 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 solid-packed 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 cylindrical-shaped body of the food materi-

al, the heat transfer is described by the classical heat conduction

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


r =R/H 0 z 1 3 + hH =0
--+ ar o =O

r = 0 0 < z 1 = 0 (33a-d)
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 Conduction-Heating 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 T-121 + T-121
minimize F(th +t) = O 10 dt + fth 10 dt (35a)

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 conduction-heating 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 real-time.

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 real-time 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 real-time when other control actions are going on.


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 on-line control of

conduction-heated 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 finite-difference

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
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
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 (2z-1) al) (37)
2 tani/2 a

where a is



mation z(n)


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


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


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


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

= 0


dz 2 tanWra/2
Lt = Lt = 0
a+ z=l a lr ao 1 + tanZna'/2

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 non-uniform 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 (steady-state or transient) desired,

range of parameters involved (particularly the Grashof number in

e .= 1
mj r-2t Z
m, m-1 ,j
m0, ?Z- -Zi2

or amj 0

ei ,0= Oi I
i,0o O1
e~0= ei

= 0-

SIr i 1 ~ I I I I I II

I 111111

j = 0 O0,j=


I t I IIIIllllll lii

tI IIII IItlll I

i ,n=
e. =
i ,n=

i ,nf


-24i ,n-I
rn-/2 r n- r-
n12n n-i

j = n

r~ Z

The grid system and the boundary conditions
for the cylindrical can



Figure 3.



.t - It -- I I ill




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

round-off errors (O'Brien et al., 1950). An uncontrolled growth of

the round-off 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 time-level 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. I-P. 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 Von-Neumann 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 Von-Neumann

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 (26-29) 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 non-linear 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-
A-t/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 space-centered 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
uOi+1- ueO
u=- i+1 i u(u < 0)

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


OR= i (uR > 0)

= i+ (UR < 0)

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 second-order 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= Ui-l,j

f f
R= i,j L ,j-I

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 i-1,j 1
An z z
Si+1- i-1
2 An

z +2- z-( UReij uLi-,j
zi+1- zi_1

UR > O, UL > 0 (a)


a(uo) 2
az zi+1- z1 (uR i+l,j uLi-1,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) (R-UR ) i+l,j+(UR+uRl-UL+ u)i,j-(uL+ i-1,j 2
t 2 z+ z i
1+1 -


In the r direction,

3(rve) ag a(rve) 1 (rve)
ar a-r 9 2r/la aE

r -1r
j+1 j-1

vRrj+ I/i j vLrj- 1/i ,j-1

rj+l rj1Rrj+ /2 ,j



=r+1 2 (vRr1Nj+ 1,j +l rJ. 1

rj+1 rj-1 Rr+ 1/2 i + vLj- ij

rj+12 rj-1vR N -/ v i Lrj. If 1,]

VR < 0 vL > 0

VR < 0 vL < 0

VR > 0 VL < 0

Therefore, combining the equations


(VRIVR )r l 1/201 j+1+{VR+IVRI )r l/(L-IvL )rr ji j(,L+ vLI )r J/i j-1




- v Lj- 1i ,j-1

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


a ao an a ae
-3 (7) = Tz Tn- a-z
1 TZ i+ 12,j
Saz/an An

a i
Tz- 1/2,j

9i+1,j 1i,j i-1,j
1 i+1 i zi zi-l
z z. i1 An

2 i__ Ei+l, ,j ei-l
z i+ zi1 i+ zi zi zi_-
i+1 i- +1 1 -


1 a ae 1 ac a ae
ar a(r) F r ar (r-




r ar

1 1
r rj+1 rj-
j 2A

S ,j+1 ij ,j ,j-1
rJI2 rj+ r 12 r. rj-
4/ rJ+1 r j-12 j -1

2. O e. .
r j(rj+1 r ) /2 j+1 rj/2 r -1



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


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 +T-F (r-r) }

At the "half" or intermediate step ("n+l/' written as "*" in short)

* n
9. .-e. .

{(e-luR)Oi+1,jUn+lu uI +lul),_(j- +lul)i *,j
Zi+l-Zi-_1 I1uR- R i U,

(V nlv n n n ,n+1+(v+i
'r(r j+l-rj


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 j-1 +1 j j j-1


At the "full" step
on+le .
0. J -e
1,j 1,3

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 2ij-1
+ rj (rj+l-rj-_1

z -z
S2 il~ji
Zi+l-zi L Zi+l- zi

* *
. .- .
1,3 i-1,j
i i-i

n+l n+1l
+ ( 2 r.-0, .
-2J j++1/ j+1 ij
r (rj+l -r ) rj+ -r.
3 j+1 3-1 j+1 3

,n+1 -n+1

-l j j-1


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 i-l ,j+bi i,j-ci i+,j = rhsi



S zi+1- Zi+l i- Zi-i1

1 + 1/ (n+U n n+1nj) + 1 1 +
S 1/ n 1 1 1
b =-- R- RI) 2 i
i T z +1-z. R R L + i-z-. --zi z-zi

C i : ,)
zi+1 1-1 '+'z +1 i

rhsi = dn +ee +fe
i1 ij-1 1 i,j+l

1/2v n+ vnjl)r -' rj
d = r- r I + 1FTr -r r-r 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 rj-rj

1 -jVR )rj rj +/2
f andR


uR = ui,j

uL = ui
L -1,j

R= Vi,j

L ,j-1

After algebraic manipulation, temperature equations at the "full"
time step

-a n+l- b l c +l = rhs
-a i,j-1 b1,j j i,j+l rhsj



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 .1-rj)
j j3+1 -1 J 3+1-1 j +1 3

rhsj = d-l,j + eOj + fi+,j

u l n+ l 1
d = +
Zi+1-zi-1 Zi+l-zi-i)(zZiz-)

1 12 n n i n 1 1 1

(UR-n 1 RI ) i
f R +
zi+l-zi-l (zi+l-Zi-1)(zi+1-zi)


uR = u j
R = ij
uL = ui.lj
L -1,j

VR = Vi j

L i ,j-1

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

-am-2 bm_2 -cm-2
-am-1 bm-1

rhsl+a1 0,j


rhsm-l+am-l m,j
1 /I


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.

bl-al cl i,1 rhs1
-a2 b2 -c2 i,2 rhs2
-a3 b3 -c3 ,i3 rhs3

-an-2 bn-2 -cn-2 0i,n-2 rhsn-2
-an-1 bn-1 in-1 rhsn-l+cn-1 i,n
n-1i,n-1 i,n


Upwind differencing of the convective terms in the vorticity equation

a(uw) a(uo)
az) = similar to (
az az

(uR-luRl)"i+l,j+ (uR+IURl-uL+IuLI)wi,j (uL+IuL)wi-l,j
zi+1- zi-1



(vR- IVRI)wi,j+l + (VR+ IVRI VL+ IVLI)wi,j (VL+ IVLI)&i,j-1
rj+1- rj-1
j+1 j-1


Upwind differencing of the diffusion terms in the vorticity equation


+2 i+ mij- i-ij
zi+1- z-1 i+1- zi zi- z i-1

a1 a(rM)3

as a a(r))
= 1' aTr ar

rj+1- rj-1"

1 a(rw)
r+ 1/ ar

,j+ 1/?

1 a(r O)
rj._ 1ar i J1/

2 rj+1i,j+l- r mi
rj+ rj-l rj+ 1/2(r+1 r.)
2 j+3+1 J

r.. j- r. l )
- J I/2, j- rj-1


Central differencing of the buoyancy term

ar r+1- rj1



ADI for the vorticity equation

a + a(u) + (v+) a2 2 + a 1 a(rw)
S+ I ar } = -GrPr-D + Pr az ar+ r 3r()

At "half" or intermediate step,

* n
13i ,j-Wi ,j

Zi+l-Z -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 j-1
j+1 :-1

n+1 n+1
1 -e P
S-GrPr2( i,j+1 i- 2 Pr
j+1 J-1 i+1 Zi-1
n n
2 Pr rj+1 ij+l- i,
rj+-r ) r j+/2rj+ r )
'j+1 j+-1 + /2 j+1 3

* *
mi+lj-ti ,j

* *

i i-i

n n
rj i ,rj-1 "ij-1
rj. /2 r-rj- )


At "full" step,

n+1 *

+ 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 +l-rj-1 )(vR-v )wij+1 vR R L L ij L L L ij- 1

n+1 n+1 *
-GrPr2 ij 1 i ,j-1 2 Pr Oi+1,j 'ij "Oi, i-1,j
Sr +1- -1 zi+-zi-1 zi+1-zi z -z i-1

n+1 n+1
2 Pr rj+oi j+1- r ij
(rj+lrj-1 r +1/2 (rj+l-r) i

n+1 n+1
rjoij -rj-1 "i ,j-1
r i-/2 r j-1r


After algebraic manipulation

-ai j+ bimj -cii+,j = rhsi (56a)


1/2 (u +jul ) Pr
ai Z i+l-Zi-l zi+lzi-ilzi-zi-)

bi E z -z (un+iui-un+tunI) + rz 1- +z 1
1 A-T zi+ i-1R R L i+- i-1 +i+1 1i

1/2 (uR- R ) Pr
c. +
1 --z -
ci z i+1-z-1 ii+1-z -1)(i+l-zi

(n+l n+1
rhs. = D2. +Ewm .+Fwn -0.5GrPr2 ij+l j
S ,j-1 1,i 3 i j+1 rj+1 -rj

1/2 (v + vjI) Pr rj-1
d = (v -I+ L J-
j+r -1 rj1/2 rj+1 rj_-1 j r-1

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 j-1 j+ -rj-1 2 rjl-r] rj

1/2 (VRVR Pr rj+1
f=1 +
rj+1-r1_ rj+1l-rj1 Jrj /2 rj+l-r )

After algebraic manipulation

-aj n+l b mn(l n+l
-a ,-1 + b i, + c ,j+ rhs (56b)


1/2 (vn+ vj ) Pr rj_1
aj= (r + r. r -r )Lr.-r. )
j+1 j-1 -1/2 j+ j j-1

(v n n n n
1 R+IVR -vL+Iv L + Pr rj 1 +
bj= A + 2(j+lrj-r (r -rj )Irj (rj+r-rj rj/2 1 rj-1r )
j+1 j-1 j+1 1

/2 (VR R Pr rj+
c.= -r +
J r -rj+l rj-r.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 1-1)

z i+1-z -z. zi+1-zi-1)(z-z
i+1 i-i LI) 1-1

1/2 (-n URI Pr
f =-I +
Zi+l-Z i-I zi+l-zi Zi+l-i)

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


b1 -c1 01,j rhsl+alw ,j

-a2 b2 -C2 "2,j rhs2
- -a3 b3 -c3 m3,j rhs3

-am-2 bm-2 -cm-2 "m-2,j rhsm-2
-am_-1 ,bm- m-l,j rhsm-l+am-l om,j


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.
b1 -cl i,1 rhsl+al wi,O

-a2 b2 -C2 "i,2 rhs2
-a3 b3 3 ,3 rhs3

-an-2 bn-2 -cn-2 "i,n-2 rhsn-2
-an-1 bn i,n-1 rhSn-l+cn-I "in


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/2-z 1(rj rj

-v t v t t u
S(vj i-1,j)(j+/2-r j 1/) (u j-uij-1)(zi+1i/2-z1/

t t t t
= 1i,j- i- ,j i,3 i',j-1 (59)
ij zi+/2-z 1/2 r.+ 2- r.
i + W zi 1/21/2


- -- Zi-/2

I zj-I
r]. l

Notations of various quantities defined on a
non-uniform grid system in cylindrical




I -1

r i./2 r] r

Figure 4.

iJ Z---+Z/2

u --


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


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,j-i-,j
I i+ Zi i-

1 ri j+l- ij 1iji,j-1
rI r )I-
rj+ 1/2-r 1/2 rj+ 1/2rj+1-r ) r.1/ -rj -1


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,j-1 + Eij+1



B (z -zi j z i i-i rj { (rj+l-rj) jr jr -j-l

+ 1 -1/22-+ 1 '-/1
C+ r-.zi.1

C [zi+ /2-zi_ 1/2 )Zi+-i )

D (rj+ 1/2-rj. 1/2 Jrj. 1/2 [r-rj-1

E =1
(r+ 1/2-r 1/2 rj+ 1/2 rj+1-r

The grid point velocities are

interpolation of the transport veloc

then achieved
ities u. j and


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
i ,j



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 .,j-1 ,f i, 1- i-1,j (64a,b)
1, r ( 1/2j 1,j j z 1/2i zi-1/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

Lt u u = 0
Ar+0 i,

Lt v v = 0

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)


j + xj = 1 (67)

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


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 vorticity-stream 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 over-relaxation 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, round-off errors can

be incurred at each mathematical operation, and simply accumulate

until final answers are obtained. When iterative techniques are

used, the presence of round-off errors at the end of any given

iteration simply results in those unknowns being somewhat poorer

estimates for the next iteration. For practical purposes, round-off

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 non-rectangular region may get

complicated, and they are not known for sure to be faster (Roache,


Another popular method to solve the unsteady formulation of the

stream function Poisson equation is the successive over-relaxation

(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 i-l ,j+Ci+,j+Dij-l+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 10-5. 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

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)

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

= 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 no-slip (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 stream-function

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 time-step 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 time-step

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,

SO,j- v,j
O,j (z l zO)

t t
UO,j- UO,j-1
( rj+ 1/-rj 12

- 1,j
r.z z12


uj= 0

_ ( ,j-Oj)
rj(z1- z0)

S0 from eqn 60a
and uj-=j 0 from eqn. 60a

At the top boundary,

SVmj- Vm-lj
m,j= (z z

t t
Smj- um,j-1
( rj/- r.j

_- m-l,j
r. (z -z z -z 1
J m m-1 m m-2


v .= 0

m- ,j z -
S m ,jm-1,j
Vm-l,j-jr rn(Z-mI)

- im-1,j
- r.(z -z )
I m m-1

utj and u t ,- 0 from eqn. 60a
mij mll=,


r z1


At the centerline,
t t t
i,0- Vi-l,0 ui,0- Ui,
i,0o- z. i+- z.1, r i/ 0

=0 (76c)


t O and v-, from eqn. 60b
1,0 1a1,0

t i 4,I- -i,0 2_ 1
io rl- O) rli


At the sidewall,
t t t
S,n Vi-l,n i,n- i,n-
i,n= (z+l/2- z l/ (rn- rnl/

= in-1 (76d)
r ni(r -r )(r r
n n n-1 n n


Vi, and vi 0 from eqn. 60b
1i ,n 1i- l,n=

t i,n i,n-1 1,n-1
1i,n-1 r/rn- r r r n -rl)
n -/2 n n-1 n_/2 n n-1

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)


n = wall vorticity value at time step n

n+1 = most recent estimate of wall vorticity at n+1
step using equation.

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

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
i j
together with a convergence criterion on the calculated Nusselt


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 10-4 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

non-dimensional 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


*n n"

ne m 1.0


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

W )

en + 1

n /2mo
wD = 0

2n1/2' f r n n4
n +3
w = f,(w, u
n+lm I
U) =O
n+1" f n 1/2

^ n+lm
) = Wr



, Q

, p 2

n+1m n+l1
uS vS

n1/2m n+1m
4 w 'a


f5 (n+1m
fo 0

S 7 j,

= n+

S,, n+
+ (1-y)w


(u,v) l

r: solid boundary
Q: interior points
: centerline
I ^ *
m,w,w: temporary

( n+1m
-f r~l!


Figure 5. Algorithm for the iterative solution
of the set of equations



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


4. Boundary vorticity values ( ab) at solid walls (top,

bottom and side) are first calculated-using 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

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)

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 under-relaxation 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 under-relaxation, 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


Modeling of Conduction Heating in Cans

The partial differential equation (PDE) governing conduction

heating of a cylindrical-shaped solid was discretized and finite

1W. E. Dunn, Mechanical and Industrial Engineering Department,
University of Illinois at Urbana-Champaign. 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

real-time for on-line 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

temperature-specified 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 conduction-heated 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

conduction-heating food model system (Niekamp et al., 1984). Eight

cans were fitted with an Ecklund-type copper-constantan thermocouple

mounted mid-way 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 head-space 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 steam-cook 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

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs