# Dynamic evolution of the particle size distribution in suspension polymerization reactors: a comparative study on Monte Carlo and sectional grid methods.

INTRODUCTIONAn important property of suspension polymerization processes is the particle size distribution (PSD), which controls key aspects of the process and affects the end-use properties of the product. Suspension polymerization processes are generally characterized by PSDs that can vary in time with respect to the mean particle size as well as to the PSD form (i.e., broadness and/or skewness of the distribution, unimodal and/or bimodal character, etc.). The quantitative calculation of the evolution of the PSD presupposes a good knowledge of the droplet/particle breakage and coalescence mechanisms. These mechanisms are coupled with the reaction kinetics, thermodynamics and other microscale phenomena, including mass and heat transfer between the different phases present in the system.

The time evolution of the PSD is commonly obtained from the solution of a population balance equation (PBE), governing the dynamic behaviour of the dispersed liquid monomer droplets that are being polymerized to solid polymer particles. The numerical solution of the dynamic PBE is a notably difficult problem due to a number of numerical complexities and model uncertainties regarding the particle breakage and coalescence mechanisms and are often poorly understood. It commonly requires the discretization of the particle volume domain into a number of discrete elements and the subsequent numerical solution of the resulting system of stiff, nonlinear differential or algebraic/differential equations (DAEs). In the open literature, a number of numerical methods have been reported on the steady-state and dynamic solutions of the general PBE. These include the full discrete method (Hidy, 1965), the method of classes (Marchal et al., 1988; Chatzi and Kiparissides, 1992), the discretized PBE (DPBE) methods (Batterham et al., 1981; Hounslow et al., 1988), the fixed and moving pivot DPBE methods (Kumar and Ramkrishna, 1996a,b), the high-order DPBE methods (Bleck, 1970; Gelbard and Seinfeld, 1980; Sastry and Gaschignard, 1981; Landgrebe and Pratsinis, 1990), the orthogonal collocation on finite elements (OCFE) (Gelbard and Seinfeld, 1979), the Galerkin method (Nicmanis and Hounslow, 1998) and the wavelet-Galerkin method (Chen et al., 1996). In the reviews of Ramkrishna (1985), Dafniotis (1996), and Kumar and Ramkrishna (1996a), the various numerical methods available for solving the general PBE are described in detail. Moreover, extensive comparative studies have been presented in the publications of Kostoglou and Karabelas (1994, 1995), Nicmanis and Hounslow (1996) and in a recent series of papers by Kiparissides and coworkers (Alexoppoulos et al., 2004; Alexopoulos and Kiparissides, 2005; Roussos et al., 2005; Meimaroglou et al., 2006). On the basis of the conclusions of these studies, the DPBE method of Litster et al. (1995), the fixed pivot method of Kumar and Ramkrishna (1996a), the Galerkin and the orthogonal collocation on finite-element methods were found to be the most accurate and stable numerical techniques for the numerical solution of the PBE.

The dynamic evolution of the PSD in a particulate process can also be obtained via stochastic Monte Carlo (MC) simulations. Spielman and Levenspiel (1965) were the first to employ an MC approach to study the effect of particle coalescence in a two-phase particulate reactive system in well-mixed reactors. Later, Shah et al. (1977) developed a general MC algorithm for time varying particulate processes. In 1981, Ramkrishna (Ramkrishna, 1981) established the precise mathematical connection between population balances and the MC approach. In MC simulations, the dynamic evolution of the PSD is inferred by the properties of a finite number of sampled particles. Based on the method employed for the determination of the sampling time step, MC simulations can be distinguished into time-driven (Domilovskii et al., 1979; Liffman, 1992; Debry et al., 2003) and event-driven ones (Garcia et al., 1987; Smith and Matsoukas, 1998; Tandon and Rosner, 1999; Kruis et al., 2000; Efendiev and Zachariah, 2002). In regard to the total number of simulated particles, MC methods can be further classified into constant number and constant volume MC methods. A more detailed description of the characteristics of the various MC formulations can be found in a number of publications (Maisels et al., 2004; Meimaroglou et al., 2006; Zhao et al., 2007).

In the present study, two numerical approaches, namely, the fixed pivot technique (FPT) and a stochastic Monte Carlo (MC) algorithm are applied for solving the general PBE, governing the PSD developments in a methyl methacrylate (MMA) free-radical batch suspension polymerization reactor, in terms of the process conditions (i.e., monomer to water volume ratio, temperature, type and concentration of stabilizer, energy input into the system, etc.) and polymerization kinetics. To the best of our knowledge, this is the first time that the two numerical methods (i.e., the FPT and MC) are applied for the calculation of the dynamic evolution of the PSD in a batch suspension polymerization reactor, using a comprehensive model taking into account all the physical and chemical phenomena in the polymerization process. The validity of both numerical methods is examined via a direct comparison of model predictions with experimental measurements on the average particle diameter and the droplet/particle size distributions for both non-reactive and reactive systems.

MODEL DEVELOPMENTS

In suspension polymerization, the monomer is commonly dispersed in the continuous aqueous phase by the combined action of surface-active agents (i.e., inorganic or/and water soluble polymers) and agitation. All the reactants (i.e., monomer, initiator(s), etc.) reside in the organic or "oil" phase. The polymerization occurs in the monomer droplets that are progressively transformed into sticky, viscous monomer-polymer particles and, finally, into rigid, spherical polymer particles in the size range of 50-500 [micro]m (Kiparissides, 1996). The polymer content, in the fully converted suspension, is typically 30-50% (w/w). Large quantities of commercially important polymers (e.g., polystyrene and its copolymers, poly(vinyl chloride), poly(methyl methacrylate), etc.) are produced by the suspension polymerization process.

One of the most important issues in the operation of a suspension polymerization reactor is the control of the final PSD (Yuan et al., 1991). The initial monomer droplet size distribution (DSD) as well as the final PSD largely depend on the type and concentrations of the surface active agents, the quality of agitation and the physical properties (e.g., densities, viscosities, and interfacial tension) of the continuous and dispersed phases. The transient droplet/particle size distribution is controlled by two dynamic processes, namely, the droplet breakage and coalescence mechanisms. The former mainly occurs in regions of high shear stress (i.e., near the agitator blades) or as a result of turbulent velocity and pressure fluctuations along the surface of a droplet. The latter is either increased or decreased by the turbulent flow field and can be assumed to be negligible for dilute dispersions, at sufficiently high concentrations of surface active agents (Chatzi and Kiparissides, 1992).

