- Email: [email protected]
Coupling population balance model and residence time distribution for pilot-scale modelling of β-lactoglobulin aggregation process
Journal of Food Engineering 177 (2016) 31e41
Contents lists available at ScienceDirect
Journal of Food Engineering journal homepage: www.elsevier.com/locate/jfoodeng
Coupling population balance model and residence time distribution for pilot-scale modelling of b-lactoglobulin aggregation process Nicolas Erabit a, b, c, Fatou Toutie Ndoye a, *, Graciela Alvarez a, Denis Flick b, c a
IRSTEA Refrigeration Process Engineering Research Unit, Antony, France AgroParisTech, UMR1145 Ing enierie Aliments Proc ed es, Massy, France c INRA, UMR1145 Ing enierie Aliments Proc ed es, Massy, France b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 1 September 2015 Received in revised form 2 November 2015 Accepted 20 December 2015 Available online 23 December 2015
This study presents a hybrid modelling approach for the prediction of the thermally-induced aggregation of whey protein solutions in a continuous process. Modelling a continuous process needs to account the fluid flow phenomena because the different fluid parcels are not subjected to the same time/temperature/shear history. The modelling approach developed combines a residence time distribution (RTD) model and a colloidal aggregation model based on Population Balance Equations (PBE). This combined model has been compared, on one side, to a Population Balance Model assuming plug flow, and on the other side, to experimental data. A 6% b-lactoglobulin solution with different CaCl2 concentrations (5.0, 6.0, 6.4, 6.6 and 7.0 mM) was used for experimental investigations. Experimental heat treatments were carried out at 87 C with a holding time of 49 s using tubular heat exchanger and holder. The integration of residence time distribution appears to be very important when the product is subject to large transformations (formation of large aggregates > 10 mm). Comparison of predictions obtained from both models to the experimental data shows that the coupled PBE/RTD model is able to predict the characteristics of the suspensions (aggregate size and viscosity) and successfully integrates the impact of calcium upon aggregation process, while the model assuming plug flow is suitable only when aggregation is low (for low calcium concentrations). © 2016 Elsevier Ltd. All rights reserved.
Keywords: Population balance model Residence time distribution Whey protein ta-Lactoglobulin be Aggregation
1. Introduction A typical dairy process involves multiple continuous operations, among which heat treatment is the most common processing step. Thermal treatment is performed with the aim of stabilizing (microbiologically and physically) or modifying the structure of the components. Considering the specific case of structuring of dairy products, whey protein (WP) aggregation process is still a challenge in term of product quality reproducibility. Thermal and mechanical operations applied to WP induce denaturation and aggregation of protein molecules leading to particles with highly sought functionalities such as foaming, emulsifying, gelling and thickening agents. These properties of WP aggregates are widely governed by the aggregates size distribution (ASD) depending on environmental factors such as temperature, pH and ionic strength. It is necessary to implement consistent process systems engineering approaches in
* Corresponding author. E-mail address: [email protected] (F.T. Ndoye). http://dx.doi.org/10.1016/j.jfoodeng.2015.12.013 0260-8774/© 2016 Elsevier Ltd. All rights reserved.
order to provide final product with the desired target functionalities according to process variables. A model-based system approach can deliver the required knowledge for process understanding, evaluation and analysis with reduced time and cost for process design. Modelling approaches commonly applied to heat induced WP aggregation are based on reaction kinetic methods (Dannenberg and Kessler, 1988; Anema and McKenna, 1996; Tolkach and Kulozik, 2007). Such an approach can be useful to provide global knowledge about the native proteins and aggregates concentrations, but it fails to depict all the physical phenomena involved in the aggregation mechanism. Moreover, the reaction kinetic based model is unable to describe the evolution of aggregates size along the process. A modelling strategy that has the potential to overcome these limitations is the application of population balance model (PBM) to the aggregates formation. There are relatively few published works in which PBM is used to describe ASD of heat induced WP aggregates (Elofsson et al., 1996; Arosio et al., 2012; Chantoiseau et al., 2012; Erabit et al., 2015). The recent work of Erabit et al. (2015) was devoted to the development of a new
32
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
Nomenclature
Latin symbols Radius of monomers (m) a0 ai Radius of aggregates in size class i (m) aj Radius of aggregates in size class j (m) ca unfolding constant (K) C Concentration of native monomers (kg m3) Cmon Concentration of monomers (kg m3) Cunf Concentration of unfolded monomers (kg m3) di Diameter of aggregates in size class i (m) Df Fractal dimension () Dv(0.5) Median volume diameter (m) E(t) Residence time distribution function (s1) E(q) Dimensionless residence time distribution function () h Proportionality factor () I Total number of size classes () k Constant of Boltzmann (m2 kg s2 K1) kij Total aggregation kernel implying size classes i and j (m3 s1) kperi,ij Perikinetic aggregation kernel implying size classes i and j (m3 s1) kortho,ij Orthokinetic aggregation kernel implying size classes i and j (m3 s1) lh Heat exchanger length (m) m, n Number of compartments in the RTD compartment model () N Number concentration of aggregates (m3) Ni Number concentration of aggregates in size class i (m3) Nj Number concentration of aggregates in size class j (m3) Nmon Number concentration of native monomers (m3) N1 Number concentration of unfolded monomers (m3) p Upper size class of aggregates involving hydrophobic interactions qi Number of monomers in aggregates of size class i () t Time (s)
modelling approach accounting physico-chemical and thermomechanical conditions to predict b-lactoglobulin (b-lg) aggregates suspensions properties such as ASD, residual native fraction and viscosity. The modelling considers the time/temperature/shear history and Population Balance Equations (PBE) based on FuchsSmoluchowski theory (Smoluchowski, 1917). The novelty of the modelling method is related to the adaptation of the stability factor values according to the specific mechanisms involved in b-lg aggregation. The model gave good results when applied to b-lg aggregation in a batch process simulator at laboratory scale, but cannot be applied directly to a continuous heat treatment process. The main reason is that phenomena related to fluid dynamics were not accounted in this model. However, it is known that WP aggregation under continuous thermal process in heat exchangers involves three aspects: fluid flow, heat transfer and product transformation (denaturation and aggregation) (Spiegel and Huss, 2002; Simmons et al., 2007; Nicolai et al., 2011). The structural transformations results from the combination of these three aspects. Moreover, these transformations induce significant rheological changes, and consequently, important modifications of the
t tr T Theat Ts Ttrans
V V_ W W1 W2 W3 z
Mean residence time (s) Fluid parcel residence time (s) Temperature (K) Maximal heating temperature (K) Saturation temperature of steam (K) Temperature of transition between unfolding limiting and aggregation limiting steps during the heat induced aggregation of b-lg (K) Volume (m3) Volumetric flow rate (m3 s1) Stability factor () Stability factor for monomers and aggregates of size classes i 3 () Stability factor for aggregates of size classes 4 i p () Stability factor for aggregates of size class pþ1 i I () Axial coordinate
Greek symbols Fraction of monomers in the unfolded state () Aggregates volume fraction () g_ Local shear rate (s1) m Dynamic viscosity of the suspension (Pa s) ms Dynamic viscosity of the continuous phase (Pa s) q Dimensionless time () q0 Breakthrough time () t Geometrical residence time in one tube (s) tth Thermal characteristic time of the heating section (s)
a f
Subscripts co cooling exp experimental ho holding i inlet o outlet th theoretical T total
fluid flow behaviour especially near the heat exchangers walls (Chantoiseau et al., 2012; Ndoye et al., 2012). Additionally, it is recognized that the different fluid parcels flowing through a continuous processing unit are not subjected to the same timetemperature-shear profiles because of the residence time distribution (RTD) (Dannenberg and Kessler, 1988). Development of a realistic model that can represent WP aggregation at continuous pilot or industrial scale needs the implementation of a modelling approach accounting fluid dynamics framework, through RTD modelling for example. A hybrid method combining PBM to fluid flow description can potentially provide the necessary supplement understanding for the WP aggregation continuous process. Several hybrid models coupling PBE to dynamical histories have been developed for particulate continuous process such as crystallization (Wynn and Hounslow, 1997; Lian et al., 2006; Arellano et al., 2013), granulation (Sen et al., 2012; Sen and Ramachandran, 2013; Lee et al., 2015), grinding (Weller et al., 2000; Bilgili and Scarlett, 2005), polymerization (Prasetya et al., 1999; Vale and McKenna, 2005), flocculation (Heath and Koh, 2003; Tanguay et al., 2014) and protein breakage (Zumaeta et al.,
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
2005, 2010; Bas¸ et al., 2009). In authors knowledge, the only study coupling fluid flow characteristics and PBM applied to heat induced WP aggregation is that of Chantoiseau et al. (2012). The authors developed a numerical approach that combines computational fluid dynamics (CFD) for solving the fluid flow and heat transfer coupled problem, with a PBE for evaluating b-lg particle aggregation inside a straight section of a tubular heat exchanger. The present study aims to model b-lg aggregation in an entire heat treatment processing unit including heating, holding and cooling sections, by coupling fluid dynamics framework to PBE. Fluid flow heterogeneities are taken into account through RTD measurements and modelling of a continuous laboratory pilot scale thermal process as performed in a previous work by Ndoye et al. (2012). This RTD model is combined to the PBM developed by Erabit et al. (2015) for the thermally-induced b-lg aggregation based on laboratory scale experiments in a batch process simulator. Such a hybrid modelling approach is compared to a PBM approach assuming plug flow. The predicted results were compared to experimental data characterizing the suspensions (ASD and residual native fraction) obtained during the aggregation process of a blg solution in a continuous heat treatment at the laboratory pilot scale. 2. Experimental approaches 2.1. b-lg heat treatment and characterization 2.1.1. Protein solution A 6% b-lg solution was prepared using industrial powder (92% of b-lg, 5.2% of moisture, 2.5% of ash and 0.3% of fat) and deionized water at 40 C under mixing. The protein concentration used is adequate to form a significant number of aggregates without gelling (Erabit et al., 2014). Different concentrations of Ca2þ were obtained by adding a CaCl2 1 M solution to the protein solution. The prepared solution was placed in a water bath at 40 C for 2 h in order to have complete rehydration of the b-lg powder. The solution was then kept at 4 C and used within a day to prevent any destabilization. Five concentrations of CaCl2 have been tested: 5.0, 6.0, 6.4, 6.6 and 7 mM. The pH of the solution was equal to 6.8 without any significant differences between the different concentrations of CaCl2. 2.1.2. Continuous thermal processing unit An ultra-high temperature (UHT) laboratory pilot scale plant (Armfield FT174X, Ringwood, UK) was used to carry out continuous heat treatments. A detailed description of the experimental setup can be found in Ndoye et al. (2012). The product was fed by a volumetric pump at a flow rate of V_ ¼ 20 L h1 in order to maintain a laminar flow. The applied heat treatment and bulk temperature profile are described on Fig. 1. A preheating at 60 C was first performed using a helical tube in a water bath. The product is then heated at T ¼ 87 C in a tubular heat exchanger made of eight concentric double tubes, using steam at saturation temperature Ts as heating media. After heating inside the heat exchanger, the product passes through a holding section composed of helical tube with adaptable volume in which heat losses lead to a slight decrease of the product temperature to Tho. The last section of the processing unit is the cooling heat exchanger similar to the heating section. The product is cooled down to a temperature Tco ¼ 4 C using cold water. Sampling was done at the outlet of the cooling section for further analysis. 2.1.3. Suspension characterization 2.1.3.1. Residual native fraction. The residual native fraction t is the fraction of the proteins which remain in the native state during
33
heat treatment but also of the proteins which unfold but do not aggregate and come back to the native state during cooling. This fraction was determined by HPLC as described by Tolkach and Kulozik (2007). An acetate buffer at pH 4.6 was used to dilute samples (1/25th) and to promote agglomeration of the denaturated proteins. Agglomerated materials was then removed with syringe filters (Ø 0.2 mm). 2.1.3.2. Granulometry. The volume size distribution of the suspension of aggregates was determined using a Laser granulometer Mastersizer ZS (Malvern Instruments Inc., Worcester, UK) for particles from 0.1 to 1000 mm. The measurement principle includes a dilution of samples before performing analysis. The required level of dilution is directly determined by the granulometer depending on the turbidity of the suspension. 2.2. Residence time distribution measurement The experimental characterization of the RTD were previously reported by Ndoye et al. (2012). The RTD measurements were performed in the UHT laboratory scale plant described above using a 6% b-lg solution with 5 mM CaCl2. Conditions permitting no WP aggregation (T ¼ 60 C) and allowing WP aggregation (T ¼ 87 C) were explored at 20 L h1 and 49 L h1. Experiments were carried out using a colorimetric method with Methylene Blue as tracer. A 10 mL sample of 6% b-lg solution coloured with MB (at concentration of 0.01 g L1 for no aggregation conditions and 0.02 g L1 for aggregation conditions) was pulse injected at the entrance of the heat exchanger or of the holding tube according to studied section. The protein samples were continuously collected at the outlet of the device and analysed by the use of a spectrophotometer (Beckman Coulter DU 730, Fullerton, USA). The age distribution function E(t) was used to characterized the RTD and mean residence time was determined for each conditions explored. 3. Modelling First a model of unfolding-aggregation of b-lg based on PBE and assuming plug flow along the heating, holding and cooling sections is described. Then the PBE model is coupled to an empirical RTD model in order to take into account that the fluid parcels of protein suspension flowing throughout the heater, the holding tube and the cooler do not have the same time-temperature-shear history. 3.1. Population balance modelling approach The complete mechanism of thermally induced aggregation of
b-lg is described as a two steps reaction composed of a reversible unfolding of monomers followed by an irreversible aggregation of unfolded proteins. The modelling framework used in this work combines the unfolding step described through an unfolding ratio a and the aggregation step described through PBE. A detailed description of the model development was recently reported by Erabit et al. (2015). The model was validated on the basis of laboratory batch experiments. This section describes briefly the main points of the modelling approach. 3.1.1. Unfolding mechanism e thermo-dependent equilibrium The unfolding phase of the thermally-induced aggregation of blg was modelled based on a temperature-dependant equilibrium between native and unfolded WP molecules as described by Tolkach and Kulozik (2007). The fraction of unfolded proteins among the protein monomers is defined as the unfolding ratio a ¼ Cunf/Cmon. This unfolding ratio is function of the temperature according to Eq. (1):
34
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
Fig. 1. Heat treatment carried out on the pilot plant: (a) schematic representation of the heat treatment sections and (b) applied bulk temperature profile.
8 1 1 > < ln a ¼ ca Ttrans T > : a¼1
for T Ttrans for T > Ttrans
(1)
where ca is the unfolding coefficient and Ttrans is the transition temperature at which unfolding or aggregation is the limiting step of the overall aggregation process of b-lg.
3.1.2. Aggregation mechanism - population balance model The population balance model was developed based on Smoluchowski PBE which describes irreversible aggregation accounting the particles collision mechanism through kernels: perikinetic aggregation for collision by Brownian motion and orthokinetic aggregation for collision by shear flow. The population balance model is described by the following equations: I1 X dNmon ¼ N1 k1j Nj dt j¼1
(2)
i1 X dNi > 1 1 X ¼ kij;j Nij Nj Ni ki;j Nj 2 j¼1 dt j1
(3)
where N1 ¼ a:Nmon is the number concentration of the unfolded monomers related to the number concentration of the total monomers (native and unfolded) Nmon by the unfolding ratio a. Aggregation kernel kij accounts both perikinetic and orthokinetic collisions, and the stability factor W represents the inverse of the probability that a collision between particles actually leads to an aggregate (Eq. (4)).
kij ¼
kperi;ij þ kortho;ij W
(4)
A novelty of the modelling strategy reported by Erabit et al. (2015) is the adaptation of the stability factor to three different mechanisms specific to b-lg aggregation. This leads to the definition of three stability factors according to the interactions involved: (i) W1 for aggregates between unfolded monomers (i ¼ 1) and small oligomers (I 3) through disulphide bonds; (ii) W2 for collisions involving hydrophobic interactions and larger aggregates (4 < i 8000) ensuing; (iii) W3 for collisions between large aggregates (i > 8000).
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
When aggregates of two different size-ranges collide, the stability factor taken into account is that of the smallest aggregate. It is assumed that the value of the stability factor does not depends on the type of collision involved (orthokinetic or perikinetic). Experimental measurements of aggregates size showed that there are few particles larger than 100 mm. This is certainly due to the difficulty for two large particles (>100 mm) to aggregate under shear flow conditions. Therefore 1/W3, the probability to form a stable aggregate when collision occurs by the third mechanism, was multiplied by a function decreasing linearly from 1 to 0 when di ¼ 2 ai increases from 50 mm to 500 mm. The fixed pivot technique was used to solve the PBE based on a discretization approach with a geometric grid defined with a progression factor equal to 2. 3.1.3. Population balance model parameters Several parameters are necessary for the resolution of the model. Some of them were estimated thanks to literature data (e.g. the radius of the monomer a0 ¼ 2.85 nm and the fractal dimension of aggregates Df ¼ 2.7). The other parameters (stability factors: W1, W2 and W3; unfolding parameters: Ttrans and ca) were fitted by minimization of the difference between laboratory-scale experimental data and predicted data obtained from the model (Erabit et al., 2015). It was found a strong dependence of the stability factors W1 and W3 to the concentration of CaCl2 (between 5.1 and 7.1 mM). The identified parameters as a function of the CaCl2 concentrations are shown in Table 1. 3.2. Coupled PBE and plug flow modelling approach The combined model considers the unfolding and the aggregation of the proteins which depend on temperature and shear rate along the different sections assuming plug flow inside the tubular heat exchanger. Hereafter the model is described for the heating section. The same approach is used for the holding tube and the cooling section. 3.2.1. Plug flow model The plug flow model assumes that all parcels of fluid flowing through the inlet section of a pipe at time t ¼ 0 have the same total residence time tr expressed as follows:
tr ¼ t ¼
V V_
(5)
where V is the volume of the pipe, V_ is the volumetric flow rate and t is the mean residence time. 3.2.1.1. Temperature estimation. The temperature evolution is estimated from experimental measurements. The inlet and outlet temperature of the product: Ti and To respectively, as well as the saturation temperature of the steam Ts, are known data. If the connecting lengths are neglected and it is assumed that the heat
Table 1 Parameters identified for modelling of unfolding and aggregation of b-lactoglobulin. Unfolding
Population balance Model
[CaCl2]
Ttrans
ca
W1
5.0 6.0 6.4 6.6 7.0
363.9 K
9064 K
3.17 2.06 1.86 1.81 1.83
mM mM mM mM mM
109 109 109 109 109
W2
W3
2512
413.6 384.8 173 41 34.6
35
transfer coefficient is the same along the heat exchanger length, then (Ts-T)/(Ts-Ti) decreases exponentially versus the position z along the heat exchanger or versus the residence time t of the fluid parcel at this position.
Ts T t ¼ exp Ts Ti tth
with
tth ¼
t r Ts Ti ln Ts To
(6)
The heater is studied in steady state condition; therefore integrating along the heat exchanger length lh (z varying from 0 to lh) is equivalent to an integration versus time (t varying from 0 to tr) by following one fluid parcel. One can imagine that such a fluid parcel occupies all the inner section of the tubular exchanger. For the holding tube and the cooling section the temperature evolution is assumed linear:
T Ti t ¼ To Ti tr
(7)
3.2.1.2. Shear rate estimation. The shear rate varies in function of the radial position in the tube: 0 on the axis, maximal near the wall. The maximal shear rate is determined from the volumetric flow rate V_ and the inner radius of the tube R. The average bulk shear rate value was approximated using the expression for a parabolic velocity profile as shown in Table 2). The heating/cooling and holding sections have different tube diameters resulting in different maximal and average shear rates (Table 2).
3.2.2. Model inputs and outputs The model is able to estimate for fluid parcels having a residence time tr: - the outlet concentration of native proteins: Co(tr) - the outlet volume fraction of aggregates: fo(tr) - the outlet aggregate size distribution: No(i, tr) representing the number of aggregates per volume unit in each pivot. The necessary input data are: - the residence time tr _ - the average shear rate g. - the inlet, outlet and steam temperatures: respectively Ti, To and Ts - the inlet concentration of native proteins: Ci - the inlet aggregate size distribution: Ni(i) It is assumed that there are no aggregates at the inlet of the heater; all the proteins are in the native form. For the holding tube, the input data about native proteins and aggregate size are the outputs from the heating section. In the same way, output data of the holding section represent inputs for the cooling section.
3.3. Coupled PBE and empirical RTD modelling approach The approach coupling PBE to actual RTD of the heat exchangers considers that the fluid parcels flowing through the heat treatment system have different residence times. As mentioned above, the RTD was experimentally measured in the heat and holding sections (knowing that the heat and the cooling sections have the same geometry, and thus similar RTD) and an empirical model was developed (Ndoye et al., 2012).
36
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
Table 2 Equations and values for maximal and average shear rate in the different section of the pilot. Step
Maximal shear rate (s1)
Diameter (mm)
_
_
Heating and cooling sections Holding section
4V g_ max ¼ pR 3
32V g_ ¼ 15pR 3
8.16
104.1
55.5
7.08
159.5
85.0
3.3.1. Empirical RTD model Experimental results obtained were accurately fitted with a new compartment model based on the generalized convection model. The compartment model was chosen in order to reproduce the physical structure of the heat exchanger (straights tubes linked by bends) and of the holding tube (helical). The developed model also takes into account the mixing effect in the bends and the curvature of the helical holding tube, due to the apparition of Dean vortices. The model assumed convective laminar flow in n tubes associated in series and separated by perfect mixing zone. The generalized convection model is used to model the RTD of a single subsystem (straight tube) and a convolution product of the local RTD was used in order to obtain the global RTD after m tubes as defined in Eq. (8)
Zq Em ðqÞ ¼
Average shear rate (s1)
Temperature mapped in Fig. 1 and RTD profiles from Ndoye et al. (2012) respectively were used to calculate the temperature profiles of the different fluid parcels. For a fluid parcel having a residence tr the unfolding model described in section 3.1.1 is applied to calculate the native protein concentration Co(tr) at the outlet of the studied device. The number of aggregates per volume unit in each pivot No(i, tr) at the outlet was obtained from the PBE as described in section 3.1.2. The fraction of fluid parcels having a residence time comprised between tr dtr/2 and tr þ dtr/2 is E(tr)dtr. Therefore the bulk outlet values become:
Z∞ Co ¼
Co ðtr ÞEðtr Þdtr
(10)
0
E1 ðq0 ÞEm1 ðq q0 Þdq0
Z∞
(8) No ðiÞ ¼
q0 n
No ði; tr ÞEðtr Þdtr
(11)
0
where E1 is the pulse response after one of the n tubes associated in series (Eq. (9)):
E1 ðqÞ ¼
1 1 1 q0 1q0 1 q0 nq nq
for
q¼
t q0 t n
4. Results and discussion
(9) 4.1. Comparison between the simulation results of the two models
q0 is the breakthrough time and t the geometrical residence time. For the heating and cooling sections, n ¼ 8 and the breakthrough time q0 were identified using experimental results. For holding section, both n and q0 were identified by minimisation between experimental and predicted data. Table 3 presents the results obtained by Ndoye et al. (2012) and used in this work. These results correspond to a mean residence time t ¼ 43.5 s in the heater and the cooler and t ¼ 49.2 s in the holding tube.
3.3.2. Coupled model of PBE and RTD Hereafter the model is described for the heater. The same approach is used for the holding tube and the cooler. It is assumed that at the outlet of the heater, all the fluid parcels are mixed together before entering the holding tube (and so on for the cooler). Therefore, for the holding tube, the input information about native proteins and aggregate size distribution are the same for all the fluid parcels and correspond to the bulk values at the outlet of the heater. It is expected that fluid parcels which remains longer in the heater (i.e. which have higher residence time) will exit with a lower fraction of native proteins and with larger aggregates.
Table 3 Data of the RTD empirical model: values of the number of tubes n and the breakthrough time q0 used for the different section of the pilot (obtained from Ndoye et al. (2012)). Step
q0
n
Heating and cooling sections Holding section
0.703 0.690
8 3
Results of simulation obtained with the approach combining PBE to plug flow model are compared to those obtained when using the model combining PBE to empirical RTD. 4.1.1. Evolution of the characteristics of the suspensions during heat treatment Fig. 2 maps the evolution of the characteristics of the suspensions (residual native fraction, volume weighted average diameter and volume fraction of aggregates) predicted by the combined PBE/ plug flow model as a function of the residence time during a 87 C and 20 L h1 pilot heat treatment for a solution with 6 mM of CaCl2. It can be seen in Fig. 2a that during the heating phase, the residual native fraction of b-lg decrease slowly from 1 to 0.8. The decrease continued more quickly during the holding step and the protein residual native fraction reached a value of about 0.5. This value remained relatively unchanged during the cooling phase. In Fig. 2b, the model predicts a slight increase of the volume weighted average diameter inside the heating section while it increased rapidly during the holding phase, correspondingly to the rapid decrease of the protein residual native fraction shown in Fig. 2a. This rise of the average diameter of aggregates continued in the cooling heat exchanger. The evolution of the aggregates volume fraction represented in Fig. 2c shows a very slight variation in the heater, followed by an obvious increase in the holder and a slowdown in the cooler. These results illustrate the significant impact of holding on the aggregation mechanism especially on the aggre~ iga et al., 2010; Erabit et al., 2014). gates size (Zún Predicted evolutions of residual native fraction, volume weighted average diameter and volume fraction of aggregates using the modelling approach combining the PBE and the empirical
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
Fig. 2. Simulated evolution, with plug flow model, of the characteristics of the suspensions during heat treatment (20 l h1, To ¼ 87 C, 6% b-lg þ 6 mM of CaCl2) Residual native fraction (a), volume weighted average diameter (b) and volume fraction of aggregates (c).
RTD are represented in Fig. 3. Simulations were performed for a b-lg solution with 6 mM of CaCl2 heated at 87 C and flowing at 20 L h1 through the pilot heat treatment unit. The model simulates aggregation in the whole heat treatment device for different fractions of fluids with various residence times (q ¼ 0.8, 1, 1.2, 1.4 and 1.6). For example q ¼ 1.6 corresponds to the fraction of fluid which spends 60% longer in the section than the fraction corresponding to the mean residence time. For each section, the suspensions
37
Fig. 3. Simulated evolution, with combined model, of the characteristics of the suspensions during heat treatment (20 l h1, To ¼ 87 C, 6% b-lg þ 6 mM of CaCl2) for different residence times - Residual native fraction (a), volume weighted average diameter (b) and volume fraction of aggregates (c).
characteristics considered are weighted using the residence time distribution model to obtain bulk values. Fig. 3a shows that the increase in residence time results in a decrease in residual native fraction of protein, and that this effect is more pronounced in the holding section as previously observed with the plug-flow model. In Fig. 3b and c, it can be seen that the slower is the fluid parcel, the higher is the aggregate size and the volume fraction respectively. Bulk values are close to that found at q ¼ 1, corresponding to the geometrical residence time. The discrepancies between results of
38
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
the different fractions of fluid are more noticeable during the holding. This can be explained by the longer residence time inside the holding tube, but especially by the important modification of the velocity profile inside the holding tube when aggregation occurs as reported by Ndoye et al. (2012). Indeed, RTD is more dispersed in the holding tube compared to heat exchangers, resulting to more differences between the fractions of fluid. 4.1.2. Evolution of the size distributions ASD at the outlet of each of the three sections of the heat treatment pilot were simulated for different fractions of fluid using the PBE/RTD model and using the PBE/plug flow model. The predicted size distributions for a b-lg solution with 6 mM of CaCl2 heated at 87 C and flowing at 20 L h1are summarized in Fig. 4. The
Fig. 4. Simulated volume size distribution during heat treatment (20 l h1, To ¼ 87 C, 6% b-lg þ 6 mM of CaCl2) for the plug flow model and for different residence times and bulk distribution of the PBE/RTD combined model e After heating section (a), holding section (b) and cooling section (c).
graphs present the volume occupied by each size class of aggregates per volume unit of suspension. For the PBE/empirical RTD model, five profiles (from q ¼ 0.8 to q ¼ 1.6) as well as the bulk profile are represented. It can be seen that the predicted distribution changed as a function of the residence time. For each of the three steps, the volume and size of aggregates are as large as the fluid parcel is slow. This result can be explained by an increase of the protein aggregation resulting in the formation of more and larger aggregates as shown on Fig. 4. After holding, the aggregates volume size increased significantly compared to the heating section. The distributions are also broader at the exit of the holding tube than for the heating section. ASD obtained at the end of the cooling section show slightly larger aggregates especially for the larger residence time. Unlike ASD for the heating and the holding sections differences between the distributions corresponding to the different residence time are less important for the cooling section. This is principally due to a lower mean residence time and a lower average temperature leading to less transformation in this section. Comparing bulk and plug flow results, it can be seen significant differences between the distributions obtained after the heating, holding and cooling sections. The graphs on Fig. 4 show that the PBE/RTD combined model (bulk) predicts larger aggregates and broader distributions compared to the model considering plug flow. This result shows that the plug flow approach can not represent satisfactorily neither the fluid flow inside the heat exchangers, nor the induced aggregation phenomena. 4.1.3. Effect of calcium concentration Ionic environment has a significant influence for b-lg aggregation. The important role of calcium on heat induced b-lg aggregation under shear was previously demonstrated (Erabit et al., 2013; Ndoye et al., 2013). Several mechanisms explaining the influence of calcium was also identified: ion-induced conformational change, electrostatic shielding, formation of protein-Ca2þ-protein complexes, screening of the protein surfaces charge (Parris et al., 1993; Simons et al., 2002; Mounsey and O'Kennedy, 2007). In order to understand how calcium impacts the dynamic of formation of aggregates in the model, it is useful to look at the size distributions. For three concentrations of CaCl2 (5.0, 6.0 and 7.0 mM), the simulated size distributions after holding are presented on Fig. 5. Comparison between heat treatment with 5.0 and 6.0 mM, shows that the models predict relatively similar aggregate size distributions but the concentrations is much higher with 6.0 mM. For the heat treatment with 7.0 mM of CaCl2, the mean size is much larger and there are large differences between the bulk distribution obtained from the PBE/empirical RTD model and the distribution given by the PBE/plug flow model. The latter shows aggregates size ranging from 0.03 mm to 20 mm with a median around 6 mm. The distribution of the PBE/empirical RTD model is very wide, ranging from 0.03 mm to about 500 mm. An analysis of the values of stability factors W1 and W3 (Table 1) can help to understand the effect of calcium chloride on the aggregation process. W1 expresses primary protein aggregation trough disulphide bonds to form small aggregates and W3 represents the stability ratio for large insoluble aggregates. It can be seen on Table 1 that the formation of large aggregates (W3) is more impacted by the calcium concentration than the formation of small aggregates (W1). Indeed, for calcium concentrations varying from 5.0 mM to 7.0 mM, the value of W3 lowers from 413.6 to 34.6 representing a decrease of 92% while W1 decreases of only 42% (from 3.17 109 to 1.83 109). This result means that it is more probable to form large aggregates at 7 mM of calcium concentration than at 5 mM, as evidenced on Fig. 5. This outcome may be explained by the fact that at low concentrations, charge shielding is the predominant effect of ionic calcium on protein interactions. Aggregation is promoted either by increasing
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
39
Fig. 6. Simulation with plug flow model and PBE/RTD combined model of the volume fraction of aggregates obtained by pilot heat treatment (20 l h1, To ¼ 87 C, 6% b-lg) with different calcium chloride concentrations.
fraction of native proteins is lowered while the average size and volume fraction of aggregates increase, meaning that calcium chloride promotes b-lg aggregation. The slope of the curve changes dramatically between 6.4 and 6.6 mM of calcium chloride. An objective of the present study is to compare the predictions of the two models in order to determine whether RTD has a significant effect on the selected characteristics of the aggregates suspension. Regarding the aggregates volume fraction, Fig. 6 shows an appreciable difference between the two models for the higher calcium chloride concentrations (6.6 and 7.0 mM), meaning that the effect of RTD on aggregation depends on the calcium concentration. This ending may be explained by the fact that few transformations occurred at the low concentrations of CaCl2 (5.0, 6.0 and 6.4, mM) as evidenced by the small size of aggregates (<1 mm) shown in Fig. 7b. Differences between plug flow and empirical RTD become important when the product is subjected to important transformations (i.e. with 6.6 and 7.0 mM of CaCl2). In these cases, considering the residence time distribution may be crucial in order to have a good representation of the aggregation phenomena. 4.2. Comparison of model prediction with experimental results
Fig. 5. Simulated volume size distributions for the plug flow model and the bulk distribution of the PBE/RTD combined model after holding (20 l h1, To ¼ 87 C, 6% blg) for different CaCl2 concentrations: 5.0 mM (a), 6.0 mM (b) and 7.0 mM (c).
electrostatic attraction between negatively charged proteins, or by decreasing electrostatic repulsion between positively charged proteins, forming essentially small aggregates. But, at high concentrations, in addition to charge shielding effects, large aggregates are formed by intermolecular ion binding to form protein-Caprotein complexes (Parris et al., 1993; Chi et al., 2003; Westerik et al., 2015),. Table 1 also shows a significant gap between values of W3 at 6.4 mM and 6.6 mM suggesting a threshold effect of calcium, probably due to the existence of a stoichiometric calcium/b-lg ratio (Mounsey and O'Kennedy, 2007). This threshold effect is evidenced in Fig. 6 which presents the aggregates volume fraction as function of the concentration of CaCl2 for both models. It can be seen that the calcium concentration affects the predicted properties of the suspensions. When calcium concentration increases, the
By combining the PBE model to the empirical RTD model, one expects a better prediction of the aggregate properties compared to the coupling with plug flow model. This section describes the comparison of experimental and simulated data on the basis of properties related to the functionalities of the suspensions: residual native fraction, viscosity of the suspension and aggregate size distributions measured with laser granulometer (range 0.1e1000 mm). 4.2.1. Functional properties of the aggregate suspensions The predicted (by both models) and experimental results of two aggregates functionality indicators (residual native fraction and median volume diameter) obtained for a heat treatment of 87 C/ 49 s and for different CaCl2 concentration are compared in Fig. 7. Both models overestimate the residual native fraction compared to experimental data (Fig. 7a) even if the general trend regarding the influence of CaCl2 is well represented. This discrepancy between experimental and modelling results may be explained by the time-temperature evolution of the different fractions of fluids inside the heat exchanger. Indeed the fraction of fluids flowing near the heat exchanger wall is submitted to a higher temperature and to a longer residence time than the others, resulting in a higher level of denaturation. On the other hand, fluid parcels moving on the tube axis undergoes low unfolding as a result of relatively lower
40
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
model and well predicted by the PBE/empirical RTD model. This result put in evidence, once again, the importance to account for the histories of the different fractions of fluids when aggregation is high (high level of calcium concentration) in order to represent satisfactorily the aggregation process. 4.2.2. Volume size distributions Using Population Balance Model make it possible to predict aggregate size distributions which are the key feature of protein aggregates functionalities. Comparisons between predicted and experimental cumulative size distributions for three concentrations of CaCl2 (5.0 mM, 6.0 mM and 7.0 mM) are presented on Fig. 8. For 5.0 and 6.0 mM of CaCl2 shown on Fig. 8a and b respectively, both models predict distributions in good agreement with experimental results and consistent with average volume size presented
Fig. 7. Experimental and simulated (plug flow model and PBE/RTD combined model) aggregate suspension properties after pilot heat treatment (20 l h1, To ¼ 87 C for 49 s) - Residual native fraction (a) and Median volume diameter (b).
temperature and shorter residence time than the other fluid fractions. Fluid trajectories between the axis and the wall have intermediates unfolding levels. Actually, samples recovered at the outlet of the cooling section constitute a blending of these different fractions of fluid, ensuing a compensation effect between the different unfolding levels. It is not surprising that the plug flow model fails to represent this effect because it considers that all fractions of fluids have the same residence time and timetemperature evolution. Contrariwise, the bulk data derived from the PBE/RTD model considers a blend of different fluid parcels. These fluid parcels are assumed to have different residence time but follow the same temperature profile (vs position). The fact that fluid parcels experience different temperature evolutions depending on distance from the heating wall is not taken into account on this approach. This could explain the observed discrepancy. These findings join those of Chantoiseau et al. (2012) who demonstrated that dynamical and thermal histories associated with fluid parcels which progress more or less quickly, far or close to the heating wall determine the level of aggregation of WP. In the region near the heat exchanger wall, the product is heated for a long time, resulting in important aggregation; by the contrary, in the axial region where the residence time is shorter, the product is less heated and weakly transformed. The authors pointed out the importance to take into account every fluid parcels trajectories instead of considering a plug flow model, in order to accurately predict the b-lg aggregation inside the heat exchanger section. Fig. 7b shows good agreement between experiments and both modelling results of median volume diameter for the three smallest concentrations of CaCl2. For the two higher concentrations, the median volume diameter is underestimated by the PBE/plug flow
Fig. 8. Experimental and simulated (plug flow model and PBE/RTD combined model) aggregate size distribution after pilot heat treatment (20 l h1, To ¼ 87 C for 49 s) with different calcium chloride concentrations: 5 mM (a), 6 mM (b) and 7 mM (c).
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
on Fig. 7b. For a calcium concentration of 7.0 mM, the plug flow model largely underestimates the aggregates size. As can be seen on Fig. 8c and 90% of aggregates predicted by the plug flow model have a size below 50 mm while this proportion is only of 30% for the experimental results. This could be attributed to the fact that formation of large aggregates in the slowest fluid parcels is not taken into account. The second model which considers RTD predicts higher concentrations of large aggregates compared to plug flow model (only about 48% of aggregates have a size below 50 mm). Fig. 8c shows that the PBE/RTD combined model predicts aggregates volume size distribution close to the experimental result. The median volume size is in good agreement with the experimental values. However, one can evidenced the fact that the experimental distribution shape is not exactly predicted. The slight dissimilarity between experimental and simulated distributions may be explained by the fact that the model does not account some other physico-chemical factors, such as pH, which can have a significant impact on aggregation process, especially on the formation of large ~ iga et al., 2010). aggregates (Zún Finally, the PBE/RTD combined model is more realistic than the plug flow model and gives better prediction for aggregates size distribution. 5. Conclusion A model for the prediction of the b-lg aggregation during the heat treatment in a continuous pilot device has been presented. It couples a RTD model and a colloidal aggregation model previously developed. This model has been compared to a plug flow model in which only the average time/temperature history is considered. The integration of RTD appears to be very important when the product is subject to much transformation (formation of large aggregates > 10 mm). The predictions were compared to the experimental data. The characteristics of the suspensions for experiments carried with the lowest CaCl2 concentrations (5.0, 6.0 and 6.4 mM) are well predicted. For higher concentrations of calcium (6.6 and 7.0 mM), the PBE/RTD combined model well predicts the average aggregates size. Further developments should interestingly integrate the influence of physico-chemical variables other than calcium concentration in order to better describe the distributions of aggregates size. The modelling approach combining PBE to RTD developed in this study is valuable for enhancing the performance of WP aggregation process. Such a hybrid model is essential in process design, evaluation and analysis. Acknowledgements This study was carried out within the framework of the Globule project (ANR-08-ALIA-08) with the help of the other partners: INRA - PIHM and GEA Process Engineering France. References Anema, S.G., McKenna, A.B., 1996. Reaction kinetics of thermal denaturation of whey proteins in heated reconstituted whole milk. J. Agric. Food Chem. 44 (2), 422e428. Arellano, M., Benkhelifa, H., Alvarez, G., Flick, D., 2013. Coupling population balance and residence time distribution for the ice crystallization modeling in a scraped surface heat exchanger. Chem. Eng. Sci. 102 (0), 502e513. Arosio, P., Beeg, M., Nicoud, L., Morbidelli, M., 2012. Time evolution of amyloid fibril length distribution described by a population balance model. Chem. Eng. Sci. 78 (0), 21e32. Bas¸, N., Catak, M., Zumaeta, N., Fitzpatrick, J.J., Cronin, K., Byrne, E.P., 2009. Population balance modelling of protein precipitates exhibiting turbulent flow through a circular pipe. Chem. Eng. Sci. 64 (18), 4051e4059. Bilgili, E., Scarlett, B., 2005. Numerical simulation of open-circuit continuous mills using a non-linear population balance framework: incorporation of non-first-
41
order effects. Chem. Eng. Technol. 28 (2), 153e159. Chantoiseau, E., Plana-Fattori, A., Doursat, C., Flick, D., 2012. Coupling fluid flow, heat transfer and thermal denaturation-aggregation of beta-lactoglobulin using an Eulerian/Lagrangian approach. J. Food Eng. 113 (2), 234e244. Chi, E., Krishnan, S., Randolph, T., Carpenter, J., 2003. Physical stability of proteins in aqueous solution: mechanism and driving forces in nonnative protein aggregation. Pharm. Res. 20 (9), 1325e1336. Dannenberg, F., Kessler, H.G., 1988. Reaction kinetics of the denaturation of whey proteins in milk. J. Food Sci. 53 (1), 258e263. Elofsson, U.M., Dejmek, P., Paulsson, M.A., 1996. Heat-induced aggregation of betalactoglobulin studied by dynamic light scattering. Int. Dairy J. 6 (4), 343e357. Erabit, N., Flick, D., Alvarez, G., 2013. Effect of calcium chloride and moderate shear on b-lactoglobulin aggregation in processing-like conditions. J. Food Eng. 115 (1), 63e72. Erabit, N., Flick, D., Alvarez, G., 2014. Formation of b-lactoglobulin aggregates during thermomechanical treatments under controlled shear and temperature conditions. J. Food Eng. 120 (0), 57e68. Erabit, N., Ndoye, F.T., Alvarez, G., Flick, D., 2015. A population balance model integrating some specificities of the b-lactoglobulin thermally-induced aggregation. J. Food Eng. 144 (0), 66e76. Heath, A.R., Koh, P.T.L., 2003. Combined population balance and CFD modelling of particle aggregation by polymeric flocculant. In: Third International Conference on CFD in the Minerals and Proocess Industries, Melbourne, Australia, pp. 339e344. Lee, K.F., Mosbach, S., Kraft, M., Wagner, W., 2015. A multi-compartment population balance model for high shear granulation. Comput. Chem. Eng. 75 (0), 1e13. Lian, G., Moore, S., Heeney, L., 2006. Population balance and computational fluid dynamics modelling of ice crystallisation in a scraped surface freezer. Chem. Eng. Sci. 61 (23), 7819e7826. Mounsey, J.S., O'Kennedy, B.T., 2007. Conditions limiting the influence of thioledisulphide interchange reactions on the heat-induced aggregation kinetics of b-lactoglobulin. Int. Dairy J. 17 (9), 1034e1042. Ndoye, F.T., Erabit, N., Alvarez, G., Flick, D., 2012. Influence of whey protein aggregation on the residence time distribution in a tubular heat exchanger and a helical holding tube during heat treatment process. J. Food Eng. 112 (3), 158e167. Ndoye, F.T., Erabit, N., Flick, D., Alvarez, G., 2013. In-line characterization of a whey protein aggregation process: aggregates size and rheological measurements. J. Food Eng. 115 (1), 73e82. Nicolai, T., Britten, M., Schmitt, C., 2011. b-Lactoglobulin and WPI aggregates: formation structure and applications. Food Hydrocoll. 25, 1945e1962. Parris, N., Anema, S.G., Singh, H., Creamer, L.K., 1993. Aggregation of whey proteins in heated sweet whey. Journal of Agricultural and Food Chemistry J. Agric. Food Chem. 41, 460e464. Prasetya, A., Liu, L., Litster, J., Watanabe, F., Mitsutani, K., Ko, G.H., 1999. Dynamic model development for residence time distribution control in high-impact polypropylene copolymer process. Chem. Eng. Sci. 54 (15e16), 3263e3271. Sen, M., Ramachandran, R., 2013. A multi-dimensional population balance model approach to continuous powder mixing processes. Adv. Powder Technol. 24 (1), 51e59. Sen, M., Singh, R., Vanarase, A., John, J., Ramachandran, R., 2012. Multi-dimensional population balance modeling and experimental validation of continuous powder mixing processes. Chem. Eng. Sci. 80 (0), 349e360. Simmons, M.J.H., Jayaraman, P., Fryer, P.J., 2007. The effect of temperature and shear rate upon the aggregation of whey protein and its implications for milk fouling. J. Food Eng. 79 (2), 517e528. Simons, J.-W.F.A., Kosters, H.A., Visschers, R.W., de Jongh, H.H.J., 2002. Role of calcium as trigger in thermal b-lactoglobulin aggregation. Arch. Biochem. Biophys. 406 (2), 143e152. Smoluchowski, M., 1917. Versuch einer mathematischen Theorie der Koagula€sungen. Z. Phys. Chem. 92, 129e168. tionskinetic kolloider, Lo Spiegel, T., Huss, M., 2002. Whey protein aggregation under shear conditions - effects of pH-value and removal of calcium. Int. J. Food Sci. Technol. 37, 559e568. Tanguay, M., Fawell, P., Adkins, S., 2014. Modelling the impact of two different flocculants on the performance of a thickener feedwell. Appl. Math. Model. 38 (17e18), 4262e4276. Tolkach, A., Kulozik, U., 2007. Reaction kinetic pathway of reversible and irreversible thermal denaturation of b -lactoglobulin. Lait 87 (4e5), 301e315. Vale, H.M., McKenna, T.F., 2005. Modeling particle size distribution in emulsion polymerization reactors. Prog. Polym. Sci. 30 (10), 1019e1048. Weller, K.R., Spencer, S.J., Gao, M., Liu, Y., 2000. Tracer studies and breakage testing in pilot-scale stirred mills. Miner. Eng. 13 (4), 429e458. Westerik, N., Scholten, E., Corredig, M., 2015. The effect of calcium on the composition and physical properties of whey protein particles prepared using emulsification. Food Chem. 177, 72e80. Wynn, E.J.W., Hounslow, M.J., 1997. Integral population balance equations for growth. Chem. Eng. Sci. 52 (5), 733e746. Zumaeta, N., Byrne, E.P., Fitzpatrick, J.J., 2010. Application of CFD and breakage modelling for predicting the size reduction of protein precipitates during transport. Chem. Eng. Res. Des. 88 (7), 911e917. Zumaeta, N., Cartland-Glover, G.M., Heffernan, S.P., Byrne, E.P., Fitzpatrick, J.J., 2005. Breakage model development and application with CFD for predicting breakage of whey protein precipitate particles. Chem. Eng. Sci. 60 (13), 3443e3452. ~ iga, R.N., Tolkach, A., Kulozik, U., Aguilera, J.M., 2010. Kinetics of formation and Zún physicochemical characterization of thermally-induced b-lactoglobulin aggregates. J. Food Sci. 75 (5), E261eE268.
Contents lists available at ScienceDirect
Journal of Food Engineering journal homepage: www.elsevier.com/locate/jfoodeng
Coupling population balance model and residence time distribution for pilot-scale modelling of b-lactoglobulin aggregation process Nicolas Erabit a, b, c, Fatou Toutie Ndoye a, *, Graciela Alvarez a, Denis Flick b, c a
IRSTEA Refrigeration Process Engineering Research Unit, Antony, France AgroParisTech, UMR1145 Ing enierie Aliments Proc ed es, Massy, France c INRA, UMR1145 Ing enierie Aliments Proc ed es, Massy, France b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 1 September 2015 Received in revised form 2 November 2015 Accepted 20 December 2015 Available online 23 December 2015
This study presents a hybrid modelling approach for the prediction of the thermally-induced aggregation of whey protein solutions in a continuous process. Modelling a continuous process needs to account the fluid flow phenomena because the different fluid parcels are not subjected to the same time/temperature/shear history. The modelling approach developed combines a residence time distribution (RTD) model and a colloidal aggregation model based on Population Balance Equations (PBE). This combined model has been compared, on one side, to a Population Balance Model assuming plug flow, and on the other side, to experimental data. A 6% b-lactoglobulin solution with different CaCl2 concentrations (5.0, 6.0, 6.4, 6.6 and 7.0 mM) was used for experimental investigations. Experimental heat treatments were carried out at 87 C with a holding time of 49 s using tubular heat exchanger and holder. The integration of residence time distribution appears to be very important when the product is subject to large transformations (formation of large aggregates > 10 mm). Comparison of predictions obtained from both models to the experimental data shows that the coupled PBE/RTD model is able to predict the characteristics of the suspensions (aggregate size and viscosity) and successfully integrates the impact of calcium upon aggregation process, while the model assuming plug flow is suitable only when aggregation is low (for low calcium concentrations). © 2016 Elsevier Ltd. All rights reserved.
Keywords: Population balance model Residence time distribution Whey protein ta-Lactoglobulin be Aggregation
1. Introduction A typical dairy process involves multiple continuous operations, among which heat treatment is the most common processing step. Thermal treatment is performed with the aim of stabilizing (microbiologically and physically) or modifying the structure of the components. Considering the specific case of structuring of dairy products, whey protein (WP) aggregation process is still a challenge in term of product quality reproducibility. Thermal and mechanical operations applied to WP induce denaturation and aggregation of protein molecules leading to particles with highly sought functionalities such as foaming, emulsifying, gelling and thickening agents. These properties of WP aggregates are widely governed by the aggregates size distribution (ASD) depending on environmental factors such as temperature, pH and ionic strength. It is necessary to implement consistent process systems engineering approaches in
* Corresponding author. E-mail address: [email protected] (F.T. Ndoye). http://dx.doi.org/10.1016/j.jfoodeng.2015.12.013 0260-8774/© 2016 Elsevier Ltd. All rights reserved.
order to provide final product with the desired target functionalities according to process variables. A model-based system approach can deliver the required knowledge for process understanding, evaluation and analysis with reduced time and cost for process design. Modelling approaches commonly applied to heat induced WP aggregation are based on reaction kinetic methods (Dannenberg and Kessler, 1988; Anema and McKenna, 1996; Tolkach and Kulozik, 2007). Such an approach can be useful to provide global knowledge about the native proteins and aggregates concentrations, but it fails to depict all the physical phenomena involved in the aggregation mechanism. Moreover, the reaction kinetic based model is unable to describe the evolution of aggregates size along the process. A modelling strategy that has the potential to overcome these limitations is the application of population balance model (PBM) to the aggregates formation. There are relatively few published works in which PBM is used to describe ASD of heat induced WP aggregates (Elofsson et al., 1996; Arosio et al., 2012; Chantoiseau et al., 2012; Erabit et al., 2015). The recent work of Erabit et al. (2015) was devoted to the development of a new
32
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
Nomenclature
Latin symbols Radius of monomers (m) a0 ai Radius of aggregates in size class i (m) aj Radius of aggregates in size class j (m) ca unfolding constant (K) C Concentration of native monomers (kg m3) Cmon Concentration of monomers (kg m3) Cunf Concentration of unfolded monomers (kg m3) di Diameter of aggregates in size class i (m) Df Fractal dimension () Dv(0.5) Median volume diameter (m) E(t) Residence time distribution function (s1) E(q) Dimensionless residence time distribution function () h Proportionality factor () I Total number of size classes () k Constant of Boltzmann (m2 kg s2 K1) kij Total aggregation kernel implying size classes i and j (m3 s1) kperi,ij Perikinetic aggregation kernel implying size classes i and j (m3 s1) kortho,ij Orthokinetic aggregation kernel implying size classes i and j (m3 s1) lh Heat exchanger length (m) m, n Number of compartments in the RTD compartment model () N Number concentration of aggregates (m3) Ni Number concentration of aggregates in size class i (m3) Nj Number concentration of aggregates in size class j (m3) Nmon Number concentration of native monomers (m3) N1 Number concentration of unfolded monomers (m3) p Upper size class of aggregates involving hydrophobic interactions qi Number of monomers in aggregates of size class i () t Time (s)
modelling approach accounting physico-chemical and thermomechanical conditions to predict b-lactoglobulin (b-lg) aggregates suspensions properties such as ASD, residual native fraction and viscosity. The modelling considers the time/temperature/shear history and Population Balance Equations (PBE) based on FuchsSmoluchowski theory (Smoluchowski, 1917). The novelty of the modelling method is related to the adaptation of the stability factor values according to the specific mechanisms involved in b-lg aggregation. The model gave good results when applied to b-lg aggregation in a batch process simulator at laboratory scale, but cannot be applied directly to a continuous heat treatment process. The main reason is that phenomena related to fluid dynamics were not accounted in this model. However, it is known that WP aggregation under continuous thermal process in heat exchangers involves three aspects: fluid flow, heat transfer and product transformation (denaturation and aggregation) (Spiegel and Huss, 2002; Simmons et al., 2007; Nicolai et al., 2011). The structural transformations results from the combination of these three aspects. Moreover, these transformations induce significant rheological changes, and consequently, important modifications of the
t tr T Theat Ts Ttrans
V V_ W W1 W2 W3 z
Mean residence time (s) Fluid parcel residence time (s) Temperature (K) Maximal heating temperature (K) Saturation temperature of steam (K) Temperature of transition between unfolding limiting and aggregation limiting steps during the heat induced aggregation of b-lg (K) Volume (m3) Volumetric flow rate (m3 s1) Stability factor () Stability factor for monomers and aggregates of size classes i 3 () Stability factor for aggregates of size classes 4 i p () Stability factor for aggregates of size class pþ1 i I () Axial coordinate
Greek symbols Fraction of monomers in the unfolded state () Aggregates volume fraction () g_ Local shear rate (s1) m Dynamic viscosity of the suspension (Pa s) ms Dynamic viscosity of the continuous phase (Pa s) q Dimensionless time () q0 Breakthrough time () t Geometrical residence time in one tube (s) tth Thermal characteristic time of the heating section (s)
a f
Subscripts co cooling exp experimental ho holding i inlet o outlet th theoretical T total
fluid flow behaviour especially near the heat exchangers walls (Chantoiseau et al., 2012; Ndoye et al., 2012). Additionally, it is recognized that the different fluid parcels flowing through a continuous processing unit are not subjected to the same timetemperature-shear profiles because of the residence time distribution (RTD) (Dannenberg and Kessler, 1988). Development of a realistic model that can represent WP aggregation at continuous pilot or industrial scale needs the implementation of a modelling approach accounting fluid dynamics framework, through RTD modelling for example. A hybrid method combining PBM to fluid flow description can potentially provide the necessary supplement understanding for the WP aggregation continuous process. Several hybrid models coupling PBE to dynamical histories have been developed for particulate continuous process such as crystallization (Wynn and Hounslow, 1997; Lian et al., 2006; Arellano et al., 2013), granulation (Sen et al., 2012; Sen and Ramachandran, 2013; Lee et al., 2015), grinding (Weller et al., 2000; Bilgili and Scarlett, 2005), polymerization (Prasetya et al., 1999; Vale and McKenna, 2005), flocculation (Heath and Koh, 2003; Tanguay et al., 2014) and protein breakage (Zumaeta et al.,
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
2005, 2010; Bas¸ et al., 2009). In authors knowledge, the only study coupling fluid flow characteristics and PBM applied to heat induced WP aggregation is that of Chantoiseau et al. (2012). The authors developed a numerical approach that combines computational fluid dynamics (CFD) for solving the fluid flow and heat transfer coupled problem, with a PBE for evaluating b-lg particle aggregation inside a straight section of a tubular heat exchanger. The present study aims to model b-lg aggregation in an entire heat treatment processing unit including heating, holding and cooling sections, by coupling fluid dynamics framework to PBE. Fluid flow heterogeneities are taken into account through RTD measurements and modelling of a continuous laboratory pilot scale thermal process as performed in a previous work by Ndoye et al. (2012). This RTD model is combined to the PBM developed by Erabit et al. (2015) for the thermally-induced b-lg aggregation based on laboratory scale experiments in a batch process simulator. Such a hybrid modelling approach is compared to a PBM approach assuming plug flow. The predicted results were compared to experimental data characterizing the suspensions (ASD and residual native fraction) obtained during the aggregation process of a blg solution in a continuous heat treatment at the laboratory pilot scale. 2. Experimental approaches 2.1. b-lg heat treatment and characterization 2.1.1. Protein solution A 6% b-lg solution was prepared using industrial powder (92% of b-lg, 5.2% of moisture, 2.5% of ash and 0.3% of fat) and deionized water at 40 C under mixing. The protein concentration used is adequate to form a significant number of aggregates without gelling (Erabit et al., 2014). Different concentrations of Ca2þ were obtained by adding a CaCl2 1 M solution to the protein solution. The prepared solution was placed in a water bath at 40 C for 2 h in order to have complete rehydration of the b-lg powder. The solution was then kept at 4 C and used within a day to prevent any destabilization. Five concentrations of CaCl2 have been tested: 5.0, 6.0, 6.4, 6.6 and 7 mM. The pH of the solution was equal to 6.8 without any significant differences between the different concentrations of CaCl2. 2.1.2. Continuous thermal processing unit An ultra-high temperature (UHT) laboratory pilot scale plant (Armfield FT174X, Ringwood, UK) was used to carry out continuous heat treatments. A detailed description of the experimental setup can be found in Ndoye et al. (2012). The product was fed by a volumetric pump at a flow rate of V_ ¼ 20 L h1 in order to maintain a laminar flow. The applied heat treatment and bulk temperature profile are described on Fig. 1. A preheating at 60 C was first performed using a helical tube in a water bath. The product is then heated at T ¼ 87 C in a tubular heat exchanger made of eight concentric double tubes, using steam at saturation temperature Ts as heating media. After heating inside the heat exchanger, the product passes through a holding section composed of helical tube with adaptable volume in which heat losses lead to a slight decrease of the product temperature to Tho. The last section of the processing unit is the cooling heat exchanger similar to the heating section. The product is cooled down to a temperature Tco ¼ 4 C using cold water. Sampling was done at the outlet of the cooling section for further analysis. 2.1.3. Suspension characterization 2.1.3.1. Residual native fraction. The residual native fraction t is the fraction of the proteins which remain in the native state during
33
heat treatment but also of the proteins which unfold but do not aggregate and come back to the native state during cooling. This fraction was determined by HPLC as described by Tolkach and Kulozik (2007). An acetate buffer at pH 4.6 was used to dilute samples (1/25th) and to promote agglomeration of the denaturated proteins. Agglomerated materials was then removed with syringe filters (Ø 0.2 mm). 2.1.3.2. Granulometry. The volume size distribution of the suspension of aggregates was determined using a Laser granulometer Mastersizer ZS (Malvern Instruments Inc., Worcester, UK) for particles from 0.1 to 1000 mm. The measurement principle includes a dilution of samples before performing analysis. The required level of dilution is directly determined by the granulometer depending on the turbidity of the suspension. 2.2. Residence time distribution measurement The experimental characterization of the RTD were previously reported by Ndoye et al. (2012). The RTD measurements were performed in the UHT laboratory scale plant described above using a 6% b-lg solution with 5 mM CaCl2. Conditions permitting no WP aggregation (T ¼ 60 C) and allowing WP aggregation (T ¼ 87 C) were explored at 20 L h1 and 49 L h1. Experiments were carried out using a colorimetric method with Methylene Blue as tracer. A 10 mL sample of 6% b-lg solution coloured with MB (at concentration of 0.01 g L1 for no aggregation conditions and 0.02 g L1 for aggregation conditions) was pulse injected at the entrance of the heat exchanger or of the holding tube according to studied section. The protein samples were continuously collected at the outlet of the device and analysed by the use of a spectrophotometer (Beckman Coulter DU 730, Fullerton, USA). The age distribution function E(t) was used to characterized the RTD and mean residence time was determined for each conditions explored. 3. Modelling First a model of unfolding-aggregation of b-lg based on PBE and assuming plug flow along the heating, holding and cooling sections is described. Then the PBE model is coupled to an empirical RTD model in order to take into account that the fluid parcels of protein suspension flowing throughout the heater, the holding tube and the cooler do not have the same time-temperature-shear history. 3.1. Population balance modelling approach The complete mechanism of thermally induced aggregation of
b-lg is described as a two steps reaction composed of a reversible unfolding of monomers followed by an irreversible aggregation of unfolded proteins. The modelling framework used in this work combines the unfolding step described through an unfolding ratio a and the aggregation step described through PBE. A detailed description of the model development was recently reported by Erabit et al. (2015). The model was validated on the basis of laboratory batch experiments. This section describes briefly the main points of the modelling approach. 3.1.1. Unfolding mechanism e thermo-dependent equilibrium The unfolding phase of the thermally-induced aggregation of blg was modelled based on a temperature-dependant equilibrium between native and unfolded WP molecules as described by Tolkach and Kulozik (2007). The fraction of unfolded proteins among the protein monomers is defined as the unfolding ratio a ¼ Cunf/Cmon. This unfolding ratio is function of the temperature according to Eq. (1):
34
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
Fig. 1. Heat treatment carried out on the pilot plant: (a) schematic representation of the heat treatment sections and (b) applied bulk temperature profile.
8 1 1 > < ln a ¼ ca Ttrans T > : a¼1
for T Ttrans for T > Ttrans
(1)
where ca is the unfolding coefficient and Ttrans is the transition temperature at which unfolding or aggregation is the limiting step of the overall aggregation process of b-lg.
3.1.2. Aggregation mechanism - population balance model The population balance model was developed based on Smoluchowski PBE which describes irreversible aggregation accounting the particles collision mechanism through kernels: perikinetic aggregation for collision by Brownian motion and orthokinetic aggregation for collision by shear flow. The population balance model is described by the following equations: I1 X dNmon ¼ N1 k1j Nj dt j¼1
(2)
i1 X dNi > 1 1 X ¼ kij;j Nij Nj Ni ki;j Nj 2 j¼1 dt j1
(3)
where N1 ¼ a:Nmon is the number concentration of the unfolded monomers related to the number concentration of the total monomers (native and unfolded) Nmon by the unfolding ratio a. Aggregation kernel kij accounts both perikinetic and orthokinetic collisions, and the stability factor W represents the inverse of the probability that a collision between particles actually leads to an aggregate (Eq. (4)).
kij ¼
kperi;ij þ kortho;ij W
(4)
A novelty of the modelling strategy reported by Erabit et al. (2015) is the adaptation of the stability factor to three different mechanisms specific to b-lg aggregation. This leads to the definition of three stability factors according to the interactions involved: (i) W1 for aggregates between unfolded monomers (i ¼ 1) and small oligomers (I 3) through disulphide bonds; (ii) W2 for collisions involving hydrophobic interactions and larger aggregates (4 < i 8000) ensuing; (iii) W3 for collisions between large aggregates (i > 8000).
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
When aggregates of two different size-ranges collide, the stability factor taken into account is that of the smallest aggregate. It is assumed that the value of the stability factor does not depends on the type of collision involved (orthokinetic or perikinetic). Experimental measurements of aggregates size showed that there are few particles larger than 100 mm. This is certainly due to the difficulty for two large particles (>100 mm) to aggregate under shear flow conditions. Therefore 1/W3, the probability to form a stable aggregate when collision occurs by the third mechanism, was multiplied by a function decreasing linearly from 1 to 0 when di ¼ 2 ai increases from 50 mm to 500 mm. The fixed pivot technique was used to solve the PBE based on a discretization approach with a geometric grid defined with a progression factor equal to 2. 3.1.3. Population balance model parameters Several parameters are necessary for the resolution of the model. Some of them were estimated thanks to literature data (e.g. the radius of the monomer a0 ¼ 2.85 nm and the fractal dimension of aggregates Df ¼ 2.7). The other parameters (stability factors: W1, W2 and W3; unfolding parameters: Ttrans and ca) were fitted by minimization of the difference between laboratory-scale experimental data and predicted data obtained from the model (Erabit et al., 2015). It was found a strong dependence of the stability factors W1 and W3 to the concentration of CaCl2 (between 5.1 and 7.1 mM). The identified parameters as a function of the CaCl2 concentrations are shown in Table 1. 3.2. Coupled PBE and plug flow modelling approach The combined model considers the unfolding and the aggregation of the proteins which depend on temperature and shear rate along the different sections assuming plug flow inside the tubular heat exchanger. Hereafter the model is described for the heating section. The same approach is used for the holding tube and the cooling section. 3.2.1. Plug flow model The plug flow model assumes that all parcels of fluid flowing through the inlet section of a pipe at time t ¼ 0 have the same total residence time tr expressed as follows:
tr ¼ t ¼
V V_
(5)
where V is the volume of the pipe, V_ is the volumetric flow rate and t is the mean residence time. 3.2.1.1. Temperature estimation. The temperature evolution is estimated from experimental measurements. The inlet and outlet temperature of the product: Ti and To respectively, as well as the saturation temperature of the steam Ts, are known data. If the connecting lengths are neglected and it is assumed that the heat
Table 1 Parameters identified for modelling of unfolding and aggregation of b-lactoglobulin. Unfolding
Population balance Model
[CaCl2]
Ttrans
ca
W1
5.0 6.0 6.4 6.6 7.0
363.9 K
9064 K
3.17 2.06 1.86 1.81 1.83
mM mM mM mM mM
109 109 109 109 109
W2
W3
2512
413.6 384.8 173 41 34.6
35
transfer coefficient is the same along the heat exchanger length, then (Ts-T)/(Ts-Ti) decreases exponentially versus the position z along the heat exchanger or versus the residence time t of the fluid parcel at this position.
Ts T t ¼ exp Ts Ti tth
with
tth ¼
t r Ts Ti ln Ts To
(6)
The heater is studied in steady state condition; therefore integrating along the heat exchanger length lh (z varying from 0 to lh) is equivalent to an integration versus time (t varying from 0 to tr) by following one fluid parcel. One can imagine that such a fluid parcel occupies all the inner section of the tubular exchanger. For the holding tube and the cooling section the temperature evolution is assumed linear:
T Ti t ¼ To Ti tr
(7)
3.2.1.2. Shear rate estimation. The shear rate varies in function of the radial position in the tube: 0 on the axis, maximal near the wall. The maximal shear rate is determined from the volumetric flow rate V_ and the inner radius of the tube R. The average bulk shear rate value was approximated using the expression for a parabolic velocity profile as shown in Table 2). The heating/cooling and holding sections have different tube diameters resulting in different maximal and average shear rates (Table 2).
3.2.2. Model inputs and outputs The model is able to estimate for fluid parcels having a residence time tr: - the outlet concentration of native proteins: Co(tr) - the outlet volume fraction of aggregates: fo(tr) - the outlet aggregate size distribution: No(i, tr) representing the number of aggregates per volume unit in each pivot. The necessary input data are: - the residence time tr _ - the average shear rate g. - the inlet, outlet and steam temperatures: respectively Ti, To and Ts - the inlet concentration of native proteins: Ci - the inlet aggregate size distribution: Ni(i) It is assumed that there are no aggregates at the inlet of the heater; all the proteins are in the native form. For the holding tube, the input data about native proteins and aggregate size are the outputs from the heating section. In the same way, output data of the holding section represent inputs for the cooling section.
3.3. Coupled PBE and empirical RTD modelling approach The approach coupling PBE to actual RTD of the heat exchangers considers that the fluid parcels flowing through the heat treatment system have different residence times. As mentioned above, the RTD was experimentally measured in the heat and holding sections (knowing that the heat and the cooling sections have the same geometry, and thus similar RTD) and an empirical model was developed (Ndoye et al., 2012).
36
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
Table 2 Equations and values for maximal and average shear rate in the different section of the pilot. Step
Maximal shear rate (s1)
Diameter (mm)
_
_
Heating and cooling sections Holding section
4V g_ max ¼ pR 3
32V g_ ¼ 15pR 3
8.16
104.1
55.5
7.08
159.5
85.0
3.3.1. Empirical RTD model Experimental results obtained were accurately fitted with a new compartment model based on the generalized convection model. The compartment model was chosen in order to reproduce the physical structure of the heat exchanger (straights tubes linked by bends) and of the holding tube (helical). The developed model also takes into account the mixing effect in the bends and the curvature of the helical holding tube, due to the apparition of Dean vortices. The model assumed convective laminar flow in n tubes associated in series and separated by perfect mixing zone. The generalized convection model is used to model the RTD of a single subsystem (straight tube) and a convolution product of the local RTD was used in order to obtain the global RTD after m tubes as defined in Eq. (8)
Zq Em ðqÞ ¼
Average shear rate (s1)
Temperature mapped in Fig. 1 and RTD profiles from Ndoye et al. (2012) respectively were used to calculate the temperature profiles of the different fluid parcels. For a fluid parcel having a residence tr the unfolding model described in section 3.1.1 is applied to calculate the native protein concentration Co(tr) at the outlet of the studied device. The number of aggregates per volume unit in each pivot No(i, tr) at the outlet was obtained from the PBE as described in section 3.1.2. The fraction of fluid parcels having a residence time comprised between tr dtr/2 and tr þ dtr/2 is E(tr)dtr. Therefore the bulk outlet values become:
Z∞ Co ¼
Co ðtr ÞEðtr Þdtr
(10)
0
E1 ðq0 ÞEm1 ðq q0 Þdq0
Z∞
(8) No ðiÞ ¼
q0 n
No ði; tr ÞEðtr Þdtr
(11)
0
where E1 is the pulse response after one of the n tubes associated in series (Eq. (9)):
E1 ðqÞ ¼
1 1 1 q0 1q0 1 q0 nq nq
for
q¼
t q0 t n
4. Results and discussion
(9) 4.1. Comparison between the simulation results of the two models
q0 is the breakthrough time and t the geometrical residence time. For the heating and cooling sections, n ¼ 8 and the breakthrough time q0 were identified using experimental results. For holding section, both n and q0 were identified by minimisation between experimental and predicted data. Table 3 presents the results obtained by Ndoye et al. (2012) and used in this work. These results correspond to a mean residence time t ¼ 43.5 s in the heater and the cooler and t ¼ 49.2 s in the holding tube.
3.3.2. Coupled model of PBE and RTD Hereafter the model is described for the heater. The same approach is used for the holding tube and the cooler. It is assumed that at the outlet of the heater, all the fluid parcels are mixed together before entering the holding tube (and so on for the cooler). Therefore, for the holding tube, the input information about native proteins and aggregate size distribution are the same for all the fluid parcels and correspond to the bulk values at the outlet of the heater. It is expected that fluid parcels which remains longer in the heater (i.e. which have higher residence time) will exit with a lower fraction of native proteins and with larger aggregates.
Table 3 Data of the RTD empirical model: values of the number of tubes n and the breakthrough time q0 used for the different section of the pilot (obtained from Ndoye et al. (2012)). Step
q0
n
Heating and cooling sections Holding section
0.703 0.690
8 3
Results of simulation obtained with the approach combining PBE to plug flow model are compared to those obtained when using the model combining PBE to empirical RTD. 4.1.1. Evolution of the characteristics of the suspensions during heat treatment Fig. 2 maps the evolution of the characteristics of the suspensions (residual native fraction, volume weighted average diameter and volume fraction of aggregates) predicted by the combined PBE/ plug flow model as a function of the residence time during a 87 C and 20 L h1 pilot heat treatment for a solution with 6 mM of CaCl2. It can be seen in Fig. 2a that during the heating phase, the residual native fraction of b-lg decrease slowly from 1 to 0.8. The decrease continued more quickly during the holding step and the protein residual native fraction reached a value of about 0.5. This value remained relatively unchanged during the cooling phase. In Fig. 2b, the model predicts a slight increase of the volume weighted average diameter inside the heating section while it increased rapidly during the holding phase, correspondingly to the rapid decrease of the protein residual native fraction shown in Fig. 2a. This rise of the average diameter of aggregates continued in the cooling heat exchanger. The evolution of the aggregates volume fraction represented in Fig. 2c shows a very slight variation in the heater, followed by an obvious increase in the holder and a slowdown in the cooler. These results illustrate the significant impact of holding on the aggregation mechanism especially on the aggre~ iga et al., 2010; Erabit et al., 2014). gates size (Zún Predicted evolutions of residual native fraction, volume weighted average diameter and volume fraction of aggregates using the modelling approach combining the PBE and the empirical
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
Fig. 2. Simulated evolution, with plug flow model, of the characteristics of the suspensions during heat treatment (20 l h1, To ¼ 87 C, 6% b-lg þ 6 mM of CaCl2) Residual native fraction (a), volume weighted average diameter (b) and volume fraction of aggregates (c).
RTD are represented in Fig. 3. Simulations were performed for a b-lg solution with 6 mM of CaCl2 heated at 87 C and flowing at 20 L h1 through the pilot heat treatment unit. The model simulates aggregation in the whole heat treatment device for different fractions of fluids with various residence times (q ¼ 0.8, 1, 1.2, 1.4 and 1.6). For example q ¼ 1.6 corresponds to the fraction of fluid which spends 60% longer in the section than the fraction corresponding to the mean residence time. For each section, the suspensions
37
Fig. 3. Simulated evolution, with combined model, of the characteristics of the suspensions during heat treatment (20 l h1, To ¼ 87 C, 6% b-lg þ 6 mM of CaCl2) for different residence times - Residual native fraction (a), volume weighted average diameter (b) and volume fraction of aggregates (c).
characteristics considered are weighted using the residence time distribution model to obtain bulk values. Fig. 3a shows that the increase in residence time results in a decrease in residual native fraction of protein, and that this effect is more pronounced in the holding section as previously observed with the plug-flow model. In Fig. 3b and c, it can be seen that the slower is the fluid parcel, the higher is the aggregate size and the volume fraction respectively. Bulk values are close to that found at q ¼ 1, corresponding to the geometrical residence time. The discrepancies between results of
38
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
the different fractions of fluid are more noticeable during the holding. This can be explained by the longer residence time inside the holding tube, but especially by the important modification of the velocity profile inside the holding tube when aggregation occurs as reported by Ndoye et al. (2012). Indeed, RTD is more dispersed in the holding tube compared to heat exchangers, resulting to more differences between the fractions of fluid. 4.1.2. Evolution of the size distributions ASD at the outlet of each of the three sections of the heat treatment pilot were simulated for different fractions of fluid using the PBE/RTD model and using the PBE/plug flow model. The predicted size distributions for a b-lg solution with 6 mM of CaCl2 heated at 87 C and flowing at 20 L h1are summarized in Fig. 4. The
Fig. 4. Simulated volume size distribution during heat treatment (20 l h1, To ¼ 87 C, 6% b-lg þ 6 mM of CaCl2) for the plug flow model and for different residence times and bulk distribution of the PBE/RTD combined model e After heating section (a), holding section (b) and cooling section (c).
graphs present the volume occupied by each size class of aggregates per volume unit of suspension. For the PBE/empirical RTD model, five profiles (from q ¼ 0.8 to q ¼ 1.6) as well as the bulk profile are represented. It can be seen that the predicted distribution changed as a function of the residence time. For each of the three steps, the volume and size of aggregates are as large as the fluid parcel is slow. This result can be explained by an increase of the protein aggregation resulting in the formation of more and larger aggregates as shown on Fig. 4. After holding, the aggregates volume size increased significantly compared to the heating section. The distributions are also broader at the exit of the holding tube than for the heating section. ASD obtained at the end of the cooling section show slightly larger aggregates especially for the larger residence time. Unlike ASD for the heating and the holding sections differences between the distributions corresponding to the different residence time are less important for the cooling section. This is principally due to a lower mean residence time and a lower average temperature leading to less transformation in this section. Comparing bulk and plug flow results, it can be seen significant differences between the distributions obtained after the heating, holding and cooling sections. The graphs on Fig. 4 show that the PBE/RTD combined model (bulk) predicts larger aggregates and broader distributions compared to the model considering plug flow. This result shows that the plug flow approach can not represent satisfactorily neither the fluid flow inside the heat exchangers, nor the induced aggregation phenomena. 4.1.3. Effect of calcium concentration Ionic environment has a significant influence for b-lg aggregation. The important role of calcium on heat induced b-lg aggregation under shear was previously demonstrated (Erabit et al., 2013; Ndoye et al., 2013). Several mechanisms explaining the influence of calcium was also identified: ion-induced conformational change, electrostatic shielding, formation of protein-Ca2þ-protein complexes, screening of the protein surfaces charge (Parris et al., 1993; Simons et al., 2002; Mounsey and O'Kennedy, 2007). In order to understand how calcium impacts the dynamic of formation of aggregates in the model, it is useful to look at the size distributions. For three concentrations of CaCl2 (5.0, 6.0 and 7.0 mM), the simulated size distributions after holding are presented on Fig. 5. Comparison between heat treatment with 5.0 and 6.0 mM, shows that the models predict relatively similar aggregate size distributions but the concentrations is much higher with 6.0 mM. For the heat treatment with 7.0 mM of CaCl2, the mean size is much larger and there are large differences between the bulk distribution obtained from the PBE/empirical RTD model and the distribution given by the PBE/plug flow model. The latter shows aggregates size ranging from 0.03 mm to 20 mm with a median around 6 mm. The distribution of the PBE/empirical RTD model is very wide, ranging from 0.03 mm to about 500 mm. An analysis of the values of stability factors W1 and W3 (Table 1) can help to understand the effect of calcium chloride on the aggregation process. W1 expresses primary protein aggregation trough disulphide bonds to form small aggregates and W3 represents the stability ratio for large insoluble aggregates. It can be seen on Table 1 that the formation of large aggregates (W3) is more impacted by the calcium concentration than the formation of small aggregates (W1). Indeed, for calcium concentrations varying from 5.0 mM to 7.0 mM, the value of W3 lowers from 413.6 to 34.6 representing a decrease of 92% while W1 decreases of only 42% (from 3.17 109 to 1.83 109). This result means that it is more probable to form large aggregates at 7 mM of calcium concentration than at 5 mM, as evidenced on Fig. 5. This outcome may be explained by the fact that at low concentrations, charge shielding is the predominant effect of ionic calcium on protein interactions. Aggregation is promoted either by increasing
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
39
Fig. 6. Simulation with plug flow model and PBE/RTD combined model of the volume fraction of aggregates obtained by pilot heat treatment (20 l h1, To ¼ 87 C, 6% b-lg) with different calcium chloride concentrations.
fraction of native proteins is lowered while the average size and volume fraction of aggregates increase, meaning that calcium chloride promotes b-lg aggregation. The slope of the curve changes dramatically between 6.4 and 6.6 mM of calcium chloride. An objective of the present study is to compare the predictions of the two models in order to determine whether RTD has a significant effect on the selected characteristics of the aggregates suspension. Regarding the aggregates volume fraction, Fig. 6 shows an appreciable difference between the two models for the higher calcium chloride concentrations (6.6 and 7.0 mM), meaning that the effect of RTD on aggregation depends on the calcium concentration. This ending may be explained by the fact that few transformations occurred at the low concentrations of CaCl2 (5.0, 6.0 and 6.4, mM) as evidenced by the small size of aggregates (<1 mm) shown in Fig. 7b. Differences between plug flow and empirical RTD become important when the product is subjected to important transformations (i.e. with 6.6 and 7.0 mM of CaCl2). In these cases, considering the residence time distribution may be crucial in order to have a good representation of the aggregation phenomena. 4.2. Comparison of model prediction with experimental results
Fig. 5. Simulated volume size distributions for the plug flow model and the bulk distribution of the PBE/RTD combined model after holding (20 l h1, To ¼ 87 C, 6% blg) for different CaCl2 concentrations: 5.0 mM (a), 6.0 mM (b) and 7.0 mM (c).
electrostatic attraction between negatively charged proteins, or by decreasing electrostatic repulsion between positively charged proteins, forming essentially small aggregates. But, at high concentrations, in addition to charge shielding effects, large aggregates are formed by intermolecular ion binding to form protein-Caprotein complexes (Parris et al., 1993; Chi et al., 2003; Westerik et al., 2015),. Table 1 also shows a significant gap between values of W3 at 6.4 mM and 6.6 mM suggesting a threshold effect of calcium, probably due to the existence of a stoichiometric calcium/b-lg ratio (Mounsey and O'Kennedy, 2007). This threshold effect is evidenced in Fig. 6 which presents the aggregates volume fraction as function of the concentration of CaCl2 for both models. It can be seen that the calcium concentration affects the predicted properties of the suspensions. When calcium concentration increases, the
By combining the PBE model to the empirical RTD model, one expects a better prediction of the aggregate properties compared to the coupling with plug flow model. This section describes the comparison of experimental and simulated data on the basis of properties related to the functionalities of the suspensions: residual native fraction, viscosity of the suspension and aggregate size distributions measured with laser granulometer (range 0.1e1000 mm). 4.2.1. Functional properties of the aggregate suspensions The predicted (by both models) and experimental results of two aggregates functionality indicators (residual native fraction and median volume diameter) obtained for a heat treatment of 87 C/ 49 s and for different CaCl2 concentration are compared in Fig. 7. Both models overestimate the residual native fraction compared to experimental data (Fig. 7a) even if the general trend regarding the influence of CaCl2 is well represented. This discrepancy between experimental and modelling results may be explained by the time-temperature evolution of the different fractions of fluids inside the heat exchanger. Indeed the fraction of fluids flowing near the heat exchanger wall is submitted to a higher temperature and to a longer residence time than the others, resulting in a higher level of denaturation. On the other hand, fluid parcels moving on the tube axis undergoes low unfolding as a result of relatively lower
40
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
model and well predicted by the PBE/empirical RTD model. This result put in evidence, once again, the importance to account for the histories of the different fractions of fluids when aggregation is high (high level of calcium concentration) in order to represent satisfactorily the aggregation process. 4.2.2. Volume size distributions Using Population Balance Model make it possible to predict aggregate size distributions which are the key feature of protein aggregates functionalities. Comparisons between predicted and experimental cumulative size distributions for three concentrations of CaCl2 (5.0 mM, 6.0 mM and 7.0 mM) are presented on Fig. 8. For 5.0 and 6.0 mM of CaCl2 shown on Fig. 8a and b respectively, both models predict distributions in good agreement with experimental results and consistent with average volume size presented
Fig. 7. Experimental and simulated (plug flow model and PBE/RTD combined model) aggregate suspension properties after pilot heat treatment (20 l h1, To ¼ 87 C for 49 s) - Residual native fraction (a) and Median volume diameter (b).
temperature and shorter residence time than the other fluid fractions. Fluid trajectories between the axis and the wall have intermediates unfolding levels. Actually, samples recovered at the outlet of the cooling section constitute a blending of these different fractions of fluid, ensuing a compensation effect between the different unfolding levels. It is not surprising that the plug flow model fails to represent this effect because it considers that all fractions of fluids have the same residence time and timetemperature evolution. Contrariwise, the bulk data derived from the PBE/RTD model considers a blend of different fluid parcels. These fluid parcels are assumed to have different residence time but follow the same temperature profile (vs position). The fact that fluid parcels experience different temperature evolutions depending on distance from the heating wall is not taken into account on this approach. This could explain the observed discrepancy. These findings join those of Chantoiseau et al. (2012) who demonstrated that dynamical and thermal histories associated with fluid parcels which progress more or less quickly, far or close to the heating wall determine the level of aggregation of WP. In the region near the heat exchanger wall, the product is heated for a long time, resulting in important aggregation; by the contrary, in the axial region where the residence time is shorter, the product is less heated and weakly transformed. The authors pointed out the importance to take into account every fluid parcels trajectories instead of considering a plug flow model, in order to accurately predict the b-lg aggregation inside the heat exchanger section. Fig. 7b shows good agreement between experiments and both modelling results of median volume diameter for the three smallest concentrations of CaCl2. For the two higher concentrations, the median volume diameter is underestimated by the PBE/plug flow
Fig. 8. Experimental and simulated (plug flow model and PBE/RTD combined model) aggregate size distribution after pilot heat treatment (20 l h1, To ¼ 87 C for 49 s) with different calcium chloride concentrations: 5 mM (a), 6 mM (b) and 7 mM (c).
N. Erabit et al. / Journal of Food Engineering 177 (2016) 31e41
on Fig. 7b. For a calcium concentration of 7.0 mM, the plug flow model largely underestimates the aggregates size. As can be seen on Fig. 8c and 90% of aggregates predicted by the plug flow model have a size below 50 mm while this proportion is only of 30% for the experimental results. This could be attributed to the fact that formation of large aggregates in the slowest fluid parcels is not taken into account. The second model which considers RTD predicts higher concentrations of large aggregates compared to plug flow model (only about 48% of aggregates have a size below 50 mm). Fig. 8c shows that the PBE/RTD combined model predicts aggregates volume size distribution close to the experimental result. The median volume size is in good agreement with the experimental values. However, one can evidenced the fact that the experimental distribution shape is not exactly predicted. The slight dissimilarity between experimental and simulated distributions may be explained by the fact that the model does not account some other physico-chemical factors, such as pH, which can have a significant impact on aggregation process, especially on the formation of large ~ iga et al., 2010). aggregates (Zún Finally, the PBE/RTD combined model is more realistic than the plug flow model and gives better prediction for aggregates size distribution. 5. Conclusion A model for the prediction of the b-lg aggregation during the heat treatment in a continuous pilot device has been presented. It couples a RTD model and a colloidal aggregation model previously developed. This model has been compared to a plug flow model in which only the average time/temperature history is considered. The integration of RTD appears to be very important when the product is subject to much transformation (formation of large aggregates > 10 mm). The predictions were compared to the experimental data. The characteristics of the suspensions for experiments carried with the lowest CaCl2 concentrations (5.0, 6.0 and 6.4 mM) are well predicted. For higher concentrations of calcium (6.6 and 7.0 mM), the PBE/RTD combined model well predicts the average aggregates size. Further developments should interestingly integrate the influence of physico-chemical variables other than calcium concentration in order to better describe the distributions of aggregates size. The modelling approach combining PBE to RTD developed in this study is valuable for enhancing the performance of WP aggregation process. Such a hybrid model is essential in process design, evaluation and analysis. Acknowledgements This study was carried out within the framework of the Globule project (ANR-08-ALIA-08) with the help of the other partners: INRA - PIHM and GEA Process Engineering France. References Anema, S.G., McKenna, A.B., 1996. Reaction kinetics of thermal denaturation of whey proteins in heated reconstituted whole milk. J. Agric. Food Chem. 44 (2), 422e428. Arellano, M., Benkhelifa, H., Alvarez, G., Flick, D., 2013. Coupling population balance and residence time distribution for the ice crystallization modeling in a scraped surface heat exchanger. Chem. Eng. Sci. 102 (0), 502e513. Arosio, P., Beeg, M., Nicoud, L., Morbidelli, M., 2012. Time evolution of amyloid fibril length distribution described by a population balance model. Chem. Eng. Sci. 78 (0), 21e32. Bas¸, N., Catak, M., Zumaeta, N., Fitzpatrick, J.J., Cronin, K., Byrne, E.P., 2009. Population balance modelling of protein precipitates exhibiting turbulent flow through a circular pipe. Chem. Eng. Sci. 64 (18), 4051e4059. Bilgili, E., Scarlett, B., 2005. Numerical simulation of open-circuit continuous mills using a non-linear population balance framework: incorporation of non-first-
41
order effects. Chem. Eng. Technol. 28 (2), 153e159. Chantoiseau, E., Plana-Fattori, A., Doursat, C., Flick, D., 2012. Coupling fluid flow, heat transfer and thermal denaturation-aggregation of beta-lactoglobulin using an Eulerian/Lagrangian approach. J. Food Eng. 113 (2), 234e244. Chi, E., Krishnan, S., Randolph, T., Carpenter, J., 2003. Physical stability of proteins in aqueous solution: mechanism and driving forces in nonnative protein aggregation. Pharm. Res. 20 (9), 1325e1336. Dannenberg, F., Kessler, H.G., 1988. Reaction kinetics of the denaturation of whey proteins in milk. J. Food Sci. 53 (1), 258e263. Elofsson, U.M., Dejmek, P., Paulsson, M.A., 1996. Heat-induced aggregation of betalactoglobulin studied by dynamic light scattering. Int. Dairy J. 6 (4), 343e357. Erabit, N., Flick, D., Alvarez, G., 2013. Effect of calcium chloride and moderate shear on b-lactoglobulin aggregation in processing-like conditions. J. Food Eng. 115 (1), 63e72. Erabit, N., Flick, D., Alvarez, G., 2014. Formation of b-lactoglobulin aggregates during thermomechanical treatments under controlled shear and temperature conditions. J. Food Eng. 120 (0), 57e68. Erabit, N., Ndoye, F.T., Alvarez, G., Flick, D., 2015. A population balance model integrating some specificities of the b-lactoglobulin thermally-induced aggregation. J. Food Eng. 144 (0), 66e76. Heath, A.R., Koh, P.T.L., 2003. Combined population balance and CFD modelling of particle aggregation by polymeric flocculant. In: Third International Conference on CFD in the Minerals and Proocess Industries, Melbourne, Australia, pp. 339e344. Lee, K.F., Mosbach, S., Kraft, M., Wagner, W., 2015. A multi-compartment population balance model for high shear granulation. Comput. Chem. Eng. 75 (0), 1e13. Lian, G., Moore, S., Heeney, L., 2006. Population balance and computational fluid dynamics modelling of ice crystallisation in a scraped surface freezer. Chem. Eng. Sci. 61 (23), 7819e7826. Mounsey, J.S., O'Kennedy, B.T., 2007. Conditions limiting the influence of thioledisulphide interchange reactions on the heat-induced aggregation kinetics of b-lactoglobulin. Int. Dairy J. 17 (9), 1034e1042. Ndoye, F.T., Erabit, N., Alvarez, G., Flick, D., 2012. Influence of whey protein aggregation on the residence time distribution in a tubular heat exchanger and a helical holding tube during heat treatment process. J. Food Eng. 112 (3), 158e167. Ndoye, F.T., Erabit, N., Flick, D., Alvarez, G., 2013. In-line characterization of a whey protein aggregation process: aggregates size and rheological measurements. J. Food Eng. 115 (1), 73e82. Nicolai, T., Britten, M., Schmitt, C., 2011. b-Lactoglobulin and WPI aggregates: formation structure and applications. Food Hydrocoll. 25, 1945e1962. Parris, N., Anema, S.G., Singh, H., Creamer, L.K., 1993. Aggregation of whey proteins in heated sweet whey. Journal of Agricultural and Food Chemistry J. Agric. Food Chem. 41, 460e464. Prasetya, A., Liu, L., Litster, J., Watanabe, F., Mitsutani, K., Ko, G.H., 1999. Dynamic model development for residence time distribution control in high-impact polypropylene copolymer process. Chem. Eng. Sci. 54 (15e16), 3263e3271. Sen, M., Ramachandran, R., 2013. A multi-dimensional population balance model approach to continuous powder mixing processes. Adv. Powder Technol. 24 (1), 51e59. Sen, M., Singh, R., Vanarase, A., John, J., Ramachandran, R., 2012. Multi-dimensional population balance modeling and experimental validation of continuous powder mixing processes. Chem. Eng. Sci. 80 (0), 349e360. Simmons, M.J.H., Jayaraman, P., Fryer, P.J., 2007. The effect of temperature and shear rate upon the aggregation of whey protein and its implications for milk fouling. J. Food Eng. 79 (2), 517e528. Simons, J.-W.F.A., Kosters, H.A., Visschers, R.W., de Jongh, H.H.J., 2002. Role of calcium as trigger in thermal b-lactoglobulin aggregation. Arch. Biochem. Biophys. 406 (2), 143e152. Smoluchowski, M., 1917. Versuch einer mathematischen Theorie der Koagula€sungen. Z. Phys. Chem. 92, 129e168. tionskinetic kolloider, Lo Spiegel, T., Huss, M., 2002. Whey protein aggregation under shear conditions - effects of pH-value and removal of calcium. Int. J. Food Sci. Technol. 37, 559e568. Tanguay, M., Fawell, P., Adkins, S., 2014. Modelling the impact of two different flocculants on the performance of a thickener feedwell. Appl. Math. Model. 38 (17e18), 4262e4276. Tolkach, A., Kulozik, U., 2007. Reaction kinetic pathway of reversible and irreversible thermal denaturation of b -lactoglobulin. Lait 87 (4e5), 301e315. Vale, H.M., McKenna, T.F., 2005. Modeling particle size distribution in emulsion polymerization reactors. Prog. Polym. Sci. 30 (10), 1019e1048. Weller, K.R., Spencer, S.J., Gao, M., Liu, Y., 2000. Tracer studies and breakage testing in pilot-scale stirred mills. Miner. Eng. 13 (4), 429e458. Westerik, N., Scholten, E., Corredig, M., 2015. The effect of calcium on the composition and physical properties of whey protein particles prepared using emulsification. Food Chem. 177, 72e80. Wynn, E.J.W., Hounslow, M.J., 1997. Integral population balance equations for growth. Chem. Eng. Sci. 52 (5), 733e746. Zumaeta, N., Byrne, E.P., Fitzpatrick, J.J., 2010. Application of CFD and breakage modelling for predicting the size reduction of protein precipitates during transport. Chem. Eng. Res. Des. 88 (7), 911e917. Zumaeta, N., Cartland-Glover, G.M., Heffernan, S.P., Byrne, E.P., Fitzpatrick, J.J., 2005. Breakage model development and application with CFD for predicting breakage of whey protein precipitate particles. Chem. Eng. Sci. 60 (13), 3443e3452. ~ iga, R.N., Tolkach, A., Kulozik, U., Aguilera, J.M., 2010. Kinetics of formation and Zún physicochemical characterization of thermally-induced b-lactoglobulin aggregates. J. Food Sci. 75 (5), E261eE268.