Brazilian Journal of Chemical Engineering MODELING PARTICLE SIZE DISTRIBUTION IN HETEROGENEOUS POLYMERIZATION SYSTEMS USING MULTIMODAL LOGNORMAL FUNCTION

This work evaluates the usage of the multimodal lognormal function to describe Particle Size Distributions (PSD) of emulsion and suspension polymerization processes, including continuous reactions with particle re-nucleation leading to complex multimodal PSDs. A global optimization algorithm, namely Particle Swarm Optimization (PSO), was used for parameter estimation of the proposed model, minimizing the objective function defined by the mean squared errors. Statistical evaluation of the results indicated that the multimodal lognormal function could describe distinctive features of different types of PSDs with accuracy and consistency.


INTRODUCTION
The particle size distribution (PSD) is one of the most important characteristics of polymer latexes/resins, since properties such as viscosity, maximum solids content, adhesion and drying time depend on the profile of this distribution (Vale and McKenna, 2005).Furthermore, the understanding of the dynamics of PSD is very important for process monitoring and final quality control of heterogeneous polymerization processes (Hosseini et al., 2013).Therefore, PSD modeling has received great attention with the proposal of models such as Population Balance Equations (PBE) (Machado et al., 2000;Araújo et al., 2001;Kiparissides et al., 2002;Immanuel et al., 2002;Coen et al., 2004).However, according to Vale and McKenna (2005), even when these difficulties in coagulation modeling are overcome, the PBEs are difficult to solve if they include kinetic and/or hydrodynamic complete models.Despite recent advances in techniques for experimental determination, difficulties in determining important variables related to the quality and productivity of polymers still remain (Machado et al., 2007).In the case of PSD, these limitations are associated with the existence of polydisperse distributions and the long times often required for sample preparation and analysis.
In this way, simpler and thus less computationally expensive alternatives for the description of complex PSDs are very useful for different types of applications as, for instance, the development of soft-sensors for on-line PSD monitoring (Clementi et al., 2013), the implementation of optimization and control strategies or the simulation of heterogeneous polymerize-tion reactions using CFD tools (Elgebrandt et al., 2006;Pohn et al., 2013).These alternatives require a model describing the PSD.
According to Limbert et al. (2001) the multimodal lognormal function has proved to be efficient and flexible in describing wide or narrow, unimodal and multimodal distributions in many particulate systems applications.
This diversity of applications can be observed in works such as Zhao et al. (2003) who investigated the dynamics of soot particle size distributions.The numerical simulation using a kinetic model proposed previously and a stochastic approach to solve aerosol dynamics equations showed a bimodal soot particle size distribution function, the consequence of the interplay between particle-particle coagulation and particle nucleation.In this same respect, Johnsson et al. (2010) claimed that soot particle distributions generally consist of a range of different particle sizes, and the lognormal size distribution has been shown to be a good approximation for various flame conditions; Hwang and Choi (2006) investigated the predictive potential of the unimodal lognormal distribution model for estimating the water retention curves, and this model was evaluated for a broader range of soil; Yuan et al. (2011) described the characteristics of atmospheric aerosols using a multilognormal distribution model with parameters estimated by Particle Swarm Optimization.The model was validated and analyzed by comparing with the measured optical properties of several bands; Pujol and Pinto (2011) approached fatigue life prediction under step-stress conditions by comparing a standard approach based on the lognormal distribution function implemented and fit to experimental data, with an approach for fatigue life prediction from neural networks.Both models were optimized by differential evolution associated to the Newton-Raphson-Method, using a maximum likelihood estimator.
Given the above, this paper aims to evaluate the use of a multimodal lognormal function to describe the behavior of particle size distributions in heterogeneous polymerization systems, using a Particle Swarm Optimization algorithm for parameter estimation.

