University of Birmingham Crystallisation in concentrated systems

Water crystallisation in concentrated systems has been studied using a combination of mathematical modelling and experimental work. Two different freezing models have been employed to describe primary and secondary mechanisms (i.e. non-seeded and seed-induced processes, respectively) in sucrose solutions up to 60% (w/w). Differential Scanning Calorimetry (DSC) has been employed to characterise the phase change of the binary water–sucrose system in primary processes, and the kinetic and thermodynamic parame-ters obtained were coupled to the heat transfer equation to obtain the product temperature distribution. A recently developed method has been also employed to measure crystal growth rates in seed-induced crystallisation systems. Simulated results for the secondary crystallisation mechanism were able to reproduce experimentally observed trends for growth rates. An evaluation of the energy consumption during freezing/crystallisation processes has been carried out to assess each mechanism performance (crystallisation will occur at temperatures approximately 20 ◦ C higher in seeded processes) considering different process conditions and product formulations (i.e. solids and air fractions).


Introduction
Freezing is one of the most important processing techniques for manufacturing in the food sector, as well as one of the most regularly used methods to preserve and extend shelflife of foodstuffs (Kiani and Sun, 2011;Zhang et al., 2015). During freezing, the available water within a system is transformed into ice progressively, as the temperature decreases and crystallisation takes place. The process of crystallisation involves two distinct mechanisms: nucleation and crystal growth. The former (nucleation) refers to the establishment of stable crystal nuclei, either in a crystal-free solution (primary nucleation) or by inclusion of previously created crystals (secondary nucleation); the latter (crystal growth) involves growth of these stable nuclei to larger crystals (Kiani and Sun, 2011). The stochastic nature of crystallisation phenomena (Zhang et al., 2015) results in a extremely complex process with many • In frozen vegetables and fruits, for example, the cell tissues could result damaged if too large crystals are formed during the freezing process, affecting both quality and consumer perception (Sun and Zheng, 2006). • In ice cream production, an enhancement of the smooth texture and cooling sensation is sought by targeting a certain crystal size and distribution (Caillet et al., 2003;Petzold and Aguilera, 2009 interconnected and stable porous structure that allows the subsequent drying and further rehydration of the lyophilised products (Ratti, 2001;Petzold and Aguilera, 2009;Mujumdar, 2014).
Freezing is also one of the most energy consuming operations employed in the food sector (Chourot et al., 2003). The energy intensity of the process is in part due to the nature of the raw materials that are being treated. Most food products and goods exhibit moisture contents typically ranging from 40% to 80%. This large water content results in high energy requirements because of the significant latent heat involved in water solidification. Therefore, processing systems with less water to crystallise -and therefore to undergo the energy demanding phase change -will significantly reduce energy consumption during the operation (Roos et al., 2013).
Some investigations have already been undertaken on the thermal study of concentrated systems, in freezing and crystal formation, although those are empirical and difficult to interpret/extrapolate (Ablett et al., 1992;Roos, 1995;Sacha and Nail, 2009). Thus, the implementation of such industrially relevant crystallisation processes still constitutes a challenge in systems characterised by a high solid content.
A good knowledge of the system behaviour, basic notions on nucleation and growth kinetics, together with the study of the effect of formulation and process conditions on the phenomenon kinetics should lead to a better understanding of crystallisation phenomena (Hartel, 2002). In this context, operational models can be employed as useful tools to help in the understanding of the system and, further, in the design and optimisation of more sustainable food manufacturing protocols, as mathematical modelling can provide an accurate description of the system behaviour and process dynamics (Cameron and Hangos, 2001), also in food systems .
A number of works are available addressing freezing in different food systems, varying from analytical solutions to numerical methods (Delgado and Sun, 2001). Among those works are (Pardo et al., 2002;Cogné et al., 2003b;Hamdami et al., 2004;Dima et al., 2014;Kaale and Eikevik, 2014;Bronfenbrener and Rabeea, 2015). Most of them focus on predicting freezing times (Plank, 1913;Salvadori et al., 1996;Heldman, 2007;Pham, 2014), while some others specifically deal with crystallisation phenomena (Hartel, 2002;Kaschiev and van Rosmalen, 2003;Martins et al., 2011) and size and crystal distribution (Miyawaki, 2001;Kiani and Sun, 2011;Lian et al., 2006;Arellano et al., 2013). In spite of all these modelling efforts, there is still a gap covering the mathematical description and prediction of water crystallisation in food systems with a high solid content.
In here, the influence of high solid concentrations has been studied, and two mathematical models (Lopez-Quiroga et al., 2015) employed to study which of the crystallisation mechanisms (i.e. primary and secondary nucleation) occur during freezing processes of concentrated systems. The first model assumes primary nucleation, and characterises the supercooling required in this system. The second one studies a seeded system (also with supercooling) in which crystal growth occurs onto an existing nucleus (i.e. secondary nucleation). Sucrose solutions with a high solid content (up to 60 %, w/w) have been chosen as representative concentrated systems. This numerical study has been supported by thermal analysis to characterise the water-sucrose system and to obtain the required latent heat and freezing point data of the different Fig. 1 -Non-isothermal DSC thermogram for the crystallisation of the 60% (w/w) sucrose solution cooled to −70 • C at 1 • C/min. The arrow points at the beginning of the phase change peak, which total area has been employed to obtain the value of the phase change enthalpy H (J/g).
formulations. A simple -yet useful -thermodynamic-based study on how processing products with a high solid content could reduce the energy consumption during the manufacturing step completes the work presented here.
The manuscript is organised as follows: in Section 2 the experimental methods employed in this work (calorimetric analysis and seeding-induced experiments) are detailed. Next, the models used to describe crystallisation phenomena and the associated energy consumption are presented in Section 3. The validation of those models as well as the results obtained for growth rates and mean crystal size are presented in Section 4, where the evaluation of the energy consumed during the freezing/crystallisation processes is discussed too. Final remarks are outlined in Section 5.

Sample preparation
Sucrose was mixed with distilled water using a magnetic stirrer under heating (reaching a maximum temperature of 50 • C) to prepare solutions of concentrations ranging from 20% (w/w) up to 60% (w/w). These solutions were then cooled at room temperature and used within the following 24 h.

Primary nucleation and growth: DSC
Differential Scanning Calorimetry (DSC) is a thermal technique widely employed in the analysis of materials undergoing phase changes (Ablett et al., 1992;Sacha and Nail, 2009;Cogné et al., 2003a;Málek, 2000). In this work, DSC was used to characterise water crystallisation kinetics and thermal properties of the sucrose aqueous solutions studied for a range of concentrations up to 60% (w/w) of solids. For each DSC experiment, approximately 5 mg (exact weight recorded) of the system under investigation were added in a standard aluminium sample pan, which was subsequently sealed. The sample and reference pans -the latter containing air -were placed in the DSC device (Perkin Elmer DSC 8000) and were cooled down and scanned from 20 • C to −70 • C at different rates (1 • C/min and 6 • C/min), then held at −70 • C for 3 min and finally heated from −70 • C to 20 • C employing the same rates. Fig. 1 shows the DSC thermogram obtained for the 60% sucrose (w/w) sample. No seed crystals  were added to the system: nucleation will occur at some supercooling level, either in the bulk or at the DSC pan wall. The figure presents the characteristic crystallisation peak, and it also reports the onset freezing temperature T onset freezing ( • C) of the process.
The phase change enthalpy, H (J/g), for each formulation was calculated from the peak area of the corresponding thermograms and used to estimate the free water in the sample, w free (%), assuming that the bound water does not freeze during the experiment and only the free liquid phase contributes to the H value (Cogné et al., 2003a): where the freezing enthalpy corresponding for a pure water system is represented by H * in Eq. (1). Table 1 summarises the phase change enthalpies, free water content, and onset freezing temperatures of all investigated systems as determined by calorimetric experiments (1 • C/min cooling rate). Data show that the energy required for the phase change, H, decreases as the solid content of the system increases, indicating a significant potential in reducing total manufacturing energy expenditure by processing concentrated systems.
From the data collected through the DSC experiments, the ice fraction ˛(T) can be obtained as a function of temperature using the following integral expression (Málek, 2000): Here (W/g) is the experimentally measured specific heat flow at temperature T ( • C) andq ( • C/s) is the cooling rate applied in the calorimetric analyses. Fig. 2 shows both the ice fraction curve recovered from the DSC data shown in Fig. 1 (60% (w/w) sample as illustrative example) together with the counterpart numerical expression fitted as a temperature function. The cooling curves have been employed in the ice fraction calculations to include supercooling effects in the system. The Least Squares method has been employed for the fitting of all the samples analysed in this study. Data listed in Table 1, as well as the fitted ice fraction function, have been further used in the mathematical models presented in Section 3.

Seed-induced experiments: secondary nucleation and growth
A series of secondary crystallisation (seeding) experiments have been conducted in a Linkam LTS 120 Peltier stage employing sucrose solution samples with a high solid content (60%, w/w). The seeding experiments consisted of placing a droplet of solution (100 L) on the Peltier stage at an initial temperature of 10 • C. Next, the droplets were cooled down at 1 • C/min, reaching final temperatures of −16 • C, −18 • C, −20 • C, −22 • C and −24 • C. As would be expected from the DSC experiments, no crystallisation occurred at these temperatures. Droplets were then held at these constant temperature conditions until the completion of the experiments. Once the desired supercooling degree was achieved, a seed (ice) of 3 L in volume, kept at the same temperature than the supercooled solution, was manually placed into the droplet sample and snapshots of the crystal growth process were recorded in real time at suitable intervals. Those snapshots were later analysed using an open source image processing software (ImageJ). The area of the expanding solid phase was measured with time and from the data obtained crystal growth rates (m/s) were calculated (Wang et al., 2014). Fig. 3 presents indicative sequential images showing the expansion of the frozen area recorded during a seeding experiment conducted at −22 • C.

Modelling of crystallisation phenomena
To help in the understanding of crystallisation processes in concentrated systems, two different models, both based on first principles, were employed in this study. The models illustrate the different possible mechanisms for crystallisation, and characterise the conditions under which they occur.

Model for primary nucleation and growth
The first model employed describes primary nucleation and subsequent crystal growth mechanism. Information provided by the DSC experiments is used to characterise the crystallisation kinetics under a range of cooling conditions and concentrations. For modelling purposes, the thermal problem defined describing the cooling process was coupled to crystallisation kinetics by means of an extra heat source representing the energy involved in the water phase change (Hu and Argyropoulos, 1996;Chégnimonhan et al., 2010). This model also considered the addition of air in the product formulation, allowing the assessment of the aeration effect on the formation and growth of crystals.

Thermal properties: sample without air
The thermophysical properties considered in the simulation have been defined by employing a mixing rule between liquid and solid state properties, weighted by fitted ice crystal fraction ˛(T) introduced in Section 2.2: where k, and C p are the heat conductivity, density and heat capacity, respectively. In Eqs.
(3)-(5) the product sample is noted by the subscript prod, while subscripts l and s represent liquid and solid states of the sucrose solution, respectively. For a sample without air, the liquid state l was considered a continuous phase formed by the water-sucrose binary system, while solid state s consisted of a solid phase (ice) dispersed in a continuous fluid phase made of dry matter (sucrose) and unfrozen (bound) water (Cogné et al., 2003a). The densities for each state were defined as: The heat capacity for the non-aerated system was described by a piecewise function dependent on the temperature value (Chen, 1985;Cornillon et al., 1995;Cogné et al., 2003a): C ps = C p (T < T f ) = 1547 + 1254(x sucrose + x bound ) In the above correlation, M is the components molecular weight, and T * f the freezing temperature of pure water. Temperatures T and T * f must be considered in Kelvin in the correlation. Continuing with the description of the sample without air, the thermal conductivities for both states were modelled as a parallel heat resistance: where Â j represents the volume fraction of each component j in the corresponding continuous phase i = l, s: The data employed in Eqs. (6)-(10) for sucrose, water and ice intrinsic properties can be found in Table 2.

Thermal properties for the aerated sample
The effects of aeration in the product were also considered in this first modelling scenario. Liquid and solid states with air were modelled using correlations commonly used in food systems, and the values found were employed in Eqs. (3) and (4) as inputs. The addition of air was not supposed to have any influence on the values of the heat capacity, so Eqs. (7) and (8) were not modified. The aerated product was considered a dispersion of the air phase (bubbles) in a continuous matrix (liquid or solid product) as defined in Section 3.1.1. According to this assumption, the Maxwell-Eucken model was chosen to define the apparent heat conductivity of the product with air (Carson, 2006;Cogné et al., 2003a): where the thermal conductivity of the air is taken as k air = 0.022 (W/m • C) and is the air volume fraction. In accordance, the density of the aerated sample -for both solid and liquid states -was obtained as follows:

Model formulation for primary nucleation and growth
As briefly introduced before, for the description of primary ice crystal formation the heat transfer mechanism has been coupled to both nucleation and crystal growth kinetics (Hu and Argyropoulos, 1996;Chégnimonhan et al., 2010). Thus, the primary crystallisation of a binary solution when convection effects can be neglected -as is the case of the sucrose solution considered -will be governed by the Fourier equation together with an extra heat source as follows: where the term Q c = ice H(∂˛/∂t) encloses the energy required for the phase change, as well as the kinetic characterisation via the fitted ice fraction function ˛(T). For the solution of this thermal problem a one-dimensional domain with length L was considered. Accordingly to this geometry description, the following boundary conditions were imposed: T c being the temperature of the cooling surface, always below the phase change temperature T f . The system is closed by the corresponding initial conditions.

Seed-induced nucleation and growth modelling
The mechanism of secondary crystallisation is defined as the growth of an ice crystal from a pre-existing nucleus: the so-called seed. This mechanism can be mathematically described by following a classical Stefan formulation for phase change problems (Crank, 1984). Consequently, the model employed considers both solid and liquid phases existing simultaneously and separated by the moving freezing front, as schematically represented in Fig. 4. In this way, and after assuming an initially small frozen domain representing the seed, the crystal growth mechanism can be simulated.
To characterise the diffusion mechanism across the sample (seed and initial liquid solution exhibit different concentrations, with no sucrose diffusion in the seed) Fick's second law has been employed: where c i is the concentration and D e i the effective diffusivity of the liquid (l) and solid (s) water states. Similarly, the heat transfer throughout the material has been described by the Fourier equation: From the mass and energy balances across the moving freezing front, the corresponding boundary conditions at r = S(t) were defined as follows (Crank, 1984;Illingworth and Golosnoy, 2005): Complementary, the continuity condition for the temperature across the moving front can be written as follows: The effect of concentration in the freezing temperature has been taken into account by the empirical expression for the freezing-point depression, T, in sucrose solutions (Roos, 1995): The following Dirichlet and Neumann boundary conditions for the temperature and concentration are considered at the edges of the sample r = 0 and r = R: In Eq. (21), T c is noting the temperature of the cooling boundary. Finally, suitable initial conditions were defined to close the system. Those considered that initially both seed and liquid phase were at the temperature of the Peltier stage (i.e. the value of the cooling boundary T c ).

Energy consumption during crystallisation processes
The total energy consumption during a crystallisation process can be calculated as the addition of two different contributions. Firstly, the energy required for changing the sample from its initial storage temperature T ini ( • C), to the subcooled and final temperature desired for the product, T fin ( • C). This first contribution is defined as Q 1 (J/m 3 ): The second contribution is the energy consumed during the phase change, which will be referred to as Q 2 (J/m 3 ), and equals the volumetric latent heat of solidification: Finally, the total energy consumption during the freezing process will be defined as the addition of the contributions: Q total = Q 1 + Q 2 .

Model validation
For validation purposes, the operational models presented in Section 3 were used to simulate the behaviour of a pure water system undergoing freezing. To simulate such scenario, the secondary crystallisation model presented in Section 3.2 was simplified to the thermal problem formed by Eqs. (16) and (18) together with the external boundary conditions for temperature imposed by Eq. (21). For the primary nucleation model described in Section 3.1, the freezing front position was obtained in a post-processing stage by means of the computed ice crystal fraction ˛(T), which for this validation case was defined as a Heaviside step function with a transition interval ıT around T * = 0 • C to represent the change in the physical state of water. The validation simulations for both models were implemented and solved by employing a commercial Finite Element (FE) software (COMSOL Multiphysics). For the seeding-induced model, the FE method was applied together with an Arbitrary Lagrangian Eulerian (ALE) algorithm (Donea et al., 2004) to track the advance of the phase change front. A fixed grid implementation was employed for the solution of the primary ice crystal formation model. The accuracy of the proposed models has been confirmed by the values of the root mean square error (RMSE) obtained when the advance of the freezing front (for both models) and the freezing front positions provided by the similarity solutions available for these freezing problems (Crank, 1984;Hu and Argyropoulos, 1996;Alexiades and Solomon, 1993) were compared. Such error has been quantified as RMSE primary = 7.32 × 10 −5 m for the primary model, and RMSE seeding = 1.21 × 10 −4 m, for the seeding model.

Modelling results for the primary nucleation and growth model
The freezing model presented in Section 3.1 for primary nucleation and growth formed by Eqs. (13) and (14) was solved for a 60% (w/w) sucrose solution employing the numerical procedure explained in Section 4.1. To ensure the completeness of the freezing process, and thus of the crystallisation phenomenon, the cooling boundary was fixed at T c =−50 • C. The model reflects the DSC experiments, as crystallisation takes place at −40.51 • C, a value 20 • C lower than the temperatures identified in the seeding experiments (see Section 2.3). This clearly demonstrates the differences between the two nucleation regimes.
The effect of including air in the product formulation was considered for three different air fractions: = [0.1, 0.2, 0.3]. The simulation results for the ice crystal fraction obtained under different aeration conditions are presented in Fig. 5(a). This figure shows how the freezing front in the system is retarded by the presence of air: i.e. the freezing front does not penetrate as far into the sample. Accordingly, samples with higher air fractions exhibit higher temperature values, as can be seen in the temperature distributions presented in Fig. 5(b). This is a consequence of the insulating effect of the air within the bulk product properties, which reduces the effective thermal conductivity and acts as an impediment for heat transfer through the sample.
The freezing front velocityṠ(t) was computed from the simulated temperature and crystal fraction distributions, and employed to obtain a quantitative estimation of the mean ice crystal size d p (m) (Miyawaki, 2001): The required values for the water diffusivity D w (m 2 /s) have been taken from (Pruppacher, 1972). Fig. 6 shows the estimated mean crystal sizes as a function of the air fraction for each simulated product formulation (from 20% up to 70% (w/w) sucrose solutions). Results show that larger sized crystals are formed in formulations with higher air fractions, and thus with delayed heat transfer through the product, which seems to be in accordance with the effect of slow cooling rates on the growth rate (Petzold and Aguilera, 2009).

Modelling results for the seed-induced crystallisation process
The Stefan model formed by Eqs. (15)-(22) was solved considering a range of cooling temperatures T c = [−16, − 18, − 20, − 22, − 24] • C for a 60% (w/w) sucrose solution by implementing the Finite Element (FE) method in an adaptive grid, as detailed in Section 4.1. To simplify the model, it was assumed that there was no sucrose present in the solid phase, so that the liquid state became more concentrated as freezing progressed. Accordingly, the thermophysical properties in the liquid phase were dependent on sucrose concentration and temperature. A constant value of the mass diffusivity in the liquid was considered in the simulations, which has been estimated from the seeding experimental data by using the expression of the mass diffusion length (Crank, 1975). As a result of the numerical simulations, the spatial distributions for both the temperature and concentration variables were obtained, together with the position of the phase change boundary at each time step. The advance of the freezing position is illustrated in Fig. 7(a), while the temperature distribution of the sample at three different times is presented in Fig. 7(b). The results correspond to a simulated seed-induced process considering T c =−24 • C. The corresponding front temperature, as well as the sucrose concentration at the freezing front are depicted, respectively, in Fig. 7(c) and (d). In these last two graphs, the effect of freeze concentration can be observed: as the solid phase is growing and the solution in getting more and more concentrated -see Fig. 7(d) -the temperature of the phase change (i.e. the front temperature) goes down (see Fig. 7(c)) obeying the equilibrium curve described by Eq. (20). Fig. 8 shows a comparison between simulated and experimentally measured growth rates for the different cooling conditions modelled is presented. Although numerical simulations correspond to a crystallisation process where equilibrium conditions were assumed for the freezing front, the simulated results obtained overall reproduced the observed trends for growth rates during the isothermal seed-induced experiments detailed in Section 2.3.

Energy consumption evaluation
The effect of increasing the sucrose content on the energy requirements during solidification is reflected in the value of Q 2 (J/m 3 ), as defined in Section 3.3. This was calculating employing the freezing enthalpy H obtained through the DSC experiments detailed in Section 2.2. Those enthalpy data (presented in Table 1) show the significant reduction of the energy requirement for the phase change as the solid content increases.
The energy requirements corresponding to the whole crystallisation process, for both primary and secondary mechanisms were calculated by Eqs. (23) and (24).
In accordance with Section 3.1 and Section 3.2, the cooling temperatures there defined (T c =−50 • C and T c =−24 • C for primary and secondary cases, respectively) were chosen as final product temperatures. For the two crystallisation scenarios devised, the initial temperature considered was T ini = 20 • C. For the evaluation of the energy consumed during primary ice crystal formation, different air fractions were considered. Table 3 shows the results obtained for the primary crystallisation process when a range of air fractions = [0.1, 0.2, 0.3] was included in the formulation of the product. The insulating effect of air is again reflected on the energy input required to freeze the sample, as the volumetric heat capacity is strongly dependent on the apparent density of the product, which decreases for the larger air fractions. As consequence, a reduction on the energy consumption was found due to the addition of air to the sample formulation.

Seed-induced vs. no-seed crystallisation
The values obtained for the case of the seed-induced crystallisation phenomenon for the 60% (w/w) sucrose sample were Q seeding 1 = 87.77 J/m 3 leading to Q seeding total = 124.90 MJ/m 3 . When the seed-induced process was compared with the primary spontaneous one (i.e. no seeding of the system), the advantages of inducing the crystal formation through the addition of small crystals were revealed, as a significantly smaller amount of energy (≈20% less energy) was required. The comparison between seeded and no-seeded processes was carried out considering non-aerated samples, being the main factor leading to this result the difference in the temperature at which crystallisation takes place, lower for non-seeded process ( 20 • C less) as experimentally observed (Wang et al., 2014).

Conclusions
Two different models describing crystallisation phenomena (both primary and secondary mechanisms) were employed to predict growth crystal rates and mean crystal sizes in a concentrated aqueous sucrose solution. The models were also used as a tool to assess energy consumption during freezing (crystallisation) processes. A classical Stefan problem formulated in terms of coupled heat and mass transfer has been employed to represent a secondary (i.e. seed-induced) crystallisation process. Numerical simulation results were able to reproduced realistic trends for growth rates. Additionally, estimate kinetic information obtained from DSC experiments was coupled with a PDE-based heat transfer model to simulate primary ice crystal formation. This model has provided realistic estimations of the mean crystal size resulting from the freezing process in both aerated and no aerated product formulations. Simulations also showed that larger crystals were formed when higher air fractions were considered in the product formulation. Finally, an evaluation of the energy requirements of freezing in concentrated systems based on the previous models was presented, too. A significant reduction on the energy demand during freezing was found in systems with reduced water content. In addition, the advantage of seeded processes over non-seeded ones (i.e. secondary over primary crystallisation) was revealed too, as well as the favourable effect of aeration in terms of energy reduction for primary crystallisation processes. Overall, this modelling approach demonstrates that process designers need to consider modification of the product formulation -i.e. increase of solid and air contents -and/or the use of seed-induced processes to achieve a significant reduction on the energy consumption of the freezing/crystallisation operations.