Determining thermal parameters in the cooling of a small-scale high-pressure freezing vessel

International Journal of Refrigeration 29 (2006) 1152e1159 www.elsevier.com/locate/ijrefrig Determining thermal parameters in the cooling of a small-scale high-pressure freezing vessel  Bérengère Guignona, Angel M. Ramosb, Juan A. Infanteb, José M. Dı́aza, Pedro D. Sanza,*,1 a Instituto del Frı́o, Centro Superior de Investigaciones Cientı́ficas (CSIC) e Ciudad Universitaria, C/ José Antonio Novais, 10, 28040 Madrid, Spain b Departamento de Matemática Aplicada, Facultad de Matemáticas e Universidad Complutense de Madrid, Pza. de Ciencias, 3, 28040 Madrid, Spain Received 5 January 2005; received in revised form 9 January 2006; accepted 30 January 2006 Available online 5 June 2006 Abstract High-pressure supported freezing processes need a more efficient refrigeration technique to be applied at industrial level. A cooling method consisting in the circulation of a refrigerant in ebullition around the product in the vessel has been tested on a lab-scale prototype built for that purpose. The cooling kinetic of a mixture of ethanol, ethylene glycol and water (a usual pressurizing medium) was followed, recording temperatures in the whole sample. A mathematical model has been developed to describe heat transfer during cooling of the sample in the vessel. The heat transfer coefficient between the refrigerant and the vessel was determined by a fitting procedure between the numerical simulation results and the experimental measurements. This model should be used to predict the cooling kinetics in other conditions (other products, larger vessels) and to optimise the process.
2006 Elsevier Ltd and IIR. All rights reserved. Keywords: Research; Freezing; Liquid; High pressure; Experiment; Evaporator; Modelling; Heat transfer Détermination des paramètres thermiques du refroidissement d’un petit récipient sous haute pression Mots clés : Recherche ; Congélation ; Liquide ; Haute pression ; Expérimentation ; Évaporateur ; Modélisation ; Transfert de chaleur * Corresponding author. Department of Engineering, Instituto del Frı́o, Centro Superior de Investigaciones Cientı́ficas (CSIC) e Ciudad Universitaria, C/ José Antonio Novais, 10, 28040 Madrid, Spain. Tel.: þ34 91 544 56 07; fax: þ34 91 549 36 27.  Ramos), E-mail addresses: angel_ramos@mat.ucm.es (A.M. psanz@if.csic.es (P.D. Sanz). 1 Member of the B1 Commission. 0140-7007/$35.00
2006 Elsevier Ltd and IIR. All rights reserved. doi:10.1016/j.ijrefrig.2006.01.007 1. Introduction Food freezing processes very often require a rapid cooling rate of the product in order to preserve its quality and safety. Nowadays, a suitable and common way for cooling in high-pressure supported freezing processes is to use a secondary refrigerant flowing through a spiral circuit around B. Guignon et al. / International Journal of Refrigeration 29 (2006) 1152e1159 1153 Nomenclature cp h k l l L n ~ n r R R < v< <1 S t T x specific heat. J kg
1  C
1 surface heat transfer coefficient. W m
2  C
1 thermal conductivity. W m
1  C
1 mid-height of the cavity. m coordinate of the entrance of the spiral circuit. m mid-height of the vessel. m number of spiral turns outward normal vector on the boundary of the domain interior radius of the vessel. m exterior radius of the vessel. m radius of the spiral circuit. m cylindrical vessel domain boundary of the domain < cylindrical cavity domain spiral domain time variable. s temperature.  C spatial variable representing the cartesian coordinates x1, x2, x3 of the vessel points. m the product [1]. A high-pressure vessel has thick steel walls to resist pressure in the range of about 700 MPa at low temperatures (down to
40  C). As the consequence of its thermal constrains, the cooling step under pressure may spend long time mainly due to the high thermal mass of this container (around 18 times more than the corresponding only to the product). For industrial applications, this is not acceptable from an energetic and economic point of view. Since food products frozen by a high-pressure process present better quality than the frozen by the conventional ones [2,3], it is a real challenge to overcome this scientific and technological problem which hinders a wider application of this emerging freezing process. Some techniques have been developed to enhance heat exchange rate [4,5]. In the case of high-pressure vessels, it is enhanced mainly by increasing heat exchange surface area: spiral circuits around the product and by using large cooling bath to immerge the vessel on it [6,7]. The first method is limited because when the number of helix turns increases the vessel wall loses its resistance and it appears more pressure losses due to frictions [8]. The second one is only conceivable at laboratory scale. The most effective and recently developed method consists in replacing the traditional cooling fluid by a refrigerant fluid which evaporates inside the heat exchanger of the vessel [9]. In this case the latent heat of vaporization involves higher heat exchange rates than only the sensible heat of usual cooling fluids. The refrigerant fluid comes in the heat exchanger in the liquid state, obtains heat to be evaporated from the sample (which is consequently cooled) and leaves the heat exchanger in gas state. In order to test the efficiency of this Greek letters a thermal diffusivity. m2 s
1

V gradient differential operator V ¼ vxv1 ; vxv2 ; vxv3 boundary of the vessel outside G1 G2 boundary of the symmetry axis G3 boundary of the spiral circuit U half a section of the vessel domain r density. kg m
3 q angle in the spiral equation. rad Superscripts app apparent Subscripts f final m metallic structure ref refrigerant s sample v vessel 0 initial cooling method, a prototype of the high-pressure equipment has been designed and constructed at laboratory scale. Because product quality and treatment uniformity (e.g. micro-organism inactivation) are directly related to temperature distribution [10,11], a numerical model is proposed to determine the evolution of temperature distribution in the sample during cooling in the cylindrical vessel. By assuming that the heat transfer coefficient is not pressure dependent the objective is to evaluate and simulate the cooling behaviour of a model sample to be frozen in this prototype. To do that, it is necessary to identify the heat transfer coefficient between the refrigerant fluid in the cooling circuit and the vessel through a mathematical modelling and experiments performed at atmospheric pressure. 2. Experiments 2.1. Prototype characteristics The prototype vessel has been made of a hollow brass cylinder whose surface was caved mechanically to form a thread. This hollow cylinder was introduced in a larger one whose internal diameter corresponds to the external diameter of the first one as shown in Fig. 1. Therefore, a closed circuit is obtained inside the wall of the vessel. To connect it with the refrigerating installation, two holes were drilled through the wall of the external cylinder: one in the lower part where the refrigerant will enter and the other one in the upper part where it will leave the vessel. The flow-rate of the refrigerant fluid (Freon 404A) was controlled with a manual expansion 1154 B. Guignon et al. / International Journal of Refrigeration 29 (2006) 1152e1159 between the heat exchange surface area and the heat capacity of the vessel was considered as the heat exchange similitude parameter. It was equal to 1.46  10
6 m2  C J
1 for the real vessel and was kept the same for the prototype vessel. The prototype had a heat capacity of 4954 J  C
1. Its surface exchange area was calculated to be about 7.23  10
3 m2 and has been obtained by means of a three turns spiral circuit with a rectangular section (3  5 mm) around the central hole (Fig. 1). 2.2. Experimental procedure Fig. 1. Technical design of the container. valve in the refrigeration installation as shown in Fig. 2. The vessel was plugged at the top and at the bottom to hold the ensemble of thermocouples and the liquid sample. Dimensions of the vessel are 1:3 sized of the ones corresponding to the real high-pressure pilot equipment (ACB GEC Alsthom, Nantes) [1]. So the prototype vessel height is 0.19 m, its interior diameter 0.034 m and its exterior diameter 0.06 m. The ratio Initially the vessel has been filled at ambient temperature with the liquid sample (volume z 90 mL) to be cooled. Three kinds of samples were used: silicon oil, ethanol and a mixture of ethanol, deionised water and ethylene glycol (20:40:40 % v/v). These liquids are chosen because they are used as pressure transmitting medium and they do not freeze in the considered temperature range. This feature is necessary in order to examine the heat transfer behaviour of the vessel at low temperature avoiding freezing. Their thermophysical properties are given in Table 1. Once the refrigerant reaches a low constant temperature (e.g.
21.3  C), the connecting valve of the vessel is opened and the refrigerant is allowed to circulate through the heat exchange spiral. The temperatures of interest is recorded every 3 s with a data acquisition system (Yokogawa Data Collector DC100, Tokyo, Japan) until the temperature at the centre of the sample is well below
17  C. Temperatures are measured by using 25 T-type thermocouples distributed along different sections of the sample. Every 2 cm along the vertical axis, five thermocouples are placed, one at the centre and four at 2 mm from the edge. A metallic structure holds the extremity of each thermocouple in a fixed and precise position Expansion R 404A Evaporator Min. 1 bar Max. 20 bar PA Sample PB Condenser Compressor Oil separator Fig. 2. Experimental refrigeration installation. 1155 B. Guignon et al. / International Journal of Refrigeration 29 (2006) 1152e1159 Table 1 Thermophysical properties of container wall and studied liquids at 0  C r (kg m
3) k (W m
1  C
1) cp (J kg
1  C
1) a (m2 s
1) a b Brass Metallic structure Silicon oila Ethanol Ethylene glycol Water Mixb 8600 110 385 3.32  10
5 8351.1 75.2 418.3 2.15  10
5 960 0.162 1510 1.12  10
7 806.3 0.182 2277 0.99  10
7 1125.6 0.265 2269 1.04  10
7 1000.1 0.594 4208.4 1.14  10
7 1011.7 0.381 3042.3 1.24  10
7 At 20  C. Ethanol: 20%, Ethylene glycol: 40%, Water: 40% (v/v); same rules that described in Section 3.2. (1 mm). Temperatures are also controlled in the mechanical cooling producer and at the entrance and at the exit of the heat exchanger. Accuracy of thermocouples is 0.1  C. Experiments are performed in triplicate. The mean of standard deviations between cooling kinetics is 0.3  C for temperature measurements compared at the same instant. 2.3. Experimental results: cooling kinetics and temperature distribution in the sample Experimental cooling kinetics for a sample of ethanole ethylene glycolewater mix are shown in Fig. 3. Cooling progresses from the surface of the sample to the centre and from the bottom (near refrigerant inlet) to the top as expected. The cooling process was considered to start when the temperature decreases more than 0.1  C from the initial stationary state. The end was considered when the temperature at the centre of the sample was two degrees above the corresponding refrigerant temperature. Three parts may be distinguished in the process: the slow cooling step at
1.2  0.1  C min
1 during its first 2 min, the cooling step at a maximal rate of
4.3  0.1  C min
1 during the following 4 min and, the final cooling step at decreasing rate (near
0.8  0.02  C min
1). In order to characterise the overall cooling process, it is defined the mean cooling rate as the slope of the timeetemperature curve at the half cooling process. It is the temperature change divided by the time change over 10 measurements. Its obtained average value is 1.0  0.1  C min
1 and the average half cooling time obtained is 11 min. When compared to the figures of 0.12  0.04  C min
1 and 90 min, respectively, from the real high-pressure vessel (i.e. with the traditional cooling system), both the mean cooling rate and the half cooling time are considerably reduced (about 8 times lower with the new cooling system). There was no significant differences found between the cooling kinetics of silicon oil, ethanol and ethanole ethylene glycolewater mix samples (data not shown). Although these liquids possess different thermophysical properties, the metallic structure contribution of the experimental set-up masked cooling rate differences due to its higher thermal diffusivity as compared to that of the liquids (Table 1). x3 L l l R Temperature (°C) 20 10 r 0 0 -l -10 -20 R x2 0 10 20 -L 30 Time (min) Fig. 3. Experimental cooling kinetics in the sample. :: mean temperature in the upper part, þ: temperature at the centre of the vesmean peripheral sel, ;: mean temperature in the lower part, temperature in the central part. Fig. 4. Container configuration for modelling. x1 -l 1156 B. Guignon et al. / International Journal of Refrigeration 29 (2006) 1152e1159 Table 2 Input data for the numerical simulation Vessel dimensions (m) R ¼ 0.050 r ¼ 0.017 L ¼ 0.095 l ¼ 0.050 R ¼ 0:030 r ¼ 0:027 Operating conditions T0 ¼ 19  C Tref ¼
21.3  C P ¼ 0.1 MPa Thermophysical properties rs ¼ 1493.3 kg m
3 rv ¼ 8600 kg m
3 cps ¼ 2079 J kg
1  C
1 cpv ¼ 385 J kg
1  C
1 ks ¼ 5.3 W m
1  C
1 kv ¼ 110 W m
1  C
1 Computing parameters tf ¼ 1800 s Dt adaptative href ¼ 2000e5000 W m
2  C
1 Mesh statistics: 58 571 tetrahedra 13 107 nodes ¼ (degrees of freedom) 3. Modelling 3.1. Mathematical description of the cooling container The container is represented by a cylinder with radius R and height 2L. Its inside cavity has a radius r and a height 2l to hold the sample to be cooled. For modelling purposes, the mass centre of the container is considered at the point x ¼ (0,0,0) and its axis along the vertical x3-coordinate axis (Fig. 4). The refrigerant circuit enters transversally through the lower part of the container at the point ðR; 0;
lÞ and exits at the corresponding upper part at the point ðR; 0; lÞ. The trajectory followed by a virtual point of the spiral circuit is then given by the curve:   l ð1Þ R cosq; R sinq;
l þ q ; q ˛½0; 2np np where n is the number of spiral turns (n ¼ 3 in our case). 3.2. Three dimensional model Let be < the domain of the whole cylinder of radius R and height 2L without considering the internal circuit previously described. At the beginning of the cooling process the initial distribution of the temperature in < is given by T0(x). However by one side the liquid sample becomes more and more viscous as it is cooled. By the other side the plates of the metallic structure holding the thermocouples avoid the sample to move freely. For those reasons the convective effect in the liquid is considered to be very small and even to be negligible for calculation purpose. The distribution of temperatures T(x,t) is given by the following heat transfer equations [12]: l ¼ 0:063 where v< is the boundary of the domain; S is the surface area of < at the boundary of the internal circuit and ~ n is the outward normal vector at the boundary of the domain. T0 is the initial temperature distribution (  C). During the cooling process, the refrigerant evaporates as soon as it enters the surface heat exchanger but at the end, the refrigerant may be considered as mainly be in the liquid state within it. Therefore, for calculation purposes the refrigerant temperature Tref (  C) is considered to be constant and equal to its boiling temperature at the experimental conditions. href is the surface heat transfer coefficient between the boundary of the container and the cooling fluid (W m
2  C
1). Since the container has been fully covered with foam, thermal insulation is supposed on the external surface of the brass vessel. Parameters r, cp and k are, respectively, the density (kg m
3), the specific heat (J kg
1  C
1) and the thermal conductivity (W m
1  C
1) of the material:  r ðTÞ if x ˛ < ðsample to be cooledÞ 1 s rðT; xÞ ¼ rv ðTÞ if x ˛ <n<1 ðvesselÞ  c ðTÞ if x ˛ < ðsample to be cooledÞ ps 1 cp ðT; xÞ ¼ cpv ðTÞ if x ˛ <n<1 ðvesselÞ  ks ðTÞ if x ˛ <1 ðsample to be cooledÞ kðT; xÞ ¼ kv ðTÞ if x ˛ <n<1 ðvesselÞ ð3Þ The thermophysical properties of the vessel (domain <n<1 ) have been considered to be independent of the temperature. The effect of the referred metallic structure holding thermocouples on the heat transfer process has been considered to be included as an additional contribution to 8
 vT > > rðT; xÞcp ðT; xÞ ðx; tÞ
V$kðT; xÞVTðx; tÞ ¼ 0 in <  0; tf > > > vt > >
 < vT on ðv<nSÞ  0; tf kðT; xÞ ðx; tÞ ¼ 0 v~ n >

 > vT > > on S  0; tf kðT; xÞ ðx; tÞ ¼ href Tref
Tðx; tÞ > > v~ n > : in < Tðx; 0Þ ¼ T0 ðxÞ ð2Þ 1157 B. Guignon et al. / International Journal of Refrigeration 29 (2006) 1152e1159  density: additive model: rapp ¼ ððxs =rs Þ þ ðxm =rm ÞÞ
1 with xi mass fraction of component i,  specific heat: additive model: capp p ¼ xs Cs þ xm Cm ,  thermal conductivity: in this case a parallel model has been used: k app ¼ 3s ks þ 3m km where 3i ¼ xi ðrapp =ri Þ is the volume fraction. As an example, in the case of the ethanoledeionised watereethylene glycol mix (20:40:40 % v/v) and a weight contribution of sample of 63.3% and of metallic part of 36.7%, calculations gave: rapp ¼ 1493.3 kg m
3, capp p ¼ 2079.3 J kg
1  C
1, kapp ¼ 5.3 W m
1  C
1. All the input data for the numerical simulation are given in Table 2. The numerical approximation of the solution of the set of Eq. (2) has been computed by Finite Element Methods by using the commercial package FEMLAB (nowadays COMSOL Multiphysics, http://www.comsol.com/ products/multiphysics), solving with Lagrange-Quadratic tetrahedral Finite Elements. where Texp(t), T(t) are the temperature of the experiment and the numerical approximation, respectively, for times t from 0 to 1800 s (see [14]). Since the 3D mathematical model is very time consuming a faster (but less reliable) 2D model has been developed 20 10 Temperature (ºC) the thermophysical properties of the sample. Apparent thermophysical properties of the sample are then calculated giving a contribution weight to the own sample and a contribution weight to the corresponding metallic part in the following way [13]: href=2000 0 -10 -20 0 5 10 ð 1800 0 Temperature (ºC) In order to choose the suitable coefficient a least square method has been used. Different surface heat transfer coefficients (href) corresponding to the boiling fluid on the spiral heat exchange circuit (Fig. 4) have been used. The final value has been chosen by minimizing, 2 Texp ðtÞ
TðtÞ dt 25 30 25 30 25 30 href=4750 0 -10 -20
20 20 10 3.3. Results: numerical simulation and optimal heat transfer coefficient 15 Time (min) ð4Þ 0 5 10 15 20 Time (min) 20 Temperature (ºC) 10 href=2750 0 -10 -20 0 5 10 15 20 Time (min) Fig. 5. The domain U and the temperature distribution for the 2D-model at time t ¼ 1800 s. Fig. 6. Cooling kinetics in the central part of the sample. e e: experimental, d: numerical simulation. 1158 B. Guignon et al. / International Journal of Refrigeration 29 (2006) 1152e1159 for getting a suitable initial guess for the value of href . By assuming that there is not vertical heat conduction, the 3D model can be simplified as: 8 vT > > > rðT; xÞcp ðT; xÞ vt ðx; tÞ
V$kðT; xÞVTðx; tÞ ¼ 0 > > > > > > kðT; xÞvT ðx; tÞ ¼ 0 > < v~ n vT kðT; xÞ ðx; tÞ ¼ 0 > > > v~ n > > 
vT > > > ðx; tÞ ¼ href Tref
Tðx; tÞ > > kðT; xÞ v~ n : Tðx; 0Þ ¼ T0 ðxÞ prototype for the href optimal value at time 1800 s. Besides, the experimental temperature values corresponding both to the centre and to the mean from four temperature edge points
 in U 0; tf
 on G1 0; tf
 on G2 0; tf
 on G3 0; tf ð5Þ in U where U is the 2D domain representing half a section of the vessel at an arbitrary height and G1 , G2 and G3 are the parts of its boundary corresponding to the outside of the vessel, the symmetry axis providing the half of the section and a part of the spiral circuit, respectively, (Fig. 5). The optimal value for href when using the 2D model was 4750 W m
2  C
1. Fig. 5 shows the 2D domain U and the temperature distribution at time t ¼ 1800 s. This value gave an idea of the range of suitable real values and was used in the 3D model as the initial guess of href when minimizing the value of expression (4). Finally, 2750 W m
2  C
1 has been obtained as the optimal value for href. Results for different coefficients href in the range 2000e5000 W m
2  C
1 are given in Fig. 6. Fig. 7 shows the 3D calculated temperature distribution of the whole for five different levels are also indicated. The agreement between calculated and experimental results is close to the experimental error. 4. Conclusions The described mathematical model appears as an interesting tool to determine the heat transfer coefficient for the cooling down of high-pressure vessels. This seems also a powerful way to optimise the different components of the vessel accounting for the heat exchange characteristics. The cooling behaviour of samples with different thermophysical properties and in larger vessels could be predicted using this model. Acknowledgements t = 1800 s Slice:Temperature 80 Experimental Max -19.38 Tc Tp 60 -18.8 -19.3 -18.9 -18.9 -19.42 This work has being carried out with the support of the Spanish ‘‘Plan Nacional de IþDþI (2004e2006). MCYT’’ through the ‘‘AGL2003-06862-C02 Project’’. References 40 20 -19.46 0 -20 -18.8 -18.8 -19.0 -18.9 -19.1 -18.2 -19.50 -40 -60 -19.54 -80 -100 50 -19.58 50 0 Min 0 Fig. 7. 3D-temperature map of the vessel at time 1800 s for href ¼ 2750 W m
2  C
1 and comparative experimental values. Tc is the temperature at the centre and Tp is the mean from four temperatures at the edge points. [1] L. Otero, A.D. Molina-Garcı́a, A.M. Ramos, P.D. Sanz, A model for a real thermal control in high-pressure treatment of foods, Biotechnol. Prog. 18 (2002) 904e908. [2] A. Le Bail, D. Chevalier, D.M. Mussa, M. Ghoul, High pressure freezing and thawing of foods: a review, Int. J. Refrigeration. 25 (2002) 504e513. [3] M.N. Martino, L. Otero, P.D. Sanz, N.E. Zaritzky, Size and location of ice crystals in pork frozen by high-pressure assisted freezing as compared to classical methods, Meat Sci. 50 (3) (1998) 303e313. [4] V. Zimparov, Energy conservation through heat transfer enhancement techniques, Int. J. Energy Res. 26 (2001) 675e696. [5] D.G. Prabhanjan, G.S.V. Raghavan, T.J. Rennie, Comparison of heat transfer rates between a straight tube heat exchanger and a helically coiled heat exchanger, Int. Commun. Heat Mass Transfer 29 (2) (2002) 185e191. B. Guignon et al. / International Journal of Refrigeration 29 (2006) 1152e1159 [6] C. Luscher, O. Schlüter, D. Knorr, High pressure-low temperature processing of food: impact on cell membranes, texture, color and visual appearance of potato tissue, Innov. Food Sci. Emerg. Technol. 6 (1) (2005) 59e71. [7] M. Thiebaud, E.M. Dumay, J.-C. Cheftel, Pressure-shift freezing of o/w emulsions: influence of fructose and sodium alginate on undercooling, nucleation, freezing kinetics and ice crystal size distribution, Food Hydrocolloids 16 (6) (2002) 527e545. [8] D.S. Austen, H.M. Solimann, Laminar flow and heat transfer in helically coiled tubes with substantial pitch, Exp. Therm. Fluid Sci. 1 (1988) 183e194. [9] J.M. Dı́az Serrano, R.M. Guerra Garcı́a, B. Guignon, P.D. Sanz Martı́nez, Procedimiento para acortar el tiempo de enfriamiento de recintos destinados al tratamiento de los alimentos por alta presión. Decreasing the cooling time in the vessel used for the high-pressure food processing. Patente P200400547; 05/03/04 Oficina Española de Patentes y Marcas. [10] A. Ardia, D. Knorr, V. Heinz, Adiabatic heat modeling for pressure build-up during high-pressure treatment in liquidfood processing. Trans IchemE, Part C, Food Bioprod. Process. 82 (C1) (2004) 89e95. 1159 [11] C. Hartmann, J.P. Schuhholz, P. Kitsubun, N. Chapleau, A. Le Bail, A. Delgado, Experimental and numerical analysis of the thermofluiddynamics in a high-pressure autoclave, Innov. Food Sci. Emerg. Technol. 5 (4) (2004) 399e411. [12] J.I. Dı́az, A.M. Ramos, L. Otero, A. Molina, P.D. Sanz, On the mathematical modelling and control of high hydrostatic pressure food processing, in: E. Balsa-Canto, J. Mora, J.R. Banga, E. Oñate (Eds.), Mathematical and Computer Techniques for Agro-Food Technologies, International Center for Numerical Methods in Engineering (CIMNE), Barcelona, 2002, pp. 170e175. [13] D.R. Heldman, Food freezing, in: D.R. Heldman, D.B. Lund (Eds.), Handbook of Food Engineering, Marcel Dekker, Inc., New York, USA, 1992, pp. 277e315. [14] R. Glowinski, A.M. Ramos, A numerical approach to the Neumann control of the CahneHilliard equation, in: R. Glowinski, H. Kawarada, J. Periaux (Eds.), Computational Methods for Control Applications, Gakuto, International Series: Mathematical Sciences and Applications, vol. 16, Gakko Tosho, Tokyo, Japan, 2002, pp. 111e155.