Claude Universite

Bernard -Lyon

M Le Professeur

Frédéric Fleury

Hamda Ben

Hadid M Le Professeur Didier

Revel M Le Professeur

Philippe Chevalier

M Fabrice

Vallée Mme

Dominique Marchand

Composantes Sante

G Rode

La Mme

C Professeure

Burillon

D Bourgeois

Vinciguerra

M X Perrot

A-M Professeure

Schott

Un système consistant d'une simulation de mécanique des fluides numérique (MFN) couplée à un modèle de bilan de population (PBM) est développé afin d'étudier l'effet des paramètres variés sur la performance d'un procédé de polymérisation en émulsion qui conduit à la production des particules de polymère dans un milieu aqueux continu.

Comme une grande gamme des produits polymériques, des latexes sont les « produits par processus » (products-by-process), et leurs propriétés sont déterminés pendant la polymérisation. LA distribution de la taille des particules (PSD) est une des plus importants paramètres qui influence la qualité finale de latex. La modélisation d'évolution du PSD est généralement réalisée par l'addition un ensemble des équations de bilan de population (PBEs) au modèle cinétique. Le PBE fournit un moyen d'étudier la contribution des différents phénomènes dans l'évolution du PSD, comme la nucléation, la croissance des particules par la polymérisation, et la coagulation des particules à cause du mouvement brownien ou le mouvement du fluide (la coagulation Perikinetic et Orthokinetic, respectivement). Afin d'évaluer l'impact du mélange non homogène et les paramètres physiques du système sur l'évolution du PSD du latex, la simulation transitoire d'écoulement à été réalisé avec l'aide d'un progiciel commercial de MFN (Fluent® 15.0) pour munir dans chaque pas du temps, les concentrations locales des espèces ioniques (pour déterminer le taux de la coagulation Perikinetic modelé par le modèle de DLVO) ainsi que certains paramètres hydrodynamiques comme le taux de dissipation de la turbulence et le taux de cisaillement (afin de déterminer le taux de la coagulation Orthokinetic).

COMPOSANTES ET DEPARTEMENTS DE SCIENCES ET TECHNOLOGIE

To assess the impact of nonhomogeneous mixing and physical parameters of the system on the evolution of the latex PSD, the transient simulation of flow was performed with the aid of a commercial CFD Package (Fluent® 15.0) to provide in each time step, the local concentrations of ionic species (to determine the rate of perikinetic coagulation modeled by DLVO model) and certain hydrodynamic parameters such as turbulence dissipation rate and shear rate (to determine the rate of orthokinetic coagulation). This information is applied simultaneously by the PBE add-on module of Fluent to calculate the PSD for the next time step; thus, a complete coupling between CFD and PBM is assured. Maximum and average values of and in both scales of 4 L and 4000 L obtained from simulations as the tank geometry is scaled from 4 L to 4000 L, maintaining a constant volumetric holding power input (P/V ratio). is set to 500 RPM at the 4 L scale and to 108 RPM at the 4000 L scale

List of Figures

Greek Letters