In general, the dynamic evolution of the PSD in a batch suspension polymerization reactor can be calculated from the solution of the governing PBE. The distribution of the droplets/particles is considered to be continuous in the volume domain and is commonly described in terms of the number density function, n(V, t). Thus, n(V, t) dV represents the number of particles per unit volume in the differential volume size range (V, V + dV). For a dynamic particulate system, undergoing simultaneous particle breakage and coalescence, the rate of change of the number density function with respect to time and volume can be described by the following non-linear integro-differential PBE (Kiparissides et al., 2004):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1)

The first term on the right-hand side (r.h.s.) of Equation (1) represents the generation of droplets in the size range (V, V + dV) due to breakage. [beta](U, V) is a daughter droplet breakage function, accounting for the probability that a droplet of volume V is formed via the breakage of a droplet of volume U. The function u(U) denotes the number of droplets formed by the breakage of a droplet of volume U and g(U) is the breakage rate of droplets of volume U. The second term on the r.h.s. of Equation (1) represents the rate of generation of droplets in the size range (V, V + dV) due to droplet coalescence. k(V, U) is the coalescence rate between two droplets of volume V and U. Finally, the third and fourth terms represent the droplet disappearance rates due to breakage and coalescence mechanisms, respectively. Equation (1) will satisfy the following initial condition at t = 0:

n(V, 0) = [n.sub.0](V) (2)

where [n.sub.0] (V) is the initial droplet size distribution of the dispersed phase.

The numerical solution of the PBE (Eq. 1) presupposes the knowledge of the breakage and coalescence rate functions. In the open literature, several forms of g(V) and k(V, U) have been proposed to describe the droplet breakage and coalescence rate functions in liquid-liquid dispersions (Narsimhan et al., 1979; Sovova, 1981; Chatzi et al., 1989; Coulaloglou and Tavlarides, 1977). More details regarding the selection of the above kernels can be found in Appendix.

In regard with the droplet/particle breakage and coalescence phenomena, the suspension polymerization can be divided into three phases (Hamielec and Tobita, 1992; Kotoulas and Kiparissides, 2006, Figure 1). During the initial low-conversion (i.e., low-viscosity) phase, the monomer droplet breakage is the dominant mechanism. As a result, the initial DSD shifts to smaller sizes. During the second sticky-phase of polymerization, the droplet breakage rate decreases while the droplet/particle coalescence becomes the dominant mechanism. Thus, the average particle size increases. In the final phase, the PSD reaches its identification point while the polymer particle size slightly decreases due to shrinkage (i.e., the polymer density is greater than the monomer one).

As can be seen in Figure 1, the second phase in the suspension polymerization of MMA coincides with the manifestation of the diffusion-controlled phenomena in the polymerization kinetics. This manifestation is marked by a sharp increase in the polymerization rate, followed by an increase in the monomer conversion and the reaction heat. This high particle polymerization rate can lead to particle overheating (i.e., the individual temperature of a particle, [T.sub.p], can exceed the temperature of the continuous phase, [T.sub.c]) due to heat-transfer limitations (see Figure 2). Particle overheating can directly affect the droplet/particle size distribution developments through the simultaneously occurring changes in the individual particles' temperatures and physical properties (i.e., densities, viscosities, and interfacial tensions). Moreover, particle overheating can considerably increase the particle agglomeration rate (i.e., via the increase of particles' stickiness). In general, particle overheating (i.e., the temperature difference ([T.sub.p] - [T.sub.c])) can be calculated by the following pseudo-steady state equation:

([T.sub.p] - [T.sub.c]) = (-[DELTA] [H.sub.r])[V.sub.p][R.sub.p] / h[A.sub.p] (3)

where (-[DELTA][H.sub.r]) denotes the heat of polymerization. [V.sub.p] is the volume of the particle and [R.sub.p] is the particle polymerization rate. Finally, h and [A.sub.p] are the heat-transfer coefficient and the particle's external surface area, respectively. In the present study, different correlations (i.e., Ranz and Marshall, 1952, etc.) were considered for calculating the particle's heat transfer coefficient.

[FIGURE 1 OMITTED]

[FIGURE 2 OMITTED]

During the second sticky-phase of polymerization, the effect of particle overheating on the droplet/particle coalescence rate was taken into account via the dependence of the [k.sub.c] term (see Eq. 36) on the polymerization rate. According to the original work of Achilias and Kiparissides (1992), the onset of the gel-effect in the free-radical polymerization of MMA usually occurs in the monomer conversion range of (25-40%) (see the inset in Figure 2). In the region A, the termination rate constant becomes diffusion-controlled since it involves the simultaneous diffusion and reaction of two "live" polymer chains. This results in a decrease of the termination rate constant, [k.sub.tc], as well as an increase of the gel-effect index, GI (GI = ([R.sub.p]/ [R.sub.p0] - 1)). The maximum value of the gel-effect index marks the crossover from region A to region B. In the region B, a decrease in the initiator efficiency takes place, due to the cage-effect, leading to an analogous decrease in the polymerization rate and, ultimately, in the gel-effect index, GI.

In accordance to the changes in the gel-effect index, GI, the droplet/particle collision frequency rate constant, [k.sub.c], was expressed, in the regions A and B, as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4)

where t and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are the time and the initial value of the collision frequency constant, respectively. The parameters [p.sub.1], [p.sub.2], [p.sub.3], and [p.sub.4], appearing in the semi-empirical equation (4), were calculated in terms of the starting and ending values of GI in the regions A and B.

The variation of the viscosity of the polymerizable monomer phase, [[mu].sub.d], was calculated using the following non-ideal mixing equation for the monomer/polymer solution (Song et al., 2003):

ln [[mu].sub.d] = (1 - [[phi].sub.p]) ln [[mu].sub.m] + [[phi].sub.p] ln [[mu].sub.p] + [[k.sub.m,p] + [l.sub.m,p](1 - 2[[phi].sub.p])] (1 - [[phi].sub.p])[[phi].sub.p] ln [[mu].sub.m,p] (5)

ln [[mu].sub.m,p] = ln [[mu].sub.p] - ln [[mu].sub.m] / 2 (6)

where [[mu].sub.p] and [[mu].sub.m] denote the viscosity of the polymer and monomer, respectively, [[psi].sub.p] is the polymer volume fraction in the polymerizing monomer droplets and [k.sub.m,p] and [l.sub.m,p] are model parameters (Song et al., 2003) (see Table 1). Additional details, regarding the calculation of the densities of the suspension system, [[rho].sub.s], continuous phase, [[rho].sub.c], and dispersed phase, [[rho].sub.d], as well as the viscosity of the continuous phase, [[mu].sub.c], can be found in the publication of Kotoulas and Kiparissides (2006).

The interfacial tension between the aqueous and the monomer/polymer dispersed phase, [sigma], is a key parameter governing the suspension stability and, consequently, the droplet/particle size distribution developments. In the present study, the variation of the interfacial tension with monomer conversion was indirectly taken into account through the variation of [k.sub.c] (see Eq. 4). The variation of Q with respect to the stabilizer (PVA) concentration was calculated by the following equation (Adamson, 1976):

