Predicting the Effect of the Homogenization Pressure on Emulsion

Mar 31, 2011 - pressure homogenizer that incorporated multiple-drop formation within two .... effect of homogenization pressure on the drop-size distr...
1 downloads 0 Views 2MB Size
ARTICLE pubs.acs.org/IECR

Predicting the Effect of the Homogenization Pressure on Emulsion Drop-Size Distributions Neha B. Raikar,† Surita R. Bhatia,† Michael F. Malone,† David Julian McClements,‡ and Michael A. Henson*,† †

Department of Chemical Engineering and ‡Department of Food Science, University of Massachusetts, Amherst, Massachusetts 01003-9303, United States ABSTRACT: We have previously developed a population balance equation (PBE) model for emulsion drop breakage in a highpressure homogenizer that incorporated multiple-drop formation within two mechanisms of turbulent drop breakage. The model was found to satisfactorily predict the effects of formulation variables on the drop-size distribution, but the model was not extensible to a range of homogenization pressures. The objective of this paper is to determine the additional model elements necessary to obtain acceptable predictions over a wide range of pressures. The most significant improvements were obtained by increasing the number of daughter drops formed upon breakage from 20 to 150 drops and by introducing a maximum stable diameter, below which drops could not break. Smaller improvements were obtained by introducing terms that describe the loss of energy available for drop breakage due to thermal heating of the sample and homogenizer and by extending the model to account for the effects of surfactant adsorption and deficiency on the interfacial tension. The simultaneous implementation of all five enhancements was shown to produce a 62% improvement over the previous model, as measured by a least-squares objective based on the difference between the measured and predicted drop-size distributions over five homogenization passes and five pressures in the range 2501250 bar. The resulting model was also validated over a range of oil and surfactant concentrations and shown to provide satisfactory predictions.

1. INTRODUCTION Emulsions constitute a broad class of structured products with a wide variety of applications including foods, petroleum recovery, cosmetics, polishes, road surfacing, and agricultural sprays. The term “structured products” is appropriate because of the presence of an internal structure at a micrometer or nanometer scale. The structure is dictated not only by the constitute ingredients and their properties but also by the processing conditions used for manufacturing. The final product quality attributes like the appearance, viscosity, partitioning rate of active ingredients, etc., are linked to the microstructure. For instance, in the case of margarine, the viscosity, texture, and mouthfeel are dictated by the size of the fat droplets and fat crystals.1 Commonly used techniques for emulsification require the application of mechanical energy to the two immiscible phases and stabilization of the newly formed interfaces by surfactants. Various types of processing equipment ranging from low-intensity devices (e.g., high-speed blenders) to high-intensity devices (e.g., high-pressure homogenizers and microfluidizers) are available for emulsion preparation. Emulsification can be a single-step or multistep process depending on the application, and the choice of equipment is based on the energy intensities required for the various steps.2 Oil-in-water emulsions are typically formed by first preparing a coarse emulsion using a low-intensity statorrotor-type device that mixes the various ingredients into a stable form. This premix is then processed with a high-intensity device, such as a high-pressure homogenizer, for a single pass or multiple passes until the target emulsion is obtained. Controlling the emulsion microstructure, in addition to the composition, is crucial to meeting the quality specifications. As a consequence, there is considerable motivation to develop robust predictive models for the homogenization process. r 2011 American Chemical Society

Population balance equation (PBE) models have been widely used to describe homogenization processes for oil-in-water emulsion preparation.37 PBE models allow prediction of the evolution of the drop-size distribution through specification of two rate processes: (1) breakage of drops due to the turbulent flow field and subsequent drop stabilization by adsorption of the surfactant at the oilwater interface and (2) coalescence of colliding drops due to insufficient free surfactant in solution. Under conditions that promote minimal coalescence such as small oil-to-surfactant ratios, the PBE model can be reduced to the breakage-only form, which has been widely used. We previously developed breakage-only PBE models based on binary drop breakage5 and multiple drop breakage8 to predict drop-size distributions of emulsions prepared with high-pressure homogenization. The models included mechanistic functions for two assumed breakage mechanisms (turbulent inertial and viscous forces) to allow prediction of the drop-size distributions over a range of formulation variables (e.g., oil and surfactant concentrations) and homogenization pressures. We used nonlinear optimization to determine unknown, adjustable parameters in the two breakage functions from measured drop-size distributions collected at a single set of formulation and homogenization conditions and then tested the extensibility of the parametrized model to different conditions. While the models proved to be extensible to different formulations,5,8 the binary breakage model yielded substantial overprediction of the breakage rate at

Received: August 30, 2010 Accepted: March 31, 2011 Revised: March 17, 2011 Published: March 31, 2011 6089

dx.doi.org/10.1021/ie101818h | Ind. Eng. Chem. Res. 2011, 50, 6089–6100

Industrial & Engineering Chemistry Research

ARTICLE

Table 1. Oil and Surfactant Concentrations Used oil

surfactant

oil

surfactant

concn (wt %)

concn (wt %)

concn (wt %)

concn (wt %)

5

1

20

4

10

2

pressures other than the single pressure used for parameter estimation.5 Our previous studies indicate that the developed PBE models5,8 lack some features required to accurately predict the effect of homogenization pressure on the drop-size distributions. Substantial increases in the emulsion temperature following multiple passes at high pressure observed in our laboratory and by other researchers9 suggest that a considerable amount of mechanical energy may be dissipated as thermal energy and not available for drop breakage. Because binary breakage models approximate multiple drop breakage as a series of binary breakage events, PBE models that include the turbulent breakage of mother drops into multiple daughter drops3,10,11 as supported by laminar-flow experiments1214 offer the potential for improved pressure predictions. Because of the limits on the mechanical energy that can be generated, there exists a pressure-dependent maximum stable diameter, below which drops are stable to breakage, and the incorporation of this feature3,11,15 offers the potential to reduce overprediction of the pressure effect. Because surfactant depletion is known to impact the resulting distributions,2,7 the explicit incorporation of a surfactant balance and the relationship between free surfactant concentration and the oilwater interfacial tension could improve predictions of the pressure effect. The objective of this paper was to determine which combination of these additional model features was necessary to obtain acceptable drop-size-distribution predictions over a wide range of pressures.

