## Abstract

High-pressure freezing processes are a novel emerging technology in food processing, offering significant improvements to the quality of frozen foods. To be able to simulate plateau times and thermal history under different conditions, in this work, we present a generalized enthalpy model of the high-pressure shift freezing process. The model includes the effects of pressure on conservation of enthalpy and incorporates the freezing point depression of non-dilute food samples. In addition, the significant heat-transfer effects of convection in the pressurizing medium are accounted for by solving the two-dimensional Navier–Stokes equations. We run the model for several numerical tests where the food sample is agar gel, and find good agreement with experimental data from the literature.

## 1. Introduction

### (a) Food freezing

Freezing is a widespread food preservation technology, as it ensures high food quality with long storage duration, and also because it has an extended implementation area (meat, fish, fruit, vegetables, dairy and egg products, etc.). Despite the benefits, freezing of foods can also cause undesirable changes in their texture and organoleptic properties, and its main drawback is the risk of food damage due to the formation of big ice crystals. Freezing can be defined as the crystallization of liquid water into its solid form (ice). Nucleation of ice is a stochastic phenomenon. However, after nucleation, crystallization always initially occurs under a known temperature defined as the initial freezing temperature. Crystal characteristics (size, location, etc.) are dependent on the freezing rate (slow freezing produces large crystals, while rapid freezing produces small ice crystals). The general purpose of food technologists working on this area has been to develop a freezing process that creates a homogeneous matrix of small ice crystals (Denys *et al.* 1997). Improvement of known freezing methods and development of new techniques are important research goals for the food industry at present. With the recent increasing impact of high-pressure technology on food processing, there has been a lot of research dealing with the potential applications of high-pressure effects on ice–water transitions, given that pressure decreases the freezing and melting point of water to a minimum of −22^{°}C at 207.5 MPa (Bridgman 1912). Therefore, several high-pressure freezing and thawing processes have been developed (Denys *et al*. 1997; Otero & Sanz 2000).

According to the path followed by the process in the phase diagram of water (Bridgman 1912), three different types of high-pressure freezing processes can be distinguished in terms of the way in which the phase transition occurs (Fernández *et al.* 2006): high-pressure assisted freezing (HPAF), high-pressure shift freezing (HPSF) and high-pressure induced freezing (HPIF). A few studies on HPIF have been carried out (Kowalczyk *et al.* 2004; Urrutia Benet *et al.* 2004), and there are numerous papers and reviews dealing with HPAF (Otero & Sanz 2003; Fernández *et al.* 2006; Norton *et al.* 2009) and HPSF (Otero & Sanz 2000, 2006; Sanz & Otero 2000). It is generally accepted that for the purposes described above, HPSF processes are the most advantageous (Fernández *et al.* 2006). Therefore, in this paper, we concentrate only on HPSF processes. In §1*b*, we describe the main features of HPSF processes and in §1*c* the state of the art is described. In §2, we present a mathematical model to describe a generic freezing process, and then couple it with a heat-transfer model, to finally give a new model describing a HPSF process. In §3, we present numerical results using known data from real experiments. In §4, we outline the final remarks.

### (b) Main features of a high-pressure shift freezing process

A sample subjected to a HPSF process (figure 1) is cooled under pressure at less than 0^{°}C and kept in the non–frozen state according to the corresponding phase diagram (points A–C in figure 1). When the desired temperature is reached in the product (point C), pressure is released. Two different processes may be distinguished depending on the pressure release rate: *rapid expansion* (in seconds) or *slow expansion* (in minutes). Phase transition occurs due to the pressure release, which induces uniform supercooling throughout the whole sample (point D for rapid expansions, point 1 for slow expansions) owing to the isostatic nature of pressure. This supercooling induces uniform formation of nuclei throughout the sample (regardless of its shape or size), then latent heat is released, raising the sample temperature to the corresponding freezing point (point E for rapid expansion, 2 for slow expansion). Freezing is then completed at constant pressure, usually at atmospheric conditions; for practical reasons of applicability, expansions are usually made from high to atmospheric pressure (i.e. rapid expansion, A–B–C–D–E in figure 1) and therefore the phase transition mainly occurs at atmospheric pressure in this kind of experiment. As a general rule (but taking into account the phase diagram of water), the higher the pressure and the lower the temperature before expansion, the more ice nuclei are formed, and hence the shorter the phase-transition time.

Various authors (Otero *et al.* 2000; Schlüter *et al.* 2004; Fernández *et al.* 2006, 2008; Alizadeh *et al.* 2009) have proved experimentally that in a HPSF process, ice nucleation occurs homogeneously throughout the whole volume of the product and not only on the surface, as they have found small granular shaped ice crystals dispersed throughout the resulting sample for several products. When comparing HPSF to classical freezing processes, important reductions of freezing times have been reported (Schlüter *et al.* 2004; Fernández *et al.* 2006; Otero & Sanz 2006). Applications of HPSF for food are still under development and the amount of available data are increasing accordingly. Most of the studies that are being carried out focus on the advantageous effects that HPSF has on the texture and structure of various products. Denys *et al.* (1997) compared the HPSF process with a conventional (at atmospheric pressure) freezing process, and concluded that the nucleation is faster and more uniform in the HPSF processes, producing a homogeneous crystallization and therefore a better product in terms of texture. Fernández *et al.* (2006) compared the HPAF and HPSF processes on gelatin gel samples and concluded that HPSF was clearly more advantageous: shorter phase-transition times and homogeneous distribution of small ice crystals throughout the sample. HPSF has been successfully applied in the processing of fruits, pork, lobster and tofu, among other products.

