Achieving Target Emulsion Drop Size Distributions Using Population

Article pubs.acs.org/IECR

Achieving Target Emulsion Drop Size Distributions Using Population Balance Equation Models of High-Pressure Homogenization Shashank N. Maindarkar,† Hans Hoogland,‡ and Michael A. Henson*,† †

Department of Chemical Engineering, University of Massachusetts, Amherst, Massachusetts 01003-9303, United States Unilever R&D, 3133 AT Vlaardingen, The Netherlands

‡

S Supporting Information *

ABSTRACT: Population balance equation (PBE) models have been used extensively to predict drop size distributions (DSDs) of dispersed phase systems. In our previous publications, we have used the PBE framework to develop increasingly sophisticated process models for oil-in-water emulsification in high-pressure homogenizers. The goal of this study was to utilize these PBE models for integrated emulsion product and process design through the formulation and solution of homogenizer optimization problems. For a specified number of homogenization passes, the initial amount of surfactant and the pressure of each pass were determined by solving a nonlinear least-squares optimization problem such that the target DSD were achieved. Three alternative objective functions that differed with respect to the distribution specification and the penalty on surfactant usage were formulated and solved for different target DSDs. The model predictions were successfully validated by performing homogenization experiments using the optimized formulation and homogenization variables.

1. INTRODUCTION An emulsion is a mixture of two or more immiscible liquids to form a dispersed phase system. Oil-in-water emulsions are very commonly encountered in consumer products, processed foods, and pharmaceutical formulations.1−4 Emulsion product properties have a complex dependence on the chemical formulation and the processing conditions. The emulsion formulation influences dispersed and continuous phase properties such as density, viscosity, and dielectric constant, while the surfactant used plays a critical role in determining interfacial properties as well as emulsion stability. Processing conditions such as mixing time, homogenization pressure, or rotor speed also have a substantial effect on emulsified product properties, including appearance, taste, mouth feel, odor, and safety. All these variables impact the drop size distribution (DSD), a key property that influences emulsion rheology, stability, texture, and appearance. A typical emulsified product requires the DSD to be maintained within acceptable limits, which includes achieving a prescribed mean drop size and maintaining acceptable variations about the mean. For most pharmaceutical and cosmetic products, the mean drop size is required to be less than 1 μm to ensure biocompatibility and stability of the product.5 Oil-in-water emulsions are generally prepared in a multistep process that depends on the application, and the choice of emulsification devices is based on the energy intensities required.4 In first step, a coarse premix is prepared using a low shear stator-rotor type device that mixes the various ingredients into a stable form. This coarse premix is then processed with an energy intensive device, such as a highpressure homogenizer or a microfluidizer. The emulsion often must be passed through the device multiple times to achieve a satisfactory DSD. Due to their scalability and relatively low cost, high-pressure homogenizers are the most common high-energy devices used in industry. A coarse emulsion is pumped through © XXXX American Chemical Society

a small orifice under very high pressure. The emulsion passes radially through the narrow gap between the piston and the valve seat at high velocity, creating a local environment of high turbulence that causes drop deformation and breakage.6 Drop breakage leads to the creation of new interfacial area that must be stabilized by the surfactant. If newly formed drops are not sufficiently covered by surfactant, the drops will coalesce to form larger drops. The inability of the surfactant to achieve adequate drop coverage can be due to either insufficient free surfactant in solution or slow adsorption kinetics.7,8 Achieving a desired DSD is a very challenging problem as the effects of the emulsion formulation (e.g., oil and surfactant concentrations) and processing conditions (e.g., the homogenization pressure and number of passes) on the DSD is complex and incompletely understood.8 Due to lack of quantitative understanding, new emulsified products are usually developed by combining a broad knowledge of previous product formulations with empirical scientific experimentation. An alternative is to use a suitable mathematical model to predict the DSD for different emulsion formulations and processing conditions. The population balance equation (PBE) modeling framework 9 is particularly well suited for this problem, as mechanistic functions describing drop breakage and coalescence can be incorporated within a fundamental number balance equation to predict the evolution of the DSD. Many alternative PBE models of high-pressure homogenizers have been developed.10−18 As part of our previous work,19,20 we have developed PBE models with mechanistic functions for Special Issue: Doraiswami Ramkrishna Festschrift Received: March 30, 2015 Revised: June 28, 2015 Accepted: August 10, 2015

A

DOI: 10.1021/acs.iecr.5b01195 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research drop breakage and drop coalescence to predict the DSD as a function of the homogenizer pass. We also have incorporated a calculated emulsion viscosity into the drop breakage and coalescence functions and developed a dynamic surface coverage model that produced good DSD predictions over range of surfactant and oil concentrations.21 Despite the availability of predictive models, the development of model-based design strategies for emulsified products has not been investigated. Such techniques have been proposed for other particulate systems including methods for achieving target polymer particle sizes in emulsion polymerization22 and target crystal size distributions in crystallization processes.23−25 In emulsified product design, the DSD can be adjusted by manipulating the emulsion formulation and/or the processing conditions. Generally the concentration of the dispersed phase (e.g., oil) is fixed based on the initial product design, but the amount of surfactant used can be varied. Due to the high cost of surfactant, the goal is to achieve the target DSD with the minimum amount of surfactant. In this work, we utilized our previously developed PBE model 21 to develop and validate model-based design strategies for achieving target DSDs in high-pressure homogenizers. We formulated three alternative objective functions to be minimized: (1) the least-squares error between the target and predicted mean sizes and polydispersities; (2) the least-squares error between the target and predicted DSDs; and (3) the total amount of surfactant used subject to the constraint that the least-squares error between the target and predicted DSDs was below a specified bound. For a given number of homogenization passes, the decision variables were the initial amount of surfactant added and the pressure of each pass. Selected modelbased designs were validated by performing homogenization experiments at the optimal conditions and comparing the measured and target DSDs.