Experimental Particle Size Distributions
In order to evaluate whether the proposed multimodal lognormal function was able to describe real PSDs, experimental PSDs of emulsion and suspension polymerization reactions were used for the parameter estimation of the multimodal lognormal function.The choice of the distributions used in this work was based on their different characteristics and on the existence of multimodalities.
The reaction techniques, monomer system, reactor type and operation mode of the experimental PSD used in this work are summarized in Table 1.Four PSDs (PSD01 to PSD04), among which there were unimodal and multimodal distributions, were obtained from the work of Araújo (1999).In short, this author investigated the effects of operational conditions (temperature, initiator and surfactant concentrations) on the PSD in an emulsion copolymerization of vinyl acetate/Veova10 in a continuous loop reactor with high solids content.The reactor pre-feeding with an aqueous (surfactant, protective colloid and Na2S2O5) and an organic phase (monomer mixture with vinyl acetate to Veova 10 weight ratio of 75/25) caused a fast increase of the particle number at the reaction beginning.As a consequence, the distributions usually show a right skewness.The higher solubility of the oligomeric radical in the aqueous phase, due to the vinyl acetate, affects the main "locus" of the polymerization, promoting re-nucleations.
Another PSD (PSD05) was obtained from the work of Zubitur and Asua (2001) which investigated the factors affecting kinetics and coagulation during the emulsion polymerization of styrene and butyl acrylate, in a semi-batch stirred tank reactor.The PSD showed some bimodality, a peak of small particles, probably resulting from secondary nucleation in the presence of seed particles (shoulder of bigger particles).Mallikarjunan et al. (2010) studied the emulsion copolymerization of vinyl acetate and butyl acrylate, conducted in a semi-batch stirred tank reactor, and another PSD (PSD06) was utilized from the work of these authors.They report a bimodal PSD determined with CHDF, since the system was fed with surfactant in excess.
In addition to the PSDs of emulsion polymerization reactions, PSDs of suspension polymerization reactions of methyl methacrylate were also evaluated.These PSDs are available in the work of Jahanzad et al. (2004), and they were included in this work with the objective of investigating the adaptability of the multimodal lognormal function to other heterogeneous polymerization techniques.In this paper, the authors carried out extensive research on the influence of process variables (temperature, inhibition, stabilizer and initiator concentrations, agitation) on the PSD in a batch stirred tank reactor.The distributions obtained by Jahanzad et al. (2004) are basically unimodal.In general, some factors promoted the growth of the mean particle size and the broadening of the distributions, such as the increase of the reaction temperature, decrease of the stirring rate, decrease of the stabilizer concentration and absence of inhibitor.
In order to normalize the distributions range, a numerical integration (trapezoid method) was used to calculate the area under each PSD.The different types of distributions (number/weight/volume density distributions) of these reference works were normalized by the division of the PSD by the calculated area.In this way, the sum of the relative frequencies was equal to one for each PSD.

Multimodal Lognormal Function
The lognormal distribution is a model with two parameters, whose density function is expressed by Equation ( 1) (Crow and Shimizu, 1988): .ln 2 where x is the independent variable, g  and g  re- present the geometric mean and geometric standard deviation of the distribution, respectively.For frequency distributions, xk values are weighted by their relative frequencies fk with k = 1, 2 ..., n.In this case, the geometric mean and geometric standard deviation can be defined by Equations ( 2) and (3), respectively (Crow and Shimizu, 1988): where However, in real situations, these distributions can be composed of more than one mode, being bimodal, trimodal or even n-modal.In such cases, the function can be described by a weighted sum of lognormal functions (Zender, 2010).Equation ( 4) describes the multimodal lognormal function: 1 MLF( , , , ) .( , , ) where ( , , , ) tion mode i; n number of modes and w i is the fraction Taking into account the number of modes of the experimental distributions used in this work, the Bimodal Lognormal Function (BLF) and Trimodal Lognormal Function (TLF) were chosen for modeling experimental PSDs.Equation ( 5) provides a representation of the function used, with the particle diameter d as independent variable: Therefore, there were always five or eight parameters to be estimated for BLF and TLF respectively, since there were three parameters to be estimated for each mode ( ln i g  , ln i g  and i w ).