### (c) Modelling of high-pressure shift freezing processes: state of the art and needs

Owing to the advantages of this method (mainly, its potential to improve the kinetics of the process and the characteristics of the ice crystals thus formed) compared with other freezing methods, apart from the published research focused on the impact of high pressure on quality aspects of the particular processed food, some modelling studies considering the temperature evolution during treatment of HPSF processes have been published to date. This is important in all high-pressure processes, as pointed out by Otero & Sanz (2003), because some industrial processes have to occur at a constant temperature or between a minimum–maximum threshold to avoid altering some properties of the food (gelification or crystalline state, protein stability, fat migration, freezing, etc.) or to ensure a uniform distribution of the pursued pressure effects (microbial or enzymatic inactivation, uniform nucleation, etc.). In particular, when studying HPSF processes, since phase changes are driven by both pressure and temperature, heat-transfer effects play a major role in this field. Therefore, the research in this area has been focused on theoretically based heat-transfer models that allow us to predict the temperature history within a product freezing under pressure.

Some notable modelling studies are those of Denys *et al.* (1997), where the authors used an explicit two-dimensional finite-difference scheme to simulate temperature profiles during high-pressure freezing of a food simile (tylose) and obtained good agreement with experimental data, but remarked that an improvement needed is the contribution of convection heat transfer by the high-pressure fluid. Sanz & Otero (2000) modelled a HPSF process of a finite cylinder using equations based on the crossed product of infinite slabs and infinite cylinders. The time needed to complete the precooling and change of phase stages is reproduced satisfactorily by the model; however, there are misalignments in the tempering phase because experimental and theoretical conditions were not the same, owing to the thermal gradients. The authors remark that for this stage, a more accurate model (based on finite elements/differences) would be needed. Kowalczyk & Delgado (2007) performed dimensional analysis of the governing equations describing high-pressure processes with forced and free convection, giving a general overview of thermo-fluid-dynamical mechanisms of these kinds of processes. Norton *et al.* (2009) developed a one-dimensional finite-difference numerical model based on the enthalpy formulation to simulate high-pressure freezing of tylose, agar gel and potatoes. This model by Norton *et al.* uses a version of the enthalpy equation derived for constant pressure systems, adapting it to their system by accounting for pressure effects on the latent heat. In addition, their model of freezing point depression of the food samples relies on the Schwartzberg equation (Schwartzberg 1976), which was derived for dilute solutions. At the high ice fractions realized in HPSF systems, the unfrozen gel between ice crystals is non-dilute, and a more general model may be required. In this work, we introduce a new enthalpy model by starting from the pressure-dependent conservation of the enthalpy equation and determining the non-dilute freezing point depression directly from experimental measurements.

## 2. Building a new model for a high-pressure shift freezing process

We propose a two-dimensional axially symmetric mathematical model for a HPSF process, derived from an enthalpy formulation based on volume fractions dependent on temperature and pressure. This model is valid for solid-type foods with a big and small filling sample versus pressurizing media ratio. Convection effects in the pressurizing fluid are taken into account when necessary.

In order to simulate a HPSF process, we first present a general heat-transfer model in a high-pressure process, and then we modify it to take into account the solidification process, by deriving a model based on the enthalpy formulation at non-constant pressure. Also, we present the equations that we use to calculate the amount of ice instantaneously produced just after expansion in a HPSF process, the amount of ice formed during the rest of the process, and the liquid volume fraction. Finally, we give the complete models used to simulate a HPSF process.

### (a) Modelling heat transfer in a general high-pressure process

When high pressure is applied in food technology, it is necessary to take into account the thermal effects that are produced by variations of temperature due to the compression/expansion that takes place in the food sample and in the pressurizing medium. The pressure evolution, *P*(*t*), is known as it is imposed by the user and the limits of the equipment. However, for the temperature evolution, *T*, it is necessary to take into account the adiabatic heating effects owing to the work of compression/expansion in the considered high-pressure device. The temperature of the processed food may change with time and with space, therefore, we need a heat-transfer model capable of predicting the temperature for the processed food. Following Infante *et al.* (2009), a heat-transfer model taking into account only conduction effects is presented when we have a solid-type food with a large filling ratio. A model also including convection effects in the pressurizing fluid, for a solid-type food with a small filling ratio is considered too (for a complete model of a liquid-type food, where convection effects have to be included also in the food, see Infante *et al.* (2009)). As the model is both time and spatially dependent, we introduce a brief description of the domain describing the high-pressure device considered in our simulations.

Usually, high-pressure experiments on food are carried out in a cylindrical pressure vessel (typically a hollow steel cylinder) that is filled with the food and the pressurizing fluid. We assume, owing to the characteristics of this kind of processes, axial symmetry, which allows the use of cylindrical coordinates, and consider a two-dimensional domain given by half a cross section (figure 2). The following four sub-domains are specified:

—

*Ω*_{F}: domain that contains the food sample;—

*Ω*_{C}: cap of the sample holder (typically rubber);—

*Ω*_{P}: domain occupied by the pressurizing medium; and—

*Ω*_{S}: domain of the steel that surrounds the rest of the domains.

The domain in the (*r*,*z*)-coordinates is the rectangle *Ω*=[0,*L*]×[0,*H*] defined by . The boundary of *Ω* is denoted by *Γ*, where we can distinguish:

—

*Γ*_{r}⊂{*L*}×[0,*H*], where the temperature is known;—

*Γ*_{up}=[0,*L*]×{*H*}, where heat transfer with the room in which equipment is located may take place; and—

