MATHEMATICAL MODELING OF HEAT TRANSFER IN THE COOLING OF FRUIT IN CLOSED CONTAINERS By TOMAS BAZAN \ 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 This dissertation is dedicated to my parents, Tomas Bazan Betancourth and Juana Maria Bolanos de Bazan, for their constant love and support. ACKNOWLEDGEMENTS I would like to express my appreciation to Dr. Robert B. Gaither, my major professor, for his support and guidance since my arrival at University of Florida. He provided valuable suggestions during the development of this study. I would also like to thank Dr. Khe Van Chau, cochairman of my dissertation committee, for his insight and support during the entire research process. I would like to acknowledge the valuable contributions to my education throughout my time here at the University of Florida made by Dr. Herbert A. Ingley, my academic advisor and professor. Dr. Ingley provided much excellent academic advice. I would like to thank Dr. C. Direlle Baird for his recommendations and guidance and Dr. William Lear for his valuable teachings and suggestions during the completion of this work. Dr. Lear also helped reinforce my background in the thermal sciences through the various courses he taught me. I would like to extend my sincere thanks to Mr. Loren Miller, without whose technical support, advice and guidance this research would not have been possible. I would like to express special gratitude to the Technological University of Panama for permitting me the leave of absence necessary to obtain this doctoral degree and to the Ful bright-LASPAU Foundation for providing both financial and administrative support during my three and one-half years in Florida. TABLE OF CONTENTS page ACKNOWLEDGEMENTS iii LIST OF TABLES vi LIST OF FIGURES vii LIST OF SYMBOLS x ABSTRACT xiii CHAPTERS 1 INTRODUCTION AND OBJECTIVES 1 Statement of the Problem 4 Research Objectives 5 2 LITERATURE REVIEW 6 Physiological Factors Relevant to the Cooling Process of Fruit 6 Mathematical Models for Heat Transfer in Bulk Cooling of Products 8 3 MATHEMATICAL MODEL FOR SIMULATION OF HEAT TRANSFER IN THE BULK COOLING OF FRUITS 12 Numerical Grid Generation 13 Model Equations 21 Numerical Algorithms 31 4 EXPERIMENTAL SETUP AND TEMPERATURE MEASUREMENT PROCEDURES 38 Description of Cooling Facility 38 Data Acquisition System 41 Reheat Facility 42 Temperature Measurement Procedures 42 i v 5 RESULTS AND DISCUSSION 48 Comparison of Experimental and Simulation Results . . 49 Results for Cubic Packing Pattern 61 Results for Random Packs 67 Effects of Room Air Velocity and Packing Arrangement 69 Effects of Uncooled Walls 80 6 CONCLUSIONS AND RECOMMENDATIONS 85 APPENDICES A MAIN PROGRAM: THREE-DIMENSIONAL TRANSIENT TEMPERATURE RESPONSE SOLUTION IN A FRUIT BIN 87 B INPUT PROGRAM: PROPERTIES AND DIMENSIONS INPUT. 3-D TRANSIENT TEMPERATURE DISTRIBUTION SOLUTION IN A FRUIT BIN 117 C INPUTS TO THE NUMERICAL MODEL 122 BIBLIOGRAPHY 123 BIOGRAPHICAL SKETCH 126 v LIST OF TABLES TABLES page 3.1 Forced convection constants 3.2 Types of position fruit surface and air nodes according to in the bin 4.1 Physical and thermal properties of Sunny tomatoes . . . 44 4.2 Measured weight and diameter of Sunny tomatoes . . . . 46 4.3 Measured contact areas vi LIST OF FIGURES FIGURE gage 1.1 Diagram of air path during room cooling 3 3.1 Top-side view of a cubic packing (CP) pattern of spheroids. . 14 3.2 Face-centered cubic (FCC) packing of spheroids 15 3.3. Top-side view of cubic volume elements to model 3.D fruits in bulk 17 3.4. Chau-Bellagha Network of fruit nodes inside a cubic volume. . 18 3.5 General flow diagram for the numerical calculation of the bin temperature response 32 3.6 Typical flow diagram of an array generating subroutine. ... 35 3.7 Flow diagram of the SLE numerical solution algorithm 36 4.1 Cooling unit 4.2 Reheat facility 43 4.3 Thermocouple locations in the tomato container 45 5.1 Fruit locations in the experimental container 50 5.2 Comparison between experimental and simulation results for fruit 15 and fruit 9 51 5.3 Comparison between experimental and simulation results for the temperature response along a diagonal axis during room cooling of a bin of tomatoes at V=18 m/min 52 5.4 Comparison between experimental and simulation results for the temperature response along a diagonal axis during room cooling of a bin of tomatoes at V=120 m/min 53 5.5 Comparison between experimental and simulation results for the temperature response along a diagonal axis during room cooling of a bin of tomatoes at V=240 m/min 54 VI 1 5.6 5.7 5.8 5.9 5.10 5.11 5.12 5.13 5.14 5.15 5.16 5.17 5.18 5.19 Comparison between experimental and simulation results for the temperature response along the central vertical axis during room cooling of a bin of tomatoes at V=18 m/min. ... 55 Comparison between experimental and simulation results for the temperature response along the central vertical axis during room cooling of a bin of tomatoes at V=120 m/min ... 56 Comparison between experimental and simulation results for the temperature response along the central vertical axis during room cooling of a bin of tomatoes at V=240 m/min ... 57 Comparison between experimental and simulation results for the temperature response of tomatoes packed in random and straight column arrangements at V= 1 20 m/min 58 Comparison between experimental and simulation results for the temperature response of tomatoes packed in random and straight column arrangements at V=240 m/min 59 Comparison of simulation results of temperature gradients development inside and at the surface of tomatoes during room cooling at V=18 m/min 62 Comparison of simulation results of temperature gradients development inside and at the surface of tomatoes during room cooling at V=240 m/min 63 Simulation of the effect of increasing fruit contact area on cooling times of tomatoes at V=1 20 m/min 64 Effect of room air velocity on the cooling time of a fast cooling fruit Effect of room air velocity on the cooling time of a slow cool ing fruit Simulation of the effect of room air velocity on the mass average temperature response of random packs of tomatoes. . . 72 Simulation of the effect of room air velocity on the mass average temperature response of straight column packs of tomatoes 7-3 Simulation of the effect of room air velocity on the mass average temperature response of tomatoes packed in a body-centered cubic (BCC) pattern 74 Comparison of mass average temperature with fastest and slowest cooling fruit temperature responses in a random pack of tomatoes 75 5.20 Simulation of mass average temperature responses of different packing arrangements of tomatoes at V=18 m/min. . . 76 5.21 Simulation of mass average temperature responses of different packing arrangements of tomatoes at V=1 20 m/min . . 77 5.22 Simulation of mass average temperature responses of different packing arrangements of tomatoes at V=240 m/min . . 78 5.23 Simulation of mass average temperature responses of containers of tomatoes with exposed and covered top and bottom walls 5.24 Simulation and test results of the fastest and slowest temperature responses in random packs of tomatoes inside containers with exposed and covered walls . 83 IX LIST OF SYMBOLS Ac: average contact area An: average conduction area between fruit nodes As: fruit surface area per unit volume of cubic element c: node thermal capacitance D: matrix right-hand side vector Df: fruit diameter D1 : container wall length g: acceleration of gravity ha: fruit surface convective coefficient hc: thermal conductance for diffusion through fruit K: specific permeability k: thermal conductivity L^: number of cubic volumes in the ^-direction NCP: number of contact points PD: fruit packing density Qr: fruit heat of respiration rfa: fruit average radius T: node temperature Tas: air temperature at convective surface Ta: temperature in the air phase Tf: fruit field temperature x room temperature time t: V: room air velocity V2: convective velocity Greek symbols aav: air isometric thermal diffusivity aap: a^r isobaric thermal diffusivity P: volume expansion coefficient A: first order forward difference operator A1 : space increment Ar: thickness of fruit spherical element At: time increment AV: fruit element volume V: first order backward difference operator 8: second order central difference operator e: porosity A: thermal conductance Aa: thermal conductance in the air-wall boundary li: dynamic viscosity direction variable p : density Pas: air density at fruit surface Pf : fruit density 0: wall identification index w: wall index operator xi TV- energy equation operator Subscripts a: air as: air at the fruit surface f: fruit h: node index for z-direction i: node index for x-direction j: node index for y-direction m: one-dimensional numbering index s: surface w: wall Superscripts n: at time n o: at initial time 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 MATHEMATICAL MODELING OF HEAT TRANSFER IN THE COOLING OF FRUIT IN CLOSED CONTAINERS By TOMAS BAZAN August, 1989 Chairman: Dr. Robert B. Gaither Cochairman: Khe Van Chau Major Department: Mechanical Engineering Postharvest technology known as room cooling to retard physiological deterioration of perishable commodities has become increasingly difficult to carry out as the procedure of tight stacking on pallets with little or no ventilation has become more common in practice. Yet, several operational and economical advantages make this method of cooling convenient for a variety of applications in this country as well as in less developed nations. There is considerable interest at the present time in developing an efficient mathematical model that contributes to the understanding of basic heat transport mechanisms relating to the process of room cooling of fruit. It is envisioned that such a model could be effectively used xi i i in the prediction of the three-dimensional temperature response and room cooling times for fruit bins stacked in commercial packing arrangements. In this study, a mathematical model was developed which accurately predicts the three-dimensional temperature response during the room cooling of a confined bin of spherical fruit. The model accounts for the contribution of heat conduction through fruit contact areas as well as convective buoyancy effects on the removal of field and internally generated heat from fruit. Experiments to validate the model were conducted with tomatoes arranged in pattern and random packs inside a closed, uninsulated container. Very good agreement between experimental and simulation results was obtained for the various room air speeds and packing arrangements used in packing houses. The results depict a more significant contribution of heat removal across contact areas for fruit located in the core of the bin than near walls. The influence of natural convection velocities was found to cause the displacement of the maximum temperature location from the center of the bin upwards as cooling advanced. Packing arrangements with high numbers of contact points exhibited faster cooling rates. Adeguate ventilation and room air velocities cut cooling significantly. Simulation of fruit containers stacked on top of each other yielded cooling times as much as 40% longer than those of containers with all walls exposed to room air. xiv CHAPTER 1 INTRODUCTION AND OBJECTIVES Fruit, as living organisms, undergo a continual process of deterioration from harvest on. The market value of this produce is consequently associated with the prompt and effective control of the variables that influence physiological decay. Of all the different environmental conditions, temperature has been identified as the single most important variable affecting the levels of tissue injuries, wilting, moisture loss and microorganism invasion which perishable commodities exhibit. This interrelation can be better understood by recognizing that harvest temperatures are often the most favorable to the growth of rot microorganisms (Mitchell et al . , 1972). On the other hand, temperature is not only a factor in the formation of ethylene, but also to the extent that the concentrations of this gas, oxygen and carbon dioxide affect the fruit. In order to effectively retard the deterioration process, the temperature of the harvested commodities must be reduced as early and quickly as possible to the desirable storage conditions in what is known as the cooling or precooling process. It is known that temperature increases of 10.0 C above optimum conditions can increase the rate of deterioration by two- to three-fold (Kader et al . , 1985). Among the several cooling techniques available, the method known as room cooling is widely used with a variety of fruit. The typical 1 2 arrangement involves stacking the harvested produce in rooms cooled by an air stream which is discharged horizontally just below the ceiling and circulated throughout the cooling room (Figure 1.1). In this arrangement, the air tends to take the path of least resistance. In order to obtain an homogeneous cooling of the produce, it is advisable to provide both for enough spacing between pallets, and air vents through the container walls. Velocities of at least 61 to 122 meters per minute are usually recommended to achieve best results (Kader et al . , 1985). The slow heat-removal rates characteristic of the room cooling method are attributed to the fact that the energy flow from the center of the bins to the container surfaces takes place largely by conduction and convective buoyancy. Furthermore, with the growing trend of handling large volumes of produce, in many cases containers have to be tightly stacked on pallets. The effect of this practice is not only to lengthen the heat diffusion path, but also to reduce the area of contact between the air streams and the container surfaces. Yet the operational simplicity, low initial cost, and advantage of cooling and storing in the same place make this method convenient for a variety of fruit. The review of literature reveals that extensive studies have been conducted on the mechanisms of heat generation, diffusion, surface convection and radiation (Chau and Bellagha, 1985), as well as the thermal properties (Baird and Gaffney, 1980) that control the heat- removal process from a single fruit. 3 Diagram of air path during room cooling. 4 Statement of the Problem Several numerical models have been implemented to analyze the method known as forced-air cooling. In this method, unlike room cooling, commodities are cooled by moving air through rather than around containers. Consequently, one-dimensional forced-bulk heat transport and surface transpiration prevail over heat diffusion as the mechanisms of energy flow (Chau et al . , 1984; Bakker-Arkema and Bickert, 1966). Instead, within the partially ventilated bulk stores typical of room cooling, three-dimensional diffusion across the air and fruit phases and natural convection at the fruit surface seem to be the dominant factors. No numerical model was found in the literature review which addresses this problem. Simulation models which could aid in understanding the heat transfer mechanisms within bulk stores of fruits could also provide valuable technical information on more efficient ways to cool these commodities. In this study, a mathematical model is developed which attempts to determine the three-dimensional temperature history during the bulk cooling of a confined bin of spherical fruits. The model accounts for the contribution of diffusional transport, contact area as well as convective buoyancy effects to the removal process of field and internally generated heat. Experimental tests to validate the numerical model were conducted with tomatoes arranged in a closed, uninsulated box placed inside a controlled -environment, walk-in cooler. Temperatures at various locations within the air phase and inside fruit were monitored. The velocity of the cooling air stream passing by the container sides was also recorded. 5 Research Objectives The main objectives of this research were to: 1. Develop a model that can provide better understanding of the contribution of the basic energy transport mechanisms related to the process of room cooling of spherical fruits. 2. Validate the numerical model with experimental data obtained from tests of bins of tomatoes conducted at different air stream velocities and packing patterns. 3. Utilize this model in the prediction of three-dimensional temperature responses and room-cooling times for bins of spherical fruits stacked in various commercial packing arrangements. CHAPTER 2 LITERATURE REVIEW Physiological Factors Relevant to the Cooling Process of Fruit The market value of fruit is closely associated with the prompt and effective control of the variables that influence physiological decay of these commodities right after harvest. Adequate temperature control in postharvesting operations has been identified as the single most important factor in containing the levels of tissue injuries, wilting, moisture loss and microorganism invasion which perishable commodities exhibit. Temperature is also an important factor in controlling the formation of ethylene and the extent of the damage that concentrations of this gas, oxygen and carbon dioxide do on the fruit. This interrelation between fruit deterioration and temperature control can be better understood by recognizing that harvest temperatures are often the most favorable to the growth of rot microorganisms (Mitchell et al . , 1972). In order to effectively retard the deterioration process, the temperature of the harvested commodities must be reduced as early and quickly as possible to the desirable storage conditions in what is known as the cooling or precooling process. It is known that temperature increases of 10.0 C above optimum conditions can significantly increase the rate of deterioration of perishable commodities (Kader et al., 1985). 6 7 The arrangement known as room cooling is widely used with a variety of fruit because of the simplicity of its design and operation. This arrangement typically involves stacking the harvested produce in rooms cooled by an air stream which is discharged horizontally just below the ceiling and circulated throughout the cooling room. Another advantage of the room-cooling arrangement is that the peak loads that its operation imposes on refrigeration systems are considerably smaller than with faster cooling methods (Mitchell et al . , 1985). In the room-cooling arrangement, the air tends to take the path of least resistance. Consequently, in order to obtain an homogeneous cooling of the produce, it is advisable to provide both for enough spacing between pallets, and air vents through the container walls. When sufficient air movement is provided to avoid the formation of warm pockets in the room, the exposed container surfaces are brought to air stream temperature in a short time (Guillou, 1960). Velocities of at least 61 to 122 meters per minute are usually recommended to achieve best results (Kader et al., 1985). Room-cooling arrangements are characterized by low cooling rates. This is attributed to the fact that heat removal from the center of the bins to the container surfaces takes place largely by conduction through contact areas and heat transport upwards by natural convection velocities. Furthermore, as a result of handling large volumes of produce, containers are often tightly stacked on pallets. Two undesireable consequences of such stacking are the lengthening of the heat diffusion path and the reduction of the area of contact between the air streams and the container surfaces. In spite of these consequences, 8 operational simplicity, low initial cost, and the advantage of cooling and storing in the same place make this method convenient. Mathematical Models for Heat Transfer in Bulk Cooling of Products Extensive studies have been conducted on the mechanisms of heat generation, diffusion and surface convection and radiation (Chau and Bellagha, 1985); as well as the thermal properties (Baird and Gaffney, 1980) that control the heat removal process from a single fruit. Much research has been devoted to modeling heat and mass transfer from bulk stores of produce. For the most part, this research has been concentrated on a process based on cooling produce by direct exposure to the air stream. This method is commonly referred to as forced-air cooling. In this forced-air cooling arrangement, unlike room cooling, one-dimensional forced-bulk heat transport and surface transpiration prevail over heat diffusion as the mechanisms of energy flow (Chau et al., 1984; Bakker-Arkema and Bickert, 1966). Rowe and Claxton (1965) studied heat and mass transfer from a single sphere to a fluid flowing through an array of spheres. A set of basic correlations for the Nusselt number of a single sphere at different Reynolds numbers was presented. The effect of the porosity of the medium was incorporated into this set of correlations. Lerew (1978) implemented a model to simulate heat and mass transfer from bulks of potatoes exposed to a cooling air stream. The model incorporated transpiration and convective effects at the surface of the product as well as the contribution of diffusion of mass and heat during the cooling process. 9 Chau and Bellagha (1985) developed an explicit finite difference method for predicting heat and mass transfer during forced-air cooling of both single and bulk tomatoes. Their model accounted for heat removal from convection and transpiration at the fruit surface, moisture loss, and the heat of respiration. A massless fruit surface node was considered in order to make the stability restrictions on the explicit formulation less strict. Talbot (1987) applied a three-dimensional porous media model to the simulation of temperature responses during the forced cooling of bulks of oranges. The model considered the effect of porosity in the three-dimensional temperature distribution throughout the medium during the cooling process. Extensive research has also been conducted on the contribution of natural convection velocities to the transport of heat from unventilated two-phase media. In most of the models developed, Darcy's law is used to simplify the momentum equation. The majority of models in this area can be categorized by the number of dimensions considered inside the medium and by whether temperature gradients are considered inside the solid phase. Beukema and Bruin (1983) developed a three-dimensional model for heat transfer from confined bins of fruit treated as a porous medium. The influence of natural convection effects in the heat-removal process and temperature distribution inside the container was thoroughly examined. Heat generation due to fruit metabolism was assumed to be uniformly distributed throughout the two-phase porous medium. The model 10 did not account for the effects of transpiration and natural convection in the heat removal process at the surface of products. Johns and Lawn (1985) presented a series solution for heat transfer between a sphere and a surrounding stationary medium enclosed in a concentric insulating sphere. The analysis revealed that a short time after the start of the heat transfer process the Nusselt number approaches a constant value which depends only on the ratio of the radii of the two spheres. Trevisan and Bejan (1985) presented a comprehensive study of the influence of natural convection inside a porous layer with both mass and heat transfer from the side walls and adiabatic conditions at the top and bottom. The analysis evaluated the effects of temperature and concentration gradients in the generation of natural circulation in the med i urn . Mural idhar et al . (1986) conducted an experimental and numerical study of natural convection flow in a porous material contained between two concentric isothermal cylinders. By using Darcy's law and the Boussinesq approximation, the Nusselt number of the variable porosity problem approximated that of the uniform porosity case. Reddy and Mulligan (1987) used macroscopic continuum mechanics to derive the governing equations for heat and mass transfer in a porous medium with internal generation. Darcy's law was used to simplify the governing equations. Mueller (1988) presented an analytical solution for the one- dimensional temperature distribution in a volumetrically heated porous 11 medium. The solution method was based on the use of the Green's function to simplify the governing equations. The assumption of one-dimensional temperature gradients during the cooling process of confined bins products is a feasible simplification in well-insulated containers (Beukema et al . , 1982). Conversely, when there are considerable rates of heat transfer at the side walls of containers, three-dimensional temperature differences are likely to develop and remain during the cooling process. Within the partially ventilated bulk stores typical of room cooling, three-dimensional diffusion across the air and fruit contact areas and heat transport by natural convection velocities are expected to be the dominant factors in the heat-removal process. No numerical model was found in the literature review which addresses the combined effect of these mechanisms in the cooling process of unventilated containers of produce. Simulation models which could aid in the understanding of the heat transfer mechanisms within bulk stores of fruits could also provide valuable technical information on more efficient ways to cool these commodities. CHAPTER 3 MATHEMATICAL MODEL FOR SIMULATION OF HEAT TRANSFER IN THE BULK COOLING OF FRUITS In this chapter a mathematical model is developed which attempts to determine the three-dimensional temperature history of the bulk cooling of a confined bin of spherical fruits. In addition to the conventional approach of diffusional transport as the controlling mechanism, the model accounts also for the contribution of convective buoyancy effects to the removal process of field and internally generated heat. The grid generation scheme used to divide the fruit and air phase spaces into differential elements, and thereafter to number and identify such elements, is discussed first. Further discussion in this section is devoted to the flexibility and advantages obtained from this scheme when dealing with the governing equations and boundary conditions. The following section deals with the derivation of the model equations for the air phase and fruit elements. Second-order central differences in an implicit format are used to express the energy balances for the interior nodes. Special operators are introduced into the energy equations for internal nodes to add flexibility for the treatment of various types of boundary conditions. This synthesizing strategy of the energy equation proves rather useful later during the formulation of the solution algorithms. 12 13 Finally, the numerical algorithms employed to solve the set of implicit energy equations for the three-dimensional temperature distribution in the fruit bin are presented. The solving procedure consists of a variation of the Gauss-Seidel Indirect Method which permits an independent treatment of different types of node equations. Additional operators are introduced with the objective of synthesizing the numerical routines and thus saving computing time. An analysis of the nature of the energy equations with regard to stability and convergence problems is considered in the discussion of the convenience of using the implicit over the explicit formulation. Numerical Grid Generation The system considered in the formulation of this model is the three-dimensional rectangular container of spherical fruits arranged in the straight column pattern shown in Figure 3.1. This pattern, also referred to as cubic packing structure (CP), is obtained by stacking spheres of adjacent layers on top of each other, yielding a theoretical packing density (PD) of about 46% and six contact points (NCP). Two other well-known packing patterns are the body-centered cubic structure (BCC) with NCP=8, and the face-centered cubic structure (FCC), shown in Figure 3.2, with NCP=12. The latter pattern is generated by positioning fruits of a given layer in the spaces formed by the fruits of the layer immediately below. Although the fee pattern could theoretically provide the highest packing densities (68 to 74%), random packing is more commonly used in the fruit industry as it facilitates automatization and handling of 17 18 Figure 3. -NON -CAPACITANCE SURFACE NODE ( h = 2 + 4k ) -INTERIOR NODE ( h=3 +4k ) •AIR PHASE NODE (h=1+4k) -CENTER NODE (h=4+4k ) k =0.1.2 (M/4 - 1) . Chau-Bel lagha Network of fruit nodes inside a cubic volume. 19 3. The interior nodes are located in the center of mass of each element. This improves accuracy (Bellagha, 1984). 4. One-dimensional heat flow in the radial direction is assumed inside the spherical fruits. Because of the very small temperature gradients expected inside the fruit, only three nodes are used in the fruit phase. The resulting grid, as illustrated in Figure 3.4, consists of one air node for the air phase and surface, and interior and center nodes for the fruit phase. An indexing scheme needs to be devised now in order to implement an algorithm to identify and locate air phase and fruit nodes throughout the two-phase system. The process is as follows: 1. The number of cubic volumes in each of the three directions can be estimated by taking the nearest integer of D1 ^ Lr 2rfa £*1,2,3 where £= numbering sequence for the x, y and z-direction respectively; rfa= fruit average radius; and Dl^= wall length in ^-direction. 2. Because of the three-dimensional, two-phase nature of the system under consideration, a set of four indexes (i,j,k,l) would normally be required in order to specify the position and phase of an element. In the system under consideration, this set can be reduced to a set of three indexes (i,j,h). This becomes possible by recognizing that the position of a given type of node, i.e., surface node, within a cubic 20 volume is the same regardless of the location of the cubic volume in the bin. The h-index in this set is introduced to specify the type and vertical position of a given node in the container by using the following formula: h=a+4k a=l ,2,3,4 k=0 ,1,2, . . .L3 where a= integer number that specifies node type following the convention illustrated in Figure 3.4, and k= number of the horizontal layer (counted from the bottom) where the element is located. Air phase and fruit nodes are located and identified by type by specifying their grid indexes (i,j,h) with respect to a set of three-dimensional cartesian coordinates centered in the north-east bottom corner of the test box (Fig 3.3). A transformation is introduced to switch from the three- dimensional indexing notation (i,j,h) to a one-dimensional index (m). This transformation becomes useful when representing and manipulating variables in intermediate matrix operations. Following an ordering sequence in which the elements are counted first in the h-direction, then in the j-direction and last in the i-direction, the transformation obtained is m=h+4L3[L2(i-l)+(j-l)] . (3.1) The second right-hand term in this equation is the total number of nodes counted, following this ordering sequence before the cubic volume where a given node is located. 5. The length of a cubic volume is estimated as A1 (s)= 2rf(s) 21 s=l,4,8. . .m. A 0-index is introduced to identify the six box walls (Fig. 3.3). Model Equations In the derivation of the model equations the following assumptions are made: 1. Heat generation due to respiration is a function of fruit element temperature. 2. Variations in fruit thermal properties during the process can be neglected. 3. Buoyancy effects due to natural convection can be represented by the Boussinesq approximation applied to Darcy flow through porous media. 4. Water vapor saturation in the air phase of confined bins occurs very early in the cooling process, making heat transport from transpiration at the fruit surface negl igible. 5. Natural convection at the fruit surface is the prevailing mechanism of heat removal from inside the product out to the air phase. The natural convection coefficient is constant over the product surface. 6. Heat from the bulk is transferred by conduction and convective buoyancy motion through the air phase and by 22 conduction across fruit contact areas. Horizontal mass transfer is negligible. 7. For the diffusion of heat through contact areas between fruit, thermal contact resistance is neglected. As an approximation, only the thermal resistance of the layer of fruit between the surface and the interior node is considered. The velocity induced by convective buoyancy in air phase of a porous medium can be expressed by Darcy's law (Beukema and Bruin, 1983): Kg Vz~ [Paz-P al]- (3.2) where subscripts 1 and 2 denote conditions in adjacent vertical layers. The specific permeability (K) is calculated using the corrected form of the Kozeny-Carmen equation presented by Fahien (1983): e3Df2 K_ 180 ( 1 - e ) 2 . Introducing the Boussinesq approximation, Pa2~Pal ~ Aa$a(Tal-Ta2) » into equation (3.2) gives the final expression for the convective velocity, V = /j _y \ yz .. v 1 al a2' ' Using a differential formulation the velocity equation becomes V n ijh ~ Kg/^a M n n ( * i jh'^i jh+4) ' (3.3) With the model assumptions previously adopted, the equation governing conservation of energy for air phase elements is 23 f jj. + lz 5Ia = aav3t “ap dz 3X2 d‘T, 32T 32T a? + a? + az2 (3.4) under initial and boundary conditions: Ta(x,y,z,0)=T°a V2(x,y,z,0)=0 (3.5) where T0= room air temperature, and M0) = thermal conductance of 0-wall, including inside and outside air films. The two terms in the left-hand side of equation (3.4) represent the energy storage rate and convective energy transport associated with the air phase element respectively. The terms in the right-hand side of equation (3.4) represent three-dimensional conductive transport of heat. For air nodes away from walls, these three terms can be expressed in implicit formulation as follows: Discretizing the left-hand terms in equation (3.4), and inserting equation (3.6) into it, the implicit representation of the energy equation for an air node away from the walls gives a 24 (6x+6y+6z) n+i . , ? L 1 i i h J Al2 ‘ ,Jh (3.7) where 6^= second order central difference operator acting on ^-direction. As an illustration of the discretization of the right-hand terms of equation (3.4) for air nodes near walls, the case of an element in contact with the East wall (0= 1) and the South wall (0=3) can be examined. Since this element is not in contact with any walls in the z-direction, the derivative in this direction can be discretized as fol 1 ows The derivatives in the directions normal to walls in contact with the element become jh+4_2Tj jh+Ti jh-4) 6 A1 (3.8) _ Aa, dx2 kaAl (3.9) 25 and 32T Aa, (T. . , . -T. * , __3 (TO"Ttjh) + ' uh' 3/ k,41 ,Jh a ^a3 k,Al AT (To-TiJh) - Vy[TiJh] AT (3.10) where A^= first order forward difference operator acting on ^-direction, and V^= first order backward difference operator acting on ^-direction. The first term in the right-hand sides of equations (3.9) and (3.10) corresponds to heat conducted from the side of the element in contact with the wall. The second term represents heat conduction in the side of the element away from the wall. For elements in contact with just one wall, two of the derivatives in the right-hand side of equation (3.4) are represented in the format of equation (3.8). The third derivative is expressed either as equation (3.8) or (3.9) depending on the wall with which the element is in contact. In the case of elements in contact with three walls, all three derivatives follow the format of equations (3.9) and (3.10). In order to achieve simplifications in numerical algorithms, equation (3.4) is expressed in a general format that permits the treatment of elements both in contact with and away from walls. This general expression is obtained by incorporating the format in equations (3.9) and (3.10) into equation (3.7), and introducing a set of special operators, giving: 26 e rTn+1 ] = (Mx+My+Mz) rTn+1 n ^ L 1jh J (VaA-V^x+V^y-V^y+riatK-riasVz) rxn+l , W [ ijh 3 , ha. ijhAs, i jh rT n+l f»l -I , 26Aa(hwf(0)) _ n+1 k L ijh+r ijh J + 0=1 j^jj I 'o" 'ijh J with initial conditions Tiih = Tfa V° = 0 i jh (3.11) where r?f(£)= direction index operator defined as ' nrtt)- ' 1 nodes away from walls normal to ^-direction 0 nodes next to wall normal to ^-direction, r/a(0)= wall index operator defined as 1 nodes next to 0-wall 0 nodes away from 0-wall, and hwf(0): wall conductance operator defined as (3.12) r)M) = (3.13) huf (0)= 1 nodes next to 0-wall , Aa(D=Aa0 0 nodes away from 0-wall Aa(0)=0. (3.14) By introducing these three operators, flexibility is added to a single general energy equation in order to make it applicable for handling both interior and boundary nodes. In particular, the hwf operator permits for each wall the simulation of different boundary 27 conditions such as surrounding air speed, adiabatic barriers, and wall conductance. This synthesizing strategy proves useful as well in simplifying the numerical algorithms and saving computer time. A procedure similar to the foregone derivation leads to the following energy conservation equation for fruit interior node: ^ rTn+1 jn -| ^nrm,h,h-l'^c) rTn+l Tn+1 °f4t ,jh‘ ,ih J ' A n,m,h,h+l A Vm hArm , , , m,n m,h,h+l [Tn+1 L i.i rj -T 1 L 1 ijh-1 1 ijh i Pf..Q"r jh+l" ijh Tfl1 ] + 2 Af (CJj-) Ac _ _n+l _ h A n+1 + 0=1 — [T0-r;Jh] +_Lff(»7flfix+My+Mz)[Tijh ] kfAV Mc kfAV kfAV ^al^x'^a6^x+^a2^y_?la3^y+,la4^z"^a5^zHTi jh ] - 0. (3.15) The fourth right-hand term in this equation represents heat conduction through wall-fruit contact areas. The next to last right-hand term represents heat conduction across contact areas between fruit away from walls. The last right-hand term corresponds to heat conduction through contact areas between fruit in contact with walls. Following the previously described model of Chau and Bellagha (1985), the energy equation of the fruit surface node is expressed as MAn,m,h+l,h-Ac) Ar n,m,h+l,h [CrO + (3.16) Similarly, for the fruit center node, the energy equation is 1 rjn+1.Tn 1 _ ^ afAt L 1 AV . Ar ~ ' m,h m,h,h-l ‘n,m,h,h-l .. n+1 Tn+1 .. AfimQ ‘■ijh-1-* ijh -I + (3.17) with initial condition 28 where Ar= radial distance between adjacent nodes, A V= volume of node, Ac= individual contact area, An= average conductive area between adjacent nodes, As= fruit surface area exposed to convection, and Qnr= heat of respiration at time step n. Implicit representations of finite difference equations can be shown to be unconditionally stable for all mesh sizes (Hornbeck, 1975). As a result, the restrictions to small space increments imposed by the numerical instability associated with the explicit method are removed. More specifically, for an explicit solution to be numerically stable, the following general criterion of compliance with the second Law of Thermodynamics, as presented in Welty (1978), must be met: At < Cm m — p|V (3-18) where Cm= thermal capacitance of node m, and \,p= thermal conductance between node m and surrounding node p. Equation (3.12) shows that the existence of very low thermal capacitance nodes in the mesh, i.e., air phase nodes in the system under consideration will also dictate the use of very small time steps for the entire network. In most explicit formulations such time steps are smaller than those required to hold truncation errors to acceptable levels. 29 Since truncation errors of the explicit and implicit methods are comparable for the same finite difference expression order, implicit schemes are more desirable when stiff instability constraints are imposed on the explicit representation. The net effect is that the computer time needed to solve the large sets of implicit equations is offset by that needed to solve the long number of step calculations required to avoid instability in the explicit scheme. A convective coefficient based on natural buoyancy effects at the product surface is used to describe heat transfer from fruit to air. The following correlation for spheres in fluids with the Prandtl number (Pr) close to unity is presented by Holman (1981): 1/4 Nu0 = 2 + 0.43 RaD where NuD= Nusselt number based on sphere diameter, and Ra0= Rayleigh number based on sphere diameter. An average natural free-convection coefficient at the inner walls can be estimated using a correlation presented by Incropera and DeWitt (1985) in the following functional form: 1/4 Nul = CjRaL . The Nusselt and Rayleigh numbers are based on the wall characteristic length. The value of the constant Cx at low Rayleigh numbers is specified by the wall position and the direction of the heat flow as ' 0.59 vertical walls Cl = 4 0.54 upper wall heated or lower wall cool ed _ 0.27 lower wall heated or upper wall cooled. 30 For the estimation of the forced convection coefficient at the outer walls, a widely employed correlation is expressed by Incropera and DeWitt (1985) in functional form as b 1/3 Nul = C2ReLPr where ReL is the Reynolds number based on the characteristic length of the wall, and the values of constants C2 and b used depend on the velocity and direction of the flow passing by the box walls (Table 3.1). Table 3.1 Forced convection constants Configuration C2 b Parallel Flow (laminar) 0.664 0.5 Parallel Flow (turbulent) 0.037 0.8 Flow normal to wall 0.228 0.731 Flow normal to wall (stagnation) 0.71 0.5 A correlation between heat of respiration and temperature for mature green tomatoes in the range of 5 and 26 °C is presented by Buffington and Sastry (1983): Qr = -20.35 + 18.11 T where Qr is in J/KG-hr and T in °C. A correlation between surface area and weight for mature green tomatoes is presented by Bellagha (1985): As = 5. 157x10’2(W)°'6875 where As is in m2 and W in kilograms. Inputs to the numerical model are listed in Appendix C. 31 Numerical Algorithms Energy conservation equations (3.5), (3.9), (3.10) and (3.11) generate the System of Linear Equations (SLE) to be solved for the three-dimensional temperature response of the fruit bin. Applying the transformation previously introduced in equation (3.1), m - h + 4L3[L2(i -l)+(j-l)] . This SLE can be expressed as M n+1 „?ApTp -° M n+1 p=i ^n.P^P = M n+1 P?1 VpTp = D 'M (3.19) where A,^ p = SLE matrix coefficients array, Dm= SLE right-hand side vector array, and M= total number of equations in the SLE ( M = 4L1L2L3) . The solution of the SLE is undertaken by using the set of algorithms shown in a general flow diagram in Figure 3.5 in a process that can be broken into two major steps. 1. The A^p coefficients and the Dm vector arrays are generated to form the SLE. This part of the process is accomplished by subroutines AIRNOD, SUFNOD, MIDNOD and CENOD, which 32 r START the bin temperature response. 33 generate the arrays for air phase and fruit surface, interior and center nodes, respectively. 2. The numerical algorithms to solve the SLE are implemented. Subroutine HERAT is invoked to handle this step of the process . Coefficient and Vector Arrays Generation As shown in Figure 3.5, bin physical properties and grid dimensions are transferred to the arrays generating subroutines from the Data and Grid Generation Block of the main program once during the overall process. On the other hand, the Codes Generation Block contains a bank of equation operators, rja, r;f, wf, and hwf, as well as the loop control indexes, MI, II, IF, JI, JF, HI, HF and WG. The use of these indexes and operators synthesizes the numerical computation of the A^ p and Dm arrays of the twenty-seven different types of air phase and fruit interior nodes (Table 3.2) into two compact algorithms contained in subroutines AIRNOD and MIDNOD, respectively. Table 3.2 Types of fruit and air phase nodes according to position in the bin. Group Number of subgroups I. Bin interior nodes (no wall contacts) 1 II. Nodes in contact with one wall 6 III. Nodes in contact with three walls 12 IV. Nodes in contact with four walls 8 Total: 27 34 The computation of the arrays for fruit center and surface nodes is done by subroutines CENOD and SUFNOD respectively. This computation is simplified by the fact that these nodes are not in contact with wall boundaries. In order to generate the time-dependent p and Dm arrays, the four previously mentioned subroutines are invoked every time step. Figure 3.6 presents a flow diagram of a typical array generating subroutine. As shown in this diagram, in every execution of the innermost loop equations, all non-zero A^ matrix coefficients in the m-row are calculated, including the Am m coefficient of the diagonal element as well as the coefficients of the nodes thermically connected to it. A set of r]a, nf, and wf operators is transferred to these loop equations for each of the twenty-seven types of nodes in the bin. This strategy provides flexibility for the computation of A^ and Dm arrays corresponding to elements in different positions throughout the bin using a common set of equations. Similarly, for each node subgroup, a set of indexes MI, II, IF, JI, JF, HI, HF and WG is sent to the three nested loops structure of the arrays subroutine. As shown in Figure 3.6, each set of indexes controls the execution of the Am p and Dm inmost loop equations over the three- dimensional range of positions of a given node subgroup. Numerical Solving Algorithms The numerical algorithms employed to solve the set of implicit energy equations governing the three-dimensional temperature response in the fruit bin are presented in the flow diagram of Figure 3.7. The solving procedure consists of a variation of the Gauss-Seidel Indirect 35 r L Figure 3.6. Typical flow diagram of an array generating subroutine. 36 Figure 3.7 37 Method which permits an independent treatment of each of the four nodes contained in a cubic volume element. For the Gauss-Seidel Iterative Method, convergence is guaranteed if the coefficient matrix \ p is diagonally dominant. Mathematically, this criterion is expressed as "f IVpl ^ I An, ml p^m ' (3.20) Nevertheless, in many cases, convergence is obtained with much weaker diagonal dominance than dictated by equation (3.14) (Hornbeck, 1975). Whenever convergence conditions are met, iterative techniques such as the Gauss-Seidel Method are widely preferred over direct solving procedures like the Gauss-Jordan Elimination for solving very large sets of equations. There are two reasons behind this preference: 1. Iterative techniques readily allow the storage of only those non-zero elements of the very sparse coefficient matrices which invariably result from very large sets of equations. Conversely, direct techniques cannot take advantage of this sparseness, and thus save storage space unless the coefficient matrix is banded, or in some cases where special programming techniques can be used. 2. With iterative techniques, roundoff errors are accumulated only in the final iteration. However, with direct techniques, roundoff errors incurred in a given step accumulate in subsequent matrix operations until the final solution is obtained. CHAPTER 4 EXPERIMENTAL SETUP AND TEMPERATURE MEASUREMENT PROCEDURES Description of Cooling Facility The environmental chamber shown in Figure 4.1 was used to conduct the various tests needed for the validation of the numerical model and the study of the effects of packing density, contact area and room air speed on the cooling of fruit. The chamber, a modification of a Bally Cooler, has inside dimensions of 1.4 x 1.8 x 2.5 m high. The chamber walls are constructed of 7.5 cm of foam insulation between aluminum sheeting panels bolted together. Figure 4.1 indicates the location of the main components of the refrigeration and air circulation systems. An induced-air Dayton model 1C791 blower runs continuously, circulating air through the chamber and the Filacell. A round section damper installed at the blower intake duct is used to control the air flowrate through the chamber. Cold water was sprayed onto a Filacell unit by a set of nozzles located at the top of the air supply duct. The purpose of the Filacell unit is to increase the rate of heat transfer between the spray water and the air. The Filacell is constructed of polypropylene filaments wrapped around a wooden frame in a position normal to the air flow. The Filacell unit is 0.61 m wide in its cross square section, and 0.91 m 38 39 COLD WATER TANK 40 deep in the direction of the air flow. The cold, wetted area formed by the adherence of water droplets to the filaments enhanced the heat transfer rate from the air stream to the water. The air leaving the Filacell unit was completely saturated. In the upper cold water loop, a pump operated continuously, feeding the nozzles with water suctioned from a 0.14 m3 mixing tank. In the lower loop, cold water was recirculated continuously by a second pump between the mixing tank and a 1 m3 storage tank installed underground. The storage tank was maintained at approximately 2.5 °C by a capillary tube connected to the compressor unit. A YSI 600 temperature probe installed at the inlet of the nozzles header controlled the temperature of the water sprayed on the Filacell unit. The signal from this probe controlled a three-way solenoid valve that regulated the flow of water from the cold storage to the mixing tank. The temperature of the water in the spray chamber was regulated by controlling the amount of cold water added to the mixing tank. This temperature was set to be equal to the desired dew point temperature inside the environmental chamber. As indicated in Figure 4.1, two additional YSI 600 probes were placed inside the chamber in order to accurately control the temperature of the air passing by the experimental container. The signals from the probes were sent to a pair of YSI model 72 proportional automatic controllers, with each controller energizing two 1000 watts cone heaters. 41 By using this system of probes and automatic temperature controllers, the room temperature at the container location was maintained within ± 0.1 °C of the set value of 8.1 °C dry-bulb and 7.2 °C dew-point temperatures selected for the tests. Figure 4.1 displays the experimental apparatus used to run room cooling tests in its position inside the environmental chamber. Room air was suctioned by a Dayton model 47525A variable volume blower through a 3.14 cm diameter PVC pipe and discharged through a wooden duct past the experimental container. A Dayton SCR controller was used to adjust the speed of the direct current motor driving the blower. The cross section area for airflow between the container and the wooden structure walls was 0.3 m2. A perforated plastic plate placed at the air inlet to the test section served to make the airflow around the container walls more uniform. An Omega HH-30 digital anemometer was used to measure air velocities at different locations in the duct. Data Acquisition System Thermocouples constructed from 36-gage, insulated copper- constantan wire were used to measure temperatures inside the experimental container. Temperatures were recorded on a Campbell Scientific 21X Micrologger microprocessor-controlled data acquisition system with installed capacity for 45 thermocouple signal inputs. The datalogger was programmed for scanning intervals varying from one to sixty minutes. Shorter intervals were programmed for the start of each run when stronger temperature variations were expected to occur 42 throughout the bin. A Kaye Icepoint apparatus was used as reference junction for the thermocouples connected to the datalogger. Data recorded during tests was sent from the datalogger to the microcomputer via modem for later processing. Reheat Facility The produce in the experimental container was brought to a selected uniform temperature initial condition before the beginning of each cooling test. The reheat unit used to carry out this preheating operation is shown in Figure 4.2. The unit, a detachahlp r nmnnnoivf n f a nvn/'nnl Inn iaaa n n .*j 14 flexible duct 43 Figure 4.2. Reheat facility. 44 Table 4.1: Physical and thermal properties of Sunny tomatoes Density 961.6 kg/nr 94.65 % 1.37 x 10~7 m2/s 4.01 KJ/KG°C 0.53 W/m °C Moisture content Thermal Diffusivity Specific Heat Thermal conductivity The experimental box was constructed of corrugated fiberboard with dimensions of 40 cm x 27.3 cm x 26.7 cm high, following the geometry of a 13.6 kg (30 lb.) tomato container (Turczyn and Anthony, 1986). For the straight column arrangement, a total of 96 fruit were placed inside the container, forming four horizontal layers (z-direction) . As illustrated by Figure 4.3, this arrangement yielded six and four vertical layers in the y- and x-directions, respectively. On the other hand, for the different tests run with random packs arrangement, the container accommodated an average of 106 fruit. Thermocouples were placed at the center of 16 tomatoes and near the surface of 8 tomatoes in key positions in the bin as shown in Figure 4.3. Ten thermocouples were also employed to measure temperatures in air voids and near walls inside the container. Two thermocouples monitored the temperature of the air flowing around the container. The fruit were weighed and placed into the experimental container. Table 4.2 presents results of experimental measurements of weight and average diameters of Sunny tomatoes. The product surface area exposed to natural convection was estimated by subtracting the measured contact area between tomatoes in a given packing arrangement from the product total surface area. 45 Figure 4.3. Thermocouple locations in the tomato container. Experimental measurements of contact areas in straight column arrangements (NCP=6) of tomatoes are presented in Table 4.3. 46 Table 4.2: Measured weight and diameter of Sunny tomatoes Sample # Weight (g) Average Diameter(cm) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 145.05 153.05 152.21 140.5 155.32 142.62 158.68 142.83 149.40 151.92 154.30 159.33 144.57 146.90 155.22 6.5 6.8 6.7 6.5 6.8 6.5 6.8 6.5 6.6 6.6 6.6 6.8 6.5 6.5 6.7 Average: 150.12 6.6 After the product was brought to a uniform predetermined temperature in the reheat unit, the container was hermetically sealed and placed in the experimental apparatus inside the environmental chamber. Cooling experiments were conducted with two packing arrangements and at three room air velocities. The two arrangements tested were the straight column pattern and random packing. The lowest room air velocity achieved, 18 m/min, corresponded to nearly still air conditions, with only the induced-air fan running. The two other air speeds used, 120 m/min and 240 m/min, are in the common range of velocities used in packinghouses. Table 4.3: Measured contact area 47 Sample # Surface 1 143 2 138 3 152 4 149 5 142 6 134 7 136 8 141 9 146 10 130 11 146 12 135 13 145 14 147 15 133 16 136 17 139 (cm2) Contact area 15.6 12.8 14.1 15.6 13.8 14.2 14.9 14.4 15.0 13.5 15.0 14.9 13.6 15.0 12.7 13.2 14.8 (cm2) % Contact area 10.9 9.3 9.3 10.5 9.7 10.6 11.0 10.2 10.3 10.4 10.3 11.0 9.4 10.2 9.5 9.7 10.6 Average: 14.3 10.2 Fourteen experiments were conducted to measure the temperature response at different positions inside the container. The scan interval for temperature measurements was set at 20 minutes for the first three hours of the process and at one hour for the remainder of the test. Tests were stopped when the temperature of the warmest fruit in the bin reached the target storage condition of about 13.5 °C. The operation of the proportional controllers maintained the room temperature within ± 0.1 °C of the target condition of 8.1 °C during the different tests conducted. CHAPTER 5 RESULTS AND DISCUSSION In order to validate the simulation model, the numerical results were compared against the experimental measurements of cubic packs of tomatoes inside a sealed fiberboard container. The containers were placed inside a cold room and subjected to three air velocities past their walls: 18 m/min, 120 m/min and 240 m/min. The tests and simulations at 18 m/min correspond to the case of containers cooled under nearly still air conditions. Additional comparisons were done against experimental measurements of random packs of tomatoes at the room air velocities of 120 m/min and 240 m/min. Numerical results for cooling processes were obtained by simulating the model in a VAX 8600 computer. The physical dimensions and thermal properties of the standard-size container of mature green tomatoes simulated by the model were presented in the previous chapter. The simulations were carried out for 26 hours, which was in all cases sufficient time for both the mass average temperature and the temperature of the warmest fruit in the bin to reach past standard storage conditions for tomatoes of about 13.5 °C (Kader et al . , 1985). The accuracy of the numerical model was examined by comparing test results with simulation results using time increments varying from 5 to 60 minutes. Improvements in accuracy averaging about 0.3 and 0.2 °C 48 49 were obtained by decreasing time increments from 60 to 20 minutes and 60 to 5 minutes, respectively. On the average, simulations used 24 seconds of computer time per hour of cooling to solve the system of 384 simultaneous equations of the experimental container of fruit, using a time step of 20 minutes. Very good agreement between experimental and simulation results at different positions in the container was obtained for the various room cooling stream velocities and packing patterns employed. The maximum deviation between measured and calculated values of temperature at any time during the different cooling processes studied was 0.5 °C. Based on its accuracy for predicting temperature responses, the model was used to simulate the cooling process of different packing patterns. Simulations were also conducted to predict the temperature responses and cooling times of fruit containers with their top and bottom walls not exposed to the room cooling air streams. Comparison of Experimental and Simulation Results Figures 5.2 through 5.10 present comparisons between experimental and simulation results for the temperature response of fruit in the six key positions of the container illustrated in Figure 5.1. As shown in these plots, there is very good agreement between experimental and simulation results for the different room air velocities and packing arrangements studied. The deviation of the simulated results from measured values did not exceed 0.5 °C during the different cooling processes studied. 50 Figure 5.1. Fruit locations in the experimental container. Temperature, 51 Figure 5.2. Comparison between experimental and simulation results for fruit 15 (bottom layer) and fruit 9 (third layer). Temperature, 52 Time, min Figure 5.3. Comparison between experimental and simulation results for the temperature response along a diagonal axis during room cooling of a bin of tomatoes at V= 18 m/min. Temperature, 53 Figure 5.4. Comparison between experimental and simulation results for the temperature response along a diagonal axis during room cooling of a bin of tomatoes at V=120 m/min. 54 Time, min Figure 5.5. Comparison between experimental and simulation results for the temperature response along a diagonal axis during room cooling of a bin of tomatoes at V=240 m/min. Tempe nature. 55 Time, min Figure 5.6. Comparison of experimental and simulation results for the temperature response along the central vertical axis during room cooling of a bin of tomatoes at V=18 m/min. Temperature 56 Time, min Figure 5.7. Comparison between experimental and simulation results for the temperature response along the central vertical axis during room cooling of a bin of tomatoes at V=120 m/min. Temperature, 57 Time, min Figure 5.8. Comparison between experimental and simulation results for the temperature response along the central vertical axis during room cooling of a bin of tomatoes at V=240 m/min. 58 Time, min Figure 5.9. Comparison between experimental and simulation results for the temperature response of tomatoes packed in random and straight column arrangements at V=120 m/min. 59 Time, min Figure 5.10. Comparison between experimental and simulation results for the temperature responses of tomatoes packed in random and straight column arrangements at V=240 m/min. 60 Figure 5.2 displays measured and calculated results of the behavior of temperature gradients inside fruit exhibiting high and low cooling rates. The closeness of the simulation results to the experimental values superimposed on this plot, indicate the validity of the model for predicting temperature gradients inside the fruit phase. Deviations between the measured and simulation results at no point during the cooling process exceeded 0.2 °C. Figures 5.3 through 5.8 allow the comparison of measured and simulation results for fruit located approximately along the diagonal and the central -vertical axes of the container. These plots reveal the different effects that diffusion and natural convection mechanisms have on the cooling rates of fruit according to their locations along these axes. At the same time, these plots illustrate the evolution through time of temperature differences among fruit inside the container. The comparison of simulation and experimental results at the three room air velocities employed verify the accuracy of the model for predicting the three-dimensional temperature distribution inside the container. The level of accuracy is consistent at the three different velocities employed. The maximum difference between measured and experimental results in the different cooling processes studied was 0.5 °C. Figures 5.9 and 5.10 display results for the tests and simulations done with random packs of fruit contrasted with those of straight column packs under the same room conditions. As expected from the difference in packing densities between these two arrangements, there is a resulting difference in their temperature responses. As observed in 61 Figures 5.9 and 5.10, the model is able to accurately predict the temperature responses of the fastest and slowest cooling fruit in the container for both packing arrangements . For the random packing arrangement, the difference between experimental and numerical results did not exceed 016 °C during the various cooling processes conducted. Based on the very good agreement between simulation and experimental results exhibited in Figures 5.2 through 5.10, the model was used in the simulation of cooling processes of random and pattern packs at different air velocities. Discussion of Results for the Cubic Packing Pattern The results in Figures 5.2 through 5.8, and Figures 5.11 through 5.12 can be used to analyze the temperature responses of fruit in the six key positions of the container illustrated in Figure 5.1. The simulation and measured values plotted on these curves correspond to the cubic packing arrangement at the three previously specified cooling stream velocities. The behavior of temperature gradients within fruit during the cooling process is illustrated in Figure 5.2 by using sets of curves representing two extreme cooling rates. The upper set of curves corresponds to results for the slow cooling rate fruit in a bin cooled at low room air stream velocity. On the other hand, the lower set was obtained from a fast cooling rate fruit in the bin cooled at high room air stream velocity. In both cases, as shown, temperature gradients between fruit surface and center are very small with differences of no more than 0.3 °C throughout the cooling process. 62 SLOT COOLING FRUIT o Center temperature of fruit 9 □Surface temperature of fruit 9 ■Air phase tancerature fflST COdLING FRUIT ^Center temperature of fruit !8 □ Surface temperature of fruit 19 iflir phase tenperaturs ■ Hal I taprature u:u o Room ai r speed; 13 n/nin u-$:o,A Hoon isiperature: 8.! C Bl,U&Ovv . D'D-DfiAA ■m. O-D-Ma., 'D:n:6:ct ■■■■ i 0 OO 120 180 240 300 360 420 480 540 600 500 720 780 840 Time, min Figure 5.11. Comparison of simulation results of temperature gradients development inside and at the surface of tomatoes during room cooling at V=18 m/min. Temperature, 63 26 25 24 oo L.-U 22 21 I'A / / 16 14 13 12 1 1 10 3 3 V I \ ■B;A ' \x Ro □v> □:a □ a □.vs Oa □:n • □:a. ■ uo.. . ■ ■ uA. ■ u,0 ■._ uAa ■ _u.0.A I . u-2a J U.O, » 1i.u'lj;o., ftxu air speeds nvnin B'b-ii^;uO ^a5 ^ra^rs: S.1 C HjjW B1-'^a " i i-V- A □ "A u-n-A U.g '“Sit. :«:24 l 0: □:n-A A aO-Q:A, ■■■■■ :m-n: — J L 0 SO 120 180 240 300 3G0 420 480 540 600 660 720 730 640 Time, min Figure 5.12. Comparison of simulation results of temperature gradients development inside and at the surface of tomatoes during room cooling at V=240 m/min. Temperature 64 24 22 2! 20 19 'ta 171- 16 *- 15 r \ 14 13 12 10 \ % % I V % \V\ % ■V ^Straight coh.im pattern ^Roidan packing ■Straight col urn pattern with 16 l increase in contact area' Room air sceed : 120 wm Room temperature: 8. ! G A Q4 storage tenperoture 4U im 'IC'AO--'. 3 SO 120 180 240 300 360 420 480 540 600 500 720 780 840 Tim, min. Figure 5.13. Simulation of the effect of increasing fruit contact area on cooling times of tomatoes at V=120 m/min. 65 The combination of small temperature gradients inside the fruit and the low natural convection coefficients at the surface of the fruit translates into low Biot numbers for the fruit phase in the container. This last fact depicts a distinctive contrast between room cooling process and forced cooling process which is characterized by high Biot numbers. Furthermore, these experimental results confirm the necessity of using only three fruit nodes in this model. Figures 5.11 and 5.12 display computer simulation results of temperature gradients both inside the fruit and at the surface of elements with slow and fast cooling rates. The upper set of curves in each of these two figures corresponds to a slow cooling rate fruit and its surrounding air phase element located towards the core of the fruit bin. As illustrated in these plots, the slow heat removal rates is due to the fact that temperature gradients both inside the fruit and at its surface remain small throughout the cooling process at these interior locations. Due to the small rates of heat removal by natural convection, the contribution of heat diffusion across fruit contact areas is more significant in locations closer to the core of the bin than near the walls. For elements with fast cooling rates, there is a stronger contribution of natural convection in the heat removal process resulting from the development of substantial temperature differences between fruit surface and surrounding air soon after the start of the cooling process. The lower set of curves in Figures 5.11 and 5.12 depicts this situation using results for a fruit element located at the corner of the 66 bottom wall. The contribution of conduction heat transport is also enhanced because of the contact areas with the cold walls. In order to examine the evolution of temperature distributions in the bin throughout the cooling process, measured and simulation results for the six key locations previously shown in Figure 5.1 are displayed in Figures 5.3 through 5.9. Figures 5.3, 5.4 and 5.5 present the results for fruit located approximately along the diagonal axis of the container. The effect of natural convection velocities on the temperature distribution along the layers in the container becomes evident by comparing the temperature differences existing between the curves for fruit 27 and 19. These fruit are located in the corners of the top and bottom walls respectively (Fig. 5.1). In a purely conductive heat removal process, the temperature responses of these two fruit would be the same. In a conductive- convective phenomenon, however, as the cooling process advances, temperature gradients, and consequently natural convection velocities, grow throughout the container. This brings about an increase in the contribution of natural convection buoyancy to heat transport from the lower towards the upper layers of the container. Figure 5.4 reveals measured temperature differences of as much as 2 °C five hours after the start of the test between these two fruit spaced three layers apart. The influence of natural convection velocities on the temperature distribution along layers is also observed in the results for fruit numbers 21 and 10 plotted in Figures 5.6, 5.7 and 5.8. Fruit 19 and 21 are located approximately along the central vertical axis of the 67 container in the second and top layers respectively (Fig. 5.1). As illustrated in these plots, fruit 21, in spite of not being in contact with any of the cold walls, cools faster than fruit 10 because of the transport of heat upwards by convective buoyancy. An interesting result of the influence of natural convection velocities previously discussed is the displacement of the location of maximum temperature from the center of the bin upwards as the cooling process advances. This upwards displacement takes place faster as the room cooling air velocity increases. The model predicts that in a standard box of tomatoes which accommodates four layers, the warmest fruit throughout the cooling process is located one layer above the center line of the container. The measured results obtained verified this fact as fruit 9 located in this position of the bin reported the warmest temperatures throughout the different tests conducted. The results in figures 5.3, 5.4 and 5.5 can also be used to assess the different impacts that variations in room air velocity have on the cooling time of fruit located near the container walls as compared to that of fruit located towards the center of the bin. An illustration of this is given by the increase in the difference between the cooling times of fruit 9 and 19, from approximately eight hours when the container is cooled in nearly still air conditions, to about ten hours when a room air speed of 120 m/min is used. Discussion of Results for Random Packs Figures 5.9 and 5.10 display the results for the experiments and simulations done with random packs of fruit contrasted with those of straight column packs with the same room conditions. In order to 68 simulate the model for random packing, the average contact area and packing density for the straight column pattern in the base program was substituted by the corresponding values for random packs. In both packing arrangements, fruit 9 and 19 are located in approximately the same positions in the container. As illustrated in Figures 5.9 and 5.10 by experimental and simulation results, fruit arranged in random packs exhibit faster cooling rates than those packed in straight column pattern. Higher rates of heat of respiration per unit volume of container are expected with the random packs arrangement because of its higher packing density. In random packs, however, there is larger contact area and consequently higher rates of heat conduction between fruit, because of the higher number of contact points, seven as compared to six for straight packs. This increase in heat removal by conduction through the fruit phase becomes the controlling effect, resulting in the faster cooling rates of random packs illustrated by the results in Figures 5.9 and 5.10. The competing effects of packing density and contact area on cooling rates can be further examined through the simulation results of mass temperature responses plotted in Figure 5.13. The fastest cooling rate curve corresponds to the simulation of a 16% contact area increase in a straight column arrangement without increasing its packing density. As shown, this brings about a 70 minutes (12.9%) reduction in cooling time in comparison to the straight column arrangement with regular contact area. In contrast, the plot for random packs reveals a reduction in the order of 30 minutes (5.6%) only, despite the fact that this arrangement also has an average contact area roughly 16% larger 69 than the straight column packs. This smaller reduction in cooling times is due to the higher packing density and the consequent higher volumetric heat generation rate in random packs. As stated above, heat conduction through contact areas is expected to be more significant for fruit located towards the center of the container than for those near walls. This fact is also demonstrated in these same figures by the stronger effect that the variation of packing arrangements has on interior fruit 9 than on fruit 19 located near the wal 1 . Effects of Room Air Velocity and Packing Arrangement The influence of room air velocity and packing arrangement on cooling times of fruit inside containers can be studied utilizing the results presented in Figures 5.14 though 5.22. Measured and calculated results for fast cooling fruit 19 at the three room air speeds used in the tests are displayed in Figure 5.14. A reduction of 2.3 hours (35%) in the time needed to reach storage temperature is obtained when the room air speed is increased from nearly still conditions (18 m/min) to 120 m/min. Doubling the room air speed from 120 m/min to 240 m/min further reduces the cooling time by one hour (21%). Smaller reductions of about 17% and 6%, respectively, are obtained for slow cooling fruit 21 with the same speed increments (Figure 5.15) since it is located further away from the walls. Based on the accuracy of the model to predict heat removal rates from fruit in containers, simulations were carried out to study the effects of room air velocity and packing arrangement on temperature response and cooling time. Temperature, 70 Figure 5.14. Effect of room air velocity on the cooling time of a fast cooling fruit. Tempe nature 71 T I ire, min Figure 5.15. Effect of room air velocity on the cooling time of a slow cooling fruit. 72 Time, min. Figure 5.16. Simulation of the effect of room air velocity on the mass average temperature response of random packs of tomatoes. Temperature. 73 24 on 22 21 20 a 19 *13 17 18 15 14 13 !2^ 1 1 r N yq \ .a, ■\\ \ \0. \'A\ \ V \ \ \D, \\ \ \ A \ Vv ■ \ " straight coium pattern 18 m/mSn roan air speed Aflt 120 m/min rcan air spged ■fit 240 i/'min roan air soeed Room ienperature: 8. \ \ X, \ A 0, x XX N[1 NX X I X "'n X 'A '-i X Storage tenperature 0.7 hrs >1,3 hrs 0 60 120 180 240 300 350 420 480 540 SCO 560 720 780 Tims, min. 840 Figure 5.17. Simulation of the effect of room air velocity on the mass average temperature response of straight column packs of tomatoes. Temperature 74 a 24& 33 22 21 20 19 15 17 16 15 14 13 12 10 \\ 3d '".tx \\°x \'A\ B,\A \ \ \ \'a ■ BCC packing pattern □fit 18 m/min room air speed Afit 120 m/min roam mr spaBd ■fit 240 m/min roan air speea _L I \ ' \\ Room taiuerature1 8-1 C A \ « ^ O.Shrs , 11 hrs \K yC / Storage temperature B\ a^/ '"'□I x.« "A ‘D^ -A. _L X '1 0 SO 120 180 240 300 360 420 480 540 5C0 660 720 780 Tim©, min. Figure 5.18. Simulation of the effect of room air velocity on the mass average temperature response of tomatoes packed in a body-centered cubic (bcc) pattern. Temperature 75 0 120 240 060 480 600 720 840 960 1080 1200 1320 1440 1550 Time, min Figure 5.19. Comparison of mass average temperature with fastest and slowest cooling fruit temperature responses in a random pack of tomatoes. Temperature,, 76 :4 \ 23 22 21 K 20 !9 - 1 !8 - !7f ia - 15 14 13 12 II 10 9 °fi T 702 of measured contact area o Straight colum pattern Random Packing ■ p pecking pattern • 52 vented container (test) \X«^ \ \ Roan air soeoa: !8 m/min Room tenpe’rature: 8.1 C / ' '•■O' r* 2.2 hrs 0 so 120 180 240 300 360 420 480 540 600 650 720 730 840 Time, min. Figure 5.20. Simulation of mass average temperature responses of different packing arrangements of tomatoes at V= 18 m/min. Temperature 77 Time, min. Figure 5.21. Simulation of mass average temperature responses of different packing arrangements of tomatoes at V= 120 m/min. Temperature 78 Time, min. Figure 5.22. Simulation of mass average temperature responses of different packing arrangements of tomatoes at V=240 m/min. 79 Figures 5.16, 5.17 and 5.18 present responses of the mass average temperature fruit in the container to different packing arrangements under varying room cooling air velocities. These figures can be used to make a general evaluation of effects of room air speed variations on different packing arrangements. Simulation results for random packs of fruit plotted in Figure 5.16 show cooling times of 10.5 hours, 8.7 hours and 7.9 hours at room air speeds of 18, 120 and 240 m/min respectively. This leads to corresponding reductions in cooling time of about 1.8 hours and 0.8 hour for the two air speed increases. A comparison of Figures 5.16, 5.17 and 5.18 shows that these cooling times and the corresponding reductions for random packs are found between the corresponding values for straight column packing and bcc pattern. It is also observed that in the range of air speeds studied, the random packing results are closer to the straight colun pattern. These results confirm the suitability of using the straight column pattern to build the base model to predict cooling rates for random packs. In packing houses, one of the most important aspects in the cooling process of produce is whether to determine cooling sufficiency by the mass average temperature or the warmest temperature in the batch. In the case of room cooling, temperature gradients generated by the end of the process are expected to continue during storage. Thus, in order to avoid deterioration, in room cooling operations it is advisable to monitor the warmest temperature in the batch (Mitchell et al . , 1985). 80 Figure 5.19 illustrates this situation using a contrast of the mass average temperature simulation results and the measured temperature responses of the fastest and slowest cooling fruit in a random pack. As shown, if the mass average temperature were used as an indicator of sufficient cooling, there would be a temperature difference of about 2.7 °C with respect to the warmest fruit going into the storage room. Figure 5.19 also provides an estimation of 24 hours for the seven- eighths cooling time of a tomato container of dimensions 40 cm x 27.3 cm x 26.7 cm high, cooled at 120 m/min room air velocity. This is by definition the time that it takes to cool produce to one-eighth of the initial produce-coolant temperature difference. In industrial practice it is known that well vented containers cool much faster than sealed ones (Kader et al . , 1985). This fact is illustrated in Figures 5.20, 5.21 and 5.22 which display mass average temperature curves of simulation results for different packing arrangements. The bottom curve in each figure corresponds to results for a random pack container with five percent vented sides. At room air speed of 120 m/min for instance, Figure 5.21 reveals a cooling time reduction in the order of 4.7 hours (54%) for a container with five percent venting area. Figure 5.22 shows that this reduction in cooling time grows rapidly as the room air speed increases. Effects of Uncooled Walls A very common practice in packing house operation consists of stacking fruit containers on top of each other on pallet boxes during the cooling process. Whenever this practice is followed, cooling times 81 are expected to rise considerably as a result of the isolation of top and bottom walls from the room cooling air stream. The numerical model was used to predict the retardation effect in cooling times of random packs caused by the isolation of top and bottom walls from the room cooling stream. Figures 5.23 and 5.24 present the results of these simulations along with numerical and test results for a fruit container with all six walls exposed to the room air streams. As illustrated in Figure 5.23, on a mass average basis, the model predicts an increase of approximately four hours in the cooling time of containers stacked on top of each other. Figure 5.23 also shows that for storage temperature criterions below 13.5 °C, the retardation in cooling times is even larger. In Figure 5.24, a similar evaluation is carried out using the temperature of the warmest fruit in the bin as the criterion for sufficient cooling. From the results plotted in this graph, it is evident that the impact of unexposed walls on retardation in cooling times is more severe for the slower cooling produce in the bin. As shown, the retardation in cooling time based on the slowest cooling fruit in the bin is approximately 4.8 hrs. The results presented in this chapter indicate the validity of the numerical model developed for predicting three-dimensional temperature responses in both pattern and random packs of produce. Both numerical and test results show an upwards displacement of the location of maximum temperature in the bin above the center line of the container due to convective buoyancy. Tempera lure, 82 24 I i 23 22 21 20 a IQ ' !8 h 17 16 L !5 14 13 12 I 0\ V 0 x, \ \ \ a., \ \ \ Random Packs ■ Covered top and bottan sides 0 Expend tip md bottom side3 Room air speed •* 120 m/min Room tenperature: 8,1 C v U \ Storage Temperature o. -5sK:-- 4 hrs i 0 60 120 180 240 300 360 420 480 540 500 560 720 750 840 Tims, min Figure 5.23. Simulation of mass average temperature responses of containers of tomatoes with exposed and covered top and bottom walls. Temperature 83 Figure 5.24. Simulation and test results of the fastest and slowest temperature responses in random packs of tomatoes inside containers with exposed and covered walls. 84 The results also reveal the effects of increasing room air velocity and packing density on achieving considerable reductions in cooling time. The differences in using the mass average temperature and the temperature of th'e warmest fruit as criterions for determining cooling sufficiency were examined. The elevated temperature gradients at the end of cooling render the latter criterion more convenient in room cool ing. Simulations carried out using the model indicate a considerable retardation in cooling time for containers stacked on top of each other. CHAPTER 6 CONCLUSIONS AND RECOMMENDATIONS The close agreement between simulation and experimental results obtained in this work indicates the validity of the numerical model developed for predicting three-dimensional temperature responses during room cooling of both random and pattern packs of fruit confined in containers. Good accuracy was obtained in the simulation of random packs when its average values of contact area and packing density substituted those of straight column pattern in the base program. Experimental and numerical results have confirmed that a pattern of very small temperature gradients exists inside fruit and that low convective heat transfer coefficients exist at the surface of each item of fruit. These have been shown to be the main factors behind the slow cooling rates that are characteristic of room cooling. Also, because these internal gradients were small, the numerical model used called for having only three nodes inside the fruit. Because of the small temperature gradients at the surface of fruit that are located in the core of the bin, the contribution of heat diffusion across fruit contact areas has been shown to be more significant at the core than near the walls. The influence of natural convection velocities was found to cause the displacement of maximum temperature from the center of the bin upwards as the cooling process advanced. 85 86 The experimental and numerical results indicate that by increasing room air velocities from 18 m/min to 120 and 240 m/min, the cooling times of tomato containers are reduced by 17.4% and 25.2%, respectively. As a result of larger contact area, random packs exhibit faster cooling rates in comparison to the straight column pattern. This trend was reflected in simulations conducted with other high-density packing patterns. Substantial temperature differences between fruit in different container locations remained throughout the cooling process. The results support determining cooling sufficiency based on the warmest temperature in the bin rather than on the mass average temperature. Experimental results obtained for tomato containers with 5% side venting areas at different room air speeds have indicated 50% reductions in cooling time when compared to that for sealed containers. The numerical model was used to predict the retardation effect on cooling times for tomato containers stacked on top of each other, a common practice in packinghouse operations. With conditions similar to those used during the experimentation, increases of up to 40% in cooling times may be predicted by the simulation results when the top and bottom walls are isolated from room air. APPENDIX A MAIN PROGRAM: THREE-DIMENSIONAL TRANSIENT TEMPERATURE RESPONSE SOLUTION IN A FRUIT BIN C NOMECLATURE C A: Matrix Coefficient C AC: Individual contact area C AAP: Air isobaric thermal C diffusivity C AAV: Air isometric thermal C diffusivity C APP: Fruit element thermal C diffusivity C RF: Fruit radius C AR: Radial distance C between fruit nodes C AS: Fruit surface area/cubic vol . C AN: Fruit nodes average C conduction area C AV: Fruit element volume C C: Matrix vector C DFP: Fruit element density C DT: Time increment KA: Air thermal conductivity KF: Fruit thermal conductivity KP: Convective velocity coefficient Li: Wall lenght NA: Air element derivative operator NF: Fruit surface node derivative operator RE: Fruit element outside radius RFA: Fruit average radius RN: Element node radius T: Single index temperature representation TF: Fruit initial temperatures TK: Triple index-time temperature representation TW: Wall temperature 87 88 C DL: Fruit diameter C DLE: Wall equivalent lenght C FAC: Contact area fraction C HA: Fruit surface C convective coefficient C HIR: Inner wall convective C coefficient C HOR: Outter wall C convective coefficient C HW: Box wall conductivity C HWA: Box wall overall C conductance C HWF: Wall-fruit surface C contact conductance C JET: Air jet direction index TO: Air jet temperature V: Room air velocity VI: Convective velocity WA: Air nodes diagonal element operator WF: Fruit surface diagonal element operator WG: Fruit surface nodes loop operators WH: Air nodes loop operator COMMON A (384, 384) ,TK(0:840,4,6,20) ,C(384) ,KF,RF(384) ,DL(384) , VI(0:30,6,6,-5:17),HW(6),JET(6),DLE(6),V(6) ,T0,AS(384) , Q,DT, LI , L2, L3,HA(4,6, 17) , HOR (6, 4, 6, 17) , AN (384,0: 17,0:17), AV(384, 0 : 17), AC(384), AR(384, 0:16,0: 17), RN(384, 0:17), RE (384,0: 17) ,TW (0:600, 4,6, 17) ,HFR(6,4,6, 17) , VO (840) DIMENSION NA(6) , NF(3) ,HWF(6) ,HWA(6) ,DLW(6) ,T(1200) ,TIME( 60 ) , DFP(600) ,APP(600) INTEGER WG ( 5 ) , WH ( 5 ) , WA ,WF,WI ,HI,HF,H,Q,F,P REAL KF, KP,KW(6),NI (0:6), NODI, NODT CHARACTER SELE*2 LOGICAL PATH 89 PARAMETER (PI=3. 1416) 0PEN(UNIT=0,FILE='TE.DAT',STATUS='0LD') OPEN (UNIT=1 , FILE=7 TEST. DAT' , STATUS='OLD' ) OPEN(UNIT=2, FILE=' TEMP. OUT', STATUS= ' UNKNOWN' ) OPEN ( UN IT=7 , FILE='TSTORY.OUT' ,ACCESS=' APPEND' ,STATUS=/ UNKNOWN' ) 10 q=o INDEX=2 IT=1 T IME ( 0 ) =0 . 0 READ(0,*) DL1 ,DL2,DL3 READ(0,*) (DLE(I) , 1=1 ,6) READ(0,*)(HW(I),I=1,6) READ(0,*) (V(I) , 1=1,6) READ(0,*) (JET (I), 1=1,6) READ(0,*) AFP,KF,RFA,FAC READ(0,*) TO,TF,DT C THE FOLLOWING ALGORITHM ESTIMATES BOX AND FRUIT ELEMENTS C DIMENSIONS AND PROPERTIES NEEDED IN SUBSEQUENT MATRIX C COEFFICIENTS SUBROUTINES. C NUMBER OF ELEMENTS IN BOX LI DIRECTION L1=IFIX(DL1/(2*RFA) ) L2=IFIX(DL2/(2*RFA) ) L3=IFIX(DL3/(2*RFA) ) N0DT=L1*L2 C FRUIT ELEMENTS DIMENSIONS P=3 DO 11 KM=4,4*L1*L2*L3,4 RF ( KM) =1 .4 90 11 CONTINUE READ(0,*)RF (4) , RF(60) , RF(96) , RF( 108) , RF( 148) , RF( 160) , RF( 168) , + RF ( 236 ) , RF(384) , RF(84) , RF( 16) ,RF(372) , RF(8) ,RF(348) DO 22 1=1, LI DO 22 J=1 , L2 DO 22 IH=2,4*L3-2,4 M=( IH+2)+4*L3*(L2*( I-1)+(J-1) ) DL(M)=2*RF(M) RN(M, IH) =RF (M) RN(M, IH+P)=0.0 RE (M, IH)=RF(M) RF(M, I H+P) =0 . 0 DO 14 F=IH+1 , IH+P-1 RE(M, F)=( IH-F+P)*RF(M)/(P-1) RE(M,F+1)=(IH-(F+1)+P)*RF(M)/(P-1) IF(F.LT. IH+P) THEN RE(M, F+2)=( ( ( IH- (F+2)+P) )/ (P-1) )*RF(M) ELSE RE(M,F+2)=0.0 ENDIF RN(M,F)=(((RE(M, F) **3 )+( RE (M, F+l )**3) )/2)**. 33333 RN(M, F+l )=(( (RE (M, F+l )**3)+(RE(M, F+2)**3) )/2)**. 33333 AR(M,F,F-1)=RN(M,F-1)-RN(M,F) AR(M,F,F+1)=RN(M,F)-RN(M,F+1) 91 AN ( M , F , F- 1 )=4*PI*RN(M, F)*RN(M, F- 1 ) AN(M, F, F+l )=4*PI*RN(M, F)*RN(M, F+l ) AV(M, F)=(4*PI/3)*( (RE(M, F)**3)- (RE(M, F+l)**3) ) 14 CONTINUE AC(M)=FAC*(4*PI*RF(M)**2)/6 AS(M)=12*(PI*DL(M)**2-6*AC(M) )/(DL(M)**3) 22 CONTINUE DO 25 IM=3 , 4*L1*L2*L3 , 4 DFP( I M ) =60 . 0 DFP ( I M+ 1 ) =60 . 0 25 CONTINUE READ (0,*) DFP (3) ,DFP(4) ,DFP(59) , DFP(60) ,DFP( 107) ,DFP(108) , + DFP(147),DFP(148),DFP(159),DFP(160),DFP(235),DFP(236), + DFP(167),DFP(168),DFP(95),DFP(96),DFP(383) , DFP (384) , + DFP (371) ,DFP(372),DFP(15) , DFP (16) , DFP (347) , DFP (348) C THIS IS THE INITIALIZATION BLOCK FOR THE COEFFICENT MATRIX C AND TEMPERATURE DISTRIBUTION. DO 88 M=1 ,4*L1*L2*L3 T (M)=TF 88 CONTINUE READ (0,*)T( 235) ,T(167) ,T(95) ,T(236) ,T( 160) ,T( 148) ,T(292) , + T(372),T(4),T(60),T(168),T(84),T(108), T(284) , + T(348) ,T(96) ,T(304) ,T(16) ,T(384) IM=1 DO 95 1=1, LI DO 95 J=1 , L2 92 DO 95 H=1 ,4*L3 TK(0, I,J,H)=T(IM) IM-IM+1 95 CONTINUE 99 TIME( IT) =T IME ( IT- 1 ) +DT DO 100 M=1 , 4*L1*L2*L3 DO 100 N=1 ,4*L1*L2*L3 A(M, N)=0 . 0 C(M)=0.0 100 CONTINUE E= . 476 N0DI=0 . 0 SUMV=0.0 DO 105 1=2, ( LI - 1 ) DO 105 J=2, (L2-1) H=4*L3-7 M=H+4*L3*(L2*(I-1)+(J-1)) KP=( ( 1 • 16E-5)*E**3*(DL(M+3)**2) )/( ( 1-E)**2) VI(Q,I,J,H)=(1 .37E+6)*KP*ABS(TK(Q, I,J,H)-TK(Q,I,J, H+4) ) SUMV=SUMV+VI (Q, I , J,H) NODI=NODI+l 105 CONTINUE VO(Q)=(NODI/ (NODT-NODI) )*SUMV C IN THE FOLLOWING BLOCK SUBROUTINES SUFNOD AND CENOD ARE CALLED TO C GENERATE THE MATRIX COEFFICIENTS OF THE ENERGY EQUATION FOR THE 93 C SURFACE AND CENTER FRUIT NODES, RESPECTIVELY. CODES AND OPERATORS C FOR BOUNDARY NODES ARE SENT TO SUBROUTINES MIDNOD AND AIRNOD TO C GENERATE MATRIX COEFFICIENTS FOR FRUIT INTERIOR AND AIR NODES, C RESPECTIVELY. CALL CENOD(AFP,DFP,APP) REWIND 1 DO 110 I M= 1 ,27 READ ( 1 , *) WF READ ( 1 ,*) (HWF (IK), I K= 1,6) READ(1 ,*) (NF (IK) , IK=1 ,3) ,NA(1) ,NA(6) , (NA(IK) , IK=2,5) READ(1 ,*) ( WG (IK) , I K= 1,5) READ( 1 ,*) (NI ( IK) , IK=0, 5) C TYPE I-S IF(IM.EQ.l) THEN MI=4*L3*(L2+l)+6 11=2 I F=L 1 - 1 J 1=2 J F=L2 - 1 H I =6 HF=4*L3-6 C TYPE I I-S ELSEIF(IM.EQ.2) THEN MI=4*L3+6 11=1 94 I F= 1 J 1=2 JF=L2 - 1 H I =6 HF=4*L3-6 ELSEIF( IM.EQ.3) THEN MI=4*L2*L3+6 11=2 I F=L 1 - 1 JI=1 JF=1 H 1=6 HF=4*L3-6 ELSEIF(IM.EQ.4)THEN MI=4*L3*(2*L2- 1 )+6 11=2 IF-L1-1 JI=L2 JF=L2 H 1=6 HF=4*L3-6 ELSEIF(IM.EQ.5) THEN MI=4*L3*(L2+l)+2 11=2 I F=L 1 - 1 J 1=2 95 JF=L2- 1 H I =2 HF=2 ELSEIF( IM. EQ.6) THEN MI=4*L3*(L2+2) -2 11=2 IF-L1-1 J 1=2 JF=L2- 1 HI=4*L3-2 HF=4*L3-2 ELSEIF(IM.EQ.7) THEN MI=4*L3*(l+(Ll-l)*L2)+6 I I=L1 IF=L1 J 1=2 JF=L2- 1 H I =6 HF=4*L3-6 ELSEIF( IM.EQ.8) THEN M I =6 11=1 I F= 1 JI = 1 JF=1 HI=6 96 HF=4*L3-6 ELSEIF(IM.EQ.9)THEN MI=4*L3+2 11=1 I F= 1 ' J 1=2 JF=L2-1 H I =2 HF=2 ELSEIF( IM. EQ. 10) THEN MI=8*L3-2 11=1 I F= 1 J 1=2 JF=L2-1 HI=4*L3-2 HF=4*L3-2 ELSEIF( IM. EQ. 11) THEN MI=4*L3*(L2-l)+6 11=1 I F= 1 JI=L2 JF=L2 H 1=6 HF=4*L3-6 ELSEIF(IM.EQ.12) THEN 97 MI=4*L2*L3+2 11=2 I F=L 1 - 1 JI=1 JF=1 ' H 1=2 HF=2 ELSEIF(IM.EQ.13)THEN MI=4*L3*(L2+1 ) -2 11=2 I F=L 1 - 1 JI=1 JF=1 HI=4*L3-2 HF=4*L3-2 ELSEIF(IM. EQ. 14) THEN MI=4*L3*(2*L2-l)+2 11=2 I F=L 1 - 1 JI=L2 JF=L2 H 1=2 HF=2 ELSEIF(IM.EQ. 15) THEN MI=8*L2*L3-2 11=2 98 I F=L 1 - 1 JI=L2 JF=L2 HI=4*L3-2 HF=4*L3-2 ELSEIF( IM. EQ. 16) THEN MI=4*L3*L2*(Ll-l)+6 I I=L1 IF-L1 JI=1 JF=1 H I =6 HF=4*L3-6 ELSEIF(IM.EQ.17)THEN MI=4*L3*L2*(L1 -1 )+4*L3+2 II=L1 I F=L 1 J 1=2 JF=L2-1 H I =2 HF=2 ELSEIF(IM.EQ. 18) THEN MI=4*L3*L2*(Ll-l)+8*L3-2 I I=L 1 IF=L1 J 1=2 99 JF=L2-1 HI=4*L3-2 HF=4*L3-2 ELSEIF( IM. EQ. 19)THEN MI=4*L3*(L2*Ll-l)+6 II=L1 IF=L1 JI=L2 JF=L2 H I =6 HF=4*L3-6 ELSE I F ( IM. EQ. 20)THEN MI=2 11 = 1 I F= 1 JI=1 JF=1 HI-2 HF=2 ELSEIF(IM. EQ.21) THEN MI=4*L3-2 11=1 I F= 1 JI=1 JF=1 HI=4*L3-2 100 HF=4*L3-2 ELSEIF(IM.EQ.22) THEN MI=4*L3*(L2-l)+2 11=1 I F= 1 JI=L2 JF= L2 HI-2 HF=2 ELSEIF(IM.EQ.23)THEN MI=4*L2*L3-2 11=1 I F= 1 JI=L2 JF=L2 HI=4*L3-2 HF=4*L3-2 ELSEIF(IM.EQ.24) THEN MI=4*(Ll-l)*L2*L3+2 II-L1 IF=L1 JI=1 JF=1 H 1=2 HF=2 ELSEIF( IM. EQ.25) THEN MI=4*(Ll-l)*L2*L3+4*L3-2 1 1 = L 1 I F=L 1 JI = 1 JF=1 ' HI=4*L3-2 HF=4*L3-2 ELSEIF ( IM. EQ . 26) THEN MI=4*L3*(L2*Ll-l)+2 I I = L1 IF-L1 JI=L2 JF=L2 H 1=2 HF=2 ELSEIF( IM. EQ.27) THEN MI=4*L3*L2*L1 -2 I I=L1 IF=L1 JI=L2 JF=L2 HI=4*L3-2 HF=4*L3-2 ENDIF 101 CALL AIRNOD(MI-l,II,IF,JI,JF,HI-l,HF-l,WF,HWF,WG,NF,NA,NI) 102 CALL MIDNOD(AFP,DFP,APP,MI+l, 1 1 , IF , J I , JF, HI+1 , HF+1 , + WF,HWF,WG,NF,NA) CALL SUFNOD(AFP, DFP, APP) 110 CONTINUE REWIND (1) MT=4*L1*L2*L3 NT=MT CALL MATRIX(*205, MT, NT, INDEX, T) Q=Q+DT M=1 DO 160 1=1, LI DO 160 J=1 , L2 DO 160 H=1 ,4*13 TK(Q, I,J,H)=T(M) M=M+1 160 CONTINUE IT=IT+1 DT=60. IF( IT. LT. 16) THEN GO TO 99 ELSE GO TO 206 ENDIF 205 WRITE(7,*) 'MATRIX HAS NO SOLUTION' 206 WRITE(7,210) (TIME(IC) , IC=0, 15) 210 FORMAT(25X,' FRUIT NODES TEMPERATURE HISTORY' ,/,9X, 'TIME(MIN) ' 103 + ,1X,17(F4.0,2X)) WRITE (7,211) 211 FORMAT (4X, 'COORDINATES' ,/,3X, ' I ' ,3X, ' J' ,3X, 'H' ,3X, 'M' ) DO 220 1=1, LI DO 220 J=1 , L2 DO 220 H=4,4*L3,4 M=H+4*L3*(L2*(I-1)+(J-1)) WRITE (7, 215) I,J,H,M,TK(0, I , J,H) ,TK(TIME(1) , I , J,H) , + (TK(IU,I,J,H), IU=TIME(2) , TIME (15) ,DT) 215 FORMAT (1X,4(I3,1X),16(1X,F5.2)) 220 CONTINUE I F ( IT . GT . 16) THEN 226 WRITE(7,230)(TIME(IC),IC=16,30) 230 FORMAT (9X, 'TIME(MIN) ' , IX, 16(F4.0,2X) ) DO 240 1=1, LI DO 240 J=1,L2 DO 240 H=4,4*L3,4 M=H+4*L3*(L2*(I-1)+(J-1) ) WRITE (7 , 235) I , J,H,M, (TK(IU, I, J,H) , I U=TIME (16) ,TIME(30) , DT) 235 FORMAT (1X,4(I3,1X),15(1X,F5.2)) 240 CONTINUE 246 WRITE(7, 250) (TIME(IC) , IC-31, IT-1) 250 FORMAT (9X, 'TIME (MIN) ' , IX, 16(F4.0, 2X) ) DO 260 1=1, LI DO 260 J=1,L2 DO 260 H=4,4*L3,4 104 M=H+4*L3*(L2*(I-1)+(J-1)) WRITE (7, 255) I,J,H,M,(TK(IU,I,J,H), IU=TIME(31) ,TIME( IT-1) ,DT) 255 FORMAT (1X,4(I3,1X),16(1X,F5.2)) 260 CONTINUE ENDIF STOP END SUBROUTINE MIDNOD( AFP, DFP, APP, MI , 1 1 , IF, JI , JF, HI , HF, + WF,HWF,WG,NF,NA) C THIS SUBROUTINE GENERATES MATRIX COEFFICIENTS FOR ENERGY EQUATION C IN FRUIT INTERMEDIATE NODES. COMMON A (384, 384) ,TK(0:840,4,6,20) , C(384) , KF , RF ( 384 ) ,DL(384) , + VI(0:30,6,6,-5:17),HW(6),JET(6),DLE(6),V(6),TO,AS(384), + Q,DT,L1,L2,L3,HA(4,6, 17) , HOR (6 , 4 , 6 , 1 7 ) ,AN(384,0: 17,0:17), + AV(384, 0: 17), AC(384), AR(384, 0:16,0: 17), RN(384, 0:17), + RE (384,0: 17) ,TW(0:600,4,6, 17) ,HFR(6,4,6, 17) , VO (840) DIMENSION NF(3) , NA(6) , HWF(6) ,HC(400) ,HCW(400) , DFP (600) , APP (600) INTEGER WG ( 5 ) , WH ( 5 ) , WA ,WF,WI,HI,HF,H,Q,F,P REAL KF,KP,KW(6) PARAMETER (PI-3.1416) DATA HRO/1 . 2 5/ , HR 1/ 1 .2 6/ M=MI DO 200 I = 1 1 , IF DO 150 J=J I , JF DO 100 H=H I , HF , 4 AFPC=AFP*(60.0/DFP(M) ) APP(M)=AFPC QR=- - 1472+ (4.326E-3) *TK (Q, I,J,H) DO 10 WI=1 ,6 HCW(M+1 )=12*(KF/ (RF(M+1 ) -RN(M+2 , H+l ) ) ) HC(M+1)=2*HCW(M+1 ) I F ( J ET ( W I ) . EQ.5)THEN HWF(WI)=0.0 ELSEIF(HWF(WI) .NE.O. )THEN HWF(WI)= 1 - 3* ( 1/HOR ( W I , I , J, H-2)+l/HCW(M+l)+l/HW(WI) )**-l ELSE HWF ( W I ) =0 . 0 ENDIF HFR(WI, I, J,H)=HWF(WI) CONTINUE A(M,H+4*L3*(L2*( I- l)+( J-l) ) )=60/(AFPC*DT)+ (144*(AN(M+1,H,H-1)-AC(M+1)))/(AV(M+1,H)*AR(M+1,H,H-1))+ (144*AN(M+1 ,H,H+1) )/ (AV(M+1 , H)*AR(M+1 , H, H+l ) ) +( ( 12*AC (M+l ) )/ (KF*AV(M+1 ,H) ) )*(WF*HC(M+1) +HWF ( 1 )+HWF(2)+HWF(3)+HWF(4)+HWF(5)+HWF(6) ) A(M, (H-l )+4*L3*(L2*( I-1)+(J-1)))= - ( 144* (AN (M+l ,H,H-1)-AC(M+1)))/(AV(M+1 , H)*AR(M+1 , H,H- 1 ) ) A(M, (H+1)+4*L3*(L2*( I-1)+(J-1) ) )= -(144*AN(M+1,H,H+1))/(AV(M+1,H)*AR(M+1,H,H+1)) I F ( H . LT . (4*L3-1) ) THEN A(M, (H+4)+4*L3*(L2*( I-1)+(J-1) ) )= - ( ( 12*HC(M+1 )*AC(M+1 ) )/ (KF*AV(M+1 , H) ) )*(NF(3)+NA(4) ) ENDIF IF(H.GT.3)THEN A(M, (fl-4)+4*L3*(L2*( I-1)+(J-1) ) )= - ( ( 12*HC(M+1 )*AC(M+1 ) )/(KF*AV(M+l ,H) ) )*(NF(3)+NA(5) ) ENDIF IF(I.LT.Ll) THEN A(M,H+4*L3*(L2*I+(J- 1 ) ) )= - ( ( 12*HC(M+1 )*AC(M+1 ) )/ (KF*AV(M+1 ,H)))*(NF( 1 )+NA( 1 ) ) ENDIF IF(I.GT.l) THEN A(M,H+4*L3*(L2*(I-2) + (J-l) ) ) = -( (12*HC(M+1)*AC(M+1) )/ (KF*AV(M+I ,H) ))*(NF(1)+NA(6) ) ENDIF IF(J.LT.L2)THEN A(M,H+4*L3*(L2*(I-1)+J))= -((12*HC(M+1)*AC(M+1) )/(KF*AV(M+l ,H) ) )*(NF(2)+NA(2) ) ENDIF I F ( J . GT . 1 ) THEN A(M,H+4*L3*(L2*( I-l)+(J-2)))= - ( ( 12*HC(M+1 )*AC(M+I ) )/ (KF*AV(M+1 , H) ) )*(NF(2)+NA(3) ) ENDIF C(M)=(60*TK(Q, I ,J,H) )/ (AFPC*DT)+(DFP(M)*QR)/KF+ ( (12*AC(M+1)*T0)/ (KF*AV(M+1 ,H) ) )*(HWF(1)+ HWF ( 2 ) +HWF ( 3 ) +HWF ( 4 ) +HWF ( 5 ) +HWF ( 6 ) ) 107 M=M+4 100 CONTINUE M=WG( l)*(M+8)+WG(2)*(M+4*L3-4) 150 CONTINUE M=WG(3)*{M+8*L3)+WG(4)*(M+4*L3*(L2-l)+8)+WG(5)*(M+4*L3*(L2-l)) 200 CONTINUE RETURN END C C SUBROUTINE CENOD(AFP,DFP,APP) THIS SUBROUTINE GENERATES MATRIX COEFFICIENTS FOR ENERGY EQUATION IN FRUIT CENTER NODES. COMMON A (384, 384) , TK(0 .*840,4,6,20) , C(384) ,KF,RF(384) ,DL(384) , VI(0:30,6,6,-5:17),HW(6),JET(6),DLE(6),V(6) , TO, AS (384) , Q , DT ,L1,L2,L3, HA (4,6, 17) ,HOR(6, 4,6, 17) , AN (384,0: 17,0:17), AV(384,0: 17), AC(384),AR(384, 0:16,0: 17), RN(384, 0:17), RE (384,0: 17) ,TW(0:600,4,6, 17) ,HFR(6,4,6, 17) , VO (840) DIMENSION DFP(600) , APP(600) INTEGER WG(5) ,WH(5) ,WA,WF,WI,HI,HF,H,Q,F,P REAL KF,KP,KW(6) PARAMETER (PI=3. 1416) F=1 M=4 DO 100 1=1, LI DO 100 J=1 , L2 DO 100 H=4,4*L3,4 108 AFPC=AFP*(60.0/DFP(M) ) APP(M)=AFPC QR=- - 1472+ ( 4 . 326E-3) *TK( Q , I , J , H ) A(M,H+4*L3*(L2*( I-1)+(J-1) ) )=60/(AFPC*DT)+ + ' ( 144*AN(M, H,H- 1 ) )/ (AV(M, H)*AR(M, H, H- 1 ) ) A(M, (H-1)+4*L3*(L2*( I-1)+(J-1) ) )= + - ( 144*AN (M, H , H- 1 ) )/ ( AV (M, H ) *AR (M , H , H- 1 ) ) C(M)=(60*TK(Q,I,J,H))/(AFPC*DT)+(DFP(M)*QR)/KF M=M+4 100 CONTINUE RETURN END SUBROUTINE SUFNOD (AFP,DFP,APP) C THIS SUBROUTINE GENERATES MATRIX COEFFICIENTS FOR ENERGY EQUATION C IN FRUIT SURFACE NODES. COMMON A (384, 384) ,TK(0:840,4,6,20) , C ( 384 ) ,KF,RF(384) ,DL(384) , + VI(0:30,6,6,-5:17),HW(6),JET(6),DLE(6),V(6) , TO, AS (384) , + Q , DT, LI , L2 , L3 , HA(4 ,6,17) ,H0R(6,4,6, 17) , AN (384,0: 17,0: 17) , + AV(384, 0 : 17), AC(384) ,AR(384, 0:16,0: 17), RN(384, 0:17), + RE (384,0: 17) , TW(0:600,4,6, 17),HFR(6,4,6,17) , VO (840) DIMENSION NF(3) , NA(6) , HWF(6) , DFP(600) , APP(600) INTEGER WG(5),WH(5),WA,WF,WI,HI,HF,H,Q,F REAL KF,KP,KW(6) PARAMETER (PI=3. 1416) DATA HR0/.85/,HRI/.86/ 109 M=2 DO 100 1=1, LI DO 100 J=1 , L2 DO 100 H=2,4*L3-2,4 A(M,H+4*L3*(L2*(I-1)+(J-1)))= + ( HA ( I , J, H- l)*AS(M+2) )/KF+ + 1 44* ( AN ( M+2 , H+ 1 , H) -AC (M+2) )/ (AR(M+2, H+l , H)*DL(M+2)**3) A(M, (H-1)+4*L3*(L2*(I-1)+(J-1)))= + -(HA(I,J,H-l)*AS(M+2))/KF A(M, (H+1)+4*L3*(L2*(I-1)+(J-1)))= + - 144* (AN (M+2 , H+l , H) -AC(M+2) )/ (AR(M+2 , H+l , H)*DL(M+2)**3) C(M)=0.0 M=M+4 100 CONTINUE RETURN END SUBROUTINE AIRNOD(MI, II, IF,JI,JF,HI,HF,WA,HWA,WH,NF,NA,NI) C THIS SUBROUTINE GENERATES MATRIX COEFFICIENTS FOR ENERGY EQUATION C IN AIR NODES. COMMON A (384, 384) ,TK(0:840,4,6,20) ,C(384) , KF , RF (384 ) ,DL(384) , + VI(0:30,6,6,-5:17),HW(6),JET(6),DLE(6),V(6) ,T0, AS(384) , + Q , DT ,L1,L2,L3, HA (4,6, 17) ,H0R(6, 4,6, 17) , AN (384,0: 17,0:17), + AV(384, 0 : 17), AC(384), AR(384, 0:16,0: 17), RN(384, 0:17), + RE (384,0: 17) ,TW(0:600,4,6, 17) ,HFR(6,4,6, 17) , VO (840) DIMENSION NF(3) ,NA(6) , HWA ( 6 ) 110 INTEGER WG(5),WH(5),WA,WF,WI,HI,HF,H,Q,F,P REAL KA , KF , KP , KW ( 6 ) , N I ( 0 : 6 ) PARAMETER (PI=3. 1416) DATA HRO,HRI,P,AAV,AAP,KA,E/l .25, 1 . 26, 2, 1 . 23 , . 79, . 0148, . 476/ M=MI DO 200 I = 1 1 , IF DO 150 J=J I , JF DO 100 H=HI,HF,4 KP=( (1 . 16E-5)*E**3*(DL(M+3)**2) )/( (1-E)**2) HA ( I, J,H)= . 4675* (ABS( (TK(Q,I,J,H+1)-TK(Q,I,J,H) )/DL(M+3) ) )**( 1 ./4) + . 3905/DL (M+3 ) +1 . 1*HRI VI(Q,I,J,H)=(1 .376E+6)*KP*ABS(TK(Q, I , J,H) -TK(Q, I , J,H+4) ) VI (Q, I, J, -3)=0.0 I F ( ( I . NE . 1 ) .AND. (I. NE. LI)) THEN IF( (J.NE. 1) .AND. (J .NE. L2) )THEN VI (Q, I, J,4*L3-3)=VI(Q, I,J,4*L3-7) ENDIF ENDIF I F ( ( I . EQ . 1 ) - OR . (I.EQ.L1)) THEN VI (Q , I , J , H) =V0 (Q) ELSEIF((J.EQ.l) .OR. (J.EQ.L2)) THEN VI (Q, I , J,H)=VO(Q) ENDIF DO 10 HI-1,6 I F ( HWA ( WI ) .NE.O. )THEN IF( JET (WI ) . EQ . 5) THEN HWA(WI)=0.0 ELSE TW(Q, I , J,H)=(TO+TK(Q, I , J,H) )/2. IF(v)ET(WI) .EQ.4)THEN I F ( WI . EQ . 4) THEN H0=.28*((TW(Q,I,J,H)-T0)/DLE(WI))**(l./4) ELSEIF(WI.EQ.5) THEN H0= .27*(TW(Q,I,J,H)-T0)**(l./3) ELSE H0=. 19*(TW(Q, I , J,H) -T0)**(1 ./3) ENDIF IF(HO.LT. .99) H0=.99 ELSE H0FP=.99+. 21*V(WI) H0FNP=3 .48*(V(WI )**. 731 )/DLE(WI)**( .269) H0FNC=1 . 1*(V(WI)**.675)/DLE(WI)**( .325) H0FNS=2. 5*(V(WI)/DLE(WI) )**( 1 ./2) HOFSH= . 3/DLE ( WI ) +1 . 4*(V(WI)/DLE(WI) )**( 1 ./2) +.6*V(WI)**.67/DLE(WI)**(l./3) I F (JET (WI ) .EQ.3)THEN H0=(H0FP+H0FNC+H0FSH)/3 . ELSEIF (JET (WI ) .EQ. 1 ) THEN HO=(HOFNP) ELSEIF (JET (WI ) .EQ.2)THEN H0=(H0FNS+H0FNC+H0FSH)/3 . 112 ENDIF ENDIF I F ( WI . EQ . 4 ) THEN HIW= . 689*(TK(Q, I,J,H)-TW(Q,I,J,H))**(l./3.) ELSEIF(WI.EQ.5)THEN HIW=.27*(TK(Q, I,J,H)-TW(Q,I,J,H))**(l./3.) ELSE HIW=. 19*(TK(Q, I , J,H) -TW(Q, I , J,H) )**(1 ./3. ) ENDIF HOR(WI , I , J,H)=HO+HRO H I R=H I W+HR I HWA(WI ) = 1 . 3* ( 1/HOR ( W I , I , J, H) + 1/HW(WI ) + l/HIR)**( - 1 ) ENDIF ELSE HWA(WI)=0.0 HOR ( W I , I,J,H)=0.0 ENDIF 10 CONTINUE A(M,H+4*L3*(L2*( I -l)+( J- 1 ) ) )=60*E/(AAV*DT) + +12*VI (Q, I , J , H )/ (AAP*DL(M+3) )+144*WA/(DL(M+3)**2) + +(HA( I , J, H)*AS(M+3) )/KA+ + (12/ (KA*DL(M+3) ) )*(HWA( 1 )+HWA(2)+HWA(3)+HWA(4)+HWA(5)+HWA(6) ) A(M, (H+1)+4*L3*(L2*(I-1)+(J-1) ))= -(HA(I,J,H)*AS(M+3))/KA IF(H.GT. 1)THEN A (M , (H-4)+4*L3*(L2*( I-1)+(J-1) ) )= + -12*VI(Q,I,J,H)*NI(0)/ (AAP*DL(M+3) ) - ( 144/DL(M+3)**2)*(NF(3)+NA(5) ) ENDIF IF(H. LT. (4*L3-3) )THEN A(M, (H+4)+4*L3*(L2*( I-1)+(J-1) ) )= - (144/DL(M+3)**2)*(NF(3)+NA(4) ) - 1 2*V I ( Q , I , J,H)*NI (1)/ (AAP*DL(M+3) ) ENDIF IF(I.LT.Ll) THEN A(M, H+4*L3*( L2*I+( J-l ) ) )= - ( I44/DL (M+3)**2)*(NF( 1 )+NA( 1 ) ) -12*VI (Q, I , J,H)*NI(2)/ (AAP*DL(M+3) ) ENDIF I F ( I . GT . 1 ) THEN A(M,H+4*L3*(L2*(I-2)+(J-l)))=-(144/DL(M+3)**2)*(NF(l)+NA(6)) -12*VI (Q, I , J,H)*NI (3)/(AAP*DL(M+3) ) ENDIF IF(J. LT. L2)THEN A(M,H+4*L3*(L2*( I-l)+J) )= - ( 1 44/DL ( M+3 ) **2 ) * ( NF ( 2 ) +NA ( 2 ) ) - 1 2*V I ( Q , I , J,H)*NI (4)/(AAP*DL(M+3) ) ENDIF IF(J.GT. 1 ) THEN A(M,H+4*L3*(L2*( I-l)+(J-2) ) )=- (144/DL(M+3)**2)*(NF(2)+NA(3)) -12*VI (Q, I, J,H)*NI (5)/(AAP*DL(M+3) ) ENDIF C(M)=(60.*E*TK(Q, I, J,H))/ (AAV*DT) +((12*T0)/(KA*DL(M+3)))*(HWA(1)+HWA(2)+HWA(3)+ HWA(4)+HWA(5)+HWA(6) ) 114 M=M+4 100 CONTINUE M=WH( 1 )*(M+8)+WH(2)*(M+4*L3-4) 150 CONTINUE M=WH(3)*(M+8*L3)+WH(4)*(M+4*L3*(L2-l)+8)+WH(5)*(M+4*L3*(L2-l ) ) 200 CONTINUE RETURN END SUBROUTINE MATRIX(*,MT,NT, INDEX, T) C SUBROUTINE MATRIX USES A VARIATION OF THE GAUSS-SEIDEL INDIRECT C METHOD TO SOLVE THE SYSTEM OF ENERGY EQUATIONS IN IMPLICIT FORMAT C FOR THE 3-DIMENSIONAL TEMPERATURE RESPONSE OF A BIN OF SPHERICAL C FRUIT. COMMON A (384, 384) ,TK(0:840,4,6,20) ,C(384) ,KF,RF(384) , DL ( 384 ) , + VI(0:30,6,6,-5:17), HW ( 6) , JET (6) ,DLE(6) ,V(6) , TO, AS (384) , + Q , DT ,L1,L2,L3, HA (4,6, 17) ,H0R(6,4,6, 17) , AN (384,0: 17,0:17), + AV(384, 0: 17), AC(384) ,AR(384, 0:16,0: 17), RN(384, 0:17), + RE (384, 0 : 17) ,TW(0 : 600,4 , 6, 17), HFR (6, 4, 6, 17), VO (840) DIMENSION T (1200) TERRMAX=0.005 TL00P=0 L00p=10000000 9 TERR=0.0 10 DO 100 IA=1 ,MT,4 SUMA=0.0 115 SUMM=0 . 0 I S= I A+l IM=IA+2 IC=IA+3 TA=f (IA) TS=T (IS) TM=T(IM) TC=T ( IC) DO 50 J=1 ,NT IF(J.NE. IA)THEN SUMA=SUMA+A( IA, J)*T (J) IF(J.NE. IM)THEN SUMM=SUMM+A( IM, J)*T (J) END IF ENDIF 50 CONTINUE T(IA)=(1 . 0/ ( A ( I A , IA)))*(C(IA)-SUMA) TERRA=ABS(TA-T ( IA) ) IF(TERRA.GT.TERR)TERR=TERRA T(IM)=(1.0/(A(IM,IM)))*(C(IM) -SUMM-A( IM,IA)*T(IA)) TERRM=ABS(TM-T ( IM) ) IF(TERRM.GT.TERR)TERR=TERRM T( IS)=( 1 .0/(A( IS, IS) ) )*(C( IS) -A( IS, IM)*T( IM) - + A( IS, IA)*T ( IA) ) TERRS=ABS(TS-T ( IS) ) IF (TERRS. GT. TERR) TERR=TERRS 116 I(IC)=(1.0/(A(IC,IC)))*(C(IC)-A(IC,IM)*T(IM)) TERRC=ABS(TC-T ( IC) ) IF(TERRC.GT.TERR)TERR=TERRC 100 CONTINUE TL00P=TL00P+1 I F (HOOP. GE. LOOP) THEN WRITE ( 7 , *) ' EXCESIVE NUMBER OF ITERATIONS' RETURN 1 ENDIF IF(TERR.GT.TERRMAX)GO TO 9 C WRITE(7,*)TL00P RETURN END APPENDIX B INPUT PROGRAM: PROPERTIES AND DIMENSIONS INPUT 3-D TRANSIENT TEMPERATURE DISTRIBUTION SOLUTION IN A FRUIT BIN C THIS PROGRAM IS DESIGNED TO ALLOW THE INPUT OF PHYSICAL PROPERTIES AND C DIMENSIONS OF A BIN OF FRUIT AND ROOM COOLING STREAM TEMPERATURE AND C VELOCITIES AROUND THE FRUIT BOX. WITH THE INFORMATION ENTERED, THIS C PROGRAM GENERATES A DATA FILE THAT CAN THEN BE USED TO RUN THE MAIN C PROGRAM IN BATCH MODE. COMMON KF,RF,HW(6),JET(6),DLE(6),V(6) DIMENSION DLW(6) INTEGER WI REAL KF,KP,KW(6) CHARACTER SELE*2 OPEN(UNIT=2, F I LE=' BFRU.DAT' , STATUS='NEW' ) REW I ND ( 2 ) C THIS IS THE BOX PROPERTIES INPUT BLOCK. AT THE PROMPT, ENTER THE C PROPERTIES FOR ONE OF THE SIDE WALLS THAT YOU SELECT AS "NORTH C WALL". THEN, AS PROMPTED, ENTER VALUES FOR THE REMAINING 3 SIDE C WALLS (SOUTH, EAST AND WEST) CONSISTENTLY WITH THE ORIENTATION OF C YOUR SELECTED NORTH WALL. SUBROUTINE AIRJET IS INVOKED TO ASSIGN C THE ORIENTATION INDEXES AND CODES THAT IDENTIFY DIMENSIONS AND C ROOM AIR VELOCITIES WITH THE SIX BOX WALLS. 10 W 1=2 50 WRITE (5,55) 117 118 55 FORMAT (7X, 'ENTER FOR NORTH WALL:',/, + 15X, ' LENGTH(DLl) -INCHES) ' + 15X, 'HEIGHT (DL3) - INCHES) ' ,/, + 15X, ' TH I CKNESS (DLW- INCHES) ' ,/, + 15X; 'THERMAL CONDUCTIVITY(KW) -BTU/HR-FT-F' ) READ* , DL 1 , DL3 , DLW ( 2 ) , KW ( 2 ) DLE(2)= DL1*(2.*DL3)/ (DL1+DL3) CALL AIRJET(WI,KW,DLW) WI = 1 WRITE ( 5 , 60) 60 FORMAT (7X, 'ENTER FOR EAST WALL:',/, + 15X, ' LENGTH (DL2) - INCHES' ,/, + 15X,'THICKNESS(DLW)-INCHES',/, + 15X, 'THERMAL CONDUCTIVITY(KW) -BTU/HR-FT-F' ) READ* , DL2 , DLW ( 1 ) , KW ( 1 ) DLE(1)=DL2*(2.*DL3)/ (DL2+DL3) CALL AIRJET(WI,KW,DLW) W I =3 WRITE (5,65) 65 FORMAT (7X, 'ENTER FOR SOUTH WALL:',/, + 15X,'THICKNESS(DLW)-INCHES',/, + 15X, 'THERMAL CONDUCTIVITY(KW) -BTU/HR-FT-F' ) READ*, DLW(3),KW(3) DLE(3)=DLE(2) CALL AIRJET(WI,KW,DLW) W I =6 119 WRITE (5,70) 70 FORMAT (7X, ' ENTER FOR WEST WALL:',/, + 15X, ' THICKNESS (DLW) - INCHES' ,/, + 15X, 'THERMAL CONDUCT I VITY( KW) -BTU/HR-FT-F' ) READ*,DLW(6) ,KW(6) DLE(6)=DLE(1) CALL AIRJET(WI,KW,DLW) W I =4 WRITE(5, 75) 75 F0RMAT(7X, 'ENTER FOR BOTTOM WALL:',/, + 15X, 'THICKNESS(DLW)-INCHES' ,/, + 15X, 'THERMAL CONDUCTIVITY(KW) -BTU/HR-FT-F' ) READ*,DLW(4),KW(4) DLE(4)=DL1*(2.*DL2)/ (DL1+DL2) CALL AIRJET(WI , KW, DLW) WI=5 WRITE(5, 77) 77 FORMAT (7X, 'ENTER FOR TOP WALL:',/, + 15X, 'THICKNESS (DLW) - INCHES' ,/ + 15X, 'THERMAL CONDUCTIVITY(KW) -BTU/HR-FT-F' ) READ*,DLW(5) ,KW(5) DLE(5)=DLE(4) CALL AIRJET (WI ,KW,DLW) C NEXT IS THE TOMATOES AVERAGE PROPERTIES INPUT BLOCK. PROPERTIES C OF INDIVIDUAL FRUITS ARE ENTERED IN THE SPECIAL DATA FILES ACCESSED C BY THE MAIN PROGRAM. 120 WRITE ( 5 , 80 ) 80 FORMAT ( 7X , ' ENTER FOR TOMATOES : ' , / , + 15X, 'THERMAL D I FFUS I V ITY ( AFP ) -FT2/HR: ' ,/, + 15X, ' DENS ITY ( DFP) -LBM/FT3 : ' ,/, + 15X, ’'THERMAL CONDUCTIVITY(KF) -BTU/HR-FT-F: ' ,/, + 15X, 'AVERAGE RADIUS ( RF) - INCHES : ' ,/, + 15X, 'AVERAGE SINGLE CONTACT AREA FRACTION (FAC) : ' ) READ* , AFP , DFP , KF , RF , FAC WRITE(5,90) 90 F0RMAT(7X, 'ENTER:',/, + 15X, 'ROOM OR AIR STREAM TEMPERATURE (TO) -F: ' ,/, + 15X, ' INITIAL TOMATOES TEMPERATURE(TF) -F: ' ,/, + 15X, 'TERMINAL STORAGE (END OF COOLING) TEMPERATURE(TS) -F: ' ,/, + 15X, 'TIME INCREMENT (EVERY HOW MANY MINUTES TEMPERATURES',/, + 15X, 'ARE WANTED) (DT>20MIN) -MIN: ' ) READ*, TO,TF,DT WRITE(2,*) DL 1 , DL2 , DL3 WRITE (2,*) (DLE(I) ,1=1,6) WRITE (2,*) (HW(I) , 1=1 ,6) WRITE (2,*) (V(I) , 1=1,6) WRITE (2,*) (JET ( I ) , 1=1,6) WRITE ( 2 , * ) AFP,DFP,KF,RF,FAC WRITE(2,*) TO,TF,DT WRITE(5,95) 95 F0RMAT(2X, 'YOU HAVE FINISHED ENTERING NECESSARY DATA.',/, + 2X, 'NEXT STEP IS TO TYPE "SUBMIT BFRUIT" AT THE',/, 121 + 2X, ' PROMPT "$" TO SOLVE TEMPERATURE DISTRIBUTION') STOP END SUBROUTINE AIRJET(WI , KW, DLW) COMMON KF,RF,HW(6),JET(6),DLE(6),V(6) DIMENSION NF(9) ,HWF(6) ,HWA(6) ,DLW(6) INTEGER WG(5),WH(5),WA,WF,WI,HI,HF,H,Q,F,P REAL KF,KP,KW(6) PARAMETER (PI=3. 1416) WRITE ( 5 , 10) 10 F0RMAT(7X, 'AIR STREAM DIRECTION( JET) : ' + 10X, 'WINDWARD(TYPE 1) , LEEWARD(TYPE 2) , PARALLEL(TYPE 3)',/, + 10X, 'STILL AIR(TYPE 4)',/, + 10X, ' BOXWALL NEXT TO WELL INSULATED WALL OR ON FLOOR(TYPE 5)') READ*, JET(WI) IF( ( JET(WI) . EQ . 4) .OR. ( JET(WI ) . EQ. 5) )THEN V(WI)=0.0 ELSE WRITE (5,20) 20 F0RMAT(7X, 'AIR STREAM VELOCITY(V) -FPS: ' ) READ*, V(WI) ENDIF HW(WI)=12.*KW(WI)/DLW(WI) RETURN END APPENDIX C INPUTS TO THE NUMERICAL MODEL Box dimensions Room air speed Time increment Initial bin temperature Storage temperature Room air temperature Number of fruit contact points Packing density Fruit thermal conductivity Fruit thermal diffusivity Fruit average density 40 cm x 27.3 cm x 26.7 cm high 18 - 120 m/min. 5-20 min. 24 - 26 °C 13.1 °C 8.1 °C 6 - 8 46 - 60% 0.53 W/m °C 1.37 x 10'7 m2/s 961.6 Kg/m3 122 BIBLIOGRAPHY Bakker-Arkema, F.W. and W.G. Bickert. 1966. A deep bed computational cooling procedure for biological products. Transactions of the ASAE 9(6): 834-836. Baird, C.D., J.J. Gaffney and D.T. Kinard. 1975. Research facility for forced-air precooling of fruits and vegetables. ASAE technical paper 72-364. Bellagha, S. 1985. Heat and mass transfer during the cooling of tomatoes individually and in bulk. M. S. Thesis. University of Florida. Beukema, K.J. and S. Bruin. 1983. Three-dimensional natural convection in a confined porous medium with internal heat generation. International Journal of Heat and Mass Transfer (26)3: 451-458. Beukema, K.J., S. Bruin and J. Schenk. 1982. Heat and mass transfer during cooling and storage of agricultural products. Chemical Engineering Science 37(2): 291-298. Buffington, D.E. and S.K. Sastry. 1983. Methodology for determining the most economic storage conditions for tomatoes. Revue Internationale du Froid, 6(4): 213-218. Carnahan, B., H.A. Luther and J. Wilkes. 1969. Applied numerical methods. John Wiley and Sons. New York. 485. Chau, K.V. and S. Bellagha. 1985. Heat and mass transfer during the cooling of tomatoes individually and in bulk. ASAE technical paper 85-6002. Chau, K.V., J.J. Gaffney and S. Bellagha. 1984. Simulation of heat and mass transfer in products with internal heat generation and transpiration. ASAE technical paper 84-6513. Fahien, R.W. 1983. Fundamentals of heat transport phenomena. McGraw- Hill. New York. 459-465. Guillou, R. 1960. Coolers for fruits and vegetables. California Agricultural Station. Bulletin 773. 1-6. 123 124 Hardenburg, R.E., A.E. Watada and C.Y. Wang. 1986. The commercial storage of fruits, vegetables, and florist and nursery stocks. Agricultural handbook 66. United States Department of Agriculture, Washington, D.C. Holman, J.P. 1981. Heat transfer (4th ed). McGraw-Hill. New York 240-247. Hornbeck, R.W. 1975. Numerical methods. Prentice-Hall. Enqlewood Cliffs, New Jersey. 104-276. Incropera, F.P. and D.P. DeWitt. 1985. Fundamentals of heat transfer. John Wiley and Sons. New York. 437. Johns, W.R. and S.J. Lawn. 1985. Unsteady-state transfer between a sphere and a surrounding stationary medium with application to arrays of spheres. International Heat and Mass Transfer Journal 28(5): 1047-1053. Kader, A., R.F. Kasmire, F.G. Mitchell, M. Reid, N. Sommer and J. Thompson. 1985. Postharvest technology of horticultural crops University of California, Davis. 3-38, 141. Lerew, L.E. 1978. Development of temperature-weight loss model for bulk stored potatoes. Ph.D. dissertation. Michigan State University. East Lansing. Mitchell, F.G. , R. Guillou and R.A. Parsons. 1972. Commercial cooling of fruits and vegetables. Agricultural publications. University of California, Davis. 3-33. Mueller, G.E. 1988. Analytical solution for transient distribution of fluid and solid temperatures of vol umetrical ly heated fixed porous nuclear debris. Nuclear Engineering Design. 106(2): 231-241. Mural idhar, K., R.A. Baunchalk and F.A. Kulachi. 1986. Natural convection in a horizontal porous annulus with a step distribution in permeability. Journal of Heat Transfer Transactions. 108(4)- 889-893. v ' Peleg, K. Produce handling, packaging and distribution. 1985. AVI Publishing Company. Westport, Connecticut. 135-161. Rowe, P.N. and K.T. Claxton. 1965. Heat and Mass Transfer from a single sphere to fluid flowing through an array. Transactions of the Institution of Chemical Engineers. 43: T321-T333. Reddy, G.B and J.C. Mulligan. 1987. Macroscopic continuum analysis of simultaneous heat and mass transfer in unsaturated porous materials containing a heat source. International Communications in Heat and Mass Transfer. 14(3): 251-263. 125 Romero, R.A. 1987. Transpiration from fruits and vegetables in storage. Department of Agricultural Engineering. University of Florida, Gainesville. Talbot, M.T. 1987. Pressure and velocity distribution for air flow through fruits packed in shipping containers using porous media flow analysis. Ph.D. dissertation. University of Florida, Gainesville, 101-129. * Thompson, J.F., Z.U.A. Warsi and C.W. Martin. 1985. Numerical grid generation. North-Holland. New York. 51-56. Trevisan, 0. and A. Bejan. 1985. Natural convection with combined heat and mass transfer buoyancy effects in a porous medium. International Journal of Heat and Mass Transfer 28(8): 1597-1611. Turczyn. M. and J.P. Anthony. 1986. MUM replacement containers. Project MUM. United States Department of Agriculture. Washington, DC. Welty, J.R. 1978. Engineering heat transfer. John Wiley and Sons. New York. 38-65. Yuge, T. 1960. Experiments on heat transfer from spheres including combined natural and forced convection. Transactions of the ASAE 82: 214-220. BIOGRAPHICAL SKETCH Tomas Bazan was born in Panama City, Panama, where he completed his elementary and high school education. In 1975, he received his Bachelor of Science in Mechanical Engineering at the University of Panama. He became a professor and was a principle actor in creating the Department of Mechanical Engineering at the Technological University of Panama before earning his Master of Science in Mechanical Engineering at the University of Wisconsin at Madison in 1982. He began doctoral studies at the University of Florida in 1986 as a Fulbright scholar while on leave of absence from the Technological University of Panama. 126 I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Robert B. Gaither, Chairman Professor of Mechanical Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. KP Khe V. Chau, Cochairman Associate Professor of Agricultural Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for~tbe deqree of Doctor of Philosophy. Herbeflr-Ar^IngTey Associate Professor Engineering Meehan i cap I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the deqree of Doctor of Philosophy. C. Direlle Baird Professor of Agricultural Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. GJJLL^- P William E. Lear Assistant Professor of Mechanical Engineering ThTS dissertation was submitted to the Graduate Faculty of the college of Engineering and to the Graduate School and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy. August, 1989 Q. • ‘ Dean, College of Engineering Dean, Graduate School