PRESSURE AND VELOCITY DISTRIBUTION FOR AIR FLOW
THROUGH FRUITS PACKED IN SHIPPING CONTAINERS
USING POROUS MEDIA FLOW ANALYSIS
By
MICHAEL THOMAS TALBOT
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF
THE UNIVERSITY OF FLORIDA
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1987
u~a: i
ACKNOWLEDGMENT
The author wishes to express his sincere appreciation
to Dr. Calvin C. Oliver who served as chairman of his
supervisory committee. His guidance, encouragement, and
assistance were invaluable throughout the course of this
study.
For the leadership, interest, and assistance, they
exhibited, appreciation is extended to the other members of
his supervisory committee: Mr. J. J. Gaffney, USDA
Agricultural Engineer; Dr. H. A. Ingley III, Dr. T. I. Shih,
Mechanical Engineering; and Dr. R. F. Matthews, Food Science
and Human Nutrition. Special thanks are also expressed for
the early contributions of Dr. R. K. Irey, an original
committee member who is presently Chairman of the Mechanical
Engineering Department at the University of Toledo, Ohio.
Appreciation is expressed to Dr. R. B. Gaither,
Chairman of the Mechanical Engineering Department, and to
the faculty members who contributed to his graduate program.
A special word of thanks is necessary for the
invaluable assistance of Dean M. M. Lockhart, Dean J. L.
Woeste and Dean G. L. Zachariah, as well as the Graduate
Council, without which pursuit of this course of study would
not have been possible.
He is indebted to the Department of Agricultural
Engineering for financial support through fulltime
employment during the parttime pursuit of this course of
study and to Dr. G. W. Isaacs, chairman, for his leadership,
encouragement, and untiring willingness to cope with
administrative problems involved. Thanks are also due to
each of the many staff and faculty members who assisted in
numerous ways, especially to Dr. C. D. Baird and Dr. K. V.
Chau who actively supported his research as project leaders.
For their cheerful sacrifices and loving patience, the
author dedicates this dissertation to his wife, Anita, his
son, Michael, and his father, James Vernon (Buck) Talbot,
who waited so eagerly and proudly but who could not be here
for the completion.
Most of all, the author wishes to acknowledge the
source of all knowledge and wisdomthe very God of the
universe through His Son the Lord Jesus Christ. To God be
the glory, for great things He hath done.
TABLE OF CONTENTS
Page
ACKNOWLEDGMENT ........................................ ii
LIST OF TABLES ........................................ vi
LIST OF FIGURES .................... ... ................ vii
ABSTRACT .......... ....................................xxvii
INTRODUCTION .......................................... 1
OBJECTIVES ............................................ 5
REVIEW OF LITERATURE ................................... 6
Background of Research .................... ........ .... 6
Fluid Flow Through Porous Media ..................... 12
PROCEDURE ........................................... 32
Commercial Finite Element Package.................... 43
Verification ........................................ 44
TwoDimensional Grain Bin ......................... 49
ThreeDimensional Grain Bin ....................... 52
Solution for Oranges Packed in Cartons .............. 57
ThreeDimensional Bulk Oranges .................... 60
ThreeDimensional Orange Carton ................... 78
Experimental Temperature Measurement ................ 88
Description of Experimental Facilities ............ 88
Basic Construction of Precooler ................... 88
Refrigeration Components and Controls ............. 89
Air Temperature and Humidity ...................... 90
Air Flow .......................................... 95
Reheat Section .................................... 97
Experimental Test Procedures ...................... 98
Heat Transfer Model ................................... 101
Selection of Model .................................. 102
Convective Heat Transfer Coefficient .............. 112
Heat Transfer Model Modifications ................. 116
ANSYS Output Data Reorganization .................. 124
Heat Transfer Program ...... ....................... 126
Temperature Response Data Reorganization ......... 128
RESULTS AND DISCUSSION ............................... 131
TwoDimensional Grain Bin............................ 131
ThreeDimensional Grain Bin........................... 139
ThreeDimensional Bulk Oranges....................... 154
ThreeDimensional Orange Carton.................... 161
Experimental Temperature Measurement................. 165
Heat Transfer Model.................................... 175
Experimental Error .................................... 179
Convective Heat Transfer Coefficient................. 181
Variable Porosity................................... 182
Approximate Model Thermocouple Location.............. 183
Experimental Versus Predicted Temperature Response... 183
SUMMARY...... ........ .................................. 197
SUGGESTIONS FOR FUTURE WORK............................ 200
APPLICATION OF RESULTS ................................. 202
CONCLUSIONS...... ........ .............................. 204
APPENDIX A METHOD FOR ACCURATELY POSITIONING
THERMOCOUPLES IN FRUITS AND VEGETABLES................. 206
APPENDIX B FORTRAN PROGRAMS AND DATA FILES............ 210
APPENDIX C PLOT INFORMATION........................... 234
APPENDIX D ADDITIONAL DISCUSSION OF RESULTS FOR TWO
AND THREEDIMENSIONAL GRAIN BIN MODELS................. 236
APPENDIX E ADDITIONAL EXPERIMENTAL TEMPERATURE RESPONSE
GRAPHS................................................. 290
APPENDIX F ADDITIONAL PREDICTED TEMPERATURE RESPONSE
GRAPHS................................................. 324
APPENDIX G ADDITIONAL PREDICTED VERSUS EXPERIMENTAL
TEMPERATURE RESPONSE GRAPHS............................ 358
APPENDIX H ADDITIONAL PREDICTED VERSUS EXPERIMENTAL
TEMPERATURE REGRESSION GRAPHS.......................... 425
BIBLIOGRAPHY..... ......... ............................. 437
BIOGRAPHICAL SKETCH. .................................. 449
LIST OF TABLES
Table Page
1. Summary of variable porosity and ANSYS input data. 85
2. Inlet and outlet air vent locations as indicated
by the numbers illustrated in Figure 11........... 87
3. Boundary conditions and air flow rates (2 flow
rates per venting arrangement) for tests
conducted with an experimental orange carton
packed with size 100 oranges ..................... 87
4. Physical and thermal properties of size 100
Valencia oranges, Gaffney and Baird (1980)........ 114
5. Pressure values (inches of water) at specific
points for the linear flow problem............... 132
6. Threedimensional pressure loss data for oranges
in bulk and orange carton sides................... 155
LIST OF FIGURES
Figure Page
1. Region divided into finite element................ 36
2. Various types of finite elements.................. 39
3. The node and element grid employed by Segerlind
(1982) ............................................. 50
4. Boundary conditions used by Segerlind (1982)...... 51
5. The node and element locations used in current
study.............................................. 53
6. Four different grainbin aeration systems studied
by Khompis et al. (1984) .......................... 55
7 Finite element model used by Khompis et al.(1984). 56
8. Tne node and element locations used in current
study.............................................. 58
9. The experimental setup used by Chau et al. (1983)
or measuring air flow resistance in orange......... 64
10. The node and element locations used to model 3D
oranges packed in bulk and in cartons ............. 71
11. The experimental orange carton..................... 80
12. The node and element locations used to model a
3D orange carton.................................. 81
13. Forcedair precooler components and air
circulation diagram, Baird et al. (1975) .......... 91
14. Modification of product bin to accommodate
experimental carton............................... 92
15. Precooler refrigeration components and controls,
Baird et al. (1975) ............................... 93
16. Fruit and thermocouple locations in experimental
carton............................................. 100
17. Diagram indicating a heat balance on the three
types of nodes used in a sphere, Baird and Gaffney
(1976) ............. ... .................... 104
18. Flow diagram for the computer program for heat
conduction in a sphere with arbitrarily varying
air temperature, Baird and Gaffney (1976)......... 107
19. Control volume for heat balance on bed, Baird and
Gaffney (1976) ................................... 109
20. Flow diagram for the computer program to solve
deep bed heat transfer equations, Baird and
Gaffney (1976) .................................. 113
21. Typical threedimensional element showing velocity
at centroid and air flow across element faces.... 118
22. Flow diagram for the computer program to solve
orange carton heat transfer equations............. 123
23. Points of comparison for linear airflow
calculations ..................................... 132
24. Plot of experimental isobaric lines from Brooker
(1969) ............................................ 133
25. 85 percent isobaric contour line from Segerlind
(1982) ............................................. 134
26. Isobaric contour lines calculated using ANSYS
for a boundary condition of 1 inch of water...... 135
27. Isobaric contour lines calculated using ANSYS
for a boundary condition of 2 inches of water.... 136
28. Isobaric contour lines calculated using ANSYS
for a boundary condition of 3 inches of water.... 137
29. Isobaric distribution in percent of total
pressure loss on a vertical cutting plane from
Khompis (1983) .................................. 141
30. Isobaric distribution in percent of total
pressure loss on horizontal cutting plane at 5
percent of bin height from Khompis (1983)......... 142
31. Isobaric distribution calculated using ANSYS in
percent of total pressure loss on the vertical
cutting plane.................................... 143
32. Isobaric distribution calculated using ANSYS in
percent of total pressure loss on the horizontal
cutting plane at 5 percent of bin height.......... 144
viii
33. Velocity distribution in meters per minute on a
vertical cutting plane from Khompis (1983) ....... 145
34. Velocity distribution in meters per minute on a
horizontal cutting plane at 5 percent of bin
height from Khompis (1983) ....................... 146
35. Velocity distribution calculated using ANSYS in
meters per minute on a verical cutting plane..... 147
36. Velocity distribution calculated using ANSYS in
meters per minute on a horizontal cutting plane
at 5 percent of bin height.... ................... 148
37. 70 percent of total pressure loss isobaric
surface from Khompis (1983), view 2.............. 149
38. 70 percent of total pressure loss isobaric
surface produced by ANSYS, view 2................ 150
39. Constant velocity surface of 2.4 meters per
minute from Khompis (1983), view 2................ 152
40. Constant velocity surface of 2.8 meters per
minute calculated using ANSYS, view 2............. 153
41. Nodal pressure printout for ANSYS solution........ 162
42. Element flow rate printout for ANSYS solution.... 163
43. Example of the data printout produced by the data
acquisition system for one experimental cooling
test............................................. 166
44. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 11,
thermocouple locations 10, 13, 66, and 69)....... 168
45. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 11,
3rd layer thermocouple locations 31, 35, 36, 44,
and 45) .......................................... 171
46. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 11,
2nd layer thermocouple locations 23, 24, 26, 27,
and 29).......................................... 172
47. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 11,
1st layer thermocouple locations 1, 9, 10, 14,
and 15)........................................... 173
48. Graphical representation of threedimensional
color graphics outpur of cooling response for
oranges within an experimental carton............. 174
49. Predicted center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 11,
3rd layer thermocouple locations 31, 35, 36, 44,
and 45) ........................................... 176
50. Predicted center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 11,
2nd layer thermocouple locations 23, 24, 26, 27,
and 29).......................................... 177
51. Predicted center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 11,
1st layer thermocouple locations 1, 9, 10, 14,
and 15)........................................... 178
52. Predicted versus experimental data for Boundary
Condition 11, 3rd layer thermocouples 31 and
35................................................. 185
53. Predicted versus experimental data for Boundary
Condition 11, 3rd layer thermocouples 44 and
45. ............................................... 186
54. Predicted versus experimental data for Boundary
Condition 11, 2nd layer thermocouples 23 and
24................................................. 187
55. Predicted versus experimental data for Boundary
Condition 11, 2nd layer thermocouples 26 and
29................................................. 188
56. Predicted versus experimental data for Boundary
Condition 11, 1st layer thermocouples 1 and 9... 189
57. Predicted versus experimental data for Boundary
Condition 11, 1st layer thermocouples 14 and
15................................................. 190
58. Regression plot of the predicted temperature
versus the experimental temperature for Boundary
Condition 11 and test times of 3, 7, and 11
hours ............................................ 191
Ai. Thermocouple positioning, part 1.................. 208
A2. Thermocouple positioning, part 2.................. 209
B1. Portion of input data file "Dimensions" which
lists for each element from left to right, the
number, volume, x, y, and z dimensions, and the
porosity ......................................... 211
B2. Portion of input data file "Data" which lists
for each element from left to right, the number,
total velocity and x, y, and z component
velocities....................................... 212
B3. Fortran program "DataReduction" which was used
to calculate convective heat transfer coefficient
and x, y, and z component mass flow rates at the
centroid.......................................... 213
B4. Portion of input/output data file "HeatMass"
which lists for each element, the number, the
convective heat transfer coefficient and x, y,
and z component mass flow rates at the centroid.. 215
B5. Fortran program "FlowWall" which was used to
calculate x, y, and z component mass flow rates
across the face of each element.................. 216
B6. Portion of input/output data file "MassFlow"
which lists for each element, the number, the x,
y, and z component mass flow rates across the
face of each element............................. 218
B7. Fortran program "TempID" which was used to
calculate mass flow rates entering each element.. 219
B8. Portion of input/output data file "MassTemp"
which lists for each element, the number, the
mass flow rates entering the x (left and right),
y (bottom and top), and z (back and front) faces
of each element along with the corresponding
temperature identification label.................. 222
B9. Fortran program "HeatDo" which was used to solve
the heat transfer model .......................... 223
B10. Input data file "DataIn1" which lists the
product thermal diffusivity, conductivity, and
radius, the specific heat of the air and product
and product density............................... 229
Bll. Input data file "DataIn2" which lists number
of product nodal points, the initial product
temperature, the time increment, total test
cooling time, data output print frequency,
entering air temperature, and number of elements. 229
B12. Portion of input/output data file "Temp
Response" which lists for each element, the
number, time (seconds),the air temperature, the
product temperature at each of the nodes, and the
product massaverage temperature................. 230
B13. Fortran program "PlotData" which was used to
reorganize the temperature response to give 78
thermocouple readings............................ 231
Dl. Numerical pressure pattern for linear air flow
through shelled corn from Brooker (1969), for 2
inches of water duct pressure.................... 237
D2. Numerical pressure pattern for nonlinear air
flow through shelled corn from Brooker (1969).
The same pattern was obtained for duct pressures
of 1, 2, and 3 inches of water.................... 238
D3. Numerical velocity distribution for nonlinear
air flow through shelled corn from Brooker
(1969) ........................................... 239
D4. Plot of experimental isobaric lines from
Brooker (1969) .................................. 239
D5. 85 percent isobaric contour line from Segerlind
(1982) ............................................ 240
D6. Isobaric contour lines calculated using ANSYS
for a boundary condition of 1 inch of water...... 242
D7. Isobaric contour lines calculated using ANSYS
for a boundary condition of 2 inches of water.... 243
D8. Isobaric contour lines calculated using ANSYS
for a boundary condition of 3 inches of water.... 244
D9. Modified boundary conditions used by Segerlind
(1982) ........................................... 246
D10. 85 percent isobaric lines from Segerlind (1982).. 246
D11. Isobaric pattern for a duct pressure of 3 inches
of water from Segerlind (1982) ................... 247
D12. Isobaric contour lines calculated using ANSYS
with a modified boundary condition of 1 inch of
water ............................................. 248
D13. Isobaric contour lines calculated using ANSYS
with a modified boundary condition of 2 inches
of water ......................................... 249
D14. Isobaric contour lines calculated using ANSYS
with a modified boundary condition of 3 inches
of water ......................................... 250
D15. Vertical and Horizontal Cuttings Planes used
by Khompis (1983) ................................ 252
D16. Location of the elements on the vertical
cutting plane from Khompis (1983) ................ 253
D17. Location of the elements on the horizontal
cutting plane from Khompis (1983) ................ 254
D18. Isobaric distribution in percent of total
pressure loss on a vertical cutting plane from
Khompis (1983) .................................. 256
D19. Isobaric distribution in percent of total
pressure loss on horizontal cutting plane at 5
percent of bin height from Khompis (1983)......... 257
D20. Isobaric distribution in percent of total
pressure loss on a horizontal cutting plane at
50 percent of bin height from Khompis (1983)..... 258
D21. 70 percent of total pressure loss isobaric
surface from Khompis (1983) ..................... 260
D22. 70 percent of total pressure loss isobaric
surface from Khompis (1983)....................... 261
D23. 70 percent of total pressure loss isobaric
surface from Khompis (1983) ..................... 262
D24. Velocity distribution in meters per minute on
a vertical cutting plane from Khompis (1983)..... 263
D25. Velocity distribution in meters per minute on
a horizontal cutting plane at 5 percent of bin
height from Khompis (1983) ....................... 264
xiii
D26. Velocity distribution in meters per minute on
a vertical cutting plane at 50 percent of bin
height from Khompis (1983)........................ 265
D27. Minimum velocity regions for the square duct
system from Khompis (1983) ....................... 266
D28. Constant velocity surface of 2.4 meters per
minute from Khompis (1983) ....................... 267
D29. Constant velocity surface of 2.4 meters per
minute from Khompis (1983)........................ 268
D30. Constant velocity surface of 2.4 meters per
minute from Khompis (1983) ....................... 269
D31. Isobaric distribution calculated using ANSYS in
percent of total pressure loss on the vertical
cutting plane.................................... 271
D32. Isobaric distribution calculated using ANSYS in
percent of total pressure loss on the horizontal
cutting plane at 5 percent of bin height.......... 272
D33. Isobaric distribution calculated using ANSYS in
percent of total pressure loss on the horizontal
cutting plane at 50 percent of bin height......... 273
D34. 70 percent of total pressure loss isobaric
surface produced by ANSYS, view 1................. 276
D35. 70 percent of total pressure loss isobaric
surface produced by ANSYS, view 2................. 277
D36. 70 percent of total pressure loss isobaric
surface produced by ANSYS, view 3................. 278
D37. Velocity distribution calculated using ANSYS in
meters per minute on a verical cutting plane..... 281
D38. Velocity distribution calculated using ANSYS in
meters per minute on a horizontal cutting plane
at 5 percent of bin height....................... 282
D39. Velocity distribution calculated using ANSYS in
meters per minute on a horizontal cutting plane
at 50 percent of bin height ...................... 283
D40. Constant velocity surface of 2.8 meters per
minute calculated using ANSYS, view 1............. 287
D41. Constant velocity surface of 2.8 meters per
minute calculated using ANSYS, view 2............. 288
D42. Constant velocity surface of 2.8 meters per
minute calculated using ANSYS, view 3............ 289
El. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 12, 3rd
layer thermocouple locations 31, 35, 36, 44, and
45)...................................... ......... 291
E2. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Conurtion 12, 2nd
layer thermocouple locations 23, 24, 26, 27, and
29)............................................... 292
E3. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 12, 1st
layer thermocouple locations 1, 9, 10, 14, and
15)... ........................................... 293
E4. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 21, 3rd
layer thermocouple locations 31, 35, 36, 44, and
45).............................................. 294
E5. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 21, 2nd
layer thermocouple locations 23, 24, 26, 27, and
29)... ........................................... 295
E6. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 21, 1st
layer thermocouple locations 1, 9, 10, 14, and
15).... ........................................... 296
E7. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 22, 3rd
layer thermocouple locations 31, 35, 36, 44, and
45) ............................................. 297
E8. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 22, 2nd
layer thermocouple locations 23, 24, 26, 27, and
29)...................................... ......... 298
E9. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 22, 1st
layer thermocouple locations 1, 9, 10, 14, and
15).... ........................................... 299
E10. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 31, 3rd
layer thermocouple locations 31, 35, 36, 44, and
45)... ........................................... 300
Ell. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 31, 2nd
layer thermocouple locations 23, 24, 26, 27, and
29)... ........................................... 301
E12. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 31, 1st
layer thermocouple locations 1, 9, 10, 14, and
15)... ........................................... 302
E13. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 32, 3rd
layer thermocouple locations 31, 35, 36, 44, and
45)... ........................................... 303
E14. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 32, 2nd
layer thermocouple locations 23, 24, 26, 27, and
29)... ........................................... 304
E15. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 32, 1st
layer thermocouple locations 1, 9, 10, 14, and
15)...................................... ......... 305
E16. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 41, 3rd
layer thermocouple locations 31, 35, 36, 44, and
45).... ........................................... 306
E17. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 41, 2nd
layer thermocouple locations 23, 24, 26, 27, and
29).... ........................................... 307
xvi
E18. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 41, 1st
layer thermocouple locations 1, 9, 10, 14, and
15)..... .................... ...................... 308
E19. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 42, 3rd
layer thermocouple locations 31, 35, 36, 44, and
45).............................................. 309
E20. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 42, 2nd
layer thermocouple locations 23, 24, 26, 27, and
29)............................................... 310
E21. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 42, 1st
layer thermocouple locations 1, 9, 10, 14, and
15).............................................. 311
E22. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 51, 3rd
layer thermocouple locations 31, 35, 36, 44, and
45).............................................. 312
E23. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 51, 2nd
layer thermocouple locations 23, 24, 26, 27, and
29)... ........................................... 313
E24. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 51, 1st
layer thermocouple locations 1, 9, 10, 14, and
15).............................................. 314
E25. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 52, 3rd
layer thermocouple locations 31, 35, 36, 44, and
45)... ........................................... 315
E26. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 52, 2nd
layer thermocouple locations 23, 24, 26, 27, and
29)... ........................................... 316
xvii
E27. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 52, Ist
layer thermocouple locations 1, 9, 10, 14, and
15)...................................... ......... 317
E28. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 61, 3rd
layer thermocouple locations 31, 35, 36, 44, and
45)... ........................................... 318
E29. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 61, 2nd
layer thermocouple locations 23, 24, 26, 27, and
29).............................................. 319
E30. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 61, 1st
layer thermocouple locations 1, 9, 10, 14, and
15)... ........................................... 320
E31. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 62, 3rd
layer thermocouple locations 31, 35, 36, 44, and
45).... ........................................... 321
E32. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 62, 2nd
layer thermocouple locations 23, 24, 26, 27, and
29)...................................... ......... 322
E33. Experimental center temperature response of
oranges in an experimental orange carton for
forcedair cooling (Boundary Condition 62, 1st
layer thermocouple locations 1, 9, 10, 14, and
15).... .............................. ............. 323
Fl. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 12, 3rd layer
thermocouple locations 31, 35, 36, 44, and 45)... 325
F2. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 12, 2nd layer
thermocouple locations 23, 24, 26, 27, and 29)... 326
xviii
F3. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 12, Ist layer
thermocouple locations 1, 9, 10, 14, and 15) ..... 327
F4. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 21, 3rd layer
thermocouple locations 31, 35, 36, 44, and 45)... 328
F5. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 21, 2nd layer
thermocouple locations 23, 24, 26, 27, and 29)... 329
F6. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 21, Ist layer
thermocouple locations 1, 9, 10, 14, and 15) ..... 330
F7. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 22, 3rd layer
thermocouple locations 31, 35, 36, 44, and 45)... 331
F8. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 22, 2nd layer
thermocouple locations 23, 24, 26, 27, and 29)... 332
F9. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 22, 1st layer
thermocouple locations 1, 9, 10, 14, and 15)..... 333
F10. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 31, 3rd layer
thermocouple locations 31, 35, 36, 44, and 45)... 334
Fll. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 31, 2nd layer
thermocouple locations 23, 24, 26, 27, and 29)... 335
F12. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 31, 1st layer
thermocouple locations 1, 9, 10, 14, and 15) ..... 336
F13. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 32, 3rd layer
thermocouple locations 31, 35, 36, 44, and 45)... 337
F14. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 32, 2nd layer
thermocouple locations 23, 24, 26, 27, and 29).. 338
F15. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 32, 1st layer
thermocouple locations 1, 9, 10, 14, and 15).... 339
F16. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 41, 3rd layer
thermocouple locations 31, 35, 36, 44, and 45).. 340
F17. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 41, 2nd layer
thermocouple locations 23, 24, 26, 27, and 29)... 341
F18. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 41, 1st layer
thermocouple locations 1, 9, 10, 14, and 15)..... 342
F19. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 42, 3rd layer
thermocouple locations 31, 35, 36, 44, and 45)... 343
F20. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 42, 2nd layer
thermocouple locations 23, 24, 26, 27, and 29)... 344
F21. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 42, Ist layer
thermocouple locations 1, 9, 10, 14, and 15)..... 345
F22. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 51, 3rd layer
thermocouple locations 31, 35, 36, 44, and 45)... 346
F23. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 51, 2nd layer
thermocouple locations 23, 24, 26, 27, and 29)... 347
F24. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 51, Ist layer
thermocouple locations 1, 9, 10, 14, and 15)..... 348
F25. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 52, 3rd layer
thermocouple locations 31, 35, 36, 44, and 45)... 349
F26. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 52, 2nd layer
thermocouple locations 23, 24, 26, 27, and 29)... 350
F27. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 52, 1st layer
thermocouple locations 1, 9, 10, 14, and 15) ..... 351
F28. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 61, 3rd layer
thermocouple locations 31, 35, 36, 44, and 45)... 352
F29. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 61, 2nd layer
thermocouple locations 23, 24, 26, 27, and 29)... 353
F30. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 61, 1st layer
thermocouple locations 1, 9, 10, 14, and 15) ..... 354
F31. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 62, 3rd layer
thermocouple locations 31, 35, 36, 44, and 45)... 355
F32. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 62, 2nd layer
thermocouple locations 23, 24, 26, 27, and 29)... 356
F33. Predicted center temperature response of oranges
in an experimental orange carton for forcedair
cooling (Boundary Condition 62, 1st layer
thermocouple locations 1, 9, 10, 14, and 15) ..... 357
Gl. Predicted versus experimental data for Boundary
Condition 12, 3rd layer thermocouples 31 and 35. 359
G2. Predicted versus experimental data for Boundary
Condition 12, 3rd layer thermocouples 44 and 45. 360
G3. Predicted versus experimental data for Boundary
Condition 12, 2nd layer thermocouples 23 and 24. 361
G4. Predicted
Condition
G5. Predicted
Condition
G6. Predicted
Condition
G7. Predicted
Condition
G8. Predicted
Condition
G9. Predicted
Condition
G10. Predicted
Condition
Gll. Predicted
Condition
G12. Predicted
Condition
G13. Predicted
Condition
G14. Predicted
Condition
G15. Predicted
Condition
G16. Predicted
Condition
G17. Predicted
Condition
G18. Predicted
Condition
G19. Predicted
Condition
G20. Predicted
Condition
G21. Predicted
Condition
versus experimental data for
12, 2nd layer thermocouples
versus experimental data for
12, 1st layer thermocouples
versus experimental data for
12, 1st layer thermocouples
versus experimental data for
21, 3rd layer thermocouples
versus experimental data for
21, 3rd layer thermocouples
versus experimental data for
21, 2nd layer thermocouples
versus experimental data for
21, 2nd layer thermocouples
versus experimental data for
21, Ist layer thermocouples
versus experimental data for
21, 1st layer thermocouples
versus experimental data for
22, 3rd layer thermocouples
versus experimental data for
22, 3rd layer thermocouples
versus experimental data for
22, 2nd layer thermocouples
versus experimental data for
22, 2nd layer thermocouples
versus experimental data for
22, 1st layer thermocouples
versus experimental data for
22, 1st layer thermocouples
versus experimental data for
31, 3rd layer thermocouples
versus experimental data for
31, 3rd layer thermocouples
versus experimental data for
31, 2nd layer thermocouples
Boundary
26 and 29. 362
Boundary
1 and 9... 363
Boundary
14 and 15. 364
Boundary
31 and 35. 365
Boundary
44 and 45. 366
Boundary
23 and 24. 367
Boundary
26 and 29. 368
Boundary
1 and 9... 369
Boundary
14 and 15. 370
Boundary
31 and 35. 371
Boundary
44 and 45. 372
Boundary
23 and 24. 373
Boundary
26 and 29. 374
Boundary
1 and 9... 375
Boundary
14 and 15. 376
Boundary
31 and 35. 377
Boundary
44 and 45. 378
Boundary
23 and 24. 379
xxii
G22. Predicted versus experimental data for Boundary
Condition 31, 2nd layer thermocouples 26 and 29. 380
G23. Predicted versus experimental data for Boundary
Condition 31, 1st layer thermocouples 1 and 9...
381
G24. Predicted versus experimental data for
Condition 31, 1st layer thermocouples
G25. Predicted versus experimental data for
Condition 32, 3rd layer thermocouples
G26. Predicted versus experimental data for
Condition 32, 3rd layer thermocouples
Boundary
14 and 15. 382
Boundary
31 and 35. 383
Boundary
44 and 45. 384
G27. Predicted versus experimental data for Boundary
Condition 32, 2nd layer thermocouples 23 and 24.
G28. Predicted versus experimental data for
Condition 32, 2nd layer thermocouples
G29. Predicted versus experimental data for
Condition 32, 1st layer thermocouples
G30. Predicted versus experimental data for
Condition 32, 1st layer thermocouples
G31. Predicted versus experimental data for
Condition 41, 3rd layer thermocouples
Boundary
26 and 29.
Boundary
1 and 9...
Boundary
14 and 15.
Boundary
31 and 35.
G32. Predicted versus experimental data for Boundary
Condition 41, 3rd layer thermocouples 44 and 45.
G33. Predicted versus experimental data for Boundary
Condition 41, 2nd layer thermocouples 23 and 24.
G34. Predicted versus experimental data for
Condition 41, 2nd layer thermocouples
G35. Predicted versus experimental data for
Condition 41, 1st layer thermocouples
G36. Predicted versus experimental data for
Condition 41, 1st layer thermocouples
G37. Predicted versus experimental data for
Condition 42, 3rd layer thermocouples
G38. Predicted versus experimental data for
Condition 42, 3rd layer thermocouples
Boundary
26 and 29.
Boundary
1 and 9...
Boundary
14 and 15.
Boundary
31 and 35.
Boundary
44 and 45.
G39. Predicted versus experimental data for Boundary
Condition 42, 2nd layer thermocouples 23 and 24.
xxiii
G40. Predicted
Condition
G41. Predicted
Condition
G42. Predicted
Condition
G43. Predicted
Condition
G44. Predicted
Condition
G45. Predicted
Condition
G46. Predicted
Condition
G47. Predicted
Condition
G48. Predicted
Condition
G49. Predicted
Condition
G50. Predicted
Condition
G51. Predicted
Condition
G52. Predicted
Condition
G53. Predicted
Condition
G54. Predicted
Condition
G55. Predicted
Condition
G56. Predicted
Condition
G57. Predicted
Condition
versus experimental data for Boundary
42, 2nd layer thermocouples 26 and 29. 398
versus experimental data for Boundary
42, 1st layer thermocouples 1 and 9... 399
versus experimental data for Boundary
42, 1st layer thermocouples 14 and 15. 400
versus experimental data for Boundary
51, 3rd layer thermocouples 31 and 35. 401
versus experimental data for Boundary
51, 3rd layer thermocouples 44 and 45. 402
versus experimental data for Boundary
51, 2nd layer thermocouples 23 and 24. 403
versus experimental data for Boundary
51, 2nd layer thermocouples 26 and 29. 404
versus experimental data for Boundary
51, 1st layer thermocouples 1 and 9... 405
versus experimental data for Boundary
51, 1st layer thermocouples 14 and 15. 406
versus experimental data for Boundary
52, 3rd layer thermocouples 31 and 35. 407
versus experimental data for Boundary
52, 3rd layer thermocouples 44 and 45. 408
versus experimental data for Boundary
52, 2nd layer thermocouples 23 and 24. 409
versus experimental data for Boundary
52, 2nd layer thermocouples 26 and 29. 410
versus experimental data for Boundary
52, 1st layer thermocouples 1 and 9... 411
versus experimental data for Boundary
52, 1st layer thermocouples 14 and 15. 412
versus experimental data for Boundary
61, 3rd layer thermocouples 31 and 35. 413
versus experimental data for Boundary
61, 3rd layer thermocouples 44 and 45. 414
versus experimental data for Boundary
61, 2nd layer thermocouples 23 and 24. 415
xxiv
G58. Predicted versus experimental data for Boundary
Condition 61, 2nd layer thermocouples 26 and 29. 416
G59. Predicted versus experimental data for Boundary
Condition 61, 1st layer thermocouples 1 and 9... 417
G60. Predicted versus experimental data for Boundary
Condition 61, 1st layer thermocouples 14 and 15. 418
G61. Predicted versus experimental data for Boundary
Condition 62, 3rd layer thermocouples 31 and 35. 419
G62. Predicted versus experimental data for Boundary
Condition 62, 3rd layer thermocouples 44 and 45. 420
G63. Predicted versus experimental data for Boundary
Condition 62, 2nd layer thermocouples 23 and 24. 421
G64. Predicted versus experimental data for Boundary
Condition 62, 2nd layer thermocouples 26 and 29. 422
G65. Predicted versus experimental data for Boundary
Condition 62, 1st layer thermocouples 1 and 9... 423
G66. Predicted versus experimental data for Boundary
Condition 62, 1st layer thermocouples 14 and 15. 424
H1. Regression plot of the predicted temperature
versus the experimental temperature for Boundary
Condition 12 and test times of 0.75, 1.75, and
2.75 hours ....................................... 426
H2. Regression plot of the predicted temperature
versus the experimental temperature for Boundary
Condition 21 and test times of 2, 5, and 8
hours ............................................ 427
H3. Regression plot of the predicted temperature
versus the experimental temperature for Boundary
Condition 22 and test times of 1, 4, and 7
hours ........................................... 428
H4. Regression plot of the predicted temperature
versus the experimental temperature for Boundary
Condition 31 and test times of 2, 4, and 6
hours . .......................................... 429
H5. Regression plot of the predicted temperature
versus the experimental temperature for Boundary
Condition 32 and test times of 0.4, 1.2, and 2.0
hours ........................................... 430
xxv
H6. Regression plot of the predicted temperature
versus the experimental temperature for Boundary
Condition 41 and test times of 1, 4, and 7
hours ........................................... 431
H7. Regression plot of the predicted temperature
versus the experimental temperature for Boundary
Condition 42 and test times of 1, 2.5, and 4
hours ............................................ 432
H8. Regression plot of the predicted temperature
versus the experimental temperature for Boundary
Condition 51 and test times of 2, 5, and 8
hours ............................................ 433
H9. Regression plot of the predicted temperature
versus the experimental temperature for Boundary
Condition 52 and test times of 1, 2.5, and 4
hours ............................................ 434
H10. Regression plot of the predicted temperature
versus the experimental temperature for Boundary
Condition 61 and test times of 2, 6, and 10
hours ............................................ 435
Hll. Regression plot of the predicted temperature
versus the experimental temperature for Boundary
Condition 62 and test times of 1, 2, and 5
hours ............................................ 436
xxvi
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
PRESSURE AND VELOCITY DISTRIBUTION FOR AIR FLOW
THROUGH FRUITS PACKED IN SHIPPING CONTAINERS
USING POROUS MEDIA FLOW ANALYSIS
By
Michael Thomas Talbot
August, 1987
Chairman: Dr. C. C. Oliver
Major Department: Mechanical Engineering
Detailed pressure and velocity fields for low air flow
rates used for forced cooling of oranges packed in
containers cannot be established experimentally. The
superficial velocity and overall pressure drop can be
determined for oranges packed in bulk and the temperature
response during cooling of oranges packed in containers can
be measured.
Flow analysis requires the consideration of at least a
second order velocity term to account for friction and
dynamic effects. The theory of fluid flow through porous
media incorporates such an analysis and has been
successfully applied to determine the pressure and velocity
fields for air flow through semiinfinite bulk agricultural
products. The objective of this study was to evaluate the
xxvii
applicability of porous media flow theory for pressure and
velocity flow field predictions for the finite boundary
conditions of air flow through oranges packed in shipping
containers.
A commercial finite element solution package was
adapted to determine the pressure and velocity distribution
and was verified by comparison with existing detailed
studies of bulk systems. Porous media input parameters and
boundary condition specifications for the current study were
obtained from studies of pressure loss through oranges in
bulk and in simulated orange cartons. An existing heat
transfer program was modified to incorporate predicted
velocity inputs and calculate temperature response.
Twelve tests were conducted using an experimental
orange carton with 88 Valencia oranges (size 100) packed in
a facecentered cubic packing arrangement using air flow
3 2 3
rates ranging from 1.5 x 10 to 2.0 x 10 m /s. The tests
involved six different venting arrangements. Temperature
responses were measured for 58 orange centers, 10 orange
surfaces, and 10 air spaces.
The measured temperature responses for the twelve tests
were compared to the predicted temperature responses and
provided indirect evaluation of the input flow predictions.
Variable porosity to account for edge effects was used to
improve initial results.
Although several areas for improvement were noted, the
porous media flow analysis was found to provide necessary
xxviii
flow field information to calculate the thermal response for
various flow boundary conditions.
xxix
INTRODUCTION
Losses of fresh fruits and vegetables and other
horticultural commodities from decay and shriveling as a
result of poor temperature management during postharvest
handling, transportation, and marketing are substantial in
the United States, as well as in other production and
consumption areas of the world.
Reduction of such losses is perhaps of more immediate
importance than increasing yields through improved
production practices. This is especially apparent when one
considers that all of the cost and energy inputs for
production, harvesting, packing, transportation, etc., are
wasted for that portion of the supply for which losses occur
near the point of consumption. For most fruits and
vegetables, the total cost for harvesting, packing,
transportation, storage, selling, etc., greatly exceeds the
cost of production. Therefore, the direct costs which can
be attributed to losses during marketing are magnified
greatly when compared with the alternative of increased
production to effect increased supply.
Respiration is the overall process by which stored
organic materials (carbohydrates, proteins, fats) are broken
into simple end products with a release of energy. Oxygen
(02) is used in this process and carbon dioxide (C02) is
produced by the commodity. The loss of the stored organic
materials from the commodity during respiration means
hastened senescence as the energy for maintaining the living
status of the commodity is exhausted, loss of food value
(energy value) for the consumer, reduced flavor quality,
especially sweetness, and loss of salable dry weight. The
rate of deterioration (perishability) of harvested
commodities is generally proportional to their respiration
rate. Product temperature is a major determinant of the
rate of respiratory activity.
Water loss can be one of the main causes of
deterioration since it results in not only direct
quantitative losses, but also causes losses in appearance
(due to wilting and shriveling), textural quality
(softening, flaccidity, limpness, loss of crispness, and
juiciness), and nutritional quality. The rate of water loss
is controlled by the vaporpressure difference between the
product and its environment, which is governed by the
temperature and relative humidity.
Ethylene (C2H4) is the simplest organic compound that
has an affect on physiological processes of plants. It is a
natural product of plant metabolism and is considered the
natural aging and ripening hormone. Generally, ethylene
production increases with maturity at harvest, physical
injuries, disease incidence, increased temperatures up to
30 C (86 F), and water stress. On the other hand, ethylene
production rates of fresh horticultural crops are reduced by
storage at the lowest safe temperature of each commodity.
Temperature affects both the rate of ethylene production and
the sensitivity of products to ethylene.
After harvest, fresh horticultural commodities are
susceptible to attack and damage by microbial pathogens.
Lower temperature affects the rate of growth and spread of
these microorganisms. Physical injuries can result from
abuses to fresh commodities at any temperature, but low
temperatures reduce the results of injuries.
Temperature is the most important environmental factor
that influences the deterioration rate of harvested
commodities. It is a well known fact that for each increase
of 10 C (18 F) above the optimum temperature, the rate of
deterioration increases two to threefold.
Thus improved cooling of fruits and vegetables before
or during shipment, along with proper temperature
maintenance throughout the marketing channels, has the
potential of greatly reducing these losses.
Nearly all fresh fruits and vegetables are now marketed
in corrugated fiberboard shipping containers. These
containers provide a barrier to proper air flow and
efficient heat transfer required for cooling. Although
cooling of fruits and vegetables before packing is much more
desirable from a heat transfer standpoint, most products are
not cooled until after packing because of other
relationships in the overall handling system. Cooling of
fruits and vegetables during transit is a particular problem
because of inherent restrictions on the size of the cooling
units and on the design of the air circulation system in
transit vehicles. Such systems are often inadequate for
cooling of fresh produce and problems of overheating, over
ripening, and losses due to decay are common. These problems
are further aggravated by increased use of palletized
shipping, introduction of new types of produce containers
and moves toward maximum utilization of vehicle volume
because of increased transportation costs. The need for
additional research on air movement as related to cooling of
fruits and vegetables in palletized shipments has been
widely recognized by industry sources.
During the last thirty years, numerous researchers and
commercial shipping companies have conducted research aimed
at improved cooling of fruits and vegetables packed in
shipping containers. Although substantial improvements have
been made, the current level of losses during marketing
indicates that much work must be done. Significant further
improvements can come only from complete design of the
containers, vent arrangements, stacking and loading
patterns, air circulation systems and refrigeration
controls, based on sound engineering principles. Before
this can be accomplished, however, basic research is needed
on the mechanisms of heat and mass transfer from the packed
product to the surrounding air and on the air flow
relationships involved. The objectives of this study follow.
OBJECTIVES
The objectives of this study were to
1. Study the basic principles related to air flow through
fruits (and vegetables) packed in fiberboard packing
containers.
2. Determine the feasibility of the theory of flow through
porous media analysis for this problem.
3. Develop a numerical model to predict air pressure and
velocity distribution for air flow through fruits (and
vegetables) packed in fiberboard packing containers.
4. Evaluate the experimental cooling response using the
mathematical model for flow in conjunction with an existing
heat and mass transfer model.
REVIEW OF LITERATURE
Background of Research
A considerable amount of information has been published
related to basic theory and experimental research on heat
and mass transfer of fruits and vegetables and on various
aspects of air cooling produce packed in shipping
containers. A number of publications (Anthony 1969, Camp
1979, Chalutz et al. 1974, Harris et al. 1969, Harris and
Harvey 1976, Hinds and Robertson 1965, Hinsch et al. 1978,
McDonald and Camp 1973, Mitchell et al. 1968, Pentzer et al.
1945, Redit and Hamer 1959, Rij 1977, Winston et al. 1959)
provide information on product cooling response during test
shipments of fruits and vegetables in refrigerated vehicles
as related to loading patterns and air flow arrangement.
Such tests give general information on whether a particular
system is better than another, but yield little of the
basic engineering information needed to design improved
transit cooling systems. It is obvious from most of these
tests that contact of the cooling air with the product
(rather than the capacity of the refrigeration system) was a
limited factor in the rate of cooling because of the rather
large cooling times in relation to the product heat loads.
Biales et al. (1973) and Kindya and Bongers (1979) conducted
shipping tests on cooling of grapefruit with outside ambient
air on nonrefrigerated ships. They showed that moderate
amounts of product cooling can be obtained in this manner in
certain regions during times of year when the outside air
temperature is relatively low. Such practice has
implications relative to cost savings and energy reduction
and points out the even greater importance of basic
engineering design of the system for optimum performance
when the source of refrigeration is limited.
Some researchers (Felsenstein and Zafrir 1975, Hale et
al. 1973, Hinsch et al. 1978, McDonald and Camp 1973,
McDonald et al. 1979, Olsen et al. 1960, Patchen 1969)
have developed new carton designs or have studied the
effect of different carton venting arrangements or air flow
directions through stacks of cartons as an aid to improved
cooling. For the most part, these were all comparative
tests and little basic engineering information on air flow
rates, air distribution and other important parameters was
obtained. Additional studies (Fountain and Hovey 1970,
Hale et al. 1969, Hale and Risse 1974, Olsen et al. 1960)
have involved packing fruits in shipping containers with
pulp or plastic tray packs. However, little information was
given to the direct influence of the tray packs on cooling
rates other than total cooling times or final temperature at
the end of test shipments.
A large amount of experimental research aimed directly
at studying cooling response of fruits and vegetables packed
in shipping containers as influenced by different variables,
has been conducted by various researchers: Chalutz et al.
(1974), Gentry and Nelson (1964), Grierson et al. (1970),
Grierson (1975), Guillou (1960), Haas and Felsenstein
(1975), Kasmire (1977), Kushman and Ballinger (1962),
Leggett and Sutton (1951), Lindsay et al. (1976), Mitchell
et al. (1972), Parson et al. (1972), Patchen and Schomer
(1971), Pentzer et al. (1945), Perry et al. (1971),
Richardson et al. (1955), Sainsbury (1951) and (1961),
Sainsbury and Schomer (1967), and Soule et al. (1969).
Results of such tests provide data on product temperature
response as related to product size, type of container,
venting and stacking arrangement and, in most cases, air
temperature and air flow rate. Most of these tests were
conducted in a valid manner and provide excellent general
information on the influence of different variables as they
affect cooling rates. Also, they provide some specific
engineering data useful in the design of cooling systems but
only for the particular product size, the particular
containers, stacking arrangements, vent designs, size of
product load, air flow rates and temperatures used in the
experiments. In most of these tests, the range of variables
studied was quite limited. In certain tests, the mass
transfer effects on product loss and evaporative cooling,
other than the overall product weight loss measurements
under the particular conditions of the test, were not
obtained.
In order to provide the necessary information for
design and analysis of alternate cooling systems,
experiments must be conducted in a much more basic manner
than those described above. The interaction of the cooling
air with the product, the individual containers, and stacks
of containers as related to the basics of heat and mass
transfer should be evaluated for different product sizes,
shapes, and packing configurations. More basic studies of
air flow patterns and pressure drops in and around cartons
and stacks of cartons are needed. A limited amount of this
basic type of research was conducted. Van Beek (1975)
studied the effective thermal conductivity of packed
spherical products and evaluated the effect of natural
convection for different heat flow directions. Meffert and
Van Vliet (1971) developed a test system for evaluation of
air circulation and product temperature in transit vehicles
and validated the results with experiments on model cargo.
Wang and Tunpun (1969) conducted basic studies on forced air
cooling of tomatoes in cartons. They presented information
on temperature response of fruits at 36 positions within two
cartons stacked in register at three different air flow
rates. They also conducted some basic studies of pressure
drop versus air flow relationships for bulk fruit and
cartons, packed and empty. Haas et al. (1976) studied
pressure drop relationships for oranges, and cartons,
packed and empty. They also studied the effect of vent hole
size, total area, position, and number on carton pressure
drop. Chau et al. (1983) studied pressure drop, as a
function of air flow and packing porosity, for four
different sizes of oranges with four different stacking
patterns. They also studied the effect the shape and number
of carton vent openings had on carton pressure drop.
Ramaker (1974) investigated the thermal resistance of
corrugated fiber board as used in shipping containers.
Some basic work has been done on analytical methods for
evaluating heat transfer in produce packed in cartons. Hicks
(1955) and Sainsbury (1961) used basic conduction equations
to evaluate theoretical cooling rates as a function of the
convective heat transfer coefficient at the container
surfaces. Meffert et al. (1971) and Meffert and Van Beek
(1975) used analytical solutions for transient heat transfer
to predict temperature response in stacks of produce
including heat generation from product respiration.
However, such solutions are applicable only to certain
limited cases. Information given in other publications
(Breakiron 1974, Gaffney 1977, Goddard 1972, Ryall and
Pentzer 1967, Stewart 1977) has pointed out the need for
additional basic research so that cooling systems may be
designed based upon sound engineering principles.
Researchers at the University of Florida have extensive
background in research related to heat and mass transfer in
fruits and vegetables. Most of the experimental work has
involved basic studies on cooling of oranges, grapefruits,
avocados and bell peppers in bulk loads (Baird and Gaffney
1974a, 1974b, and 1976; Chau et al. 1983, 1984; Gaffney
1977; Gaffney and Baird 1977, 1980). Certain methodologies
and techniques for measurement of product temperature, air
temperature, flow rates and humidity have been developed and
refined during the course of this research (Baird et al.
1975; Gaffney et al. 1980, 1985a, 1985b). In addition to
the experimental work, a computer based mathematical model
has been developed to analyze heat and mass transfer during
cooling of fruits and vegetables in bulk loads (Baird and
Gaffney 1974a, 1976; Eshleman et al. 1976; Chau et al.
1984). Currently experimental work in this area involves
measurement of thermal and physical properties and
transpiration coefficients of oranges and grapefruits and
theoretical studies of heat and mass transfer during surface
moisture drying of citrus.
This study was a subcomponent of a project concerning
heat and mass transfer relationships of fruit and vegetables
packed in fiberboard shipping containers. To understand and
model the heat transfer during cooling and storage of fresh
fruits and vegetables packed in shipping containers, the
pressure and velocity field characteristics within the
container must be established. Both distributions are
difficult to establish experimentally. There are many
interrelated variables involved in air cooling of fruits
and vegetables. These include thermal properties, physical
properties, size, and shape of the product and temperature,
flow rate and relative humidity of the cooling air. When
cooling products in containers, additional variables of
importance are container size, shape and wall thickness,
venting arrangements, stacking arrangements, product packing
configurations, and the direction of air flow.
The actual process that occurs between the individual
fruits (or vegetables) and the flowing air is complex and
not well understood. Previous researchers have studied the
cooling, heating and drying of semiinfinite systems of bulk
piled agricultural products such as fruits, grains,
vegetables, nuts, and root crops using the theory of fluid
flow through a porous media to determine the pressure and
velocity fields. The validity of porous media flow theory
for finite boundary conditions of air flow through fresh
fruits and vegetables packed in shipping containers will be
evaluated. The theory of porous media flow will be
considered before looking at previous work related to other
agricultural products.
Fluid Flow Through Porous Media
In addition to the interest received from a purely
scientific point of view, the theory of flow through porous
media has become important in many fields of scientific and
engineering applications.
Such diversified fields as ground water hydrology,
petroleum engineering, civil engineering, agricultural
engineering, chemical engineering, soil mechanics (physics),
water purification, industrial filtration, ceramic
engineering, powder metallurgy, medicine and bioengineering,
all rely heavily upon porous media flow theory as
fundamental to their individual applications.
Perhaps the most important of these areas of technology
that depend on the properties of porous media is hydrology,
which relates to water movement in earth and sand
structures, like dams, flow to walls from waterbearing
formation, intrusion of sea water in coastal areas, filter
beds for purification of drinking water and sewage, etc.,
and petroleum engineering which is primarily concerned with
petroleum and natural gas production, exploration, well
drilling, and logging, etc. In the introductory chapter of
the book by Muskat (1946), R.D. Wyckoff points out that
despite the great similarity of the physical systems and the
processes in these two fields of technology, each has
produced distinct technical literature and terminology.
Therefore, it is no surprise that as each branch of science
and engineering has addressed specific problems, each has
contributed to the vast amount of literature available
related to porous media flow theory. Unfortunately, this
literature is often difficult to correlate.
Porous media was defined by Collins (1961) as a solid
containing pores or voids of sufficient number, either
connected or nonconnected, dispersed within it in either a
regular or random manner. The interconnected pore space is
termed effective pore space. Fluid can flow through a
porous material only through the effective pore space.
Muskat (1946) represented an ideal porous media as a body of
unconsolidated sand in which all the voids are
interconnected, which means total pore space is actually
effective pore space.
By the above definition, a system of agricultural
products, in bulk piles or in storage structures, which are
of nonuniform shape and size, such as fruits, vegetables,
grains, beans, peanuts and root crops, can be considered as
an ideal porous media. A homogenous fluid, such as air or
water, used as a medium for cooling, heating or drying, can
flow through an assembly of interconnected channels of
varying sizes and shapes in this system of agricultural
products.
The actual process that occurs between individual
product particles and the flowing fluid is complex and not
well understood. For this reason, the only feasible method
of analysis is a macroscopic approach, which implies that
the values of variables under consideration are averages of
indeterminate instantaneous values.
Muskat (1946), Scheidegger (1960), and others indicated
that the earliest recorded investigation of the flow through
porous media was that of Darcy in 1856. Darcy was interested
in the flow characteristics of sand filters. Because of the
analytical difficulties in describing this type of flow,
Darcy had to resort to an experimental study of the problem.
His classic experiments led to the real foundation of the
quantitative theory of the flow of homogeneous fluids
through a porous media. He ran experiments in water flowing
through a vertical pipe filled with sand. From his
investigations, Darcy concluded that the flow rate was
proportional to the pressure loss and the crosssectional
area and inversely proportional to the length of the flow
path. The famous Darcy's Law is written as follows:
(1) Q = k A Ah/ L
Dividing both sides by the crosssectional area, A, yields
(2) V = k Ah/ L
s
where
Q = flow rate in volume per unit time
k = permeability or hydraulic conductivity
A = the crosssectional area
Ah = pressure loss or head
AL = the length of the flow path
V = the superficial velocity
s
The study of the validity of Darcy's Law has been the
subject of numerous investigations. These investigations
have been of essentially two types: those with the
objectives of either verifying Equation 1 or establishing an
appropriate modification of this equation and those
concerned with the nature of the constant, k, as determined
by the properties of the sand or porous media. The bases of
nearly all engineering calculations for porous media flow
problems have originated from Darcy's Law and/or purely
empirical determinations.
Later investigators found that Darcy's Law is limited
in its application only to very low velocity (creep) flow
and becomes invalid when inertial forces become effective.
Since then, several related theories and approaches have
been developed to approximate the flow of fluid through
porous media. Scheidegger (1960) recorded six different
theories, namely 1. heuristic correlation (curve fitting);
2. dimensional analysis; 3. capillaric model; 4. tube bundle
theory (sometimes known as Kozeny Theory); 5. drag theory;
and 6. statistical theory. Bird et al. (1960) recognized
two general theories, namely tube bundle and drag theories.
The most acceptable one and the one which has been pursued
by other workers is the tube bundle theory. This theory
assumes that the packed column is a bundle of tangled tubes
of irregular crosssection and that the established
principles of flow through a single tube are applicable.
Using the tube bundle approach in introducing the
parameter's porosity, hydraulic radius and mean particle
diameter into the HagenPoiseuille equation, Bird et al.
(1960) presented the following equation for laminar flow
through porous media:
2 2 3
(3) AP/AL = 150v V (1c) /D C
s p
where
= porosity
= void volume/total volume
p = fluid viscosity
D = mean diameter of particles of the porous material
P
AP = pressure loss
and the rest of the symbols are as previously defined.
This is the BlakeKozeny equation (first derived by
Kozeny and modified by Blake), which is valid for modified
Reynolds number equivalent to V D p/w(lE) less than 10.
This equation is identical to Darcy's equation where the
permeability, k, is expressed in terms of the characteristic
properties of the material and the fluid.
Similarly, for turbulent flow through a porous media,
the relationship is
2 3
(4) AP/AL = 3.50 (pV /2)(1e)/D
s p
where is the fluid density and the other symbols are as
previously defined. The kinetic energy term appears inside
the first parenthesis. This implies that the pressure drop
is proportional to the loss of kinetic energy which becomes
appreciable when the fluid velocity is large. This is the
BurkePlummer equation (first derived by Burke and Plummer),
which is valid for modified Reynolds number greater than
1,000.
Based on the theory of Reynolds for resistance to fluid
flow, Ergun (1952) illustrated that a pressure drop through
a porous media is caused by the simultaneous viscous and
kinetic energy losses. Pursuing this, he developed a
general equation of fluid flow through porous media, which
combines the BlakeKozeny equation for laminar flow and the
BurkePlummer equation for turbulent flow. This equation
(the Ergun equation) is
2 2 3 2 3
(5) AP/AL = 150P V (1E) /D E + 1.75 P V (lE)/D E
s p s p
The first term on the right side represents the viscous
energy loss and the second term on the right side represents
the kinetic energy loss.
In terms of dimensionless groups, the Ergun equation
can be written as follows:
3 2
(6) APD e / LpV (1e) = 150w(lE)/pV D + 1.75
p s sp
In this form, one can see that if the velocity becomes
very large, the first term on the right becomes negligible
and the Ergun equation reduces to the BurkePlummer
equation. On the other hand, if the velocity is very small,
the second term on the right is negligible compared to the
first and the Ergun equation reduces to the BlakeKozeny
equation.
In threedimensional form, the Ergun equation, Equation
(5), may be written as follows:
2
(7a) 9P/9x = f U + f U
1 2
and
2
(7b) 3P/ay = f V + f V
1 2
and
2
(7c) 9P/az = f W + f W
1 2
where
x,y,z = Rectangular coordinates
DP/3x = Pressure gradient decreasing along positive x
direction
DP/3y = Pressure gradient decreasing along positive y
direction
BP/8z = Pressure gradient decreasing along positive z
direction
2 2 3
f = 150l(lE) /D
1 p
3
f = 1.75p(lE)/D
2 p
U,V,W = Velocity in x,y, and z directions, respectively.
BakkerArkema et al. (1969), after an extensive review
of static pressure/air flow relationships in packed beds of
granular biological products, adopted the Ergun equation in
their experiments on onedimensional cooling of cherry pits.
They indicated that the Ergun equation was the most general
among the available equations.
Stanek and Szekely (1972, 1973) applied the Ergun
equation in twodimensional flow in both isothermal and non
isothermal packed beds. They found that nonuniform bed
porosity could cause flow maldistribution.
Lai and Haque (1976) applied the Ergun equation in
analyzing the threedimensional flow of air through non
uniform grain beds by adding a third coordinate axis, z, as
indicated in Equation 7. The equations were first written
in vector notation as follows:
(8) (grad P) = V (f + f IVI)
1 2
where
grad = gradient of a scalar function
V = Velocity in vector form
IVI = Magnitude of the velocity
and the rest of the symbols are as previously defined.
This, together with the continuity equation, was simplified
by transforming the equation into dimensionless forms and by
using a stream function. The resulting nonlinear partial
differential equations were solved numerically by a line by
line method. Nonuniform porosity was indicated to cause
different velocity and pressure distribution.
Lai (1980) used Ergun's equation to calculate the
threedimensional airflow and pressure distributions
through a circular grain bin using the method of lines,
coupled with a sophisticated ordinary differential equation
integrator. This method is similar to the finite difference
method.
Hague et al. (1980) modified Ergun's equation to
calculate the pressure and velocity fields for airflow
through conicaltop grain beds containing corn and fines
distributed nonuniformly. A twodimensional (cylindrical
coordinate) finite difference technique was used.
Geertsma (1974) reported it necessary to apply a more
general flow law for the prediction or analysis of the
production behavior of gas wells. These wells produce high
singlephase fluid flow rates through the surrounding
reservoir rocks. He indicated the appropriate formula was
given in 1901 by Forchheimer and is written as
(9) (grad P) = ajV + 8p JV V
where
a = coefficient of viscous flow resistance
= 1/k (approximately)
S= coefficient of inertial flow resistance.
This equation is noted to be of the same form as Equation 8.
Geertsma (1974) and Fircozabadi and Katz (1979)
discussed the widely differing views used to describe the
mechanisms that consume energy at more than a linear rate
with velocity. The 8 term has been given various names
depending on the users interpretation of the flow mechanism.
For example Fircozabadi and Katz (1979) indicated 8 should
be called "velocity coefficient" while Geertsma (1974)
defined B as indicated above as the "inertial coefficient."
Also noted was the difficulty relating the 8 term with the
porosity, permeability and other measurable porous media
characteristics.
Geertsma (1974) reported that the validity of Equation
9 could be checked by plotting the ratio pV/p as a function
of (grad P)/pV. Hie indicated that this plot must yield a
straignr line. Tne value of coefficient a will then be
determined by the intercept of the straight line plotted and
the (grad P)/pV axis while coefficient ( is found from the
slope of the line with respect to the xaxis. This
procedure is applicable for liquids only, since density p is
nearly independent of the pressure P.
For a gas, Equation 9 cannot be applied directly.
However, the equation is applicable for a small distance dx
in the flow direction. By introducing the mass flow rate
G = pV of the gas and using the ideal gas law, p = P/c in
which
c = RT/M
where
R = universal gas constant,
T = absolute temperature and
M = molecular weight,
Equation 9 can be rewritten as
2
(10) p dP = apGc+BG c/dx.
By performing integration between x = L and x = 0, Equation
10 becomes
2 2
(11) (P P )/2LcGy = a + B (G/p).
L 0
Therefore for gases, Geertsma (1974) reported that the
validity of Equation 11 could be checked by plotting G/p as
a function of the term on the left side of Equation 11. He
indicated that this plot must yield a straight line. The
value of coefficient a will then be determined by the
intercept of the straight line plotted and the axis of the
term on the left side of Equation 11, while coefficient P is
found from the slope of the line with respect to the xaxis.
The Agricultural Engineering Department, in addition to
other departments of the University of Florida, leases a
university version of the ANSYS*** general purpose finite
element computer program offered by Swanson Analysis Systems
(SASI) (DeSalvo and Swanson 1983).
***The information given herein is supplied with the
understanding that no discrimination is intended and no
endorsement by the Institute of Food and Agricultural
Sciences of the University of Florida is implied. The
listing of specific trade names herein does not constitute
endorsement of these products in preference to other
products which have equivalent capabilities.
During system verification and preparation for student
laboratory use, many of the practical structural, thermal
and flow engineering analysis techniques were tested. An
option exists that allows both 2D and 3D isoparametric
thermal solid elements to be used to model nonlinear
steadystate fluid flow through a porous medium (DeSalvo and
Swanson 1983).
The porous media flow problem is formulated in a
manner identical to that used for the thermal analysis,
requiring only a change of variables to use thermal analysis
to obtain a solution. Pressure is the variable rather than
temperature. The momentum equation is simplified to
(12) (grad P) = Reff V
where
grad = gradient of a scalar function
P = pressure
V = seepage velocity vector
and
(13) Reff = i/K + Spl1v
where
p = gas viscosity
K = absolute permeability of porous medias
B = viscoinertial parameter
p = density.
The components of Equations 13 are
(14) aP/3x = Reff Vx
(15) 2P/2y = Reff Vy
(16) aP/az = Reff Vz.
The continuity equation is
(17) div (pV) = 0
where div = divergence of a vector field, or
(18) 3(pVx)/3x + 3(pVy)/3y + a(pVz)/3z = 0.
Combining Equations 14, 15, 16, and 18 yields
(19) a(k3P/ax)/ax + 3(k3P/ay)/3y + a(kRP/az)/az = 0
where k =p /Reff. Equation 19 is nonlinear because Reff is
a function of velocity. The coefficients of permeability,
k, are (kx, ky, kz) internally calculated for each
coordinate direction as
(20) k = Kp/(p+Kpp6VI)
To handle the flow of power law fluids, an exponent can be
applied to IVI as an additional analysis input variable.
Combining Equation 12 and 13 yields
(grad P) = IV/K + BpIVj V,
~
which is the same as Equation 9. In fact, personal inquiry
to SASI revealed that the porous media option was based on
the work of Geertsma (1974) and Fircozabadi and Katz (1979).
Another modification of Darcy's Law that has received
considerable attention is a simplified, semiempirical
approach to account for nonlinear flow patterns of air
through porous media. This is accomplished by assuming that
the air velocity is proportional to the pressure gradient
raised to a power. Shedd (1953) introduced the following
equation in the study of air flow through grain storage:
B
(22) V = A(DP/3n)
where
V = interstitial fluid velocity
DP/an = pressure gradient along any direction
A,B = experimentally determined constants.
Equation 22 is sometimes rearranged in terms of
b
(23) 3P/an = aV ,
where
a,b = experimentally determined constants
and referred to as the Ramsin Equation.
One major difference between Equation 22 and the Ergun
Equation, Equation 6, is the absence of porosity. It is
assumed to be lumped with other characteristic properties of
the fluid and porous medium into constants a and b.
Shedd (1953) stated that A and B were not constant over
a wide variation in pressure gradients. Equation 22 is a
straight line on a loglog plot; however most of the
experimental data plots did not produce straight lines.
By generalizing information given by Shedd (1953),
Brooker (1961) obtained the twodimensional velocity
component equations
2 2 (B1)/2
(24) Vx = A [(SP/9x) + (OP/3y) ] P/9x
and
2 2 (Bl)/2
(25) Vy = A [(DP/3x) + (2P/Dy) ] 3P/3y
where Vx and Vy are the velocity components and the rest of
the symbols are as previously defined.
Brooker (1961) substituted Equations 24 and 25 into the
continuity equation, Equation 18 assuming constant density
and considering only the x and y directions, to obtain the
partial differential equation
2 2 2 2 2 2
(26) [(SP/3x) + (8P/8y) ] [ P/Dx + 3 P/y ]
2 2 2 2
2m[(8P/3x) 3 P/3x + 2(DP/x) (3P/3y)(3 P/3x8y)
2 2 2
+ (3P/3y) 3 P/iy ] = 0
where
m = (1B)/2.
The finite difference method was used to solve this
differential equation. The numerical solution agreed
favorably with experimental pressure values measured in a
rectangular wheat bin.
Brooker (1969) proposed that the loglog plots of
velocity versus pressure gradient be modeled by a series of
straight line segments using Equation 22 and reported
several values of A and B for a range of pressure gradients.
The finite difference method was used to solve Equation 26
with the value of B adjusted at each node corresponding to
the value of the gradient P/3n. Brooker (1969) used this
model to calculate pressure and velocity distributions for
twodimensional nonlinear airflow through a rectangular corn
bin. The numerical results did not agree well with
experimental findings.
Jindal and Thompson (1972) applied Brooker's (1969)
model to calculate the twodimensional pressure distribution
for air flow through conicalshaped piles of grain sorghum.
Pierce and Thompson (1974) also applied Brooker's
(1969) model by modifying the Jindal and Thompson (1972)
model for the prediction of airflow characteristics in
conical piles of corn or sorghum.
Marchant (1976a) applied Shedd's equation and used the
finite difference method to establish the total airflow
through large round hay bales. Threedimensional airflow
was discussed briefly. Marchant (1976b) also used the
finite element method to solve linear airflow problems in
two dimensions. The experiments of Brooker (1958) and
Borrowman and Boyce (1966) were used in verifying the model.
Segerlind (1982) developed a model for twodimensional
nonlinear airflow and used the finite element technique to
solve for the pressure distribution in rectangular grain
bins. The model was based on a generalization of Shedd's
equation and calculated constant pressure lines through
shelled corn similar to the lines measured by Brooker (1969).
Smith (1982) applied the two dimensional finite element
method used by Marchant (1976b) to estimate the pressure and
velocity distributions in a threedimensional field of a
rectangular heap of grain with an onfloor monoduct system.
A comparison was made between experimental values of
pressure and velocity of Marchant (1976a) and the calculated
results using the finite element method. Smith (1982) found
that the pressure could be calculated more accurately than
velocity but that both could be calculated with reasonable
accuracy in practical problems.
Khompis (1983) applied Shedd's equation and used the
finite element method to solve the threedimensional
nonlinear pressure and velocity distribution within
cylindrical grain storage. Four types of perforated floor
systems were studied using corn as the grain media.
Extensive graphics were demonstrated to display the three
dimensional pressure and velocity distributions.
From the above discussion it is evident that the theory
of porous media flow has been used extensively in the study
of air flow through agricultural commodities. This previous
work can be generally classified into three areas of study.
The first area of study is a development of equations
for prediction of pressure drop as a result of air flow
through bulk fruits and vegetables. In addition to the
previous work cited, several workers applied the Shedd or
Ergun equation to determine the relationship of pressure
drop to velocity of air flow through the various
commodities. Patterson et al. (1971) studied cherry pits,
shell corn, and navy beans using both the Shedd and Ergun
equations. Steele (1974) studied peanuts using Shedd's
equation. Akritidis and Siatras (1979) used Shedd's
equation in their study of pumpkin seeds. Neale and Messer
studied root and bulb vegetables (1976) and leafy vegetables
(1978), also using Shedd's equation. Gaffney and Baird
(1977) studied bell peppers using Shedd's equation. Wilhelm
et al. studied snap beans (1978) and lima beans and southern
peas (1981) using Shedd's equation. Abrams and Fish (1978)
applied Shedd's equation in their study of sweet potatoes.
Rumsey (1981) used the Ergun equation in the study of
English walnuts.
The second area of study is the prediction of pressure
and velocity field distributions for airflow through bulk
agricultural products. These are sufficiently outlined
above. It is noted that each of these individual studies
has several common aspects. They all use some modification
of Darcy's approach to porous media flow. They all involve
a large bulk of porous media (semiinfinite). The air
entrance into the porous media was through a line or area
source. The air exit from the porous media was through a
large exit area. Finally, these various studies used
several modeling analysis techniques which were specifically
developed for the particular situation being studied. This
factor is emphasized as a result of difficulty in using the
various models for other than the specific situation for
which they were intended.
The third area of study and most limited is the
development of equations to predict the pressure drop for
airflow through fruits and vegetables packed in shipping
containers.
As mentioned above, Wang and Tunpun (1969) studied the
pressure drop versus airflow relationships for tomatoes in
bulk and incartons. Haas et al. (1976) used a modified
Shedd equation to predict equations relating pressure drop
to air velocity for oranges inbulk and incartons. Chau
et al. (1983) used both the Shedd and Ergun equation to
predict the pressure drop as a function of airflow for
oranges inbulk and incartons.
No studies exist which predict the pressure and
velocity fields for air flow through fruits or vegetables
packed in shipping containers.
PROCEDURE
In order to better understand the heat transfer and
the product temperature change which occurs when forced air
cooling is used to reduce the temperature of fruits and
vegetables packed in fiberboard containers, the velocity and
pressure fields within the carton must be analyzed.
Intuitively it was anticipated that the velocity
characteristic of the air flowing through packed cartons of
fruits and vegetables would vary in three dimensions. As
mentioned previously, direct measurement of the pressure and
velocity within and throughout the packed carton is not
feasible with current instrumentation. The measurement of
temperature is the only direct method which can be used as
an indication of the velocity pattern of the air passing
through the carton as the product is cooled.
Extensive temperature response experiments were
conducted with oranges packed in an experimental packing
container of a size similar to regular packing cartons. The
number and location of vent holes could be changed and
several flow rates were used. Eighty 36gage thermocouples
were used to measure air temperature entering and leaving
the box, air temperature at various locations inside the
box, and surface and interior temperatures of oranges at
various locations within the box.
The temperature response data were analyzed graphically
and by computer simulation which indicated changes in
temperature within the carton through various cutting planes
by use of color variation for various temperature ranges.
From this temperature response data, no symmetrical
patterns could be identified. This leads to the assumption
that the airflow within the orange carton varies in three
dimensions.
In order to evaluate the velocity and pressure fields,
an analysis technique had to be selected. The two principal
numerical solution methods employed by previous researchers
are the finite difference and the finite element techniques.
Each method has advantages and disadvantages.
The flow of air through a porous media may be described
in terms of partial differential equations. In real
situations, the geometry of the system to be modeled is
usually not simple and these equations cannot be solved by
analytical solution. In these cases the use of some
numerical solution is necessary. Brooker (1961) and others
proposed a finite difference method of solution. In this
method the porous medium is divided into a regular grid.
The partial differential equations of static pressure at
each grid point are approximated by truncated Taylor Series
expansion in terms of the static pressure at the surrounding
points. In this way algebraic equations can be formed for
each grid point therefore setting up a system of equations
to be solved simultaneously. While the finite difference
method can give good results, it has several disadvantages:
1. The method is not easily adapted to regions which have
boundaries not lying on the grid lines.
2. It is difficult to write a computer program to take
account of a general geometric shape.
3. A grid size which gives acceptable accuracy in areas
with rapidly varying pressure will be far too small in the
areas of slight variation, thus leading to a solution of a
large number of superfluous equations.
4. It is difficult to apply the method to media which are
not homogeneous and isotropic.
5. For each equation relating air velocity to pressure
gradient, a completely different partial differential
equation has to be solved thus requiring a new set of
equations to be set up in each case.
The finite element method is especially suited for the
study of airflow through fruit and vegetable crops because,
as well as overcoming the general points 1 to 3, it can also
deal easily with 4 and 5 which are very pertinent.
In the past ten years, much work has been done to
develop the numerical solution of partial differential
equations using the finite element method. The finite
element approach has its origins in the field of solid
mechanics (elasticity, plasticity, and structural analysis),
but its application for use with a much wider range of
problems was soon realized. Finite element methods have
replaced finite difference methods in many areas of solid
mechanics and are making inroads into fluid mechanics, heat
transfer, and other fields. More recently, considerable
interest has been shown in its application to nonlinear
problems, such as the problem under consideration.
The finite element method is a technique that
approximates a continuous quantity, such as a partial
differential equation, by a discrete model composed of a set
of piecewise continuous functions. This discretization is
done on a small but finite portion of the component. As its
name suggests, the finite element method consists of
dividing the region into a number of smaller elements which
can be of any shape (Figure 1). Nodal points are
distributed around the boundary of the element and
occasionally inside it. These nodes are arbitrarily
numbered. The overall model is comprised of a finite number
of nodes. At each node the value of the variable (i.e.,
pressure, flow, temperature, deflections) for a differential
equation is either known or unknown. Between the nodes are
subdomains called elements. Mathematics converts the
partial differential equations into matrix equations that
describe each element. The individual elements (in matrix
form) are joined together to form an overall model. The
results of the finite element simulation output are given in
terms of the variable at each nodal point. After the known
variables are obtained, internal quantities can be
Figure 1. Region divided into finite element.
estimated. An example structural problem would solve for
the deflection on a loaded component and then evaluate the
stresses inside each element.
With the finite element approach, the partial
differential equations describing the desired quantity (such
as pressure, flow, temperature, displacement, deflections)
in the continuum often are not dealt with directly.
Instead, the continuum is divided into a number of finite
elements, which are assumed to be joined at a discrete
number of points along their boundaries. A functional form
is then chosen to represent the variation of desired
quantity over each element in terms of the values of this
quantity at the discrete boundary points of the element. By
using the physical properties of the continuum and the
appropriate physical laws (usually involving some sort of
minimization principal), a set of simultaneous equations in
the unknown quantities at the element boundary points can be
obtained. This set of equations is in general quite large,
but the matrix is banded.
For those situations where the finite element
technology has been developed, there are three primary
advantages of the finite element approach over the finite
difference methods. These are
1. Irregularly shaped regions can be handled easily,
without the special treatment usually required by finite
difference methods.
2. The size of the finite element can easily be varied over
the region, permitting the use of small elements where
strong variations occur and large elements where only
slight variations are expected. With finite difference
methods, at least in their conventional form, the use of
several different mesh sizes can cause bookkeeping
difficulty.
3. For comparable accuracy, the finite elements can usually
be considerably larger than the mesh elements of a finite
difference grid. As a result, when elliptic problems are
involved the band matrix referred to earlier is usually
small enough to be solved directly without recourse to the
iterative methods which are usually necessary in finite
difference methods.
The main advantage of the finite element method is that
a general computer program can be written to create the
overall model matrix. For different applications (e.g.,
different shape geometry) input is easily changed and a new
solution obtained with the same computer program. Other
important advantages include the following:
1. Material properties can be different between each
element.
2. Irregular shaped geometries can be modeled.
3. The size of elements can be varied.
4. Mixed boundary conditions can be handled (e.g., loads
and temperatures).
The number of nodes in an element can vary with the
addition of internal nodes that provide a means to curve the
element or specify a nonlinear variable. Many types of
elements are available (Figure 2).
POINT (MASS)
LINE
(SPRING, BEAM, SPAR, GAP)

AREA
(THIN SHELL, 2D SOLID,
AXISYMMETRIC SOLID)
VOLUME
(3D SOLID)
THICK SHELL
Figure 2. Various types of finite elements.
N2
E: 40
The two methods (finite element and finite difference)
are similar in that an algebraic equation is set up for each
nodal point, but the way in which they are set up differs.
Compared with finite difference programs, for example, they
are easier to use and require fewer computer resources. For
modeling nonlinear material properties, finite element
programs require less than 20 iterations compared with the
200 or 300 required for finite difference programs. With
recently developed finite element software, the calculated
distributions can be shown graphically on a computer screen,
which makes it relatively easy to see areas of high and low
pressure distribution. Working with a geometrical model of
a particular size, an engineer can vary the input
specifications such as material properties and inlet and
outlet locations and see the effect of these variations on
the resulting pressure distributions. The program also
calculates performance parameters. With this information,
the engineer can upgrade the design geometry and change the
material requirements and modify the energy losses, and
thereby achieve the best design performance with the least
amount of material or achieve the optimum cooling. The
procedure is quick and relatively inexpensive when compared
with conditional procedures such as building and testing
prototypes.
After review of the above literature and analysis of
the temperature response data, application of a 3D finite
element nonlinear porous media flow analysis for the problem
under consideration was thought feasible.
As mentioned above, previous modeling involved a porous
media which was of much larger scale and without the more
significant boundary conditions for the situation of air
flow through a packing container. Also the finite element
models developed by previous workers were developed
specifically for the problem which they were addressing.
Although the 3D finite element porous media flow analysis
appeared feasible, in the early stages of this study it was
not deemed appropriate to devote a considerable amount of
time in developing a specific model for this situation
before feasibility of the approach was confirmed.
Finite element programs are available and are becoming
the most utilized engineering tool for design. These large
and powerful tools have been written by engineering graduate
students and professors with advanced technical knowledge in
their areas of expertise and therefore might pose a problem
to the majority of engineers who have little background in
finite elements and other dynamic simulation computer
programs. This leads the average engineer relying on the
program written by others and thus using the analysis
techniques as a black box: entering input and receiving
output. The degree of success in using a finite element
program is related to the assumptions and models developed
by the program originator. The user of the finite element
program must realize that the actual results from the finite
element model are a function of engineering judgement used
to select the element type and size, and a function of the
computer program with its associated mathematical
assumptions.
Also mentioned above was the discovery that the SASI
finite element analysis package, DeSalvo and Swanson (1983),
contained an option which allows the modeling of nonlinear
steady state fluid flow through a porous medium. This
solution technique was studied and appeared to be applicable
to the case at hand. A limited review was unable to
identify other commercial solution packages which offered
the needed analysis option and capabilities.
If the SASI package or a similar general finite element
computer program was readily accessible to the workers
interested in porous media flow analysis applied to
agricultural crops, a tool would be available which would
allow for a uniform solution technique and easy comparison
of the results of various workers. Therefore, in addition
to evaluating the feasibility of the finite element porous
media flow analysis, this study is also evaluating the
feasibility of applying a commercial general computer
program. If it is established that the analysis technique
is feasible, as well as, the general computer program, then
the industry which provides such software can be encouraged
and solicited to provide the needed software for widespread
adoption and application.
Commercial Finite Element Package
The ANSYS finite element computer program is based on
classical engineering concepts and documented finite element
and numerical analysis techniques (DeSalvo and Swanson
1983). This selfcontained general purpose program uses
efficient solution techniques to solve a large class of
engineering problems. It is user oriented, command driven,
has extensive graphic capabilities, and is documented,
benchmarked, and verified.
A typical analysis consists of three phases: Pre
Processing (Analysis Definition), Solution, and Post
Processing (Interpretation of Results).
The PreProcessing phase is very important since the
accuracy of the solution depends directly upon the degree of
accuracy of the problem description. Input data prepared in
the analysis definition would include the model description,
boundary conditions, and the analysis type and options.
The model description involves creating the desired
geometry, selections) from the element library,
specification of geometric (real) constants describing
properties of elements (e.g., area, moment of inertia, and
height of a beam), and identification of material properties
(e.g., viscosity, conductivity, and density). The user must
ensure dimensional homogenity.
The analysis is performed in the solution phase. For
the nonlinear porous media flow case, this involves the
solution of the matrix equation
(27) [K] {P} = {Q}
where
[K] = the transmissivity matrix, defined by Reff in
equation 13.
{P} = pressure vector (unknown)
{Q} = mass flow rate vector
and the calculations of the pressure and mass flow
distributions.
The solution results are evaluated in the third phase.
The user should determine if the objective of the analysis
was met. Tools to use in making this decision include
graphics (contour maps, distorted shape plots, and graphs of
one variable versus another); scans of results; and
combinations of results.
Therefore, the use of ANSYS is straightforward, but
careful engineering judgment must be exercised during the
analysis definition phase in selecting the geometry and
elements, as well as specifying the input parameters and
boundary conditions.
Verification
In order to verify the ANSYS calculated solutions, the
solutions must be compared with known theoretical solutions,
experimental results, and/or with other calculated
solutions. The extensive verification problem manual,
DeSalvo (1982), provided by SASI does contain a ground water
flow verification problem. However, it does not contain a
porous media flow problem of the type under consideration
with which to compare the ANSYS output to a theoretical
solution. Therefore, before attempting to solve the problem
under consideration, ANSYS was used to solve similar porous
media flow problems which have been solved by other workers.
ANSYS was used to solve the twodimensional rectangular
grain bin problem solved by Segerlind (1982) and one of the
threedimensional grain bin problems solved by Khompis et
al. (1984). As noted above, ANSYS is based on an Ergun
Equation analysis while the work of Segerlind (1982) and
Khompis et al. (1984) is based on a Shedd Equation analysis.
Patterson (1969) determined the coefficients for pressure
drops through grain beds using both the Shedd Equation and
Ergun Equation. It is assumed that appropriate coefficients
can be selected for the Ergun Equation, that correspond to
the coefficients used in the Shedd Equation by Segerlind
(1982) and Khompis et al. (1984). It is further assumed,
that ANSYS provides a satisfactory solution if the
calculated results are within reasonable accuracy of the
results determined by Segerlind (1982) and Khompis et al.
(1984).
The input coefficients needed to use ANSYS are
indicated in Equation 13. The gas viscosity and the density
are for air and are taken from standard data tables. The
absolute permeability of the porous media and the visco
inertial parameter are derived from experimental results.
The work of Patterson (1969) is used to obtain values for
these last two parameters.
Patterson (1969) used a modified Ergun Equation
2 2 3 2 3
(28) AP/AL = Ke [150 PV (1E) /D + 1.75 p V (lE)/D ]
s p s p
where
Ke = modified Ergun product constant.
Equation 28 is identical to Equation 5 with the
addition of the product multiplication factor, Ke.
Comparing Equation 28 and Equations 12 and 13, the absolute
permeability of the porous media, K, and the viscoinertial
parameter, 0, needed for ANSYS input can be equated as
follows:
2 2 3
(29) 1/K = Ke [150 (1E) /D E ]
p
and
3
(30) A = Ke [1.75 (lE)/D E ].
p
Patterson (1969) reported all the parameters required
to solve Equations 29 and 30, for shelled corn and various
other biological materials. For shelled corn at a moisture
content of 16.01 %, a temperature of 87 F, and a normal
fill, Patterson (1969) reported the following values:
D = 0.032 feet; E = 0.43; and Ke = 4.5
The conditions of the corn are similar to that used in
the work of Segerlind (1982) and Khompis et al. (1984).
Both studies used coefficients which were reported in the
research of Shedd (1953). This research was based on normal
fill, shelled yellow dent corn at a moisture content of
12.4 %. However, Shedd (1953) reported only a 1 % variation
in the pressure drop for the same corn at a moisture content
of 15.8 %.
Solving Equations 29 and 30 using the values given
above results in the following:
6 2 7 2
1/K = 2.69 x 10 ft ; or, K = 3.71 x 10 ft
3 1
and 8 = 1.76 x 10 ft The following viscosity and
density for the air were used:
5 2 3
S= 1.24 x 10 Ibm/ftsec ; and p = 7.25 x 10 lbm/ft
The English Engineering system of units was used by
Segerlind (1982) and the SI system of units was used by
Khompis et al. (1984) and Chau et al. (1983). Conversion
factors were used to modify the units of the viscosity and
density to insure dimensional homogenity. When the desired
units for the pressure are inches of water and the desired
units for the velocity are feet per minute, the unit
conversion results in the following values for viscosity and
density for input to ANSYS:
48
9
p/g = 1.24 x 10 in. watermin;
c
7 2 2
p/g = 1.20 x 10 in. watermin /ft ;
c
where the Newton's Law gravitational constant is
2
g = 32.2 ftlbm/lbfsec
c
When the desired units for the pressure are Pascals and
the desired units for the velocity are meters per second,
the unit conversion results in the following values:
5 3
p = 1.846 x 10 kg/msec; and p = 1.1614 kg/m.
After conversion of the units, the other input
parameters for ANSYS become
7 2
1/K = 2.9 x 10 m ; or,
3 1
and B = 5.77 x 10 m
8 2
K = 3.45 x 10 m
Now all the necessary input data to use the ANSYS
finite element program are known.
TwoDimensional Grain Bin
The physical problem selected for solution by Segerlind
(1982) was the same rectangular grain bin that was analyzed
by Brooker (1961) and (1969). The bin was 121.92 cm wide (48
inches) and 243.8 cm high (96 inches). The duct was 20.32 cm
high and 40.64 cm wide (8 inches by 16 inches).
Segerlind (1982) utilized the six node quadratic
triangular element in his finite element method solution.
The bin was divided into 32 elements with 85 nodes. The
horizontal node spacings was 10.16 cm (4 inches) while the
vertical spacing was 5.08 cm (2 inches) in the region of the
inlet duct. The upper 76.2 cm (30 inches) was modeled with
two elements because all of Brooker's (1969) results
indicated that a linear pressure gradient existed in this
region. The element and node locations used by Segerlind
(1982) are shown in Figure 3 and the boundary conditions
are shown in Figure 4. Segerlind (1982) noted that the
known pressure values at the duct inlet and at the free
surface are easily incorporated into all finite element
programs. The noflow or impermeable boundary condition,
3P/3n = 0, is automatically enforced on the boundary when
the pressure values are not specified. This differs
significantly from the finite difference technique used by
Brooker (1969) where extra rows of nodes had to be added
outside the region to satisfy the no flow boundary
condition.
 24" 