always same regardless of extra surfactant solution added in the second step. The coarse emulsion obtained from the two step procedure was passed through a high-pressure homogenizer (Emulsiflex C-3, Avestin Inc.) five times to reduce the drop size. For model validation, experiments were performed at six different homogenization pressures (100, 200, 300, 400, 600, and 800 bar). After each homogenization pass, approximately 2 mL of emulsion was sampled, and the DSD was measured using static light scattering (Mastersizer 2000, Malvern Instruments). For validation of model-based designs, experiments were performed with the optimal surfactant concentrations and homogenization pressures predicted by the model. 2.3. Dynamic Simulation and Parameter Estimation. The integral expressions in the PBE (see Appendix, eq A.1) were numerically approximated as summations using the fixed pivot technique 26 with the drop volume discretized into 100 node points equally spaced in a logarithmic scale. Spatial discretization produced 100 nonlinear ordinary differential equations (ODEs) in time with the volume percent distribution at each node point as the dependent variables. These 100 ODEs along with two additional ODEs for the free surfactant concentration (eq A.9) and the surface coverage (eq A.8) were integrated with the Matlab code ode15s. The premix distribution N(v,0) was used as the initial condition for the PBE, while eqs A.10 and A.11 were used to generate initial conditions for the other two ODEs. A homogenization pass was defined to be one dimensionless time unit. Bulk properties of the emulsion system used in the model equations are listed in Table S1. The homogenizer model had five adjustable parameters, K1− K5, in the breakage and coalescence functions. These parameters were estimated by nonlinear optimization using measured drop volume distributions collected for five homogenization passes at the base case conditions of 50 wt% oil, 1 wt% surfactant, and 400 bar pressure following the procedure detailed in our previous publications.19−21 The ODE system obtained from spatial discretization was temporally discretized using orthogonal collocation with 15 finite elements and two internal collocation points per element to obtain a set of nonlinear algebraic equations. The algebraic equations were posed as equality constraints, and continuity conditions were imposed across the finite elements. The parameter estimates were determined by minimizing the following objective function,

2. MATERIALS AND METHODS 2.1. Materials. Oil-in-water emulsions were prepared using 50 wt% sunflower oil (Sigma) as the dispersed phase, 1 wt% Pluronic F68 (PF68, Sigma) as the surfactant, and nanopure water as the continuous phase. For model-based design validation, the optimal amount of surfactant (0.5−2.5 wt%) predicted by the model was used. 2.2. Emulsion Preparation and Characterization. Emulsions were prepared in two steps. In the first step, a coarse pre-emulsion was prepared by mixing the ingredients in a rotor-stator device (Ultra-Turrax Model T25, Rose Scientific Ltd.) at 13000 rpm for 20 min. In the second step, the coarse emulsion was passed through the high-pressure homogenizer. To ensure that the same premix DSD was obtained for any x wt % of surfactant used (0.5 ≤ x ≤ 2.5 wt%), the coarse emulsion with 50 g of oil, x g of PF68, and (50 − x) g of water was prepared in two steps. In the first step, 0.5 g of PF68 was dissolved in 39.5 g of water and mixed with 50 g of oil using the Ultra-Turrax blender to form the first coarse emulsion. In the second step, another surfactant solution was prepared with (x − 0.5) g of PF68 in (10 − x + 0.5) g of water. This second surfactant solution was then added to the first coarse emulsion and mixed for a minute using a magnetic stirrer at low stirring speed of 200 RPM such that the DSD of the first coarse emulsion was not altered. Since the coarse emulsion prepared in first step has fixed amount of oil (50 g), surfactant (0.5 g), and water (39.5 g), the DSD of the final coarse emulsion is

Np

Ψ=

∑ i=1

100

∑ j = 1 [N (vj , i) − Ne(vj , i)]2 100

∑ j = 1 [Ne(vj , i)]2

(1)

where Ne(vj,i) is the measured value of the drop volume distribution at drop volume vj and homogenizer pass i, N(vj,i) is the corresponding predicted value obtained from the discretized model, and Np is the number of passes. This objective function ensured that the denominator would remain nonzero even if the number of drops in a particular volume class became zero. The parameter estimation problem was formulated in AMPL (AMPL Optimization Inc., Albuquerque, NM, USA)27 and solved with the nonlinear programming code CONOPT (ARKI Consulting & Development A/S, Bagsvaerd, Denmark).28 B


Article

Industrial & Engineering Chemistry Research

Figure 1. Empirical functions for the dependencies of the (a) parameter K1 and (b) the energy dissipation rate ε on the homogenization pressure.

