- Email: [email protected]
Transient thermal behaviour of phase-change processes in solid foods with variable thermal properties
Journal of Food Engineering 54 (2002) 331–336 www.elsevier.com/locate/jfoodeng
Transient thermal behaviour of phase-change processes in solid foods with variable thermal properties F. Alhama *, C.F. Gonz alez Fern andez Department of Applied Physics, Technical University of Cartagena, Campus Muralla del Mar, 30203 Cartagena, Spain Received 11 August 2000; received in revised form 29 June 2001; accepted 19 November 2001
Abstract In this work we study the transient temperature fields and heat fluxes for simple or multilayer foods with planar, cylindrical or spherical geometry in one-phase or two-phase (phase-change) processes. The temperature dependencies of the thermal properties conductivity and/or specific heat are assumed to be specified as arbitrary and continuous functions. These dependencies are normally very pronounced in these products in the phase-change temperature interval. Interpolation tables or piece-wise functions, which are widely used in foods, may also be considered. Simple boundary conditions such as isothermal, convective and radiative, or any kind of compatible combination of lineal and/or non-lineal conditions may also be assumed with no special requirements. The problem, expressed by way of only one governing equation for both phases is solved numerically by the Network Simulation Method. This method, which has been recently applied to a great variety of non-lineal transport problems, only requires the use of finite-difference schemes for the spatial variable. The time remains as a continuous variable. Application to the freezing of white fish is presented. Ó 2002 Elsevier Science Ltd. All rights reserved.
1. Introduction As regards the freezing or unfreezing processes in foods, the dominant mechanism of heat transfer is the conduction. The solid–liquid phase change occurs progressively in a finite and (approximately) well-defined interval of temperatures. As a consequence of this, the two phases are not clearly separated and both co-exist in a finite portion of the food at a given instant. In this phase-change temperature interval the thermal characteristics – conductivity (k) and specific heat (ce ) may be well approximated as arbitrary and continuous or piecewise functions dependent on temperature, T. Generally, the specific heat exhibits a sharp peak at the initial phase-change temperature, whereas the conductivity has a more regular variation. Expressions for ce ðT Þ and k(T) are frequently found in a tabulated form in the research literature (Polley, Snyder, & Kotnour, 1980; Rao & Rizvi, 1986; Sweat, 1974). Both the mentioned kinds of dependencies together with the boundary lineal and/or non-lineal conditions, such as natural convection or radiation, make these kinds of problems strongly non-lineal and prevent an *
Corresponding author. Fax: +968-32-53-37. E-mail address: [email protected] (F. Alhama).
analytical solution of the conduction equation. Some numerical methods, based generally on finite differences and finite element techniques, allow a solution by using an only governing equation to be applicable for both phases. Each one of these methods has their inherent advantages and disadvantages. On one hand, the use of temperature as dependent variable (named apparent specific heat formulation) is limited by the stability criteria since the time increment must be high enough to prevent oscillations caused by the sharp peak of the specific heat (DeBaerdemaeker, Singh, & Segerlind, 1977; Cleland & Earle, 1979; Lin, 1972). On the other hand, the use of the enthalpy variable (named enthalpy formulation) requires to know the enthalpy dependencies of both the conductivity and the temperature (Mannapperuna & Singh, 1988; Voller & Cross, 1985). Besides the interest in knowing the transitory time of a particular thermal engineering process, the solution to the transient thermal fields ensures not only that the critical points reach the required temperature but also that it determines the design parameters for the energysaving systems within the agro-alimentary industry and predicts the quality loss in stored products due to temperature fluctuations. Numerical solution to these phase-change processes obtained by the Network Simulation Method (NSM) is
0260-8774/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 2 6 0 - 8 7 7 4 ( 0 1 ) 0 0 2 1 9 - 9
332
F. Alhama, C.F. Gonzalez Fernandez / Journal of Food Engineering 54 (2002) 331–336
Nomenclature specific heat ðJ kg1 K1 Þ capacitor voltage-control voltage source control current source heat transfer coefficient ðW m2 K1 Þ enthalpy ðJ kg1 Þ heat flux ðW m2 Þ thermal conductivity ðW m1 K1 Þ characteristic length (also thickness of the plate) m and n 1; 2; 3; . . . N number of volume elements Nu Nusselt number S surface (m2 ) t time (s) T temperature (K) Ta surrounding fluid temperature T0 initial temperature x space co-ordinate (m)
ce C E G h H j k L
Greek symbols a; b natural number
presented in this paper, including the simpler processes of heating and cooling without phase-change (one-phase problem). 1-D, planar, cylindrical and spherical geometries are assumed and an application to the freezing of white fish is presented. Common simple boundary conditions such as isothermal, convection and radiation as well as other kinds, including a compatible combination of lineal and/or non-lineal conditions may also be assumed by the numerical method with no special requirements. In relation to the use of non-lineal natural convection, standard non-lineal formulas for the value of the heat transfer coefficient or Nusselt number may be used. The NSM is a numerical technique for the solution of lineal and non-lineal problems, which are defined by a mathematical model. It has already been applied to lineal and non-lineal problems both in heat transfer (Alhama, L opez S anchez, & Gonz alez-Fernandez, 1997; Gonz alez-Fern andez, Alhama, L opez Sanchez, & Horno, 1998) and in other areas as well (Horno, Garcıa Hern andez, & Gonz alez-Fern andez, 1993). As in the Lines Method, only the spatial variable is discretised and the resulting differential equation with time a continuous variable is used for the design of the network model. Assuming that the electrical variables voltage and current are equivalent to the heat variables temperature and heat flux, respectively, a network model for each volume element is designed in such a way that its electrical equations are formally equivalent to the spatial discretised equations. The whole network model, in-
q Dx Dce r C
density ðkg m3 Þ thickness of the volume element (m) defined in Eq. (7) gradient operator volume of the volume element
Subscripts 1; 2; 3 associated to specific devices in the volume element C defined in Eq. (8a) ce associated to specific heat i index (and also central point) of the volume element i D; i þ D end points of the volume element k associated to conductivity mean mean value s associated to the surface a 1; 2; . . . ; n b 1; 2; . . . ; m c defined in Eq. (8b)
cluding the devices associated with the boundary conditions, is solved by the commercial computer code PSPICE (PSPICE, 1994). Among the advantages of the NSM it is important to mention: (i) very few and simple devices that are necessary to simulate complicated non-lineal expressions; as a consequence, a limited number of programming rules are used for writing the circuit file. (ii) non-linearities associated either with the boundary conditions, with the temperature dependencies of the thermal characteristics or with certain special phase-change conditions, are implemented easily in the model by means of simple electrical devices or auxiliary circuits (Alhama, 1999). (iii) an only network is necessary to analyse the whole process, including the transitory and stationary responses. (iv) the user does not need to manipulate the discretised equations nor to have in mind considerations about the convergence or the stability of the solution; so that no mathematical manipulations are necessary by the user because the software selected to solve the circuits does this work (Nagel, 1975). (v) NSM permits direct visualisation and evaluation of the main variables of the transport phenomena heat fluxes and temperatures as well as other physical aspects associated with the process. The circuit resolution software selected in this work, PSPICE, has an extensive library which deals with non-lineal electric devices that are suitable for the solution of this kind of problems.
F. Alhama, C.F. Gonzalez Fernandez / Journal of Food Engineering 54 (2002) 331–336
The errors of the NSM are studied by Alhama (Alhama, 1999) for both lineal and non-lineal problems including the phase change in pure materials. These errors are acceptable in engineering using a number of volume elements N P 50 for 1-D geometry; average relative deviations between theoretical and simulated temperatures are reduced by about 0.1–1% depending on N.
333
temperature. Other kinds of boundary conditions can be easily assumed. Equations (4) and (5) summarise the discrete values of k and ce into the temperature interval to which the food is subjected in the process. All kinds of real dependencies may be approximated by these equations (interpolation tables) with negligible error by taking high natural values for a and b. For different temperatures from those ones on the table, values for k and ce are found out by simple lineal interpolation.
2. The mathematical model 3. The network model A food product of 1-D, planar, cylindrical or spherical geometry whose initial temperature is T0 is subjected to a thermal phase-change process of freezing or unfreezing, under the application of defined boundary conditions. Natural and forced convection as well as radiation, or isothermal boundary conditions, may be applied to the outer surface of the food. Variations of the thermal conductivity and specific heat or enthalpy in relation to the temperature are specified in a tabulated form by n and m values, respectively, which are associated to n and m discrete temperatures inside the work temperature interval. Continuous arbitrary functions may also be assumed for these dependencies as simpler cases. Making use of Fourier law, the governing partial differential equation for both phases is obtained by performing a heat balance over a small region (volume element) of the food, the form of this region depends on the co-ordinate system. The complete problem statement can be written as qce oT =ot ¼ rðkrT Þ; hðTs Ta Þ ¼ koT =oxjs ; Tðt¼0Þ ¼ T0 ; kðTa Þ ¼ ka ; ce ðTb Þ ¼ ce;b ;
1 < b < m;
½SiþD kiþD ðTi TiþD Þ=ðDxi =2Þ ;
ð6Þ
ð4Þ ð5Þ
ce;i; ¼ ce;i;mean þ Dce;i :
ð2Þ ð3Þ
1 < a < n;
Ci qce;i ðdTi =dtÞ ¼ ½SiD kiD ðTiD Ti Þ=ðDxi =2Þ
where SiD and SiþD are the boundary spherical surfaces of the crown and C the volume. Two values for the conductivity in the volume element are defined; one on the left, kiD , whose value obtained from Eq. (4) is going to be a function of the temperature at that point, TiD ; other one on the right, kiþD (the same but using the temperature TiþD ). This is the conductivity ki D ¼ kðTi D Þ that is evaluated by means of the table given in Eq. (4) by simple lineal interpolation. On the other hand the specific heat for the cell, ce;i , is evaluated according to the temperature in the centre of the cell, ce;i ¼ ce;i ðTi Þ, using Eq. (5). To design the network we use the equivalence ‘‘temperature voltage’’ and ‘‘heat flux electric current’’. Now in order to implement the initial conditions, for which a condenser is necessary, the values of the specific heat are separated into two terms in the form
ð1Þ h ¼ hðTs Þ;
Spherical co-ordinates are used for the design of the network model. For 1-D geometry, Eq. (1) applied to an elemental volume of spherical crown form (which is defined by a thickness Dxi and radii xiD ; xi and xiþD , see nomenclature in Fig. 1) may be written as
where Eq. (1) is the governing equation in vector notation, which is applied to the whole medium; q is the density, k is the thermal conductivity, ce is the specific heat and r is the gradient operator. Eq. (2) is the natural or forced convective boundary condition. As regards the heat transfer coefficient, h, or its dimensionless form Nu ¼ hL=k, with L the characteristic length of the product, standard expressions dependent on temperature for natural convection (such as Yuge, 1960 or Churchill, 1983), or constant values in forced convection may be used; Ta is the temperature of the surrounding fluid. In the end, Eq. (3) is the initial condition. The set of equations (1)–(3) makes up the so-called apparent specific heat formulation. The use of dimensionless form for these equations is not possible because of the variations of the thermal properties with the
ð7Þ
Defining the currents ji;C ¼ Ci qce;i;mean ðdTi =dtÞ ¼ Ci ðdTi =dtÞ;
ð8aÞ
ji;c ¼ Ci qðDce;i ÞðdTi =dtÞ;
ð8bÞ
jiD ¼ SiD kiD ðTiD Ti Þ=ðDxi =2Þ;
ð8cÞ
jiþD ¼ SiþD kiþD ðTi TiþD Þ=ðDxi =2Þ
ð8dÞ
Eq. (6) can be written in the form ji;C þ ji;c ¼ jiD jiþD ; which is a kind of Kirchhoff’s law for these currents. The network model for the volume element or cell is now designed. It is integrated by only four devices and three simple special sources, named voltage-control voltagesource (E), which provide the values of kiD ; kiþD and ce;i . These special sources are Ek;1 ; Ek;2 ; Ece , Fig. 2. In this
334
F. Alhama, C.F. Gonzalez Fernandez / Journal of Food Engineering 54 (2002) 331–336
Fig. 1. Nomenclature and geometry of the volume element.
way, the output voltage of Ek;1 is the value of kiD associated to input temperature TiD , according to Eq. (4). The same is applied to the auxiliary source Ek;2 , which provides kiþD , and to Ece , which provides the value of ce;i associated to the temperature Ti . Eq. (8a) defines a condenser of capacity Ci ¼ Ci qce;i;mean and Eqs. (8b), (8c) and (8d) are implemented by a non-lineal source named control current generator, G. This active device is a current source whose output is evaluated as an arbitrary function of one or more voltages or currents at different points of the network. In particular, Gi;1 is controlled by the voltages Ei;k;1 ; TiD and Ti , according to the expression ‘‘ðSiD 2=Dxi Þ Ek;1 ðTiD Ti Þ’’ and Gi;2 by the voltages Ei;k;2 ; TiþD and Ti , in the form ‘‘ðSiþD 2=Dxi ÞEk;2 ðTi TiþD Þ’’. In the end, Gi;3 is controlled by the voltage Ei;c and the current of Ci according to the function ‘‘ji;C ððEi;ce ce;i;mean Þ= ce;i;mean Þ’’.
Fig. 2 shows all the devices of the network model; variables between brackets and in cursive mean control variables of the adjacent non-lineal sources. Voltages and current source must be connected with the adequate polarity in order to satisfy Eq. (6). Geometrical parameters in each cell (radii xiD ; xi and xiþD , surfaces SiD and SiþD and volume, Ci ) are different for the cylindrical and spherical co-ordinates and constant values for planar co-ordinates. Once the network of the volume element has been designed a number of them (N) are connected in series to conform the whole region. For N P 50 errors are reduced to values less than 0.5% in 1-D geometries (Alhama, 1999). Boundary condition is implemented by a voltagecontrolled current-source in the cases of convection and/or radiation, and by a constant voltage source in the case of isothermal condition. As regards the initial condition, this is applied to the network by fixing the initial voltage of the condenser to the value of T0 . Once the whole network is completed, it must be written in the way of a program that could be read by the computer code Pspice. Only a few rules are necessary for this purpose since only three different devices are needed in the network. It is also possible to change the parameters of the problem by means of simple programming so that this can be adapted to different geometrical or thermal values. The solution is obtained for all the temperatures and fluxes, at the time, by running the program in Pspice. A few seconds of computation time are needed using a PC (Pentium II, 500).
4. Numerical results and discussion
Fig. 2. Network model of the volume element.
An application is made to the process of freezing of white fish by forced convection. For the storage and
F. Alhama, C.F. Gonzalez Fernandez / Journal of Food Engineering 54 (2002) 331–336
preservation of the white fish, temperatures about 20 °C must be reached; the initial temperatures being about 5 °C. With this aim the fish is taken into the freezing unit multi-plate with temperature operation of 40 °C, undergoing forced convection on both sides of the plate. Once the critical point of the fish reaches the required stored value (20 °C) the food is taken to the frozen-food vessel. To simulate this problem a planar 1-D geometry is used, being only necessarily half a plate due to the symmetry of the problem. The parameters are: L ðthickness of the slabÞ ¼ 0:1 m; q ¼ 1054 kg m3 ; T0 ¼ 5 °C; Ta ¼ 40 °C; N ¼ 40. The temperature dependencies of k; ce and enthalpy (H), taken from the literature (FAO, 1985), are shown in Tables 1 and 2. As regards the heat transfer coefficient typical values around h ¼ 600 W m2 K1 , reported by Gruda (Gruda & Potolsky, 1974), have been adopted. Fig. 3 shows the transient temperature at five regular positions, x ¼ 0, L=8; L=4; 3L=8 and L=2, taking x ¼ 0 at the surface of the plate. The time required to reach 20 °C at the centre of the plate is about 2h. At this moment the temperature on the surface of the fish is 38 °C. The effect of the near regions of the surface which is undergoing the phase change over the most internal regions is observed in the figure, as well as the
335
Fig. 3. Transient temperatures at five regular positions.
Table 1 Thermal conductivity of the white fish Temperature (K)
Conductivity (W m1 K1 )
278 273 272 243
0.43031 0.43031 1.30256 1.87243 Fig. 4. (a) Instantaneous heat flux. (b) Integrated heat flux.
Table 2 Specific heat and enthalpy of the white fish Temperature (K)
Enthalpy (J kg1 )
Specific heat (W m1 K1 )
233 253 263 265 267 269 270 271 272 273 275 277 279 283 293 313
0 41985.6 74217.8 83678.1 96319.9 117417 136882 176482 297876 322992 330275 337601 344885 359535 396205 469878
1841.84 2595.32 4227.86 5316.22 7744.1 15111.4 26539.2 65636.5 102724 4144.14 3641.82 3641.82 3641.82 3683.68 3683.68 3683.68
increase of the time required for freezing the same piece of fish as process advances. Instantaneous heat flux crossing the surface is shown in Fig. 4(a). An inflexion point is observed about t ¼ 2h in the curve. The instantaneous heat flux maintains a high value (almost constant) around t ¼ 2h, in which all the fish has changed the phase. A rapid decreasing takes place after this time due to the values of k and ce in the solid state. Integrated heat flux as a function of time, shown in Fig. 4(b), may be compared with the enthalpy step between two temperatures, Table 2. In this way and for the complete process, the enthalpy step between 5 and 40 °C is 341.24 kJ kg1 , nearly the same value as that of the total integrated flux in the simulation 342.43 (error ¼ 0.34%). When the centre of the fish reaches 20 °C the total integrated flux is 316.46; a higher value
336
F. Alhama, C.F. Gonzalez Fernandez / Journal of Food Engineering 54 (2002) 331–336
as we could expect than the enthalpy step between 5 and 20 °C because most of the fish is below 20 °C. Nevertheless the last difference is very small (5.74%) because the heat eliminated from the fish below 20 °C is negligible compared to the latent heat of phase change needed for the complete freezing.
5. Conclusions The network simulation method has been applied for the first time, to the authors’ knowledge, in the way of apparent specific heat formulation to solve the transient problem of freezing foods, whichever be the boundary conditions. The application of the method is simple since very few electrical devices are needed for the network design of the volume element. This, in turn, requires the use of very few programming rules to make the file to be read by the circuit resolution code. The user does not require mathematical manipulations of the finite-difference equations since the computer code specially designed to this end does this work. The only mathematical step is to convert the governing equation into a spatial finite-difference equation. An application is made to the freezing of white fish under forced convection, a severe non-lineal problem due to the strongly temperature dependencies of the thermal properties. Specific heat has a sharply peak near the initial freezing front. Any kind of dependencies for k and ce , and boundary conditions (including natural convection and radiation) may be considered without special requirements. Errors produced by comparing the integrated heat flux with the enthalpy step for two given temperatures are negligible. Computer time for this application is 10–12 s in a PC (Pentium II, 500). References Alhama, F. (1999). Transient thermal responses in non-lineal heat conduction processes. Ph.D. Thesis. University of Murcia, Spain, 234 pp.
Alhama, F., L opez Sanchez, J., & Gonzalez-Fernandez, C. F. (1997). Heat conduction through a multilayered wall with variable boundary conditions. Energy, 22, 797–803. DeBaerdemaeker, J., Singh, R. P., & Segerlind, L. J. (1977). Uses of finite element method in food processing. Journal of Food Engineering, 1, 37. Churchill, S. W. (1983). Free convection around inmerser bodies. In E. U. Schl€ under (Ed.-in-Chief), Heat exchange design handbook, Section 2.5.7. New York: Hemisphere Publishing. Cleland, A. C., & Earle, R. L. (1979). A comparison of methods for predicting the frezzing time of cylindrical and spherical foodstuffs. Journal of Food Science, 44, 964. FAO (1985). Freezing and refrigerates storage in fisheries. Fisheries Technical Paper 340 (ISSN 0429-9345). Gonzalez-Fernandez, C. F., Alhama, F., L opez Sanchez, J. F., & Horno, J. (1998). Application of the network method to heat conduction processes with polynomial and potential-exponentially varying thermal properties. Numerical Heat Transfer A, Applications, 33, 549–559. Gruda, Z., & Potolsky, J. (1974). Zarazanie zywnosci. Wydawnictwa Naukowo-Techniczne Ed. Warszawa. Horno, J., Garcıa Hernandez, M. T., & Gonzalez-Fernandez, C. F. (1993). Digital simulation of electrochemical processes by network approach. Journal of Electroanalytical Chemistry, 352, 83–97. Lin, S. (1972). Similarity solution of freezing or melting process with variable thermal properties. Transactions CSME, 1, 51. Mannapperuna, J. D., & Singh, R. P. (1988). Prediction of freezing and thawing times of foods using a numerical method based on enthalpy formulation. Journal of Food Science, 53(2), 626– 630. Nagel, L. W. (1975). SPICE2: A computer program to simulate semiconductor circuits. Memo. No. UCB/ERL M520. Electronic Research Laboratory, Univ. de California, Berkeley, CA 94720. Polley, S. L., Snyder, O. P., & Kotnour, P. (1980). A compilation of thermal properties of foods. Food Technology, 76–94. PSPICE 6.0 (1994). Microsim Corporation, 20 Fairbanks, Irvine, California 92718. Rao, M. A., & Rizvi, S. S. H. (1986). Engineering properties of foods. Food-Analysis. Series: Food science and technology. New York: Marcel Dekker. Sweat, V. E. (1974). Experimental values of thermal conductivity of selected fruits and vegetables. Journal of Food Science (39), 1080. Voller, V., & Cross, M. (1985). Application of control volume enthalpy methods in the solution of Stefan problems. In Lewis, Morgan, & Johnson (Eds.), Computational techniques in heat transfer. Swansea, UK: Pineeridge Press. Yuge, T. (1960). Experiments on heat transfer from spheres including combined natural and forced convection. Journal of Heat Transfer, 82, 214–220.
Transient thermal behaviour of phase-change processes in solid foods with variable thermal properties F. Alhama *, C.F. Gonz alez Fern andez Department of Applied Physics, Technical University of Cartagena, Campus Muralla del Mar, 30203 Cartagena, Spain Received 11 August 2000; received in revised form 29 June 2001; accepted 19 November 2001
Abstract In this work we study the transient temperature fields and heat fluxes for simple or multilayer foods with planar, cylindrical or spherical geometry in one-phase or two-phase (phase-change) processes. The temperature dependencies of the thermal properties conductivity and/or specific heat are assumed to be specified as arbitrary and continuous functions. These dependencies are normally very pronounced in these products in the phase-change temperature interval. Interpolation tables or piece-wise functions, which are widely used in foods, may also be considered. Simple boundary conditions such as isothermal, convective and radiative, or any kind of compatible combination of lineal and/or non-lineal conditions may also be assumed with no special requirements. The problem, expressed by way of only one governing equation for both phases is solved numerically by the Network Simulation Method. This method, which has been recently applied to a great variety of non-lineal transport problems, only requires the use of finite-difference schemes for the spatial variable. The time remains as a continuous variable. Application to the freezing of white fish is presented. Ó 2002 Elsevier Science Ltd. All rights reserved.
1. Introduction As regards the freezing or unfreezing processes in foods, the dominant mechanism of heat transfer is the conduction. The solid–liquid phase change occurs progressively in a finite and (approximately) well-defined interval of temperatures. As a consequence of this, the two phases are not clearly separated and both co-exist in a finite portion of the food at a given instant. In this phase-change temperature interval the thermal characteristics – conductivity (k) and specific heat (ce ) may be well approximated as arbitrary and continuous or piecewise functions dependent on temperature, T. Generally, the specific heat exhibits a sharp peak at the initial phase-change temperature, whereas the conductivity has a more regular variation. Expressions for ce ðT Þ and k(T) are frequently found in a tabulated form in the research literature (Polley, Snyder, & Kotnour, 1980; Rao & Rizvi, 1986; Sweat, 1974). Both the mentioned kinds of dependencies together with the boundary lineal and/or non-lineal conditions, such as natural convection or radiation, make these kinds of problems strongly non-lineal and prevent an *
Corresponding author. Fax: +968-32-53-37. E-mail address: [email protected] (F. Alhama).
analytical solution of the conduction equation. Some numerical methods, based generally on finite differences and finite element techniques, allow a solution by using an only governing equation to be applicable for both phases. Each one of these methods has their inherent advantages and disadvantages. On one hand, the use of temperature as dependent variable (named apparent specific heat formulation) is limited by the stability criteria since the time increment must be high enough to prevent oscillations caused by the sharp peak of the specific heat (DeBaerdemaeker, Singh, & Segerlind, 1977; Cleland & Earle, 1979; Lin, 1972). On the other hand, the use of the enthalpy variable (named enthalpy formulation) requires to know the enthalpy dependencies of both the conductivity and the temperature (Mannapperuna & Singh, 1988; Voller & Cross, 1985). Besides the interest in knowing the transitory time of a particular thermal engineering process, the solution to the transient thermal fields ensures not only that the critical points reach the required temperature but also that it determines the design parameters for the energysaving systems within the agro-alimentary industry and predicts the quality loss in stored products due to temperature fluctuations. Numerical solution to these phase-change processes obtained by the Network Simulation Method (NSM) is
0260-8774/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 2 6 0 - 8 7 7 4 ( 0 1 ) 0 0 2 1 9 - 9
332
F. Alhama, C.F. Gonzalez Fernandez / Journal of Food Engineering 54 (2002) 331–336
Nomenclature specific heat ðJ kg1 K1 Þ capacitor voltage-control voltage source control current source heat transfer coefficient ðW m2 K1 Þ enthalpy ðJ kg1 Þ heat flux ðW m2 Þ thermal conductivity ðW m1 K1 Þ characteristic length (also thickness of the plate) m and n 1; 2; 3; . . . N number of volume elements Nu Nusselt number S surface (m2 ) t time (s) T temperature (K) Ta surrounding fluid temperature T0 initial temperature x space co-ordinate (m)
ce C E G h H j k L
Greek symbols a; b natural number
presented in this paper, including the simpler processes of heating and cooling without phase-change (one-phase problem). 1-D, planar, cylindrical and spherical geometries are assumed and an application to the freezing of white fish is presented. Common simple boundary conditions such as isothermal, convection and radiation as well as other kinds, including a compatible combination of lineal and/or non-lineal conditions may also be assumed by the numerical method with no special requirements. In relation to the use of non-lineal natural convection, standard non-lineal formulas for the value of the heat transfer coefficient or Nusselt number may be used. The NSM is a numerical technique for the solution of lineal and non-lineal problems, which are defined by a mathematical model. It has already been applied to lineal and non-lineal problems both in heat transfer (Alhama, L opez S anchez, & Gonz alez-Fernandez, 1997; Gonz alez-Fern andez, Alhama, L opez Sanchez, & Horno, 1998) and in other areas as well (Horno, Garcıa Hern andez, & Gonz alez-Fern andez, 1993). As in the Lines Method, only the spatial variable is discretised and the resulting differential equation with time a continuous variable is used for the design of the network model. Assuming that the electrical variables voltage and current are equivalent to the heat variables temperature and heat flux, respectively, a network model for each volume element is designed in such a way that its electrical equations are formally equivalent to the spatial discretised equations. The whole network model, in-
q Dx Dce r C
density ðkg m3 Þ thickness of the volume element (m) defined in Eq. (7) gradient operator volume of the volume element
Subscripts 1; 2; 3 associated to specific devices in the volume element C defined in Eq. (8a) ce associated to specific heat i index (and also central point) of the volume element i D; i þ D end points of the volume element k associated to conductivity mean mean value s associated to the surface a 1; 2; . . . ; n b 1; 2; . . . ; m c defined in Eq. (8b)
cluding the devices associated with the boundary conditions, is solved by the commercial computer code PSPICE (PSPICE, 1994). Among the advantages of the NSM it is important to mention: (i) very few and simple devices that are necessary to simulate complicated non-lineal expressions; as a consequence, a limited number of programming rules are used for writing the circuit file. (ii) non-linearities associated either with the boundary conditions, with the temperature dependencies of the thermal characteristics or with certain special phase-change conditions, are implemented easily in the model by means of simple electrical devices or auxiliary circuits (Alhama, 1999). (iii) an only network is necessary to analyse the whole process, including the transitory and stationary responses. (iv) the user does not need to manipulate the discretised equations nor to have in mind considerations about the convergence or the stability of the solution; so that no mathematical manipulations are necessary by the user because the software selected to solve the circuits does this work (Nagel, 1975). (v) NSM permits direct visualisation and evaluation of the main variables of the transport phenomena heat fluxes and temperatures as well as other physical aspects associated with the process. The circuit resolution software selected in this work, PSPICE, has an extensive library which deals with non-lineal electric devices that are suitable for the solution of this kind of problems.
F. Alhama, C.F. Gonzalez Fernandez / Journal of Food Engineering 54 (2002) 331–336
The errors of the NSM are studied by Alhama (Alhama, 1999) for both lineal and non-lineal problems including the phase change in pure materials. These errors are acceptable in engineering using a number of volume elements N P 50 for 1-D geometry; average relative deviations between theoretical and simulated temperatures are reduced by about 0.1–1% depending on N.
333
temperature. Other kinds of boundary conditions can be easily assumed. Equations (4) and (5) summarise the discrete values of k and ce into the temperature interval to which the food is subjected in the process. All kinds of real dependencies may be approximated by these equations (interpolation tables) with negligible error by taking high natural values for a and b. For different temperatures from those ones on the table, values for k and ce are found out by simple lineal interpolation.
2. The mathematical model 3. The network model A food product of 1-D, planar, cylindrical or spherical geometry whose initial temperature is T0 is subjected to a thermal phase-change process of freezing or unfreezing, under the application of defined boundary conditions. Natural and forced convection as well as radiation, or isothermal boundary conditions, may be applied to the outer surface of the food. Variations of the thermal conductivity and specific heat or enthalpy in relation to the temperature are specified in a tabulated form by n and m values, respectively, which are associated to n and m discrete temperatures inside the work temperature interval. Continuous arbitrary functions may also be assumed for these dependencies as simpler cases. Making use of Fourier law, the governing partial differential equation for both phases is obtained by performing a heat balance over a small region (volume element) of the food, the form of this region depends on the co-ordinate system. The complete problem statement can be written as qce oT =ot ¼ rðkrT Þ; hðTs Ta Þ ¼ koT =oxjs ; Tðt¼0Þ ¼ T0 ; kðTa Þ ¼ ka ; ce ðTb Þ ¼ ce;b ;
1 < b < m;
½SiþD kiþD ðTi TiþD Þ=ðDxi =2Þ ;
ð6Þ
ð4Þ ð5Þ
ce;i; ¼ ce;i;mean þ Dce;i :
ð2Þ ð3Þ
1 < a < n;
Ci qce;i ðdTi =dtÞ ¼ ½SiD kiD ðTiD Ti Þ=ðDxi =2Þ
where SiD and SiþD are the boundary spherical surfaces of the crown and C the volume. Two values for the conductivity in the volume element are defined; one on the left, kiD , whose value obtained from Eq. (4) is going to be a function of the temperature at that point, TiD ; other one on the right, kiþD (the same but using the temperature TiþD ). This is the conductivity ki D ¼ kðTi D Þ that is evaluated by means of the table given in Eq. (4) by simple lineal interpolation. On the other hand the specific heat for the cell, ce;i , is evaluated according to the temperature in the centre of the cell, ce;i ¼ ce;i ðTi Þ, using Eq. (5). To design the network we use the equivalence ‘‘temperature voltage’’ and ‘‘heat flux electric current’’. Now in order to implement the initial conditions, for which a condenser is necessary, the values of the specific heat are separated into two terms in the form
ð1Þ h ¼ hðTs Þ;
Spherical co-ordinates are used for the design of the network model. For 1-D geometry, Eq. (1) applied to an elemental volume of spherical crown form (which is defined by a thickness Dxi and radii xiD ; xi and xiþD , see nomenclature in Fig. 1) may be written as
where Eq. (1) is the governing equation in vector notation, which is applied to the whole medium; q is the density, k is the thermal conductivity, ce is the specific heat and r is the gradient operator. Eq. (2) is the natural or forced convective boundary condition. As regards the heat transfer coefficient, h, or its dimensionless form Nu ¼ hL=k, with L the characteristic length of the product, standard expressions dependent on temperature for natural convection (such as Yuge, 1960 or Churchill, 1983), or constant values in forced convection may be used; Ta is the temperature of the surrounding fluid. In the end, Eq. (3) is the initial condition. The set of equations (1)–(3) makes up the so-called apparent specific heat formulation. The use of dimensionless form for these equations is not possible because of the variations of the thermal properties with the
ð7Þ
Defining the currents ji;C ¼ Ci qce;i;mean ðdTi =dtÞ ¼ Ci ðdTi =dtÞ;
ð8aÞ
ji;c ¼ Ci qðDce;i ÞðdTi =dtÞ;
ð8bÞ
jiD ¼ SiD kiD ðTiD Ti Þ=ðDxi =2Þ;
ð8cÞ
jiþD ¼ SiþD kiþD ðTi TiþD Þ=ðDxi =2Þ
ð8dÞ
Eq. (6) can be written in the form ji;C þ ji;c ¼ jiD jiþD ; which is a kind of Kirchhoff’s law for these currents. The network model for the volume element or cell is now designed. It is integrated by only four devices and three simple special sources, named voltage-control voltagesource (E), which provide the values of kiD ; kiþD and ce;i . These special sources are Ek;1 ; Ek;2 ; Ece , Fig. 2. In this
334
F. Alhama, C.F. Gonzalez Fernandez / Journal of Food Engineering 54 (2002) 331–336
Fig. 1. Nomenclature and geometry of the volume element.
way, the output voltage of Ek;1 is the value of kiD associated to input temperature TiD , according to Eq. (4). The same is applied to the auxiliary source Ek;2 , which provides kiþD , and to Ece , which provides the value of ce;i associated to the temperature Ti . Eq. (8a) defines a condenser of capacity Ci ¼ Ci qce;i;mean and Eqs. (8b), (8c) and (8d) are implemented by a non-lineal source named control current generator, G. This active device is a current source whose output is evaluated as an arbitrary function of one or more voltages or currents at different points of the network. In particular, Gi;1 is controlled by the voltages Ei;k;1 ; TiD and Ti , according to the expression ‘‘ðSiD 2=Dxi Þ Ek;1 ðTiD Ti Þ’’ and Gi;2 by the voltages Ei;k;2 ; TiþD and Ti , in the form ‘‘ðSiþD 2=Dxi ÞEk;2 ðTi TiþD Þ’’. In the end, Gi;3 is controlled by the voltage Ei;c and the current of Ci according to the function ‘‘ji;C ððEi;ce ce;i;mean Þ= ce;i;mean Þ’’.
Fig. 2 shows all the devices of the network model; variables between brackets and in cursive mean control variables of the adjacent non-lineal sources. Voltages and current source must be connected with the adequate polarity in order to satisfy Eq. (6). Geometrical parameters in each cell (radii xiD ; xi and xiþD , surfaces SiD and SiþD and volume, Ci ) are different for the cylindrical and spherical co-ordinates and constant values for planar co-ordinates. Once the network of the volume element has been designed a number of them (N) are connected in series to conform the whole region. For N P 50 errors are reduced to values less than 0.5% in 1-D geometries (Alhama, 1999). Boundary condition is implemented by a voltagecontrolled current-source in the cases of convection and/or radiation, and by a constant voltage source in the case of isothermal condition. As regards the initial condition, this is applied to the network by fixing the initial voltage of the condenser to the value of T0 . Once the whole network is completed, it must be written in the way of a program that could be read by the computer code Pspice. Only a few rules are necessary for this purpose since only three different devices are needed in the network. It is also possible to change the parameters of the problem by means of simple programming so that this can be adapted to different geometrical or thermal values. The solution is obtained for all the temperatures and fluxes, at the time, by running the program in Pspice. A few seconds of computation time are needed using a PC (Pentium II, 500).
4. Numerical results and discussion
Fig. 2. Network model of the volume element.
An application is made to the process of freezing of white fish by forced convection. For the storage and
F. Alhama, C.F. Gonzalez Fernandez / Journal of Food Engineering 54 (2002) 331–336
preservation of the white fish, temperatures about 20 °C must be reached; the initial temperatures being about 5 °C. With this aim the fish is taken into the freezing unit multi-plate with temperature operation of 40 °C, undergoing forced convection on both sides of the plate. Once the critical point of the fish reaches the required stored value (20 °C) the food is taken to the frozen-food vessel. To simulate this problem a planar 1-D geometry is used, being only necessarily half a plate due to the symmetry of the problem. The parameters are: L ðthickness of the slabÞ ¼ 0:1 m; q ¼ 1054 kg m3 ; T0 ¼ 5 °C; Ta ¼ 40 °C; N ¼ 40. The temperature dependencies of k; ce and enthalpy (H), taken from the literature (FAO, 1985), are shown in Tables 1 and 2. As regards the heat transfer coefficient typical values around h ¼ 600 W m2 K1 , reported by Gruda (Gruda & Potolsky, 1974), have been adopted. Fig. 3 shows the transient temperature at five regular positions, x ¼ 0, L=8; L=4; 3L=8 and L=2, taking x ¼ 0 at the surface of the plate. The time required to reach 20 °C at the centre of the plate is about 2h. At this moment the temperature on the surface of the fish is 38 °C. The effect of the near regions of the surface which is undergoing the phase change over the most internal regions is observed in the figure, as well as the
335
Fig. 3. Transient temperatures at five regular positions.
Table 1 Thermal conductivity of the white fish Temperature (K)
Conductivity (W m1 K1 )
278 273 272 243
0.43031 0.43031 1.30256 1.87243 Fig. 4. (a) Instantaneous heat flux. (b) Integrated heat flux.
Table 2 Specific heat and enthalpy of the white fish Temperature (K)
Enthalpy (J kg1 )
Specific heat (W m1 K1 )
233 253 263 265 267 269 270 271 272 273 275 277 279 283 293 313
0 41985.6 74217.8 83678.1 96319.9 117417 136882 176482 297876 322992 330275 337601 344885 359535 396205 469878
1841.84 2595.32 4227.86 5316.22 7744.1 15111.4 26539.2 65636.5 102724 4144.14 3641.82 3641.82 3641.82 3683.68 3683.68 3683.68
increase of the time required for freezing the same piece of fish as process advances. Instantaneous heat flux crossing the surface is shown in Fig. 4(a). An inflexion point is observed about t ¼ 2h in the curve. The instantaneous heat flux maintains a high value (almost constant) around t ¼ 2h, in which all the fish has changed the phase. A rapid decreasing takes place after this time due to the values of k and ce in the solid state. Integrated heat flux as a function of time, shown in Fig. 4(b), may be compared with the enthalpy step between two temperatures, Table 2. In this way and for the complete process, the enthalpy step between 5 and 40 °C is 341.24 kJ kg1 , nearly the same value as that of the total integrated flux in the simulation 342.43 (error ¼ 0.34%). When the centre of the fish reaches 20 °C the total integrated flux is 316.46; a higher value
336
F. Alhama, C.F. Gonzalez Fernandez / Journal of Food Engineering 54 (2002) 331–336
as we could expect than the enthalpy step between 5 and 20 °C because most of the fish is below 20 °C. Nevertheless the last difference is very small (5.74%) because the heat eliminated from the fish below 20 °C is negligible compared to the latent heat of phase change needed for the complete freezing.
5. Conclusions The network simulation method has been applied for the first time, to the authors’ knowledge, in the way of apparent specific heat formulation to solve the transient problem of freezing foods, whichever be the boundary conditions. The application of the method is simple since very few electrical devices are needed for the network design of the volume element. This, in turn, requires the use of very few programming rules to make the file to be read by the circuit resolution code. The user does not require mathematical manipulations of the finite-difference equations since the computer code specially designed to this end does this work. The only mathematical step is to convert the governing equation into a spatial finite-difference equation. An application is made to the freezing of white fish under forced convection, a severe non-lineal problem due to the strongly temperature dependencies of the thermal properties. Specific heat has a sharply peak near the initial freezing front. Any kind of dependencies for k and ce , and boundary conditions (including natural convection and radiation) may be considered without special requirements. Errors produced by comparing the integrated heat flux with the enthalpy step for two given temperatures are negligible. Computer time for this application is 10–12 s in a PC (Pentium II, 500). References Alhama, F. (1999). Transient thermal responses in non-lineal heat conduction processes. Ph.D. Thesis. University of Murcia, Spain, 234 pp.
Alhama, F., L opez Sanchez, J., & Gonzalez-Fernandez, C. F. (1997). Heat conduction through a multilayered wall with variable boundary conditions. Energy, 22, 797–803. DeBaerdemaeker, J., Singh, R. P., & Segerlind, L. J. (1977). Uses of finite element method in food processing. Journal of Food Engineering, 1, 37. Churchill, S. W. (1983). Free convection around inmerser bodies. In E. U. Schl€ under (Ed.-in-Chief), Heat exchange design handbook, Section 2.5.7. New York: Hemisphere Publishing. Cleland, A. C., & Earle, R. L. (1979). A comparison of methods for predicting the frezzing time of cylindrical and spherical foodstuffs. Journal of Food Science, 44, 964. FAO (1985). Freezing and refrigerates storage in fisheries. Fisheries Technical Paper 340 (ISSN 0429-9345). Gonzalez-Fernandez, C. F., Alhama, F., L opez Sanchez, J. F., & Horno, J. (1998). Application of the network method to heat conduction processes with polynomial and potential-exponentially varying thermal properties. Numerical Heat Transfer A, Applications, 33, 549–559. Gruda, Z., & Potolsky, J. (1974). Zarazanie zywnosci. Wydawnictwa Naukowo-Techniczne Ed. Warszawa. Horno, J., Garcıa Hernandez, M. T., & Gonzalez-Fernandez, C. F. (1993). Digital simulation of electrochemical processes by network approach. Journal of Electroanalytical Chemistry, 352, 83–97. Lin, S. (1972). Similarity solution of freezing or melting process with variable thermal properties. Transactions CSME, 1, 51. Mannapperuna, J. D., & Singh, R. P. (1988). Prediction of freezing and thawing times of foods using a numerical method based on enthalpy formulation. Journal of Food Science, 53(2), 626– 630. Nagel, L. W. (1975). SPICE2: A computer program to simulate semiconductor circuits. Memo. No. UCB/ERL M520. Electronic Research Laboratory, Univ. de California, Berkeley, CA 94720. Polley, S. L., Snyder, O. P., & Kotnour, P. (1980). A compilation of thermal properties of foods. Food Technology, 76–94. PSPICE 6.0 (1994). Microsim Corporation, 20 Fairbanks, Irvine, California 92718. Rao, M. A., & Rizvi, S. S. H. (1986). Engineering properties of foods. Food-Analysis. Series: Food science and technology. New York: Marcel Dekker. Sweat, V. E. (1974). Experimental values of thermal conductivity of selected fruits and vegetables. Journal of Food Science (39), 1080. Voller, V., & Cross, M. (1985). Application of control volume enthalpy methods in the solution of Stefan problems. In Lewis, Morgan, & Johnson (Eds.), Computational techniques in heat transfer. Swansea, UK: Pineeridge Press. Yuge, T. (1960). Experiments on heat transfer from spheres including combined natural and forced convection. Journal of Heat Transfer, 82, 214–220.