Parameter Estimation
The parameter estimation procedure of the MLFs was supported by Particle Swarm Optimization (Kennnedy and Eberhart, 1995).This algorithm uses its intelligent strategy to randomly minimize the objective function defined by the sum squared errors [SSE], according to Equation ( 6): where exp k f is the experimental frequency value in k; ( ) k MLF d is the frequency calculated by the multi- modal lognormal function in k and n is the number of experimental points.
In short, according to Yang (2014), the PSO strategy considers that a particle i in an interaction k, moves through the search space with two attributes: position (X) and velocity (V) vectors, represented by equations ( 7) and (8), respectively: where, k i P is the best previous position along the d th dimension of particle i in iteration k th , while k g P is the best previous position among all the particles along the d th dimension in iteration k th .According to Blum and Merkle (2008) the velocity vector is determined by three terms: the "momentum", the "cognitive", and the "social" term.The "momentum" term w.vi represents the previous velocity which is used to carry the particle in the direction it has travelled so far; the "cognitive" term c1r1(Pl k -Xi k ), represents the tendency of the particle to return to the best position it has visited so far; the "social" term c2r2(Pg k -Xi k ), represents the tendency of the particle to be attracted towards the position of the best position found by the entire swarm.Parameters c1 and c2 are acceleration coefficients, while r1 and r2 are two independent random numbers uniformly distributed in [0,1].
The procedure proposed by Jiao et al. (2008), which considers the implementation of the dynamic inertial weight (Equation ( 9)) as a way to control the particle momentum, weighing the contribution of the previous velocity and thus, basically controlling how much memory of the previous flight direction will influence the new velocity, was used in this work.In this way, algorithm parameters follow the references of Jiao et al. (2008) and Van den Bergh (2002) who investigated the convergence criteria of the algorithm.
As this is a strategy that performs a random search, it is necessary to set upper and lower limits of the parameters to be adjusted.So, the limits of the parameters ln g  were defined from the limits of the diameter variable on a logarithm scale; in turn the limits of parameters ln g  were defined based on previous works on modeling of particle size distribution using the lognormal function (Limpert et  al., 2001).Table 2 shows the employed search ranges, and the initial parameter estimates generated randomly from a uniform distribution, on the PSO.Two stopping criteria of the algorithm were used: the maximum number of iterations (8000) or evaluation of Equation (10): where g Y and z Y are, respectively, the best result found and the best z last results.The parameter estimation strategy was implemented in Fortran 90 language and parameter estimations were carried out in triplicate.Reported optimal parameters are arithmetic mean values of the three estimations.Standard deviations of each parameter are also reported, representing the algorithm variability due to the random feature of the search.Therefore, it is important to point out that standard Brazilian Journal of Chemical Engineering Vol. 33, No. 03, pp. 469 -478, July -September, 2016 deviation values have no relationship with the confidence interval of the optimal parameter.
Statistical interpretation of the results was performed through the correlation coefficient [r] between experimental and calculated frequency distributions.In addition, the Kolmogorov-Smirnov test was used to compare experimental frequency distributions with estimated frequency distributions.The Kolmogorov-Smirnov test and Chi-square test are widely used in the scientific literature and the first one is more efficient for continuous frequency distributions (Evren and Tuna, 2012).