*Γ*_{0}=*Γ*\{*Γ*_{r}∪*Γ*_{up}}, that has zero heat flux, either by axial symmetry or by isolation of the equipment.

We use star notation ([ ]*) to denote the three-dimensional domains generated by rotating all the domains explained above along the axis of symmetry ({0}×(0,*H*)). These spatial considerations are typical for experimental high-pressure equipments, like the pilot unit used for our numerical experiments presented in §3. For the mathematical modelling, two significantly different cases can be studied (see Infante *et al.* 2009): solid- and liquid-type foods. In this paper, we deal only with the case of solid-type foods. First, we study solid-type foods with a large filling ratio, and therefore a model that takes into account only conduction effects (presented in §2*a*(i)) can give quite precise results (Otero *et al.* 2007; Infante *et al.* 2009). We will also study solid-type foods with a small filling ratio, in which the model includes the convection effects (presented in §2*a*(ii)) in the pressurizing medium.

#### (i) Heat transfer by conduction

When solid-type foods are considered, the starting point is the heat-conduction equation for temperature *T* (K),
2.1where *ρ* is the density (kg m^{−3}), *C*_{p} the specific heat (J kg^{−1} K^{−1}), *k* the thermal conductivity (W m^{−1} K^{−1}) and *t*_{f} the final time (s). The right-hand side of equation (2.1) is the heat production due to the change of pressure *P*=*P*(*t*) (Pa) applied by the equipment (chosen by the user within the machine limitations) and *β* is the thermal expansion coefficient, which is given by
This term results from the following relation (for derivation, see Knoerzer *et al.* 2007), which is valid for isentropic processes:
2.2where Δ*T* denotes the temperature change due to the pressure change denotes the volume and *M* denotes the mass.

Equation (2.1) has to be completed with appropriate boundary and initial conditions depending on the high-pressure machine and the problem we want to solve. We use the same conditions as in Otero *et al.* (2007) for a pilot unit (ACB GEC Alsthom, Nantes, France) located at the Instituto del Frío, Consejo Superior de Investigaciones Cientificas, Spain, which we use in §3 for numerical experiments,
2.3where **n** is the outward unit normal vector on the boundary of the domain, *T*_{0} the initial temperature, *T*_{cool} the cooling temperature that is constant on *Γ**_{r} (cooling the food sample), *T*_{env} the environment temperature (constant) and *h* (W m^{−2} K^{−1}) the heat-transfer coefficient.

By using cylindrical coordinates and taking into account axial symmetry, systems (2.1) and (2.3) may be rewritten as the following two-dimensional problem: 2.4

This model is suitable when the filling ratio of the food sample inside the vessel is much higher than that in the pressurizing medium, since convection effects due to the pressurizing fluid can be neglected. This has been confirmed to be true in Otero *et al.* (2007), by validation with several comparisons between numerical and experimental results. Otero *et al.* (2007) showed that when the filling ratio of the food inside the vessel is not much higher than in the pressurizing medium, the solution of this model differs a lot from the experimental results. Therefore, they improve the model by including convection effects in the pressurizing medium. We present this model in §2*a*(ii).

Another case when (2.4) may be used, even if the filling ratio of the food sample is not much higher than in the pressurizing medium, is for highly viscous fluids—which could be the case at low temperatures—where convection effects may be negligible. This would reduce the computing time without any significant loss in prediction.

#### (ii) Heat transfer by conduction and convection

The non-homogeneous temperature distribution induces a non-homogeneous density distribution in the pressurizing medium and consequently a buoyancy fluid motion, i.e. free convection. In order to take into account the fact that the fluid motion influences the temperature distribution, a non-isothermal flow model is considered. We assume that the fluid velocity, **u** (m s^{−1}), satifies the Navier–Stokes equations for compressible Newtonian fluid under Stokes' assumption (Aris 1989). The resulting system, with appropriate point, boundary and initial conditions, is (Infante *et al.* 2009)
2.5where **g** is the gravity vector (m s^{−2}), *η* the dynamic viscosity (Pa s), *p*=*p*(*x*,*t*) the pressure generated by the mass transfer inside the fluid, and *P*+*p* the total pressure (Pa) in the pressurizing medium *Ω**_{P}. **A**_{1} is a corner point of , which is the boundary of *Ω**_{P} (figure 2).

As in §2*a*(i) for the conductive heat-transfer model, system (2.5) can also be rewritten as an equivalent two-dimensional problem by using cylindrical coordinates (system not shown in this paper). The numerical experiments considered in this paper were carried out using the two-dimensional version of the corresponding equations.

### (b) Modelling a solidification process using the enthalpy formulation

In general, the difficulties of any solidification process are to control the position of the solid/liquid interface and to deal with the release of latent heat, which evolves over a very small temperature range. Crank (1987) classifies the numerical methods for solving the ‘moving-boundary’ problem into *front-tracking* methods and *fixed-grid* methods. The former methods may be used when there is a distinct phase change and a smooth continuous front (Voller *et al.* 1990). However, as the solid/liquid interface becomes less distinct (e.g. if it does not move smoothly or monotonically with time), it may sometimes be difficult or even impossible to track the moving boundary directly. It may have sharp peaks, or it may even disappear. Therefore, the possibility of reformulating the problem in such a way that there is no need to track the position of the solid/liquid interface, but instead it is bound up in a new form of the equations, which applies over the whole of a fixed domain, is an attractive one. These are the so-called *fixed-grid* methods in which the position of the moving boundary appears, *a posteriori*, as a part of the solution.