coagulation rate coefficient between particles i and j (m 3 part -1 s -1 ) coagulation rate coefficient between particles i and j in the absence of particle interaction (m 3 part -1 s -1 ) rate of turbulent energy dissipation (m 2 s -

Introduction

Synthetic polymers can be denoted as the materials of the twentieth century. The enormous growth of synthetic polymers is due to the vast range of their applications.

Polymer latexes are an important class of this kind of materials consisting of colloidallystable dispersions of nanoscale polymer particles dispersed in a continuous medium.

The typical polymer latex recipe is comprised of monomer, water, surfactant, initiator and buffer and most polymer latex products are produced via emulsion polymerization process (EP). The need for large amounts of synthetic rubber that arose as a result of World War II was the motivation to develop the emulsion polymerization process technique. This technique was found to offer a number of advantages over bulk and solution polymerization such as the improvement of product properties, higher rates of productivity and controllability. In general, there are three major phenomena participating in the final latex properties: nucleation, growth and coagulation of particles. Careful manipulation of each of these steps will result in tailoring the latex properties in order to meet the rapid growing applications of this important class of materials.

The particle size distribution (PSD) is one of the most important characteristics of a polymer latex determining its rheological properties, maximum solid content, adhesion, drying time, etc. This is why having a fine control on this property is one of the most necessary aspects of its manufacturing process. A means of tracking changes to the particle size distribution during the course of polymerization is necessary when modeling EP processes. This will be done through the construction of a framework consisting a set of population balance model (PBM) added to the polymerization kinetic model. Using this framework, researchers have developed new manufacturing strategies that offer improved control over the particle size distribution.

Changing certain reaction parameters such as characteristic mixing times (processrelated parameter) with respect to characteristic reaction or coagulation times (chemistry/physico chemistry) will likely have an impact on the PSD. Considering the importance of understanding the impact of changing characteristic times in chemical processes, Computational Fluid Dynamics CFD simulations can serve perfectly. CFD simulations can provide information about flow characteristics that are very difficult to obtain from experiments alone. They can be used to generate detailed flow fields for a wide range of reactor configurations and operating conditions and in particular, quantify the mixing in these different configurations. When the information abstracted from CFD analysis (such as mixing characteristics, reactant gradients, or the shear rate distribution) are incorporated into a framework containing a descriptive process model, a more precise prediction of the product quality can be achieved giving a tool to evaluate the influence of reactor scale and configuration on the evolution of latex properties over the course of a reaction.

Speaking of the construction of a framework, a number of works have been done to couple CFD and PBE in order to predict some features such as PSD in different systems.

Due to the high computational requirements of the simultaneous coupling, most of these works employed a simplified sequential approach, where the CFD simulations are first realized then used in a PBM using a different software such as Matlab or Fortran.

The objective of this work is to develop a simultaneous two-way coupled CFD/PBM computational framework to efficiently predict the particle size distribution of a specific polymer latex called poly(vinylidene fluoride) PVDF. This fluoropolymer is widely employed in industrial applications demanding materials with high performance thermoplastic characteristics, such as architectural coating industry, manufacturing of fitting, valves and pumps in chemical industry, and as insulation material in wire and cable industry. In industrial scale PVDF is produced by suspension and emulsion polymerization which will result in different products in terms of properties. The choice of the polymerization method depends on the area of application. The focus of this work is the EP process of vinylidene fluoride VDF. The framework is used to investigate the impact of changes in key process length and time scales on the properties of latex products. This was accomplished by using the "Discrete Method" as the PBM solution method which is the most precise at the price of being more computationally expensive than other solution methods that exist for the solution of PBM. To the best of our knowledge, this is the first time a sequential two-way coupling approach has been used

to study the effect of changing characteristic times on the latex quality.

The remainder of the thesis is covered as follows:

Chapter 2 is a bibliographic review covering emulsion polymerization fundamentals and modeling, computational fluid dynamics and fluid mixing. Besides, it includes the modellings used in this work.

Chapter 3 discusses in detail the CFD simulation of the single phase flow in the stirred tank EP reactor. The aim of this chapter is to use the appropriate flow models through the comparison of CFD simulation results with the experimentally measured ones, in order to achieve a validated single phase flow which will be coupled to the PBM to complete the simultaneous CFD/PBM computational framework in Chapters 4 and 5.

Chapter 4 deals with the modeling of perikinetic coagulation phenomenon and investigates its effect on the evolution of PSD through the use of the developed simultaneous framework. CFD simulation is used to obtain the local distribution of the species participating in the coagulation phenomenon. This data is used in each time step in the PBM in order to calculate the particle diameter which will be used in the next time step in the flow equations. In this way a two-way coupling is assured. The mentioned framework is used to investigate the effect of vessel scale on the perikinetic coagulation of the latex.

Chapter 5 is built on the CFD/PBM framework with the objective of studying the simultaneous perikinetic and orthokinetic coagulation of the particles, in order to study the effect of the operating conditions and configuration of the reactor on the evolution of PSD. CFD simulations is used to produce the hydrodynamic data playing role in the shear-induced coagulation; this information is then used by the PBM in order to calculate the particle diameter. The effect of reactor scale-up through maintaining geometric similarity and volumetric power input is investigated on the simultaneous coagulation of polymer particles.

Chapter 6 is a summary of the main conclusions and a list of recommended guidelines for future investigation.

Chapter 2

Literature Review and Review of Modeling

2-1-Polymer Lattices

2-1-1-Early History

Synthetic polymers can be denoted as the materials of the twentieth century. Since World War II the production volume of polymers has increased by a factor of 50 to a current value of more than 311 million tons annually [1]. The enormous growth of synthetic polymers is due to the fact that they are lightweight materials, act as insulators for electricity and heat, cover a wide range of properties from soft packaging materials to fibers stronger than steel on a kilo per kilo basis, and allow for relatively easy processing. Moreover, parts with complex shapes can be made at low cost and at high speed by shaping polymers or monomers in the liquid state [2]. The major application areas of polymers can be divided into five categories: automotive industry, building and construction, electrical and electronic industry, packaging, and agriculture.

Polymers are composed of very large molecules, each of which includes a large number of repeating structural units. The oldest group of polymers consists of natural polymers, among which natural rubber occurs in the form of a latex that is defined as a viscid milky juice. By far the most important natural latex is the one obtained from the rubber tree Hevea brasiliensis which is usually called "natural latex". The need for large amounts of synthetic rubber arose as a result of World War II which led to a world-wide natural rubber shortage and a motivation to develop the emulsion polymerization (EP) process technique. Since the natural rubber is in the form of latex, it is obvious that the hope was mimicking of it, at least as a starting point.

Emulsion polymerization and polymer latex technology started due what war necessitated, remained as active fields of research afterwards. In addition to improved product properties, emulsion polymerization was found to offer a number of additional advantages over bulk and solution polymerization, mainly higher rates of productivity and controllability.

In academia, the Harkins [3] and Smith-Ewart [4] theories were the most prominent and important ones served the understanding of the mechanistic and, subsequently, kinetic mechanisms. Today, a wide range of synthetic latexes, with rubbery and non-rubbery polymers, are produced using emulsion polymerization. The range of applications includes paints, coatings (including paper finishing), adhesives and carpet backing with annual production exceeding 20 million tons [2].

2-1-2-Emulsion Polymerization Fundamentals

A typical formulation for emulsion polymerization contains monomer, water, surfactant, initiator and additives such as buffers and chain transfer agents. Usually emulsion polymerization is carried out in stirred-tank reactors, which commonly operate in a semi-continuous mode, although both batch and continuous operations are also used. If we consider a batch emulsion polymerization process the monomer is dispersed in water using emulsifiers or surfactants, which normally is of anionic or nonionic type, or a combination of both. The surfactant partitions between the surface of the monomer droplets and the aqueous phase and forms emulsified monomer droplets. Since in most recipes the amount of surfactant is more than the amount needed to saturate the aqueous phase (which is called critical micelle concentration or CMC) and completely cover the monomer droplets, it forms micelles. Since the interior of the micelles is hydrophobic, the micelles will be swollen with monomer. Besides, a small fraction of the total monomer in the system will remain in the aqueous phase. Polymerization is commonly initiated using a water-soluble initiator, although oil-soluble initiators may be used.

As mentioned earlier, the original emulsion polymerization reaction theory was proposed by Harkins, and it divided the EP process into three intervals: During Interval I particle nucleation occurs. Interval II is the particle growth stage, and Interval III is the stage where the most of polymerization occurs and as a result the remaining monomer in the droplets consumes.

Since the first radicals, which are formed by addition of a water-soluble initiator are too hydrophilic to enter the organic phases of the systems, they react with the monomer dissolved in the aqueous phase and form oligoradicals. After adding some monomer units, they become hydrophobic enough to be able to enter the micelles (entry into the monomer droplets is not likely because their total surface area is about three orders of magnitude smaller than that of the micelles). Since the concentration of monomer in the micelle is high, the oligoradicals that have entered the micelle grow quickly and form polymer chains. The process of formation of polymer particles by entry of radicals into micelles is called micellar nucleation [2]. That part of oligoradicals that do not enter into micelles undergo homogeneous nucleation. After adding some monomer units they become too hydrophobic and precipitate. Later, by adsorption of the emulsifier present in the system onto their surface they become stabilized. Then, monomer diffuses into the new polymer particles. In general, homogeneous nucleation is predominant for monomers of relatively high water solubility and heterogeneous nucleation is predominant for water-insoluble monomers [2].

During stage I, monomer droplets, monomer swollen micelles, and monomer swollen polymer particles coexist in the reactor. The monomer is consumed by polymerization inside the polymer particles that compete efficiently for radicals. If an external source of monomer is available, monomer consumed in the particles is replaced by monomer that diffuses from the monomer droplets through the aqueous phase. So the size of the particles increases and that of the monomer droplets decreases. On the other hand, the number of micelles decreases because they become polymer particles or disappear to contribute to the stabilization of polymer particles. After some time, all of the micelles disappear. This is considered to be the end of the nucleation. Unless coagulation occurs, the number of particles remains constant during the rest of the process [2].

At the end of interval I the system is composed of monomer droplets and polymer particles. The monomer is still being consumed by polymerization in the polymer particles and then being replaced by the monomer diffusing from the monomer droplets through the aqueous phase. Monomer partitions between different phases of the system according to the thermodynamic equilibrium and the concentration of the monomer in the polymer particles reaches a maximum value in the presence of monomer droplets.

Because of the polymerization the polymer particles grow in size, and after some time the monomer droplets disappear. This is the end of Interval II. In Interval III, the monomer concentration in the polymer particles decreases continuously. The final product is a waterborne, concentrated dispersion of tiny polymer particles called latex.

However, as mentioned in chapter 1, the overall goal of the joint research project of which this thesis is a part is to study the emulsion polymerization of vinylidene fluoride (VDF) by means of CFD/PBM coupling approach. It should be noted that since the VDF monomer is present in supercritical state in the overhead gas phase from which it is transferred to the polymerization loci via the continuous aqueous phase, the emulsion polymerization of VDF does not follow the typical stages of the classical emulsion polymerization process that was described previously. In other words, even if the pressure were to drop in a batch process of this type, monomer would continue to be transferred to the polymer particles. If we were to use the Harkin's model process, the VDF process will never leave Interval II. Furthermore, for many reasons of product quality and environmental concerns, the amount of surfactant used in these processes is often below the CMC, so the nucleation process will take place without micelles. As will be seen later on in this document, we will concentrate on modelling the stability of particles, rather than their rates of production (mostly for reasons linked to the highly time consuming nature of the complex situations we will undertake). For this reason, we will not focus on the differences between the particle nucleation that takes place in the system of interest and that described by the idealized Harkin's model -the latter is simply presented to give the uninitiated reader a sense of how reactions of this type usually take place.

2-1-3-Stability of Polymer Latexes

Without sufficient stabilization (e.g. by surfactant, charged species on particle surface etc.), thermodynamics will drive the polymer latexes to minimize the surface area between the oil and water phases. Thus, unstable systems will see a certain amount of aggregation and coalescence of the particles that leads to a reduction in Gibbs free energy of the system which accompanies the loss of interface between two phases (the particles and the aqueous phases). This in turn results in a positive interfacial free energy between the two phases. The intermolecular forces of attraction cause molecules in the condensed state to cohere together [5].

There are three main mechanisms affecting the relative motion of the two particles and that eventually govern coagulation in emulsion polymerization. The particles in a polymer latex are continually undergoing Brownian motion which brings them into close proximity with each other (diffusive forces). Velocity gradients in the flow field at the microscopic scale may enhance the motion of particles relative to one another (convective forces). Furthermore, interparticular forces act as the third mechanism affecting the relative motion of the two particles. Therefore, there is both a thermodynamic tendency for the particles to aggregate and coalesce, and also a mechanistic pathway by which these processes can occur [5]. The contributions of these mechanisms to the overall coagulation rate depend both on the polymerization recipe and the mixing conditions, namely, on particle size, shear rate, and colloidal stability.

Most polymer latexes are kinetically stable, and this stability is a consequence of existing barriers between the particles arising from the balance between the various attractive and repulsive forces which are operative between two particles when they approach each other closely. The higher the potential energy barrier, the more stable the latex.

The attractive potential is due to van der Waals forces and is essentially related to the composition and shape of the polymer particles and the composition of the medium in which they are dispersed. These forces arise from electric dipoles in the atoms which may be permanent, induced by other molecules or simply a result of fluctuations in the atoms' electron clouds. The forces of repulsion which act between particles can be grouped into two general types: "electrostatic forces", which depend on the presence of an electric charge at the surface of the particles; and "steric forces", which arise from the presence of hydrophilic molecules bound to the surface of the particles. Further information can be found in Vale & McKenna [6] and references cited therein. In this work we are concerned only with electrostatic repulsion due to the formulation of the product of interest.

Although creating a repulsion force stabilizes colloidal systems, it is possible that in systems subjected to shear the energy barrier between particles is not sufficiently high to prevent coagulation. In fact collisions of particles are strongly affected by flow pattern, which might vary significantly in the course of the process, especially for turbulent flows.

2-2-Emulsion Polymerization Modeling 2-2-1-Introduction

Models of Emulsion Polymerization can be classified into two levels according to the way they treat the particle size [6].

1-In level-one models the assumption is that all particles have the same average volume (mono-dispersed particles). Under circumstances where the PSD is relatively monodisperse and monomodal, this is a useful and practical simplification which is used in some cases and often provide useful information.

2-Level-two models are closer to reality since they account for the particle size distribution by means of a transport equation: the population balance equation. Of course there will still be different levels of complexity for this type of Level-two model, for example the way of handling the differential equation to consider the particle size change.

In general the number and the type of equations also represents a big challenge in modelling the PSD of EP process. The interested reader can consult a review on the uses and limitations of PBEs in modeling emulsion polymerization [6].

2-2-2-Overview of Emulsion Polymerization Kinetics

Emulsion polymerization follows the same basic mechanisms as in free-radical polymerization with the main difference due to the confinement of radicals within particles, the phenomena so called compartmentalization, which is particular to emulsion polymerization. For this reason, the particles containing one radical should be treated as being quite distinct from those with two or more radicals. Accordingly different rate equations should be used to describe the kinetics of particles containing different number of radicals. Moreover, emulsion polymerizations exhibit events that have no counterpart in bulk or solution polymerizations, such as phase transfer processes, where radicals move between the particle and the water phases in two ways, radical entry and radical exit.

The basic rate equation for homogeneous batch free-radical polymerization in bulk and solution polymerization is normally defined based on the rate of consumption of monomer:

= [ ] = [ ][ ] (2.1)
where is the propagation rate coefficient, [ ] is the concentration of monomer and [ ] is the total radical concentration.

Since in emulsion polymerization the locus of polymerization is within the latex particles, [ ] is replaced by [ ] and the total radical concentration is . So equation (2.1) is modified as:

= [ ] (2.2)
where is the average number of radicals per particle, is the total particle concentration per unit volume of aqueous phase, and is Avogadro's number.

Generally, three existing phenomena in the production of polymer latexes can effect : nucleation, by which new polymer particles are formed; growth of these newly formed particles due to polymerization; and coagulation of the particles. So a means of tracking changes to this parameter during the course of polymerization is necessary when modeling emulsion polymerization processes.

2-2-3-Modeling Particle Coagulation

The production of many polymers often involves stages in which the monomer, the polymer, or the intermediate product is dispersed into small particles. This fact that the system is in the form of a dispersion brings along some difficulties which are mainly related to the colloidal stability and thus the process of coagulation and fouling [7].

Coagulation is a major cost to the latex manufacturing industry. This phenomenon can have a negative effect on the product quality and causes large costs due to product loss and reactor down time for cleaning, because the coagulum reduces the efficiency of heat exchange. On the other hand, coagulum formation is one of the major problems when scaling up a reaction from laboratory to industrial size, as the agitation regime in the reactor changes.

One of the main issues in the field of emulsion polymerization modeling has been the absence of a semi-phenomenological coagulation model that can predict the rate of coagulation in both the presence and absence of shear rates, the main reason of which being high mathematical complexity that it can add to the systems of equations and also the lack of experimental data.

For modeling purposes, the coagulation between particles is assumed to occur as a result of binary collision events, although the assumption that all particles will coagulate through a binary collision mechanism is valid only when working with dilute systems.

However, due to the complexity of formulating and solving the equations for multi-body collisions, binary collision models have been used at all latex concentrations [8].

Due to the complexity of determining the contribution of different mechanisms affecting the relative motion of the two particles, two cases are typically considered in the literature, where particle aggregation is dominated by Brownian diffusion (perikinetic coagulation) or by transport due to fluid motion (orthokinetic coagulation) [9]. Even though in reality both cases coexist and affect the coagulation rate simultaneously.

Therefore, it will be useful to have an estimate of the coagulation rate in the ranges between the two extremes of pure perikinetic coagulation and Pure orthokinetic coagulation.

Von Smulochowski's [10] approach to solve the convection-diffusion equation which describes the distribution of particles around a reference particle, serves as a basis for the quantitative description of aggregation kinetics. Nearly a century ago he solved this equation using these boundary conditions: at the reference particle surface the concentration is set to zero, since it is assumed that the particles that come in contact with the reference particle will undergo irreversible aggregation. At an infinite distance from the reference particle, the particle concentration is set at its bulk value. He ignored the interparticle forces and hydrodynamic interactions which arise between particle pairs as they approach one another (essentially a form of viscous friction) and considered that the motion of particles relative to one another arises only from their Brownian motion.

Under the assumption that all collisions between particles are of a binary nature, the change in particle concentration ( ) with time is:

= (2.3)
where is the rapid coagulation rate coefficient.

Fuchs [11] modified the work of Smulochowski by taking the interparticle forces into account as well: the diffusion-convection equation is modified to include a convective term that captures the net effect of the attractive and repulsive forces that may act on a particle as they diffuse towards one another; but the hydrodynamic interactions are still ignored. Considering all of the three effects will result in the addition of a second convective term such that the net convective rate is a sum of the interaction and shearinduced convective terms: [11] . ( . .

( + ) ) = 0 (2.4)
This is the basis of most modern coagulation models. Here is the diffusion tensor (modeled as Brownian diffusion and often modified with terms to account for hydrodynamic interactions), is the convective velocity and is the velocity induced by interparticle interactions.

2-2-3-1-Perikinetic Coagulation

For the simplest case, where aggregation of non-interacting particles occurs by Brownian diffusion only (usually referred to as the fast aggregation limit), the estimates of the coagulation rate in EP are frequently obtained on the basis of the DLVO theory of colloid stability. DLVO theory framework, despite its limitations, is currently the best model available to model the interaction velocity vector . Under the DLVO framework, may be computed as follows:

= .

(2.5)

where is the gradient in the interaction potential and is the diffusion tensor which is expressed as a function of the diffusion coefficient and may be modified to account for hydrodynamic interactions:

= ( , ) 0 0 0 ( , ) 0 0 0 ( , ) (2.6)
where ( , ) and ( , ) are the radial and angular (polar) hydrodynamic interaction functions. When two particles approach one another, fluid must be squeezed out from the gap between the particles; this action creates a viscous drag on the particles.

The mutual diffusion coefficient, is expressed according to the Stokes-Einstein equation as:

= 6 1 + 1 (2.7)
In this theory particle stability is dependent on the total particle potential energy of interaction, . As stated previously, the total potential energy of interaction is the net sum of all attractive and repulsive contributions, including van der Waals' attraction and electrostatic repulsion. The reader may refer to Ottewill [12] for details on calculating each of these contributions and the overall perikinetic coagulation rate coefficient.

Treated within the DLVO framework, the total interaction potential energy between particles is assumed to be the sum of the attractive and the repulsive energies between particles: [12] = + (2.8)

is the overlap potential between the diffuse double layers of the particles. According to Hamaker, the energy between two interacting particles, with radii and , can be described as: [13] = 6 2

+ + 2 + + (2.9)
where is the center-to-center separation.

The most widely used expression for computing the electrostatic repulsive energy potential is the Hogg-Healy-Fürstenau expression: [14]

= ( + ) 4( + ) 2 + 1 + exp ( ) 1 exp ( ) + [1 exp ( 2 )] (2.10)
which is a function of the zeta potential ( ) of each particle, the Debye-Huckel parameter ( ), and the distance between the surfaces of the two particles ( ). In this equation is the permittivity constant of water which is a function of the permittivity constant of vacuum ( ) and water ( ):

= 4 (2.11)
The distance between the surfaces of the two particles is calculated as a function of the particle radius of each particle and the center-to-center separation by the following equation:

= + (2.
12)

The Debye-Huckel parameter is related to the ionic force ( ) by the equation:

= 8 . (2.13)
where is the electron charge, is the Boltzman constant, and is the temperature.

The ionic force included in this equation is given by:

= (2.14)
where and are the concentration and the valenceof the ionic species, respectively.

The potentials are determined as a function of thesurface potential and the Stern layer thickness :

= 2 exp( ) + 1 exp( ) 1 (2.15) = + ln exp( ) + 1 exp( ) 1 (2.16) = 2 (2.17)
where is the counter-ion valence.

The appropriate description of the surface potential depends on the product ; for values of this product lower than 1, spherical surfaces may be approximated to plate surfaces and is computed based on following equation:

[15] = 4 (1 + ) (2.18)
For values of the product higher than 1, spherical geometry has to be assumed, and then surface potentials are calculated by: [15]

= 2 2 (2.19)
Finally, the coagulation rate between two particles with the radii and ( ) is dependent on the Fuch's stability ratio ( ): [15] = 2 ( + ) (2.20)

= = 2 3 + (2.21)
The total surface charge on the particles is obtained by adding the charges due to added surfactant and generated charges:

[16] = + (2.22)
The generated charges due to initiator decomposition is very low compared to the charges due to adsorption of surfactant on the particles' surface, which is given by: [16]

= (2.23)
where is the area actually occupied by a surfactant molecule.

The concentration of surfactant in the aqueous phase [ ] is dependent on the concentration of added surfactant per unit volume of the continuous phase [ ] and the total amount of surfactant adsorbed onto the particles' surfaces: [16]

[ ] = [ ] (2.24)
where is the total area of particles per unit volume of the reactor.

Since the amount of surfactant is very low in the particular case of interest here, it is assumed that the whole amount of added surfactant is adsorbed onto the particles' surface, therefore:

= [ ] (2.25)
Combining equations (2.22) and (2.23), the surface charge density on the surface of particles can be obtained:

= = [ ] (2.26)
While the DLVO framework is widely used in modeling emulsion polymerization coagulation, the theory has quantitative limitations with respect to the prediction of the coagulation rate of colloids. Behren et al. [17] found that while DLVO theory was quantitatively applicable for weakly charged particles (i.e. surface charge less than 3 mC/m2) and low ionic strengths (< 10 mM), discrepancies emerged between theory and experiments at high ionic strengths. An analysis of the interaction energy profiles suggested that such deviations were related to the position of the energy barrier. At high ionic strengths the energy barrier appears at a separation distance less than 1 nm and non-DLVO forces that are not included in the model become influential. These non-DLVO forces arise from surface heterogeneities, the discrete nature of charges and the finite size of ions and water molecules. Since higher ionic strengths and higher surface charges are typically encountered in emulsion polymerization, the DLVO interaction model should be used with caution. Another important aspect has to do with the fact that the DLVO theory supposes dilute dispersions, for which the influence of the surrounding particles on the pair of interacting particles can be neglected. In concentrated systems (e.g. high solid content latexes), the surrounding particles cause an effective reduction in the total potential energy of the interacting particles leading to a lower stability ratio [18][19][20]. The concentration effect may lead to a higher rate of coagulation that would be predicted by DLVO theory.

The surface charge density of the particles is given by the sum of the charges due to ionic end groups, adsorbed surfactant, and in situ generated surfactant, among all the adsorbed surfactant has typically the most important contribution. In the development of DLVO-based models, it is assumed that every particle in the system has the same surfactant coverage which is invariably calculated by assuming thermodynamic equilibrium between the aqueous phase and the surface of the particles. However, in reality the precursor particles formed from homogeneous nucleation may undergo coagulation rapidly, such that the surfactant does not have sufficient time to diffuse to the newly-created surface and equilibrate. Obviously, this is not a limitation of the DLVO theory itself, but a consequence of applying an equilibrium hypothesis to a dynamic process such as emulsion polymerization. Despite all this, DLVO theory remains the only basis for non-empirical modeling of particle coagulation. It should be noted that it is possible to fit certain model parameters to obtain coagulation rate coefficients that correctly predict experimental results. In Chapter 4, results of a detailed study performed on perikinetic coagulation phenomenon of a PVDF EP will be presented.

2-2-3-2-Orthokinetic Coagulation

While good mixing of the ingredients is important to ensure reactor homogeneity, it is also the cause of shear-induced coagulation. This poses a significant performance tradeoff for the manufacture of latex by emulsion polymerization. Many researchers neglect orthokinetic coagulation when modeling emulsion polymerization, as perikinetic coagulation is thought to be the dominant mechanism at small particle sizes (~ 50 nm) [21]. This simplification is often justified based on experimental results, [22,23] however theoretical studies of particle coagulation suggest that hydrodynamics can play an important role when ionic strength is low and the double layer is thick [9]. The development of models for this effect would improve the production process and scaleup [24].

2-2-3-2-1-Laminar Flow

Smoluchowski [10] with the assumption of monodispersed particles and binary coagulation developed the following expression for determining the orthokinetic coagulation rate coefficient for laminar shear flow:

= 4 3 ( + ) (2.27)
The velocity gradient used in this equation is that for a flow exhibiting a simplified two-dimensional form of pure-shear strain where only one component of the relative velocity is considered [25]. In order to generalize the work of Smoluchowski, Camp and Stein [26] replaced the velocity gradient of Eq. (2.27) with a local parameter ( ) referred to as the absolute velocity gradient (1/s) and defined as:

= = + + + + + (2.28)
where is the viscous energy dissipation rate per unit volume and is the dynamic viscosity (Pa.s), and so proposed the following equation for the orthokinetic coagulation rate:

= 4 3 ( + ) (2.29)

2-2-3-2-2-Turbulent Flow

For turbulent flow a global parameter rather than a local one as in laminar case ( ), referred to as the root mean-square velocity gradient, was defined by Camp and Stein and expressed as:

= (2.30)
In this equation is the mean value of the work input to the tank per unit time and unit volume. The validity of the Camp and Stein approach has been questioned [START_REF] Pedocchi | Camp and Stein's velocity gradient formalization[END_REF][START_REF] Clark | Critique of Camp and Stein's RMS velocity gradient[END_REF][START_REF] Saatçi | Discussion of "Critique of Camp and Stein's RMS Velocity Gradient[END_REF][START_REF] Kramer | Influence of strain-rate on coagulation kinetics[END_REF][START_REF] Graber | Discussion of "Influence of strain-rate on coagulation kinetics[END_REF].

Pedocchi and Piedra-Cueva [START_REF] Pedocchi | Camp and Stein's velocity gradient formalization[END_REF] concluded that for laminar flows, Camp and Stein approach is valid, provided that the following expression for the viscous dissipation function is utilized:

= 2 + 2 + 2 + + + + + + (2.31)
A similar expression, Eq. (2.32), was derived by Saffman and Turner [START_REF] Saffman | On the collision of drops in turbulent clouds[END_REF] and also presented by Spielman [START_REF] Spielman | Hydrodynamic aspects of flocculation[END_REF] based on homogeneous isotropic turbulence with particles smaller than the Kolmogorov microscale ( ). This equation differs from the Eq. (2.29) not only in the value of the constant used, 1.333 versus 1.294, but also in the definition of the rate of energy dissipation ( ), where rather than is used:

= 8 15 ( + ) (2.32)
Some workers [START_REF] Lowry | Mechanical Coagulation in Emulsion Polymerizations[END_REF][START_REF] Lowry | Kinetics of Agitation-Induced Coagulation of High-Solid Latexes[END_REF][START_REF] Chern | Shear-induced coagulation kinetics of semibatch seeded emulsion polymerization[END_REF][START_REF] Mayer | Dynamic modeling of limited particle coagulation in emulsion polymerization[END_REF] have addressed the issue of orthokinetic (shear-induced) coagulation. However, the shear-induced coagulation as used in these models makes use of an averaged shear strain rate over the entire reactor as a function of the impeller type and/or the power usage, so it is in fact a simplistic method, because it does not take into account the fact that the shear-induced coagulation rate will vary spatially in the reactor.

Nevertheless this method can provide an estimate of the effect of shear on the PSD.

Elgebrandt et al. [24] investigated the difference in the predictions between the two methods, one which uses the average shear rate and the other that uses the local shear rates in order to calculate the rate of simultaneous (perikinetic/orthokinetic) coagulation and found out that using local shear rates produces more precise and useful results.

2-2-3-3-Simultaneous Perikinetic/Orthokinetic coagulation

As mentioned previously, in reality perikinetic and orthokinetic coagulations coexist and effect the PSD evolution simultaneously. More precisely, since Brownian motion is always present, regardless of whether the shear flow exists or not, no pure orthokinetic coagulation can actually exist. When the shear flow does exist, coagulation should be simultaneously orthokinetic and perikinetic [38].

The rate of aggregation between two particles can be obtained from the solution of Equation (2.4) by computing the total flux of particles towards the reference particle (refer to [9] or [39] for more details regarding the numerical techniques employed):

= (

) .

(2.33)

where the integral runs over the entire collision surface and is the unit vector normal to the surface and pointing inwards. Note that the integral can be evaluated on any closed surface surrounding the reference particles, since application of the divergence theorem to Equation (2.4) implies that the integral of the total flux is conserved over any closed surface surrounding the reference particle [39]. The boundary conditions of the convection-diffusion equation (Equation 2.4) are as follows: = 0 at = + and = at = .

Batchelor [START_REF] Batchelor | Mass Transfer from a Particle Suspended in Fluid with a Steady Linear Ambient Velocity Distribution[END_REF] indicated that a turbulent flow can be effectively represented by assuming an axisymmetrical extensional flow arising from isotropic turbulence as the latex particles are smaller than the Kolmogorov scale (the length scale at which the directionality of the turbulent eddies is lost, typically in the order of a few microns).

Therefore the components of the relative particle velocity in spherical coordinates are given by Batchelor and Green: [START_REF] Batchelor | The hydrodynamic interaction of two small freelymoving spheres in a linear flow field[END_REF]

, = 2 (1 ( , ))(3 cos( ) 1) (2.34) , = 2 (1 ( , )) cos( ) ( ) (2.35)
where is the radial distance from the center of the reference particle and is the colatitude (the difference between the latitude and 90°) and is the local shear rate. = is the dimensionless particle size ratio between the interacting particle pairs. ( , )

and ( , ) are the two hydrodynamic functions defined by Batchelor [START_REF] Batchelor | The hydrodynamic interaction of two small freelymoving spheres in a linear flow field[END_REF], which account for hydrodynamic interactions between a pair of particles in an axisymmetric linear flow field.

In order to facilitate its numerical solution, Equation (2.4) has been rewritten in dimensionless form and in spherical coordinates: [39] 1 ( ( , ) . ( , )(3 cos( ) 1) ( , )

(2.36)

+ 1 sin ( ) sin( ) ( , ) ( . ( , ) cos( ) sin ( )) = 0
where = is the dimensionless distance, = is the normalized particle concentration, and = is the normalized interaction potential. In the case of extensional flow ( , ) and ( , ) are defined as:

( , ) = 2 (1 ( , )) (2.37) ( , ) = 3 2 (1 ( , )) (2.38)
The Peclet number, , which provides a measure of the relative importance of shear versus diffusion, is defined in the case of extensional flow as:

= 3 + 2 (2.39)
Melis et al. [9] tried to investigate theoretically the changes that arise in the magnitude of as the ratio of fluid convection to particle diffusion changes (i.e. when varies).

He tried to develop an analytical expression to fit the simulation data by using the numerical solution to Equation (2.36). He found out that under unstable conditions (low electrostatic forces), an assumption of the perikinetic and orthokinetic contributions being additive in nature works properly and the rigorous model (Equation (2.36)) results could be fit with the following semi-empirical expression:

= + = 8 + . (2.40)
where is an adjustable parameter.

Nevertheless, Melis found that Equation (2.40) provided a poor fit for more stable systems. This indicates that for slowly aggregating systems, the diffusive and the convective aggregation mechanisms interact with each other, and, therefore, the additivity assumption is not valid [9].

Lattuada and Morbidelli [39] derived an analytical expression that reproduces the results of the numerical simulation by neglecting the angular dependence in the velocity profile functions in Equation (2.36) and so using a simplified velocity profile. Their model requires the use of a boundary layer approximation and finally the following expression is derived for . Expressed in terms of a stability ratio, :

= 1 2 exp( (2 + ) . ( ) + 2 + . ( 2 ) ( , ) (2.41)
where the boundary layer thickness and is an adjustable parameter.

= 2 ( + ) (2.42)
where is an adjustable pre-factor that is used to fit the simplified expression's predictions to the numerical data generated from the rigorous solution.

( ) is the Heaviside function and is the diffusion-limited rate of aggregation:

= 2 3 + (2.43)
The presence of two adjustable parameters in Equation (2.41) may limit its applicability.

If the complex velocity function is replaced by a symmetrical velocity gradient

( , =
), the following expression can be ob tained (refer to section 4.6.1.3 of Blackley [5] for the derivation), which is the expression used in this study in order to compute the simultaneous orthokinetic/perikinetic coagulation phenomena:

= = 1 2 exp 8 3 ( , ) (2.44)
While inclusion of hydrodynamic interaction ( ( , )) brings about better accord with experiment (e.g., Peula et al. [START_REF] Peula | Comparative study on the colloidal stability mechanisms of sulfonate latexes[END_REF]), this effect is not included in this study, because our objective is to establish basic methodology rather than attempt precise accord with experiment. In Chapter 5, results of a detailed study performed on simultaneous perikinetic/orthokinetic coagulation phenomenon of a PVDF EP will be presented.

2-2-4-Modeling the Particle Size Distribution (PSD)

The particle size distribution is one of the most important characteristics of a polymer latex determining its rheological properties, maximum solid content, adhesion, drying time, etc. The modeling of the rheological properties of the latex depends on knowledge of the PSD. High solid content latexes are an excellent example of a product requiring an accurate control of the PSD. From the perspective of economic productivity and improved latex properties, there is a strong interest in increasing the solids content of latex formulations [START_REF] Guyot | High solid content latexes[END_REF], mainly for two reasons: (i) the production and transport of concentrated latexes is more economical, an especially-important consideration for manufacturers of commodity-grade latexes; (ii) high solid content latexes are uniquelysuited to a wide range of applications (for example, fast-drying paints and coatings) [START_REF] Pohn | Scale-up of Latex Reactors and Coagulators: A Combined CFD-PBE Approach[END_REF].

One advantage inherent to synthetic polymer latexes is the ability to control the PSD and so the properties of the latex product by manipulating the manufacturing process.

Recently published results suggest that it is possible to prepare low viscosity, stable latexes with solid loading above 70 vol.% if the PSD is well-engineered [START_REF] Pohn | Scale-up of Latex Reactors and Coagulators: A Combined CFD-PBE Approach[END_REF]. A means of tracking changes to the particle size distribution during the course of polymerization is necessary when modeling emulsion polymerization processes that are sensitive to changes in the PSD. That will be accomplished by level-two models through the addition of a set of population balance equations (PBEs) to the kinetic model.

We may consider an open system where the particles are distributed according to their size, , and position, , and let the domains of and be represented by (internal coordinates) and (external coordinates). In addition, we may assume that there exists an average number density function, ( , , ) such that the average number of particles in the infinitesimal volume would be ( , , )

.The population balance for this density function of particles with the properties ( , ) can be shown to give: [START_REF] Ramkrishna | Population Balances Theory and Applications to Particulate Systems in Engineering[END_REF] ( , , ) =

where ( , , ) is the rate of particle growth, ( , , ) is the rate of change of the external coordinates (i.e. particle motion in the space), and ( , , ) is the net rate of generation of particles in the infinitesimal volume . The last term on the RHS of Equation (2.45) may include a variety of phenomena, namely particle formation and depletion due to coagulation. It may also account for nucleation, but nucleation is usually treated through the boundary condition [6].

The equation presented above hold true irrespectively of the variable chosen as a measure of the particle size (radius, volume, mass, etc.). However, the derivation of the PBE is somewhat simpler when done in terms of the volume of the particles [6]. This is because both the expression for the rate of particle growth and the coagulation kernels are easier to write in terms of these variables. The relationship between the two density functions is given by: [6] ( , , ) = ( , , ) = 4 ( , , )

For a review on the uses and limitations of PBEs in modeling emulsion polymerization, the reader is referred to the comprehensive review written by Vale and McKenna [6]. integrating. The resulting algebraic system of equations is solved to yield the unknown coefficients of the piecewise approximating functions [START_REF] Bathe | Finite Element Procedures[END_REF]. In the finite volume method, the conservation equations are integrated over all control volumes of the domain, and the resulting integral expressions are discretized to form a system of algebraic equations, which are solved using an iterative method. Over the years, advances in research have addressed both methods' weaknesses. For process engineers, the choice of software is likely to be strongly-influenced by software compatibility and the availability of both technical support and application-specific support in the research community.

2-3-Simulation of Reactor Performance Using Computational Fluid Dynamics

2-3-2-Grid Generation

Grid generation strategies can be classified as Cartesian, structured, unstructured, hybrid and gridless [START_REF]Grid generation: Past, present, and future[END_REF]. In Cartesian gridding a network of grid lines containing uniform spacing is placed in a rectangle (2-D) and a rectangular box (3-D) application. The boundary conditions implementation is established by cutting interior geometrical entities with grid lines [START_REF]Grid generation: Past, present, and future[END_REF]. While Cartesian gridding is the simplest way to discretize a given field, the approach is generally not suited for solving complex computational fluid flow problems with complex geometries.

The structured grid is represented by a network of curvilinear coordinate lines such that a one-to-one mapping can be established between physical and computational space [START_REF]Grid generation: Past, present, and future[END_REF]. Since the curvilinear grid points conform to the solid surfaces/boundaries, this method provides the most economical (in terms of computational expenses) and accurate way for specifying boundary conditions. Thus the structured gridding has the advantage of being computationally efficient, requiring fewer cells to compute the solution to a desired level of accuracy. The main drawback to structured grid generation is that the method is not automated. That is for complicated geometrical configurations, the physical region should be divided into sub-regions and within each a structured grid should be generated. The transfer of solution information at the block interface is very critical for successful simulation [START_REF]Grid generation: Past, present, and future[END_REF].

Unstructured grids are composed of triangles (2-D) and tetrahedrons (3-D). The grid information is presented by a set of coordinates (nodes) and the connectivity between the nodes. The connectivity table specifies the connections and appropriate neighborhood information between nodes and cells [START_REF]Grid generation: Past, present, and future[END_REF]. The unstructured grids offer a greater geometric flexibility, but at the same time they require a higher number of elements relative to structured grids. On the other hand the automatic generation software may generate poor quality elements (skewed) in regions with high aspect ratios.

These skewed cells may potentially introduce errors in the numerical solution.

The usual practice in Hybrid gridding is to generate structured grids near solid components where high-aspect ratio cells are required, and to fill in the remaining region with an unstructured grid. Doing this, the method profits the advantages of both structured/unstructured gridding methods letting one the generation of meshes of high quality with a large degree of automation, but the algorithms are often changeable and may require the user to be experienced with the pre-processing software before they function properly.

There is another approach called gridless approach in which the numerical treatment of governing equations is performed without requiring any type of explicit connectivity links between points. The only structure required is a distribution of cloudless points within the computational domain and the discretized numerical scheme is developed based on the points registered in the neighborhood of each point. The development of this technology is in its infancy and has a long way to go before its utilization in practical industrial applications [START_REF]Grid generation: Past, present, and future[END_REF]. Anyhow the approach holds much promise in modeling the unsteady motions of rigid bodies, where mesh quality would tend to degrade during large rigid movements of a body's boundaries [START_REF] Wang | A study of gridless method with dynamic clouds of points for solving unsteady CFD problems in aerodynamics[END_REF].

2-3-3-Simulation of Flow Field Using ANSYS Fluent

ANSYS CFD solvers are based on the finite volume method. They solve conservation equations for mass (or continuity equation) and momentum. The general form of the mass conservation equation which is valid for incompressible as well as compressible fluids reads as follows:

+ ( ) = 0 (2.47)
where is the density of the medium fluid, is the velocity vector and is a source term.

The conservation of momentum in an inertial reference frame reads as:

( ) + ( ) = + . ( ) + + (2.48)
where is the static pressure, is the stress tensor, and and are the gravitational body force and external body forces, respectively [START_REF]ANSYS Fluent 17.1, Theory Guide[END_REF].

For flows involving species mixing or reactions, a species conservation equation is solved. Additional transport equations are also solved when the flow is turbulent.

2-3-3-1-Turbulent Modeling

The The method of Large Eddy Simulation involves space filtering of the unsteady Navier-Stokes equations prior to the computations, which passes the larger eddies and rejects the smaller eddies. Large eddies are directly resolved, but eddies smaller than the mesh are modeled. This method is less expensive in terms of computational expenses than DNS, but still not practical for most practical applications.

In Direct Numerical Simulations the full unsteady Navier-Stokes equations are solved on spatial grids that are sufficiently fine that they can resolve the Kolmogorov length scales at which energy dissipation takes place and with time steps sufficiently small to resolve the period of the fastest fluctuations. These calculations are highly costly in terms of computing resources, so the method is not used for industrial flow computations.

The vast majority of turbulent flow computations has been based on the Reynoldsaveraged Navier-Stokes (RANS) equations. In Reynolds averaging, the solution variables in the instantaneous (exact) Navier-Stokes equations are decomposed into a steady mean value and a fluctuating component ( ) superimposed on it [START_REF] Versteeg | An introduction to computational fluid dynamics: The finite volume method[END_REF]:

( ) = + ( ) (2.49)
Substituting expressions of this form for the velocity into the instantaneous continuity and momentum equations, and taking a time (or ensemble) average yields the ensembleaveraged mass and momentum equations: where the shear stresses are important [START_REF] Ranade | Computational flow modeling for chemical reactor engineering[END_REF]. This model has already been utilized elsewhere [START_REF] Raikar | Experimental studies and population balance Equation models for breakage prediction of emulsion drop size distributions[END_REF][START_REF] Wang | Population Balance Model for Influence of Bubble Coalescence and Breakup Models[END_REF][START_REF] Zucca | Implementation of the population balance Equation in CFD codes for modelling soot formation in turbulent flames[END_REF] and has been shown to provide acceptable agreement with experimental results, particularly with regards to the overall flow field [START_REF] Aubin | Modeling turbulent flow in stirred tanks with CFD: the influence of the modeling approach, turbulence model and numerical scheme[END_REF][START_REF] Montante | CFDsimulations and experimental validation of homogenisation curves and mixing time in stirred Newtonian and pseudoplastic liquids[END_REF]. In the derivation of the Model, the assumption is that the flow is fully isotropic turbulent, and the effects of molecular viscosity are negligible. model is based on model transport equations for the turbulence kinetic energy ( ) and its dissipation rate ( ). The basic definitions of turbulent kinetic energy and its dissipation rate are as follows [START_REF] Launder | The Numerical Computation of Turbulent Flows[END_REF]:

+ ( ) = 0 (2.50) ( ) + = + + 2 3 + ( ) (2
= 1 2 + + (2.53) = (2.54)
where is the turbulent or fluctuating velocity and [START_REF] Bird | Transport phenomena[END_REF].

, the turbulent (or eddy) viscosity is computed by combining k and :

= (2.55)
where is a constant.

The transport equations for turbulent kinetic energy and its dissipation rate are as follows:

( ) + ( ) = + + + (2.56) ( ) + ( ) = + + ( + ) (2.57)
where represents the generation of turbulent kinetic energy due to buoyancy and is the generation of turbulent kinetic energy due to the mean velocity gradients and is computed as:

= (2.58)
where is the component of the flow velocity parallel to the gravitational vector and is the component of the flow velocity perpendicular to the gravitational vector. Considering the computational expenses and the preciseness and reliability of prementioned models, a decision was made to use the RNG Model to model the turbulence.

2-3-3-2-Transition Model

In this study Transition Model was used to address the transition of the boundary layer from laminar to turbulent regime in the case with 100 RPM agitation rate, which is the default model for modeling the transitional regime, and considering the type of flow in our study (confined flow) and the computational expenses of other options, it was preferred. Transition Model is a three-equation eddy-viscosity type that includes three transport equations: for turbulent kinetic energy ( ), laminar kinetic energy ( ), and the inverse turbulent time scale ( ):

= + + + + (2.59) = + (2.60)
The reader is referred to [START_REF]ANSYS Fluent 17.1, Theory Guide[END_REF] for the comprehensive details of the transport equations, the definition of parameters and model constants.

2-3-3-3-Multiphase Modeling

Currently there are two approaches for the numerical calculation of multiphase flows:

the Euler-Lagrange approach and the Euler-Euler approach. In Euler-Lagrange approach, the fluid phase is treated as a continuum by solving the Navier-Stokes equations, while the dispersed phase is solved by tracking a large number of particles, bubbles, or droplets through the calculated flow field. This approach requires that the dispersed second phase occupies a low volume fraction. In the Euler-Euler approach, the different phases are treated mathematically as interpenetrating continua. Since the volume of a phase cannot be occupied by the other phases, the concept of phasic volume fraction is introduced. These volume fractions are assumed to be continuous functions of space and time and their sum is equal to one. Conservation equations for each phase are derived to obtain a set of equations, which have similar structure for all phases. These equations are closed by providing constitutive relations that are obtained from empirical information. Since in this study we aim to use the PBM add-on module of Fluent and this module necessitates using the Euler-Euler approach, inwhat follows we discuss this option for modelling the multiphase flows.

There In the two last methods, which are of interest for our application, the phases are treated as interpenetrating continua. Mixture Model solves momentum equation for the mixture and prescribes relative velocities to describe the dispersed phases (but it can also be used without relative velocities for the dispersed phases to model homogeneous multiphase flow). The Eulerian model solves a set of momentum and continuity equations for each phase. Then the coupling is achieved through the pressure and interphase exchange coefficients.

"Mixture Model" was used in this study as the multiphase modeling approach. It is a good substitute for the full Eulerian multiphase model (it can perform a full multiphase simulation while solving a smaller number of variables than the Eulerian approach [START_REF]ANSYS Fluent 17.1, Theory Guide[END_REF])

and it is less computationally expensive compared to the Eulerian approach. In the Mixture Model, the momentum, continuity, and energy equations for the mixture, the volume fraction equations for the secondary phases, and algebraic expressions for the relative velocities are getting solved simultaneously.

The continuity equation for the mixture is given by: [50]

( ) + . ( ) = 0 (2.61)
where and are the mass averaged velocity and the density of the mixture respectively.

The momentum equation for a mixture containing phases reads as follows: [START_REF]ANSYS Fluent 17.1, Theory Guide[END_REF] ( ) + . ( )

+ . + + + + . , , (2.62)
where is the number of phases, is a body force, is the viscosity of the mixture and , is the drift velocity of the secondary phase .

Since in our case we deal only with two phases (one secondary and one continuous phase ), from now on, the equations are related to a two-phase case. The drift velocity (of secondary phase ) and the slip velocity (or relative velocity of phase relative to phase ) are related as follows: [START_REF]ANSYS Fluent 17.1, Theory Guide[END_REF] , = (1 )

(2.63)
where is the volume fraction of the secondary phase .

In Mixture Model, an algebraic formulation is used for the slip velocity which is based on the assumption of reaching to a local equilibrium between the phases over a short spatial length scale. In this work, the equation proposed by Manninen et al. [START_REF] Manninen | On the mixture model for multiphase flow[END_REF] was used for computing the slip velocity:

= (2.64)
where is the particle density, is the secondary phase particle's acceleration, and is the particle relaxation time given by:

= 18 (2.65)
here is the diameter of the particles of secondary phase and is the drag function.

The default drag function is taken from Schiller and Naumann [START_REF] Schiller | A Drag Coefficient Correlation[END_REF]:

= 1 + 0.15 . 1000 0.0183 > 1000 (2.66)
This drag function can be employed for the solid spherical particles or for the fluid particles that are sufficiently small and may be considered spherical such as our case in this study. Furthermore, according to Loth [START_REF] Loth | Quasi-steady shape and drag of deformable bubbles and drops[END_REF], empirical Schiller-Naumann expression is appropriate as long as the particle deformation is negligible. A few researchers have employed this drag expression and achieved good results [START_REF] Schütz | Modeling fluid behavior and droplet interactions during liquid-liquid separation in hydrocyclones[END_REF][START_REF] Yarin | Turbulence intensity in dilute two-phase fows-3: The particles turbulence interaction in dilute two-phase flow[END_REF][START_REF] Silva | Implementation and analysis of numerical solution of the population balance Equation in CFD packages[END_REF].

The acceleration is of the form:

= ( . )
(2.67)

2-3-3-4-Species Transport

The dynamic distribution of any species getting distributed in the reactor domain was calculated by solving a Reynolds-averaged time-dependent scalar ( ) transport equation, assuming a constant fluid density:

+ = (2.68)
where is the turbulent viscosity and is the turbulence Schmidt number.

2-3-4-Mixing

Mixing is usually defined as the reduction in inhomogeneity of concentration, phase or temperature to achieve a desired process result. Mixing is said to be occurred in three The study of mixing dates back many years before the first journal publications, and the idea of "well mixed" is easily discarded as intuitively obvious [START_REF] Kukukova | Measuring the scale of segregation in mixing data[END_REF]. One of the first attempts to quantify the amount of mixedness was the work of Danckwerts in 1952 [START_REF] Danckwerts | The definition and measurement of some characteristics of mixtures[END_REF] who defined the intensity of segregation for a binary mixture of liquids. While this parameter can be used to quantify the simplest mixing problems, it cannot quantify a number of important mixing processes. During the 1970s, Chemineer published the Chemscaleas, a qualitative description of the intensity of mixing in a tank and this concept was widely used for process design for many years. From the 1960s to the 1980s, Bourne, Villermaux and others developed more refined ideas about macromixing, mesomixing, and micromixing, but again the definitions are somewhat indirect. Concurrently, Corrsin (1957,1964), Toor (1969), and Brodkey (1964) all investigated the impact of turbulence on mixing through measurements of concentration fluctuations at a point, sometimes calling this the segregation, with the idea that as the variance in concentration drops to zero, the fluid approaches perfect homogeneity. In the early 1990s, Ottino examined laminar chaotic mixing with a fresh analytical perspective using chaos theory, and later in 2004 computational fluid dynamics promised numerical solutions to many complex problems as appeared in the work of E.

M. Marshall and A. Bakker [START_REF] Marshall | Computational fluid mixing[END_REF]. The most complete work on quantifying the degree of mixing has been done by Kukukovaet al. [START_REF] Kukukova | Measuring the scale of segregation in mixing data[END_REF] who proposed that the state of segregation can be characterized using three variables: (i) intensity of segregation, which is the variance in concentration throughout the vessel; (ii) scale of segregation, which is the distribution of length scales; (iii) exposure, which is the rate of change in segregation.

The intensity of segregation can be computed using a wide range of metrics, all of which are functions of either the variance or standard deviation in tracer concentration.

Blending time can be defined as the time taken to reduce the maximum variation in particle volume fraction by 95% from its initial value. The definition of this parameter will be discussed in chapter 3.

2-4-Coupling CFD and PBM

A wide range of synthetic latexes are produced using emulsion polymerization. To cite an often used adage, latexes are products by process. In other words, their properties are influenced both by the chemistry chosen to make them and by the way their production process is handled. One of the main properties of such latexes is the particle size distribution (PSD), and it is clear that changing certain reaction parameters such as characteristic mixing times (process-related parameter) with respect to characteristic reaction or coagulation times (chemistry/physico chemistry) will likely have an impact on the PSD. Considering the importance of understanding the impact of changing characteristic times in chemical processes, and in particular trying to quantify the mixing in different reactor configurations, CFD simulations can serve perfectly. Computational Fluid Dynamics (CFD) simulations can be used for this purpose since they can be used to generate detailed flow fields for a wide range of reactor configurations and operating conditions and provide information about flow characteristics such as its velocity, mixing quality, etc. (which might vary spatially specially in big reactors). This kind of information are very difficult, if not impossible, to obtain from experiments alone.

Researchers have used CFD simulations to evaluate the mixing performance of different reactor configurations [START_REF] Marshall | Computational fluid mixing[END_REF] and identify regions within the reactor where the latex will be subjected to high rates of shear [24].

As it was mentioned previously in Level-two models, one uses tools such as Population Balance Model (PBM) as the most detailed process model to keep track of the changes in the PSD during the emulsion polymerization process under a given set of process conditions (velocity, concentration, temperature, etc.) and accounts (to a certain extent)

for the impact of different particle sizes on kinetics, coagulation, viscosity, etc. When the information abstracted from CFD analysis (such as mixing characteristics, reactant gradients, or the shear rate distribution) are incorporated into a framework containing a descriptive process model, a more precise prediction of the product quality can be achieved. It gives a tool to evaluate the influence of reactor scale and configuration on the evolution of latex properties over the course of a reaction. On the other hand the information obtained from the process model, such as PBM can be used to update the flow simulations in each time step.

Some researchers coupled CFD and PBE to predict certain features (e.g. particle size) in different systems, including slurries, emulsions, and gas-liquid systems. On the other hand, very few studies employed the coupling approach to calculate the droplet/particle size distribution. Heath and Koh [START_REF] Heath | Combined population balance and CFD modeling of particle and aggregation by polymeric flocculant[END_REF] modeled the aggregation of solid particles in a slurry system by coupling CFD with PBE. The distribution of particle size was not taken into account, and only the evolution of the mean particle size with time was considered.

Agterof et al. [START_REF] Agterof | Prediction of emulsion particles sizes using a computational fluid dynamics approach[END_REF] employed this approach in a sequential way to calculate the variation size. In emulsion polymerization, Pohn et al. [START_REF] Pohn | Scale-Up of Emulsion Polymerization Reactors Part I -Development of a Model Framework[END_REF][START_REF] Pohn | Scale-Up of Emulsion Polymerization Reactors Part II -Simulations and Interpretations[END_REF][START_REF] Pohn | Development of a Computational Framework to Model the Scale-up of High-Solid-Content Polymer Latex Reactors[END_REF] used CFD to predict flow properties which were saved in multi-zones. Then a multi-zonal population balance containing nucleation, growth and coagulation was coupled to the saved flow on each zone. They used this sequential (one-way coupling) method to study the scale up of a semi-batch styrene emulsion polymerization reactor and the effect of mixing on the latex PSD. Very recently Roudsari et al. [START_REF] Roudsari | Experimental and CFD Investigation of the Mixing of MMA Emulsion Polymerization in a Stirred Tank Reactor[END_REF] used a two-way coupling approach in order to investigate the impact of different process parameters on the evolution of PSD in an emulsion polymerization process. They considered nucleation and growth as the means of PSD change but neglected coagulation phenomena. To conclude, although certain attempts have been made to use PBM/CFD coupling approach for emulsion polymerization, few have really looked at the problem in detail.

Due to the high computational requirements of the simultaneous coupling, most of these works employed a sequential approach where the CFD simulations are first realized then used in a PBM using different software, which brings a significant simplification into the solution [START_REF] Agterof | Prediction of emulsion particles sizes using a computational fluid dynamics approach[END_REF][START_REF] Srilatha | Relation between hydrodynamics and drop size distributions in pump-mix mixer[END_REF][START_REF] Pohn | Scale-Up of Emulsion Polymerization Reactors Part I -Development of a Model Framework[END_REF][START_REF] Pohn | Scale-Up of Emulsion Polymerization Reactors Part II -Simulations and Interpretations[END_REF][START_REF] Pohn | Development of a Computational Framework to Model the Scale-up of High-Solid-Content Polymer Latex Reactors[END_REF][START_REF] Silva | Implementation and analysis of numerical solution of the population balance Equation in CFD packages[END_REF][START_REF] Silva | Development and implementation of a polydispersed multiphase flow model in Open FOAM[END_REF][START_REF] Silva | Comparison of the accuracy and performance of quadrature-based methods for population balance problems with simultaneous breakage and aggregation[END_REF][START_REF] Bayraktar | Numerical aspects and implementation of population balance Equations coupled with turbulent fluid dynamics[END_REF]. In general, two-way coupling is a more complete approach than one-way coupling (where only information from CFD is given to PBE)

as it allows taking into account the effect of changes in the dispersed phase properties on the flow dynamics. One way of reducing the high computation time related to twoway coupling of CFD with PBE consists of using the method of moments, where different versions can be found in the literature: Standard Method of Moments (SMM) [START_REF] Randolph | Theory of Particulate Processes. Analysis and Techniques of Continuous Crystallization[END_REF], the Quadrature Method of Moments (QMOM) [START_REF] Mcgraw | Description of Aerosol Dynamics by the Quadrature Method of Moments[END_REF][START_REF] Marchisio | Quadrature Method of Moments for Aggregation-Breakage Processes[END_REF], or the Direct Quadrature Method of Moments (DQMOM) [START_REF] Fox | Optimal moment sets for multivariate direct quadrature method of moments[END_REF][START_REF] Marchisio | Quadrature method of moments for population-balance Equations[END_REF]. In the current work the Discrete Method will be used (also known as the classes or sectional method) [START_REF] Kumar | On the solution of population balance Equations by discretization-I. A fixed pivot technique[END_REF][START_REF] Hounslow | A Discretized Population Balance for Nucleation, Growth, and Aggregation[END_REF][START_REF] Litster | Adjustable Discretization Population Balance for Growth and Aggregation[END_REF] as it allows for a better precision of the PBS.

The objective of this work is to investigate a simultaneous two-way coupling computational framework consisting of a computational fluid dynamics simulation coupled to a population balance model to predict the impact of changing characteristic times and the effect of scale-up on the particle size distribution of a polymer latex. We will attempt to achieve these objectives using Ansys Fluent for modelling the flow field while population balance add-on module of Fluent will be used to solve the PBE. The RNG will be employed for modelling the turbulence and the mixture model will be used for multiphase modeling. For the solution of population balance equation the discrete method will be applied.

3-1-Introduction

This chapter explains the employed strategy pursued in order to model and validate the flow field inside the reactor which is composed of a stirred tank with a series of four pitched blade impellers. As described in chapter 2, the hydrodynamics inside the reactor can play an important role in the emulsion polymerization process and effect the achieved product quality. Therefore, modeling and validating the hydrodynamics inside the reactor is essential. For performing this task, Ansys package was employed to create the geometry, to mesh it, and finally to model the flow field inside the reactor. The details of the steps pursued will be presented in the following sections.

3-2-Flow in a Multiple Impeller Stirred Tank 3-2-1-Description of Reactor Geometry

The considered reactor is a cylindrical stirred tank with a truncated-cone bottom. It is equipped with a central agitation shaft upon which is mounted a series of four pitched blade impellers with a diameter of 5 cm. The tank has a total capacity of four liters with the inner diameter and the height equal to 10 and 50 centimeters respectively. Therefore, the aspect ratio of the tank when it is filled is 5; however most of the experiments and simulations were done when 3/4 portion of tank was filled, which gives an aspect ratio of around 3.8. The reactor is metallic and so the operating pressure of it can go up to 80 bars. A schematic view of the mixing system is illustrated in Figure .1.

Part Symbol Relative Dimension

3-2-3-Outline of CFD Simulation Method

Fluent package (version 15.0) was used to solve all transport equations. The transport equations were integrated using control volume method. The symmetry boundary condition was used at the liquid level. First-order upwind discretization scheme was used to calculate the face fluxes in momentum and phase transport equations. The PRESTO scheme was employed for the pressure discretization. The velocity-pressure coupling was solved using SIMPLE algorithm.

The simulations were done in two consecutive steps. First the single-phase flow was solved without solution of the species transport equation to reach a steady state condition. Reaching steady state was judged by monitoring the velocity magnitude on the planes made at the regions of high turbulence in the domain. At this moment the second step consisting of the tracer injection was started. The injection of tracer was performed using the same amount as in experiments, i.e. 5

, by patching the tracer to a sphere that was defined numerically in the domain. It should be noted that sodium thiosulfate or iodine solution do not change the properties of water, so in the simulations water was used as the substance of continuous phase and the buoyant tracer. The simulations were done only for the case of injection at the bottom of the vessel, as it is

shown in Figure 3.3 schematically. was checked by monitoring the residuals of all transport equations, the area weighted average of velocity magnitude, as well as the area-weighted average of (turbulenct kinetic energy) and (turbulent energy dissipation rate) on afore-mentioned planes.

Transient simulations were done in two different agitation rates, 100 RPM and 500 RPM. Based on the rule of thumb that a reactor is operating in the transitional regime when 10 < < 10 and fully-turbulent regime when 10 , the reactor is working in transitional regime at 100 RPM agitation rate and in turbulent regime at 500 RPM agitation rate.

3-2-3-1-Experimental Validation of Flow Simulations

A set of discoloration experiments was performed to help validate the single-phase simulations and the mesh. Tap water colored by sodium thiosulfate was used as the continuous phase. The vessel was filled with 3 liters of this aqueous solution, and then a iodine solution (potassium iodide plus iodine) was injected as a neutrally buoyant discoloring tracer through a plastic tube in the tank in three different positions. The experiments were done for different agitation rates, varying from 100 RPM to 600 RPM.

Since the main vessel is made of metal, in order to have visual observation, it was replaced by a glass vessel with the same dimensions. In each case, the flow field was allowed to become formed (i.e. agitation for 1 minute minimum before performing the injection). After the injection, the medium loses its color when the iodine and thiosulfate reaction. Mixing was judged to be complete when the aqueous phase completely changed color, and the time it took for this to occur was the experimentally measured mixing time. Finally the measured mixing time was compared to the calculated one (by simulations) through the afore-mentioned UDF. We admit that this criteria is somewhat subjective, but as we show later it can help in validation of single flow results and the hydrodynamic models chosen. Figure 3.4 shows the experimental setup in the case of 500 RPM and in its half-way of discoloration for illustration purposes. As it can be seen, it takes some time for the decoloring agent to get distributed completely, and this suggests that we indeed need a tool in order to take into account the gradient in local distribution of any species mixing in the reactor domain.

3-2-3-2 Optimization of the Mesh

As mentioned earlier, the goal of this part is validating the quality of the mesh and optimizing (i.e. minimizing) the number of required cells to be employed for the rest of the simulations in subsequent chapters. The optimal number of cells of the mesh was calculated as follows: Four different meshes with increasing level of refinement were considered. The number of cells in the meshes were, in increasing order of refinement, 743090, 1374497, 3166363 and 4452251 cells. The same set of simulations as described in last section (3-2-3) were done using these different meshes. When the additional cells did not change the calculated area-averaged velocity magnitude near the impellers and the calculated mixing time by more than 5%, the flow field was considered gridindependent. The mesh with the minimum required number of cells was selected as the optimized mesh to perform the rest of simulations.

3-3-Mathematical Model 3-3-1-Turbulence and Transition Model

In this study, in the case of turbulent regime (i.e. agitation rate = 500 RPM) RNG Model was used for turbulence modeling which is a modified version of standard model and is refined to take the swirling and rapidly strained nature of the flow under study into account. A detailed description of this model is presented in section 2-3-3-2 of Chapter 2. It should be noted that the standard Model was also tested but it was found that the RNG Model gives better results in the sense of compatibility with experimental results. In the case of transitional regime (i.e. agitation rate = 100 RPM),

Transition Model was utilized (refer to section 2-3-3-3 of Chapter 2 for the equations).

3-3-2-Species Transport

The dynamic distribution of the tracer inside the reactor was calculated by activating the species transport equation in Ansys Fluent which solves a Reynolds-averaged timedependent scalar ( ) transport equation as explained in section 2-3-3-4 of Chapter 2.

The turbulence Schmidt number ( ) was set as its default value to 0.7. is 0.7 as an empirical constant and is relatively insensitive to the molecular fluid properties. Hence, there was no need to alter the default value [1]. The diffusivity of the tracer, was set to 1.0 .

3-3-3-Mixing Time

In this work we proposed the following expression to calculate the variation in tracer concentration:

( ) = (3.1)
where is the local mass fraction of tracer that evolves until reaching to a state of uniform concentration throughout the entire vessel. The local concentration is taken from Fluent in each time step and in each cell, and ( ) gets calculated through a user defined function. By definition, for a completely homogeneous system ( ) = 1 , and the mixing time is considered to be the time required for ( ) to reach to 0.95.

3-4-Results and Discussions

3-4-1-Flow Simulations and Experimental Validation in the 4 L Tank and 100 RPM

The flow simulations and their experimental validation was done at two different agitation rates: 100 RPM and 500 RPM. These two agitation rates were chosen since they represent an upper and lower limit used in polymerization experiments [2].

3-4-1-1-Formation of Steady State Flow Field before the Injection of Tracer

Figure 3.5 shows the temporal variation of area-weighted average of velocity magnitude on the plane where the upper impeller is located. As it can be seen, the value of velocity magnitude reached a plateau at approximately 20 seconds. This also confirms that the solution is very well converged. However, the formed flow at 25 second was chosen for the rest of the simulations (injection of tracer) to ensure being at steady state.

Tracer mass fraction

3-4-2-Flow Simulations and Experimental Validation at 4 L Tank and 500 RPM

The mixing time analysis and mesh optimization were also performed at an agitation rate of 500 RPM [3].

3-4-2-1-Formation of Steady State Flow Field before the Injection of Tracer

Figure 3.12 shows the temporal variation of area-weighted average of velocity magnitude on the plane where the upper impeller is located. As it can be seen, the value of velocity magnitude reached a plateau at approximately 5 seconds. This also confirms that the solution is very well converged. Comparing this graph with Figure 3.5, it can be seen that at higher agitation rate the flow reaches steady state faster. The formed flow at 10 second was chosen for the rest of the simulations (injection of tracer) to ensure being at steady state. Figures 3. 13 As mentioned previously, the same sets of simulations have been repeated for the other levels of mesh refinement and it was found out that by increasing the cell number from 743090 to 4452251 cells we do not gain a significant improvement in mixing time prediction. In fact increasing the mesh refinement level did not changed the calculated mixing time by more than 5%, so a decision was made to continue the rest of simulations with the lightest mesh (with the less cell cumber) which contains 743090 cells.

3-4-3-Flow Simulations at 4 L Tank with a New Impeller Configuration and 500 RPM

It should be noted that the reactor configuration, and in particular the type and arrangement of the impellers, used in the simulations presented above was based on a reactor used in a related study in our laboratories [2]. Since the characteristic mixing time was poor (found to be 53 s in the reference simulations), a decision was made to replace the last pitch blade impeller with a large hydrofoil agitator mid-way through the accompanying study. In fact pitched-blade turbine is an axial-flow impeller, which generates flow along the impeller axis; but with the previous design of the impeller, the generated axial flow was not strong enough and so the quality of mixing was not excellent. We therefore did a rapid evaluation of mixing time on the new impeller configuration as well. A schematic view of this impeller is shown in Figure 3

U(t) time (s) 3-4-4-Flow Simulations at 4000 L Tank and 108 RPM

In order to investigate the impact of reactor scale-up on mixing times, CFD simulations of flow were run on a bigger scale geometrically similar reactor with 4000 L capacity corresponding to tank diameter (T) of 1 m and original agitator configuration. Based on the objective of maintaining a constant volumetric power input (note a common rule of thumb for scaling up dispersions of particles or for gas-liquid mixing is to keep P/V constant [4]), the impeller rotational speed is lowered when the vessel is scaled-up, based on the following relationships [5]:

(3.2) , = , , , (3.3)
It should be noted that Equation (3.2) holds true for turbulent conditions only, where inertial forces dominate. Applying Equation (3.3), the impeller rotational speed is calculated equal to (approximately) 108 RPM. At this case also the flow regime is fully turbulent, so the RNG Model was chosen as in the 4 L reactor.

3-4-4-1-Formation of Steady State Flow Field before the Injection of Tracer

Figure 3.22 shows the temporal variation of area-weighted average of velocity magnitude on the plane where the upper impeller is located. As it can be seen, the value of velocity magnitude reached a plateau at approximately 60 seconds. This also confirms that the solution is very well converged. However, the formed flow at 75 second was chosen for the rest of the simulations (injection of tracer) to ensure being at steady state.

3-5-Conclusions

Ansys Fluent was used in this Chapter to model the flow field inside the stirred tank.

The obtained results on mixing time are in good compatibility with the experimentally measured mixing time, which is a strong proof of the validity of the approach, the developed mesh, and the models used in order to solve the hydrodynamics of the reactor.

Considering the flow regime in the reactor which necessitates using a complex viscous model such as the "RNG Model" and "Transition " which inherits a certain level of impreciseness, this mesh and the chosen hydrodynamic models were considered valid for the rest of simulations.

4-1-Introduction

As was explained previously, in general there are three phenomena that effect the PSD of a latex:

Nucleation, by which new polymer particles are formed;

Growth of these newly formed particles due to polymerization;

Coagulation of the particles.

The objective of this chapter is to investigate a simultaneous two-way coupling computational framework consisting of a computational fluid dynamics simulation coupled to a population balance model to predict the impact of changing time scales on the particle size distribution of a polymer latex.

In order to simplify the problem, and to show the potential importance of changing the different time scales for mixing, coagulation etc. in the reactor, only the coagulation of particles was considered in the model of the evolution of the particle distribution in the current chapter, and was modeled using the DLVO approach. Thus, the PBE considered only contains a coagulation term with no reaction. It should be mentioned that that the growth rate is relatively slow and its effect on the flow dynamics during the simulation time (few seconds) may be neglected. However, nucleation is a rapid phenomenon and it will be interesting to include this in future studies.

The production of many polymers often involves stages in which the monomer, polymer, or the intermediate products are dispersed into small particles. This adds a certain number of challenges that are mainly related to the colloidal stability of the latex [1]. Unwanted particle coagulation is a major cost to the latex manufacturing industry mainly because of its negative effect on the product quality, product loss and reactor down time. On the other hand, the ability to perform controlled coagulation in other circumstances is desirable, and it is important to be able to predict the impact of mixing and species concentrations in this particular case.

A typical emulsion polymerization formulation consists of four main substances that exist in the aqueous medium: the monomer; polymer particles; the initiator that decomposes and creates the radicals; and the surfactant used to create and stabilize polymer particles. In addition, one can find buffers and other additives in smaller quantities. The effect of the ionic species, i.e. initiator and surfactant, on latex stability is often modeled through DLVO model that predicts the coagulation rate between particles. This model was explained in detail in section 2-2-3-1 of Chapter 2.

In the DLVO model, the particle coagulation rate is dependent on the local concentration of ionic species, which in turn depends on the amount of species injected, and on the quality of mixing. A limited range of tools exist in order to quantify the mixing in different reactor configuration. As was explained in Chapter 2, CFD simulations can serve perfectly for this purpose, since they can generate detailed flow fields for a wide range of reactor configurations and operating conditions and provide useful information about flow characteristics which are not easy to obtain experimentally. Contrarily, the quality and speed of mixing might be effected by changes predicted by PBE, as the relative velocity of the two phases is influenced by particle relaxation time which is dependent on the PSD. This justifies the two-way coupling considered in this work.

Different scenarios were simulated in order to study the effect of agitation rate, reactor scale, impeller type, and distribution of ionic species on the final PSD.

4-2-CFD-PBM Coupling Framework

The geometry and its meshing procedure of the reactor applied in this study which is a cylindrical stirred tank with a truncated-cone bottom equipped with a series of four pitched blade impellers was described in Chapter 3 in detail.

The Fluent package (version 15.0) was used for solving all transport equations in order to compute the single-phase flow simulation to obtain a fully formed steady state flow field before coupling the PBM to it. The symmetry boundary condition was used at the liquid level. A first-order upwind discretization scheme was used to calculate the face fluxes in the momentum and phase transport equations. The PRESTO scheme was employed for the pressure discretization. The velocity-pressure coupling was solved using SIMPLE algorithm. RNG Model was used for turbulence modeling. The time step was selected equal to 0.01s. Mixture model was used for multiphase modeling which was described in detail in section 2-3-3-3 of Chapter 2.

The number of attempts using the simultaneous two-way approach is very limited in the literature for the main reason of computational expense of this method which is very high. This is why in these few attempts the bin number was chosen to be very low (e.g.

7 in [2]). Since a two-way coupling approach in which the equations of flow and the PBE are solved simultaneously is used in this study, and as will be shown, later, the chosen number of bins is high (31 bins), the computational time is quite high (calculation time of 6 hours for simulating 1 second of reality). This is why in general the simulations were not continued for long periods of time. In order to decrease the calculation times and isolate the more important parameters to look at, a part of simulations (were the distribution of species was considered uniform) were done on a 2D mesh with 950 cells.

4-2-1-Uniform PBM Simulations

Before simulations using the real reactor geometry, optimization of the PBE parameters (i.e. discretization method, number of bins) was done using a uniform system in which the ionic species existing in the recipe are uniformly distributed in the whole domain.

As mentioned previously, in order to accelerate the calculations, this set of simulations was done on a simplified mesh consisting a 2-dimensional surface with 950 cells and initializing it with the average velocity magnitude computed in real geometry of the reactor. This system will be referred to as, "uniform simulation". A two phase system was simulated with the latex as the continuous phase and the polymer particles as the secondary particulate phase.

These simulations were used to identify the key process parameters and to guide the choice of the more complete simulations shown below. Preliminary calculations revealed that the evolution of the PSD was found to most sensitive to the surface charge density of polymer particles ( ) calculated with equation (2.26).

4-2-2-Coupling PBM to CFD to Model Latex Particle Size Distribution

The simulations were done in two consecutive steps: First a single-phase flow was defined, constituting the continuous aqueous phase. The equations of flow were solved on this case until we reached a steady state condition (fully developed flow). The establishment of steady state was judged by monitoring the area-weighted averaged velocity magnitude and the area-weighted average of (turbulent kinetic energy) and

(turbulent energy dissipation rate) on the planes made at the regions of high turbulence in the domain; when this value reached a plateau and did not change anymore, the flow was assumed to have reached steady state. On the other hand this can be considered as a means of convergence check.

Then, the second step was commenced by coupling the population balance equation to the flow field found in the first step by activating the Fluent PBM add-on module and patching the polymer particles defined as the dispersed particulate solid phase in the reactor domain. It should be noted that in the second step, the equations of flow were still being solved alongside with the population balance equation. So the framework utilized in this study is not sequential, rather it is a simultaneous framework in which all of the transport equations are being solved simultaneously with the PBM.

4-3-3-Population Balance Model

As coagulation was considered as the only means of particle diameter change in this chapter, the second term on the RHS of the equation (2.45) (related to particle growth)

A user defined function (UDF) was coded (in C language) to calculate the coagulation rate coefficient and was hooked to Fluent to be used in the PBM. In every time step, the value of is updated to calculate a new value for the surface charge density and coagulation rate coefficient (look at Equation (2.26)).

4-3-Results and Discussions

In this section, first of all, the simulation results of 3D single-phase flow are presented and discussed (so CFD simulations alone, without coupling with PBM). The second part is dedicated to the results of the preliminary set of simulations where a uniform case (thus instantaneous distribution of ionic species) was considered to optimize the PBM discretization method and evaluate the influence of major parameters on coagulation.

The third part shows the CFD-PBM coupling simulation results, where a more realistic modelling was considered to study the effect of mixing on the evolution of particle size distribution.

4-3-1-Single-Phase Flow Simulation

The considered solid content of the latex was SC = 15 weight percent (w/w), so as a first approximation the viscosity of the aqueous phase can be assumed to be equal to the viscosity of water. The density of polymer particles is set equal to 1.1 g/cm 3 .

The agitation rate was chosen equal to 500 RPM, which is the agitation rate experimentally employed in this reactor [4]. The Reynolds number is therefore: Besides, another algorithm was also tested where the discretization was evenly-spaced, and since it showed almost no difference, the geometrical algorithm was chosen which is the default case in Ansys Fluent.

4-3-2-2-Uniform PBM Simulations -Effect of Key Parameters

In this set of simulations, it was assumed that the ionic species, i.e. the initiator and the surfactant, are distributed uniformly in the reactor domain. This means that at the instant when the PBM is coupled to the flow equations, the concentration of these species is initially uniform over the domain at the same amount. As it was mentioned previously and it will be shown here, the parameter with a paramount effect on PSD was found to be , the surface charge density of polymer particles, through preliminary calculations.

Our simulations will therefore focus on this parameter, and the impact that changing relative mixing times has on the evolution of the PSD.

Based on a review of the patent literature, the initiator concentration was chosen equal to 0.4 moles/m3 and that of the surfactant was set equal to 4 moles/m3 for the surfactant (a very low level).

The initial diameter of the latex particles was set at 65 nm, and 21 bins were chosen for the PSD. The discrete method (which is a finite difference method) was used to solve the population balance equation which is more precise than the other PBM solution methods (at the price of being more computationally expensive). The time step was set equal to 0.01 s.

Figure 4.5 shows the evolution of PSD at different times with a uniform distribution of surfactant (i.e. instantaneously mixed). It is clear that with this amount of surfactant in the system, the particles coagulate from the very beginning of the simulation because the surface coverage is very low. After 1 second, two peaks can be seen in the particle size distribution, the first one is at particle diameter equal to 111 nm and the second one at particle diameter of 130 nm. After few seconds, the first peak disappears gradually and a third pick at particle diameter equal to 164 nm appears.

4-3-3-Results of Coupling of 3D Two-Phase CFD with PBM

A second case-study was considered in which the original latex were just sufficiently covered with surfactant to be stable (i.e. there is enough surfactant coverage at time 0 to stabilize the particles). We imposed a constant growth rate on the particles (5 10 -12 m 3 /s), and due to this growth, the surface area increases in such a way that the particles would become destabilized if no additional surfactant were added to the reactor. Since these simulations are quite time consuming, it was decided that we would simulate our system for a very short period of time, starting from the conditions of destabilization.

Separate calculations (not shown here for the sake of brevity) showed that we need to maintain a surface charge density of approximately 0.6 C/m 2 to ensure latex stabilization. For the simulation conditions discussed below, this means that we need a surfactant concentration of 30 mol/m 3 to avoid coagulation if the surfactant were perfectly and instantaneously mixed. However, mixing is not instantaneous; rather the injected surfactant is dispersed throughout the reactor as a function of time (as we saw in Chapter 3). For this reason, not all of the particles receive the same amout of surfactant at the same time. This means that even if we theoretically have enough surfactant to avoid coagulation, the fact that we will have variability in the local concentration of surfactant means that coagulation will occur. The idea behind these simulations is to show that the approach presented in this work allows us to understand the impact of different process changes (e.g. intensity of mixing, length scales etc.) on the variation of surfactant concentration (through generating information on the local concentration of surfactant by CFD), and thus the impact of "real world" mixing on product quality.

4-3-3-1-Effect of Surfactant Injection

Based on the scenario presented above, for particles of size 65nm and surface charge density of 0.6 C/m 2 , the results of injecting of 30 mol/m 3 surfactant (enough to stabilize the particles if mixing were instantaneous) are shown in Figures 4.7

4-3-3-2-Effect of Surfactant Injection Position

If one supposes that we wish to improve upon the scenario presented above, it is obvious that we need to improve (i.e. reduce) the characteristic time for mixing in the reactor.

There are different ways of doing this, but principally one can either stir faster of reduce the characteristic distance over which mass transfer takes pace. One way to reduce the distance over which mass transfer occurs in the reactor is to inject the surfactant at three different spots to the reactor (the total amount remains same as in the previous scenario) Therefore the number of bins chosen was set to 31 bins instead of 21 bins in order to take this normal distribution into account. While the change in the number and PSD of particles is not exceptional, it can be seen that the enhancement in surfactant distribution in the case of injection at three spots prevents a certain fraction of the particles from coagulating. The number of small particles is higher when the surfactant is injected at three spots than the case where the surfactant is injected in one spot at the top of the reactor, and the peak at 165 nm is less pronounced.

4-3-3-4-Effect of Stirring Rate

In order to study the effect of agitation rate on surfactant distribution and thus stabilization of the particles, a simulation was run with the agitation rate of 2000 RPM. This is of course not physically realistic as similar conditions in a real experiment would provoke orthokinetic coagulation, but rather serves to show the impact of decreasing the characteristic mixing times (another way to improve the mixing time would of course be to change the type and number of agitators -we saw that adding one hydrofoil improved the time -but this would make it difficult to directly compare the different simulations if we change the physical set up) . The latex was assumed to have a size distribution with an initial mean diameter of = 65 (so not all particles are put in the bin 65 nm) The surfactant was added only in one spot at the top of the reactor. Figure 4.12 shows the average particle diameter contours at the same instants and same injection position for agitations at 500 RPM (Figure 4.12-a) and 2000 RPM (Figure 4.12-b), respectively. It can be clearly seen that when the mixing time decreases due to increasing the agitation rate from 500 RPM to 2000 RPM, the particles stability is improved, which is due to their capture of the surfactant more quickly.

The calculated mixing time in these cases was equal to 13 s when the agitation rate is 2000 RPM, whereas it was 53 s for the agitation rate of 500 RPM.

4-3-3-5-Effect of Impeller Type

It should be noted that the reactor configuration, and in particular the type and arrangement of the impellers, used in the simulations presented above was based on a reactor used in a related study in our laboratories [5]. Since we noticed that the characteristic mixing time was very poor (found to be 53 s in our reference simulations), we decided to replace the last pitch blade impeller with a large hydrofoil-type agitation.

A schematic view of this impeller was shown in Figure 3 The same series of simulations presented above was done for this stirrer configuration as well. The surfactant was injected through three different spots to the reactor as in As it can be seen, changing the impeller does not appear to have a significant impact on the distribution of the surfactant and thus on the stability of the latex in the particular situations simulated here. This shows that the characteristic mixing time is not sufficient in the system at hand. It seems that in order to see a meaningful enhancement in particle stabilization by surfactant addition, the decrease in mixing time should be higher than the value this impeller configuration provides.

-0.5 0 0.5 The CFD/PBM coupling framework is used to investigate the impact of reactor scaleup on the evolution of PSD due to Brownian coagulation during emulsion polymerization process.

The main objective of scale-up is to design a large scale system that will achieve the same environment quality as in a laboratory one in order to achieve the same product quality which is usually best at the laboratory scale because of excellent mixing and heat transfer and the presence of the academically trained chemists who put a high premium on yield and selectivity. The same product quality can be achieved upon scale-up by employing many bench chemists working in parallel, but the economics will be poor [6].

Scaling-up can be practiced by maintaining some selected parameters constant, the most common of which is / (power per unit volume). The emergence of affordable computational power has led to the development of emulsion polymerization modeling frameworks that can be used to numerically investigate process scale-up, thereby potentially reducing process development times and costs [7,8]. The information abstracted from CFD analysis when input into a process model, can be used to evaluate the influence of reactor scale and configuration on the evolution of latex properties over the course of a reaction.

In order to investigate the impact of reactor scale-up, CFD simulations were run on a bigger scale geometrically similar reactor with 4000 L capacity corresponding to tank diameter ( ) of 1 m. Based on the objective of maintaining a constant volumetric power input, the impeller rotational speed is lowered when the vessel is scaled-up, based on the following relationships [24] (c.f. Chapter 3, Section 3-4-4):

(4.2) , = , , , (4.3)
It should be noted that Equation (4.2) holds true for turbulent conditions only, where inertial forces dominate. Applying Equation (4.3), the impeller rotational speed is calculated equal to (approximately) 108 RPM. At this case also the flow regime is fully turbulent, so the RNG-Model was chosen as in the 4 L reactor.

Figure 4.18 shows the temporal variation of area-weighted average of velocity magnitude on the plane where the upper impeller is located. As can be seen, the value of velocity magnitude reached a plateau at approximately 25 seconds. However, the value of 35 second was chosen for the rest of the simulations (to be coupled with the PBM) to ensure being at steady state. As it can be seen from the peaks, in bigger reactor particles are coagulating more which is reasonable considering the higher mixing time of 37s in 4000 L reactor comparing to 12 s in 4 L reactor.

5-1-Introduction

As was mentioned in chapter 2, in reality perikinetic and orthokinetic coagulation both coexist and effect the PSD evolution simultaneously. To be more precise, since Brownian motion is always present, regardless of whether the shear flow exists or not, no purely orthokinetic coagulation can actually exist. When the shear flow does exist, coagulation should be simultaneously orthokinetic and perikinetic [1]. Of course there might be conditions when orthokinetic will be dominant. [START_REF] Swift | The coagulation of hydrosols by brownian motion and laminar shear flow[END_REF] [2] analyzed the kinetics of simultaneous perikinetic and orthokinetic coagulation by assuming that the contributions are additive [2]. This approximation was easy to use, so it was adopted and discussed in subsequent studies [3][4][5]. Nevertheless, this assumption was criticized and refused by some other researchers. Van de Ven and Mason [6][7][8] tried to find asymptotic solutions of the diffusion-convection Fuchs-Smoluchowski equation, Equation (2.4), by developing a perturbation analysis of the relevant pair probability equation with respect to a parameter representing the ratio between the convective and the diffusive terms, that is, the Peclet number . It was found that shear enhances the aggregation rate and that this enhancement is proportional to the square root of the Peclet number. They concluded that the perikinetic and orthokinetic contributions are not directly additive. This analysis was limited to very small values of . Feke and Schowalter [9][10], continuing the approach proposed by Van de Ven and Mason, pointed out that the additive assumption is inconsistent with the asymptotic corrections for solid spheres with very small and very high levels of shear rate. The last technique that has been used to tackle the problem is the full numerical solution of Fuchs-Smoluchowski convection-diffusion equation. Zinchenko and Davis [11,12] studied the collision rates of spherical particles and fluid droplets by using the Fokker-Plank equation for a broad range of Peclet numbers and concluded that the additive approximation is usually questionable and a modified approximation provides better results. The semi-empirical coagulation models that have been validated experimentally require extensive parameter-fitting [13,14], therefore, their applicability to a wider range of processes and conditions is limited. Melis et al. [15] have solved Fuchs-Smoluchowski equation (Equation (2.4)) for the case of extensional flow including all interactions. They compared the results of a thorough numerical solutions to the convection-diffusion equation with the predictions of the limiting models where either orthokinetic or perikinetic coagulation dominates. Even though their approach provided useful insight, attempts to develop an analytical expression that matched the numerical data over a full range of shear rates was generally unsuccessful. However, Lattuada and Morbidelli [16] developed such an analytical expression and obtained a simplified model which matches the more thorough numerical simulations, even in the moderate shear rate range. Such an expression could feasibly be incorporated into an emulsion polymerization population balance model, although it is still more resource intensive than DLVO-only coagulation kernels. Very recently Sheng-Hua et al.

[1] examined the case of simultaneous perikinetic and orthokinetic coagulations and after refusing the additive assumption they examined the kernels proposed by Van de Ven and Mason [7] and Zinchenko and Davis [12] for small and large Peclet numbers.

They concluded that these two kernels can successfully describe the coagulation feature but only for the early stages of coagulation.

The objectives of this chapter are as follows: i) determining via CFD simulation, if shear-induced coagulation is expected to be important in the system under our study, and ii) if shear-induced coagulation is expected to become more of an issue as a reactor is scaled-up (maintaining geometric similarity and a constant power-per-unit-volume input).

Both Hong [17] and Elgebrandt et al. [18] used CFD simulation to investigate the distribution of turbulent shear inside laboratory-scale stirred tanks, but neither group extended their set of simulations to numerically investigate the shear rate distribution in lager vessels. By the way Pohn [19] used this approach to study the simultaneous perikinetic-orthokinetic coagulation phenomenon during the scale-up of a polymerization reactor.

5-2-Outline of CFD Simulation Method

The simulations in this chapter were performed on the same reactors (4L and 4000L) explained in chapter 3 and applied in chapter 4.

The simulations were done in two consecutive steps. First the single phase flow simulation was done on the reactor until reaching to steady state condition (fully developed flow) the establishment of which was judged by monitoring the areaweighted average of velocity magnitude and turbulent kinetic energy on the planes made at the regions of high turbulence in the domain; when the graphs of these two values versus time show a plateau the flow is assumed to have reached steady state. The second step was started by activating the add-on module of PBM and patching the particles into the domain on the formed flow solved and saved in the first step. It should be noted that that the equations of flow are still being solved in the second step along with the Population Balance Equation, so a two-way coupling between flow and PBM is assured.

The Fluent package (version 15.0) was used for solving all transport equations as described previously. A first-order upwind discretization scheme was used to calculate the face fluxes in the momentum and phase transport equations. The PRESTO scheme was employed for the pressure discretization. The velocity-pressure coupling was solved using SIMPLE algorithm. In this study, RNG Model was used for turbulence modeling. The time step was selected equal to 0.01s. Mixture model was used for multiphase modeling which was described in detail in section 2-3-3-3 of Chapter 2.

The population balance equation which is used to predict the evolution of PSD due to coagulation is explained in detail in chapter 4. While in that chapter the contribution of orthokinetic coagulation was ignored in the PSD evolution, the objective of this chapter is to study the importance of this contribution.

5-4-Results and Discussions

In the simulation that will be presented here, the latex particles were assumed to have an initial mean diameter of = (normal distribution) and the number of bins was set to 31. The considered solid content of the latex was SC = 15% wt. Therefore, the viscosity of the aqueous phase can be assumed to be equal to the viscosity of water.

The density of polymer particles is equal to 1.1 g/cm 3 . The agitation rate of 500 RPM was chosen in the case of 4L reactor (which is the agitation rate experimentally employed in this reactor [20]) and results in = ; therefore the flow regime is fully turbulent.

5-4-1-Case Study 1: Simulation of Pure Orthokinetic Coagulation

A first set of simulations was done to see the effect of pure orthokinetic coagulation, i.e.

in the absence of stabilizing substances like surfactant to suppress perikinetic contribution in coagulating the particles. For this case the orthokinetic kernels proposed by Saffman and Turner [21] was considered as described in Chapter 2, section 2-2-3-2-2. They proposed an expression for calculating orthokinetic coagulation coefficient based on homogeneous isotropic turbulence with particles smaller than the Kolmogorov microscale, the situation that fits the application under our study. This expression was also presented by Spielman [22] and its applicability in producing correct results was confirmed by Mei and Hu: [23] = ( + )

where is the local rate of energy dissipation and is the kinematic viscosity of the medium. As it was explained in Chapter 2, given Elgebrandt's conclusion [18], in this work we use the local shear rate rather than an average one. For calculating the coagulation rate using this equation, is calculated using the equations of flow in Fluent and fed to the kernel coded in C, and hooked to Fluent as a UDF, at each time step. As it was mentioned previously, no pure orthokinetic coagulation can exist [1].

Therefore another set of simulations was done using the simultaneous coagulation kernel represented by Equation (5.2) (refer to Equation (2.44) of Chapter 2, section 2-2-3-3) assuming ( , ) = , which is called "Combined kernel" henceforth. As it was mentioned previously in Chapter 2, the Peclet number provides a measure of the relative importance of shear versus diffusion and is defined as: [24] and defined above (500 RPM, 4 moles/m 3 surfactant) As it can be seen the contribution of orthokinetic coagulation is not important since the predicted PSD's are not that different from each other. This result is in agreement with the experimental results done on the same reactor at the same agitation rate in which the coagulation phenomena was seen to be not effected by the shear generated due to agitation. This suggests that the perikinetic coagulation is dominant in the system under the conditions of our study ( = . × ), and the stabilization brought by the Melis uses a parameter called the Extension rate which is a function of the rate of energy dissipation to study the behavior of the coagulation rate as a function of particle size. For a locally homogeneous and isotropic turbulent flow field, the mean value of is given by [25]:

‹ › = (5.6)
where is a skewness factor for the distribution of and is usually set to 0.6. He concludes that at low values of the value of decreases continuously as particle size increases. For our system, assuming constant kinematic viscosity, the average value of is equal to 44 1/s and it has the maximum value of 1095 1/s. Melis presents the graphs of the evolution of as a function of particle size at different rates of extension ( ) for two different systems with = . × as a system with a moderately thick double layer and = . × as system with a thin double layer. He concludes that for the system with a moderately thick double layer, if the extension rate is moderate, the coagulation rate decreases as the particle size increases whereas at larger extension rates the curves exhibit a minimum that can be regarded as a critical particle size beyond which the coagulation rate is favored for large particles. He suggests that this can explain the formation of few and very large pieces of coagulum, thus leading to some kind of "runaway" of the coagulation phenomenon often reported in industrial applications. On the other hand for the system with a thin double layer, he shows that the curve of versus particle diameter exhibits a monotonously decreasing behavior at all values of . The system under our study is placed in the between of the cases studied by Melis ( = . × ). Therefore in order to examine the behavior of the beyond this it does not provide more insight towards the study of shear induced coagulation. The change in the local shear distribution with reactor scale is shown in Figure 5.11. As it can be seen, the CFD simulations suggest that the maximum shear rate arising from the turbulent eddies will decrease rather than increase with reactor scale. From this result it is expected that the contribution of orthokinetic coagulation in the coagulation phenomenon decreases by increasing the scale.

5-5-Conclusions

In this chapter a coagulation kernel proposed by Blackley (Equation 5.2) was used to model the simultaneous perikinetic and orthokinetic coagulation phenomenon. CFD simulations showed that at the conditions simulated in this study, the contribution of shear-induced (orthokinetic) coagulation is not so important in the simultaneous coagulation phenomena. This result perfectly matches with the suggestion made by Melis [15] about existing a threshold for a parameter defined as modified peclet number ( ) equal to 0.1 above which the hydrodynamics of the system starts to play an important role in the coagulation phenomenon. The calculated values by simulations remain lower than this value at least for the segments of time simulated in this chapter.

Besides CFD simulations showed that the shear rate arising from the turbulent eddies decrease instead of increase with reactor scale-up maintaining geometric similarity and volumetric power input constant. This is why the contribution of shear-induced coagulation is even less pronounced in the bigger scale through using the previously explained Combined kernel.

6-1-Overall Conclusion

In this work a framework was developed in order to couple the computational fluid dynamics (CFD) simulation of a stirred tank reactor with the population balance were in good agreement with the experimental ones which justifies the meshing procedure as well as the chosen models in performing the flow simulations. These results suggest that the mixing in the 4L experimental set-up is actually very slow, and that there might be dead zones at the bottom of the reactor. It is suggested that in future studies with this reactor, the users reflect upon a choice of an improved agitation system before trying to do experiments.

In Chapter 4 the evolution of latex PSD due to "Perikinetic coagulation" (as a result of "Brownian motion" of polymer particles) was modeled using the well-known DLVO model. The Perikinetic coagulation rate coefficient was modeled using a kernel based on DLVO model which was coded in C and hooked into the Fluent PBM add-on module.

The local concentrations of ionic species, which are dependent on the time and length scale of mixing, are obtained from the solution of species transport equations, along with the rest of conservation equations in each time step. This data is fed simultaneously into the coagulation kernel in order to calculate the particle size distribution through the solution of PBM. Population balance equation was solved using the discrete method, and the particle size coordinate was discretized into 31 bins (which was chosen based on discretization study). The calculated particle size was fed simultaneously to the drag term, needed to calculate the relative velocity of the continuous and secondary (particulate) phases. The effect of different parameters (such as surface charge density on the polymer particles, the reactor performance and configuration, and the injection policy) was studied on the evolution of PSD through performing different case studies.

It was concluded that the main physical parameter affecting the Perikinetic coagulation phenomenon and so PSD evolution is the surface coverage of polymer particles, which is determined by the amount of surfactant injected to the system. In terms of the impact of the physical system on coagulation, the time scale of mixing has an important effect on the evolution of the PSD, because it determines the distribution of ionic species playing role in the Perikinetic coagulation of the particles. That's why by increasing the agitation rate, and so decreasing the mixing time in the reactor, the stabilization of particles gets enhanced. It was shown that another way of enhancing the particle stabilization is choosing another policy for the injection of surfactant into the system.

By increasing the reactor scale, the particles receive the surfactant in longer times and so the rate of Perikinetic coagulation is shown to be higher in bigger reactors.

Clearly a more indepth analysis of emulsion polymerization would be of more use.

Because of the extremely long computational times required for detailed solutions, we were constrained to use short time frames and to eliminate nucleation and particle growth from our simulations. One would expect that including nucleation in particular would have a great impact on the final form of the PSD, especially for large reactors when we can have very high local concentrations of surfactant (and any other species fed to the reactor).

Since the emulsion polymerization process under study is performed in a stirred tank reactor, a complete coagulation model should take into account the effect of flow motion as well. "Orthokinetic coagulation" is a coagulation phenomenon caused as a result of fluid motion. In Chapter 5 the simultaneous Perikinetic/Orthokinetic coagulation was studied. It was shown that for the system and the operating conditions under study, the contribution of shear-induced coagulation is not high. This result is in agreement with the experimental results done on the same reactor at the same agitation rate in which the coagulation phenomena was seen to be not affected by the shear generated due to agitation. The same comments made above concerning the simulations in Chapter 4 hold here. Having access to a computation platform that allows one to run more detailed, longer simulations would provide greater insight.

6-2-Recommendations for Future Work

In the hydrodynamic section of this work, it is recommended that other turbulent models and other meshing approaches such as multiple reference frame and sliding mesh be tested and compared to experimental results. Performing more quantitative experimental measurements can be useful in validation of the flow simulations more confidently.

More sophisticated DLVO model can be used in modeling the Perikinetic coagulation phenomenon in order to see the enhancement of model predictions. Likewise, employment of more complete Orthokinetic kernels which consider all hydrodynamic effects into account may be considered.

As mentioned above, an increase in calculation power and so simulation of the coagulation for longer segments of time will help to obtain more detailed and realistic predictions made by the proposed framework.

Nucleation phenomenon should be added to the framework components in order to have a more complete image of the polymerization process. It should be noted that considering nucleation will add certain number of differential equations to the system and so it will increase the calculation time even more.

Other solution methods for PBM, such as standard method of moments or quadrature method of moments can be tested and compared to the results of discrete method used in this work.

Faculté

  Techniques des Activités Physiques et Sportives Observatoire des Sciences de l'Univers de Lyon Polytech Lyon Ecole Supérieure de Chimie Physique Electronique Institut Universitaire de Technologie de Lyon 1 Ecole Supérieure du Professorat et de l'Education Institut de Science Financière et d'Assurances Directeur : M. F. DE MARCHI Directeur : M. le Professeur F. THEVENARD Directeur : Mme C. FELIX Directeur : M. Hassan HAMMOURI Directeur : M. le Professeur S. AKKOUCHE Directeur : M. le Professeur G. TOMANOV Directeur : M. le Professeur H. BEN HADID Directeur : M. le Professeur J-C PLENET Directeur : M. Y.VANPOULLE Directeur : M. B. GUIDERDONI Directeur : M. le Professeur E.PERRIN Directeur : M. G. PIGNAULT Directeur : M. le Professeur C. VITON Directeur : M. le Professeur A. MOUGNIOTTE Directeur : M. N. LEBOISNE To my family quality of the latexes is the particle size distribution (PSD). Modeling the evolution of PSD is usually accomplished through the addition of a set of population balance equations (PBEs) to the kinetic model. PBE provides a means of considering the contribution of different phenomena in the PSD evolution, being nucleation, growth of polymer particles by polymerization, and coagulation of particles due to Brownian or fluid motion (Perikinetic and Orthokinetic coagulation, respectively).

Figure 3 .

 3 Figure 3.1. Schematic representation of the reactor geometry and the key design dimensions Figure 3.2. Schematic view of the reactor mesh Figure 3.3. Injection of tracer at t = 0 s Figure 3.4. Schematic view of experimental setup at its half-way of discoloration Figure 3.5. Temporal variation of area-weighted average of velocity magnitude (m/s) on the plane where the upper impeller is located; Agitation rate = 100 RPM

1 .Figure 3 . 6 .Figure 3 . 7 .Figure 3 . 9 .Figure 3 . 12 .Figure 3 . 13 .Figure 3 . 14 .Figure 3 . 16 .Figure 3 . 19 .Figure 3 . 21 .Figure 3 . 22 .Figure 3 . 23 .

 1363739312313314316319321322323 Figure 3.1. Schematic representation of the reactor geometry and the key design dimensions Figure 3.2. Schematic view of the reactor mesh Figure 3.3. Injection of tracer at t = 0 s Figure 3.4. Schematic view of experimental setup at its half-way of discoloration Figure 3.5. Temporal variation of area-weighted average of velocity magnitude (m/s) on the plane where the upper impeller is located; Agitation rate = 100 RPM Figure 3.6. Temporal variation of area-weighted average of turbulence kinetic energy, (m 2 s -2 ) on the plane where the upper impeller is located; Agitation rate = 100 RPM Figure 3.7. Temporal variation of area-weighted average of rate of turbulent energy dissipation, (m 2 s -3 ) on the plane where the upper impeller is located; Agitation rate = 100 RPM Figure 3.8. (a) Vertical & (b) horizontal cross-sectional slices showing the change in the velocity magnitude profile (m/s) at 20 s; Agitation rate = 100 RPM Figure 3.9. Evolution of tracer mass fraction distribution at different instants of time & its comparison with video shots; Agitation rate = 100 RPM; Injection at the bottom of the reactor Figure 3.10. Variation of tracer concentration with time; Agitation rate = 100 RPM; Injection at the bottom of the reactor; Simulated with 1st refinement level of mesh Figure 3.11. Evolution of tracer mass fraction distribution at t = 3 minutes 58 seconds ~ ( ) = 0.95 Agitation rate = 100 RPM Figure 3.12. Temporal variation of area-weighted average of velocity magnitude (m/s) on the plane where the upper impeller is located; Agitation rate = 500 RPM Figure 3.13. Temporal variation of area-weighted average of turbulence kinetic energy, (m 2 s -2 ) on the plane where the upper impeller is located; Agitation rate = 500 RPM Figure 3.14. Temporal variation of area-weighted average of rate of turbulent energy dissipation, (m 2 s -3 ) on the plane where the upper impeller is located; Agitation rate = 500 RPM Figure 3.15. (a) Vertical & (b) horizontal cross-sectional slices showing the change in the velocity magnitude profile (m/s) at 10 s; Agitation rate = 500 RPM Figure 3.16. Evolution of tracer mass fraction distribution at different instants of time & its comparison with video shots; Agitation rate = 500 RPM; Injection at the bottom of the reactor Figure 3.17. Variation of tracer concentration with time; Agitation rate = 500 RPM; Injection at the bottom of the reactor; Simulated with 1st refinement level of mesh Figure 3.18. Evolution of tracer mass fraction distribution at t = 40 seconds ~ ( ) = 0.95; Agitation rate = 500 RPM Figure 3.19. Schematic representation of the new impeller and the key design dimensions Figure 3.20. Temporal variation of area-weighted average of velocity magnitude (m/s) on the plane where the upper impeller is located; New impeller configuration -Agitation rate = 500 RPM Figure 3.21. Variation of tracer concentration with time; New impeller configuration -Agitation rate = 500 RPM; Injection at the bottom of the reactor Figure 3.22. Temporal variation of area-weighted average of velocity magnitude (m/s) on the plane where the upper impeller is located; 4000 L -Agitation rate = 108 RPM Figure 3.23. Variation of tracer concentration with time; 4000 L -Agitation rate = 108 RPM; Injection at the bottom of the reactor

Figure 4 . 1 .Figure 4 . 5 .Figure 4 . 7 .Figure 4 . 8 .Figure 4 . 9 .Figure 4 . 10 .Figure 4 . 11 .Figure 4 . 12 .

 4145474849410411412 Figure 4.1. Temporal variation of area-weighted average of velocity magnitude (m/s) on the plane where the upper impeller is located; Agitation rate = 500 RPM Figure 4.2. Temporal variation of area-weighted average of turbulence kinetic energy, (m 2 s -2 ) on the plane where the upper impeller is located; Agitation rate = 500 RPM Figure 4.3. Temporal variation of area-weighted average of rate of turbulent energy dissipation, (m 2 s -3 ) on the plane where the upper impeller is located; Agitation rate = 500 RPM

Figure 4 . 16 .Figure 4 . 18 .Figure 4 .Figure 4 . 20 .Figure 5 . 1 .Figure 5 . 3 .Figure 5 . 5 .Figure 5 . 6 .Figure 5 . 7 .Figure 5 . 8 .Figure 5 . 10 .Figure 5 . 11 .Figure 5 . 12 .

 4164184420515355565758510511512 Figure 4.16. Particle Size Distribution at 5 s for two different impeller configurations Figure 4.17. Contours of average particle diameter at 5 s (a) new impeller configuration (b) original impeller configuration Figure 4.18. Temporal variation of area-weighted average of velocity magnitude (m/s) on the plane where the upper impeller is locatedFigure 4.19. (a) Vertical, and (b) horizontal cross-sectional slices showing the change in the velocity magnitude profile (m/s) at 35 s Figure 4.20. Particle Size Distribution at 5 s for two different reactor scales; Injection at three spots Fiure 4.21. Contours of average particle diameter at 0.5 s (a) 4L & 500 RPM. (b) 4000 L & 108 RPM; Injection at three spots Figure 5.1. PSD evolution predicted using the kernel proposed by Saffman & Turner [21]; Agitation rate = 500 RPM Figure 5.2. (a) Contours of diameter at 5s, predicted using the kernel proposed by Saffman & Turner [21] , (b) Contours of energy dissipation rate ( ); Agitation rate = 500 RPM Figure 5.3. PSD evolution predicted using the Combined kernel; Agitation rate = 500 RPM Figure 5.4. Predicted PSD at 5 s, Pure perikinetic coagulation vs. Combined Coagulation; Agitation rate = 500 RPM Figure 5.5. Vertical (a) and horizontal (b) cross-sectional slices showing the evolution in the localized shear rate profile, and diameter at 5 s (c) ; predicted using the Combined kernel Figure 5.6. Vertical cross-sectional slice showing the change in the Peclet number profile, at 5 s; predicted using the Combined kernel Figure 5.7. (a) The variation of Maximum Peclet Number and (b) Average Peclet Number by time; predicted using the Combined kernel Figure 5.8. Histogram of shear rate in the 4L reactor domain Figure 5.9. (a) The evolution of Maximum Modified Peclet Number and (b) Average Modified Peclet Number by time; predicted using the Combined kernel

  the particles occupied by a surfactant molecule (m 2 ) total area of particles per unit volume of the reactor (m 2 m -3 ) secondary phase particle's acceleration (m s -2 ) turbulent model constant local mass fraction of tracer concentration of the ionic species (mol m 3 ) blade clearance (m) particle diameter (m) diffusivity of the tracer (m 2 s -1 ) impeller diameter (m) diffusion tensor (m 2 s -1 ) charge on a single electron (C) body force (N) ( , , ) average number density function absolute velocity gradient (m/s) gravitational acceleration (m/s 2 ) generation rate of turbulent kinetic energy (Kg m -1 s -3 ) generation rate of turbulent kinetic energy due to buoyancy (Kg m -1 s -3 ) -surface distance between particles (m) [ ] concentration of monomer (mol m -3 ) [ ] concentration of monomer in the particle phase (mol m -3 ) average number of radicals per particle total particle concentration per unit volume of aqueous phase (part m -3 ) Avogadro's number (mol -1 ) agitation rate (RPM) or (1/s) pressure (Pa) power consumed (J s -1 ) particle Peclet number ratio factor for discretization of volume coordinate in PBE total rate of polymerization (mol m -3 s -1 ) total radical concentration (mol m -3 ) center-to-center distance between particles (m) Reynolds number radius of particle i (m) [ ] concentration of surfactant in the aqueous phase (mol m -3 ) xv [ ] concentration of added surfactant per unit volume of the continuous phase (mol m -direction (m/s) velocity in y direction (m/s) velocity in z direction (m/s) velocity vector in Equations ( ) fluctuating component of the velocity (m/s) steady mean velocity (m/s) mass averaged velocity of the mixture DLVO interaction potential (interparticle interactions (m s -1 ) relative velocity of phase relative to phase (m/s) mass averaged velocity (m/s) xvi , drift velocity of the secondary phase (m/s) Fuchs' stability ratio between two particles with the radii and blade width (m) ( , , ) vector of radial growth rates (m s -1 ) liquid height (m) valence of the ionic species

3 )

 3 kinematic viscosity (m 2 s -1 ) density of continuous phase (kg m -3 ) density of the mixture (kg m -3 ) particle density (kg m -3 ) viscosity of the mixture (kg m -1 s -1 ) turbulent viscosity (N m -2 s) viscosity of aqueous medium (kg m -1 s -1 ) or (Pa.s) volume fraction of the secondary phase particle relaxation time turbulence Schmidt number permittivity of water (s 2 C 2 m -1 kg -1 ) xvii vacuum permittivity (s 2 C 2 m -1 kg -1 ) relative permittivity of water zeta potential (V) particle surface potential (V) Stern layer thickness (m) total particle surface charge density (C m -2 ) particle surface charge density due to adsorption of surfactant on the particles' surface (C m -2 ) generated surface charge density due to initiator decomposition (C m -2 ) shear rate (s -1 ) dimensionless particle size ratio boundary layer thickness (m) dimensionless particle center-to-center distance stress tensor (N m -2 ) particle relaxation time (s) inverse Debye length (m -1 ) viscous energy dissipation rate per unit volume (kg m -2 s -3 ) volume fraction of the secondary phase Chapter 1

2- 3 - 1 -

 31 IntroductionComputational fluid dynamics or CFD is the science of predicting fluid flow, heat transfer, mass transfer, chemical reactions, and related phenomena by solving the mathematical equations which govern these processes using a numerical algorithm. The first industry to make extensive use of CFD was the aerospace industry, which from the 1960s onwards integrated CFD techniques into the development and manufacture of aircraft and jet engines. The availability of affordable high-performance computing hardware and the introduction of user-friendly interfaces have led to a recent upsurge of interest, and CFD has entered into the wider industrial community since the 1990s[START_REF] Versteeg | An introduction to computational fluid dynamics: The finite volume method[END_REF]. Today, CFD is used extensively in the field of chemical engineering to study a wide range of phenomena such as heat transfer, polymeric fluid flow, combustion analysis and separation and mixing.All CFD codes should contain three main elements: (i) a pre-processor, (ii) a solver and (iii) a post-processor. Pre-processing is the stage where the geometry is defined, which in fact is the computational domain, and in the next step is divided into a number of smaller, non-overlapping sub-domains (grid or mesh), then the problem is specified by selection of the physical and chemical phenomena, definition of fluid properties, and specification of appropriate boundary and initial conditions. In general, there are three distinct streams of numerical solution techniques: finite difference, finite volume and finite element and spectral techniques (an extension of finite element), each of which has its strengths and weaknesses. When the solution has been completed, a postprocessor can be used to display the domain geometry and grid, view the solution field, produce graphics and animation, and export data for further processing or as an input to an external process model. The most common discretization techniques are the finite element method and the finite volume method. Both methods begin with the integral form of the conservation equation. The finite element method uses simple piecewise functions (linear or quadratic) to describe variations in the unknown flow variables. The piecewise functions are substituted into the exact conservation equations and the residuals of the resulting approximation are minimized by multiplying them by a set of weighing functions and

  Reynolds number of a flow gives an estimation of the relative importance of inertia and viscous forces in that flow. At values below the so-called critical Reynolds number , the flow is smooth and adjacent layers of fluid slide past each other in an orderly fashion. If the applied boundary conditions do not change with time the flow is steady. This regime is called Laminar flow. At values of the Reynolds number above a complicated series of events takes place which eventually leads to a radical change of the flow character. The velocity and all other flow properties vary in a random and chaotic way. This regime is called Turbulent flow. Visualization of turbulence would reveal the presence of rotational flow structures called eddies or vortices. The largest eddies acquire energy through interaction with the mean flow and they, in turn, interact with smaller eddies, transferring energy to them before breaking-up. These smaller eddies transfer their energy to even smaller eddies. This so-called energy cascade continues until the eddies are small enough such that their kinetic energy can be effectively dissipated through viscous action. The energy cascade converts kinetic energy from the mean flow into heat at the molecular level. A substantial amount of research has been dedicated to the development of numerical methods to capture the important effects due to turbulence which can be grouped into the following three categories: Large eddy simulation (LES), Direct numerical simulation (DNS), and Reynolds-averaged Navier-Stokes equations (RANS).

  are three different Euler-Euler multiphase models available in Ansys Fluent: i) VOF Model, ii) the Mixture Model, and iii) the Eulerian Model. The VOF Model is a surface-tracking technique designed for two or more immiscible fluids where the position of the interface between the fluids is of interest. In this model, a single set of momentum equations is shared by the fluids, and the volume fraction of each of the fluids in each computational cell is tracked throughout the domain.

  scales: (i) Macromixing, which is driven by the largest scales of motion in the fluid and is characterized by the blend time in batch systems; (ii) Mesomixing, which is on a scale smaller than the bulk circulation, but larger than the micromixing scales; and (iii) Micromixing, which is the mixing on the smallest scales of motion (the Kolomogrov scale) and at the final scales of molecular diffusivity[START_REF] Marshall | Computational fluid mixing[END_REF]. A strict distinction should be considered between different mixing mechanisms according to flow regime; laminar mixing is mechanistically different from turbulent mixing. The flow for a given geometry is a continuous development from very low Reynold numbers (laminar operation) to very high ones (fully turbulent operation). In low Reynolds numbers, viscosity is the dominant phenomena affecting the flow structure, while in high Reynolds numbers changes in viscosity have no effect on process results, and the inertial forces dominate.