Figure 2. Drop volume distributions (blue, first pass; green, second pass; red, third pass; cyan, fourth pass; violet, fifth pass) at homogenization pressures of (a) 200 bar (Ψ = 0.0585), (b) 400 bar (Ψ = 0.0310), (c) 600 bar (Ψ = 0.0323), and (d) 800 bar (Ψ = 0.0630) obtained using parameters K2−K5 estimated at 400 bar pressure and fitted equations for the parameter K1 and the energy dissipation rate ε.

3. RESULTS AND DISCUSSION

concentrations and oil fractions using a single set of optimized parameters. However, model extensibility for the homogenization pressure was not investigated in that study. Therefore, we performed homogenization experiments at six different pressures (100−800 bar) using fixed amounts of oil and surfactant (50 wt% oil, 1 wt% PF68). Parameter estimation was performed with DSDs measured at 400 bar pressure, and the extensibility of the model with these optimized parameters was

3.1. Model Extensibility for Homogenization Pressure. To be successfully used for emulsified product design, the homogenizer model must be extensible over a wide range of emulsion formulations and processing conditions. In our previous work,29 we demonstrated that the model could generate satisfactory predictions over a range of surfactant C


Article

Industrial & Engineering Chemistry Research investigated at the other five pressures. The predicted DSDs deviated substantially from the measured DSDs as indicated by the relatively large Ψ values obtained at the five pressures (Table S2). Based on available computational fluid dynamics results,30 we hypothesized that the dependence of the energy dissipation rate on homogenization pressure was more complex than the simple linear form used (eq A.7). Furthermore, we anticipated that the pre-exponential factor in the breakage frequency function g1 (eq A.2) should depend on the energy dissipation rate ε such as the breakage frequency function g2 (eq A.3). To address these model limitations, we re-estimated the parameters K1 and ε at each pressure while keeping the parameters K2−K5 fixed at their previous values. We found that the dependencies of K1 and ε on pressure could be accurately fit with a power law equation and a quadratic equation, respectively (Figure 1). These empirical equations allowed the model to generate accurate DSD predictions over the entire range of pressures with a single set of parameters (Figure 2, Table S2). 3.2. Model-Based Emulsified Product Design Methodology. The goal of emulsified product design was to obtain the desired DSD properties while minimizing the usage of surfactant. Three alternative objective functions were defined to achieve this goal. The first objective function Ψ1 was defined as the normalized root-mean-square error between the target and model-predicted Sauter mean diameters d32 and polydispersities PD, Ψ1 =

(d32,tar − d32)2 2 d32,tar

+

W penalizes the surfactant term and therefore provides a tradeoff between the applied pressures and surfactant usage. To solve each problem, the homogenizer model was spatially and temporally discretized to generate a nonlinear algebraic equation system that was posed as a set of equality constraints. The decision variables of the optimization problems were the (i) the initial (total) surfactant concentration and (ii) the homogenization pressure at each pass. As number of passes is a discrete rather than continuous variable, we solved the optimization problems for different numbers of passes to avoid a considerably more complex mixed-integer nonlinear programming problem. Given a specified objective and a fixed number of homogenization passes, the problem of determining the optimal amount of surfactant and the optimal pressure of each pass was solved within AMPL using the nonlinear optimization code CONOPT. Three different design problems were formulated with different objective functions (see Table 1). The pressure of each pass was constrained to be between Table 1. Overall Objective Functions for Different Design Problems

achieve the target d32 and PD while minimizing the surfactant used achieve the target DSD while minimizing the surfactant used minimize the surfactant used with the DSD error maintained below a critical value

(PDtar − PD)2 PD2tar

(2)

Ψ2 =

Ψ2 + WΨ3 Ψ3, subject to Ψ2 < 0.01

Table 2. Decision Variables variable

dimensions

units

initial guess

lower bound

upper bound

Cs,in P

scalar vector of length equal to the number of passes

wt% bar

1 400

0.5 100

2.5 800

(3)

3.3. Problem 1: Achieve Target Mean Diameter and Polydispersity. First we considered the problem of achieving the target mean drop size (d32,tar) and the target polydispersity (PDtar) using the least amount of surfactant. Two different targets were considered: (1) d32,tar = 0.75 μm, PDtar = 2, and (2) d32,tar = 0.5 μm, PDtar = 1.5. For 1−5 homogenization passes, the optimization problem was solved to find the optimal surfactant amount and the optimal pressure of each pass by minimizing the objective function Ψ1 + WΨ3 with the weighting factor W = 0.1. The model-predicted mean drop sizes, polydispersities, and objective function values are listed in Table 3 for the first set of targets. As the number of passes was increased, the predicted d32 and PD approached their target values. The objective function decreased rapidly for the first two passes, with additional passes producing small improvements. The solutions for one and two passes were characterized by pressures at the upper bound of 800 bar, suggesting that additional passes were necessary to achieve the desired DSD properties. The amount of surfactant used decreased with the number of passes such that four passes required 15.6% less surfactant than a single pass.

where ntar is the target DSD and n is the model-predicted DSD. To minimize the amount of surfactant used, a third objective function Ψ3 was formulated as follows, ⎛ Cs,in − 0.5 ⎞2 Ψ3 = ⎜ ⎟ ⎝ ⎠ 0.5