96"
8"
Figure 3. The node and element grid employed by Segerlind
(1982).
P= 0
p= 0
aP
aP o / ax
ax
aP
ayESSURE SPECIFIED
PRESSURE SPECIFIED
Figure 4. Boundary conditions used by Segerlind (1982).
The ANSYS analysis option to be used did not contain
an element type identical to the type used by Segerlind
(1982). A four node 2D isoparametric thermal solid element
was used. The element option which allowed modeling of
nonlinear steadystate fluid flow through a porous media
was used. This element could be formed into a triangular
(three node) shape by the duplication of two of the four
nodes. The node and element locations used are shown in
Figure 5. The model consisted of 85 nodes and 71 elements.
All but nine nodes were placed identical to the locations
used in Segerlind's (1982) model. Over twice as many
elements were required because smaller elements were used
which did not contain the midpoint nodes like the elements
used by Segerlind (1982). Although this resulted in
additional computer solution time, the results were as
refined or more refined when compared to Segerlind's (1982)
solution. The boundary conditions were specified in a
manner similar to Segerlind's (1982) model.
ThreeDimensional Grain Bin
Khompis et al. (1984) extended the procedures of
Segerlind (1982) to account for three dimensional pressure
and velocity fields in cylindrical grain storage. Four
different aeration systems were studied. The four systems
were differentiated by the size and and shape of the
perforated floor. The system to be considered in the
present study consisted of the square perforated type floor
Figure 5. The node and element locations used in current
study.
for a bin 12.8 meters high and 18.3 meters in diameter. The
perforated floor was 9.15 x 9.15 meters square. Specific
details of these sizes and shapes are shown in Figure 6.
Flush floor ducts were assumed.
Khompis et al. (1984) used a twenty node, quadratic,
three dimensional element for the finite element model. The
cylindrical grain storage investigated were symmetrical and
the finite element analysis was performed using only half of
the grain bin. The limitation of computer central memory
and desire for reasonable computing restricted the finite
element model to five layers with 12 elements per layer.
The 60 elements had 406 nodes. The height of the elements
at each layer depended on the distance from the perforated
floor and height of the grain storage. The shortest
elements were located in the layer in contact with the floor
and had a height of about 5 percent of the total grain
storage height. The element grid is shown in Figure 7.
The boundary conditions consisted of a pressure of 500
Pascals (2 inches of water) specified at the square
perforated floor and atmospheric pressure at the free
surface on the top of the grain bin. Again, the noflow or
impermeable boundary condition, 2P/8n = 0, is automatically
enforced on the boundary when the pressure values are not
specified.
Again, the identical type element which Khompis et al.
(1984) used was not available for use with the ANSYS solution
package. An eight node 3D isoparametric thermal solid
FULL PERFORATED FLOOR
SQUARE DUCT
Y DUCT
Figure 6. Four different grainbin aeration systems studied
by Khompis et al. (1984).
STRAIGHT DUCT
x
Crosssection of the elements
in the XY plane.
Y
Division into the elements
of a half cylindrical
grain storage.
Figure 7. Finite element model used by Khompis et al.(1984).
element, which allows modeling of nonlinear steadystate
flow through a porous media, was used. This element could
be formed into prism or tetrahedron shaped elements by
defining duplicate node numbers. Although ANSYS did have
twenty node elements like those used by Khompis et al.
(1984), these elements did not have the porous media flow
option capability. The node and element locations used with
ANSYS are shown in Figure 8. The model consisted of 741
nodes and 552 elements.
The comparison of results with those of Segerlind
(1982) and Khompis et al. (1984) received detailed coverage
in the Results and Discussion section, below. However, the
results supported the use of ANSYS to study the problem
under consideration. The procedures for this application
follow.
Solution for Oranges Packed in Containers
The use of the ANSYS porous media analysis to study the
pressure and velocity distributions as air flows through a
container of oranges presented several questions that were
not significant problems for previous workers. A primary
concern was the overall scale of the porous media used in
this case. This scale was finite when compared to the semi
infinite cases studied by others. Because of the small
dimension of the packing container, boundary (wall) effects
could be significant. The wall contact with the fruit or
vegetable presented two possible difficulties. The drag
E
N
30 #.m am( SIATIC P~U~ftAWInfl .ATrmUW
Figure 8. The node and element locations used in current
study.
caused as the air passes the wall was one consideration.
The second concern, which has been reported by other
workers (Pillai 1977, Ridgway and Tarbuck 1968, Stanke and
Eckert 1979) relates to the variance of the voidage or
porosity adjacent to the walls when compared to the central
portion of the porous media.
Other researchers have reported the significant
pressure drops produced by the air inlets of packing
containers. This point was another possible significant
difficulty when comparing the pressure drop of the porous
media (fruit or vegetable within the cartons) to that of the
pressure drop across the inlet(s) and exit(s).
Another consideration in terms of fruits and vegetables
was the compaction possible during packing and subsequent
shrinkage with time due to physiological changes and
moisture loss. The porosity may change, for example, if the
product become more compact due to shrinkage or handling.
The porosity next to the inside top of the carton could
increase allowing more air to pass through this area.
The impact of some or all of these points of interest
is addressed in the Results and Discussion section.
Before verifying the applicability of the ANSYS finite
element program modeling air flow through a carton of fruit,
the work of Chau et al. (1983) was modeled in order to
better understand how the model of a commercial orange
carton should be developed.
ThreeDimensional Bulk Citrus
As previously mentioned, Chau et al. (1983) used both
the Shedd and Ergun equations to predict pressure drop as a
function of air flow for oranges in bulk and packed in
cartons. This work involved the determination of pressure
drop though empty cartons (for a single carton wall and for
both walls of an empty carton) and for oranges packed in the
simulated carton. Fruit stacking pattern (bed porosity),
fruit size, vent hole shape, vent hole number, and air flow
rate were some of the variables studied.
The experimental equipment employed by Chau et al.
(1983) is shown in Figure 9. The crosssection of the
product bin was 55 cm x 55 cm. The actual crosssection
for each test varied slightly depending on fruit size and
stacking pattern.
There were 16 tests involving oranges in bulk, which
consisted of four orange sizes and four stacking patterns.
Since Chau et al. (1983) presented more data for Florida
size 80 (8.1 cm) waxed Valencia oranges stacked in a square
staggered pattern than other combinations of size and
stacking pattern, the present study considered only these
data. In the square staggered pattern each layer had a
square pattern (every 4 adjacent fruit forming a square)
with the next layer offset so each fruit on the subsequent
layer rested on 4 fruit from the previous layer. The air
flow rate was varied for each test to cover a range of
approach velocities from 0.1 to 1.0 meters per second.
Chau et al. (1983) fitted the experimental data using
the Ramsin Equation (Equation 22) in terms of pressure
change per bed depth of fruit:
b
(31) AP/h = a V
For the size 80 oranges stacked in a square staggered
pattern, with a bulk depth of 0.33 meters and a porosity of
0.435, the constants in Equation 31 were found to be
a = 278; and b = 1.937.
Chau et al. (1983) presented a variation of Ergun's
equations as
2 2 3 2 3
(32) AP/h = K P V (1E) /D E g + K P V (lE)/D E g
1 s p c 2 s p c
where
h = AL in Equation 5
K = 150 in Equation 5
1
K = 1.75 in Equation 5
2
g = Newton's Law gravitational constant.
c
This work did not employ a modified Ergun product
coefficient, Ke, like that used by Patterson (1969).
Therefore, eliminating Ke and substituting K and K
1 2
Equations 29 and 30 can be rewritten
2 2 3
(33) 1/K = K [ (1E) /D E ]
1 p
and
3
(34) = K [ (1e)/D e ].
2 p
Chau et al. (1983) determined K and K by fitting
1 2
the experimental data and reported all the parameters
required to solve Equations 33 and 34 for bulk packed
oranges. For size 80 oranges arranged in a square staggered
stacking pattern, the following values were reported:
D = 0.081 meters; E = 0.435; K = 890; and K = 2.73.
P 1 2
Solving Equations 33 and 34 using the values given
above resulted in the following:
5 2 6 2
1/K = 5.2608 x 10 m ; or, K = 1.9009 x 10 m
1
and 8 = 231.3 m
Since the SI system of units was used by Chau et al.
(1983), the desired units for the pressure are Pascals and
the desired units for the velocity are meters per second.
The inlet air condition used by Chau et al. (1983) are
approximately the same as those used by Segerlind (1982) and
Khompis et al. (1984). The air viscosity and density
specified above as input for the ANSYS finite element
program to solve the problem reported by Khompis et al.
(1984) were used again,
5 3
V = 1.846 x 10 kg/msec; and p = 1.1614 kg/m .
To determine the pressure drop across rectangular and
round vent holes, Chau et al. (1983) conducted tests using
the experimental setup shown in Figure 9, for one layer
of corrugated fiberboard (simulating one side of a carton)
with 4, 6, and 8 vent holes. The tests were repeated with 2
layers of fiberboard separated by 30 centimeters (simulating
2 sides of a carton). Each fiberboard side had 4, 6, or 8
vent holes placed uniformly in the fiberboard crosssection.
The crosssectional area of each vent hole was 12.94 square
centimeters. The thickness of the fiberboard used to make
the cartons was 3.2 millimeters.
The experimental data were normalized on the basis of
actual air velocities through the vent holes (calculated as
total flow rate divided by total vent area), and the data
points for the tests with 4, 6, and 8 vent holes fitted the
same straight line on a loglog plot.
The equations of best fit for the experimental data
were presented by Chau et al. (1983) as pressure drop
across one side of carton with rectangular vent holes
2
(35) AP = 1.59 V ;
PRESSURE
TAPS
R
BLOWER
FLOW ELEMENT7
I,
DUCT
EXTENSION
: FRUIT
AI R
L AIR
STRAIGHTENER
L SLIDING
DAMPER
Figure 9. The experimental setup used by Chau et al. (1983)
for measuring air flow resistance in orange.
1.
I
AI
pressure drop across two side of carton with rectangular
vent holes
2
(36) Ap = 2.78 V ;
1
pressure drop across one side of carton with round vent
holes
2
(37) AP = 1.76 V ; and
i
pressure drop across two side of carton with round vent
holes
2
(38) AP = 2.89 V
i
Chau et al. (1983) noted that the pressure drop across
two sides of the carton was less than twice the pressure
drop across one side of the carton alone. The explanation
for this difference given by Chau et al. (1983) was that the
loss for one carton side is similar to an exit loss. For
two sides of the carton, the loss is more like a duct
enlargement loss as the air travels between the sides of the
carton. The enlargement loss is smaller than an abrupt exit
loss.
The explanation for the slightly higher pressure drop
for the round vent holes was attributed to roughness of the
edges of the holes. The rectangular holes were machine
stamped while the round holes were cut using a hole saw.
It is interesting to note that Kaminski (1986a)
presented the following equation for static pressure
differences:
2
(39) AP = C p V / 2g
s D m c
where
Ap = Static pressure difference across a given
s configuration
C = Static pressure coefficient based on highest the
D highest average velocity for a given configuration
V = Highest average velocity
m
and the other terms are as previously defined.
For turbulent flow, Reynolds Number greater than
10,000, through a squareedged orifice the static pressure
coefficient was reported as C = 2.80.
D
If Equations 35 and 37 are rearranged into the form of
Equation 39, the coefficients of static pressure are
respectively, 2.74 and 3.03. Therefore, the pressure loss
across a single fiberboard can be approximated by the
pressure loss across an orifice or orifice plate.
No such comparison of the pressure loss across two
fiberboards, separated by a fixed distance, could be made
with the typically presented expansion and contraction
configurations reported in various fluid flow text and data
books.
Chau et al. (1983) experimentally evaluated the
cumulative effect of air flow resistance of the vent holes
and the oranges. Measurements were made of the pressure
drop across oranges stacked on fiberboard with 4 rectangular
vent holes which simulated the fruit with one side of a
carton. The pressure drop was also measured across oranges
stacked between 2 fiberboards with 4 rectangular vent
holes of the same size and location on each fiberboard.
All the tests were repeated using 4 different stacking
patterns and 4 different orange sizes.
Chau et al. (1983) fitted the experimental data using
a modification of Equation 31 which did not normalize the
pressure drop to a unit depth of fruit
b
(40) AP = a' V
This was necessary because the measured pressure drop was a
combination of the pressure change across the oranges, which
was a function of fruit depth, and the pressure change across
the vent holes, which was not a function of fruit depth.
For the size 80 oranges stacked in a square staggered
pattern, with a bulk depth of 0.33 meters and a porosity of
0.430, the constants in Equation 40 for oranges with one
carton side were found to be
4
a' = 1.85 x 10 and b = 1.98, while for oranges with
two carton sides the constants were found to be
4
a' = 2.46 x 10
and b = 1.95.
Chau et al. (1983) noted that most of the pressure loss
was caused by the resistance of the vent holes. The
pressure drop across a carton of fruit was always greater
than the sum of the pressure drops across each side of the
carton plus the pressure drop across the fruit. Wang and
Tunpun (1969) reported similar results.
According to Chau et al. (1983), the pressure drop
across the sides of the carton was expected to be higher
than the sum of the pressure drops across the sides of the
carton and the fruit. The fruit disturbed the flow past the
vent holes since there was always an amount of obstruction
by the fruit. Further, the air flow pattern through the
oranges in a carton with vent holes was much different from
the uniform flow through bulk loads of oranges even for the
same approach velocity. Near the area of the vent holes for
flow through oranges in cartons, high localized velocities
occurred. Since the pressure loss through the vent holes
was proportional to the velocity squared, approximately, the
high localized velocities produced a higher overall pressure
drop across the vent holes and through the fruit near the
vent holes.
All the necessary input data to use the ANSYS finite
element program to model the experimental work of Chau et
al. (1983) are known.
Extensive modeling using twodimensional elements,
like that employed previously in modeling the work of
Segerlind (1982) and Khompis (1983), was carried out but
was not be reported here since the work of Chau et al.
(1983) was threedimensional. The twodimensional work was
very helpful initially in providing necessary direction
without the additional computer time required when working
with threedimensions. The twodimensional modeling results
were essentially the same as the threedimensional modeling,
which is discussed in detail below.
Modeling the experimental work of Chau et al. (1983)
involved more than one application of the capabilities of
the ANSYS finite element solution package. First, the
flow through the oranges in bulk was easily modeled as
illustrated previously for the grain bin problems of
Segerlind (1982) and Khompis (1983). The boundary
conditions at the crosssection of the entrance and exit
were easily applied for this uniform flow situation. The
second application dealt with the reduced boundary
conditions as a result of the entrance and exit vents of the
oranges packaged in cartons. The third application involved
modeling the orificetype loss through the vents of the
carton walls. Finally, the combination of the second and
third applications was applied to study the combination
losses which were reported by Chau et al. (1983) and Wang
and Tunpun (1969).
The physical problem selected for solution was the
experimental setup shown in Figure 9 and described above.
An eight node 3D isoparametric thermal solid element,
which allowed modeling of nonlinear steadystate flow
through a porous media, was used to model the pressure loss
through oranges in bulk and oranges in cartons, without
consideration of the pressure loss through the vent holes.
This element could be formed into prism or tetrahedron
shaped elements by defining duplicate node numbers.
However, only 3D rectangular elements were needed. The
node and element locations are shown in Figure 10. The
model consisted of 320 nodes and 196 elements. The nodes
were placed every 8.1 cm. The boundary conditions were
specified at the inlets and outlets corresponding to data
reported by Chau et al. (1983). For the situation with
oranges in bulk, the pressure or flow rate across the inlet
and outlet crosssections was specified. For the situation
with oranges in cartons, the pressure or flow was specified
at nodal locations corresponding to the four inlet and exit
vent hole locations of the experimental setup used by Chau
et al. (1983). The inlet and outlet locations were point
sources and the velocity and pressure varied in the x,
y and zaxis. As noted in previous models, the noflow
or impermeable boundary conditions are automatically
enforced on the boundary when no pressure values were
specified.
In order to use ANSYS to model the pressure loss
through the vent holes in the fiberboard, a new element type
was used. The element was the hydraulic conductance element
(DeSalvo and Swanson 1983). This uniaxial two node element
had the ability to transmit flow between its nodal points.
I 4
4
56.7cm
3o AIR rNou RSIST4NC rat amReS IN LUK
2 3
Figure 10. The node and element locations used to model 3D
oranges packed in bulk and in cartons.