Figure 3 . 1 .

 31 Figure 3.1. Schematic representation of the reactor geometry and the key design dimensions

Figure 3 . 2 .

 32 Figure 3.2. Schematic view of the reactor mesh

Figure 3 . 3 .

 33 Figure 3.3. Injection of tracer at t = 0 s

Figure 3 . 4 .

 34 Figure 3.4. Schematic view of experimental setup at its half-way of discoloration

Figure 3 . 5 .

 35 Figure 3.5. Temporal variation of area-weighted average of velocity magnitude (m/s) on the plane where the upper impeller is located Agitation rate = 100 RPM

Figure 3 . 6 .Figure 3 . 7 .Figure 3 . 8 . 3 - 4 - 1 - 2 -Figure 3 .

 36373834123 Figure 3.6. Temporal variation of area-weighted average of turbulence kinetic energy, (m 2 s -2 ) on the plane where the upper impeller is located Agitation rate = 100 RPM

Figure 3 . 10 .

 310 Figure 3.10. Variation of tracer concentration with time Agitation rate = 100 RPM; Injection at the bottom of the reactor Simulated with 1 st refinement level of mesh

Figure 3 .

 3 Figure 3.11 shows the distribution of tracer mass fraction at t = 3 minutes 58 seconds

Figure 3 . 11 .

 311 Figure 3.11. Evolution of tracer mass fraction distribution at t = 3 minutes 58 seconds ~ ( ) = 0.95 Agitation rate = 100 RPM

  and 3.14 show the temporal variation of area-weighted average of , turbulence kinetic energy and , turbulent energy dissipation rate on the plane where the upper impeller is located. These figures assure one more time the convergence of flow simulations.