Ψ1 + WΨ3

100 and 800 bar, and the amount of surfactant used to prepare the premix emulsion was constrained to be between 0.5 and 2.5 wt% (Table 2).

where the subscript “tar” represents target values and the other variables represent model-predicted values. The polydispersity PD was calculated from measured or predicted DSDs as ratio of d43 to d32. While these target values are relatively simple to specify for a particular problem, their attainment does not ensure that a satisfactory DSD will be achieved. Therefore, the second objective Ψ2 involves the full distribution: n ∑ j = 1 [ntar(vj) − n(vj)]2 n ∑ j = 1 [ntar(vj)]2

objective function (Ψ)

problem

(4)

where Cs,in is the total amount of surfactant used in wt% and the lower bound of 0.5 wt% was specified to ensure that the homogenized emulsions would be stable. Three alternative optimization problems were formulated by combining the individual objective functions as shown in Table 2, where W is an adjustable weighting factor. The first problem involved attainment of the target d32 and PD while minimizing the amount of surfactant used. The second problem was formulated similarly, except that the target emulsion was characterized by the entire DSD. In the third problem, the amount of surfactant used was minimized subject to the constraint that the DSD squared error Ψ2 was maintained below a critical value known from previous experience to produce satisfactory attainment of the target DSD. The weight D


Article


Table 3. Model-Predicted and Experimental Results for the Target d32 = 0.75 μm and PD = 2, Using a Weighting Factor W = 0.1 pressure (bar) total surfactant (wt%)

first pass

second pass

third pass

fourth pass

0.716 0.651 0.634 0.604

800 800 800 434

− 800 485 100

− − 100 194

Model Prediction − − − 217

0.604

434

100

194

Experimental Result 217

d32 (μm)

PD

Ψ1

Ψ3

Ψ

1.353 0.767 0.745 0.750

3.129 2.103 2.005 1.986

0.982 0.056 0.0077 0.0071

0.187 0.091 0.072 0.043

1.001 0.065 0.015 0.011

0.811

2.620

0.321

0.043

0.325

Table 4. Model-Predicted and Experimental Results for the Target d32 = 0.5 μm and PD = 1.5, Using a Weighting Factor W = 0.1 pressure (bar) total surfactant (wt%)

first pass

second pass

third pass

fourth pass

fifth pass

0.790 0.803 0.903 0.921 0.908

800 800 800 800 795

− 800 800 800 800

− − 800 800 800

Model Prediction − − − − − − 758 − 799 634

0.921 0.908

800 795

800 800

800 800

Experimental Result 758 − 799 634

d32 (μm)

PD

Ψ1

Ψ3

Ψ

1.337 0.726 0.558 0.506 0.501

3.146 2.124 1.680 1.539 1.502

2.003 0.614 0.163 0.029 0.0024

0.336 0.367 0.650 0.708 0.665

2.036 0.651 0.228 0.100 0.069

0.689 0.667

1.365 1.347

0.389 0.349

0.708 0.665

0.460 0.416

Table 5. Model-Predicted and Experimental Results for a Target DSD Generated from a Log-Normal Distribution with μ = 0.6 μm and σ2 = 0.5, Using a Weighting Factor of W = 0.01 pressure (bar) total surfactant (wt%)

first pass

second pass

third pass

fourth pass

1.464 1.352 1.050 0.977 0.963

762 796 790 605 100

− 634 635 475 296

Model Prediction − − − − 625 − 363 473 350 383

0.977

605

475

Experimental Result 363 473

fifth pass

Ψ2

Ψ3

Ψ2 + WΨ3

− − − − 526

0.577 0.092 0.005 0.0016 0.0024

3.720 2.904 1.210 0.909 0.857

0.614 0.121 0.017 0.0107 0.0109

−

0.0127

0.909

0.0218

3.4. Problem 2: Achieve Target Drop Size Distribution. Direct specification of a target DSD provides additional flexibility beyond that possible by specifying just a few DSD moments (e.g., mean, polydispersity). However, specification of a meaningful target DSD is problem specific and requires careful consideration. Our experience is that Gaussian-like DSDs are approximately achievable when the drop size is represented in logarithmic scale. Therefore, we generated three different target DSDs from the log-normal distribution by specifying the mean μ = 0.5, 0.6, and 1.0 μm and the variance σ2 = 0.5. For 1−5 homogenization passes, the objective function Ψ2 + WΨ3 with W = 0.01 was minimized to determine the optimal surfactant amount and the optimal pressure at each pass. First we considered a target log-normal distribution with μ = 0.6 μm and σ2 = 0.5. The objective function decreased rapidly for 1−3 passes and appeared to reach a minimum value at four passes (Table 5). Because the objective function weight W = 0.01 placed increasingly more emphasis on minimizing

To validate these model predictions, we performed homogenization experiments with the optimal surfactant amount and the optimal pressures calculated for four passes. While the target values were approximately achieved, the model overpredicted drop breakage such that the actual d32 was slightly above the target value (last row of Table 3). A larger error was observed in the PD value, suggesting that the model underpredicted variability in the drop size. Similar results were obtained for the second set of targets (Table 4). In this case, large improvements in the predicted d32 and PD were achieved for the first four passes. With the exception of five passes, the amount of surfactant used increased with the number of passes, presumably due to the large amount of surface area created when very small drops were produced. We performed homogenization experiments for four and five passes to validate these model predictions. As before, the model overpredicted drop breakage such that the actual d32 values was larger than the target value. Interestingly, the actual PD values were slightly smaller than the target value. E