2. MATERIALS AND METHODS 2.1. Experimental Methods. Materials. We used soybean oil (Spectrum Organic) as the dispersed phase, nanopure water as the continuous phase, and nonionic surfactant Pluronic F-68 (Sigma Aldrich) as the emulsifier. Nanopure water was prepared using a standard system (Nanopure Infinity, Barnstead) purified to 18 mΩ cm1. The base-case emulsion consisted of 5 wt % oil and 1 wt % surfactant, with the remainder water. A relatively small oil-to-surfactant ratio of 5 wt %/wt % was used to reduce the possibility of coalescence. Emulsions were also prepared at different concentrations while keeping the ratio of the oil-to-surfactant added constant at the base-case value of 5 wt %/wt %, although the actual oil-to-surfactant ratio in solution might be slightly different because of the surfactant solubility in water. The different oil and surfactant concentrations used are listed in Table 1. Emulsion Preparation. Emulsions were produced using a twostep process. First, a coarse premix was prepared by mixing the different ingredients in a statorrotor device (Ultra-Turrax model T25, Rose Scientific Ltd.) at 16 000 rpm for 15 min. A standard premix was made at 40 wt % oil and 8 wt % surfactant and then diluted to the concentrations in Table 1 to ensure similar initial drop-size distributions for all experiments. About 100 mL of the premix was then processed in a high-pressure homogenizer (Emulsiflex C-3, Avestin Inc.) to prepare the fine

emulsion. The homogenizer had a volumetric flow rate of approximately 3 L h1. Multiple passes were performed by collecting and reprocessing the emulsion from the previous pass. After each pass, the emulsion was cooled to room temperature to minimize the effects of heating on the emulsion properties. Samples from the premix and all five passes were subjected to particle size analysis (see below) with approximately 5 mL of the sample from each pass used for the analysis. Experiments were performed for five homogenization pressures ranging from 250 to 1250 bar in steps of 250 bar. Emulsion Characterization. Drop-size distributions of the premix and the five passes were measured using a static lightscattering device (Malvern Mastersizer S, Malvern). Previously, we found such drop-size distributions to be very reproducible,5,8 with errors much less than the difference between model predictions and experimental data. Oilwater interfacial tensions were measured by drop-shape analysis (tensiometer model DSA-10, Kruss Instruments) at 25 C. The densities of the emulsion samples were measured using a densitometer (model DE50, Mettler Toledo). The viscosities of the constituent phases were measured using a stress-controlled rheometer (TA-Ar 2000, Texas Instruments). 2.2. Theory. PBE Model. The PBE model was formulated for a pure breakage process under the assumption that the coalescence was negligible because of the small oil-to-surfactant ratio used.5,8 In this case, the PBE can be written as1618 Z ¥ ∂nðv, tÞ ¼  gðvÞ nðv, tÞ þ βðv, v0 Þ gðv0 Þ nðv0 , tÞ dv0 ð1Þ ∂t 0 where n(v,t) dv is the number of drops in the volume range [v, v þ dv] per unit volume of the dispersion at time t, g(v) is the breakage rate representing the fraction of drops of volume v breaking per unit time, and β(v,v0 ) is the daughter drop distribution function representing the probability of forming a daughter drop of volume v from breakage of a mother drop of volume v0 . The high-pressure homogenizer was modeled as a well-mixed batch system, with each pass corresponding to one dimensionless time unit. The PBE (eq 1) was reformulated in terms of the volume percent distribution np(v,t) to allow a direct comparison with measured distributions from the particle-size analyzer Z ¥ 0 ∂np ðv, tÞ gðv Þ βðv, v0 Þ np ðv0 , tÞ 0 ¼  gðvÞ np ðv, tÞ þ v dv ∂t v0 0 ð2Þ The initial condition np(v,0) was taken to be the measured dropvolume distribution of the premix. Following our previous work,5,8 the breakage rate was specified as the sum of two individual breakage functions that accounted for drop breakage due to inertial [g1(v)] and viscous [g2(v)] forces under turbulent conditions: g(v) = g1(v) þ g2(v). The turbulent inertial breakage rate g1(v) was based on a previously proposed function18 that was modified for highpressure homogenization through the appropriate specification of the energy dissipation rate ɛ (see below) " # K2 σð1 þ jÞ2 2=9 1=3 g1 ðvÞ ¼ K1 v E exp  ð3Þ Fd v5=9 E2=3 where σ is the oilwater interfacial tension, φ is the dispersedphase volume fraction, Fd is the density of the dispersed phase, and 6090

dx.doi.org/10.1021/ie101818h |Ind. Eng. Chem. Res. 2011, 50, 6089–6100

Industrial & Engineering Chemistry Research

ARTICLE

K1 and K2 are adjustable constants to be determined by parameter estimation. The second breakage rate g2(v) was derived considering breakage of drops due to turbulent viscous forces5 !  1=2  1=2 2 EFd K4 σ2 λ g2 ðvÞ ¼ K3 exp ð4Þ π ηd ηc v2=3 EFd where ηd and ηc are the viscosities of the dispersed and continuous phases, respectively, λ = ηc/ηd, and K3 and K4 are adjustable constants to be determined by parameter estimation. In our previous studies,5,8 we utilized a simple expression for the energy dissipation rate ɛ. In this study, we adopted the following more detailed description3,15 with the potential for improving predictions of drop-size distributions over a range of homogenization pressures E¼

ΔPQ Vdiss

ð5Þ