Figure 3 . 12 .Figure 3 . 13 .Figure 3 . 14 .Figure 3 . 3 - 4 - 2 - 2 -Figure 3 .Figure 3 . 16 .

 312313314334223316 Figure 3.12. Temporal variation of area-weighted average of velocity magnitude (m/s) on the plane where the upper impeller is located Agitation rate = 500 RPM

Figure 3 . 17 .

 317 Figure 3.17. Variation of tracer concentration with time Agitation rate = 500 RPM; Injection at the bottom of the reactor Simulated with 1 st refinement level of mesh

Figure 3 .

 3 Figure 3.18 shows the distribution of tracer mass fraction at t = 40 seconds where the

Figure 3 . 18 .

 318 Figure 3.18. Evolution of tracer mass fraction distribution at t = 40 seconds ~ ( ) = 0.95 Agitation rate = 500 RPM

Figure 3 . 19 . 3 - 4 - 3 - 1 -Figure 3 .

 31934313 Figure 3.19. Schematic representation of the new impeller and the key design dimensions

Figure 3 . 20 . 2 -Figure 3 .

 32023 Figure 3.20. Temporal variation of area-weighted average of velocity magnitude (m/s) on the plane where the upper impeller is located New impeller configuration -Agitation rate = 500 RPM

Figure 3 . 21 .

 321 Figure 3.21. Variation of tracer concentration with time New impeller configuration -Agitation rate = 500 RPM Injection at the bottom of the reactor