Article

Industrial & Engineering Chemistry Research surfactant usage once the DSD error Ψ2 was small, Ψ2 showed an increasing trend after four passes. Compared to a single pass, four passes produced a 99.7% decrease in the DSD error and a 33% decrease in surfactant usage. While five passes produced a slight decrease in surfactant usage, the overall objective function increased because the pressure was constrained by the lower bound for the first pass. Experimental implementation of the optimal solution for four passes produced good agreement with the target DSD as only small errors were observed near the tails of the distribution (Figure 3).

Figure 4. Model-predicted and experimental distributions for a target DSD generated from a log-normal distribution with μ = 1.0 μm and σ2 = 0.5 using a weighting factor of W = 0.01.

optimal solution for four passes showed larger deviations from the predicted and target DSDs (Figure 5). The results suggest that the model overpredicts breakage and/or underpredicts coalescence of drops in the 0.1 μm size range.


For the target log-normal distribution with a larger mean μ = 1.0 μm and the same variance σ2 = 0.5, the objective function decreased slowly after two passes because Ψ2 was already small and the surfactant concentration approached the minimum value of 0.5 wt% (Table 6). Compared to a single pass, three passes produced a 48% reduction in surfactant usage. When implemented experimentally, the three pass solution showed only slight deviations from the predicted DSD at small and large drop sizes and good agreement with the target DSD (Figure 4). A target log-normal distribution with a smaller mean μ = 0.5 μm and the same variance σ2 = 0.5 required more passes to generate a small objective function, with negligible improvement beyond four passes (Table S3). The optimal solution for four passes was dominated by surfactant usage since the DSD error Ψ2 was reduced to a very small value and a minimum amount of surfactant was required to stabilize the large interfacial area created. Experimental implementation of the


3.5. Problem 3: Minimize Surfactant Usage. The third formulation of the emulsified product design problem involved minimization of surfactant usage via the objective function Ψ3 subject to the constant the DSD error Ψ2 was below some threshold (see Table 2). This problem formulation offered the

Table 6. Model-Predicted and Experimental Results for a Target DSD Generated from a Log-Normal Distribution with μ = 1.0 μm and σ2 = 0.5, Using a Weighting Factor of W = 0.01 pressure (bar) total surfactant (wt%)

first pass

second pass

1.128 0.631 0.544 0.534

800 796 724 315

− 800 661 270

Model Prediction − − 580 390

0.544

724

661

Experimental Result 580

third pass

F

fourth pass

Ψ2

Ψ3

Ψ2 + WΨ3

− − − 470

0.260 0.011 0.0073 0.0072

1.579 0.068 0.0079 0.0045

0.276 0.011 0.0073 0.0072

−

0.0096

0.0079

0.0097


Article


effect of W on the tradeoff between attainment of the target DSD and minimization of surfactant usage, we solved second problem where the objective function was Ψ2 + WΨ3 over a range of W values. More specifically, W was varied in the range 0.001−0.1, and the optimal design problem was solved for four passes with a target DSD generated from a log-normal distribution with μ = 0.5 μm and σ2 = 0.5. Small W values emphasized DSD attainment at the expense of relatively high surfactant usage (Table S4, Figure 6). On the other hand, large W values substantially reduced surfactant usage but resulted in large DSD errors. While the appropriate balance between these two competing objectives is subjective and problem specific, we preferred W values that produced Ψ2 ≤ 0.01 such that a high degree of DSD attainment was achieved. According to this criterion, the previously used value W = 0.01 was most appropriate for our design problem.

advantage that surfactant costs could be minimized directly, but the inequality constraint Ψ2 < 0.01 introduced the possibility of problem infeasibility. The design problem was solved with two target DSDs generated from a log-normal distribution: (1) μ = 0.6 μm, σ2 = 0.5, and (2) μ = 0.5 μm, σ2 = 0.5. For the first target DSD, a minimum of three passes was required to achieve feasibility (Table 7). An additional pass Table 7. Model-Predicted Results for the Minimization of Surfactant Usage with Two Target DSDs Generated from a Log-Normal Distribution pressure (bar) total surfactant (wt%)

first pass

second pass

third pass

fourth pass

fifth pass

Ψ2

Ψ3

0.975 0.894

Target DSD Mean = 0.6, Variance = 0.5 800 631 599 − − 0.0097 614 477 393 473 − 0.0098

0.902 0.621

1.135 1.126

Target DSD Mean = 0.5, Variance = 0.5 676 576 557 566 − 0.0098 100 463 500 506 562 0.0099

1.613 1.569