The essential feature of the fixed-grid methods is that the latent heat evolution is accounted for in the governing energy equation by defining either a total enthalpy, an apparent specific heat or a heat source term (Voller *et al.* 1990). Consequently, the numerical solution can be carried out on a space grid that remains fixed throughout the calculations. Another advantage of these methods is that the numerical treatment of the phase change can be achieved through simple modifications of existing heat-transfer numerical methods and/or software. To model HPSF processes, we use a combination of such methods, adapted to the case of non-constant pressure.

For a pure non-convecting material obeying Fourier's law of heat conduction, conservation of enthalpy can be written as (Bird *et al.* 1960)
2.6where *e* is the enthalpy per unit volume (J m^{−3}). If the pressure is constant, then equation (2.6) reduces to the enthalpy equation used in Voller *et al.* (1990). In the enthalpy model of a HPSF process developed by Norton *et al.* (2009), the pressure term on the right-hand side of (2.6) was neglected.

In the general case of non-constant pressure, the enthalpy is a function of both temperature and pressure, and from thermodynamics, we may derive the relation (Elliott & Lira 2005)
2.7Integrating (2.7) from some reference temperature *T*_{ref} and reference pressure *P*_{ref} to the current temperature and pressure yields an expression for the enthalpy,
2.8

In the case of a mixed-phase region composed of ice crystals (which we will refer to here as the solid part) and the solid and liquid parts of the food that has not frozen yet (which we will refer to here as the liquid part), we may extend previous enthalpy models based on phase fractions (Voller *et al.* 1990) by defining a ‘mixture enthalpy’ as
2.9where *g*_{l} is the volume fraction of the liquid part, *g*_{b} the volume fraction of solids in the food (therefore, *g*_{l}−*g*_{b} is the volume fraction of water in the food) and *λ* the latent heat of freezing of water (J kg^{−1}). We consider that *g*_{l,ref}=*g*_{b}, that *T*_{ref} is the temperature (K) at which all the latent heat has been released (typically *T*_{ref}=−40^{°}C=233.15 K, see Schwartzberg 1976), and that *P*_{ref}=*P*_{atm} (Pa), the atmospheric pressure. Similarly, a ‘mixture conductivity’ is defined as
2.10In the previous equations, subscripts [ ]_{s} and [ ]_{l} refer to the solid and liquid phases, respectively. Taking the total derivative of (2.9) gives
2.11where *C*_{vol}=(1−*g*_{l})*ρ*_{s}*C*_{ps}+*g*_{l}*ρ*_{l}*C*_{pl}, and *β*_{vol}=(1−*g*_{l})*β*_{s}+*g*_{l}*β*_{l}. Inserting (2.11) into the enthalpy equation (2.6) leads to
2.12

For a HPSF process with rapid expansion (see ABCDE in figure 1), typically the pressure profile is as follows:
2.13where *t*_{p1} is the time at which the maximum pressure, *P*_{max}, has been reached; *t*_{p2} is the time at which we start releasing pressure; and *t*_{p3} is the time at which all the pressure has been released down to atmospheric pressure, *P*_{atm}≈0.1 MPa.

Experiments suggest (Otero & Sanz 2006) that as the pressure is rapidly released (*t*∈(*t*_{p2},*t*_{p3})), nucleation is at first delayed (leaving the system in a metastable state with respect to the formation of ice) until the pressure nears atmospheric when ice suddenly nucleates and forms in some small time interval. Therefore, we assume that there is no ice while the pressure is changing (until *t*=*t*_{p3}), so during this stage, equation (2.12) becomes
2.14Once the pressure has been released, ice suddenly nucleates and forms in some small time interval (*ε*), where we assume that the change in ice fraction owing to the drop in pressure is some known function of time, hence (2.12) now becomes
2.15
where now , since *P*=*P*_{atm}. After this point and until the end of the process, we assume that the rest of the ice is computed as a function of temperature (an explicit formula for *g*_{l} is given in §2*c*(iv)), which leads to
2.16

As will be shown in §3, the model given by (2.12) (equivalent, under the previous assumptions, to (2.14)–(2.16)) yields satisfactory agreement between theory and experiment. All the parameters involved in the model can be easily found in the literature for many food similes, except for the liquid fraction, *g*_{l}. In §2*c*, we explain different characteristics of a HPSF that will enable us to derive a formula for *g*_{l} for a gel food simile as a function of time and temperature. We note, though, that a more rigorous approach would be to combine with equation (2.12) a stochastic nucleation law for *g*_{l}(*T*,*P*) quantifying the ice formed at each instance as the pressure is released. As no such expression for *g*_{l}(*T*,*P*) is currently available, we leave this to future work.

### (c) Deriving an expression for the volume fraction *g*_{l} for a gel food simile

#### (i) Supercooling reached after expansion

Otero & Sanz (2006) define the extent of supercooling, Δ*T*_{sc} (^{°}C), as the difference between the lowest temperature at the sample centre just before nucleation, *T*_{min}, and the sample freezing temperature, *T*_{F}, at pressure *P*_{N}, where the nucleation takes place. The extent of supercooling is a crucial factor in the dynamics of a freezing process. In conventional freezing at atmospheric pressure, it is generally admitted that supercooling and nucleation only occurs at the surface of the sample. Otero & Sanz (2006) show that in HPSF experiments, upon pressure release, a metastable state was reached throughout the sample before nucleation. As can be seen in figure 1, after expansion from (*P*_{max},*T*_{prev}), the pressure/temperature coordinates of the sample move to (*P*_{N},*T*_{min}), indicating extensive supercooling. Therefore, the extent of supercooling only depends on the minimum temperature reached after expansion and on the pressure at which nucleation occurs. The minimum temperature, *T*_{min}, reached at nucleation pressure, *P*_{N}, after a rapid pressure release can be estimated according to (2.2), and is accounted for in our heat-transfer model. However, the prediction of *T*_{min} using (2.2) implies foreknowledge of the nucleation pressure, *P*_{N}, and this is, in general, very difficult owing to the stochastic nature of the nucleation phenomenon. In §2*c*(ii), we explain how to avoid this problem.