where ΔP is the applied pressure, Q is the volumetric flow rate, and Vdiss is the valve gap volume calculated as π Vdiss ¼ ðDo 2  Di 2 Þhgap ð6Þ 4 Table 2. Values of the Constants K1K4 for the Base Case K1 K2

2.62  108 1

1.75  10

K3

1.19  106

K4

3.80  103

The following equation relating the applied pressure to the valve gap hgap has been reported:19 !2 !2 Fc Q Fc Q ΔP ¼  4 πDi hgap 2 πDo hgap ! 5Fc  hgap 3

ηc Fc

!3=5 

Q 2π

7=5 " 2=5  2=5 # Di Do  2 2 ð7Þ

An equation for the valve gap was obtained after retaining only the last term of eq 7 and rearranging in terms of hgap. We observed that, for our geometry, similar values of hgap were obtained with or without the previous terms. Hence, only the last term was retained to obtain an explicit equation for hgap (eq 8). 

hgap ¼

5Fc ΔP

1=3

ηc Fc

!1=5 

Q 2π

7=15 " 2=5  2=5 #1=3 Di Do  2 2

ð8Þ The daughter drop distribution function β(v,v0 ) in the PBE (eq 2) was specified to allow breakage of a single mother drop into more than two daughter drops. To this end, the function was considered to be the power law product form of the generalized

Figure 1. Base-case modeling results. Experimental and predicted drop-volume distributions for the (P) premix, (1) first pass, (3) third pass, and (5) fifth pass. Ψp is the objective function value for the particular pressure set: (a) 250 bar pressure set; (b) 750 bar pressure set; (c) 1250 bar pressure set. (d) Pass-by-pass objective function Ψi as a function of the number of passes for the different pressures. 6091

dx.doi.org/10.1021/ie101818h |Ind. Eng. Chem. Res. 2011, 50, 6089–6100

Industrial & Engineering Chemistry Research

ARTICLE

Figure 2. Modeling results when heating loss was added to the base-case model. Experimental and predicted drop-volume distributions for the (P) premix, (1) first pass, (3) third pass, and (5) fifth pass. Ψp is the objective function value for the particular pressure set. (a) Temperature increase averaged over the five passes (ΔT) and effective homogenization pressure as a function of the pressure applied to the homogenizer: (b) 250 bar pressure set; (c) 750 bar pressure set; (d) 1250 bar pressure set.

HillNg distribution,2022  q  1   p v v r1 βðv, v0 Þ ¼ 1  , Bðq, rÞ v0 v0

r ¼ qðp  1Þ ð9Þ

where p g 2 is the number of daughter drops formed, which can be considered as an adjustable parameter, B(q,r) is the β function, and q > 0 is a parameter that determines the shape of the distribution function. For q = 1, the following simplified equation that represents a uniform distribution for p daughter drops is obtained:   v p2 0 ð10Þ βðv, v Þ ¼ pðp  1Þ 1  0 v All of the results reported in this paper are based on the uniform distribution function with p = 20 unless otherwise stated (eq 10). We also considered a PBE model for combined drop breakage and coalescence to investigate the possible significance of the coalescence phenomenon in our experimental system. Expressed in terms of the number distribution n(v,t), the PBE is written as Z ¥ dnðv, tÞ ¼  gðvÞ nðv, tÞ  βðv, v0 Þ gðv0 Þ nðv0 , tÞ dv0 dt 0 Z ¥  nðv, tÞ Cðv, v0 Þ nðv0 , tÞ dv0 0 Z 1 ¥ þ Cðv  v0 , v0 Þ nðv0 , tÞnðv  v0 , tÞ dv0 ð11Þ 2 0

where C(v,v0 ) is the coalescence frequency, which is a product of the collision frequency h(v,v0 ) and the coalescence efficiency λ(v,v0 ). These functions were described by18 E1=3 2=3 ðv þ v02=3 Þðv2=9 þ v02=9 Þ1=2 1þj 2 !4 3 1=3 01=3 Cμ F E v v 5 λðv, v0 Þ ¼ exp4  K6 c 2 c σ v1=3 þ v01=3

hðv, v0 Þ ¼ K5

ð12Þ

ð13Þ

where K5 and K6 are adjustable constants to be determined by parameter estimation. The PBE models for breakage only (eq 2) and combined breakage and coalescence (eq 11) were solved numerically by discretizing the model equations in volume space with 100 node points using the fixed pivot method,23 which conserves both mass and volume. The resulting system of 100 nonlinear ordinary differential equations describing the time evolution of the drop-volume distribution at each node point was solved using the Matlab integration code ode45. Heating Loss. The substantial increases in the emulsion temperature observed following multiple homogenization passes at high pressure9 suggest that a considerable amount of mechanical energy may be dissipated as thermal energy and may not be available for drop breakage. Energy is lost both to the liquid sample and to the homogenization equipment. This effect was neglected in our previous models.5,8 The thermal energy gained EH was estimated as EH ¼ Fs ½jCpd þ ð1  jÞCpc ΔT þ Fc Cpc ΔTins 6092

ð14Þ

dx.doi.org/10.1021/ie101818h |Ind. Eng. Chem. Res. 2011, 50, 6089–6100

Industrial & Engineering Chemistry Research

ARTICLE

surfactant concentration in the bulk solution c, both of which can be obtained experimentally from measurements of the oilwater interfacial tension for various amounts of surfactant Γ¼ 

1 dσ aRT d lnðcÞ

ð17Þ

where a is 1 for nonionic surfactants and 2 for ionic surfactants. Thus, Γ can be estimated from a plot of the first-order derivative of σ versus ln(c).2 The Gibbs adsorption isotherm also allows expression of the surface pressure Π in terms of Γ Z c Π ¼ σ o=w  σ ¼ aRT ΓðcÞ d ln c ð18Þ 0

The Langmuir adsorption isotherm relates Γ to the bulk concentration c c c1=2 Γ ¼ ð19Þ c Γ¥ 1þ c1=2