4. CONCLUSIONS A model-based framework for the design of emulsified products in high-pressure homogenizers was developed and experimentally evaluated. The framework was based on a population balance equation model that included mechanistic descriptions of drop breakage and coalescence, surfactant adsorption and desorption kinetics, and the variation of emulsion viscosity with dispersion parameters. Model extensibility over a wide range of homogenization pressures was achieved by introducing empirical relations for the dependence of two model parameters on pressure. Three alternatives of the product design problem were formulated by defining appropriate objective functions: (1) attainment of a target mean size and polydispersity while minimizing surfactant usage; (2) attainment of an entire target drop size distribution (DSD) while minimizing surfactant usage; and (3) minimization of surfactant usage with the constraint that the predicted DSD must be sufficiently close to the target DSD. The optimal design problems were modeled in AMPL and solved with the nonlinear programming code CONOPT to generate the amount of surfactant required and the pressure of each homogenization pass. Our approach required the number of passes to be specified a priori to avoid the complexities of solving a mixed-integer nonlinear programming problem. While we found that acceptable solutions could be obtained with each problem formulation, the second formulation based

offered the advantage of reducing surfactant usage by an additional 8%, which could offer substantial economic benefits for large-scale production. Compared to the alternative objective function Ψ2 + WΨ3, this problem formulation reduced surfactant usage by 7% for three passes and by 8% for four passes while maintaining a small DSD error (see Tables 5 and 7). The second target DSD with smaller mean required a minimum of four passes to establish feasibility, and an additional pass offered very little additional reduction in surfactant usage. Compared to the alternative objective function, this problem formulation only reduced surfactant usage by 4% for four passes due to relatively large amount of surfactant required for stabilization of the small drops (see Tables S3 and 7). We did not perform experimental validation of these predictions. 3.6. Effect of Objective Function Weighting. The previous results were generated with fixed objective function weighting factor: W = 0.1 for the first problem formulation, and W = 0.01 for the second problem formulation. To explore the

Figure 6. Effect of the objective function weighing factor W on (a) the objective functions Ψ2 (in green) and Ψ3 (in blue), and (b) the modelpredicted DSD. G


Article


The breakage frequency function g(v) was assumed to be the sum of two functions representing different breakage mechanisms for turbulent inertia and turbulent viscous shear such that g(v) = g1(v) + g2(v) . Drop breakage due to turbulent inertia was described as

on simultaneous minimization of the DSD error and surfactant usage was preferred for several reasons. The first problem formulation offered the advantage that the target mean size and polydispersity were easily specified, but attainment of these two DSD moments does not ensure that the resulting DSD will be acceptable. The third problem formulation offered the advantage of direct minimization of surfactant costs, but imposition of a hard constraint on the DSD error could result in problem infeasibilities. Therefore, we found the second problem formulation to be most flexible and effective. Despite errors in modeling the highly complex homogenization process, we found that experimental implementation of the optimal solutions resulting from model-based design produced DSDs that provided good agreement with target DSDs. The largest source of error between model-predicted and target DSDs occurred at the tails of the distribution, presumably due to overprediction of drop breakage or underprediction of drop coalescence. Experimental attainment of target DSDs was most difficult at small drop sizes, where model errors seemed to be most prevalent. Despite these limitations, our study represents an important first step in developing model-based design strategies for emulsified products.

■

g1(v) = K1

⎛ σ(1 + ϕ)2 ⎞⎟ 1 exp⎜⎜ −K 2 ⎟ 1+ϕ ρd v 5/9ε 2/3 ⎠ ⎝

(A.2)

where ϕ and ρd are the dispersed phase volume fraction and density, respectively, ε is the energy dissipation rate, and K1 and K2 are adjustable parameters. Drop breakage due to turbulent shear was described as g2(v) = K3 ⎛ ⎜ exp⎜ − ⎜ ⎝

(ερc /ηem)1/2 (1 + ϕ)2

⎛ σ(1 + ϕ)2 ⎞⎟ exp⎜⎜ −K4 ⎟ (ερc ηem)1/2 v1/3 ⎠ ⎝

1/3 ⎞

( π6 v) λK

⎟ ⎟ ⎟ ⎠

(A.3)

where ρc is the continuous phase density, ηem is the calculated emulsion viscosity, σ is the interfacial tension, λK = (ηc3/ερc3)1/4 is the length scale of the smallest size eddies (Kolmogorov length scale), and K3 and K4 are adjustable parameters. The daughter drop distribution function β(v,v′) was constructed such that drops underwent multiple breakage events during a single pass, resulting in the formation of 80 drops per pass.19 The coalescence frequency C(v,v′) was dependent on the drop collision frequency h(v,v′), the coalescence efficiency Λ(v,v′), and the surfactant surface coverage of drops Γ as

APPENDIX

Population Balance Equation

We utilized our most recent PBE model of high-pressure homogenization 21 to develop model-based design strategies for achieving target DSDs. Below the PBE model is briefly presented with a focus on the salient points for this study. The interested reader is directed to our previous studies for more details on the model formulation as well as the numerical methods used for model solution and parameter estimation.19−21 As the light-scattering device used in this study directly measured drop volume distributions, the PBE model was formulated in terms of drop volume fraction rather than drop number,

2 ⎛ Γ ⎞ C(v , v′, Γ) = ⎜1 − ⎟ h(v , v′)Λ(v , v′) Γ∞ ⎠ ⎝

∂N (v , t ) = − g ( v ) N (v , t ) ∂t ∞ β(v , v′)g (v′)N (v′ , t ) dv′ +v v v′ ∞ C(v , v′)N (v′ , t )V t dv′ − N (v , t ) 0 v′ v C(v − v′ , v′)N (v − v′ , t )N (v′ , t )V v t dv′ + 2 0 v′(v − v′)