[sigma] = [[sigma].sub.0] - [K.sub.[sigma]] [K.sub.A] [C.sub.PVA] / 1 + [K.sub.A] [C.sub.PVA] (7)

where [K.sub.[sigma]] and [K.sub.A] are model parameters and [[sigma].sub.0] denotes the interfacial tension between the pure water phase and the dispersed monomer phase, in the absence of the stabilizer (see Table 1).

THE MONTE CARLO METHOD

The stochastic Monte Carlo (MC) method is based on the principle that the dynamic evolution of an extremely large population of particles/droplets, [N.sup.p](t), (e.g., [10.sup.20]) can be followed by tracking down the relevant events (i.e., coalescence and breakage) occurring in a smaller population of sample particles/droplets, [N.sub.s](t), (e.g., [10.sup.4]). It is commonly assumed that the effects of the droplet mechanisms on the properties of the sample population are representative of the corresponding effects on the properties of the total droplet population. The MC algorithm employed in the present study is schematically depicted in Figure 3. Initially, the size of the sample population, [V.sub.ms](0), is defined via the division of the initial monomer volume, [V.sub.mp](0), by a factor f. Subsequently, monomer droplets are randomly created to form the sample population in such a way so that the particle array at time zero closely represents the initial particle size distribution and the total volume of the generated droplets is equal to [V.sub.ms](0). This procedure involves the random assignment of droplet diameters, [D.sub.i], and respective volumes, [V.sub.i] ([V.sub.i] = [pi][N.sup.3.sub.i]/6), to every single droplet in the sample population. It should be noted that for every random procedure described in the present work, a random number generator from a uniform distribution in the range [0,1] was employed. In the present study, the initial droplet-diameter distribution was assumed to follow a normal distribution:

P([D.sub.i]) = (1 / [sigma][square root of 2[pi]]) exp (- [([D.sub.i] - [mu]).sup.2] / 2[[sigma].sub.2]), i = 1,2, ..., [N.sub.s] (8)

where the parameters [mu] and [sigma] denote the mean and the standard deviation of the distribution, respectively. The initial sample population was generated according to the approach of Box and Muller (1958). A typical sample population usually contains about [10.sup.3]-[10.sup.5] droplets in order to ensure an accurate representation of the initial distribution and, at the same time, to minimize the computational requirements.

Once the sample population has been formed, the MC algorithm is initiated and the effects of droplet coalescence and breakage mechanisms on the dynamic evolution of the droplet population are stochastically simulated in a series of consecutive variable-duration time steps. At the beginning of every time step, two consecutive decisions need to be made on the basis of the calculated droplet/particle mechanism rates, using a random number generator. Initially, the type of event (i.e., coalescence or breakage) that will take place, in the next infinitesimal time interval, is determined and, subsequently, a specific droplet or pair of droplets that will undergo droplet breakage or coalescence, are chosen from the sample population. The occurrence of a droplet/particle coalescence or breakage event is decided in accordance to the following probabilities:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (9)

where [R.sub.i] denotes the total rate ([s.sup.-1]) of each event (i.e., coalescence or breakage) in the sample population. For the selection of a specific droplet/particle (i.e., for breakage) or a pair of droplets/particles (i.e., for coalescence), the following acceptance-rejection procedure was employed (Garcia et al., 1987):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (10)

where [R.sup.*.sub.i] denotes the rate of droplet breakage (i=1) or droplet coalescence (i=2) and [rn.sub.j], is a randomly generated number. Based on extensive numerical simulations (Kotoulas, 2006), it was found that the maximum values of [g.sub.max](V) and [k.sub.max](V, U) could be related to the maximum and minimum values of [V.sub.max] and [V.sub.min]

[g.sub.max](V) = g([V.sub.max]), [k.sub.max](V, U) = k([V.sub.max], [V.sub.min]) (11)

[FIGURE 3 OMITTED]

Finally, the time interval, [DELTA]t, required for the occurrence of a discrete event was given by the reciprocal of the sum of the rates of the two events, [summation.sup.2.sub.i=1] [R.sub.i]:

[DELTA]t = 1/[2.summation over (i = 1)] [R.sub.i] (12)

The main drawback of the MC method is often the long simulation times. One of the new features of the proposed algorithm, that resulted in a significant decrease of the simulation time, was related to the calculation of the summation terms in Equations (9) and (12). Note that the commonly employed approach requires the calculation of a double summation term over the total number of droplets in the sample population:

[2.summation over (i = 1)][R.sub.i] = [[N.sub.s].summation over (j=1)] (g([V.sub.j]) + [[N.sub.s].summation over (l=j+1) k([V.sub.j], [V.sub.i])) (13)

In the present study, the following stochastic approach was employed for the above calculation that led to a significant decrease in the total simulation time:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (14)

The symbols [V.sub.l], [V.sub.m], and [V.sub.n] denote the volumes of the droplets that are randomly selected from the sample population within each summation and [q.sub.1] and [q.sub.2] are the percentage parameters. The identification of parameters [q.sub.1] and [q.sub.2] was the outcome of extensive numerical simulations. It was found that the values of these parameters had to vary in the range 5-25 for an accurate calculation of the coalescence and breakage rate terms. The above simplification resulted in a reduction of the computational effort approximately by 80-96%. Typical simulation times in a 2.21 GHz processor (AMD Athlon[TM] 64 X2, Dual Core Processor 4200 +) were of the order of 10-500 min, depending on the specific process conditions. Finally, it should be noted that the sample droplet population was kept within a range of predetermined limits, [[f.sub.A][N.sub.s](0), [f.sub.N][N.sub.s](0)], in order to preserve statistical accuracy of the solution and maintain the computational time at acceptable levels. This was accomplished by a refreshing procedure. Accordingly, the population of sample droplets/particles was re-instated, at frequent time intervals, to its initial value in such a way so that the information gathered from the precedent events was preserved. A more detailed description of the refreshing procedure can be found in a recent publication by Meimaroglou et al. (2006).

For comparison and graphical representation of the results, the droplet-diameter domain was discretized into a number of discrete elements (i.e., (D(i), D(i + 1)), i = 1,2, ..., [NE.sub.d]). The number diameter density function, n(D, t), was then calculated by the following equation:

[n.sub.i](D,t) = [N.sub.i](t) / D(i + 1) - D(i), i = 1,2, ..., [NE.sub.d] (15)

where [N.sub.i](t) is the number of sample droplets/particles in the size range [D(i), D(i + 1)] at time t. Thus, [n.sub.i] (D, t) dD denotes the total number of droplets/particles per unit volume, in the diameter range (D(i), D(i + 1)). Similarly, the respective volume probability diameter density function was calculated by the following equation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16)