#### (ii) Modelling the amount of ice formed instantaneously after expansion

Otero & Sanz (2006) explain that the pressure release in HPSF processes can be divided, ideally, into two different phases (figure 1). In the first phase, expansion takes place under metastable conditions. A percentage of water is instantaneously frozen when *P*_{N} is reached. The latent heat raises the sample temperature to the corresponding freezing point. The amount of ice formed at *P*_{N} can be calculated from the following equation (other equations are given in Otero & Sanz 2006):
2.17where *m*_{i} is the mass fraction of ice (defined as mass of ice divided by initial mass of water) formed after expansion, (J kg^{−1} ^{°}C^{−1}) is the specific heat capacity of ice at the nucleation pressure (taken as the mean value of the specific heat capacity at the minimum temperature reached after expansion and the specific heat capacity at the corresponding melting point), (J kg^{−1} ^{°}C^{−1}) is the specific heat capacity of water at the nucleation pressure (taken as the mean value of the specific heat capacity at the minimum temperature reached after expansion and the specific heat capacity at the corresponding melting point), Δ*T*_{sc} (^{°}C) is the extent of supercooling, and *λ*(*P*_{N}) (J kg^{−1}) is the latent heat at nucleation pressure *P*_{N}.

Latent heat of freezing of water (J kg^{−1}) as a function of pressure *P* (MPa) can be estimated by the following equation (for more details, see Otero & Sanz 2006):
2.18Therefore, we can calculate the percentage of instantaneously frozen water, *m*_{i}, by using equations (2.17) and (2.18), taking into account the extent of supercooling attained and the latent heat released at nucleation pressure. However, this is just a theoretical evaluation and is not useful for modelling purposes because experimental data of nucleation pressure and temperature during expansion (*P*_{N},*T*_{min}) are needed and, as was pointed out in §2*c*(i), it is in general very difficult to predict the nucleation pressure. To overcome this problem, Otero & Sanz (2006) proposed a simplified method for HPSF processes with rapid expansions: to assume that nucleation occurs at atmospheric pressure. On this basis, by using experimental pressure and temperature values immediately prior to expansion (*P*_{max},*T*_{prev}), instead of (*P*_{N},*T*_{min}), the amount of ice formed after expansion can be determined. They do as follows: the minimum temperature after expansion, *T*_{min}, is calculated using (2.2); the corresponding supercooling is attained as Δ*T*_{sc}=*T*_{F}−*T*_{min}; equation (2.18) is used to calculate the latent heat released at atmospheric pressure, and finally we introduced these values in (2.17) to determine the percentage of ice instantaneously formed in HPSF experiments with rapid expansions. We use this simplification in our simulations.

#### (iii) Modelling the mass fraction of ice of the rest of the process

Once we have calculated the amount of ice formed instantaneously after expansion, we need an expression to calculate how the rest of the ice is formed as a function of temperature. As it occurs at atmospheric pressure, we do not have to worry about taking into account pressure effects this time. After expansion, an amount *m*_{i} of ice is instantaneously formed, which can be calculated, as we explained in §2*c*(ii). The temperature is raised to the corresponding freezing point at atmospheric pressure due to the release of latent heat.

*Remark 2.1.* The mass fraction of ice *f*_{s} as a function of temperature can be calculated with the following expression:
2.19where *x*_{0} is the initial mass fraction of solids in the food and *x*(*T*) is the mass fraction of solids in the food (without including the ice generated after the initial reference instant) at temperature *T*.

