[B - IHIIHIIIIHlllllllulllllflllllllfIIHllllllllllllllllllllll 293 10376 5156 This is to certify that the thesis entitled A NUMERICAL PREDICTION MODEL FOR FOOD FREEZING USING FINITE ELEMENT METHODS presented by Hadi Karia Purwadaria has been accepted towards fulfillment of the requirements for _Eh.._D_.__degree in Wral Engineering 9. 6.4% Major professor Date tC. 28’ 14.19 0-7 639 M .1535 e o r I “flip "‘ 2 .‘JAN Ure'flg'ég -. ‘4 '1" " z .142 W ‘6')!“ a 0 3 7 WW1 5351‘ e 3 . 34 S I l OVERDUE FINES: 25¢ per day per item RETUMIIG LIIRARY MTERIALS: Place in book return tore-0 charge from circulation records A NUMERICAL PREDICTION MODEL FOR FOOD FREEZING USING FINITE ELEMENT METHODS BY Hadi Karia Purwadaria A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Agricultural Engineering 1980 ABSTRACT A NUMERICAL PREDICTION MODEL FOR FOOD FREEZING USING FINITE ELEMENT METHODS BY Hadi Karia Purwadaria The rate of freezing is one of the most important factors in designing an efficient freezing process for foods in order to achieve good product quality and to avoid excessive energy consumption. Significant improvements have been achieved in the area of freezing process simula- tion, however, the phenomena of phase—change, the influence of product thermal properties, the importance of product geometry and the effect of freezing environment on the freezing rate are not fully understood. The objective of this investigation was to develop a numerical simulation model using the finite element method to predict freezing rate in anomalous food product geometries while accounting for the non-linear temperature dependent product properties and various boundary condi- tions. To verify the model, experimental tests were con- ducted for elliptical and trapezoidal product shapes using ground beef as the food product. The experiments Hadi Karia Purwadaria were conducted in a wind tunnel placed in a low temperature room and temperature measurement was recorded for 24 node locations within the critical cross-section of the product during the freezing process. The finite element computer simulation used to predict the food product freezing rate of anomalous shape has been developed and verified by experimental data. The results illustrate the capability of the simulation model to incorporate various boundary conditions and various product geometries. Closer approximation to the experi- mental data was obtained by using the prediction incor- porating a boundary condition with the surface heat transfer coefficient varying as a function of location. More efficient freezing times are predicted by utilizing an approach based on area average enthalpy as compared to the conventional method based on the slowest freezing point location. Time steps in the range from one to three minutes do not influence the stability of the finite ele- ment scheme. While geometric size has significant influ- ence on the rate of freezing, the influence on initial product temperature in the range from 14.0 to 22.0°C interval is negligible. Approved 1 Major Professor Approved ACKNOWLEDGMENTS The author is deeply indebted to Dr. Dennis R. Heldman, Professor of the Agricultural Engineering Depart- ment and Department of Food Science and Nutrition, for suggestion of the research tOpic and his continuous support, interest, and patient guidance throughout the course of this study. Sincere appreciation is extended to the guidance committee members: Dr. Larry J. Segerlind, Agricultural Engineering Department; Dr. K. Jayaraman, Chemical Engin- eering Department; and Dr. Lawrence E. Dawson, Department of Food Science and Nutrition; and to the external exam- iner, Dr. James F. Steffe, Agricultural Engineering Department. Their review and suggestions to the manu- script have been the most helpful. Special thanks go to Dr. Haruhiko Murase, Research Associate in Agricultural Engineering Department, for his helpful suggestions in the laboratory experiments. The author also wishes to express his gratitude to Dr. Andi Hakim Nasoetion, President of Bogor Agricul- tural University, for his encouraging support and to the Government of Indonesia which provided the opportunity for the author to continue his study. ii The contribution from Dr. Merle L. Esmay in lab-- oratory equipment and the supportive cooperation of all fellow students, especially Rong-ching Hsieh and Fred W. Hall are fully appreciated. iii TABLE OF CONTENTS LIST OF TABLES . . . . . . . . . . . LIST OF FIGURES . . . . . . . . . . LIST OF APPENDICES . . . . . . . . . . LIST OF SYMBOLS . . . . . . . . 1. INTRODUCTION . . . . . . . . . 2. LITERATURE REVIEW . . . . . . . . . 2.1. Analytical Solution for Phase-Change Problems . . . . . . . . . . Numerical Methods Using Finite Differ- ences for Solving Phase-Change Problems The Implementation of Finite Element Models to Phase-Change Problems . Factors Influencing the Rate of Food Freezing . The Stability of the Finite Element Method in Comparison to the Analytical Solution and the Finite Difference Method . . . . . . . . . . . THEORETICAL CONSIDERATIONS . . . . . . 3.1. 3.2. 3.3. Temperature Dependent Physical Prop- erties . . . . . . . . . 3.1.1. Unfrozen Water Con ent . . . .2 Thermal Conductivity . . . . 3. Product Density . .4. Enthalpy and Apparent Specific Heat . . . . Governing Equations, Initial and Bound- ary Conditions . . . . . . Finite Element Formulation . . urge) wrap: 3.3.1. Development of the Model for Two— Dimensional Heat Transfer in Elliptical Geometry . . . . iv Page vi vii xi xii 12 15 23 25 25 26 27 28 28 3O 33 33 Page 3.2.2. Development of the Model for Axi- symmetrical Heat Transfer in Trapezoidal Geometry . . . . 42 3.4. Computer Implementation and Finite Element Grid . . . . . . . . . 45 4. EXPERIMENTAL . . . . . . . . . . . 55 4.1. Equipment . . . . . . . 55 4. 2. Material for Food Product . . . . . 63 4.3. Procedures . . . . . . . . . . 64 5. RESULTS AND DISCUSSION . . . . . . . . 67 5.1. Comparison of Numerical Simulation Model with Finite Difference Method During Food Freezing . . . 67 5.2. Verification of Computer Simulation with Experimental Data . . 71 5.2.1. Freezing of Product with Ellipti- cal Geometry . . . 71 5.2.2. Freezing of Product with Trape- zoidal Geometry . . . . . . 80 5.3. Influence of Area Average Enthalpy on Freezing Time . . . . . . . . . 87 5.4. Sensitivity Analysis . . . . . . . 90 5.4.1. Geometric Size . . . 90 5.4.2. Initial Product Temperature . . 95 5.4.3. Time Step . . . . 96 5.5. Application of the Numerical Model in the Food Freezing Process . . . . . 100 6. CONCLUSIONS . . . . . . . . . . . 107 7. RECOMMENDATIONS FOR FURTHER STUDY . . . . 109 APPENDICES . . . . . . . . . . . . . 110 BIBLIOGRAPHY . . . . . . . . . . . . 144 LIST OF TABLES Page The "Student t-test" for Experimental and Predicted Temperature in Codfish Freezing . 70 The Results of "Student t-test" to Evaluate Agreement Between Experimental and Pre- dicted Temperature at the Slowest Freezing Point Location . . . . . . . . . . 76 Experimental and Predicted Freezing Time for Various Shapes and Air Speeds . . . 89 Moisture Content and Density of Ground Beef at Various Experimental Treatments . . . 126 Physical Properties of Ground Beef . . . 127 Input Parameters for Computer Program . . 128 Recorded Temperature Data (°C) at One-Hour Time Increments During Freezing Process . 129 Mathematical Models Describing Surface Temperature as a Function of Time and Location at v = 7.9 m/s . . . . . . . 143 vi Figure 2. L») LIST OF FIGURES The relationship of thermal properties of food product and temperature during freez- ing according to Comini et a1. (1974) . The relationship of thermal properties of food product and temperature during freez- ing according to Tarnawski (1976) . The relationship of thermal properties of food product and temperature during freez- ing according to Heldman and Gorby (1974) Simplex two-dimensional element and area coordinates . . . . . . . . . . Simplex triangular axisymmetric element Flow diagram of Main Program . Flow diagram of subroutines SETMAT and TRANSIENT . . . . . . . . . Flow diagram of subroutines READ l and PROP 1 O O O O I O O O O 0 Flow diagram of subroutines HVAR and TVAR The two-dimensional elliptical finite ele- ment grid with simplex triangular elements . . . . . . . . . . The axisymmetrical trapezoidal finite ele- ment grid with simplex triangular elements . . . . . . . . . . Schematic diagram of the experimental set- ups for freezing process (top view) . . Elliptical geometric product with thermo- couple junctions and wire structure . . vii Page l7 17 36 43 49 50 51 52 53 56 57 Page Nodal locations of thermocouple junctions on the cross-sectional area in elliptical product . . . . . . . . . . . . 59 Trapezoidal geometric product with thermo- couple junctions and wire structure . . . 6O Nodal locations of thermocouple junctions on the cross-sectional area in trapezoidal product . . . . . . . . . . . . 61 The two-dimensional finite element grid for slab geometry with simplex triangular elements . .. . . . . . . . . . . 68 Temperature history of codfish fillets during freezing . . . . . . . . . 69 The time-temperature history of elliptical product during freezing at air speed 7.4 m/s . . . . . . . . . . . . 72 The time—temperature history of elliptical product during freezing at air speed 11.3 m/s . . . . . . . . . . . . 73 The time-temperature history of elliptical product during freezing at air speed 1502 m/S O O O O O O O O O I O I 74 The isothermal fields inside elliptical geometry product after 2.5 freezing hours at air speed 11.3 m/s . . . . . . . 77 The isothermal fields inside elliptical geometry product after 2.5 freezing hours at air speed 7.4 m/s . . . . . . . . 78 The isothermal fields inside elliptical geometry product after 2.5 hours freezing at air speed 15.2 m/s . . . . . . . 79 The time-temperature history of trapezoidal product during freezing at air speed 7.4 m/s . . . . . . . . . . . . 81 viii Figure Page The time-temperature history of trapezoi- dal product during freezing at air speed 11.3 m/s . . . . . . . . . . . 82 The time-temperature history of trapezoi- dal product during freezing at air speed 15.2 m/s . . . . . . . . . . . 83 The isothermal fields inside trapezoidal product after 1.0 freezing hour at air speed 11.3 m/s . . . . . . . . . 84 The isothermal fields inside trapezoidal product after 1.0 freezing hour at air speed 7.4 m/s . . . . . . . . . . 85 The isothermal fields inside trapezoidal product after 1.0 freezing hour at air speed 15.2 m/s . . . . . . . . . 86 The isothermal fields inside elliptical product at the time the freezing process terminated based on area average enthalpy (v = 15.2 m/s) . . . . . . . . . 91 The isothermal fields inside trapezoidal product at the time the freezing process terminated based on area average enthalpy (v = 15.2 m/s) . . . . . . . . . 92 The influence of geometrical size on the freezing rate of elliptical product . . 93 The influence of geometrical size on the freezing rate of trapezoidal product . . 94 The influence of initial temperature on the freezing rate of elliptical product . 96 The isothermal fields for various initial temperatures inthe elliptical product after 2.5 freezing hours . . . . . . 97 The influence of initial temperature on the freezing rate of trapezoidal product . 98 The isothermal fields for various initial temperatures bathe trapezoidal product after 1.0 freezing hour . . . . . . 99 ix Figure Page The effect of time step used in numerical scheme on the freezing rate of elliptical product . . . . . . . . . . . . 101 The isothermal fields for two different time steps in the elliptical product after 2.5 freezing hours . . . . . . . . 102 The effect of time step used in numerical scheme on the freezing rate of trapezoidal product . . . . . . . . . . . . 103 The isothermal fields for two different time steps in the trapezoidal product after 1.0 freezing hour . . . . . . 104 The air flow pattern around a circular obstacle . . . . . . . . . . . 136 The air flow pattern around an elliptical obstacle . . . . . . . . . . . 137 Local surface heat transfer coefficient for circular cylinder and elliptical cylinder . . . . . . . . . . . 138 The surface temperature on trapezoidal ground beef during air-blast freezing at v = 7.9 m/s . . . . . . . . . . 142 LIST OF APPENDICES Appendix Page A. Computer Program TWODFR and AXISFR . . . 111 B. Physical Properties of Ground Beef, Com- puter Input Parameters and Variables, and Experimental Data . . . . . . . . . 125 C. The Transformation of Local Surface Heat Transfer Coefficient from the Circular Cylinder to the Elliptical Cylinder . . . 130 D. The Data Fitting of Local Surface Tempera- ture at Trapezoidal Product Geometry Using Least Square Regression . . . . . . . 139 xi a(T) a, b, c [B] [C] CPA CPI CPS CPW [D] DI DS Dw EMS EMW {F} {9} LIST OF SYMBOLS Area of a triangular finite element, In2 Coefficient of temperature balance at tempera— ture T, m2/h Parameters defined in equation (3.34), a in m2 and b & c in m Gradient matrix defined in equation (3.39) Capacitance matrix defined in equation (3.41) Volumetric specific heat capacity, J/m3 K Apparent specific heat of product, J/kg K Specific heat of ice, J/kg K Specific heat of solid, J/kg K Specific heat of water, J/kg K Specific heat of product, J/kg K Matrix defined in equation (3.28) Density of ice, kg/m3 Density of solid, kg/m3 Density of water, kg/m3 Effective mass of solids, kg solids/kg product Effective mass of water, kg water/kg product Force vector defined in equation (3.41) Matrix defined in equation (3.27) xii :3" LT. 5| IK ICP IWC [K] KPI KPS KPW k(T) XX YY 22 t1 W WW LW Enthalpy, J/kg 2 Local surface heat transfer coefficient, W/m K Average surface heat transfer coefficient, W/m2 K Specific enthalpy: J/m3 Initial thermal conductivity, W/m K Initial specific heat of product, J/kg K Initial water content, kg water/kg product Stiffness matrix defined in equation (3.41) Thermal Thermal Thermal Thermal Thermal conductivity conductivity conductivity conductivity conductivity T, W/m K Thermal Thermal W/m K Thermal Thermal Thermal conductivity conductivity conductivity conductivity conductivity of ice, W/m K of solids, W/m K of water, W/m K of product, W/m K of product at temperature of continuous phase, W/m K of discontinuous phase, in x direction, W/m K in y direction, W/m K in z direction, W/m K Area coordinates of a triangular element defined in equation (3.36) Latent heat of fusion of solvent A, J/kg Latent heat of fusion of water, J/kg xiii :3) PD *3 DF *3 F3 t-J IF *3 IN Latent heat of fusion of product, J/kg Molecular weight of solvent A, kg/kg-mole Molecular weight of solvent B, kg/kg-mole Molecular weight of solids, kg/kg-mole Molecular weight of water, kg/kg-mole Volume fraction of discontinuous phase, dimensionless Effective mass of solvent A, kg/kg product Effective mass of solvent B, kg/kg product Mass fraction of water in the product, kg water/ kg product Shape function defined in equation (3.33) Number of nodal volume Product density, kg/m3 Parameter defined in equation (3.9) Heat generation inside the body Gas constant, 8314 J/kg-mole K Coordinates of triangular nodes for axisymmetric element, m Surface area, m2 Temperature, °C Depressed freezing point, K Freezing point temperature, K Initial freezing point, K Initial product temperature, K xiv At UFWC VFD VFI VFS VFW <| WC Product surface temperature, °C Temperature at the half thickness of the slab, °C Temperature at the first nodal point from the product surface, °C Ambient temperature, °C Time, second Time step, second Unfreezable water content, kg water/kg product Volume, m3 Volume fraction of discontinuous phase, dimensionless Volume fraction of ice, dimensionless Volume fraction of solids, dimensionless Volume fraction of water, dimensionless Air speed, m/s Specific volume, m3/kg Unfrozen water content, kg water/kg product Mole fraction of solvent A, dimensionless Mole fraction of water, dimensionless Coordinates of triangular nodes for two- dimensional element, m Space increment in finite difference scheme, m Density, kg/m3 Latent heat effect, J/m3 XV 1 Heat fraction of water, J/m3 w x Chi function Subscripts f Condition after freezing i, j, k Nodes of a triangular finite element p Condition before freezing u Define nodal point in finite difference scheme, x = u A x Superscripts m Define time level in finite difference scheme e A triangular finite element xvi l . INTRODUCTION Frozen food is one of the most important food products in the United States. The total pack and sales value of frozen foods in the United States reported by product category in 1966 was 6,284 million kilograms and 6,244 million dollars (Tressler et al., 1968) respectively. The design of food freezing processes is based primarily on refrigeration requirements for freezing and rate of product freezing. Appropriate methods to estimate the rate of freezing is important to achieve optimum design which results in good product quality and in an efficient process to avoid excessive energy consumption. Even though there have been significant improvements in the prediction of food freezing rates, the phenomena of phase- change, the influence of product thermal properties, the importance of product geometry and the effect of freezing environment on the freezing rate are not fully under- stood. The most recognized exact solutions to predict the freezing time in the food freezing process are Plank's equation and Newmann's solution (Bakal and Hayakawa, 1973) or the solution recently developed by Golovkin et a1. (1973). All of these are limited to special boundary 1 conditions, constant product thermal properties and to geometrically regular shapes, i.e., infinite slab, infinite cylinder, and sphere. Numerical solutions utilizing finite differences methods for temperature dependent thermal properties and regular geometric shapes have been discussed by Bonacina et a1. (1973), Charm et a1. (1972), Cleland and Earle (1977), Fleming (1973), Joshi and Tao (1974), Heldman (1974b), Heldman and Gorby (1974b), Hsieh et a1. (1977), Lescano and Heldman (1973), and Tarnawski (1976). The finite difference approxima- tions are acceptable, but the simulation method lacks flexibility to incorporate more complex geometry of food products and boundary conditions. Comini and Bonacina (1974), De Baerdemaeker et a1. (1977), Rebellato et a1. (1978) and Singh and Segerlind (1974) suggested Unaimplementation of finite element methods for estimating freezing time of food to reduce complex geometry and boundary condition problems. While the finite element analysis has proved excellent in accommodating linear and non-linear heat conduction in the freezing process, no investigations have verified the simulation results experimentally. It should be empha- sized that previous investigations have assumed that the product thermal properties are constant or a linear func- tion of temperature. Numerical techniques utilized by Heldman and Gorby (1974b), Hsieh et a1. (1977) and Lescano and Heldman (1973) to account for variation of product thermal properties during freezing have provided accurate predictions of freezing time. During the freezing process when air is used as refrigerant, convective surface heat transfer coefficients become an important factor in prediction of freezing rate. Most of the investigations conducted have assumed constant surface heat transfer coefficient as the boundary condition for a given food product. Katinas et a1. (1976) and Zdanavichyus et a1. (1977) published an experimental result suggesting that the convective heat transfer coef- ficient varies sinusoidally along surface of a cylinder. The primary objective of this research was to develop a numerical solution model using finite element methods to predict the rate of freezing of food products with anomalous shapes. Specific objectives were as follows: 1. To develop a computer program utilizing the finite element method to simulate the freez- ing process of elliptical and trapezoidal food products subjected to various boundary conditions. 2. To incorporate temperature dependent thermal properties during phase-change into the computer algorithm for both two-dimensional and axisymmetric heat transfer problems. To investigate the influence of various bound- ary conditions on the freezing time: con- stant surface heat transfer coefficient, heat transfer coefficient as a function of location on the surface of food product and variable surface temperature during freezing. To incorporate a method to estimate an optimum freezing time based on area average enthalpy in the finite element analysis and to compare this estimate to a conventional method based on the slowest freezing point location. To conduct experimental measurement of the freezing rate of elliptical and trapezoidal shape food product in the laboratory, using air-blast freezing method and ground beef for the freezing material, in order to verify the numerical solution. 2 . LITERATURE REVIEW 2.1. Analytical Solution for Phase-Change Problems Most analytical solutions used to estimate the freezing time in phase-change problems have been based on either solving heat balance equations (Plank's equation and Tanaka and Nishimoto's formula) or solving Fourier's equation of unsteady-state heat conduction (Newmann's solution, Tao's chart, and Tien's approach). All approaches have limitations of assuming constant product thermal properties and assuming regular geometrical shapes, i.e., infinite slab, infinite cylinder, and sphere (Bakal and Hayakawa, 1973; Carslaw and Jaeger, 1959). Bakal and Hayakawa (1973) indicated further that all the above methods used either a single temperature or a spe- cified range of temperatures, during phase-change. Slavin (1964) pointed out the inaccuracy of Plank's formula for calculation of freezing times for food, and Charm and Slavin (1962) reported 40 to 80 percent differences between Newmann's equation and experimental data in freezing time for cod fillets. Cho and Sunderland (1974) attempted to improve the exact solution by assuming thermal conductivity to vary linearly with temperature. The analysis applied to both melting and solidification of semi—infinite bodies but used the fusion temperature as a fixed temperature while phase-change occurred and did not account for the variability of product thermal properties other than thermal conductivity. likhailov (1976) developed an exact solution for freezing of a humid porous body, thus solving for moisture distribution as well as temperature distribu- tion. The analytical method is applied to Stefan-like problems which consider the occurrence of phase-change at a single temperature. Riley and Duck (1977) used the heat-balance integral method for the Stefan problems in freezing of a three-dimensional cuboid with all thermal properties of the product assumed constant. The authors also mentioned the unresolved question of accuracy even though some criteria have been established for semi- infinite region by Langford (1973). Golovkin et a1. (1973) suggested mathematical models for freezing of meat in two-sided slab, cylinder, and sphere geometry. The Stefan assumptions on the phase interface in the integral form were applied to obtain more accurate solution than if the differential form was used. Hayakawa and Bakal (1974) proposed formulas to predict transient temperatures in food during freezing and thawing. Phase changes are assumed to occur over a range of temperatures and the geometry of food is an infin- ite slab with insulation on one side. During freezing, the material is observed to move through an unfrozen state, partly frozen state and a frozen state. Freezing processes are divided into several periods: (a) precooling, (b)first phase-change, where a partly frozen zone moves along the direction of heat transfer, (c) intermediate phase-change, where a partly frozen zone exists throughout the body until the surface body temperature reaches the final freezing point, (d) second phase-change, where a frozen zone moves along the slab thickness, and (e) tem- pering when the body is completely frozen. The experiment conducted to verify the formula indicated that the mathe- matical model was in good agreement for all periods of the freezing process except intermediate phase—change. Difficulty was also encountered in determining the final freezing point during the intermediate phase-change. 2.2. Numerical Methods Using Finite Differ- ences for Solving Phase-Change Problems Many researchers have explored numerical techni- ques in order to get more accurate solutions for phase— change problems than obtained from analytical methods. Bonacina and Comini (1973a) developed a numerical solu- tion using an implicit finite difference scheme suggested by Lees (1966) which involves three time levels. 31:” k1; - T111 + I «21+: ‘23:"? + - é kit: - T1: 3 ‘23:” [ #1.; mix. - T? + T1; - T11> - k?_% (T? - T?_l + T?-1 - TT:i)] + c? T?—1 (2.1) The scheme was unconditionally stable and convergent and the applied boundary conditions were the first kind--pre- scribed surface temperature, and the fourth kind--variable surface temperature. Bonacina et a1. (1973) checked the numerical method against the analytical solution for the one-dimensional freezing problem (Luikov, 1968) and found the agreement to be within 3 percent. Using the same method for two-dimensional heat transfer, Bonacina and Comini (1973b) investigated the second and third kind of boundary conditions which were constant heat flux and linear heat transfer at the surface, respectively. How- ever, the analysis and the experiment to verify the method were conducted only for heating and cooling processes. Cleland and Earle (1977b) discussed the above numerical solutions thoroughly and suggested the use of the third kind of boundary condition for food freezing. _ _ 91 h (Too - TS) — k (dX for t > o (2.2) x=0 Instead of using the finite difference boundary condition as proposed by Bonacina et a1. (1973), the boundary con- dition was derived from a heat balance over the surface space increment which extended a distance of 0.5Ax from the surface. +— dT _ _ __§ T1? (T1 - Ts) — hm?S Tm) + C (Ts)dt (2.3) k l 2 The numerical scheme for one-dimensional heat conduction was proved to be in good agreement with experimental data in freezing of mashed potato and minced lean beef. Goodrich (1978) outlined a numerical procedure to solve one-dimensional phase-change problems with a defined moving boundary condition at a fixed temperature which was the product freezing point. The central difference scheme was utilized in the numerical technique and the thermal properties of product were considered to vary linearly with temperature. Hashemi and Sliepcevich (1967) presented a numeri- cal solution for one- and two-dimensional temperature distribution in an isotropic medium where phase-change occurred in finite temperature intervals. The procedure utilized predictor-cOrrector and implicit finite differ- ence methods in solving Neumann's solution for one- dimensional heat transfer and incorporated the 10 alternating direction method with the predictor-corrector formula for two-dimensional heat transfer. Shamsundar and Sparrow (1975) employed an enthalpy model to analyze multidimensional conduction phase-change g? pI dV = Jr k grad T - 8 dA (2.4) V A The model was approximated by using the implicit finite difference method, but no experimental work was conducted to confirm the simulation. Joshi and Tao (1974) utilized the finite differ- ence method to solve the problem of axisymmetrical freez- ing of food products. The method implemented the first and third kind of boundary conditions and used forward- difference for the time derivative and central-difference for the space derivative. The product thermal properties varied with temperature and fraction of frozen water except for product density which was assumed constant. The verification of the numerical method in rectangular beef freezing experiments gave satisfactory results. Tarnawski (1976) proposed a mathematical model to solve one-dimensional heat and mass transfer during food freezing using the third kind of boundary condition. The mathematical model took into account the discontinuity and nonlinearity of product thermal properties and was ll approximated by finite difference methods. The simulation results did not use the mass transfer potential as described in the model, and were not verified by experi- mental data. Lescano and Heldman (1973) developed a mathemati- cal model to predict the thermal properties of a food product based on variable composition of water and ice within the product as temperature changed during freezing. A numerical scheme was later outlined to solve the one- dimensional symmetric heat transfer problem with a bound- ary condition of the third kind using the Crank-Nicholson formula for finite difference analysis. The application of the computer simulation yielded good agreement with experimental data in freezing of slab codfish. Heldman and Gorby (1974a).improved the prediction model for vari- able product thermal properties by implementing the Kopel- man equation (1966) to describe the relationship of thermal conductivity with product composition which was changing as temperature decreased during freezing. The improved mathematical model along with finite difference methods were utilized successfully to solve one- dimensional transient heat transfer in ice cream freezing. Numerical solutions, using finite difference method, to simulate the freezing process for spherical geometric food products, were developed incorporating the above 12 prediction model for variable product thermal properties (Heldman and Gorby, 1975). The finite difference equa- tions were derived by both forward and pure implicit methods and applied to IQF (Individual Quick Freezing) of cherries with acceptable results. Hsieh et a1. (1977) modified the above computer simulation techniques to pre- dict the freezing times and temperature history for dif— ferent fruits and vegetables. 2.3. The Implementation of Finite Element Models to Phase-Change Problems The application of finite element methods in solv- ing heat conduction problems has been discussed by Zienkiewicz and Cheung (1965), Visser (1965), Wilson and Nickell (1966), and Richardson and Shum (1969). The analysis has included problems with steady and unsteady state heat transfer, linear and nonlinear boundary condi— tions and nonstationary temperature distribution. More examples of the finite element models applied to transient heat conduction can be found in references such as Zienkiewicz (1971), Desai and Abel (1972), and Seger- lind (1976). Emery and Carson (1971), and Bruch and Zyvoloski (1974) discussed the accuracy and efficiency of the finite element method and illustrated acceptable com— parison to exact solution and finite difference method for both linear and nonlinear two-dimensional heat l3 conduction problems. Comini and Lewis (1976) developed a numerical solution using finite element methods for two— dimensional and axisymmetrical problems involving heat and mass transfer in porous media. The simulation was in agreement with analytical solution in drying of a geometri- cally slab material. Singh and Segerlind (1974) applied the finite element models to describe time-dependent axisymmetric problems in heating of a cylindrical food can containing homogeneous material and to simulate heat- ing of a chicken leg composed of four different materials. De Baerdemaeker et a1. (1977) discussed the application of finite element analysis to pear cooling and to rec- tangular beef steak frying. Bonnerot and Jamet (1974) introduced the imple- mentation of the finite element method for the one- dimensional Stefan problem to determine the position of the free boundary of phase-change. The quadrilateral elements were used and the temperatures were calculated for all element nodes as well as additional nodes along the moving boundary at each time step. The expanding grid used to track the free boundary is only useful for small boundary motions and for cases in which the tem- peratures on one side of the boundary are always zero (Wellford Jr. and Ayer, 1977). The latter authors pro- posed a fixed grid of standard space-time finite elements l4 and discontinuous interpolation to define the finite ele- ment model on the special elements which were the quadri- lateral elements crossed by the free boundary at any particular step in time. Krutz (1976) developed a finite element computer model for phase-change from solid to liquid in determining the time-temperature history of a welded joint. Thermal conductivity and specific heat were considered as a linear function of temperature and the model included radiation, convection and heat flux on the surface. Comini et a1. (1974) applied the finite element method to freezing analysis with nonlinear boundary con- ditions. The physical properties were considered to vary linearly with temperature in addition to a jump within a small temperature interval (2AT) at the freezing point. The nonlinear boundary condition took into account imposed heat flux and rates of heat flow per unit area due to convection and radiation on the surface. Simple triangular elements were used and the three level scheme suggested by Lees (1966), as discussed previously in Sec- tion 2.1, was introduced for time-stepping instead of Crank-Nicholson algorithm. The freezing simulation pro- gram was used to predict the position of frozen boundary in slab form and for soil freezing. Comini and Bonacina (1974) presented the appli— cation of the above method in food freezing to overcome 15 the lack of flexibility of finite difference method in solving the irregular geometrically shape problems. The thermal properties were calculated based on the decreas- ing mass fraction of water during freezing but the method to detect the fraction of water during the freezing process was not given. The latter method will be dis- cussed in the next section. Bonacina et a1. (1974) com- pared results of the above finite element method with experimental data in freezing of Tylose samples which have been modeled after lean beef. The heat conduction problem was selected to be one-dimensional with boundary condition of the first kind. The error between measured and calculated temperature at the center and surface of the slab was found to be less than 2 percent. Rebellato et al. (1978) used this simulation program to solve the two-dimensional heat conduction problem utilizing a second-order quadrilateral element grid in estimating the freezing rate of lamb carcass and beef side. However, no experimental data has been reported so far to verify the results of this two-dimensional irregular geometry food freezing problem. 2.4. Factors Influencing the Rate of Food Freezing The rate of food freezing is influenced by sev- eral factors including temperature of surrounding medium, 16 size and shape of frozen product, thermal properties of the product and surface heat transfer coefficient. It has been widely known that the lower the surrounding tem- perature, the higher the rate of freezing. Tarnawski (1976) and Hsieh et a1. (1977) showed that the relation- ship between freezing time and surrounding temperature was nonlinear. Size and shape have always been important factors to consider in the analysis of freezing as pre- viously discussed for the analytical solution, the finite difference method and the finite element method. Hsieh et a1. (1977) investigated the influence of product diameter on freezing time for various fruits and vege- tables. The freezing time increased linearly as the product diameter became larger. The determination of the transient temperature field and the rate of freezing for food products using the assumption that phase-change exists at a constant temperature is not accurate. Several investigators have proposed formulas to estimate product thermal properties as freezing occurs over a temperature range or during the whole process. Comini et a1. (1974) suggested that phase- change could be assumed to occur at a temperature range from T1 (initial temperature) ;-1°C until a certain value of Tf (final temperature) where there was Tp (peak N temperature) = - 3°C in between (Figure 2.1.). The 17 k m o o 1: oc J? T Figure 2.1. The relationship of thermal properties of food product and temperature during freezing according to Comini et a1. (1974). k,c:, a T fr T Figure 2.2. The relationship of thermal properties of food product and temperature during freez- ing according to Tarnawski (1976). 18 final temperature is estimated as the value that gives the best fit between calculated results and experimental data. The formulas to calculate thermal conductivity and specific heat capacity above and below freezing are given as follows cp = mc-CPW + (1 - mc)-CPS (2.5) kp = mc-{W + (l - mc)-KS (2.6) cf = mC°CPI + (l - mc)-CPS (2.7) kf = mcoKI + (l - mc)-KS (2.8) The value of heat capacity at Tp can be obtained from evaluating the latent heat effect which is the area of heat capacity versus temperature as Ti’ Tp, and T are f known A = mC°A (2.9) The disadvantages of this method are that phase-change is assumed to occur only over a short temperature range, the final temperature has to be chosen arbitrarily if experi- mental data do not exist and the relationship between thermal properties and temperature is actually nonlinear. Tarnawski (1976) presented nonlinear function to describe the relationship between physical parameters of 19 a food product and temperature where the function was dis— continuous and nondifferential at the freezing point (Figure 2.2). The model to compute thermal conductivity and specific heat is not published and the calculation for beef are given as follows for 248 K i T i TF a(T) [-11032.6509 - 261.42196 (Z) - 18252.4705 (2)2 - 26.133707 (2)3 - 133.87762 (z)4 - 7.31961337 (2)5 9 - 0.11645322 (2)6] x 10' mz/h (2.10) k(T) [~124634.1825 - 744465.0959 (Z) — 160251.468 (2)2 - 17695.612 (2)3 - 102.877816 (2)4 - 29.874884 (2)5 - 0.343013138 (2)6] x 10‘6 W/m K (2.11) for Tf 3 T 3 303.16 K a(T) = 0.00042 - 0.000001 (2) m2/h (2.12) k(T) = 0.476079324 - 0.0004026324 (Z) W/m K (2.13) where Z = T - 273.16 Since changes in product thermal prOperties dur- ing freezing are due to continuous depression of freezing 20 point and thus continuous changes in unfrozen water con- tent, the best approach is the method proposed by Held- man (1974a), and Heldman and Gorby (1974b). The unfrozen water content can be detected at any given time assuming food product as a mixture consists of water (solvent A), and ice together with food solids (solute B) M X _ A A mB mA 7 MB (1 - xA) (2'14) where XA = exp [(LA-MA/R1)(1/TIF - l/TDF)] (2.15) Thermal conductivity is obtained from the Kopelman equa- tion (1967) k = kc (l - Q)/[l - Q (l - M)] (2.16) _ 2 _ Q — M (l kd/kc) (2.17) Enthalpy and specific heat are computed from equations (Lescano and Heldman, 1973) H = EMS-CPS (T + 40) + wc-L + wc-CPW (T + 40) + MI-CPI (T + 40) - UFWC°L (2.18) c = AH/AT (2.19) P 21 The implementation of the above equations (2.14) - (2.19) are further discussed in the next chapter. Figure 2.3 illustrates the nonlinear relationship between thermal properties of food product with temperature as calculated using equations (2.14) - (2.19). The influence of surface heat transfer coefficient on freezing time has been investigated by several research— ers. Heldman (1974b) compared the surface heat transfer coefficient versus freezing time curves for lean beef obtained from various analysis, Charm (1971), Lescano and Heldman (1973), analysis graphical method, and modified Planck's equation (1958). The results indicated that the freezing time decreased significantly as the heat transfer coefficient increased to 25 W/m2 K. The same result was confirmed by Tarnawski (1976) for beef freezing. Hsieh et a1. (1977) found that the freezing time could be reduced significantly as the surface heat transfer increased to 40 W/m2 K for freezing of various fruits and vegetables. All the previous investigations were carried out for uniform surface heat transfer coefficient. The influ- ence of variable local heat transfer coefficient on the surface of food product during freezing has not been published. Katinas et a1. (1976) and Zdanavichyus et a1. (1977) presented the local heat transfer coefficient as a 22 kr pr Cp Figure 2.3. The relationship of thermal properties of food product and temperature during freez- ing according to Heldman and Gorby (1974). 23 function of location on the surface of a cylinder. The function behaved sinusoidally because the degree of tur- bulency of an inflowing air stream around the circular cylinder. 2.5. The Stability of the Finite Element Method in Comparison to the Analyti- cal Solution and the Finite Differ- ence Method Emery and Carson (1971) evaluated the use of the finite element method in the computation of temperature for linear and nonlinear two—dimensional problems and compared it with the finite difference method. The finite element method applied was utilizing linear, quad- ratic, cubic and special cubic elements. For the finite difference method, three different kinds of time step- ping schemes were analyzed, i.e., explicit, Crank- Nicholson, and Lex-Wendroff (1967). The authors concluded that the finite difference method required less core memory in the computer and gave faster execution time especially for variable thermal properties problems which needed computation at each time step. The finite ele- ment method had the advantages of solving heat conduction problems for arbitrary geometry and was more accurate. Furthermore, there were advantages associated with the ease of inputting the required data and the capability of altering the basic accuracy of the method. 24 Bruch and Zyvoloski (1974) discussed the implica- tion of the finite element method to solve linear and non- linear two-dimensional heat conduction problems. Rectangular prisms inaaspace-time domain were used as the finite elements and the implicit method was applied for the time-stepping scheme. The finite element solutions compared favorably with the results from analytical solu- tions and finite difference methods. It was found to be stable and convergent to the exact solution. Yalamanchili and Chu (1973) analyzed the stability and oscillation characteristics of the finite difference method and the finite element method with and without use of Galerkin's method of weighted residuals. The stability criteria were established by utilizing the‘general sta- bility, von Neumann formulas, and Dusinberre concepts (1961). For transient two-dimensional heat conduction in solids, the results showed that the region of favorable stability and oscillation characteristics were found to be significantly larger for the finite element method than for the finite difference method. The use of Galerkin method improved the degree of stability and reduced the oscillation. 3. THEORETICAL CONSIDERATIONS 3.1. Temperature Dependent Physical Properties The change in thermal prOperties of a food product during freezing is due to continuous freezing point depres- sion caused by a reduction in unfrozen water content. The development of mathematical equations to predict the thermal properties of food based on the freezing point depression has been described by Heldman (1974a) and Heldman and Gorby (1974a). This method has been success- fully utilized to predict the thermal properties of food product during freezing (Heldman, 1974b; Heldman and Gorby, 1974b; Heldman and Gorby, 1975; Hsieh et al., 1977). Assumptions regarding the temperature dependent physical properties are as follows: 1. The food product is homogeneous and iso- tropic. 2. The thermal properties are constant above the initial freezing point. 3. The food product consists of solids, water, and ice during the freezing process. While the thermal properties of food product vary nonlinearly according to temperature, thermal properties of product solids remain constant. 25 26 4. Below the initial freezing point, the food product is assumed to be an ideal binary system with continuous and discontinuous phases. Since the frozen food consist of three phases, it is reduced to one binary system during the first step and to another system during the second step. 3.1.1. Unfrozen Water Content The relationship between the mole fraction of solvent and the freezing point depression in a solution is described in equations (2.14) and (2.15). Heldman (1974a) proposed that the unfrozen water content at a given temperature during food freezing can be predicted using the above equations assuming the liquid water is the solvent and solids is the solute. The molecular weight of product solids above freezing point is calcu- lated assuming product solids as the solute and water as solvent. Thus, equations (2.14) and (2.15) become xw = exp [(LW-MW/R)(1/TIN - l/TIF)] (3.1) _ EMS-XW°MW MS ’ BMW (1 - XW) (3'2) where EMS = 1 - Iwc (3.3) 27 BMW in equation (3.2) is modified taking into account the small amount of unfreezable water content at low tempera- ture. This approach has been thoroughly discussed by Lescano and Heldman (1973). The effective mass of water becomes BMW = IWC - UFWC (3.4) The total unfrozen water content at any given temperature below the freezing point can be computed as follows XW = exp [(Lw-MW/Rlxl/TIF - l/TDP)] (3.5) _ BMs-xw-Mw EMW — Ms (1 _ xw) (3.6) WC = BMW + UFWC (3.7) 3.1.2. Thermal Conductivity Kopelman (1967) derived mathematical models to predict thermal conductivity in food products for both isotropic and anisotropic systems. The model for an isotropic system in a two component product is _ 1 - Q k _ kc (l _ Q<1 _ M)) (3.8) M2 (1 - kd/kc) (3.9) 10 ll 28 The first step in the use of the Kopelman equation has water considered as the continuous phase for the water-ice system. Then, the water-ice mixture is treated as a continuous phase while the product solids is taken as discontinuous phase for the second use of the equation. 3.1.3. Product Density The product density decreases according to the prOportional changes of the mixture and can be expressed in the following manner (Heldman and Gorby, 1974a) V = l/PD = WC/DW + MI/DI + EMS/DS (3.10) For the computer program, the density of solids before freezing must be obtained by substituting initial product density, initial water content and MI = 0 into equation (3.10) DS = l/(l/IPD - WC/DW) (3.11) Then, the product density at any given time is solved by using equation (3.10) and unfrozen water content as com- puted according to Section 3.1.1. 3.1.4. Enthalpy and Apparent Specific Heat The enthalpy of the food product can be obtained based on the specific heat of product components and the 29 unfreezable water content. Utilizing a reference tempera- ture of -40°C, Lescano and Heldman (1973) expressed the enthalpy as H = BMs-CPs-(T + 40) + wc-LW + WC°CPW-(T + 40) + MI'CPI'(T + 40) - UFWC°LW (3.12) The multiplication of water content by its latent heat of fusion, WC-LW, accounts for the heat released during phase-change inside the food product. The specific heat of solids can be determined by solving the equation ICP = IWC'CPW + MS°CPS (3.13) while CPI is obtained from Dickerson equation (Dickerson, 1969) CPI = A + B°T (3.14) where A = 1.9507941 ED ll 0.00206153 for T as absolute temperature WC in equation (3.12) is calculated by solving equations (3.5) - (3.7), while MI is obtained as follows MI = 1 - WC - MS (3.15) Assuming the relationship between enthalpy with temperature is a continuous function, the apparent specific 30 heat of a food product can be expressed as differential change in enthalpy CPA = AH/AT (3.16) For this study, AT = 0.003°C was used in the numerical solution. 3.2. Governing Equations, Initial and Boundary7Conditions In food freezing, heat transfer occurs primarily by conduction. This research dealt with two—dimensional heat transfer in eliptical and trapezoidal shapes. The governing differential equation for heat conduction in isotropic bodies is known as the Fourier heat conduction equation (Carslaw and Jaeger, 1959). For two-dimensional heat transfer, the equation is as follows 1-3 8 = 21[ 8T] + 4: [k :2] + Q' (3.17) 9p08t: 8x xx 8x 8y yy 8y The body is at uniform temperature initially T = T0 at t = 0 (3.18) The boundary conditions are 3T 8T _ kxx[F§] + kyy[§§] + h (T - Tm) — 0 (3.19) at the convective surface and for t > 0 and 3) 1% o; ‘< 31 = 0 at the insulated surface and for y=0 any given time (3.20) Further assumptions regarding the modes of heat transfer are listed below: 1. Heat transfer from the freezing medium (air) to the product occurs by convection and heat moves by conduction within the product. Energy transports occur only in x and y direction. The rate of heat transfer within the food is uniform along the x and y direction. Thus, kxx = kyy’ The surface heat transfer coefficient is a function of location along the surface; however, it remains constant at a given posi- tion during the freezing process. The surrounding temperature and velocity of freezing medium are constant and uniform. Water vapor transport from the product to the air is negligible. Thus, mass transfer within the product and on the product surface are neglected. 32 The governing differential equation, initial and boundary conditions for axisymmetrical heat transfer in trapezoidal bodies are expressed in equations (3.21) - (3.24) Cppg—E— = 11:- krrgg + 831? [kn-‘3"; + 3% [kzzg—E] + Q' (3.21) Initial condition T = T0 at t = 0 (3.22) Boundary conditions krr[-:—:] + kzz [35%] + h (T - Too) = o (3.23) along the convective surface and for t > 0 and g; = 0.at the insulated surface for any r=0 given time (3.24) The assumption made for axisymmetric heattransfer are nearly the same as for two-dimensional, except for surface heat transfer coefficient and surface temperature: 1. Heat transfer from the freezing medium to the product occurs by convection and by con- duction within the product. 2. Energy transports occur only in r and 2 directions. 33 3. The rate of heat transfer within the food is uniform along the r and 2 directions. Thus, k = k . rr 22 4. The surface temperature is a function of time and location along the surface. 5. Velocity of freezing medium is stable and uniform. 6. Mass transfer within the product and on the product surface are neglected. 3.3. Finite Element Formulation 3.3.1. Development of the Model for Two-Dimensional Heat Transfer in Ellipti- cal Geometry In the finite element method for field problems such as heat conduction, the integral of a function is minimized using the calculus of variations (Segerlind, 1976). The governing equation for two-dimensional heat transfer (3.17) and its boundary conditions (3.19) and (3.20) can be formulated as follows _1 312 312 . 3T 'X - J[; [kXX(8X + kyy(8y) 2Q T + 2cppT 5E) dv + {—121 (T-Toc)2ds (3-25) where T0° is ambient temperature while V denotes the total volume of the body and s is surface area. 34 The advantage of calculation of temperature depen- dent thermal properties using freezing point depression as described in 3.1. is that it has taken change of phase into account. Specific heat of product is derived as differential change in enthalpy which is a function of latent heat, thus the specific heat is also a function of latent heat. Since latent heat has been incorporated into thermal properties of product (k, cp, and p), the term for heat generation inside the body, Q', can be eliminated. _ _1_ 312 M2 31-: X - .1; 2 [kXX(ax + kyy (8y + 2cppT at] dV +- Jr g-(T - Tw)2 ds (3.26) S Defining two matrices {911‘ = {—33 93'- 8X 3y (3.27) k 0 and [0] = XX (3.28) 0 k yy Equation (3.26) can be rewritten as X ll NI l-' ( {g} T [D] {g})dv + jf CppT 33 dV v + [g [T2 — 261200 + T3,) «is (3.29) S 35 Since the function for T is defined over individual sub- (e) regions called elements, T , equation (3.29) can be transformed to a sum of the integrals over the total num— ber of elements, E. B X = Z jr %({g(e)}T [D(e)] {9(9)})dv =1 V(e) + c 6—3-2 dV + 2(T(e)T(e) - 2T(e)T V e S(e) + Ti) ds B or x = z x(e) (3.30) e=l To minimize the function which is equating the summation (e) of derivatives of x with respect to T with zero, it is necessary to express the equation in terms of nodal values of temperature {T}. Utilizing two—dimensional simplex elements (Figure 3.1.a), then T(e) = [N(e)] {T} (3.31) (e) or T N.T. + N.T. + N T (3.32) 11 j] k K where N8 is a shape function, and i, j, k are denoting nodes of each triangular element whose equations are given by 36 .moumcflpuooo comm paw acoEmHm Hmcoflmcoeflpuozu onQEHm .H.m muooflm 1.x :me x "1 h x x a H H A» 5 x 31.x: d—-—— _ _ . _ _ _ _ _ 37 _ 1 Ni — 2A (a + bix + ciy) l N. = —— a. + b.x + c. 3.33 3 2A” 3 3y) ( ) N = ;L (a + b x + c y) k 2A k k k a1 = xjyk - xkyj bl = y3 - yk C1 = xk - X3 aj = Xkyi - x yk bj = yk - yl C3 = xi - xk (3.34) ak = xlyj - ijl bk = y1 - y] Ck = X] - X1 F— —H 1 X1 Y1 l A = — 1 x. . 3.35 2 3 Y3 < ) l xk yk L _ Introducing the area coordinates L L L (Figure 3.1.b) 1' 2' 3 as values indicating the area, equation (3.33) can be expressed as follows _ _ 1 L1 — Ni — 2A (ai + bix + ciy) L=N=—l-(a+bx+cy) (336) 2 3 2A3 3' j ' L = N = 1L (a + b x + c y) 3 k 2A k k k Substituting equations (3.31) and (3.33) into equations (3.27) and (3.30) 33(9) 7311“?) 80(9) BN (9’ 8X 1 —l ---2(- T1 { (e)} 8x 3x 8x 9 — = T. . Em) 8Ni(e) 631%) 3Nk(e) 3 (3.37) By a? a, 3;— T). Or {g(e)} = [3(9)] {T} (3.38) X(E) = j %{T}T [B(e)]T[D(e)] [B(e)]{T} dV V(e) +fcp0[N(e)] {T} [N(e)] dV 2.191} V(e) +f aim}T [N(e)]T [N(e)]{T} ds (e) S —f h "I. [N(e)]{T}ds (e) S +/ 3T: ds (3.38) (e) S where [B(e)] is called the gradient matrix and is defined as bl bj bk (e) _ l [B 1 - 2A (3.39) f1 03' Ck- The minimization of x becomes 39 E (e) a _ 3x _ T 3%- _ 87—10—— fua] [D] [B]{T} dV e-l v T a(T} +jf cpo [N] [N] dV at v T + h [N] [N]{T} ds h Too [N]T = 0 (3.40) J\_‘F‘m\“‘fi Equation (3.40), in compact form, can be formulated as follows [C] {T} + [K] {T} = {F} (3.41) where [C] = 2 p [N] dV, capacitance matrix, fvcp [K] = z. :03 8.”.ommm5m :92 - 18.2 6.... 658.00.. 2.03:2? 8 :0 oz... Lo. 3...; - 3...: 2632...... 9:2. 33:030 .n H H) a eve: :02. he Ede—3:0 0:: .3355 5205.09 . o 0:: .o 3. you 3 39.: 0:330:95 Sou .~ m03~e> >umv:=cn me Auk. easy-309:0» vonauumvud 0.3 .31: :oquosvo contented on... vcauuomo: 0:32. .93... :0 330 3.0330 9.0.5 0:7... 1 A .uu Ba x ..n .«n oueuaoaeu 2.. _ "cognac... vouch _, cc: cc:c.....c.,_c: .mmucuuuun acofiouo Dana-.0460 A»: + N: 925: cyan acoleuu £085: H ope: 500.5: 0:330 g — can» :o 320 ”20530 acocquuon L 88... 3:9..ch no .395: on» :o 3 on 3:23.. no ~ 03...... 05 co 2 8 — H nouuuoaoku :onauuunnzn “6:35.. £095.: etc: £083: €060; cuaaxoen an 95:50: «Eu woodman... To: Us... .oocuoeo: ESE; 6:9. 05 ocfiuonrucoou .3 Aoc + a: H H Lo .5330: 5330 cu gum :eu :ozeuoua «0 3:03.30 .0 ~03: 05 :o a 8 1.09.5: 05 :0 ca)— 51 READ Input initial values R, LW, Dw, DI, TO, KI, CPW, IWC, TF, UFWC, ICP, IK, MW CALL SUBROUTINE PARAM 1 Calculate the initial parameters necessary DS, KS, MS, CPS DS is solved from equation (3.11) KS is solved from equations (3.8) and (3.9) using water-solid system KS = [KW (VFSZ/3-Q)]/VFSZ/3 Q = (Kw - IK)/m3xm 0:9 E0 m3 EU 03 EU Ln .x.m musm3m EU EU EU EU EU EU m.m o.m m.h 4 . EXPERIMENTAL 4.1. Equipment An air-blast wind tunnel located in a low— temperature room was used to freeze the food product as illustrated in Figure 4.1. The length of the wind tunnel was 3.8 m with circular cross section of 46 cm I.D. The fan drew air into the tunnel using a three-phase electric motor of 3.73 KW; 440/220 Volt. The air speed was con- trolled by a circular opening with baffles which were regu- lated by a lever and positioned in front of the fan. The internal size of the freezing room was 6.7 m x 1.8 m x 2.4 m. Two evaporators provided the refrigeration effect to reduce the temperatureixlthe room to as low as -34.4°C with a deviation of :1.8°C when the fan was not operating. The food product, supported by a styrofoam plate was placed at the Center of the wind tunnel. The styro- foam support plate had dimensions of 46 cm x 35 cm x 5 cm, and functioned as an insulator to avoid heat losses through the bottom and edge of the food products (Fig- ure 4.2). The boundary conditions in equations (3.17) and (3.21) were satisfied by the design of the plate. A sharp leading edge was designed on the side of the 55 56 .ABwH> QOpv mmmooum mcflummum How masumm HmucwEfluwmxw may mo Ewummflo oflumesom .H.v musmflm 1"‘t Boom mcflumwum Hommoq mumo mxsam . _ H _ u “ i|.L j cumuumm 30am “fin mmnfiz Al/ umcmucmflmuum wamsoofisumna cam :oflaosm but HI 3 f/ VI! Ma. - . A w a n P n — _ Ta ”HM" . . . . .w" " mumllm..a 3"... 1+ .3" ..w mmw , n ..... “+5; m.” .M—\. m V Hos/com. . poswvoum 1 J\ on“: // p on um>mA// mumam Ecomonmum mcfitmasmom l\ 57 .wuduosuum muflB can mcofluocom mHQSOOOEuwcu saws uoscoum Ufluuweowm Hmuflumflaam .N.v wusmflm 58 styrofoam support plate facing the air flow to reduce turbulence created by air flow over the support plate and product. Unsheathed fine copper-constantan thermocouple wires (beaded by Omega Engineering, Inc.) were utilized to sense the temperature of the product. The thermo— couple had precision of i 0.417°C in the temperature range from -59.444°C to 93.333°C, and the diameter of each wire was 0.0125 cm. The wires were sheathed by teflon tubes with 0.055 cm I.D. to avoid direct contact with product. Copper-constantan thermocouple extension wires (gauge 20; 0.08 cm diameter) with polyvinyl insulation were used to connect the sensing thermocouple wire to the temperature recording instrumentation. The precision of the extension wires was the same as the sensing wires. The connectors were heavy-duty copper constantan miniature thermocouple connectors designed with fine gauge thermo- couple wires. Two different types of product shape were investi- gated in this study, a two-dimensional elliptical geome— try (Figure 4.2) and an axisymmetric trapezoidal geometry (Figure 4.4). The nodal locations measured by the thermo- couple junctions are illustrated in Figure 4.3 for the elliptical shape and Figure 4.5 for the trapezoidal shape. The teflon sheathed thermocouple wires were partially UN IT 'N CH Food Thermocouple 23 junctions Styrofoam plate I—ak—I—ak—Hl L2 k—z+2—k—3—*t.e+i—K 4 (a) Overall locations C’ (b) Locations in cross-sections AA, BB and CC (Figure 4.2) Figure 4.3. Nodal locations of thermocouple junctions on the cross-sectional area in elliptical product. 60 .ousuoouum oufl3 pom wcofluocsn wamsooosuosu :ufl3 uCDuCLQ oHuquomm ampMONwmmuB . v . v 6.8on 61 .uoopoum Hmcflonmmuu cw moum Hmcofluomm Immouo on» :0 mcoHuocsn mHQDOOOEHwnu mo mcofiumooH Hmpoz .m.v wusmflm A¢.v musmflmc mm can << wcofluommlmmouo ca mcoHumooq ADV wcoflumooH Hamum>o Amv oumHm Emowou>um rm~.~.*m~.~.+m~.~.+m~.~+l mlrl ml. poopoum 000% vm mcofluocsm mamsooosumne E z. :2: quapzsq 62 imbedded in the styrofoam plate and partially supported by a 0.05 cm diameter copper wire structure inside the product. Since the teflon sheathed wires and the c0pper wire structure had different thermal properties from the product, there was concern that the 24 nodal measurements located in one cross-section might influence heat conduc- tion within the product and produce significant error. To eliminate this possibility, the measured locations were installed in three different cross-sections (Figures 4.3 and 4.5) assuming that heat transfer occurred uniformly along the length of the elliptical body and along the circumference of the trapezoidal body. Three other thermocouple junctions were positioned to detect the air temperature at three locations near the product. The temperature measurements were recorded by the Fluke Model 2240A Data Logger with scanning speed capacity of 2.5 channels per second or about 11 seconds for one cycle containing 27 nodal readings. The output was printed in degrees centigrade for each nodal location and for each time step. The system accuracy was i 0.5°C in the temperature range of -l30°C to 0°C, and i 0.4°C in the temperature range from 0°C to 400°C when copper- constantan thermocouple wires were employed. 63 4.2. Material for Food Product Ground beef was used as the food product in the freezing process experiments for two considerations: (1) it was easily molded and shaped into the desired geometry and (2) it responded to the assumption of product homogeneity since its fibers and tissues had been broken during the grinding process. Ground beef was purchased as commercial lean ground beef from a local grocery store. Thermal conductivity values and initial freezing point of ground beef were obtained from literature refer- ences. Density of ground beef was determined by weighing the product before and after freezing. By calculating the volume of the product shapes and dividing the weight by the volume, the density of unfrozen and frozen meat was established (Table B.l). Specific heat of ground beef was calculated using Dickerson's formula (1965) and the moisture content of meat (MC) cp= 4.185 (0.4 + 0.006 MC) (4.1) where<fi9is .in J/g K and MC is .in percent. Charm (1971) suggested more elaborate equations to compute the specific heat of a food product Cp =4.184 (0.5 X + 0.3 XS + 1.0 Xm) (4.2) f 64 and c: = 4.184 (0.34 X + 0.37 X + 0.4 X P c p f + 0.2 x + 1.0 X) (4.3) a m where X was the weight fraction and subscripts a, c, p, f, s, and m indicated ash, carbohydrate, protein, fat, solid and moisture content of food, respectively. However, the comparison of the results from equations (4.1), (4.2), and (4.3) for specific heat of beef product illustrated only i 2 percent deviation (Heldman, 1975). In this investigation, equation (4.1) was used after the moisture content of lean ground beef was deter- mined experimentally according to AOAC procedure (AOAC, 1975). Three samples of ground beef, each 2 g, from each freezing treatments, were placed in aluminum dishes and dried for 4 hours in a Precision Scientific Model 625-A oven at 125°C temperature. The dishes were covered with lids to avoid contact with moist air when removed from the oven and cooled inside a desiccator prior to being weighed. The moisture content of samples for each treat- ment are listed in Table 3.1, and other physical proper- ties of ground beef are presented in Table 8.2. 4.3. Procedures The freezing experiments were conducted for three magnitudes of air velocity for both elliptical and 65 trapezoidal shapes. Each experiment was repeated three times. Air velocities at the center of the wind tunnel had been monitored by Chavarria (1978) utilizing a micro- manometer and a Pitot tube. The air speeds used in this investigation are presented in Table B.1 along with the air temperatures. The air temperature fluctuated during the freezing process due to the effect of turbulence within the closed circuit air pattern inside the freezing room and motor heat dissipation. The magnitude of fluctua- tion varied from : 0.8°C to i 2.0°C for an air temperature of -20.0°C (Table 3.1). The data from the experiments were needed to verify the results of computer simulation program described in Chapter 3. The operation of each experiment included the following steps: 1. Setting the freezing room temperature at -20°C at least 24 hours before the experi- ment. 2. Tempering the ground beef to a uniform ini— tial temperature of 18°C to satisfy the initial condition for the governing heat transfer equations. 3. Molding the ground beef into either ellip- tical or trapezoidal shape on the styrofoam support plate. 66 Programing the Fluke Data Logger and supplying the input parameters, number of channels, and scanning intervals. Placing the styrofoam support plate and the product at the center of the wind tunnel. Connecting the sensing thermocouple wires with the extension thermocouple wires by matching the nodal number with the channel number of the Fluke Data Logger. Starting the electric motor and the fan to initiate the freezing process and, at the same time, starting the Fluke Data Logger to monitor the temperature history. Removing the ground beef model from the wind tunnel and refrigerating room after freezing process. 5. RESULTS AND DISCUSSION 5.1. Comparison of Numerical Simulation Model with Finite Difference Method During Food Freezing The numerical simulation was compared with experi- mental data and results from the finite difference analy- sis used by Lescano (1973) to describe freezing of codfish fillets in a flat plate geometry. Since Lescano's inves— tigation was conducted for one-dimensional heat transfer and the finite element program was for two-dimensional, some modifications were made to accommodate the comparison. A rectangular geometry of codfish fillets with a ratio of length to thickness being 12 was used (Figure 5.1) and insulated boundaries were applied in all sides except the product surface. This resulted in heat conduction along x-direction being negligible compared to the y-direction, which was the product thickness and a one-dimensional heat transfer problem could be assumed in the product. Using the same product parameters, the predicted temperature history using finite element method was com- pared to results from Lescano's prediction and experi- mental data in Figure 5.2. In general, the results indicated that the finite element prediction was in good 67 .mucoEme amasmcmauu xmameflm SDHB muumeomm anm you beam unmEmHm muflsflm Hmcoflmcmefiolo3u one .H.m musmflm :8 0m 50 ON :8 0H Eu 0 68 EU EU 69 .mcflummum mcflusp mumHHHw LmHmUoo mo wuoumwc musumquEoB .m.m musmflm monocfifi .oEHB ow om 0v Om om OH o Anhma u d d1! W a - OMI mom.o 03H m\E m > 8 HH.O DEMO Dom.hNI B a xNE\3 M.NOH L Uo®.mH .B ONI .OcmomQvao;qu mocwumwwwp muwcflm IIIII pozuws ucwEme muHCHm “mama .ocmommav mama Hancosfiumdxm . . . ON 30 aznnexedmel 70 agreement with both experimental data and the finite difference method. The "Student t-test" (Netter and Wasserman, 1974) was used to quantify the differences among the experimental data, the predicted temperature by finite difference and predicted results from the finite element simulation. The results (Table 5.1) indicate that the finite element method provides more favorable agreement with experimental data than the finite difference method. The finite element method will predict the same temperature as experimental data at the confidence level of 95 percent. TABLE 5.1.--The "Student t-test" for Experimental and Predicted Temperature in Codfish Freezing Degree of t-tab1e* Treatment Calculated t Freedom (a=0.05) Experimental vs finite element method ' 0.48 24 2.064 Experimental vs finite difference method -5.25+ 24 2.064 Finite element vs finite difference method -2.29+ 24 2.064 *t for two-tail test (Neter and Wasserman, 1974). 71 It has been observed that the finite element method gives more accurate results than finite difference method and converges to an exact solution in heat conduc- tion problem (Emery and Carson, 1971; Bruch and Zyvoloski, 1974). Hence, the results predicted in this study (Fig- ure 5.2 and Table 5.1) are similar to the previous publi- cation. 5.2. Verification of Computer Simulation with Experimental Data 5.2.1. Freezing of Product with Elliptical Geometry The results of computer simulation using the finite element method to predict temperature history of elliptical geometry product are presented in Figures 5.3- 5.5 at various air velocities; 7.4, 11.3, and 15.2 m/s. The predicted temperatures were compared to data obtained experimentally from 3 replications at each air velocity. Experimental data were tabulated in Table B.4. Tempera- ture histories are illustrated at three node locations including the center of the ellips (node 1), 4 cm above the center (node 2) and at the surface adjacent to the product center (node 3). The standard deviations of the temperature measurements at one hour time step were com— puted and presented with the experimental curves. The results indicate that the standard deviations varied from i 0.l°C to i 3.5°C. Higher deviations occurred at the 72 .m\E ¢.h woman “Hm um mcflmomnm mcfluoc uosooum Havaumwaam mo wuoumwn musumuwQEmulmEflu one .m.m ousmflm muson .mEMB OH m m b m m v m m H o 1 . j a u . u d 4 CM! ;.K kw (AMW7J. mm moo 1 ON: 7’ 1"! I/lflli-lufi .. I.’ I ’1” Fl, // I, I I I . l CHI . o _ , _ OH ounces H u< m.v wusqwm on Hmmwm /r/ /M m\E «.5 > noucmo w>onm monuusm " mm /z .1 Uoomu 8H Hmucwo m>onm 80 v n ma II Uoma fie Houcoo " m /, mumc HmucmEHummxm mo :ofiumfl>mp pumpcmum H ON M NE\3 0.05 : .coHuoHpmum Housmeov lllll n Hmoom..cofluofiomum Housmsoo I.Il dump amucmeummxm om ’aanexadmam Do 73 OH mason m .msne .m\Em.HH ooomm “Hm um mnfluoouw mnfluso poopoum Hooflumwaao mo Snowman ouououomfiouloeflu one .v.m ousmflm ouonwE H u< mxs m.HH > UOONI 8? UomH MB m.v onsoflh on nowom uoucoo o>ono oommusm Houcoo o>ono Eu v Hounoo ouop Hounosauomxo mo nofluofi>o© ouopnmum M E\3 n.Ho.m .nofluoflooum nousmeoo lllll N n HmooH .cofluoflooum Mouamfioo muwc Hounosfinomxm m m N H m H ' mm oboz ’ ma owoz m ocoz OMI ON! OHI OH Om 'eznqexadma; Do .m\E N.mH UQQLm nHm no mcHNooum DCHHSU noononm HooHuQHHHo mo xuoumHn onoumncoeonuoEHu onfi .m.m onqum mnson .oEHB m w b w m w m N H o Cm: 74 1 H H _ H 1 H H H ouscHE H u< m.v ouson op pomom / OH m\E m.mH > noncoo o>ono ooowusm " mm UoONw 8e uoucoo o>ono So v ” mH oomH fie uwucwo ” m ON anon Hmucoswuoaxc wo :oHumH>o© buoccmum H u NE\enm.mvH m .coHuoHUouQ housmeou IIIII n HmcoH .coHuoHnouQ ucuocEOU .I II must HoucoEHuomxm om 3° aanexadmaL 75 initial freezing point of food product and decreased when the flat plateau of the freezing curve ended and the slope began to increase. In general, the computer prediction was in close agreement with experimental data. By introducing variable local surface heat transfer coefficients, a better simula- tion was obtained than when the surface heat transfer coefficient was constant along the surface of the product with elliptical shape. The discrepancy between the simu- lated and experimental curves can be attributed to fluc- tuating air temperature (i 2.0°C), as indicated in Table B.l., and inconsistency of product density which might have occurred during the shaping of product into the elliptical geometry. Measured surface temperatures (node 3) were higher than predicted values due to place- ment of thermocouple junctions on the product surface. It is anticipated that the temperature sensors will give higher temperature reading if their locations are not exactly on the product surface but several millimeters under. A statistical test was conducted using the t- distribution to evaluate the agreement between experimental and predicted temperature history at the center of the elliptical shape. The results in Table 5.2 indicate a good agreement between the experimental data and computer simulation. The finite element simulation predicts the 76 TABLE 5.2.-~The Results of"Student t-test" to Evaluate Agreement Between Experimental and Predicted Temperature at the Slowest Freezing Point Location ' * General Shape Air Speed Calculated Degree of t Table m/s t Freedom (d=0.05) Elliptical 7.4 -l.63 16 2.12 11.3 0.89 15 2.13 15.2 0.95 15 2.13 Trapezoidal 7.4 -3.59+ 10 2.23 11.3 -4.49+ 10 2.23 15.2 -1.44 10 2.23 *t for two-tail test (Neter and Wasserman, 1977). same temperature as the experimental data at the confi- dence level of 95 percent. The isothermal fields at 2.5 hours after the freezing process started are illustrated in Figures 5.6 to 5.8. The shape of the isothermal fields from the experiment are closer to the slowest freezing point for the portion of the product facing the cold air stream, indicating a higher convective heat transfer on the sur- face of the upstream portion compared to the surface of the downstream portion of the elliptical product. Com- parison of experimental data and the computer predicted isothermal fields using average surface heat transfer 77 .m\E m.HH ooomm uHm no mason mnHNooum m.m Houwo uosnoua wuuoeoom HmoHumHHHo oonnH mpHon HoEuonuomH one .o.m ousmHm muscns H no m\E m.HH > x~E\3 0.2. n conuoHpoud poundsou I II Uoo~n 8e n HoooH .coHuoHUouQ nousaeou lllll UomH He came HmucoEHuooxm HH< cHOU 78 .m\E v.h oooam HHo um muson mnHNooum m.~ “ovum uosooum wuuofioom HooHuaHHHo oonnH moHon HoEuonuOmH one «uscns H u< x ~53 ES m m\E v.5 > Deon: 8e n HmooH .coHuoHpoum Housdeou 9.2 He mama HmucoeHuoLxm .h.m musmflm .1 79 .m\E N.mH poomm uHm um mnHNoon muoon m.m nouwm uospoud muuoEoom HooHumHHHo oonnH moHon HoEuonuOmH onell.m.m ouquh ouscHE H u< x e\z m.~vn m N > m\E N.mH Uoo~n SF n HmUOH .coHuoHcouQ housnsou I ll UowH He camp HoucoEHuomxm U ./ (w\ 80 coefficient and local surface heat transfer coefficient revealed that local heat transfer coefficient was in favor over average heat transfer coefficient. This was supported by the result of mathematical analysis using Katinas' data (Katinas et al., 1976) as described in Appendix C. 5.2.2. Freezing of Product with Trapezoidal Geometry The temperature history, during the freezing proc- ess of trapezoidal geometry product measured in the laboratory was compared to the computer prediction in Figures 5.9 to 5.11. The results of the numerical model was in close agreement with experimental; the standard deviation varied from i 0.1°C to i 3.l°C. More favorable results were observed from utilizing local surface tem- perature (Appendix D) than uniform surface temperature. This implies that the heat transfer coefficient varies as a function of location on the product surface. The isothermal fields presented in Figures 5.12 to 5.14 indicate that the heat transfer depends on the local surface temperature. The numerical model incorporating the surface temperature as a function of location and time gave more exact results to experimental data compared to the model using uniform surface temperature. The t-tests for trapezoidal (Table 5.2) indicate that the predicted temperature at the slowest freezing 81 mnHuso noncOHQ HmcHoNommuu mo euoumHn ououmuomeouloEHu one .w\E v.5 coomm uHo um mCHNooum muson .oEHe e o m v m N H q . q H H Ht I ’l”” I I // I I/ / Al / I . I / ouchE H an m\E v.5 > UOONI 8%. A muop HmucoEHquxo no :oHuoH>oo oumccoum EuOWHcs we .coHuoHcouQ nousmsou nu .xvuuume .coHuoHpouQ nouameoo mumv HoucoEHquxm m.v ouqum Cu UomH .e nowou nonfidc Havoc uom .m.m ousmHm Om: 'aaneladmal Do ON Om 82 .m\E m.HH pooam nHm no mnHNoOnw mcHusc uooooum HopHoNomouu mo enoumHn ousuonomEouloEHu one mason .oEHe .OH.m ousonm n o v m m H o H H H H H Cm: m 6 oz H c om: % / i l I / / / / / / / l / a / m x a / / w x m J a m ouscHE H u< m\E m.HH > // ooo~- as m.v unsung on “com“ /,I UomH He hog: Havoc NO...— /. ON mump kucoEHquxo mo :oHuoH>ou pumccmum H EMOHHSJ me coHuoHpoum houses—ow nnnnn Ha .wi ume .coHuoHvouQ nousQECU l.l| numb HmucoEHuoaxm Om 83 .m\E N.mH Toomm uHo no UCHNoonw ocHuoo poocoum HopHoNoQouu mo >uoumHn cuoumuerouuoEHu one .HH.m onque muson .oEHe h o m v m N H o Om: Om: OH: 1 a .m a J F 3 o n J O m CH ouscHE H m\E ~.mH > UooNI fie m.v onoon Ou nouou UomH .e muonfisc Hobo: non ow mumc HmuccEHquxo mo coHumH>op cumccoum H 4 ELOWHCD me .ccHuchonL ncuzoeoo cccccc Au .xvu u we .coHuoHcona nouscscu l.lll moan HmuccEHuoaxm 0m 84 .m\E m.HH poomm nHm um noon mnHNoouw o.H Houmm nooooum HooHoNomoHu oonnH mpHon HoEHonHOmH one .mH.m ousmHm EMOMHCD me .noHuoHpoum nousmEOU lllll M E 3 . Nova W.M% .m H» .xv m u we :oH oH ohm no smeo .I.II UoONI 89 .H .mu U U UomH He muop HoucoEHHomxm \ MHn CHOU AWHU 85 .m\E «.5 coomm HHo um Hson mnHNooum o.H Houmo uoocoum HmcHoNomouu oonnH moHon HoEuonHOmH one .MH.m ousmHm ochHE H u< Doom: 8e Hu.xvm n e .coHuoHooum Honsmfioo .I.II m\s v.e > oomH He mama HmucmsHumdxm HH< GHOU AU 86 .m\E m.mH pooam HHo um Moon mnHNoon o.H Houmo uosooum HmoHoNommuu opHmnH mpHon HoEHonMOmH one .vH.m ouomHm ouncHs H pa ooomu we H» .xvm n we .c0HpoHooum Hopoaaou.l.n.u mxe «.3 > 003 .e 38 Egg! 87 point location, for air Speeds of 7.4 and 11.3 m/s, do not give results as good as for 15.2 m/s air speed. The discrepancies can be attributed to the method used to average the thermal product properties of a triangular finite element after the initial freezing point was achieved. While this averaging method may provide insig- nificant error for a small size triangle, the error will increase as the size of the triangular element becomes larger. Thus, the predicted curves tend to be decreasing slower than the experimental curves beyond the thermal arrest time. 5.3. Influence of Area Average Enthalpy on Freezing Time Freezing time can be predicted conventionally by expressing the temperature history at the slowest freezing point location in the product or by using the total heat content (enthalpy) of the food product. The conventional method for predicting freezing time represents the length of process required to achieve a temperature equal to the desired storage temperature at the slowest freezing point location. An alternate criteria involves time required for the food product enthalpy to be reduced to an enthalpy equivalent to the storage temperature (Gorby, 1974). When the freezing process is stopped at this point, the product 88 temperature will equilibrate uniformly to the storage temperature with the heat content from the portion with higher temperature transferring to the portion of the product with lower temperature. The mass average enthalpy, H can be expressed as follows m Jf Hdm (5.1) m m n BlH Implementing the above equation into the finite element method, an area average enthalpy was used as a criteria to determine the freezing time (e) (5.2) The predicted freezing times were compared for both conventional method and area average enthalpy (Table 5.3) at various product geometries and air speeds. The storage temperature used for these comparisons was -l7.0°C. The average enthalpy method predicts lower freezing times, ranging from 10.7 to 24.7 percent com- pared to conventional method. These results may encourage the use of the average enthalpy method to predict the freezing time in order to achieve a more efficient process and to reduce energy consumption. 89 .erHI ousuouomfiou ommuoum« me.m m.v m.mH no.4 m.a m.HH mm.v o.m 4.5 HmUHonamue om.m m.» «.mH me.m m.e m.HH mm.m H.m v.5 HmoflumHHnm wonuwz smHmnuam Begum: ommuo>< moan no pommm HonoHuno>nou no pommm m\E omonm Houonow 666mm “Ha mason .oEHe mcHNooum pouoHoon «mooomm HHd onm momonm msoHuo> MOM oEHe mnHNooum couoHcoum new HounoEHuomxmll.m.m mqmV amHmnuno ommuo>o mono no oomon oouonHEuou mmoooum mnHNooum man man» man um uoououm HmoHuQHHHm monmcH mcanm HmeHmnuomH was UOONI n HmooH ouscHE H u< m\E ~.mH > 0.2 .mH.m musmHn 92 .Hm\E m.mHnu>v emHonuco ommuo>m mono no woman pouonHEHou mmoooum mnHNoonm onu oEHu onu um uosnoum HoUHoNomouu ochnH mpHon HoEuonuowH one .mH.m onsmHm n HoooH ouanE H u< UoONI 8e m\E m.mH > UomH .e UomHI uHfi CHOU AU 93 .uoooonm HmOHumHHHo mo oumu mnHNoouu on» no oNHm HMUHHuoEooq mo oonosHmnH one .eH.m ousmHm mason .oEHe vH nH 0H a Q h m w v n N H o 1"! d J 4 4 d 1 1 d d 4 CHI ‘MN'II. Il‘l'l'M-N L ONI z/uff. / I I H * oHa . _ o W. ouacHs H no n.v ousmHm Ou nouom OH m\s ~.mH > noncoo o>onu ouuuusm u mm UoONu 8? noucoo o>ond EU v " mH UooH He woucoo “ m . om oNHn HoucoaHuomxo x m uuuuuuu ouHm HmucoEHuoaxo «Hm HancosHuomxo x m.H ll l.ll On 'aznaezadma¢ 3o 94 .uoscoum HooHoNonuu mo ouou mnHNooum onu no oNHm HmOHHuoEoom mo oonoonnH one muson .oEHe m b m _ m v m N .mH.m ousmHm H 1 H H H q H ouanE H u< m\E N.mH > UoONI 8e m.v ousmHm ou OomH He Homo“ nonfion Hoooz oNHm HounoEHuomxo m. lllllllll oNHm HounoEHuomxm oNHm HounoEHuomxo m.H III I III OM! ON on 'eanexadme; Do 95 useful information in the selection of optimum size and conditions for efficient freezing time. 5.4.2. Initial Product Temperature A sensitivity analysis was conducted to determine the effect of initial product temperature on the freezing time. Initial temperatures of 22°C and 14°C were compared to 18°C, which was initial temperature used to generate the experimental results. The results, presented in Figures 5.19 uD5.22, indicate that initial product tem- perature does not influence freezing time significantly. The isothermal field after 2.5 hours illustrates a devia- tion in the range of : 2.0°C for elliptical shape and : l.5°C for trapezoidal shape. The results, presented in Figures 5.19 u:5.22, indicate that the freezing time is about the same for various initial product temperatures. The freezing rate curves converge to 7.5 hours of freez- ing time for elliptical product and to 4.0 hours of freezing time for trapezoidal product. In general, it could be concluded that initial product temperature has small influence on freezing time. 5.4.3 Time Step Time increment, At, has been known as an important factor in the stability of numerical analysis. Two .Hoopoum HMUHHQHHHo mo oumu mnHNooum on» no ousuouomEou HMHHHCH mo oonoonnH one .mH.m ouomHm muaon .oEHe 96 w h m m v m N H o d H d a (q d u d OMI IIIII Il.// omI I I , I. I . mm 6602 / . I / ,/ r 2- I m .d a 1 o w / /_ . m6 ouanm on Homom / / «w ouanE H u< Hounoo o>ono oommusm "MN I // H / OH m\E N.mH > Hounoo o>ono So v "mH /, / UooNI 8e uounoo "m [I U vH HE u o H I . om UomH .e H UoNN .e I I III Om 97 .muson mnHNooum m.~ Houmm Hoopoua HMUHumHHHo on» nH monouwuomsou HoHanH msoHuo> MOM moHoHu HoEuonHOmH one ouscHs H ac oovH e IIIIIIIII mxs «.mH > oomH “a Uoomn 8e Uomw He I. II I II .om.m ouomHm °C Temperature, 98 30‘ ——-—--—-T. 22°C 1 T. 18°C 1 20 H ------- T. 14°C I l I M O -30 Nodal number refer T0° -20°C to Figure 4.5 v 15.2 m/s At 1 minute ~ ,_ ~‘ Node 15 1 l l Figure 5.21. 3 4 5 6 Time, hours The influence of initial temperature on the freezing rate of trapezoidal product. 99 .Hson mnHNooum o.H Houmo pooooum HonHoNomoHu onu nH monouoHoQEou HmHanH mooHHo> MOM moHon HoEHonuomH one .NN.m ouomHm H M#DCHE H U< UOVH HE lllllllll m\E mH > UomH e 8 H UOONI 9 UONN E.III.II||| qu CHOU AH“. 100 different types of time step, one minute and three minutes, were supplied to the simulation program. Figures 5.23 to 5.26 indicate that the numerical solution remains stable, at least in the range of At equal to one to three minutes. The discrepancies shown are small and insignificant; 0.0 to l.0°C for the trapezoidal shape and 0.0 to 0.5°C for the elliptical shape. Isothermal fields indicate that discrepancies exist between the time step of one minute and three minutes on the product portion facing the cold air stream (Figure 5.26). This phenomena can be attributed to the fact that thermal properties of product were computed at each time step in the simuation. Large time step will provide inaccuracy which becomes more significant for higher temperature difference between the product and the cold air occurred at the upstream portion than the down- stream portion. 5.5. Application of the Numerical Model in the Food Freezing Process The finite element simulation model has been veri- fied by the experimental data as described in previous sections. The application of the model in food freezing is not restricted to trapezoidal and elliptical geometry but to other anomalous shapes as well. For practical value, an observation of size parameter influence on 101 .uoopoum HmoHuQHHHo mo ouou mnHNooum on» no oEonom HonuoEs: CH poms mono oEHu mo noowmo one .mN.m ousmHm mason .oEHe m b O m v m N H O MN oUOZ mH ovoz m.v ousmHm on Homom Houcoo o>ono oomwusm u mN m\E N.mH uounoo o>ono Eu w " mH UoONI 8e uounoo u m UomH e oHDCHE H u< mouanE m u< I I I OM! ONI OHI OH ON Om 102 .muoon mnHNoouw m.N “ovum nonpoum HmoHumHHHo on» CH mmoum oEHu unouommHO 03» How moHon HoEHonuOmH one .vN.m ousmHm m\E ~.mH > 0.8.. 8e 32:... H u< UomH He moaacHE m an IIIII 103 30 ----- At 3 minutes Ti 18°C At 1 minute Too -20°C 20 Nodal numbers refer v 15.2 m/s to Figure 4.5 °C Temperature, _30 l I j l 1 O l 2 3 4 5 6 Time, hours Figure 5.25. The effect of time step used in numerical scheme on the freezing rate of trapezoidal product. 104 .uson OCHNooum o.H noumm posooum HooHoNomouu onu CH mmoum oEHu unouowpr 03H MOM mpHon HmEuonCOmH one .om.m ousmHm m\s m.mH > ooomI 8e ouscHs H u< oomH e mouanE m u< IIIII qu CHOU AWIIIU 105 freezing time can be conducted for any given product shape. Then a chart illustrating the relationships among freezing time, size and shape can be developed and a con- venient method to select freezing times for desired prod- uct, shape, and size would result. Requirements that must be satisfied when using the numerical model include: 1. The initial product thermal properties must be available either from literature or experimental measurements. These properties include thermal conductiv- ity, specific heat, and product density of unfrozen prod- uct and freezing point, water content and unfreezable water content. 2. Average heat transfer coefficient should be known or measured accordingly over a range of air veloci- ties. 3. A finite element grid should be developed for the given product configuration. Again careful judgment must be made since the coarseness of the grid might affect the accuracy of the estimated freezing process. Smaller grid will use more core memory in the computer and produce long execution time and in turn high computer cost. For example, the elliptical grid in this study used 360 to 390 seconds execution time and the trapezoidal grid 140 to 190 seconds for total computation at the CDC 6500 computer. 4. 106 Computer facilities must be close at hand to run the computation. The finite element analysis has advantages and limitations when used to simulate the freezing process. The advantages include: 1. Easiness in changing the input boundary conditions and parameters including product properties and freezing environment. Readiness to accommodate various product shapes and sizes. The limitations to applying the simulation model approach are: Modifications are needed for a specific anomalous shapes and for non-homogeneous materials. The model is valid for food freezing utilizing air velocities in laminar region with Reynolds Numbers smaller than 2.1 x 105. Required more computer time than finite dif- ference method, thus higher computer cost. The method used in averaging nodal product thermal properties for a triangular element causes some error in the simulation scheme as previously discussed especially in large grids. 6. CONCLUSIONS 1. A computer simulation utilizing the finite element method to predict the freezing rate for food products with elliptical and trapezoidal shapes has been developed and verified by experimental data. 2. The computer simulation using the finite ele- ment method has the ability to accommodate various bound- ary conditions as well as the non—linearity of product thermal properties during phase-change, and different product geometries. 3. The utilization of local convective heat trans- fer coefficients over the product surface as boundary condition in the computer model provides better simula- tion of the actual freezing process as compared to an average surface heat transfer coefficient. 4. The predicted freezing times based on area- average enthalpy method are 10.7 to 24.7 percent lower than the predicted times based on the conventional method --the slowest freezing point location. 107 108 5. Geometric size affects the freezing time sig- nificantly, while the influence of initial product tem- perature over the range of 14.0 to 22.0°C can be considered negligible. 6. The stability of the finite element algorithm is not influenced by the magnitude of time step at At from one to three minutes. 7. RECOMMENDATIONS FOR FURTHER STUDY 1. To develop a chart illustrating the relation- ship between freezing time and various product shapes and sizes. This chart would have practical value in selecting the optimum size of a given product shape for an efficient freezing process. 2. To modify the computer program in order to accommodate the ability to simulate the freezing process for non-homogeneous food products. 3. To determine the ratio of local to average surface convective heat transfer coefficient for product shape other than flat plate, circle, and ellips. 4. To utilize higher order elements and a non— consistent capacitance matrix to investigate the possi- bility of getting more accurate and stable results than using a simplex triangular element and a consistent capacitance matrix. 109 APPENDICES 110 APPENDIX A COMPUTER PROGRAM TWODFR AND AXISFR lll TWODFR FORTRAN LISTING FOR ELLIPTICAL PRODUCT --TWO DIMENSIONAL FREEZING 75:00: IhPt‘Ing-FE:1=CUYPJT, YAFU. .- \h u 1700 'rfinr‘owy v-«A- If C)‘ ’CUC‘U“ A! I nt-Jnanmgm “3.02400 0 8202222010. P‘C‘XC‘C‘CCO II. a 3. xsxxz‘vc, § r/bz 'R. 1.1530“, CI otagKigKI.CPI,IF.; WC Uri-1C5 KSgHS;CFSo 1 Y 10 220 . uvvvp ary In... F-OH ("--CC:fl Xfif‘OO-‘Al"! D! 15 r O dal I’ " I u ... u 9.4: up. (¢(((((C(( (J (' 7M0. 20 " IO!’ M as :32: CA ‘ 1., .1 ' 30 Lmvbv-oa NPPI\;’N$~,!I.. .IH", K‘x,‘YY,H'::I'F’ 079.1%“; QUO‘? 35 N u 26"!” W." ’CMPI S‘Cr. ‘ ’7‘. 0’. Men: vwou *W I. “IN‘UW'JU.IU N M " N‘ "O D'DIWMH/ lint!" .40 10 lqnflall *4" H o (' ma t‘ahJO Naatm- r ‘ ‘0. HO CV! (‘VO .0 .' LA. 0. n OI POINTERS Aw; INITIBLLZHI;JN or THL catuqu vECYtha ( 4. I LC -| “F OlP’h- - 9|:P'hJ" L (a(J'. N W) t' n- vuu “‘1 nu LW‘I’IO I!“ n 3 I 45) '0' 2 :ztrtes -c 0 I»: :'r v 1.:- c 0» Cu U 11 8" - - 7W zh: .IX ace-will 51K!) =3630-lll91 “FLU; 5‘.~HY‘1”6A’.hl‘ H” 07;".79 IiTx C {‘0 3 ptlgo LZATIUN; -9;U’2Xo 30: ..T£- I” 1’ 11C 6. o ‘70-1I1‘.:7Hh:5 N! :5 N: 1'" (( :*F3-9’1£’12?‘“Y' '- ,‘S..,1 ‘- I ) V, (f FP,"} yr A. L 2) ox,~nv(:).o‘ubhl\3’tb'14 “NP " NVIZ'J 7.7". O N to“! ’13:” “.hr" o-ounu 1"."- "‘048 2A.. .11 IIHE SIEP = =1o.5/ “-51?“ : EL . r%." .FP; a ( ‘. ( U NF- C 4 3 yr ,3 v. :r... p O ELHL fk-~T 0F EsEthT DAT; I O 2 VU' (f “I RE?’ )),r(1),x(2),¥(2l,l(3),1(31.:SZIE é:(1),7(1),x(2),1(2).H3),Y(3)o 151.}! "-”V’ l U tt-ONU'II 2H 137-4 40>. «A «JUV0* ,f' O 'H 'l N C M c. NC “N H O-‘I " rrr1‘r I 3. " epur- XH ( 6. I Hfi mu!“ l-r) "C AANX I In M!- ’U ‘|" m V V 75 b;. yKAXp KYY 301.7th ,NFoJGJ'HJuCFHUT 9N. ybkaVl'ouE (~50 ) i4(JGC‘OI),A(JGF01),a(JTnox3,;(1)1~P, N5~,N;- Co CA .1 01K; AM HgT‘NF’ JbS‘yJL'H. "VE'JL~;,Lv’ -a-4m< 9,,qu 11L2 11.3 SIDE(3) . ‘ H57 (A,CP.HJ ,KXX. K77 ,h,T1NF,NP, JbbflvJGLHOCT’NEpNhgl‘Bh, ,NS(1).VStZ).N$(3).x(1),V(1),X‘ZD.Y(Z),X(3).V(3). ox(3)'v(1)oX(1y-V(2)-x<2)'1(1;-x(3)'v¢z)-x¢1)'V(3)) 3(2)och(s)'P;J«3))/3. :13 aHSH-‘Uvol I730‘19K'15Flo7fi;“GgU‘UCgCS vKSflflSgCFS. - - 3 ‘ or E ’ ’ 3 . 3 ‘ H s E . ‘l b 3 Q I. A 3 / ( n, H .1 c J E , l\ 9 3 C I! 0‘ L O 3 J n 1 (I o T I c P s c 1| 9 I .1 7.5), 33 Y. o .a I. \I I 000 [P 3 U 3 3 T 000 \l. .L I '- l\ l\ \I R 111 \l \l) h ’ V. 8 J 2 .L O 9 9 1 32 A \r ‘ 9 ‘ \I .1 D F. 111 7.x- 1 (1\ 13 .0 \l 9 O 2 .U.\ 111 0 i 1 JJ Y 1‘ O.) 3 1 E I PU Urns 03 a) K; 9. Eu. fiJ 1‘ 3 J X 9..- 333 3D. . 0 o . s In ( ‘5 l‘ l\ 2 u I. I. G. n. 1 )0 on 1L U KN) J .l 9 7K 111 GK 1 2.! Mn 10 To ’ 1.. 05 L q ”5. 111 T.‘ 3 *1 0|. 6) O O \I \\r( 1 s ‘7 a L. \I o .P 000 l I 1‘ .I. l... a. 0 *1 7:01. 9) l. 0. U .1 L20 NH 333 CE. 1 KJ .c 5(l 23.( duo. 1):) FJ PT v El\1 14 an. n.” )1) 2L .9? nu O s .Ulfi :3... T1E1L31 SK N ’E HEP—(2222211111), 123 LT. \l” l.— 0 \IP ‘ UHF. KO 0 E OJ3A‘1nD H I (IE/J 0.0 o o o o 0 03122313 .LEL .L o o I3 s: l 3 1. F. 1:55 .ILX 37(1VJHI lv. 1.: at. E); 10033.50(o\(1\((( Ouh." F 222 J :0 (.1 F.7- oo o hX . .LXTgllrr 1V : 07L N fiz3333333¥YY11yzY (13‘ F0 3331 2‘ 1 J1 at. 1).; u \er EN CrUFK 001‘ 90.1II/III- . o . o .0 LLL 0 . . n = 31A 0 1 vs K‘ .L 01E! JC‘ N ,PVNHNK,KLL; I L1,..." 10 'b.\1\11\l\l)\l\l\l\l\r\l) 1), h‘h ,J 1))“ opp 511:?» l h (J F” 1“ h o: :30 . .JUQD JP .dl ZJ 1751.. 11(1231232313122 231VVV NH 123 .IOCKPZ no .1 ..E: . 0L 210...£33 0)) I o T b . .1 o. ./ 0.112 o 2 = U E 11 ‘(((r\(((((.\(( (((TT-l .L U (((11bfi .. A: . :3 = K = U) L . . H15 U l 1 .CJJ10 UNSS.DSC.~D NKab a 00(1)No.11‘oleuXYv-Y'YYKXN‘ 5.35: t : .91 JJJON up):\lK\l 1:1VXI‘ NE KEH 6 N11 0 O O: a .U 'N \NN U1». .0 0.9.0 0N1 IDIN (.1. : s .. .. = 8 : 8 : 3 : :( .\ \u‘)’, olol citations: N)!‘ H.JN\IV.AK : .U K1135 7.: .. 0.1.11.1) .KJ...E::..EH1H o0 : .. 31.’\J‘E.lo"hdalv)la|l-lulc’lu”)1 : z : 317.3 5.0 ((‘SJ‘ .(V‘1‘1‘ ‘(t- : 3 v.2? (: ...>ol.1J1l..-|l|11 DNHHHHHH 'H.A: 1.1.. H A RJNH $.éfi1231231231qcs h. 0123((.\ LN .l LJ‘J 1 J‘JNY H 7.0 11 : “3&5: HH1‘ UX .1. .. . .uK DEF HLL.L£U(0 mCK.oO((((((((((((R2EELJJJ UU FFFL:FA .JLF ht .J uLY n « Fhfi‘ bLtKF Lb L... s 000 Ugltvkc. SSPDKF‘ D *XX‘V‘VVDBBLLCAO NNN'IQIT CF .1o1C11rVCKPIL’U‘LKPLKL «Hui 1‘55 LU?.¢£K7UK 1 ‘ 1 1 s -f. 5 n 5 LS 7.81 1 2 1 UR 1 9 5 1 D 1 o c h 001 1 0 o S.‘ 1 O 1 c .3 r“ ,v - 007. 3 D G d . U D \U 1 1 .o .u 113 3 3 .u v 1 a u by a s . v 1 5 0 p: «U .7 c 5 0 ‘4 C r: r; 5. Rd E. . v p: 1 1 2 7. 3 3 h h 2 .2 fl. C 7 7 O O 90 95 130 105 120 125 130 135 lIL4 u~ C£.CLgllION OF IHL IONV-C710N R;Lll;£ 0U~h11f1gs U0! 00101:; 2 1r¢xsxut(1) ).LE. 0) 6070c h=H315F(hEL) J=2 tOL‘.) n=Jo1_ Lo=>2n7((xtfl)- X(J))"2°(Y(H)-Y(J))"2) HL=4'LG EF‘J)8§F(J)OFL'YIN‘/2o %F(1)=LF(H)OHL'S NF/d‘ :H(J,J)=£;H(J.J OnLI L9H(J,HJ=L-H(J,:)OHs/Lo EaH(H.J)=LbH(J;1 ."10 L5H‘q'fl’gl-”(H’~).HLl 30 .. INaLRIZON 0F LLLHLNT PHJP.i71E5 INTO THE GLOBI. S11“NESS It T‘I’ 6 Lb 7.13193 I,=1:( J- JD‘Z’N;*LI A(JS)=A(JE)OLF(ID OU17J=1 3 JJ=vS¢gf JJ‘JJ’- 91 _ 1F¢Jg)17 17,16 . 1o J:=J,s:o1JJ-1)°upoz- ‘ Jt‘JaCH‘(JJ‘1)'N”II v AJ:3| HtTiZX ~. A'TCNEH)=P'T(U-D)-2'F CU. A‘J§’=M‘J5).EJH‘ ‘J’020 0’07:E;" (IpJ) £(Jb[=~(Jb)02o a? Ebn(.oJ) OEJH (3.4) 17 ELNTzth ‘_ 7 LHIINUL UV. u CJ.CU5‘TIOH “‘DY UUT’UT OF IQEA .__aka;GE Eh.HALP . IF ((K-ZIO) 711,712,712 12 ‘uaF-Sfll Pho" W1:1“|:" _ b111 FLK‘ 1-1 ,‘Ar.;. £V.¥Lci nghu.FY'g710.q) ‘711 Gehrlhuc by .. DLH’BD DECtfiPLbE: A -HIJ 09°;fi IxIAh5.E «Alnlx va Lés._gLHPbL(L(J5-P013g\’.Nod) FLTJ! :hu 10 p \h 20 30 01 I" ‘n 50 CS 70 50 11.5 C 4 SH JL'vflpN .JENJQ" to (IV (y 00—0“). \-Q'N :‘U? r 1% 5 4‘5an '7" 0:? w 00"” ‘ Q rt ['5 ”it W¥¥X¢WTUOCODH’ r§ 331 nl(23).18(2 )glp 9 E Z < : hp )NBH) GHH‘NPgNaN’ p‘NP 1’ ti ("Fgl’ 3 9 5 évl 3) 2 ’ 5.92) 2‘712UOORU xmwmmuuyz=1u¢JJ: :fTs;«znc or pacscnzsa: Bourosfiv VILUiS 1n Bevan; via 22 1F“~0Eu00) 007315 0L13lz1’K'I g§.°:V£,) 16 .79H1(.)3X;(J11) 15 Cghrguuc_ >LT=30TO.T AL. HU-‘bL‘GH.‘"X:'T".gh‘g~fi'NCL’ L51=1 N‘ . _. 9 nCI.13=TK$2.1)02.'P(-91i .‘-DE;0HPJSE§ F ANL 53-V£3 Fb§ 'EH’ERETUiE ‘ CAL. S-VbD£6LH,Y-..xl.NF,N6~.«CHID) ZTTaEJETrlus or THE FRLSCHluED aouao;av vn.ucs "' 1f(<~.EQ.O) 507321 Dkzoigl' 'I J3-’:V(-) - zc X-({.1)8bCYPn.(1) 21 LG" sub; KEI§~G 2 Du-710 .K=1 HF _ NL-tE (2 100;) x-(1LK.1) 1.9: FLn4ST(£i7.5) ‘210 Lc4r1~u£ . CUTPLT 3F int RLSJ;?J 1F 355K311N7KD'2dTa).NE.KKI 30 I; 6011 _- 0% ,7, .E%"‘P _ :25 ) (1)8(x- 1)-32.)‘9.I:. HL.YE (b1 7! :.t,(; ’Fti).1=1.h~) LLJJ‘T=-L UNTon°ISO§ 5511 LCLIXNUE , KLCJ\T=3'NFOZ’N9’N=4 DL 501 ;‘11K50UVT At.’=0.c _261 LchTINUL . SE! Tn? NEE £-Ln£ut HA7 .;55 F0\ . In; hi1? 7.H£ SYEF va LA.L SETfiAI¢£gQRiOgKXYgKYY,H.Tih‘,NF.JLSN,JGLH,CT,hL,Ah,hg-,JEuS, IML’NIa‘-VI‘1'V~S) DO .30 -’"=1.'1“ X.(IIH.1)=~(-X') Yh(I‘H,1)8C.C N'.K;=IIHONP'2 F‘Ll1,1)=I-(HLKJ) .30 Luntzwu, LL1=3'ltr91 LC2=3'thN:h'NPo1 UL. 531 95:1.“64 0L b32 tbglg'ip GN1(QE,-5):H‘5J?) 6Lm(.b .E)‘~(LU1’ Lc:=.oéo1 Lul=-0101 ~32 Cyu7LhU; 931 LLI‘YZNUL RLTJxN Ehu '53” (KH’3HH£‘QT‘91IgN‘oNU"N;L3379N¢tgleipNNpC5HQp 0 "$10.03: 9'3";,K~.;?"TF’; “CQU:H:905"5|HS'CP50 3,3loil7otgi‘1.13.319E17ot 35‘ 9:313} lIL6 SEE Iapihpll C ’ 5},"9‘5‘11.’ F1 0 1 i HETA1‘hH) ca.UU..rxou or NEH ruar.. UHH a K.a.r.: M CL. .L thh SCDRFPFNI 1 2 5 5 .u «a . a u 1 s 10 VETh(TH,fHZ) gététfié UUU G c D . 0 £5 2 .3 . Av \l O 3 3 O o O I! .4 3 H O .1 O 0‘ av U H .9 9. a I. . O F.) 6 5.0 o .30 c 0. L») 1.9 90 O H 30 ‘v- 0. 2‘ CC 0 O OH U 2 ’9' :1 2‘ H. O 0 vii O 2 (6 .91 08 H. Q. ‘E n w ‘3 o o O a El: 3 o 15 03 7.0 . . 0. . it 5C .3 .H 2. .vf... o. .11 :Jav If. on... O H .U U 3 xvvl a 30 H1 .r. 0. To .3 E I! .E . g 5 V ll 2 £1 ’0 I is .\ 4.. .. .Q. 5. 5 ZN OI .55 o 0.: .9 H 7‘ )7. 3.1 as l . ll. "Inf”. .1 7 w «LL )ol’i/O’ 1 O 5 7.! N(2 09% .522 F IL LHhhaoo no. .f I! \I 0 1H ‘TQIVI 0. '0) t. 5 K 77 1EOJ1C .25 H 1 K N 1.35V V9..Ht...o 3E: .111 0 9 07" .Hc’vladO’C wL£.UJ NO.” I. .0 £20.!"th BLCH ' ‘ NVVN Y h534 (H.031LHYCITi ‘N‘N ‘ ~1ltL:l“?-I:‘Yl:rlrlHJ =.... L 1.. z .9 F3 .vF‘L N 10.9 FhK s U FH.Uh.th. .. LP. a. PLLILN INKL C :TGTLTlHEGTh7VGLPfiE a 1 1 37 A .t 13 .h h .b 5 55 c 7 7.7 7 5 5 8. O .U .u a v v v r: 0 C. 0. c. G 1 2 2 3 3 ~ 22 \I - 0 o 7L I O \l . \I) .5 . ‘12.... II .a... .1 App‘ H (4‘.- T O .1 1| 2)., .0 . I II .V O I 2 .U \I 1”. O 5 7374. a 5 .Ki) 5 O h a? 4\ 0 111.0 . .6 (l‘ D 2 0 II. . .‘ 1...): i In 1.18 \t O 33‘ \I .I o..16 .Iv " ‘ ‘afl’. Z O... o ((8 h .1 2 O O O I: I: o ‘ c.c..flt . H b)” 7..., 3 3. vl O "h... n.- u D c o... at 18.17....) ((0 “1 P A. .h!\..\0h1 6 o . O I) 1.. Y: 015N191, it): o I. 1 o .L 31.. ..AZ 33. 21 [2 V ,b/rvs 1: .PF‘ (6 VII .. 320.20‘ (uhh ID(/ 2 /. )\.l. ) - \ll‘l“ 6‘2 -0 x ‘ ' f... 12 n.1,. . l\ . .v L0 E D 9r)PR)kF)X())Ttu 1 L6 NP L; b111A‘P1 911.. V5. 0 £3 . . 1 o.£.l~l...((( . P... .d 631. 0 7X3 ~23 .rd 3 X...AH #850 . .EEZ JP. 0 o 1. I09‘ OLRE.1((O 01V .HN 3 .13 o .41.).3‘z: 3.! 5(13‘ 1E . OIQ I‘D: 11.0311 8 : 8 z ‘5: .. a .. ol.‘ : V: : .. J .04.... z : : = 21122 5 53322:. .1112]..- . UUL.AL.fif.L-. .t. .::L....l. .LhrL’th. .n SCVAEPTEAkhflPXaAbePV.VTT L 21 6b f. 'n 10 (:4 | H 1n ‘7‘ (1 u 1" 60 11.7 AXISFR FORTRAN LISTING FOR TRAPEZOIDAL PRODUCT --AXISYMMETRICAL FREEZING PnLSian AX-SFR(1NPUTthTRUTgPUNJH,TIPEOO=thuTgl‘Ffié:=OUTPUTgTAPEt 12=FJ~C§ 11111 11r£2 .1:53,11=£¢31 8:n£us 9n n§g§1 x13!,v«31.1s;c£131 canon? 15).- .E120: Cun13~ 5:1351) COHQON A(1) - . . . - ‘Ekpsggo EggF/bAoQ.nH5gU-.DI.73,KI.KR,CP~,IFHHC.UFI‘..CS.K>.HS,CPS. tv - acfi10u§vv1.vv¢121 ‘ EA. KXX,KYY,LG ."8Erlnzrgon or La: ‘gutnug =A§1nsrsas . u — nun ER 6 0190.. £n=£a.ruazs. u: - nunasa or . c.cncnt; ~u~_- 6.x; d1011 xxx - .ou.u;11v1vv 1n 1 DIPLSTZLNv . xvv - CuhuugT-VIYY ;~ v ozkacrxcn. n - couvtcrxcn cocrr..1£nr. . IINE ‘.5§UIu TEMPELLTJRE ' LR"?§§'EC?FEL "£37;0£“5;3; tsuarzcn a = :N00 3 .c c Y!~z=1?ziaf1ou: béFan =RIN¥th ' zo=§1 31:25 ‘ xx:1 .”' «>17 3= Tn: 117*: :aao an; tn: co~1a01 p1a1n37£a3 v... RlVL 3 “1-710.. E‘IPERICJUK; "' REL)(XZ.3) 7:1-5 3 wa'h‘téZO‘n-r’ ' HLADtI ,2) up u£.~uu,uzr.1ura.xxx,xvv,n,t:nr.07.:a»a 2 FLA‘IT (523,6’1 . 1h71=10 REugho ; DL 310 -HO=1.NP GLAUESbn. NLLYE(Z 10u5) GRAVi 1305 F.&1LY 1:17.51 “:10 CLNIINU£ focA.Lu.1txon or POLNTERJ :1: 1u17111121115~ or was :31un~ v::7:a:& "' JTN=~F JoF=?.NP Jgu=JQFONP Jb.'!=Jb'.H0N°’h‘E- 4-1!3: Jbzro1zp.h.‘“ Lat-5.. SE. r L (a(JLV: )) D.131=1.J£L. 13 A(.!=0.J .j'ourput as 7111: AhD 011. nanoz~5s UV. Hr.1":('.0gh) T:T._,KX),KYV.h 0'1": 8 FL'.‘:T(1H1glI/o'1) Z‘Jv-flc/1xo 1KAA -’- F3.~giO’g"‘"' : oF3.-:1A'."",_ 1NVE;YIUH CL'fiFF 3 alRt'l"12H:Lt/;c; *il‘; ‘9‘“. ~/I1‘.17"~"L. l!.-E N» ZFCLitbxl‘rhkl1),6 3N?"(1)96XQ¢HX(2)’6I9~“Y(2)9519~Hll3)9313~"7‘() m.-!E(;O,62) $fi"».'g3?gf11T91“7\ 62 FLFN£Y(21'21H :‘lNO “Hub/21.13172“ 37E: =,'1C.‘/ %?"19th L? lTE'.i~T-.NS =.1~I21,3on1tE—.T;L~; fiftnffh pKith; :, 1‘! 1.; ..~1h3u7 IND ELHL PR4h? CF ELEHEhT saw. "' H-~:~D 1 L 77 L KK:1;NE RL$3(12,1) REL UL,X(1’,Y11).X(2’,Y(2),A(3)9'(3).1$I)E 1 Fux‘-T‘~.3 br16.-, :2 _. _ Hr..' (1.31.152) Nng"4%X(1),V(1’.l(2)p'(2)’l(3),V‘31913.J: 300? FL'.“Y (~13’DF101", -1) 0L”: 0.0 77 cunlxuui LAB- REAL 1 .- . kLaJ «1 (vv1.-).-:=1.:2) ~1 FLKQIT 11215.0) _ LA»- SETHATtA.Ci10,th,KYY,H.IINF,NF,JuS".JG;M,nghE,AN,M*—,J[ND: Ist'NIT IHTP..KK) C‘LL 'klNS‘cT(‘-1Jt~SH°1)'.(JGC19111.1.)55’1’ 9H JTNO1, ‘(1’2NF9"D“9I\~~ 15555117.IIIR,NI‘,.Hhu,Kxx.KYY.H.Y1NF.JG.~.1.J1’.H.N£.J{h3.a‘~) 118 I I h 3 I I \I B l\ S ) ) N E F 3 3 N YA I v I ~ 3 S I O E I I 3 1 ~ I .5 l.\ l.\ I 3 ( 4A x T 1 o. I . Flu F ‘u I \I I ‘5 0 2 7.. l1 I I l\ l‘ . a I .v I V- o G 3 H I I 3 o J I F ) II I .U H l\ I l\ l\ \I 1 3 H C X X 3 I to s u I o ( \l J L '1 D \l J \l P \I F 1 0‘ P 1 N 3 7. I v- E O A K I I ' ’ . .P \I .5“ F 3 H \I II J I. )0 h ‘ P 1 2 .I 1 . 0’ .1 H .u l\ l\ L J .l) T C o x l a .1 IJ " I 3 ‘ I \I p I pp“ 9 ) C 9 3 2 H a 1. Lo 7 3 J A l\ 1 E 2 p- ‘3) Y I. D K 5 Y .I ( T R... x c P I N C 5 OJ 0 O (I . ' ' a“ . I .- 75 ’ .o 3. J T I .1. x I I T I 1 I 000 Olp 3 3. l 3 3 I 2 1| TI DOG» )0 .U / (I ‘ l! 1| .1 l\ I x a 111 ) N)) h ) O V. I a J D .5 2 I .L I I I 1 032 ‘ I ir‘ a I x I I N I I O a E 111 2" 1 .511 13 O. H ) I 1. pk I z 1 3 .U.P 111 0 I 1 TJJ \- 1‘ 7..) H 3 ) 0 E ) .. ( / LU 009 05 3 LK.. P 6: IJ .UIII‘ 3 I J -& I Y I Po... 313 3:? I .116 b . In ‘1‘ IK—J (1 Pd 1‘ I\ z ‘ \l A I I I I I. U). A 1.5. T8 ‘Kh) J) H S o. ) 3 TR 111 CK 1 F2) h 10 (o I. a 1.113 a G L N 1 3 ( NE 111 TA 3 (1 T a) a u ) DITOoIl‘nD. 5‘ .. \l . ofi. I O I‘ ! up DOME I I IHIJ‘ .11 I1. s. 0.50. A70! . IDIR 2 o v ) .20 x 0 in 333 3.. 1 KJ E v(/ 2.<( HH) ll) 3 I 1.4 FT. 5 £111~ wank. O ) ))) CE GP 3 59C 014 .795 0'. OE) Ls’I‘o.‘ (K N IE NEF222222)))))I\I 2 123 I), .IU 0 ‘IP ‘ 0h.“ K. . L IJ3A‘31b V I )5J 06 o o o o 0 03122313 1 LE 1. u o o 03 ) I 3 10 E 1 Li TU) taT((VJ(H.o V7 15 .( E) . 00000C1(((((( X NHH F. 222 l :a Y() R? 00 o hi :AXvAIZLF .IY .. 37 1 .h ..L3333333YYY>1XY 9 .. ((1 F0 3331 2:1 ) 7.41 hi” ))1 ~ )LK LN 001K 001A 10.111/llll. . . . . .0 )L.. LLL D . . . .. 31 h 5 1 C .1K( E 51kt J01. ~ INNNNHNHPV I £1Ev 10 I~)))’)\IIII)))’ 1a))..l “thoh s III)” 0P1. 5310:?) H '(J ‘H 3““ I: :30 .sOUUUUJP. .Vx ZJ 19:151.. 1)(12312323131£Z ( a331VVV NH 123 7 J .v..r1fl0$x n .. 2.. . UE 11¢ r.[33 U)) o 0 TC. .1 2.1.9.4 a/XZ . 2 : U E :1 ((((((((((((( X.n(((.IT.1 0C (((1),b.r. : A: . :3: K .. UTX L L . H1AU t 1 .LJ ~10 UNSSSSSSNPNKIS o DD(T)~011(TXIXYYYYYYIXXX (1555...... 11 JJJON .P) ..)K) ) : )NCXK NE KEN. N11 o a. o: .. 3 ofiKNNNNO..U 933 0 V1 0.1.1N (Ea. : .. .. .. = 8 .. : .. a .. :1 3.5.-11%).!) f-I. ololvlCl‘T N)VKM:JN)\1JK : O ((5.1: : .11. .v.) .‘DE...E:..1::I1I. .0: :01.931.:.ol.07334)))))))))))): 1R: : :123 1.. (((3J( .(N‘1‘f‘!())=3 1r 1.. ...>T1J1‘(l\11 uhHHHHHHH 3H 5.01:“ Avthd £1K123123123123 w o.~..123((( .N .1 .J(J 1 J(JH.nYH TO 11...“:r:= HH1I. .LE....-. . 1.. .1 .LKCEF hr 5 LL0‘ LECL:.C(((((((((((( LZ FLEEJJJ LU FFF 0...? a .J. .Fh L .JL .v LYL A FhréL LLKF LLLF $JEDCDDC..CRCR5$ DRF‘CRDKIFXXXYYYBBBCECL. ANNNIIY CF 11. .bTILLKPILbCKPLCKC HA 1&550091ELUK 1 c 1 1 . H... 5 0 Q ‘5 751 1 7.. 1 UR 1 9 5 1 0 1 0 CI 001 1 D 0 SI 1 0 1 0 :a 0 .a a C 00 3 U 0 . .c A. 0 1 1 Au .3 113 3 3 .a 3 1 u a u u y u u .C 1 5 «L C o 5 .6 :c 94 C. AJ Ca . . r: p . r: 0 p: .U 1 1 7. 2 3 3 ~ ~ : : t c 7 7 O o 9 119 ITILS CONVLCTZCN RELLTEC OUAhT CALCU.AIION CF Tfii C-L U C 30 r. TI [- H to S E N t: 9? I T s o I- 8 AJ L in O o .0 I )2 2 2 0H1 1 E O 61’ o I H O ’anz L 9| I I. hi h I II I. I . av J "2)... I 91 l\ (9)” I N VI ‘IH. “ 1 o 0J1) It I I“) x s u. JXOH o E ll (1)1 o I Y XOJX 3 T 1 O )1. I .‘ 6 I O 0.!) I up. 0 2 ZCOJ J D. .I . (cl .1 l\ U 0 O OF3X I. R G ) €N11 1 F \l "7‘. I . I J .T)\I)) II a 1 .11J‘HH N 0 x l J I I I I .‘h E 1. LLJJJH N L 3) HH1111 E o ”H OOHHHH . \l 1 )).~ .0... C 70‘] )X JHLEL'L 21.‘ 3‘ 11:32: F 1L1 0‘ FF)\I\I\I fly 10?. QTI.ULEJHJH =10 ERk:8IIII N 15511 03. ))JJI‘~ a OISIJS4JI‘1111 .1 11. .J1: ..11~..HHH .l LF::FG LFF 37.1. t. DIJHILHEELELE E .3 o N 1591c .v rb 6 av v r: 0 5 9 0 0 Id 1 1x I \l -I** 11JJ 6 I11 1)61 11. .76 O 8. 2FA1J. .1MHH I1”: :1 . )sc .)0 )J.>JJ.J.U r~z.,7 ‘JJJJ ....J1...1=: U151 UJ F56 .UIJ‘UJJIJJ I 6 1 110 115 D)-2' a a 0 ELL OUTPUY UP “REL NTEALPY F (Kl-13S) 711 ( ¢ E 2: H 712 ALh=SH1/;A 1E 120 125 130 .712)712 7 O 5) & hVEf‘SE ENThéLPY'.‘10. Pn-~t alil‘kLn _ FLnHAT (' ,‘An: 111 I11 CLMTINUL .2! R TillhfiuE HAT' DEULHPOLEE & :NID UPPE . acqce (&(J N01),N9.N6H) JiN DLH’BJ 1~0 6;.) L :0 s A11. CPE 10 15 20 25 30 80 60 75 00 90 95 100 105 ,‘S'"3’CFSO $USRJU7-NE-1H"\SNT (“H.GHH FnggllnggNadghsL’UT9N1791NTK.NH.CO!FG. :531’=1-~. s 90:52.. 2 01:57.2 10:32. 3 x;;1.3 s x~=.32 : 3Pd=1. PK1V: 1000 L-;RI ”SaJH’JI’TOIK19"9CPN REA) 2000 fr 15:‘0F~. 120 ;?0 1x Ph-N! 1001,1i..4..urut.126. 25.1K CH5; P‘P.‘H1 . _ - . _ 1000 FLEH&T('1','INE \0n00n MODEL ran tn:aa.. CQNDU 0110.1v n.; Sikh . 5 'SPEL F [0‘1'0' '7": FoLLauxue an.;La ZONSIANIS n-v; BEEN . 02:10 1N‘.00' 'CENVEdsxon :0 oesaees L35:.uTL = ',£12.5/ 3 ‘0'.‘GAS consinut = - £12.5/ 3 '0',’nE&T er $05101 f5i ungza = ~.£12../ Z 001,002n;=¥v 8r HATER s - 112.51 '0'.'JL~L._Y r 10- = '.£12.5/ r '0'.'FHELZ.hu Pox“: F02 PuiE utTER = 2.512.5/ 0 'O’g'THLhHAL :cwouzrzvzrv or 1:: = '.L12.s/ a ‘0','THE£§LLASDNDULTLJIYY_JF unrza =,- £12.51 1 '0' 'bPEL1F1L aLaz or uAtzi = °.E12.>r}/0 10;: scan;¥('o‘ 'THE FOLLouxus 92:00:! 01:1 n15 :ecu 2:13 1~~r 1 001.05ugiz-uu 2010: 0r Pquucw a ~.L12.s/ g '0-,'.n T “3 natan CONTENT = 0,212.5/ 00'.-unfhtt aa,£ ants: ccyt£~1 s 0,212.5, s '0',¢;u27;~L PnoauL? DENS1IY = 0 212.5: . 5 '0'.' u 314. srsc:r;: Hfizf‘OF 9250001 = '.£12.5/ a ' 0 ' ..~L liciflnL :0uou.7101tv or 91000.! = 'pE12.S/ll) 2000 Foaqa}(tr1c.0) RLTJxfl Can. 9‘0 SUoQDUT-NE PKOP1 (1'CP‘QKOP0. H1, "‘0 Tnxs SUBROUTINE CALCULAYES FiGPE‘IY 01.055 0? a azvzn TEHPERATUEE l&§";°~/PKUFIBZ’:.'HS.O.H 0901": KIgKflpCF-HFJHCNFIC.DS,KS.H;.CFS II Kn K1 K3 H “E §31 38¢ 5 ’ ' ' I gozozo is. o b” o “c ’ EAL» VE,E?rL,E§. fis,fi-,51.v?;.4r.,VFut bfiL..THr.EEL(V‘-ao Sde‘K) gc~;1'g‘cpuo(1.0-1uL» .=s L C 1 Ca . 0EPHEL-(FL n; ta H3 JFh; 1 T ,u: H; 32:. 3‘“5%¢“‘1“§ §£.03.5505é>. “' ' ’ b- .R a We '1- NF F wry: F?r?vf;‘vru1 ’ ' 1’ °' "v " LAL. 'Hht1L1VFJoK;,KJ,KS, LA-. THKEEU1VF99$39KL; K) L3g,_§§é€:s‘f:do~P'1‘1b;1|‘C;fl.9~:’UF:IC::,:1) ALL .5 .' u; H H .1“: at; . ~11) 05L. EnrnA.«LPs'upl,§fi.Iic. fi..u:JtFu$,15fi crastu1-h2>/0.0b£ c LUuIANU’. kEIJ - 'r. H L U~°I~G A M u . g /:g;é9gpicfifibg5: EUEE. gga‘gn§3§§§ 885505$§U;99L;33‘0§’7 . v - ; N" ;; ; = o .5 1 ’ 3 :ggiagszEOEch ;7:t3‘0; 3cL33§ g°-,£:§E§? I.“.’ q , _~n 3 Hu.£ J.lh IE1 «T F ’-"i = ’ . 5 33;o':PL.IF;. near St 50L105°= ~?czifs5}/: .212 D 10 m Q" 1C ‘7‘ \H 1224 01 E 0EN;3-(1PJ,1IS.DI.DS) spa 10~ 1111. N (“TV ! MP 5/1P0-1u»I0u1/11.0-I~;11 t -0npu 80-" I" E PZJ<</OS) ['17H§-(CP-,L3H0H50I'JI”AO';OU:“§ITO”) NU! I“? C PTTOG’DM 1'1"? «11:1:- . ~11 .0 '0 o 1 r0 5:019:11: A: -. a ' $u;i3UT.NI 0Een535(L1.b;,Hs.ro.ns.urdc.zuc.t.«(1121 RL“. Hui ‘HL fl. "fl nr=1s.0i§3u' ' ”TINIW‘MX "f'. 'fif'fi‘!‘ "‘0 n gun-1'1" "f‘l Al. 11".)- 2OC’tI‘ '0 7090.0)9H.'(HSOCPU'(YO~J.O))onz';P;'(IO~G.a1-JF~:'H: APPENDIX B PHYSICAL PROPERTIES OF GROUND BEEF, COMPUTER INPUT PARAMETERS AND VARIABLES, AND EXPERIMENTAL DATA 125 meo mm.pmoa AANn+nm+mmvspv o\H HOGOE chflowmmmuu mo OEDHO> “mm m m50 mv.nmm~ n a 6‘ 2.99: m) u 1602 $03330 mo 9:38, \\\\‘1 .mE9H0> >3 cmofl>fl© unmfim3 Scum cwcflmuno was xuflmcmon .ucmEummuu 50mm uwm mcoflumofiamou wwucu Eoum ovmuw>e omo.H boo.H m.mmaa m.thH mo.ow o.a+o.om1 m.HH N>B mmo.H who.a m.omHH >.mnHH H5.0> m.HHo.omn ¢.m H>B HMUHONwmmMB omo.H vmo.H n.mmmm m.moom mm.mo m.a~o.om1 m.mH m>m vmm.o vao.H H.mm>~ m.momm mm.oc o.mmo.om1 m.HH m>m moo.H wvo.H v.mmmm m.mmmm om.on o.m+o.om1 v.n H>m Hmoflumflflam cmwoum cmmouwca camoum cmNouwco w Uo .mpsu m\E .ucmucou ImmeEme >DHUOHm> ucoEummuB m .u: mEo\m-fi>MHMJmQ mama ummz ousumwoz ufl< 0fl< mmucoEumwue HmucmEHuwam msoflum> um mmmm ccsouu mo xuflmcoo Ucw ucmucou musumHOZ11.H.m mqde 126 127 TABLE B.2.--Physical Properties of Ground Beef Thermal Specific Density,b Initial Conductivity,a Heatb k /m3 Freezing Point,a W/m K J/Kg K 9 °c 0.4 3430 1060 -l.5 aValues obtained from literature for mince meat (Screnfors, 1974). bValues determined by experiment. 128 TABLE B.3.--Input Parameters for Computer Program 1. Product Properties Parameters IWCa (ground beef) UFWC = 0.01 IPD = 1060 kg/m3 ICP = 3430 J/kg K IK = 0.4 W/m K TIF = -l.5°C 2. Physical Properties Parameters DI = 920 kg/m3 gislce, Water, and DW = 1000 kg/m3 KPI = 2.32 W/m K KPW = 0.55 W/m K Lw = 335.103J/kg CPW = 4180 J/kg k TWF = 0°C R1 = 8314 J/kg-mole K 3. Freezing Medium Variables v, m/s h, W/mzK Properties 7.4 61.7 11.3 70.6 15.2 142.5 T = -l4°C, —18°C, -22°C 4. Finite Element Parametersb NP:number of nodal Properties temperature NE:number of element NBW:number of bandwith NIT:number of itera- tions Variable Dthime step = l min— ute and 3 min- utes aSee Table B.1. bSupplied by the GRID program. .... ::_.=.>;: cuccccbu "m uctsaeh;;Ecu no ma~c> ccmhc>< "I. u .Lo:_ zgzec»:;ELu gusto»; unqu.c~0 .33..»933 b: ....e 3:0: .5... .2: 3:... 2c To 32.: 3 .85... .a:;E~mo»u :11; ha; r::_.c;.~;:t tangy EOLE accuc>< a 129 w I m I m I w I m I m 2 m I x Z a I 0.0 0.0NI 0.0 0.0NI 0.0 0.0NI 0.0 0.0NI mbwlflfiwww c— 0.0 0.0NI 0.0 0.0nl c.~ 0.0dl c.~ m..~l n.~ n.— I v. o.c o.c~.. ~.~ o.v~1 En m4 I m.~ NJ ~.~ v.2 v wim— 0.0 0.0NI 0.0 0.0NI 0.0 0.0NI v.0 0.0—I N.~ 0.3—I :— o.o 0.0N1 a; CALI m4 n.m~1 h.~ v.0 I c.— we I 3. 6.: m.o~1 m4 ~.v~1 3..” A; 1 v.~ o.~ m._ 9.». v n4. 0.0 0.0NI 0.0 0.0NI 0.0 0.0NI 0.0 0.0NI 5.0 o.o~I e._ ~.C_1 m— 0.0 0.0NI 0.0 0.0u1 0.0 O.0~1 N.~ 0.0—1 m.— v.-l o.~ c.m I @— 0.0 ¢.G~1 ~.0 G.0~1 m.~ 0.v~1 0.N 0.~ I Q.N 0.~ :._ ~.v_ v v.h an‘ONchuP o.o 04L: 0.0 0.4:: 0.0 m.o~1 3.0 08:1 v.0 ~.o~1 NJ 04:1 v.~ 957. c.~ 9.21 mm 0.0 0.c~1 0.0 N.o~1 0.0 h.0~1 0.0 o.h~1 ®.N h.m~1 v.N 0.0—1 m.N m.m I m.~ m.~ :~ 0.0 0.0-1 0.0 h.m~1 G.~ m.m~1 c.m ~.o I ~.m Q.N I —.N m.~ I m.~ c.m 0.~ ©.m— m N.m~ 0.0 0.0N1 0.0 0.®~1 0.0 O.O~I 0.0 0.0a1 n.0 0.0nl m.~ d.0~1 0.N @.b~l ~.N v.m~1 mm 0.0 o.c~1 0.0 Q.0~I N.0 m.m~1 0.0 m.m~1 0.N «.mfil 0.~ N.@ I v.~ v.m I m.~ 0.N ¢~ ._.0 $57 ~.~ ¢.o~1 o.~ o.-1 in c.o~1 ~.n m4 1 m.~ m.o 1 V4. o.v v.~ Q: r I: 0.0 0.0NI 0.0 o.o~1 0.0 h.a~1 0.0 m.¢—1 0.0 m.0_1 0.0 0.0»1 ~.0 n.m~1 €.~ v.h~I €.N 0.m~1 FM 0.0 o.a~1 0.0 0.3~I 0.0 N.9~1 N.0 0.0—1 h._ m.v~1 m.~ m.-1 m.~ N.h 1 v.~ 0.N I N.“ c.p 0~ 0.0 0.0—1 ~.0 h...~n 0.0 m..:1 m._ m.~—I A... 0.? I h.~ ~.— 1 m.N m.c N.~ 0.x. NI— _..: n v.h AQUwhnmmqfiw a a h c r v m N d 0;;E:z m\E .> >huvsczc ,1111111111 I! I- 1- 1.11 I--I a .ct:2 u:;E.1¢.F aware»; .1: SEC. 1mm;;:_; 7:.N;;_; tzdnsc r.:;£0u;:_ ;E_% _:::1;:: .1 .;o. 1.1: 1t:acLL;E;& 6:2»:uchI.v.m m4m um mcflNmoum ummHQIHflw mcfluso moon ocsoum HMUHONmmmuu co musumquEmu mommusm ch H.o mucous 143 ON... H .H. X UCM mmmuflm HHM .HOnm m.m A U m.wa I uI u E x pcm mopflm Ham Mom m.H A u Hw.o + uv.vm I u x bv.H I my x om.m n B m opflm vo.o + uo~.o + u x mn.m I Nu x Hm.H n B m.mH v x v Ha ma.m I uoa.om I u x mm.o I m» x m>.H n 9 Ha v x x o "N moflm mm.m I unm.m I u x m¢.v I my x mm.m n B H opflm m.H v u v o om- n was ma Havoc 00m o A 0 ma n B x pcm mooflm Ham Mom 0 n u 06 .musumquEmB SO .x u: .wEflB m\E m.b u > um coflumooq cam mEHB mo cofluocsm m mm musumquEme momwusm mcflnfiuommo mambo: HmoflumemcumZII.H.o mqm