Figure 3 . 22 . 3 - 4 - 4 - 2 -Figure 3 .Figure 3 . 23 .

 32234423323 Figure 3.22. Temporal variation of area-weighted average of velocity magnitude (m/s) on the plane where the upper impeller is located 4000 L -Agitation rate = 108 RPM

Figure 4 . 1 .Figure 4 . 2 .Figure 4 . 3 . 4 - 3 - 2 - 4 - 3 - 2 - 1 -

 4142434324321 Figure 4.1. Temporal variation of area-weighted average of velocity magnitude (m/s) on the plane where the upper impeller is located Agitation rate = 500 RPM

Figure 4 . 5 .Figure 4 . 6 .

 4546 Figure 4.5. Evolution of PSD in different instants of time

  and 4.8. The stirring rate was chosen equal to 500 RPM as mentioned previously[4]. The initial conditions in this set of simulations are similar to the previous set, i.e. the initial diameter of the latex particles and the solid content of the latex were chosen equal to 65 nm (all particles are put in this size) and 15 % wt. The number of bins used in the PBE calculations was set equal to 21. It should be noted that for solving the distribution of surfactant in the reactor domain, the transport equation (Equation2.68) of this species was activated in the simulations with injection.

Figure 4 .

 4 Figure 4.7 shows the contours of average particle diameter on a plane at the center of