This can easily be seen by taking *m* as the initial mass of solids, *m*_{w0} as the initial mass of water and *m*_{w}(*T*) as the mass of water in the food at temperature *T*, then
and, therefore,
Next, we derive an equation for *x*(*T*). In Rahman *et al.* (2010), the following extended Clausius–Clapeyron equation is presented to calculate the freezing point depression (*T*_{w}−*T*) of gel, as a function of mass fraction of solids (*x*):
2.20In (2.20), *T*_{w} is the freezing point of water at atmospheric pressure (i.e. 0^{°}C), *T* is the freezing point of the food (^{°}C), *α*_{w} is the molar freezing point constant of water (1860 kg K mol^{−1}), *γ*_{w} is the molecular weight of water (18 kg mol^{−1}), *E* is the molecular weight ratio of water and solids (*γ*_{w}/*γ*_{s}) and *B* is the ratio of unfrozen water to total solids. *E* and *B* are model parameters that Rahman *et al.* (2010) estimate using Statistical Analysis Software non-linear regression.

From (2.20), we can calculate the freezing point depression depending on the mass fraction of solids, but we want to calculate the inverse, i.e. given a certain freezing point temperature *T*, what is the mass fraction of solids at that temperature. From (2.20), it is easy to find an expression for *x* as a function of *T*,
2.21Substituting (2.21) into (2.19), we finally obtain the equation for the mass fraction of ice as a function of temperature,
2.22

#### (iv) Expression for the volume fractions

In our model, we work with volume fractions instead of mass fractions. The relation between them is given in the following remark.

*Remark 2.2.* In a mixture of two phases, the volume fractions (*g*_{l},*g*_{s}) as functions of the mass fractions (*f*_{l},*f*_{s}) are
2.23The density of the mixture can be written as *ρ*=*ρ*_{l}*g*_{l}+*ρ*_{s}*g*_{s} and the relationships between the mass and volume fractions in liquid and solid phases are
As *ρ*=*ρ*_{l}*g*_{l}/*f*_{l} and also *ρ*=*ρ*_{l}*g*_{l}+*ρ*_{s}*g*_{s}=*g*_{l}(*ρ*_{l}−*ρ*_{s})+*ρ*_{s}, after some straight-forward calculations, we have that

As we said in §2*b*, the liquid and solid volume fractions for our model depend on temperature and time. For a classical freezing process, they are considered to be only dependent on temperature, but in a HPSF process we have information *a priori* about when the sample starts to freeze. We know that the sample remains unfrozen (even at subzero temperatures) until the pressure has been completely released (i.e. at *t*=*t*_{p3}), and at that point, there is a percentage of ice instantaneously formed, *m*_{i}, the mass of which we calculate using (2.17) and (2.18). After this point, the rest of the ice, *f*_{s}(*T*), is computed as a function of temperature following (2.22). We consider that the initial mass of solids in the food is the one before the expansion plus the mass of ice instantaneously formed after releasing pressure. Therefore, we take *x*_{0}=*f*_{b}+(*f*_{l}−*f*_{b})*m*_{i}, where *f*_{b} is the mass fraction of solids in the food and *f*_{l} is the mass fraction of liquid in the food. Finally, we have that
2.24where *ε* is the small time interval in which the mass fraction *m*_{i} of ice is instantaneously formed after the pressure release.

### (d) Resulting full model for a high-pressure shift freezing process

We focus on two different situations: HPSF of a big solid-type food (i.e. when the filling ratio of the food sample inside the vessel is much higher than the one of the pressurizing medium) and a small solid-type food (i.e. when the filling ratio of the food is not much higher than the one of the medium).

#### (i) For a big solid-type food

We use system (2.4), but we replace the first equation with (2.12), resulting in
2.25where *g*_{l} is given by (2.24), *P* is given by (2.13), and *C*_{vol}, δ*e*, *k* and *β*_{vol} are defined in §2*b*.

#### (ii) For a small solid-type food

We use system (2.5), but we change the first equation for (2.12) plus the convective term of (2.5), resulting in
2.26We point out that in the cases considered in §2*d*(i,ii) the pressurizing medium used in HPSF processes does not freeze (therefore, *g*_{l}=1 in *Ω*_{P}). For both cases, in order to reduce computational complexity, following the work done in Infante *et al.* (2009) for non-freezing processes, we assume that the thermophysical properties of the food sample are constant (we set them to their mean value in the range of temperature and pressure considered in the process, as done in Infante *et al.* 2009) but different in the unfrozen and frozen states. The thermophysical properties of the steel and the rubber cap remain constant during the whole process (*g*_{l}=1 in those domains). The thermophysical properties of the pressurizing fluid are assumed to be constant for the big sample case, and for the small sample case, where the convection effects are taken into account, we base the model on the Boussinesq approximation. That is, we assume that the coefficients *C*_{p}, *k*, *α* and *η* are constant, and *ρ* is also chosen as a constant value everywhere except in the gravitational force *ρ***g** that appears in the second equation of system (2.26).

## 3. Numerical tests

We consider for the numerical tests the size of the pilot unit (ACB GEC Alsthom, Nantes, France) that was used in Infante *et al.* (2009) and Otero *et al.* (2007). The dimensions of the machine are *L*=0.09 m, *H*=0.654 m, *L*_{2}=0.05 m, *H*_{1}=0.222 m and *H*_{5}=0.472 m (figure 2). We simulate in this section the two cases described in §2*d*. The size and location of the sample and the rubber cap are given by *H*_{3}=0.404 m and *H*_{4}=0.439 m in both cases; *L*_{1}=0.045 m and *H*_{2}=*H*_{1} in the big sample case, and *L*_{1}=0.02 m and *H*_{2}=0.294 m in the small sample case (figure 2). The numerical tests we present are computed in cylindrical coordinates assuming axial symmetry. We use the finite-element method solver COMSOL Multiphysics 3.5a. Velocity and pressure spatial discretization is based on P2-P1 Lagrange finite elements satisfying the Ladyzhenskaya, Babuska and Brezzi stability condition. The time integration is performed using the variable-step-variable-order backward-differentiation-formula-based strategy implemented in the platform. The nonlinear systems are solved using the unsymmetric multi-frontal method for sparse linear systems combined with the stabilization technique Galerkin least squares.

For both cases, we consider agar gel as the solid-type food sample (this gel is a solid food simile that contains 99 per cent water and therefore its properties are taken as those of water). The thermophysical properties for the agar gel, in both unfrozen and frozen states, are, respectively, *ρ*_{Fl}=997 kg m^{−3}, *C*_{p}_{Fl}=4179 J kg^{−1} K^{−1}, *k*_{Fl}=0.613 W m^{−1} K^{−1}, *β*_{Fl}=3.351×10^{−4} K^{−1}, *ρ*_{Fs}=918 kg m^{−3}, *C*_{p}_{Fs}=2052 J kg^{−1} K^{−1}, *k*_{Fs}=2.31 W m^{−1} K^{−1} and *β*_{Fs}=7.97×10^{−4} K^{−1}. The parameters for equation (2.22) are taken from Rahman *et al.* (2010) to be *E*=0.026, *B*=0.050, *α*_{w}=1860 kg K mol^{−1} and *γ*_{w}=18 kg mol^{−1}. Parameters *B* and *E* given in Rahman *et al.* (2010) were estimated for bovine gel, and we are simulating the freezing of agar gel. We did not find in the literature any fitting parameters for agar gel, and that is why we used the ones for bovine gel, as they give good results. For *x*_{0} in (2.22), we need the mass fraction of bounded solids in the food *f*_{b}, which we can calculate using (2.23) with *g*_{b}=0.01, as agar gel contains 99 per cent water.

The thermophysical properties of the steel and rubber cap are *ρ*_{S}=7833 kg m^{−3}, *C*_{p}_{S}=465 J kg^{−1} K^{−1} and *k*_{S}=55 W m^{−1} K^{−1} (for steel) and *ρ*_{C}=1110 kg m^{−3}, *C*_{p}_{C}=1884 J kg^{−1} K^{−1} and *k*_{C}=0.173 W m^{−1} K^{−1} (for rubber). The environment temperature and the heat-transfer coefficient used in all the tests are *T*_{env}=19.3^{°}C and *h*=28 W m^{−2} K^{−1}, respectively. All of these data have been obtained from Fernández *et al.* (2006) and Otero & Sanz (2003). For the two cases, different pressurizing fluids are considered, and for each case, different temperature and pressure conditions are considered. We have considered these exact conditions and pressurizing medium in order to be able to compare our results with experimental published data (Otero & Sanz 2003; Fernández *et al.* 2006).

For the big sample, following the experiments described in Otero & Sanz (2000) and Sanz & Otero (2000), the pressurizing fluid is considered to be a mixture of ethylene glycol and water (75/25, v/v), which has a very low freezing point and therefore does not freeze during the process (Otero & Sanz 2003). The thermophysical properties of this fluid are taken to be constant (as explained in §2*d*(ii), we take their mean values in the range of pressure and temperatures of the experiments) and using data from Guignon *et al.* (2010), we obtain *ρ*_{P}=1127.91 kg m^{−3}, *C*_{p}_{P}=2972.6 J kg^{−1} K^{−1}, *k*_{P}=0.345 W m^{−1} K^{−1} and *β*_{P}=5.655×10^{−4} K^{−1}. As explained in Otero & Sanz (2000) and Sanz & Otero (2000), before the experiment started, the high-pressure pilot unit was tempered to the final subzero freezing temperature (*T*_{prev}) in order to avoid heat loss during freezing of the agar gel. The pressurizing fluid was also kept at the final freezing temperature, i.e. *T*_{0S}=*T*_{0P}=*T*_{prev}. The initial temperature of the food sample and the rubber cap are *T*_{0F}=*T*_{0C}=2^{°}C. Several HPSF processes were carried out at different final subzero temperatures and maximum pressures, from which we have chosen *T*_{prev}=−18^{°}C/*P*_{max}=180 MPa and *T*_{prev}=−21^{°}C/*P*_{max}=210 MPa to simulate. Pressure was applied at a rate of 2.5 MPa s^{−1}. The temperature of the cooling medium, *T*_{cool}, was 0.5^{°}C lower than the subzero final freezing temperature, *T*_{prev}, for each experiment. Sanz & Otero (2000) determine the initial freezing temperature of agar gel from the freezing plateaus of all the experimental processes and on the basis of the recorded data, *T*_{F}=−0.3^{°}C was considered.

For the small sample, following the experiments described in Otero & Sanz (2006), we consider a mixture of ethylene glycol, water and ethanol (40/40/20, v/v/v), as the pressurizing medium. This also has a very low freezing point and therefore does not freeze during the process (Fernández *et al.* 2006). As we explained in §2*d*, all the thermophysical properties are considered to be constants, except for *ρ* in the gravitational force term. The constant values have been taken from Guignon *et al.* (2006) as *ρ*_{P}=1011.77 kg m^{−3}, *C*_{p}_{P}=3042.3 J kg^{−1} K^{−1}, *k*_{P}=0.381 W m^{−1} K^{−1} and *β*_{P}=6.219×10^{−4} K^{−1}. The viscosity has been taken as *η*_{P}=0.02 Pa s. For the density of the fluid as a function of temperature and pressure, we follow Guignon *et al.* (2010), where the volumetric properties for binary mixtures of pressure-transmitting fluids are given, and also equations for calculating these properties for other mixtures from those of their pure components. In this case, as explained in Otero & Sanz (2006), again before the experiment started, the high-pressure pilot unit was tempered to the final subzero freezing temperature (*T*_{prev}) in order to avoid heat loss during freezing of the agar gel. The pressurizing fluid was also kept at the final freezing temperature, i.e. *T*_{0S}=*T*_{0P}=*T*_{prev}. The initial temperature of the food sample and the rubber cap are *T*_{0F}=*T*_{0C}=5^{°}C. Several HPSF experiments were performed at different final subzero temperatures and at various pressures, from which we have chosen to simulate the following: *T*_{prev}=−8^{°}C/*P*_{max}=120 MPa and *T*_{prev}=−20^{°}C/*P*_{max}=210 MPa. The temperature of the cooling medium, *T*_{cool}, was 0.5^{°}C lower than the subzero final freezing temperature, *T*_{prev}, for each experiment. Pressure was also applied at a rate of 2.5 MPa s^{−1}. Otero & Sanz (2006) determine the initial freezing temperature of agar gel from the freezing plateaus of all the experimental processes and on the basis of the recorded data, in this case, *T*_{F}=−0.1^{°}C was considered.

For both cases, as we know the rate at which pressure is applied and the maximum pressure we have to reach, *t*_{p1}, i.e. the time in which we are increasing pressure, is calculated as *P*_{max}/*P*_{rate}. For the sake of simplicity (we could have also used models (2.4) or (2.5) to compute it), we have obtained time *t*_{p2} (the time when pressure is at its maximum and the temperature of the sample has reached the desired freezing temperature, *T*_{prev}) from the experimental data (Otero *et al.* 2000; Otero & Sanz 2006; Sanz & Otero 2000). Finally, *t*_{p3}=*t*_{p2}+2, as we are considering the case of rapid pressure release, that is considered to be in 1–2 s.

### (a) Results

We have simulated the experiments described in §3 using the models explained in §2. The results of solving (2.25) for a HPSF process of a big sample of agar gel are shown in figure 3. These experiments with agar gel and the different temperature and pressure profiles are described in Otero *et al.* (2000) and Sanz & Otero (2000). We do not have the experimental data to compare our simulations with (graphically), so instead we compared the phase-transition times of our simulations with the ones published in Otero & Sanz (2000) and Sanz & Otero (2000). Following Otero & Sanz (2006), we have calculated the phase-transition times (plateau times shown in table 1) as the time span between nucleation and reaching a temperature 5^{°}C below the corresponding initial freezing point at the centre of the sample. We compare the predicted plateau times of our model to the experimental data, and also to the plateau times predicted by other HPSF models (table 1). For the big sample, we compare them with the theoretical times calculated by Otero & Sanz (2000); for the small sample , we compare them with those predicted in Norton *et al.* (2009). The times in table 1 are not given as such in Norton *et al.* (2009), where instead a ‘reduction in plateau time’ (%) was given; this refers to how long the plateau time has been reduced when compared with atmospheric pressure freezing (APF) process. The APF times are given in Otero & Sanz (2006).

In figure 4, the results of solving equation (2.26) for a HPSF process of a small sample of agar gel are shown compared with experimental data, and agree very well. These experiments with agar gel and the different temperature and pressure profiles are described in Otero & Sanz (2006). In that paper, the authors remark that the free convection in the pressure fluid can acquire relative importance, as the pressure medium occupies 88.7 per cent of the total vessel volume. In our model, we include the convective effects, and we simulated exactly the same experiments without including convection effects, i.e. with only conduction effects, and the results (not shown here) were very different and did not agree with the experimental data. Without including convection effects, the small sample takes much longer to cool.

Our model predictions capture the trend in measured plateau times with applied temperature and pressure, and tend to underpredict the data slightly (figures 3 and 4; table 1), although in only one case of the four, the simulated time is not within experimental uncertainty, but the remaining difference is of the order of 0.5 min, thus it should be considered acceptable for a process that takes about 1 h. In any case, we discuss possible reasons for this underprediction in the conclusions. The other published models tend to overpredict the plateau times slightly, but also give quite satisfactory results. However, we remark that Otero & Sanz (2000) calculated the plateau times taking into account the amount of ice instantaneously produced after expansion, and the time required to freeze the sample at atmospheric pressure, but did not obtain them as a result from a heat-transfer model. Norton *et al.* (2009), for the small sample case, proposed a one-dimensional heat-transfer model for the food sample, but they did not model the heat transfer in the pressurizing medium. They assumed a boundary condition of the third class at the surface between the food and the pressurizing fluid, and had to choose the surface heat-transfer coefficient to best fit the experimental curve, i.e. they had a fitting parameter. Also, they did not validate the precooling stage, so their model is validated only from the temperature at which the pressure was released.

## 4. Conclusions

The models described in this paper provide a useful tool to simulate the temperature profile at all points inside both big and small solid food samples going through a HPSF process. The model for solidification is based on the enthalpy formulation at non-constant pressure and volume phase fractions. The models are two-dimensional axisymmetric, and can therefore predict the temperature distribution at any point inside the food sample, not just at any radial component, which happens with one-dimensional models, but also at any height. Convection effects in the pressurizing fluid are considered for the small sample case (and shown to be important).

For the small sample case, the numerical tests agree with the experimental data very well, especially at the centre of the sample. We emphasize that no fitting parameters were used in the simulations. At the surface, there is not an accurate match between the data and the model. However, we consider that the modelling efforts are good enough, given the relatively large experimental errors. It is important to point out that locating a thermocouple exactly at the surface is very difficult, and often generates significant measurement errors from the neighbouring material temperature, which is colder than the food sample. For the big sample case, there is no experimental data to plot against the simulated temperature profile, but the results look very similar to the temperature profiles plotted in Sanz & Otero (2000). When comparing the model plateau times to the published experimental ones, we find that they are, in general, shorter, but in the range of the latter. In all cases, the prediction errors are of the order of minutes, which is not too bad given that the processing time is of the order of hours.

Future improvements, in case experimental measurements become more precise, could potentially be obtained by including the temperature and pressure dependence of material properties, anisotropy effects in the frozen region, and extending our two-dimensional convection simulations to three dimensions.

## Acknowledgements

This work was carried out owing to the financial support of the Spanish ‘Ministry of Science and Innovation’ under projects MTM2008-04621 and MTM2011-22658; the research group MOMAT (ref. 910480) supported by ‘Banco Santander’ and ‘Universidad Complutense de Madrid’; and the ‘Comunidad de Madrid’ and ‘European Social Fund’ through project S2009/PPQ-1551. This publication was also based on the work supported in part by award no. KUK-C1-013-04, made by King Abdullah University of Science and Technology (KAUST). Finally, we would like to thank Dr Pedro Sanz and Dr Laura Otero for providing us with experimental data.

- Received October 13, 2011.
- Accepted April 5, 2012.

- This journal is © 2012 The Royal Society