Figure 3. (a) Measured oilwater interfacial tension (dots) and quadratic fit (line) as a function of the bulk surfactant concentration. (b) Normalized surface load as a function of the bulk surfactant concentration.

where Fs is the sample density measured using the densitometer, Cpd and Cpc are the heat capacities of the oil and water, respectively, ΔT is the total temperature rise of the sample measured from the inlet to the outlet of the homogenizer over all of the passes, and ΔTins is the temperature rise of the instrument measured by passing water at room temperature through the homogenizer and measuring the outlet water temperature after each pass. When EH is known, the effective pressure available for drop disruption can be calculated as follows: ΔPeff ¼ ΔP  EH

ð15Þ

Variable Interfacial Tension. The interfacial tension σ increases as the surfactant is depleted from the solution, and less surfactant is available for adsorption, thereby reducing the stability of newly formed drops. For simplicity, our previous homogenizer models5,8 treated σ as a constant determined from oilwater interfacial tension measurements. To account for this effect, we considered pass-by-pass recalculation of σ for incorporation into the breakage rate functions (eqs 3 and 4). The reduction in the interfacial tension due to surfactant is referred to as the surface pressure Π2 Π ¼ σo=w  σ

ð16Þ

where σo/w is the interfacial tension in the absence of surfactant, i.e., pure oilwater. The Gibb’s adsorption isotherm was used to relate the surface load (Γ), which is the amount of surfactant present at the interface, to the interfacial tension and the

where Γ¥ is the surface load when the interface is fully covered with surfactant found from the adsorption isotherm as the value of Γ at the highest concentration used and c1/2 is the surfactant concentration when Γ = 1/2Γ¥. When eq 19 is substituted into eq 18, an expression relating the interfacial tension to the bulk concentration is obtained. ! c σ ¼ σ o=w  aRTΓ¥ ln 1 þ ð20Þ c1=2 The adsorbed surfactant concentration was calculated as2 cads ¼

6Γj d32

ð21Þ

where the Sauter mean diameter d32 was calculated from the appropriate moments of the predicted drop-volume distribution np(v). The bulk surfactant concentration c was calculated by subtracting the concentration of the surfactant adsorbed cads from the initial concentration of that added cinit. c ¼ cinit  cads

ð22Þ

Initially, σ was set equal to the value corresponding to cinit. For each pass, the adsorbed surfactant concentration cads was calculated from eq 21 and the bulk surfactant concentration c was calculated from eq 22. The updated c value was used to recalculate the interfacial tension σ using eq 20. The updated σ value was used in the breakage rate functions (eqs 3 and 4), with increasing σ values resulting from surfactant depletion, leading to decreased values of the breakage rate g(v). Critical Drop Diameter for Breakage. Because of the limits on the mechanical energy that can be generated, there exists a pressure-dependent maximum stable diameter, below which drops are stable to breakage. Our previous models5,8 have neglected this phenomenon under the assumption that the breakage rate would be sufficiently small at small drop volumes to effectively incorporate this effect. The maximum stable diameter dmax for turbulent inertial breakup is15 dmax ¼ c1 6093

σ3=5 E2=5 Fc 1=5

ð23Þ

dx.doi.org/10.1021/ie101818h |Ind. Eng. Chem. Res. 2011, 50, 6089–6100

Industrial & Engineering Chemistry Research

ARTICLE

Figure 4. Modeling results when both heating loss and variable surface tension were added to the base-case model. Experimental and predicted dropvolume distributions for the (P) premix, (1) first pass, (3) third pass, and (5) fifth pass. Ψp is the objective function value for the particular pressure set. (a) Calculated interfacial tension increase as a function of the pass number and the homogenizer pressure: (b) 250 bar pressure set; (c) 750 bar pressure set; (d) 1250 bar pressure set.

where c1 is an order unity proportionality constant estimated from experimental data. Below dmax, drops cannot be broken, but they can be formed by breakup of larger drops. Because measured drop-volume distributions changed only slightly between the fourth and fifth passes, the Sauter mean diameter d32 for the fifth pass was considered to be a steady-state value that represented the smallest achievable diameter for a particular experiment. The proportionality constant c1 in eq 23 was estimated by performing linear regression with dmax equal to the fifth pass d32. The data available for regression were the energy dissipation rate ɛ at the five equally spaced pressures in the range 2501250 bar and the fifth pass d32 measured for each experiment. The breakage rate was modified as g(v) = 0 for d < dmax. Complete Surfactant Depletion. If there remains no free surfactant in solution because all of the surfactant added has adsorbed onto drop interfaces, new drops formed by breakage cannot be stabilized and they will immediately recoalesce. Our previous models5,8 have neglected this effect under the assumption that there would always be sufficient surfactant for stabilization of newly formed drops. We calculated the free surfactant concentration c using eq 22. The breakage rate g(v) was set to zero for any pass where c e 0. While simple, this method does have the disadvantage of introducing discontinuities into the model through g(v). Parameter Estimation. The constants K1K4 in the breakage rate functions (eqs 3 and 4) and the constants K5 and K6 in the coalescence frequency functions (eqs 11 and 12) were treated as adjustable parameters to be estimated from homogenization

experiments. Available experimental data for parameter estimation included emulsion ingredient properties (φ, σ, Fd, Fc, ηd, and ηc) that were measured a priori and drop-size distributions for the premix np(v,0) and the processed emulsion obtained after the ith homogenizer pass np(v,i). The PBE models (eqs 2 and 11) with the various enhancements discussed above were spatially discretized by finite differences with 100 node points. The resulting set of nonlinear ordinary differential equations was reduced to a large system of nonlinear algebraic equations suitable for incorporation into nonlinear optimization codes through temporal discretization with orthogonal collocation on finite elements. For this purpose, we used Radau collocation, where each pass corresponded to 1 finite element and 2 internal collocation points were employed within each finite element. We determined that the number of node points, finite elements, and collocation provided an acceptable compromise between the solution accuracy and computational efficiency. The parameter estimation problem was posed as a constrained minimization problem with the following objective function and equality constraints corresponding to the discretized model equations and continuity conditions across the finite elements Np