Figure 4 . 7 -

 47 a as a reference case, shows that if a perfect mixing system exists, in which the injected surfactant can uniformly be distributed (in 0 second), then the particles are prevented from coagulation and remain perfectly stable. Figure4.7-b shows that if no surfactant is injected, massive particle coagulation is observed.

Figure 4 . 7 -

 47 c shows that by injecting the needed amount of surfactant in a finite drop at the top of the domain (through patching the surfactant to a defined sphere), particles close to the injection spot are prevented from coagulating, but those farther away are not.

Figure 4 .

 4 Figure 4.8 shows the particle size distribution of the three mentioned cases after 5 s.

Figure 4 .

 4 Figure 4.8-a shows that despite the injection of the required amount of surfactant,

Figure 4 . 7 .

 47 Contours of average particle diameter at 5 s (a) Reference case, instantaneous surfactant distribution; (b) no surfactant injection; (c) surfactant injection at top of the reactor

Figure 4 . 8 .

 48 Figure 4.8. Evolution of the PSD of a latex as a function of time at 5 s for three different cases: (a) Reference case, instantaneous surfactant distribution; (b) no surfactant injection; (c) surfactant injection at top of the reactor

Figure 4 . 9 .Figure 4 .

 494 Figure 4.9-d that injection of the surfactant at three spots enhances the stabilization of