THE FIXED PIVOT TECHNIQUE

In the present study, the FTP of Kumar and Ramkrishna (1996a) was employed for solving the continuous PBE (Eq. 1). According to the original developments of Kumar and Ramkrishna (1996a), the droplet/particle volume domain was initially discretized into a number of NE elements (i.e., ([V.sub.i], [V.sub.i+1]), i=1, 2, ..., [NE.sub.v] and the overall particle population (i.e., [N.sub.i], i= 1, 2, ..., [NE.sub.v]) was assigned to selected discrete points of the domain, also known as grid points (i.e., [x.sub.i], i = 1, 2, ..., [NE.sub.v]. It should be noted that breakage and coalescence events may lead to the formation of droplet/particle sizes that do not coincide with the selected grid points. In this case, the concentrations of the generated droplets/particles are assigned to two neighbouring grid points in such a way so that two selected properties of the population (i.e., total number, mass, etc., corresponding to the leading moments of the distribution) are exactly preserved:

a (V, [x.sub.i])[f.sub.1]([x.sub.i]) + b(V,[x.sub.i+1])[f.sub.1][x.sub.i+1] = [f.sub.1](V) (17)

a (V, [x.sub.i])[f.sub.2]([x.sub.i]) + b(V,[x.sub.i+1])[f.sub.2][x.sub.i+1] = [f.sub.2](V) (18)

Assuming that the number density function remains constant in the discrete volume interval ([V.sub.i] to [V.sub.i+1]), one can define a particle number distribution, [N.sub.i] (t), corresponding to the "ith" element of the discrete volume domain:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (19)

By integrating Equation (1) over the discrete volume interval ([V.sub.i], [V.sub.i+1]), the following set of discretized equations can be derived (Kotoulas and Kiparissides, 2006):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (20)

where [n.sub.b] ([x.sub.i], [x.sub.k]) denotes the fraction of droplets/particles of size [x.sub.i], resulting from the breakage of a droplet/particle of size [x.sub.k].

To preserve the number and the volume of particles in the size range ([V.sub.i], [V.sub.i+1]), the fraction [n.sub.b] ([x.sub.i], [x.sub.k]) must satisfy the following equation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (21)

The respective expression for [n.sub.c] ([x.sub.i], V), accounting for the fraction of droplets/particles of size V (=[x.sub.j] + [x.sub.k]) formed via the coalescence of two droplets/particles of volumes [x.sub.j] and [x.sub.k], will be given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (22)

Accordingly, from the calculated values of [N.sub.i] (t), one can easily calculate the values of the average number density function, [[bar.n].sub.i] (D, t):

[[bar.n].sub.i] (V, t) = [N.sub.i] (t) / ([V.sub.i+1] - [V.sub.i]) (23)

It is often desirable to know the number diameter density function, n (D, t). By noting that n (V, t) dV = n (D, t) dD, one can easily calculate the average number density function, [[bar.n].sub.i] (D, t), in terms of [N.sub.i] (t):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (24)

In many cases, experimental measurements on PSD are given in terms of number or volume fractions, from which it is not easy to derive the actual particle number distribution, [N.sub.i] (t), or the number volume density function, n(V, t). Thus, it is preferable to directly compare available experimental measurements on number and volume fraction distributions with respective simulation results obtained from the numerical solution of the PBE. Accordingly, one can define the following number, A(D, t), and volume, [A.sub.v] (D, t), probability density functions:

A(D,t) = n(D,t) / [N.sub.t](t) [approximately equal to] f[N.sub.i] / [D.sub.i+1] - [D.sub.i] (25)

[A.sub.v](D,t) = ([pi][D.sup.3] / 6)n(D,t) / [V.sub.t](t) [approximately equal to] f[V.sub.i]/[D.sub.i+1] - [D.sub.i] (26)

where A (D, t) dD and [A.sub.v] (D, t) dD represent the number, f[N.sub.i], and volume, f[V.sub.i], fractions of particles in the size range (D, D + dD), respectively. [N.sub.t] (t) and [V.sub.t] (t) denote the respective total number and volume of particles per unit volume of the reaction medium. It is apparent that the number and volume probability density functions will satisfy the following normalization conditions:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (27)

Very often, experimental measurements on some average particle diameters are only available. In general, the average particle diameter, [D.sub.qp], can be calculated, in terms of the number probability density function, using the following equation (Chatzi and Kiparissides, 1992):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (28)

where q and p are characteristic exponents that define the desired average particle diameter (e.g., mean Sauter diameter, [D.sub.32], average number diameter, [D.sub.10], average volume diameter, [D.sub.30], etc.)

RESULTS AND DISCUSSION

Extensive simulations were carried out, using the two numerical methods (i.e., FPT and MC), for both non-reactive liquid-liquid dispersions (i.e., methyl methacrylate (MMA)-water) and reactive (i.e., MMA free-radical suspension polymerization) systems. In order to validate the MC and FPT numerical results, model predictions were directly compared with experimental data reported in the open literature (Jahanzad et al., 2005). The experiments were carried out in a 1 L jacketed glass reactor (internal diameter 10 cm) equipped with four baffles. A four-bladed flat turbine impeller with a diameter of 5 cm was used for agitation. The polymerization temperature was controlled within [+ or -] 0.5[degrees]C of the specified temperature. More details, regarding the realization of the experiments can be found in the original work of Jahanzad et al. (2005). The PSD model developed in the present work was solved together with a comprehensive dynamic kinetic model, describing the polymerization rate and monomer conversion of MMA in a batch suspension polymerization reactor. A detailed description of the kinetic model can be found in a recent publication by Meimaroglou et al. (2007). The physical and transport parameters, appearing in the PSD model of the MMA suspension polymerization system, are reported in Table 1.

In Figures 4 and 5, the effect of agitation speed (i.e., 220, 300, and 500rpm) on the Sauter mean diameter, [D.sub.32], is depicted for a non-reactive liquid-liquid dispersion (i.e., MMA-water) and the free-radical suspension polymerization of MMA, respectively. It is evident that both numerical methods display a very good agreement with the experimental measurements for both reactive and non-reactive systems. In both cases, the initial liquid-liquid dispersion undergoes significant droplet breakage, leading to a decrease in the mean size of the droplets. In reactive systems (i.e., suspension polymerization), the initial phase is followed by a second sticky-phase where the viscosity of the dispersed phase increases and particle coalescence becomes the dominant mechanism, resulting in a large increase in the average particle size. At the final phase, the effects of the two mechanisms are diminished, leading to an almost time-invariant PSD. Note that an increase in the agitation speed results in a subsequent increase of the droplet breakage rate (see Eq. 33), leading to an increase in the number of small droplets.

In Figure G, the time evolution of the DSD, in terms of the volume density function, is depicted for a non-reactive liquid-liquid dispersion at 500 rpm. Apparently, there is a very good agreement between the DSDs, calculated by the MC and FP methods, and the experimentally measured distributions at different times. It can also be seen that the DSD becomes narrower as the agitation time increases. At the same time, the maximum value of the distribution moves to smaller droplet sizes, an observation consistent with the results of Figure 4.

[FIGURE 4 OMITTED]

[FIGURE 5 OMITTED]

In Figure 7, experimentally measured PSDs are compared with model predictions obtained by the two numerical methods (i.e., FPT and MC), for the free-radical suspension polymerization of MMA, at different times (i.e., 1, 20, and 70 min). As can be seen, the simulation results are in very good agreement with the experimental measurements. The initial liquid-liquid dispersion undergoes significant droplet breakage that shifts the distribution to smaller droplet sizes. However, in the following second sticky-phase, the particle coalescence mechanism becomes dominant and thus, the distribution moves to larger particle sizes.

Figure 8 illustrates the effect of the surfactant concentration on the Sauter mean diameter, [D.sub.32]. Experimental measurements on [D.sub.32] are compared with model predictions obtained by the two numerical methods (i.e., FPT and MC) for four different surfactant concentrations (i.e., [C.sub.PVA] = 0.5, 1,4 and 10 g/L). Both non-reactive liquid-liquid (lower curves) and reactive (upper curves) dispersions are considered. It is apparent that both numerical methods are in very good agreement with the experimental results. Furthermore, it can be easily observed that droplets having a higher stabilizer coverage are less susceptible to coalescence. This results in a decrease in the mean particle size as the PVA concentration increases.

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

It should be noted that in most cases studied in the present work, the FPT exhibited an instantaneous "undershoot" in the prediction of [D.sub.32] at the onset of the gel-effect. This can be attributed to the inherent difficulty associated with the numerical solution of the stiff system of DAEs, generated by the FPT discretization. Note that, at this specific point, the stiffness of the numerical solution considerably increases due to the sharp increase in the polymerization rate (i.e., large change in the particle growth rate) due to the appearance of the gel-effect.

In Figures 9 and 10, model predictions and experimental measurements are depicted for both non-reactive (lower curves) and reactive (upper curves) dispersions, for two initial monomer volume fractions (i.e., ([[PHI].sub.m] = 0.05 % and 0.2 %), respectively. As can be seen, as the monomer volume fraction increases, the final mean droplet/particle size increases due to the higher concentration of monomer droplets and thus, the higher droplet collision frequency.

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

In Figures 11 and 12, the effect of temperature on the Sauter mean diameter, [D.sub.32] is shown. Experimental measurements on [D.sub.32] are compared with model predictions obtained by the two numerical methods (i.e., FPT and MC) for both non-reactive (lower curves) and reactive (upper curves) dispersions. Note that, as the polymerization temperature increases (i.e., polymerization rate increases), the duration of the droplet breakage phase decreases (i.e., the polymerization moves to the second sticky-phase at an earlier time). As a result, larger droplet/particles are produced at higher polymerization temperatures due to the shorter breakage time.

[FIGURE 10 METHOD]

Finally, in Figure 13, the calculated PSDs by the two numerical methods are depicted at various polymerization times. It is apparent that both methods display an excellent agreement throughout the polymerization process.

CONCLUSIONS

In the present work, a comprehensive population balance model, coupled with a system of differential equations governing the conservation of the various molecular species present in the system, was employed for the calculation of the DSD/PSD in non-reactive and reactive dispersion systems. Two different numerical methods were developed for the solution of the governing PBEs, namely, the Fixed Pivot technique and an efficient stochastic Monte Carlo method.

[FIGURE 11 OMITTED]

[FIGURE 12 OMITTED]

[FIGURE 13 OMITTED]

Numerical simulations were carried out under various process conditions (i.e., agitation speed, temperature, surfactant concentration, and initial monomer volume fraction). The predictive capabilities of the model and the validity of the two numerical methods were demonstrated via the successful simulation of experimental measurements on the DSD/PSD and the average droplet/particle diameter for both non-reactive liquid-liquid dispersions and the free-radical suspension polymerization of MMA. The numerical efficiency and accuracy of the two numerical methods were tested via numerous simulations. The FP numerical method was found to be computationally less demanding (i.e., <30 min in a 2.4 GHz processor, Intel[R] Core[TM] 2 Duo Processor E6600) but its accuracy was highly dependent on the discretization of the particle size domain. On the other hand, the MC algorithm proved to be very easy to implement but its accuracy was highly dependent on the size of the sample particle population. Typical MC simulation times in a 2.21 GHz processor (AMD Athlon[TM] 64 X2, Dual Core Processor 4200+) were of the order of 10-500 min, depending on the specific process conditions.

NOMENCLATURE [a.sub.b],[a.sub.c] model parameters A surface area ([m.sup.2]) A (D, t), [A.sub.v] (D, t) number and volume probability density functions ([m.sup.-1]) [C.sub.PVA] concentration of surface-active agent (kg/[m.sup.3]) D diameter (m) [D.sub.10] average number diameter (m) [D.sub.30] average volume diameter (m) [D.sub.32] mean Sauter diameter (m) [E.sub.p] polymer elasticity modulus (kg/m [s.sup.2]) f sampling factor in the MC simulation [f.sub.A] parameter for sample particle introduction in the MC simulation [f.sub.N] parameter for sample particle removal in the MC simulation g(V) breakage rate ([s.sup.-1]) h heat-transfer coefficient between the particle and the continuous phase (J/s [m.sup.2] K) k(V, U) coalescence rate ([m.sup.3]/s) [MATHEMATICAL EXPRESSION model parameters NOT REPRODUCIBLE IN ASCII] [k.sub.f] the thermal conductivity of the continuous phase (J/s m K) [K.sub.[sigma]] [K.sub.A] model parameters L macroscale of turbulence (m) [n] intrinsic viscosity ([m.sup.3] /kg) NE total number of discrete elements in the droplet/particle size domain [N.sub.i] number of particles having volume equal to [x.sub.i] per reactor unit volume ([m.sup.-3]) [N.sub.p] (t) number of droplets in the MC simulation total population at time t. [N.sub.s] (t) number of droplets in the MC simulation sample population at time t. n(V, t) number density function ([m.sup.-6]) P event occurrence probability in MC simulation Pr Prandtl number R total rate of event in MC simulation ([s.sup.-1]) r volume ratio of daughter over the satellite drops R* individual rate of event in MC simulation, for given droplet/pair of droplets ([s.sup.-1]) Re Reynolds number [R.sup.max] maximum individual rate of event in MC simulation ([s.sup.-1]) [MATHEMATICAL EXPRESSION polymerization rate at degree of NOT REPRODUCIBLE IN AS conversion equal to zero (mol/[m.sup.3] s) [R.sub.p] polymerization rate (mol/[m.sup.3] s) [S.sub.Nsa] model parameter T temperature (K) t time (s) [??][([D.sub.u]).sup.2] mean square of the relative velocity between two points separated by a distance D (m/s) u (V) number of droplets formed by a breakage of a droplet of volume V V, U, x volumes ([m.sup.3]) [V.sub.mp] (t) volume of droplets in the MC simulation sample population at time t ([m.sup.3]) We Weber number Greek Symbols [[alpha].sub.b][[alpha].sub.c] model parameters [beta](V, U) daughter droplets probability function (1/[m.sup.3]) [DELTA][H.sub.r] heat of polymerization (J/mol) [DELTA]t time interval in MC simulation (s) [epsilon] average energy dissipation rate per unit mass ([m.sup.2]/[s.sup.3]) [eta] microscale of turbulence (m) [[lambda].sub.b] [[lambda].sub.c] breakage and coalescence efficiencies [mu] viscosity (kg/m s) [rho] density (kg/[m.sup.3]) [sigma] interfacial tension (kg/[s.sup.2]) [[sigma].sub.0] interfacial tension between the pure water and the dispersed phase (kg/[s.sup.2]) [psi] volume fraction x monomer conversion [[omega].sub.b], [[omega].sub.c] breakage and coalescence frequencies ([s.sup.1]) Subscripts c continuous phase d dispersed phase m monomer p polymer s suspension system w water Abbreviations DAE differential-algebraic equation DPBE discretized population balance equation DSD droplet size distribution FPT fixed pivot technique MC Monte Carlo MMA methyl methacrylate OCFE orthogonal collocation on finite elements PBE population balance equation PSD particle size distribution PVA polyvinyl alcohol)

APPENDIX: CALCULATION OF BREAKAGE AND COALESCENCE RATES

According to the original work of Alvarez et al. (1994) and the proposed modifications of Maggioris et al. (2000), the breakage and coalescence rates can be expressed in terms of the breakage, [[omega].sub.b], and collision, [[omega].sub.c], frequencies and the respective, Maxwellian efficiencies, [[lambda].sub.b] and [[lambda].sub.d]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (29)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (30)

where [[lambda].sub.b] and [[lambda].sub.c] denote the corresponding ratios of required to available energy for an event to occur.

In the present study, the breakage of a drop exposed to a turbulent flow field was supposed to occur as result of energy transfer from an eddy to a drop having a diameter equal to the eddy wave length, [D.sub.u]. Eddy fluctuations with wavelengths smaller (larger) than the drop diameter [D.sub.u] produce an oscillatory (rigid body) motion of the drop that does not lead to breakage (Alvarez et al., 1994). The frequency term, [[omega].sub.b] (V), was assumed to be equal to the inverse fluctuation time period, corresponding to the time required for a drop to reach its mean drop displacement:

[[omega].sub.b](V) = [??][D.sub.u] / [D.sub.u] (31)

where [??][([D.sub.u]).sup.2] is the mean square of the relative velocity between two points separated by a distance [D.sub.u], or the mean square fluctuation velocity of drops of diameter [D.sub.u].

For drops in the inertial subrange of turbulence (i.e., [eta] < [D.sub.u] [much less than] L), the energy spectrum will be independent of the kinematic viscosity, [v.sub.c], and the mean fluctuation velocity will be solely determined by the rate of energy dissipation (Hinze, 1959):

[??][([D.sub.u]).sup.2] = [k.sub.b] [([[epsilon].sub.s][D.sub.u]).sup.2/3] (32)

where [k.sub.b] is a model parameter and [[epsilon].sub.s] is the average energy dissipation rate for the dispersion.

For droplets in the viscous dissipation range (i.e., [D.sub.u] < [eta]), the inertial forces are of the same order of magnitude as the viscous shear forces and the mean square of the relative velocity between two points separated by a distance [D.sub.u], will be given by (Shinnar, 1961):

[??][([D.sub.u]).sup.2] = [k.sub.b] [D.sup.2.sub.u] ([[epsilon].sub.s] / [v.sub.c]) (33)

According to Alvarez et al. (1994), for an effective drop breakage to occur, the drop surface energy and drop viscoelastic resistance must be overcome. Considering that the flow within a drop can be described as one-dimensional simple-shear flow of a Maxwell fluid, the breakage efficiency can be expressed as follows:

[[lambda].sub.b] = [[alpha].sub.b][OMEGA]([D.sub.u]) (34)

where [[alpha].sub.b] is a model parameter and [OMEGA] ([D.sub.u]) is the ratio of required to available energy for a drop of diameter [D.sub.u] to break:

[OMEGA]([D.sub.u]) = 2 / Re (1 + ReVe) + [C.sub.ds] / We (35)

where Re and We denote the droplet Reynolds and Weber numbers, respectively. The dimensionless quantity Ve accounts for drop viscoelastisity whereas the scalar quantity [C.sub.ds] can be expressed in terms of the numbers and volumes of daughter and satellite drops (Chatzi and Kiparissides, 1992).

Assuming that the droplet collision mechanism in a local isotropic flow field is analogous to collisions between molecules as in the kinetic theory of gases, the collision frequency between two drops with volumes V and U can be expressed as (Maggioris et al., 2000):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (36)

In general, the monomer drops will behave like deformable drops at the beginning of polymerization while, at high monomer conversions, they will behave like rigid polymer particles. Thus, the coalescence efficiency over the whole monomer conversion range can be written as:

exp{-[[lambda].sub.c](V,U)}=(1 - x)exp{-[[lambda].sup.(a).sub.c](V,U)} + x exp{-[[lambda].sup.(b).sub.c](V,U)} (37)

where [chi] is the fractional monomer conversion. The terms [[lambda].sup.(a).sub.c] and [[lambda].sup.(b).sub.c] can be calculated from the following expressions for deformable drops and rigid spheres, respectively (Kotoulas and Kiparissides, 2006):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (38)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (39)

where [[alpha].sub.c] is a model parameter, [[mu].sub.c] and [[rho].sub.c] are the viscosity and density of the continuous phase, respectively, and [sigma] is the interfacial tension between the dispersed and aqueous phases, respectively. More details regarding the breakage and coalescence rates can be found in the publications of Maggioris et al. (2000) and Kotoulas and Kiparissides (2006).

Manuscript received March 7, 2008; revised manuscript received May 30, 2008; accepted for publication June 4, 2008.

REFERENCES

Achilias, D. S. and C. Kiparissides, "Development of a General Mathematical Framework for Modeling Diffusion-Controlled Free-Radical Polymerization Reactions," Macromolecules 25(14), 3739-3750 (1992).

Adamson, A. W., "Physical Chemistry of Surfaces," Wiley, New York (1976).

Alexopoulos, A. I. and C. Kiparissides, "Part II: Dynamic Evolution of the Particle Size Distribution in Particulate Processes Undergoing Simultaneous Particle Nucleation, Growth and Aggregation," Chem. Eng. Sci. 60, 4157-4169 (2005).

Alexopoulos, A. H., A. I. Roussos and C. Kiparissides, "Part I: Dynamic Evolution of the Particle Size Distribution in Particulate Processes Undergoing Combined Particle Growth and Aggregation," Chem. Eng. Sci. 59, 5751-5769 (2004).

Alvarez, J., J. Alvarez and M. Hemansez, "A Population Balance Approach for the Description of the Particle Size Distribution in Suspension Polymerization Reactors," Chem. Eng. Sci. 49, 99-113 (1994).

Batterham, R. J., J. S. Hall and G. Barton, "Pelletizing Kinetics and Simulation for Full-Scale Balling Circuits," In "Proc. 3rd Int. Symp. Agg., Nur.," West Germany, A136 (1981).

Bleck, R., "A Fast, Approximate Method for Integrating the Stochastic Coalescence Equation," J. Geoph. Res. 75, 5165-5171 (1970).

Box, G. E. P. and M. Muller, "A Note on the Generation of Random Normal Deviates," Ann. Math. Stat. 29, 610-611 (1958).

Brandrup, J. and E. H. Immergut, "Polymer Handbook," 3rd ed. John Wiley & Sons, New York (1989).

Chatzi, E. G. and C. Kiparissides, "Dynamic Simulation of Bimodal Drop Size Distributions in Low-Coalescence Batch Dispersion Systems," Chem. Eng. Sci. 47, 445-456 (1992).

Chatzi, E. G., A. D. Gavrielides and C. Kiparissides, "Generalized Model for Prediction of the Steady-State Drop Size Distribution in Batch Stirred Vessel," Ind. Eng. Chem. Res. 28, 1704-1711 (1989).

Chen, M. Q., C. Hwang and Y P. Shih, "A Wavelet-Galerkin Method for Solving Population Balance Equations," Comp. Chem. Eng. 20(2), 131-145 (1996).

Coulaloglou, C. A. and L. L. Tavlarides, "Description of Interaction Processes in Agitated Liquid-Liquid Dispersions," Chem. Eng. Sci. 32, 1289-1297 (1977).

Dafniotis, P., "Modeling of Emulsion Copolymerization Reactors Operating Below the Critical Micelle Concentration," Ph.D. Thesis, University of Wisconsin, Madison, WI (1996).

Debry, E., B. Sportisse and B. Jourdain, "A Stochastic Approach for the Numerical Simulation of the General Dynamics Equation for Aerosols," J. Comp. Phys. 184, 649-669 (2003).

Domilovskii, E. R., A. A. Lushnikov and V N. Piskunov, "A Monte Carlo Simulation of Coagulation Processes," Izvestiya Akademii Nauk SSSR, Fizika Atmsfery I Okeana 15, 194-201 (1979).

Efendiev, Y and M. R. Zachariah, "Hydrid Monte Carlo Method for Simulation of Two-Component Aerosol Coagulation and Phase Segregation," J. Coll. Interf. Sci. 249, 30-43 (2002).

Garcia, A. L., C. van de Broek, M. Aertsens and R. Serneels, "A Monte Carlo Simulation of Coagulation," Physica 143 (3), 535-546 (1987).

Gelbard, F. and J. H. Seinfeld, "Exact Solution of the General Dynamic Equation for Aerosol Growth by Condensation," J. Coll. Interf. Sci. 68(1), 173-183 (1979).

Gelbard, F. and J. H. Seinfeld, "Simulation of Multicomponent Aerosol Dynamics," J. Coll. Interf. Sci. 78(2), 485-501 (1980).

Hamielec, A. E. and H. Tobita, "Polymerization processes," In: "Ulmann's Encycl. Ind. Chem.," Vol. A21. CVH Publishers, New York (1992), pp. 305-428.

Hidy, G. M., "On the Theory of the Coagulation of Noninteracting Particles in Brownian Motion," J. Coll. Sci. 20, 123-144 (1965).

Hinze, J. O., "Turbulence," McGraw-Hill, New York (1959).

Hounslow, M. J., R. L. Ryall and V R. Marshall, "Discretized Population Balance for Nucleation, Growth, and Aggregation," AIChE J. 34(11), 1821-1832 (1988).

Jahanzad, F., S. Sajjadi and B. W Brooks, "Comparative Study of Particle Size in Suspension Polymerization and Corresponding Monomer-Water Dispersion," Ind. Eng. Chem. Res. 44, 4112-4119 (2005).

Kiparissides, C., "Polymerization Reactor Modeling: A Review of Recent Developments and Future Directions," Chem. Eng. Sci. 51, 1637-1659 (1996).

Kiparissides, C., A. Alexopoulos, A. Roussos, G. Dompazis and C. Kotoulas, "Population Balance Modeling of Particulate Polymerization Processes," Ind. Eng. Chem. Res. 43, 7290-7302 (2004).

Kostoglou, M. and A. J. Karabelas, "Evaluation of Zero Order Methods for Simulating Particle Coagulation," J. Coll. Interf. Sci. 163, 420-431 (1994).

Kostoglou, M. and A. J. Karabelas, "Evaluation of Numerical Methods for Simulating an Evolving Particle Size Distribution in Growth Processes," Chem. Eng. Commun. 136, 177-199 (1995).

Kotoulas, C., "Prediction of the Particle Size Distribution in Suspension Polymerization Reactions," PhD Thesis, Aristotle University of Thessaloniki (2006).

Kotoulas, C. and C. Kiparissides, "A Generalized Population Balance Model for the Prediction of Particle Size Distribution in Suspension Polymerization Reactors," Chem. Eng. Sci. 61, 332-346 (2006).

Kruis, F. E., A. Maisels and H. Fissan, "Direct Simulation Monte Carlo Method for Particle Coagulation and Aggregation," AIChE J. 46(9), 1735-1742 (2000).

Kumar, S. and D. Ramkrishna, "On the Solution of Population Balance Equations by Discretizations. I. A Fixed Pivot Technique," Chem. Eng. Sci. 51(8), 1311-1332 (1996a).

Kumar, S. and D. Ramkrishna, "On the Solution of Population Balance Equations by Discretizations. II. A Moving Pivot Technique," Chem. Eng. Sci. 51(8), 1333-1342 (1996b).

Landgrebe, J. D. and S. E. Pratsinis, "A Discrete-Sectional Model for Powder Production by Gas-Phase Chemical Reaction and Aerosol Coagulation in the Free-Molecular Regime," J. Coll. Interf. Sci. 139(1), 63-86 (1990).

Liffman, K., "A Direct Simulation Monte Carlo Method for Cluster Coagulation," J. Comp. Phys. 100, 116-127 (1992).

Litster, J. D., D. J. Smit and M. J. Hounslow, "Adjustable Discretized Population Balance for Growth and Aggregation," AIChE J. 41, 591-603 (1995).

Maggioris, D., A. Goulas, A. H. Alexopoulos, E. G. Chatzi and C. Kiparissides, "Prediction of Particle Size Distribution in Suspension Polymerization Reactors: Effect of Turbulence Nonhomogeneity," Chem. Eng. Sci. 55, 4611-4627 (2000).

Maisels, A., F. E. Kruis and H. Fissan, "Direct Simulation Monte Carlo for Simultaneous Nucleation, Coagulation, and Surface Growth in Dispersed Systems," Chem. Eng. Sci. 59, 2231-2239 (2004).

Marchal, P., R. David, J. P. Klein and J. Villermaux, "Crystallization and Precipitation Engineerings. I. An Efficient Method for Solving Population Balance in Crystallization with Agglomeration," Chem. Eng. Sci. 43(1), 59-67 (1988).

Meimaroglou, D., A. I. Roussos and C. Kiparissides, "Part IV: Dynamic Evolution of the Particle Size Distribution in Particulate Processes. A Comparative Study between Monte Carlo and the Generalized Method of Moments," Chem. Eng. Sci. 61, 5620-5635 (2006).

Meimaroglou, D., A. Krallis, V Saliakas and C. Kiparissides, "Prediction of the Bivariate Molecular Weight--Long Chain Branching Distribution in Highly Branched Polymerization Systems Using Monte Carlo and Sectional Grid Methods," Macromolecules 40, 2224-2234 (2007).

Narsimhan, G., G. Gupta and D. Ramkrishna, "A Model for Translational Breakage Probability of Droplets in Agitated Lean Liquid-Liquid Dispersions", Chem. Eng. Sci. 34, 257-265 (1979).

Nicmanis, M. and M. J. Hounslow, "A Finite Element Analysis of the Steady-State Population Balance Equation for Particulate Systems: Aggregation and Growth," Chem. Eng. Sci. 20, S261-S266 (1996).

Nicmanis, M. and M. J. Hounslow, "Finite-Element Methods for Steady-State Population Balance Equations," AIChE J. 44, 2258-2272 (1998).

Ramkrishna, D., "Analysis of Population Balance IV The Precise Connection Between Monte Carlo and Population Balances," Chem. Eng. Sci. 36, 1203-1209 (1981).

Ramkrishna, D., "The Status of Population Balances," Rev. Chem. Eng. 3(1), 49-95 (1985).

Ranz, W E. and W R. Marshall, "Evaporation from Drops," Chem. Eng. Prog. 48, 173-180 (1952).

Roussos, A. I., A. H. Alexopoulos and C. Kiparissides, "Part III: Dynamic Evolution of the Particle Size Distribution in Batch and Continuous Particulate Processes: A Galerkin on Finite Elements Approach," Chem. Eng. Sci. 60, 6998-7010 (2005).

Sastry, K. V. S. and P. Gaschignard, "Discretization Procedure for the Coalescence Equation of Particulate Processes," Ind. Eng. Chem. Fund. 20, 355-361 (1981).

Shah, B. H. J. D. Borwanker and D. Ramkrishna, "Simulation of Particulate Systems Using the Concept of the Interval of Quiescence," AIChE J. 23, 897-904 (1977).

Shinnar, R., "On the Behaviour of Liquid Dispersions in Mixing Vessels," J. Fluid Mech. 10, 259-275 (1961).

Smith, M. and T. Matsoukas, "Constant-Number Monte Carlo Simulation of Population Balances," Chem. Eng. Sci. 53 (9), 1777-1786 (1998).

Song, Y., M. Mathias, D. Tremblay and C. C. Chen, "Liquid Viscosity Model for Polymer Solutions and Mixtures," Ind. Eng. Chem. Res. 42, 2415-2422 (2003).

Sovova, H., "Breakage and Coalescence of Drops in a Batch Stirred Vessel. II. Comparison of Model and Experiments," Chem. Eng. Sci. 36, 1567-1573 (1981).

Spielman, L. A. and O. Levenspiel, "A Monte Carlo Treatment of Reacting and Coalescing Dispersed Phase Systems," Chem. Eng. Sci. 20, 247-254 (1965).

Tandon, P. and D. E. Rosner, "Monte Carlo Simulation of Particle Aggregation and Simultaneous Restructuring," J. Coll. Interf. Sci. 213, 273-286 (1999).

Yuan, H. G., G. Kalfas and W H. Ray, "Suspension Polymerization," JMS-Rev Macromol. Chem. Phys. C31, 215-299 (1991).

Zhao, H., A. Maisels, T. Matsoukas and C. Zheng, "Analysis of Four Monte Carlo Methods for the Solution of Population Balances in Dispersed Systems," Powder Techn. 173, 38-50 (2007).

Vassilis Saliakas, Costas Kotoulas, Dimitris Meimaroglou and Costas Kiparissides * Department of Chemical Engineering, Aristotle University of Thessaloniki and Chemical Process Engineering Research Institute, P.O. Box 472, Thessaloniki, Greece 54124

* Author to whom correspondence may be addressed. E-mail address: cypress@alexandros.cperi.certh.gr

Table 1. Physical, transport properties and PSD model parameters for MMA/PMMA system [[rho].sub.m] = 968-1.225 (T - 273.15) (kg/[m.sup.3]) [[rho].sub.m] = 1212-84.5 (T - 273.15) (kg/[m.sup.3]) [[micro].sub.m] = [[10.sup.[453.25(1/T-1/254.92)]].sup.-3] (kg/m s) [n] = 9.8 x [10.sup.-5] [M.sup.0.64.sub.w] ([m.sup.3]/kg) [E.sub.p] = 3.3 x [10.sup.11] (kg/m [s.sup.2]) r = 35, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], [k.sub.b] = 40, [a.sub.b] = 41, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] = 6 x 10-7, [a.sub.c] = 1 x [10.sup.9], [a'.sub.c] = 3 x [10.sup.3], [k.sub.m,p] = 0.8, [l.sub.m,p] = -3, [K.sub.[sigma]] = 0.0104, [K.sub.A] = 10.56 Achilias and Kiparissides Achilias and Kiparissides (1992) Achilias and Kiparissides (1992) Achilias and Kiparissides (1992) Achilias and Kiparissides (1992) Brandrup and Immergut (1989) This study