h(v , v′) = K5

(A.4)

ε1/3 (v 2/3 + v′2/3 )(v 2/9 + v′2/9 )1/2 (1 − ϕ) (A.5)

∫

⎡ −η ρ ε ⎛ 1/3 1/3 ⎞4 ⎤ v v′ ⎟⎥ Λ(v , v′) = exp⎢ 2 em c 3 ⎜ 1/3 ⎢⎣ σ (1 + ϕ) ⎝ v + v′1/3 ⎠ ⎥⎦

∫

∫

(A.6)

where Γ∞ is the maximum surfactant surface coverage and K5 is an adjustable parameter. The breakage and coalescence frequency functions both depend on the energy dissipation rate, which was modeled as14,15

(A.1)

where v is the volume of the drop; N(v,t) dv is the volume fraction of drops with volume in the range [v,v + dv] per unit volume of dispersion at time t; Vt is the conserved total volume of the drops per unit volume of dispersion; g(v) is the breakage frequency which represents the fraction of drops of volume v breaking per unit time; β(v,v′) is the daughter drop distribution function which represents the probability of forming a daughter drop of size v from breakage of a mother drop of size v′; and C(v,v′) is the coalescence frequency which represents the rate at which drops of size v and size v′ coalesce. We assumed that drop breakage and coalescence occurred only in the homogenizer valve gap, and that these mechanisms were uniform throughout the gap region.

ε=

ΔPQ Vdiss

(A.7)

where ΔP is the homogenization pressure, Q is the volumetric flow rate, and Vdiss is the valve gap volume, which depends on valve gap distance.31 The emulsion viscosity ηem was calculated using an empirical equation 32 that accounts for the effects of ϕ, ηd, and ηc.21,29 Surfactant Adsorption Model

The dynamic surface coverage model had the form H


Article

Industrial & Engineering Chemistry Research dΓ = K a(Γ∞ − Γ)Cs − Kd Γ − dt ∞ ⎛ 1 ∂N ⎞ Γ ⎜ ⎟ dv ; ∞ ∫ (N /d) dv 0 ⎝ d ∂t ⎠

∫

d=

(3) Israelachvili, J. The science and applications of emulsions - An overview. Colloids Surf., A 1994, 91, 1−8. (4) McClements, D. J. Food Emulsions: Principles, Practice, and Techniques; CRC Press: Boca Raton, FL, 2005. (5) Marti-Mestres, G.; Nielloud, F. Emulsions in health care applications−an overview. J. Dispersion Sci. Technol. 2002, 23 (1−3), 419−439. (6) Walstra, P. What happens during homogenization. Netherlands Milk Dairy J. 1980, 34 (2), 135. (7) Becher, P. Encyclopedia of Emulsion Technology, Vol. I; Marcel Dekker: New York, 1983. (8) Walstra, P. Principles of emulsion formation. Chem. Eng. Sci. 1993, 48 (2), 333−349. (9) Ramkrishna, D. Population Balances: Theory and Applications to Particulate Processes in Engineering; Academic Press: New York, 2000. (10) Soon, S.; Harbidge, J.; Titchener-Hooker, N.; Shamlou, P. Prediction of drop breakage in an ultra high velocity jet homogenizer. J. Chem. Eng. Jpn. 2001, 34 (5), 640−646. (11) Maguire, L.; Zhang, H.; Shamlou, P. Preparation of small unilamellar vesicles (SUV) and biophysical characterization of their complexes with poly-L-lysine-condensed plasmid DNA. Biotechnol. Appl. Biochem. 2003, 37, 73−81. (12) Kelly, W.; Muske, K. Optimal operation of high-pressure homogenization for intracellular product recovery. Bioprocess Biosyst. Eng. 2004, 27 (1), 25−37. (13) Tcholakova, S.; Vankova, N.; Denkov, N.; Danner, T. Emulsification in turbulent flow: 3. Daughter drop-size distribution. J. Colloid Interface Sci. 2007, 310 (2), 570−589. (14) Vankova, N.; Tcholakova, S.; Denkov, N.; Vulchev, V.; Danner, T. Emulsification in turbulent flow 2. Breakage rate constants. J. Colloid Interface Sci. 2007, 313 (2), 612−629. (15) Vankova, N.; Tcholakova, S.; Denkov, N.; Ivanov, I.; Vulchev, V.; Danner, T. Emulsification in turbulent flow - 1. Mean and maximum drop diameters in inertial and viscous regimes. J. Colloid Interface Sci. 2007, 312 (2), 363−380. (16) Hakansson, A.; Tragardh, C.; Bergenstahl, B. Dynamic simulation of emulsion formation in a high pressure homogenizer. Chem. Eng. Sci. 2009, 64 (12), 2915−2925. (17) Hakansson, A.; Tragardh, C.; Bergenstahl, B. Studying the effects of adsorption, recoalescence and fragmentation in a high pressure homogenizer using a dynamic simulation model. Food Hydrocolloids 2009, 23 (4), 1177−1183. (18) Håkansson, A.; Innings, F.; Trägårdh, C.; Bergenståhl, B. A highpressure homogenization emulsification model−improved emulsifier transport and hydrodynamic coupling. Chem. Eng. Sci. 2013, 91, 44− 53. (19) Maindarkar, S.; Raikar, N.; Henson, M. Incorporating emulsion drop coalescence into population balance equation models of high pressure homogenization. Colloids Surf., A 2012, 396, 63−73. (20) Maindarkar, S. N.; Bongers, P.; Henson, M. A. Predicting the effects of surfactant coverage on drop size distributions of homogenized emulsions. Chem. Eng. Sci. 2013, 89, 102−114. (21) Maindarkar, S.; Hoogland, H.; Henson, M. Predicting the combined effects of oil and surfactant concentrations on the drop size distributions of homogenized emulsions. Colloids Surf., A 2015, 467, 18−30. (22) Crowley, T. J.; Meadows, E. S.; Kostoulas, E.; Doyle Iii, F. J. Control of particle size distribution described by a population balance model of semibatch emulsion polymerization. J. Process Control 2000, 10 (5), 419−432. (23) Nagy, Z.; Aamir, E.; Rielly, C. Internal fines removal using population balance model based control of crystal size distribution under dissolution, growth and nucleation mechanisms. Cryst. Growth Des. 2011, 11 (6), 2205−2219. (24) Aamir, E.; Nagy, Z. K.; Rielly, C. D. Optimal seed recipe design for crystal size distribution control for batch cooling crystallisation processes. Chem. Eng. Sci. 2010, 65 (11), 3602−3614.