Figure 4 . 10 . 4 - 3 - 3 - 3 -

 4104333 Figure 4.10. Particle Size Distribution at 5 s for two different injection policies: Injection at top & Injection at three spots

Figure 4 .Number

 4 11 shows the PSD evolution after 5 s for these two different initializations. As it can be seen in this figure, by choosing the more realistic initial condition for the initial particle size distribution, the predicted form of PSD is smoother. This can be considered as an enhancement in the prediction of PSD evolution, -08 1.1E-07 1.6E-07 2.1E-07 2.6E-07 3.1E-07 3.6E-07 fact that having a normal distribution as the initial particles state is closer to the reality of latex systems. Coagulation is slower with a normal distribution, but the formation of larger agglomerates is a little more pronounced.

Figure 4 . 11 .

 411 Figure 4.11. Particle Size Distribution at 5 s for two different initialization conditions Injection at top

Figure 4 . 12 .

 412 Figure 4.12. Contours of average particle diameter at 5s (a) reference mixing conditions (500 RPM). (b) Rapid mixing (2000 RPM) Normal distribution initialization for both cases

Figure 4 .

 4 Figure 4.13 shows the evolution of the PSD with different stirring rates. It can be seen

Figure 4 . 13 .

 413 Figure 4.13. Particle Size Distribution at 5 s Normal distribution initialization for both cases

Figure 4 .

 4 Figure 4.14 shows that by choosing an injection at three spots, one can keep the realistic