Ψ¼

N

n

∑∑∑

p¼1 i¼1 j¼1

½^np ðp, vj , iÞ  np ðp, vj , iÞ2 np ðp, vj , iÞ2

ð24Þ

where np(p,vj,i) is the measured value of the drop-volume distribution corresponding to drop volume vj and the ith 6094

dx.doi.org/10.1021/ie101818h |Ind. Eng. Chem. Res. 2011, 50, 6089–6100

Industrial & Engineering Chemistry Research

ARTICLE

Figure 5. Critical diameter fits to fifth pass d32 data. (a) Experimental mean diameter d32 versus pressure at 5 wt % oil and 1 wt % surfactant. (b) Experimental mean diameter d32 and linear regression equation versus σ3/5/ɛ2/5Fc1/5 in eq 23 at 5 wt % oil and 1 wt % surfactant. (c) Experimental mean diameter d32 versus pressure at 10 wt % oil and 2 wt % surfactant. (b) Experimental mean diameter d32 and linear regression equation versus σ3/5/ɛ2/5Fc1/5 in eq 23 at 10 wt % oil and 2 wt % surfactant.

homogenizer pass at the pth pressure, n^p(p,vj,i) is corresponding predicted value from the PBE model, n is the total number of node points, N is the number of passes, and Np is the number of pressure sets used. We also optimized both pass-by-pass objective functions over all of the pressure sets Ψi and individual pressure objective functions over all of the passes Ψp. The various optimization problems were formulated in AMPL24 and solved using the nonlinear solver CONOPT. Our optimization method provides point estimates of the parameters without confidence intervals. Values of the estimated parameter are not provided for each case because they provide little insight into the underlying physicochemical behavior. Because of the potential discontinuities in the breakage rate g(v) when accounting for complete surfactant depletion and critical drop diameter effect, this version of the model was not suitable for optimization and parameters were obtained from another optimization run as discussed below.

3. RESULTS AND DISCUSSION 3.1. Base-Case Model. We developed a base-case model without any of the model enhancements listed above to establish a baseline by which to assess the predictive improvements afforded by the various enhancements. The base-case PBE model (eq 1) used the simple energy dissipation rate (eq 5) and the multiple drop breakage function (eq 10) with p = 20 drops to be consistent with our previous work.8 PBE model parameters K1K4 were estimated from drop-volume distribution data