RESULTS AND DISCUSSION
The first set of PSDs, PSD01 to PSD04, corresponds to the evolution during a continuous vinyl acetate/Veova 10 emulsion copolymerization reaction conducted in a tubular loop reactor.Due to successive nucleations this reaction presents broad and often multimodal PSDs (Rawlings and Ray, 1988).This feature generates consequences for the determination of the search interval of the parameter geometric mean.Besides, the positively skewed modes affect strongly the geometric standard deviation, minimizing it in this case.
Figure 1 compares the experimental PSDs with and PSDs predicted using a BLF.It may be observed that the BLF could fit nicely the broad monomodal and positively skewed distribution PSD01, as well as different bimodal distributions with different mode locations and relative fractions (PSD02 and PSD04).As expected, since the proposed model provides the generation of two modes, it was not able to represent all the three modes of PSD03.Nevertheless, even in this case the BLF was able to predict correctly the location, width and relative intensity of the two bigger particle populations, while the less significant mode of the distribution was ignored.
When a trimodal lognormal function (TLF) is used to represent the evolution of the PSD during this continuous emulsion polymerization reaction, Figure 2, the TLF fits perfectly the experimental distributions, no matter if mono-, bi-or trimodal, showing better agreement than with BLF for all PSDs.This result is due to the higher number of parameters that ensures the representation of trimodal distributions and favors data fitting when the distribution is simpler and, thus, with overlapping modes, with very close parameters, favoring the parameter estimation procedure.The higher number of parameters, though, can lead to a high correlation among the TLF parameters.In order to verify the adaptability of the TLF to other heterogeneous polymerization techniques, this function was used to predict the PSDs of methyl methacrylate of suspension polymerization reactions (PSD07 and PSD08, Jahanzad et al., 2004).Again the TLF predicted correctly both mono-as well as multimodal PSDs, despite the different order of magnitude of the particle diameters and forms of the PSDs, when compared to those of the previously evaluated emulsion polymerization reactions.
Table 3 and 4 show the average and standard deviation values for the parameters set for BLF and TLF, respectively.These tables also show the results of the objective function value (SSE) for each case.The error in the parameter estimation was quantified from the results obtained in minimization triplicates.It can be verified that all standard deviation and SSE values are low, confirming the ability of the algorithm to reproduce consistently the PSDs, thus, corroborating PSO robustness.Finally, a comparison between standard deviations presented in Table 3 and 4 shows that they are higher when using TLF, suggesting a greater correlation between the parameters of TLF, since in this function there exists a greater number of possible parameter combinations that result from the minimization of the function.
It can be pointed out that scale differences in the particle diameters between emulsion polymerization [nm] and suspension polymerization [µm] did not affect the proposed model for PSDs, and also that the TFL model was very sensitive to small variations in the distribution like the presence of a very small population of small particles in PSD08 in Figure 4.  (Jahanzad, 2004) and predicted PSDs of batch suspension methyl methacrylate polymerization reaction in a stirred tank reactor using a trimodal lognormal function (TLF).

Brazilian Journal of Chemical Engineering
Evaluation of the correlation coefficient (r) was consistent with the results shown in Figures 1 to Figure 4, with r ≥ 0.98 in all cases using BLF and r ≥ 0.99 in all cases using TLF, evidencing a strong correlation between experimental and predicted frequency distributions.The Kolmogorov-Smirnov test and the established critical values resulted in a dcritical value equal to 0.210 for a significance level (α) of 5% and for a number of observed points (n) equal to 40, used in this work.As the largest value dmáx evaluation for PSDs predicted is equal to 0.0147, the hypothesis that the data refer to the same distribution should be accepted, increasing the confidence in the performed parameter estimation.
Finally, it can be pointed out that all variability in the evaluated PSDs could be perfectly assimilated by TLF, and, except for the smaller third population of large particles of PSD03, also by the BLF.The proposal did not present inconsistencies even for the different particle diameter ranges commonly provided by different polymerizations techniques.This difficulty was overcome by using normalized relative frequencies.

CONCLUSIONS
A multimodal lognormal function was proposed to describe the particle size distributions in heterogeneous polymerization systems, with parameters estimated by the global optimization algorithm Particle Swarm Optimization.Distributions with different characteristics, unimodal or multimodal, wide or narrow, but all originating from polymerization reactions in emulsion or suspension were modeled.
High correlations between experimental and predicted frequency distributions using bimodal and trimodal lognormal functions and the Kolmogorov-Smirnov test were obtained for all particle size distributions, confirming that the proposed function can describe successfully particle size distributions with different characteristics.
Even when different particle diameter ranges were used, the proposed model could perfectly represent the distribution profile.Furthermore, the accuracy of the parameter estimation was evidenced by low standard deviation values, attesting to the ability of Particle Swarm Optimization to reproduce consistent results.
Comparing the trimodal lognormal function with the bimodal lognormal function, the use of the former is justified only when the investigated particle size distributions are very polydisperse with three relevant modes.

ACKNOWLEDGMENTS
The authors thank the financial support from CNPq -Conselho Nacional de Desenvolvimento Científico e Tecnológico.

Figure 3
Figure3compares the predictions of a TLF with the final experimental PSDs of semi-batch styrene/butyl acrylate (PSD05,Zubitur and Asua, 2001) and vinyl acetate/butyl acrylate (PSD06,Mallikarjunan et al., 2010) emulsion copolymerization reactions and the TLF was able to predict correctly both bimodal PSDs, shoulder of bigger particles and two separate populations.In order to verify the adaptability of the TLF to other heterogeneous polymerization techniques, this function was used to predict the PSDs of methyl methacrylate of suspension polymerization reactions