Figure 4 . 14 .

 414 Figure 4.14. Particle Size Distribution at 5 s Normal distribution initialization for both cases

. 18 .

 18 Simply replacing the last pitched blade agitator with the new one leads to a reduction in the mixing time from 53 to 33 seconds (calculated from simulations) versus 40 seconds and 25 measured experimentally (both cases stirred at 500 RPM).

Figures 4 .

 4 Figures 4.15-a and 4.15-b show the contours of velocity magnitude at t = 10 seconds on

Figure 4 .

 4 Figure 4.9-b. The latex was assumed to have size distribution with an initial mean

Figure 4 .

 4 Figure 4.16 shows the PSD in both systems with two different impeller configurations

Figure 4 . 16 .

 416 Figure 4.16. Particle Size Distribution at 5 s for two different impeller configurations

NumberFigure 4 . 17 . 4 - 3 - 3 - 6 -

 4174336 -08 1.1E-07 1.6E-07 2.1E-07 2.6E-07 3.1E-07 3.6E-07 Contours of average particle diameter at 5 s (a) new impeller configuration (b) original impeller configuration Effect of Length Scale

Figure 4 . 18 .Figure 4 .

 4184 Figure 4.18. Temporal variation of area-weighted average of velocity magnitude (m/s) on the plane where the upper impeller is located

Figure 4 . 20 .

 420 Figure 4.20. Particle Size Distribution at 5 s for two different reactor scales Injection at three spots

Figure 5 .

 5 Figure 5.1 shows the PSD results of this simulation set and Figure 5.2 shows the

Figure 5 . 1 .

 51 Figure 5.1. PSD evolution predicted using the kernel proposed by Saffman & Turner [21]Agitation rate = 500 RPM

5 - 4 - 2 -Case Study 2 :

 5422 Figure 5.2. (a) Contours of diameter at 5s, predicted using the kernel proposed by Saffman & Turner [21] , (b) Contours of energy dissipation rate ( ); Agitation rate = 500 RPM Simulation of Combined Perikinetic/Orthokinetic Coagulation

Figure 5 . 3 .

 53 Figure 5.3. PSD evolution predicted using the Combined kernel Agitation rate = 500 RPM

Figure 5 . 5 .

 55 Figure 5.6 shows the local distribution inside the reactor using Equation (5.3) and

Figure 5 . 6 .

 56 Figure 5.6. Vertical cross-sectional slice showing the change in the Peclet number profile, at 5 s; predicted using the Combined kernel

Figure 5 . 7 .

 57 Figure 5.7. (a) The variation of Maximum Peclet Number and (b) Average Peclet Number by time; predicted using the Combined kernel

Figure 5 . 10 .

 510 Vertical cross-sectional slices showing the evolution in velocity magnitude profile (m/s) as the tank geometry is scaled from 4 L to 4000 L, maintaining a constant volumetric holding power input ( / ratio). is set to 500 RPM at the 4 L scale and to 108 RPM at the 4000 L scale. (a) 4 liters reactor and (b) 4000 liters reactor

Figure 5 . 11 .Figure 5 .

 5115 Figure 5.12 shows the histogram of shear rate in the reactor domain for the 4000 L

Figure 5 . 13 .

 513 Vertical cross-sectional slice showing the change in the Peclet number profile, at 5 s; predicted using the Combined kernel as the tank geometry is scaled from 4 L to 4000 L, maintaining a constant volumetric holding power input ( / ratio). is set to 500 RPM at the 4 L scale and to 108 RPM at the 4000 L scale. (a) 4 liters reactor and (b) 4000 liters reactor The maximum and average values of and in both scales of 4L and 4000 L are shown in Table 5.1 at two instants of time, 0 s and 5 s. It can be seen that the values of and are much lower in the 4000 L reactor which confirms the fact that the contribution of shear-induced coagulation in the coagulation phenomenon is lowered by scaling the reactor up. Besides by comparing the value of in 4000 L with the threshold suggested by Melis [15], = 0.1, upon which the contribution of orthokinetic coagulation becomes significant, it can be expected that at 4000 L reactor with the stirring rate of 108 RPM this contribution is very low at least for the periods of time simulated in this study.

  equation (PBE) to study the effect of various parameters on the performance of an emulsion polymerization process which leads to the production of a fine dispersion of polymer particles in a continuous aqueous medium. To assess the impact of imperfect mixing and the physical parameters of the system on the evolution of the latex PSD as one of the main parameters determining the final product quality, the transient simulation of flow was performed with the aid of a commercial CFD Package (Fluent® 15.0) to provide in each time step, the local concentrations of ionic species and certain hydrodynamic parameters such as turbulence dissipation rate and shear rate. This local information, which is difficult (if not impossible) to measure experimentally, is applied simultaneously by the PBE add-on module of Fluent to calculate the PSD for the next time step; thus a complete coupling between CFD and PBM is assured.In Chapter 3, after creating the 3D geometry of the reactor (which consists of a central agitation shaft upon which is mounted a series of four pitched blade impellers) using Ansys Design Modeler 15.0, Single Reference Frame (SRF) approach was applied to discretize the 3D computational domain via unstructured tetrahedral cells using Ansys Meshing 15.0. The conservation equations of flow were solved using Ansys Fluent 15.0. Mixture model was used as the multiphase model. "RNG " model was applied in the case of fully turbulent flow (associated with the agitation rate = 500 RPM) while "Transition " in the case of transitional regim (associated with the agitation rate = 100 RPM). Mixing time was measured through using a user defined function, coded in C and hooked into the Fluent. The calculated mixing time was compared to the experimentally measured mixing time obtained in a similar setup. The simulation results

Table 3 .

 3 1. Comparison of calculated mixing time with experimental results; Agitation rate = 100 RPM; Injection at the bottom of the reactor; Simulations done with 1 st refinement level of mesh

Table 3 .

 3 2. Comparison of calculated mixing time with experimental results; Agitation rate = 500 RPM; Injection at the bottom of the reactor; Simulations done with 1 st refinement level of mesh

Table 5 .

 5 1.

  .51) These two equations are called Reynolds-averaged Navier-Stokes (RANS) equations.

	They have the same general form as the instantaneous Navier-Stokes equations
	(Equations (2.47) and (2.48)), with the velocities and other solution variables now
	representing ensemble-averaged (or time-averaged) values. Additional terms which are
	called Reynolds stresses,	now appear that represent the effects of turbulence.
	=	+	2 3	2 3		(2.52)
	Eddy Viscosity Models include various types: i) One-Equation Models, like Spalart-
	Allmaras Model that solves a modeled transport equation for	, ii) Two-Equation
	Models, that solve two separate transport equations for turbulent dissipation rate and
	turbulent kinetic energy to calculate the turbulent viscosity	=	2	in which
	is a constant, iii) Standard (and SST)		Models and iv) --Transition
	Model.
	Among two-equation models, the standard		model has become the workhorse of

These terms must be modeled in order to close equation. There are two ways to close this equation: i) Eddy Viscosity Models (EVM), and ii) Reynolds-Stress Models (RSM).

In Eddy Viscosity Models which are mostly used for simple turbulent shear flows, Reynolds stresses are modeled using an eddy (or turbulent) viscosity, :

practical engineering flow calculations. This model owes its popularity to its robustness, economy, and reasonable accuracy for a wide range of turbulent flows. It is a validated turbulence model that has been used extensively in industrial engineering confined flows

  Stress Model (RSM) closes the Reynolds-averaged Navier-Stokes equations by solving transport equations for the Reynolds stresses, together with an equation for the dissipation rate. Since the RSM accounts for the effects of streamline curvature,

	the RNG	model more accurate and reliable for a wider class of flows than the
	standard	model. Both the Realizable and RNG	models have shown
	substantial improvements over the standard	model where the flow features include
	strong streamline curvature, vortices, and rotation. So it may be the best option for our
	application [50].
	Reynolds-
	The model constants	,	, , and have the following default values [57]:
	= 1.44,	= 1.92,		= 0.09,	= 1 and = 1.3
	"RNG	" and "Realizable	"models are two other options in the category of
	two-equation models. The RNG	turbulence model (which is derived from the
	instantaneous Navier-Stokes equations, using a mathematical technique called
	"renormalization group" (RNG) methods) has an additional term in its	equation
	compared to standard		model that improves the accuracy for rapidly strained flows.
	Also, the effect of swirl on turbulence is included in the RNG model which enhances
	the accuracy for swirling flows. While the standard	model is a high-Reynolds

number model, the RNG theory provides an analytically derived differential formula for effective viscosity that accounts for low-Reynolds number effects. These features make swirl, rotation, and rapid changes in strain rate in a more rigorous manner than oneequation and two-equation models, it has greater potential to give accurate predictions for complex flows, but it should be noted that it is more computationally expensive, because the numbers of equations solved in this model is higher.

Table 3 .

 3

	completely based on the experiments, the simulations predict ( ) = 0.907. This means
	that the simulations slightly overestimate the mixing time. The mixing time is calculated
	equal to approximately 3 minutes and 58 seconds from simulations whereas from the
	experiments it's measured equal to 3 minutes.

1 compares the mixing time obtained from simulations with the mixing time measured experimentally. Every time the discoloring border passes a specific level in the reactor (e.g. the level of each impeller), the time is noted from the videos (second row); then the value of ( ) is read at this time instant from the simulations (fourth row). As it can be seen from this table, when the discoloring agent gets distributed

Table 3 .

 3 1. Comparison of calculated mixing time with experimental resultsAgitation rate = 100 RPM; Injection at the bottom of the reactor, Simulations done with 1 st refinement level of mesh

	Experiments:	1 st impeller ~ 0.2 7s	2 nd impeller ~ 0.4 27 s	3 rd impeller ~ 0.6 4 th impeller ~ 0.8 whole reactor ~ 1 58 s 113 s 180 s
	Simulations:
	U(t)	0.34	0.51	0.66	0.81	0.907

Table 3 .

 3

2 compares the mixing time obtained from simulations with the mixing time measured experimentally which is made like Table

3

.1. As it can be seen from this table, when the discoloring agent gets distributed completely based on the experiments, the simulations predict ( ) = 0.904. This means that the simulations overestimate the mixing time. The mixing time is calculated equal to approximately 53 seconds from simulations whereas from the experiments it's measured equal to 40 seconds.

Table 3 . 2 .

 32 Comparison of calculated mixing time with experimental resultsAgitation rate = 500 RPM; Injection at the bottom of the reactor Simulations done with 1 st refinement level of mesh

	Experiments:	2 nd impeller ~ 0.4 3 rd impeller ~ 0.6 4 th impeller ~ 0.8 whole reactor ~ 1 3 s 7 s 24 s 40 s
	Simulations:
	U(t)	0.37	0.50	0.79	0.904

Table 5 .

 5 1. Maximum and average values of and in both scales of 4 L and 4000 L obtained from simulations as the tank geometry is scaled from 4 L to 4000 L, maintaining a constant volumetric holding power input (P/V ratio). is set to 500 RPM at the 4 L scale and to 108 RPM at the 4000 L scale

			4L	4000L	4L	4000L
	0 s	0.333	0.002	0.078	0.00056 4.81e-2 0.0003 1.12e-2 7.64e-5
	5 s	3.78	0.025	0.891	0.006	2.26e-1 0.0014 5.29e-2 0.0003

of the mean droplet size with time in an oil/water emulsion but size distribution was not obtained. Srilatha et al. [START_REF] Srilatha | Relation between hydrodynamics and drop size distributions in pump-mix mixer[END_REF] also employed this approach in a sequential manner and considered binary droplets' breakage and coalescence to calculate the size distribution in different zones of the mixing tank in two different emulsion systems. Kerdouss et al. [START_REF] Kerdouss | Two-phase mass transfer coefficient prediction in stirred vessel with a CFD model[END_REF] employed this approach on a gas-liquid (air in water) dispersion. They used discrete method for population balance modeling and Luo's model for breakage and coalescence.

They did not present any results on the bubble size distribution and only the contours of mass transfer coefficient and Sauter mean diameter were presented. Selma et al. [START_REF] Selma | Simulation of bubbly flows: Comparison between direct quadrature method of moments (DQMOM) and method of classes (CM)[END_REF] considered a gas-liquid system as well and coupled PBE with CFD to calculate the Sauter mean diameter considering the coalescence and breakage phenomena.

This approach has been used in a number of polymerization systems. In suspension polymerization, different approaches were proposed to predict the effect of mixing on the droplet size distribution. For instance, Maggioris et al. [START_REF] Maggioris | Prediction of particle size distribution in suspension polymerization reactors: effect of turbulence nonhomogeneity[END_REF] developed a twocompartment population balance model for taking into account the large spatial variations of the local turbulent kinetic energy in a high holdup suspension polymerization system, in order to predict the evolution of droplet sizes due to coalescence and breakage. They observed a reduction in the drop size when the stirring rate was increased in their CFD simulations Vivaldo-Lima et al. [START_REF] Vivaldo-Lima | Calculation of the particle size distribution in suspension polymerization using a compartment-mixing model[END_REF] used compartmentmixing model coupled to a PBM including breakup and coalescence terms to effectively estimate the PSD for a suspension polymerization system. The values for the rate of energy dissipation of each compartment were estimated from CFD. Alexopoulos et al. [START_REF] Alexopoulos | CFD analysis of turbulence nonhomogeneity in mixing vessels: A two-compartment model[END_REF] examined the same approach in an emulsion polymerization system to predict the effect of the turbulent kinetic energy and dissipation rate on the evolution of the droplet

Chapter 4 Modeling Perikinetic Coagulation in a Multiple Impeller Stirred Tank

is eliminated. To find an expression for the last term on the RHS, a particle coagulation kernel is needed. If we let account for particle formation due to coagulation, and we assume binary aggregation to be the dominant aggregation mechanism, it is possible to develop fairly simple expressions for the coagulation term: [3] ( , , ) = ( , ; , ) ( , , ) (

A particle of size is created when two particles with a total volume of undergo aggregation. A particle of size is consumed when it undergoes aggregation with any other particle in the system. The coagulation rate coefficient, ( , ; , ), represents the collisions between particles of size and .

When the equations of flow are coupled to the population balance equation, the rate of change of the external coordinates of the particles, ( , , ), is obtained from the solution of flow equations in each time step and is fed into the PBM to serve as the coupling parameter in the direction of CFD to PBM. Another important information from CFD used in the PBM in this chapter is the local concentrations of ionic species required to calculate . Note also that the shear rate calculated by CFD can be useful in calculating the coagulation rate when shear induced coagulation is taking place in the system which is the subject of Chapter 5. The particle diameter calculated by PBE is fed in each time step to the particle relaxation time (Equation 2.65) to finally calculate the relative velocity of two phases (Equation 2.64) for the flow equations.

= = 21141

where is the agitation rate and is the impeller diameter. Based on the rule of thumb that a reactor is operating in fully-turbulent regime when 10 , the considered stirring rate leads to a turbulent regime, which justifies the use of the RNG Model.

Figure 4.1 shows the temporal variation of area-weighted average of velocity magnitude on the plane where the upper impeller is located. As it can be seen, the value of velocity magnitude reached a plateau at approximately 5 seconds. However, the value of 10 seconds was chosen for the rest of the simulations (to be coupled with the PBM) to ensure being at steady state. As it can be seen the number of particles is higher with the higher amount of surfactant (as expected), and that at 6 / there is less of a tendency to form a peak at particle sizes around 165 nm.

4-4-Conclusion

The CFD/PBM coupling approach demonstrated in this work has the potential to capture information that is often difficult to isolate experimentally. By having the local concentration distribution of all existing species in the system in each instant of time, the information provided by the CFD simulation, and using this information in DLVO model that has been chosen for modeling the particle coagulation of a latex, the detailed evolution of its particle size distribution can be obtained in any instant of time. Being able to have such information is interesting from two aspects. On the one hand, it As it can be seen in Figure 5.3 by the passage of time the particles are getting more coagulated so the second peak is growing while the first one decreases.

existing surfactant in the system is sufficient to inhibit the effect of shear in coagulating the particles. As expected, the rate of shear is higher around the impellers. But since, as just mentioned, perikinetic coagulation is the dominant mechanism of coagulation, the results of the diameter contours for the pure orthokinetic case (Figure 5.2) and contours of diameters in Figure 5.5-b are different, and those in the combined case do not follow the distribution of shear rate inside the reactor.

Percentage of cell volume

Shear rate (1/s) According to Melis [15] the Peclet number does not identify the situation where convection starts to play a role, but modified Peclet number ( ) does:

At time = 0 second the maximum value of is equal to 0.0477 and the average value of is equal to 2.9×10 -4 . Figures 5.9-a and 5.9-b show the variation of the maximum and average values of peclet number during time, respectively. According to Melis [15] the threshold value for is 0.1; it means that > . is considered as a general criterion for the coagulation rate to be dependent on system fluid dynamics. Based on this threshold, the reactor fluid dynamics are not likely to play a role in the coagulation rate which justifies the results seen in Figure 5.4. As it is mentioned, the value of increases by the passage of time due to the increase in particles' size. This means that at longer times, the contribution of fluid dynamics may become more important in the coagulation rate. coagulation rate as a function of particle size and verify the probable existence of a critical particle size, the simulations need to be continued towards higher times. If the curve of versus particle diameter does not exhibit a minimum, it means that the value of coagulation rate will stay uneffected by the fluid dynamics of the system.

Unfortunately due to the high calculation times in our system (because of the large number of cells and high bin number), continuation of the simulations towards higher times is not likely to be practical.

5-4-2-Case Study 3: Simulation of Combined Perikinetic/Orthokinetic Coagulation in 4000L Reactor

The same set of simulations was run on a bigger scale geometrically similar reactor with 4000 L capacity corresponding to tank diameter ( ) of 1 m. Based on the objective of maintaining a constant volumetric power input, the impeller rotational speed is lowered to 108 RPM when the vessel is scaled-up, based on the following relationships (c.f. Chapter 3, Section 3-4-4): [26] (5.7)

The reader is referred to Chapter 4, section 4-3-3-6 to find out the procedure followed in order to scale the reactor up.