collected for 5 wt % oil, 1 wt % surfactant, five homogenization pressures (250, 500, 750, 1000, and 1250 bar), and five passes at each pressure. The pass-by-pass objective function Ψi, individual pressure objective function Ψp, and total objective function Ψ (eq 24) were used as error measures to assess deviations of the predicted distributions from the experimental data. The values of the parameters are given in Table 2. Modeling results for the base case are shown in Figure 1ac for three of the five pressures studied. Each pressure set includes the experimental data and model prediction for the first, third, and fifth passes as well as the objective function Ψp for that pressure. Relatively large objective function values were obtained at the two extreme pressures because of underprediction of the breakage at lower pressures and overprediction at higher pressures. When summed over all five pressures, a total objective function Ψ of 4.46 was obtained. Figure 1d shows trends in the pass-bypass objective function Ψi with respect to the number of passes. At 250 bar, the prediction error was higher for the initial passes and decreased for later passes, which is consistent with the underprediction observed at this pressure. At higher pressures, the largest errors were observed for intermediate passes. To examine the possible role of unmodeled coalescence in degrading predictions of the breakage-only PBE model, we repeated the analysis using the combined breakagecoalescence PBE model (eq 11 by estimating the model parameters K1K6 from the same drop-volume-distribution data. Because limitations in the problem size that AMPL could handle, we utilized only a single pressure set for each estimation run. For all five 6095

dx.doi.org/10.1021/ie101818h |Ind. Eng. Chem. Res. 2011, 50, 6089–6100

Industrial & Engineering Chemistry Research

ARTICLE

Figure 6. Modeling results when heating loss, variable surface tension, and critical diameter were added to the base-case model. Experimental and predicted drop-volume distributions for the (P) premix, (1) first pass, (3) third pass, and (5) fifth pass. Ψp is the objective function value for the particular pressure set: (a) 250 bar pressure set; (b) 750 bar pressure set; (c) 1250 bar pressure set. (d) Combinations of pressures and pass numbers where the critical diameter limit was reached at 750 and 1250 bar. The limit was not reached at lower pressures.

pressures, the optimizer returned K6 = 0, which corresponds to a zero coalescence frequency. Because this solution was obtained over a wide range of initial guesses, the results suggest that adding coalescence to the PBE model has little value and that any coalescence that occurs can be adequately captured by the breakage-only PBE model. These computational results are consistent with the experimental results from our previous paper,5 in which we demonstrated negligible coalescence at 0.5 wt % oil and 0.1 wt % surfactant, conditions that also yielded an oil-to-surfactant ratio of 5. Therefore, the breakage-only PBE model was used for all further tests. 3.2. Heating Effects. We have observed substantial increases in the emulsion temperature for multiple passes at high pressure. While in this study the emulsion was cooled to room temperature after each pass to minimize the effects of heating on emulsion properties, these results suggest that a considerable amount of mechanical energy may be dissipated as thermal energy and may not be available for drop breakage. Figure 2a shows the temperature rise ΔT averaged over the fives passes as a function of the pressure, where the error bars represent the standard deviation. As expected, the ΔT increase was larger at higher pressures, with 1250 bar producing about a 15 C increase per pass. The average ΔT of the sample and Tins of the homogenizer were used to calculate the energy gained by the sample EH at each pressure (eq 14). Figure 2a also shows the effective pressure actually available for drop breakage calculated for each pressure by eq 15. Generally, the effective pressure was less than 50% of the applied pressure, with the largest percentage differences observed at

higher pressures. Using these effective pressure values in the breakage rate functions (eqs 3 and 4), the constants K1K4 were reestimated to minimize the least-squares error measure (eq 24). Modeling results with the heating effect included through the effective pressure are shown in Figure 2bd for three of the five pressures included in the estimation data set. Each pressure set includes the experimental data and the model prediction for the first, third, and fifth passes as well as the objective function Ψp for that pressure. As compared to the base case, the largest reduction in Ψp was obtained at 250 bar (42%) because the new model parameters reduced the amount of breakage underprediction at low pressures. However, only marginal improvements were observed at higher pressures, resulting in a reduction in the total objective function Ψ to 3.89 (13%). Therefore, inclusion of the mechanical energy loss through thermal heating of the emulsion and homogenizer only had a significant effect on model predictions of the drop-volume distribution at low pressures. 3.3. Variable Interfacial Tension. We recalculated the oil water interfacial tension σ on a pass-by-pass basis to account for depletion of the surfactant available for adsorption on drop interfaces. The interfacial tension σ was measured as a function of the bulk surfactant concentration c (Figure 3a). These data were fit to a quadratic equation to allow calculation of the derivative needed in eq 17. The normalized surface load Γ/Γ¥ (Figure 3b) was used to estimate c1/2. Using this information, an updated σ value was calculated for each pass and used in the breakage rate functions as described in section 2.2. The constants K1K4 were reestimated with a revised model containing both 6096

dx.doi.org/10.1021/ie101818h |Ind. Eng. Chem. Res. 2011, 50, 6089–6100

Industrial & Engineering Chemistry Research

Figure 7. Modeling results when heating loss, variable surface tension, critical diameter, and surfactant depletion were added to the base-case model. The formulation used was 5 wt % oil and 0.1 wt % surfactant, and the pressure used was 250 bar. (a) Pass numbers where the critical diameter limit and the surfactant limit were reached. (b) Experimental and predicted drop-volume distributions for the (P) premix, (1) first pass, (3) third pass, and (5) fifth pass.

the interfacial tension effect and the heating effect to minimize the least-squares error measure (eq 24). The largest increases in the interfacial tension were calculated at high pressures (Figure 4a). Simulation results with the revised model are shown in Figure 4bd for three of the five pressures included in the estimation data set. Compared to the revised model with only the heating effect included, only marginal improvements were obtained by adding the surface tension effect. The maximum improvement in an individual pressure objective function Ψp was 6.5% at 250 bar, and the total objective function Ψ was only decreased by 6% from the heating effect case. Therefore, extension of the model to include increasing surface tension due to depletion of the bulk surfactant did not have a significant effect on predictions of the drop-volume distribution. 3.4. Critical Drop Diameter for Breakage. We utilized eq 23 to account for the existence of a pressure-dependent maximum stable diameter, below which drops were stable to breakage due to mechanical energy limitations of the homogenizer. Parts a and c of Figure 5 show the results of the linear regression used to estimate the proportionality constant c1 for two different formulations. Very similar fifth pass d32 values were obtained for both formulations, with d32 decreasing monotonically with pressure as expected. Reasonably good fits of fifth pass d32 values versus σ3/5/ɛ2/5Fc1/5 values were obtained for both formulations, supporting the assumed linear relationship in eq 23.

ARTICLE

Moreover, the two formulations yielded virtually identical values of the proportionality constant, suggesting that only a single data set was needed for estimation and that c1 = 0.37 could be used over a range of oil and surfactant concentrations. The breakage rate constants K1K4 obtained for the model with the heating and interfacial tension effects only were used in this case because the critical diameter limit was active only for certain conditions and reestimation proved difficult because of the discontinuity introduced by setting g(v) = 0 for d32 < dmax. Figure 6d shows the combinations of pressures and pass numbers for which the ratio d32/dmax dropped below zero, and breakage was stopped for the revised model that included the heating effect, the interfacial tension effect, and the critical diameter limit. No results are shown for lower pressures because the critical diameter limit was never active. The limit had the greatest impact at 1250 bar, where breakage was stopped at a finite element between the second and third passes. At 750 bar, breakage was stopped between the third and fourth passes. Predicted drop-size distributions obtained with the revised model are shown in Figure 6ac. The critical diameter limit had no effect at 250 bar, so the individual pressure objective function Ψp was identical with that obtained for the model with just the heating and interfacial tension effects included. The Ψp value was reduced by 67% at the 1250 bar because overprediction was avoided, while the limit offered no advantage at 750 bar. The total objective function Ψ was decreased by 22% compared to the revised model without the critical diameter limit. Therefore, we concluded that inclusion of the critical limit was advantageous primarily at high pressures, where the breakage rate would be overpredicted otherwise. 3.5. Complete Surfactant Depletion. We included a surfactant balance (eqs 21 and 22) in the model and set the breakage rate g(v) = 0 when the balances predicted that no free surfactant remained for stabilization of newly formed drops. The breakage rate constants K1K4 obtained for the model with the heating and interfacial tension effects only were used because reestimation proved difficult because of the discontinuity introduced by setting g(v) = 0 for c e 0. The revised model used for simulation included the heating effect, variable interfacial tension, critical breakage diameter, and surfactant depletion effect. For the three formulations listed in Table 1, we found that complete surfactant depletion only occurred when the critical diameter limit had already been reached because of the relatively large oil-to-surfactant ratio. To explore conditions under which complete surfactant depletion might have occurred independent of the critical diameter limit, we performed experiments at a higher oil-to-surfactant ratio (5 wt % oil and 0.1 wt % surfactant) and a lower pressure (250 bar pressure). This combination produced conditions under which complete surfactant depletion occurred at a finite element between the third and fourth passes, while the critical diameter limit was never reached (Figure 7b). Predicted drop-size distributions obtained with the revised model are shown in Figure 7a. Compared to an equivalent model with the surfactant depletion effect eliminated (not shown), individual pressure objective function Ψp was increased slightly, although the predictions were not degraded significantly. On the basis of these results, we concluded that inclusion of complete surfactant depletion could be important at low pressures and relatively high oil-to-surfactant ratios that might also promote coalescence, which was unmodeled in this study. 3.6. Effect of the Number of Daughter Drops. We have shown previously that the inclusion of multiple drop breakages in 6097

dx.doi.org/10.1021/ie101818h |Ind. Eng. Chem. Res. 2011, 50, 6089–6100

Industrial & Engineering Chemistry Research

ARTICLE

Figure 8. Modeling results when the number of daughter drops formed p was optimized and heating loss, variable surface tension, critical diameter, and surfactant depletion were added to the base-case model. Experimental and predicted drop-volume distributions for the (P) premix, (1) first pass, (3) third pass, and (5) fifth pass. Ψp is the objective function value for the particular pressure set. (a) Total objective function Ψ as a function of the different numbers of daughter drops formed upon breakage: (b) 250 bar pressure set for p = 150; (c) 750 bar pressure set for p = 150; (d) 1250 bar pressure set for p = 150.

the PBE model generated more accurate drop-volume distribution predictions than an equivalent model with binary drop breakage.8 In that study, we assumed that 20 daughter drops were formed upon breakage of a mother drop. Here we performed optimization of the number of daughter drops formed in an effort to improve pressure predictions. The number of daughter drops p was varied between 20 and 500, and the breakage rate constants K1K4 were reestimated for each case. The model used for this purpose was the base-case model and included none of the proposed effects. Figure 8a shows the total objective function Ψ as a function of p. A minimum was observed at p = 150 drops, with both smaller and larger p values producing larger Ψ. While experimental studies on the number of daughter drops formed under turbulent flow are lacking, our results were consistent with other modeling studies in which the formation of about 100 drops have been reported.3,11 Predicted drop-volume distributions obtained from the revised model with p = 150 and all four other effects included are shown in Figure 8bd for the optimized breakage parameters obtained with this p value. Compared to the equivalent model with p = 20, the individual pressure objective function Ψp was reduced by 40% at 250 bar and 57% at 750 bar, while no significant change was observed at 1250 bar. The total objective function Ψ was decreased by 62%, suggesting that including an appropriate number of daughter drops formed is important for an accurate prediction of the impact of pressure on the drop-volume distribution. 3.7. Extensibility of the Revised Model to New Formulations. The primary objective of this study was to determine the

Figure 9. Total objective function values obtained with various model enhancements for three emulsion formulations. The breakage rate parameters were obtained from data collected for the base-case emulsion of 5 wt % oil and 1 wt % surfactant. The extensibility of the model to the other two formulations listed in Table 1 was evaluated. The model enhancements were implemented sequentially such that “heating loss” represents the base case with inclusion of the heating effect, “variable interfacial tension” represents the base case with inclusion of the heating and variable interfacial tension effects, etc.

enhancements necessary to improve the pressure predictions of our previously developed PBE model.8 To be useful for homogenizer simulation and ultimately emulsified product design, the 6098

dx.doi.org/10.1021/ie101818h |Ind. Eng. Chem. Res. 2011, 50, 6089–6100

Industrial & Engineering Chemistry Research

ARTICLE

Table 3. Total Objective Function Ψ for 5 wt % Oil and 1 wt % Surfactant effect

objective function

% improvement

base case

4.46

heating interfacial tension

3.89 3.68

13 17

critical diameter

2.85

36

number of daughter drops

1.69

62

model must be extensible to new formulations not included in the data set used for parameter estimation. As a modest test of this capability, the revised model, which included all of the aforementioned enhancements, was developed at 5 wt % oil and 1 wt % surfactant as described above and tested for two other formulations with the same oil-to-surfactant ratio (Table 1). Figure 9 summarizes the improvements achieved with the various enhancements, as measured by the total objective function Ψ for the three formulations. Because it offered no advantage for these formulations, the complete surfactant depletion effect was not considered in this analysis. As before, the effects were implemented sequentially such that “heating loss” represents the base case with inclusion of the heating effect, “variable interfacial tension” represents the base case with inclusion of the heating and variable interfacial tension effects, etc. At the basecase conditions of 5 wt % oil and 1 wt % surfactant, simultaneous implementation of all enhancements produced a 62% improvement over the base-case model (Table 3). For all three formulations, the most significant improvements were obtained by increasing the number of daughter drops formed upon breakage from 20 to 150 drops and by introducing a critical diameter, below which drops could not break. The other two effects considered here, namely, the introduction terms that describe heat loss due to thermal heating of the sample and homogenizer and extension of the model to account for the effects of surfactant adsorption on the interfacial tension, offered only marginal improvements over the base-case model. Therefore, we concluded that the specification of an appropriate number of daughter drops and the inclusion of a pressure-dependent critical diameter for breakage were important for improving the ability of the PBE model to generate accurate drop-volume distribution predictions over a wide range of pressures.

’ AUTHOR INFORMATION Corresponding Author

*Phone: 413-545-3481. Fax: 413-545-1647. E-mail: henson@ ecs.umass.edu.

’ ACKNOWLEDGMENT We acknowledge financial support from the National Science Foundation (Grants CBET-0730795 and DMI-0531171) and Unilever Foods. ’ NOTATION a parameter for solute type c aqueous-phase surfactant concentration, kg m3 c1/2 half-life aqueous-phase surfactant concentration, kg m3 cads adsorbed surfactant concentration at the drop interface, kg m3

C(v,v0 ) Cp Cs Di Do d dmax EH g(v) g1(v) g2(v) G Gb h(v,v0 ) hgap K1, K2, K3, K4 K5 K6 K1, K2, K3, K4 n n(v,t) np(v,t) ^np(p,vj,i) N Np ΔP p Q q R r Τ ΔΤ t tB v Vdiss Vtot

coalescence frequency of drops of size v and v0 , s1 specific heat capacity, J kg1 K1 surfactant concentration, kg m3 inlet diameter of the valve, m exit diameter of the valve, m drop diameter, m maximum stable drop diameter, m thermal energy, N m2 breakage rate of drops of size v, s1 turbulent inertial breakage rate of drops of size v, s1 turbulent viscous breakage rate of drops of size v, s1 shear rate, s1 critical shear rate, s1 collision frequency of drops of size v and v0 , s1 valve gap, m adjustable parameters in breakage rate function g(v) adjustable parameters in the collision frequency function h(v,v0 ) adjustable parameters in the coalescence effciency function λ(v,v0 ) adjustable parameters in the breakage rate function g(v) number of node points number of drops with size v per unit volume of the dispersion at time t volume percent distribution of drops of size v at time t predicted volume percent distribution at pressure p for drops of size vj at the ith pass number of passes number of pressures applied pressure, Pa parameter of eq 10 that represents the number of daughter drops volumetric flow rate, m3 s1 parameter of eq 10 gas constant, J K1 mol1 parameter of eq 10 absolute temperature, K temperature rise, K time, s breakage time, s volume of drops, m3 dissipation volume, m3 total volume of the drops, m3

Greek Symbols

β(v,v0 ) ɛ φ ΨT ηc ηd Fd Fc σ Γ Γ¥ 6099

probability density function for daughter drops local energy dissipation per unit volume, kg m1 s3 volume fraction of the dispersed phase total objective function defined in eq 24 viscosity of the continuous phase, Pa s viscosity of the dispersed phase, Pa s density of the dispersed phase, kg m3 density of the continuous phase, kg m3 interfacial tension, N m1 surface load of the surfactant, kg m2 plateau value of the surface load of the surfactant, kg m2 dx.doi.org/10.1021/ie101818h |Ind. Eng. Chem. Res. 2011, 50, 6089–6100

Industrial & Engineering Chemistry Research λ λ(v,v0 ) ν(v0 ) Π

ARTICLE

viscosity ratio coalescence efficiency of drops of size v and v0 number of daughter drops formed by breakage of a drop of size v0 surface pressure, N m1

Subscripts

c d ins o/w oil water

continuous phase dispersed phase instrument temperature rise oilwater oil water

’ REFERENCES (1) Becher, P. Emulsions: Theory and Practice; Oxford University Press: New York, 2001. (2) McClements, D. J. Food Emulsions: Principles, Practice, and Techniques; CRC Press: Boca Raton, FL, 2005. (3) Vankova, N.; Tcholakova, S.; Denkov, N.; Vulchev, V.; Danner, T. J. Colloid Interface Sci. 2007, 313 (2), 612–629. (4) Soon, S.; Harbidge, J.; Titchener-Hooker, N.; Shamlou, P. J. Chem. Eng. Jpn. 2001, 34 (5), 640–646. (5) Raikar, N.; Bhatia, S.; Malone, M.; Henson, M. Chem. Eng. Sci. 2009, 64 (10), 2433–2447. (6) Ramkrishna, D. Population Balances: Theory and Applications to Particulate Processes in Engineering; Academic Press: New York, 2000. (7) Hakansson, A.; Tragardh, C.; Bergenstahl, B. Chem. Eng. Sci. 2009, 64 (12), 2915–2925. (8) Raikar, N.; Bhatia, S.; Malone, M.; Almeida-Rivera, C.; Bongers, P.; McClements, D.; Henson, M. Colloids Surf., A 2010, 361, 1. (9) Cortes-Munoz, M.; Chevalier-Lucia, D.; Dumay, E. Food Hydrocolloids 2009, 23 (3), 640–654. (10) Tcholakova, S.; Vankova, N.; Denkov, N.; Danner, T. J. Colloid Interface Sci. 2007, 310 (2), 570–589. (11) Wieringa, J.; Vandieren, F.; Janssen, J.; Agterof, W. Chem. Eng. Res. Des. 1996, 74, 554–562. (12) Rallison, J. Annu. Rev. Fluid Mech. 1984, 16, 45–66. (13) Stone, H. Annu. Rev. Fluid Mech. 1994, 26, 65–102. (14) Zhao, X.; Goveas, J. Langmuir 2001, 17 (13), 37883791. (15) Vankova, N.; Tcholakova, S.; Denkov, N.; Ivanov, I.; Vulchev, V.; Danner, T. J. Colloid Interface Sci. 2007, 312 (2), 363–380. (16) Canu, P. Ind. Eng. Chem. Res. 2005, 44 (8), 2649–2658. (17) Chen, Z.; Pruss, J.; Warnecke, H. Chem. Eng. Sci. 1998, 53 (5), 1059–1066. (18) Coulaloglou, C.; Tavlarides, L. Chem. Eng. Sci. 1977, 32 (11), 1289–1297. (19) Phipps, L. J. Phys. D: Appl. Phys. 1975, 8 (4), 448–462. (20) Diemer, R.; Olson, J. Chem. Eng. Sci. 2002, 57 (19), 4187–4198. (21) Hill, P.; Ng, K. AIChE J. 1996, 42 (6), 1600–1611. (22) Zaccone, A.; Gabler, A.; Maass, S.; Marchisio, D.; Kraume, M. Chem. Eng. Sci. 2007, 62 (22), 6297–6307. (23) Kumar, S.; Ramkrishna, D. Chem. Eng. Sci. 1996, 51 (8), 1311– 1332. (24) Fourer, R.; Gay, D. M.; Kernighan, B. W. AMPL: A Modeling Language for Mathematical Programming; Brooks/Cole Publishing Company: Pacific Grove, CA, 2003.

6100

dx.doi.org/10.1021/ie101818h |Ind. Eng. Chem. Res. 2011, 50, 6089–6100