⎛ 6v ⎞1/3 ⎜ ⎟ ⎝π⎠

0

(A.8)

where Ka and Kd are the surfactant adsorption and desorption rate constants, respectively, Cs is the concentration of free surfactant, and d is the emulsion drop diameter. The first two terms on the right-hand side of the equation account for the change in surface coverage due to surfactant adsorption and desorption, while the last term accounts for the change in coverage due to variations in total interfacial area resulting from drop breakage and coalescence. The following mass balance equation was used to calculate the free surfactant concentration, dCs 6ϕ d ⎛ Γ ⎞ =− ⎜ ⎟ dt 1 − ϕ dt ⎝ d32 ⎠

(A.9)

where d32 is the Sauter mean diameter of the emulsion droplets. The initial free surfactant concentration Csinit was calculated using the following mass balance equation, Cs,init = Cs,input −

6ϕ Γinit 1 − ϕ d32,pre

(A.10)

where d32,pre is Sauter mean diameter of the coarse premix emulsion and Cs,input is the total surfactant concentration used to prepare the emulsion. The Langmuir isotherm equation was used to calculate the initial value of the surface coverage Γinit, Γinit =

Γ∞Cs,init Cs,init + C1/2

(A.11)

where C1/2 is the free surfactant concentration corresponding to half of the maximum surface coverage.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.5b01195. Tables S1−S4, listing bulk physical properties of emulsion system, objective function values at six different homogenization pressures, results to achieve target DSD, and effect of objective function weighing factor (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*Phone: (413) 545-3481. Fax: (413) 545-1647. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS We acknowledge financial support from Unilever Foods and National Science Foundation Grant DMI-0531171. REFERENCES

(1) Becher, P. Emulsions: Theory and Practice; Oxford University Press: New York, 2001. (2) Chappat, M. Some applications of emulsions. Colloids Surf., A 1994, 91, 57−77. I


Article

Industrial & Engineering Chemistry Research (25) Majumder, A.; Nagy, Z. K. Prediction and control of crystal shape distribution in the presence of crystal growth modifiers. Chem. Eng. Sci. 2013, 101, 593−602. (26) Kumar, S.; Ramkrishna, D. On the solution of population balance equations by discretization 1. A fixed pivot technique. Chem. Eng. Sci. 1996, 51 (8), 1311−1332. (27) Fourer, R.; Gay, D. M.; Kernighan, B. W. AMPL: A Modeling Language for Mathematical Programming; Brooks/Cole Publishing Co.: Pacific Grove, CA, 2003. (28) Drud, A. S. CONOPT - A large-scale GRG code. ORSA J. Comput. 1994, 6 (2), 207−216. (29) Maindarkar, S.; Dubbelboer, A.; Meuldijk, J.; Hoogland, H.; Henson, M. Prediction of emulsion drop size distributions in colloid mills. Chem. Eng. Sci. 2014, 118, 114−125. (30) Dubbelboer, A.; Janssen, J.; Hoogland, H.; Mudaliar, A.; Maindarkar, S.; Zondervan, E.; Meuldijk, J. Population balances combined with computational fluid dynamics: A modeling approach for dispersive mixing in a high pressure homogenizer. Chem. Eng. Sci. 2014, 117, 376−388. (31) Raikar, N.; Bhatia, S.; Malone, M.; Henson, M. Experimental studies and population balance equation models for breakage prediction of emulsion drop size distributions. Chem. Eng. Sci. 2009, 64 (10), 2433−2447. (32) Jansen, K.; Agterof, W.; Mellema, J. Viscosity of surfactant stabilized emulsions. J. Rheol. 2001, 45, 1359−1371.

J


Achieving Target Emulsion Drop Size Distributions Using Population

Recommend Documents