*Article* **Mathematical Modeling of the Production of Elastomers by Emulsion Polymerization in Trains of Continuous Reactors**

#### **Enrique Saldívar-Guerra 1,\*, Ramiro Infante-Martínez <sup>1</sup> and José María Islas-Manzur <sup>2</sup>**


Received: 12 October 2020; Accepted: 13 November 2020; Published: 20 November 2020

**Abstract:** A mechanistic model is proposed to describe the emulsion polymerization processes for the production of styrene–butadiene rubber (SBR) and acrylonitrile–butadiene rubber (NBR) elastomers in trains of continuous stirred tank reactors (CSTRs). A single model was used to describe both processes by choosing the proper physicochemical parameters of each system. Most of these parameters were taken from literature sources or estimated a priori; only one parameter (the entry rate coefficient) was used as an adjustable value to reproduce the kinetics (mainly conversion), and another parameter (the transfer to polymer rate coefficient) was used to fit the molecular weight distribution (MWD) experimental values from plant data. A 0-1-2 model for the number of particles and for the moments of the MWD was used to represent with more fidelity the compartmentalization effects. The model was based on approaches used in previous emulsion polymerization models published in the literature, with the premise of reaching a compromise between the level of detail, complexity, and practical value. The model outputs along the reactor train included conversion, remaining monomer composition, instantaneous and accumulated copolymer composition, the number of latex particles and particle diameter, polymerization rate, the average number of radicals per particle, average molecular weights, and the number of branches per chain.

**Keywords:** emulsion polymerization; styrene–butadiene rubber; nitrile rubber; mathematical modeling

#### **1. Introduction**

Two of the most important rubber products (copolymers) from the economical point of view—styrene–butadiene rubber (SBR) and acrylonitrile–butadiene rubber or nitrile rubber (NBR)—are produced via emulsion polymerization, mainly via the continuous cold process (5–15 ◦C) in trains of continuous stirred tank reactors (CSTRs). Some types of styrene–butadiene copolymer rubbers can also be produced via living anionic polymerization in solution, with a better control of the composition microstructure along the polymer chain, but the subject of this paper is the random copolymer produced via radical polymerization in emulsion. The world production of SBR was estimated to be around 8.3 × 10<sup>6</sup> ton/year in 2020 (based on [1]), and that of NBR was estimated to be 1.56 × 10<sup>6</sup> ton/year by 2023 [2], together representing roughly 3% of the total world polymer production (assuming a total world production of 300 × 10<sup>6</sup> ton/year).

The main use of SBR is in tires and automotive applications, while NBR is used in seals, gaskets, hoses, etc., where chemical resistance, higher strength, and lower gas permeability are needed.

As mentioned above, a large portion of these copolymers is produced worldwide via continuous cold processes that, due to the lower temperatures used, avoid undesired crosslinking during polymerization. The processes for the production of SBR and NBR have many similarities, basically differing in the formulation of the emulsion components fed to the process. Both processes are performed by free-radical emulsion polymerization in a train of CSTRs at low temperatures, and they use a redox system of initiators to maintain a flux of radicals at these low temperatures. Additionally, they both are designed to reach maximum monomer conversions between 70 and 75% to avoid crosslinking due to branching reactions and to limit the composition drift. Side streams of the monomer that are consumed faster can also be included in the process to correct the composition drift in both the SBR and the NBR processes.

Concerning the cold SBR process, early models were published by the Hamielec and MacGregor group in the 1980s. Both steady-state (Kanetakis et al. [3]) and dynamic models (Broadhead et al. [4] and Penlidis et al. [5]) were presented, capable of calculating the monomer conversion, the copolymer composition, the molecular weight averages, and the long-chain branching frequency, as well as the particle size distribution (PSD), in each reactor. In the steady-state model [3], the PSD was calculated based on the known residence time distribution of a CSTR using a convolution integral to estimate the PSD in reactors downstream of the train. Only micellar nucleation was included, considering radical entry competition between micelles and already formed particles (Smith–Ewart theory); Smith–Ewart case II kinetics were assumed with 0.5 as the average number of radicals per particle (e*n*) [6]. The molecular weight averages were calculated by using steady-state expressions derived by applying the method of moments. The dynamic model [4,5] had similar bases to the steady-state one, although it contains some differences in certain respects. Naturally, ordinary differential equations (ODEs) for the components of the redox initiation system, monomers, and converted monomers, surfactant, number of particles, and molecular weight distribution (MWD) moments per reactor were explicitly included. The treatment of the MWD equations was made by using the pseudo-homopolymer approach [7] (also known as the pseudo-kinetic [8] or apparent rate constants [9] methods). The evolution of the PSD was calculated in an approximated way via the discretization of the population by particle age (generations) and the use of several simplifying assumptions. In this way, only one differential equation was solved for the volume growth of a particle of the first generation, and this information was used to estimate the growth of the particles of any other generation. Interestingly, this model included the energy balances of the reactor and the cooling jacket, as well as an equation for a temperature controller. Smith–Ewart case II kinetics were also used in this model during the nucleation period and at low conversions, while at higher conversions, the extended fraction solution of the Smith–Ewart recurrence equation, proposed by Ugelstad et al. [10], was utilized to estimate <sup>e</sup>*n*. These models were used in several applications: to determine operation policies aimed at improving the productivity of the reactor train or the quality of the SBR produced, to manipulate the PSD, to enhance the operation stability (avoiding sustained oscillations for example), and to design start-up procedures to minimize the quantity of the off-spec product [4,5]; the dynamic model adapted for semibatch operation was later used to design midcourse correction policies for the control of molecular weight and cross-linking density in semibatch SBR reactors [11].

Years later, Gugliotta et al. [12] published a model based on the Broadhead et al. model with added detail in some respects. The main differences in Guggliota et al. model with respect to [4] were:



Similarly to the Broadhead et al. model [4], the Gugliotta et al. model [12] was applied to design policies for increasing production and transient optimization [15–17]. More recently, it was also used for the calibration of a soft sensor for the monitoring and control of an SBR production process using seven CSTRs in series. The model replaced actual plant data by simulating the operation of the process in an extended variable space [18].

Zubov et al. presented a detailed dynamic model for the SBR process in a train of CSTRs based on population balances and the method of moments for the description of the MWD, besides mass balances for other species [19]. Several non-conventional features are used in the model. Interestingly, they use an empirical relation for the average number of radicals in particles (e*n*) as a function of the initiator concentration; as a consequence, no radical entry or radical desorption terms are used. Additionally, the transport of the monomer from the droplets to the particles is modeled by using mass transfer terms, where the concentration gradients are calculated using partition coefficients for the estimation of the equilibrium concentrations. The model is validated with industrial data for conversion, molecular weight, and Mooney viscosity for the steady-state and transient case, and its usefulness is illustrated by its application to generating a grade transition procedure to minimize the production of off-spec product. More recently, Mustafina et al. developed a rather simple model to describe the copolymer composition and other characteristics of SBR, and they used it to fit the corresponding kinetic parameters to represent intrinsic viscosity, conversion, and molecular weight data [20]. Apparently, the used data came from a train of CSTRs, but a batch model was implemented to approximately represent the process. The used MWD model was too simple and did not take branching reactions into account.

Regarding the NBR case, the apparently first mathematical models for batch and semibatch industrial processes were published, almost simultaneously, by Dubé et al. [21] and by Vega et al. [22]. The Dubé et al. model was built to represent batch and semibatch reactors while considering redox initiation, the thermodynamic partitioning of monomers among phases using the partition coefficients model, and the pseudo-kinetic rate constants method for the copolymerization kinetics;e*<sup>n</sup>* was calculated by using the Ugelstad [10] expression. The model was capable of predicting conversion, the number of particles per L of water *N*p, the rate of polymerization *R*p, the copolymer composition, the number and weight average molecular weights *M*<sup>n</sup> and *M*<sup>w</sup> (by the method of moments), the branching averages, and the average particle diameter *D*p.

The model of Vega et al. [22] was also written for batch and semibatch reactors, and it was an extension of the previous model built by this group for the SBR process [12] (which was itself based on the Broadhead et al. model), with some modifications to capture the mechanistic features generated by the presence of the more water-soluble acrylonitrile (AN) monomer. The modifications of the model included the polymerization of AN in the aqueous phase, radical desorption, and homogeneous nucleation. However, some of these modifications were only applied to the polymerization rate but not to the MWD calculations, making these two parts of the model somewhat inconsistent. The latter model was further extended in [23] to describe in more detail the bivariate chain length and branching distribution in the batch case. It was also used to simulate plant data for testing an open-loop estimator based on calorimetric measurements to monitor NBR properties in a semibatch process [24]. Notice that none of the previous cited works addressed continuous NBR production processes; it was only in the works of Minari et al. [25,26] that the Vega et al. model was extended to the case of a train of CSTRs for NBR production to study the effect of intermediate feed streams of chain transfer agents (CTAs) and AN on the polymer quality [25] and to reduce transients during grade transitions [26]. More recently, the model of Vega et al. was used to evaluate a closed-loop control strategy to produce the NBR of uniform copolymer composition in a semibatch reactor [27].

Another model for NBR production in a train of CSTRs (also used for batch and semibatch operations) was proposed by Washington et al. [28,29]. This model was based on that of Dubé et al. [21] with many similarities: redox initiation, the partitioning of monomers described by partition coefficients, the use of pseudo-kinetic rate constants for the reaction kinetics, and the calculation of <sup>e</sup>*<sup>n</sup>* by using the extended fraction expression [10]. It also used an average particle size and consisted of differential equations (dynamic model). As in the case of other models for continuous reactors, it was used to study start-up procedures and the control of polymer properties (copolymer composition, molecular weight, and particle size).

More recently, the Washington et al. model was used to study, in more detail, the effect of the monomer partitioning model, the radical desorption, and the monomer soluble impurities of several polymer properties in the batch and the CSTR train case. [30] The model was also applied to study dynamic aspects of a train of eight CSTRs and the effect of different operation policies on the off-spec material produced during transient operation. [31] Scott et al. also used this model to demonstrate the application of the Bayesian experimental design technique to maximize the information that can be experimentally obtained for an NBR process carried out in a train of CSTRs [32].

In summary, there are few models for trains of CSTRs for the production of SBR ([3,4,12,19,20]) and NBR ([25,26,29]). Many of the models ([12,25,26]) stemmed from the same model of Broadhead et al. [4] with some adaptations, and others ([3,29]) are not very different either in terms of assumptions and structure.

Given the similarities of the SBR and NBR processes, it was the goal of this work to develop and present a single model that is flexible enough to represent both processes with relatively small variations and to demonstrate that the same mathematical structure could be applied with minor adjustments to both processes. Some previous models have addressed only one of the two polymerization systems at a time (not necessarily for the continuous process but for batch or semibatch modes), although, in some cases, they used essentially the same model with adaptations. We believe that it is instructive to treat the two processes with the same model simultaneously in order to emphasize their similar underlying structures. Other goals of this work were to introduce improvements in some key aspects of the models (to be specified below) and to reach a good balance between the level of detail, the ease of implementation, and the predictive power of the model. To do this, a different model is proposed that is intended to study the effect of the configuration of the CSTR's train in the production (rate and quality) of both SBR and NBR elastomers. Configuration refers to the number of reactors in the train (in the range of 7–16), as well as the possible presence and position (reactor number) of side feed streams of one of the monomers or CTAs in the reactor train.

The goal of the study was to compare steady-state operations; nonetheless, a dynamic model was written in the form of ordinary differential equations (ODEs) since, in our experience, the numerical solution of this kind of system is more robust than that of the corresponding steady-state non-linear algebraic equation system. This selection was based on the experience of many years of our group working on the mathematical modeling of industrial emulsion polymerization processes.

#### **2. Mathematical Model**

In this section, the different components of the used mathematical model are presented and briefly discussed. Since some of these components were taken or adapted from existing models—only the main features of them are provided, and their details can be found in the given references and in the Supplementary Materials (SM).

The philosophy of the present work during the definition of the different components of the model was to use only the level of complexity deemed necessary to have a reliable description of each of the involved phenomena. In some cases, as in the calculation of <sup>e</sup>*n*, which was based on a 0-1-2 particle model adapted from one by Vale and McKenna [33], this choice resulted in an increased level of detail compared with previous models. Additionally, the present model uses a more detailed calculation of the leading moments (0–2) of the MWD which, for consistency with the particle model, were also calculated using a 0-1-2 model based on a subset of the population balance model of Saldívar et al. [9] The model is not just a collection of existing pieces from other models, but the equations were re-derived and/or adapted as needed. In other cases, our choice was a simpler model than those in

other works; this was the case, for example, of the monomer partitioning model, which is quite simple but still effective. A more detailed discussion of these differences is provided in the specific sections of the paper.

The used kinetic scheme is given in Table 1. Notice that the medium part of the table shows, in detail, the terminal model on which the pseudo-homopolymer approach was based (for propagation and chain transfer to monomers and CTAs). Diffusion-controlled termination using a single coefficient was used. The lower part of the table shows the kinetics in the pseudo-homopolymer form, which is particularly useful for the MWD (moments) calculations. Notice that the internal and terminal double bond reactions only occur on the butadiene units of the dead polymer, and this is reflected in the multiplication of the corresponding rate constant by *F c* 1,*mo* , which is the accumulated molar fraction of butadiene in the polymer.


**Table 1.** Kinetics in the aqueous and particle phase.

<sup>1</sup> A pseudo-homopolymer approach with apparent rate constants is used to represent the copolymerization kinetics based on the terminal model. <sup>2</sup> A single value *k<sup>t</sup>* = *ktc* + *ktd*was used for the termination step, calculated as *k<sup>t</sup>* = √ *kt*11*kt*<sup>22</sup> while assuming diffusion-controlled termination, where *ktjj*, *j* = 1, 2, are the homotermination constants for monomers 1 and 2. <sup>3</sup> The transfer-to-polymer and double bond polymerization reactions are only considered for the MWD calculations.

The symbols used in Table 1 for the aqueous phase are: *Iw*, initiator; *Y r* 1 and *Y o* 1 , reducing agents in the reduced and oxidized states, respectively; *Y*2, secondary reducing agent; *Rw*, primary radicals in the aqueous phase; *Pw*, polymeric radicals in the aqueous phase; and *Dw*, the dead polymer in the aqueous phase. For the particle phase, the nomenclature is: *P<sup>i</sup>* , polymeric radicals of type *i* terminal unit; *M<sup>j</sup>* , type *j* monomer; *D*, the dead polymer; *P <sup>r</sup>* and *D<sup>r</sup>* , length (r) of live and dead polymers, respectively; *M*, total monomer (*M* = *M*<sup>1</sup> + *M*2); and *T*, chain-transfer agent. Notice that *P <sup>r</sup>* = *P r* 1 + *P r* 2 , where *P r i* are r-length polymeric radicals of type *i* (*i* = 1,2).

The following equations are valid for any reactor in the train, except where noted.

#### *2.1. Monomer Mass Balances*

The remaining mass of monomer *M<sup>m</sup> i* (*i* = 1,2)

$$\frac{dM\_i^m}{dt} = -F\_{\rm in} w\_{\rm Mi, in} - \frac{M\_i^m F\_{\rm out}}{M\_T} - F\_{i, \rm mu} R\_p V\_{\rm w} M\_{\rm wi} + F\_{\rm ri, Mi} \tag{1}$$

Monomer bound in the polymer chain *H<sup>m</sup> i* (*i* = 1,2)

$$\frac{dH\_i^m}{dt} = F\_{in} w\_{Hi,in} - \frac{H\_i^m F\_{out}}{M\_T} + F\_{i,ma} R\_p V\_w M\_{wi} \tag{2}$$

where *Mwi* is the molecular weight of monomer *i*; *ws*, *in* is the inlet mass fraction of species s entering the reactor; *Fi*,*ma* is the mass fraction of monomer *i* instantaneously formed in the reactor according to the Mayo–Lewis equation in mass form; *Fin* and *Fout* are the inlet and outlet total mass flows entering and exiting the reactor, respectively; *M<sup>T</sup>* is the total mass in the reactor; R*<sup>p</sup>* is the reaction rate, mol/(L s); and *V<sup>w</sup>* is the water volume in the reactor. Notice that the reaction rate is defined per liter of the aqueous phase. *Fri*,*Mi* are side feed streams (mass flow) of monomer *i* to reactor *ri*.

The Mayo–Lewis equation in molar units is:

$$F\_{i,no} = \frac{r\_1 f\_1^2 + f\_1 f\_2}{r\_1 f\_1^2 + 2f\_1 f\_2 + r\_2 f\_2^2} \tag{3}$$

where *Fi*,*mo* (*i* = 1,2) is the mole fraction of *i*-monomer units instantaneously incorporated in the copolymer; *f* <sup>1</sup> (*i* = 1,2) is the mole fraction of monomer *i* remaining in the reaction site; and *r*<sup>i</sup> (*i* = 1,2) comprises the reactivity ratios of the copolymer system. The two forms (molar and mass) of the equation are related by a simple function *f* <sup>M</sup> that converts the molar fraction into the mass fraction:

$$F\_{i,ma} = f\_{\rm M}(F\_{i,mo}, f\_{1\prime}, M\_{\rm ui}) \tag{4}$$

where *Mwi* is the molecular weight of monomer *i* (*i* = 1,2).

Notice that for practical reasons, mass, instead of molar units, is used for some of the balances, since the recipes are handled on a mass basis in industrial practice.

#### *2.2. Rate of Polymerization*

The polymerization rate can be calculated with the classical expression for emulsion polymerization:

$$R\_p = \frac{\overline{k}\_p \overline{n} \Big[M\_p] \mathbf{N}\_p}{N\_a} \tag{5}$$

where [*Mp*] is the monomer concentration in particles, mol/L; *n* is the average number of radicals per particle; <sup>N</sup>*<sup>A</sup>* is the Avogadro number, mol−<sup>1</sup> ; <sup>N</sup>*<sup>p</sup>* is the number of particles per L of water, L−<sup>1</sup> ; and *k<sup>p</sup>* is the apparent rate constant for propagation, L (mol s)−<sup>1</sup> , calculated as:

$$\overline{k\_p} = \sum\_{i=1}^{2} \sum\_{j=1}^{2} p\_i \mathcal{D}\_j k\_{p\_{ij}} \tag{6}$$

where *p<sup>i</sup>* is the probability of type *i* radical in particle (*i* = 1,2) and ∅*<sup>i</sup>* is the mole fraction of monomer *i* in particle (*i* = 1,2). The probabilities *p<sup>i</sup>* can be calculated as shown in [34].

#### *2.3. Population Balance Equations for Particles and Calculation of Rpand n*

Defining <sup>F</sup>*n*,*ri* as the number of particles with *<sup>n</sup>* radicals (*<sup>n</sup>* <sup>=</sup> 0,1,2) per liter of water, L−<sup>1</sup> , in the reactor *r<sup>i</sup>* in the train, population balance equations (PBEs) based on those postulated by Vale and McKenna [33] for a 0-1-2 system can be obtained. These equations are different from those of Vale and McKenna in two aspects: they do not consider the full-size distribution (partial derivatives with respect to size) but an average particle size instead that can increase with time and/or reactor. Another difference is that these equations include the inflow and outflow terms not present in the batch model of Vale and McKenna:

$$\frac{dV\_w F\_{0,ri}}{dt} = -\rho F\_{0,ri} + k\_{\rm des} F\_{1,ri} + \frac{2k\_t'}{\nu} F\_{2,ri} + Q\_{in,ri} F\_{0,ri-1} - Q\_{out,ri} F\_{0,ri} \tag{7}$$

$$\frac{dV\_w \mathbf{F}\_{1, \text{ri}}}{dt} = \rho \mathbf{F}\_{0, \text{ri}} - (\rho + k\_{\text{des}}) \mathbf{F}\_{1, \text{ri}} + (\rho + 2k\_{\text{des}}) \mathbf{F}\_{2, \text{ri}} + \rho\_{\text{mi}} M\_{\text{ic}} + Q\_{\text{in}, \text{ri}} \mathbf{F}\_{1, \text{ri}-1} - Q\_{\text{out}, \text{ri}} \mathbf{F}\_{1, \text{ri}} \tag{8}$$

$$\frac{dV\_w F\_{2,ri}}{dt} = \rho F\_{1,ri} - (\rho + 2k\_{\rm des} + \frac{2k\_t'}{\nu})F\_{2,ri} + Q\_{in,ri}F\_{2,ri-1} - Q\_{out,r}F\_{2,ri} \tag{9}$$

where ρ, *kdes*, and *k* 0 *t* are the radical entry, radical exit, and termination rate coefficients, respectively (see details and units below); *v* is the volume of an average particle, L; *Mic* is the micelle concentration L −1 ; *Qin*,*ri* and *Qout*,*ri* are the total in and out volumetric flows to and from reactor *ri*, respectively; and *k* 0 *<sup>t</sup>* = *kt N<sup>A</sup>* . Notice that the radical entry and exit coefficients are particle size-dependent as detailed below. The micellar nucleation mechanism is included in Equation (8); in general, the radicals in the aqueous phase may enter micelles or particles in a competitive way.

A dimensionless version of these equations, used for the numerical solution of the model, can be found in the Supplementary Materials.

From the solution of Equations (7)–(9) it is possible to calculate *n* and N*<sup>p</sup>* from the following expressions:

$$\overline{n} = \frac{\sum\_{n=0}^{2} n \text{F}\_n}{\sum\_{n=0}^{2} \text{F}\_n} \tag{10}$$

$$\mathbf{N}\_p = \sum\_{n=0}^{2} \mathbf{F}\_n \tag{11}$$

#### *2.4. Radical Entry and Exit Coe*ffi*cients*

The used specific expressions were taken from [9], but they were based on concepts previously proposed and generally accepted [35–37]. See the Supplementary Materials for specific details; Supplementary Material Equations (S1a)–(S1o).

#### *2.5. Mass Balances of Species*

In all the balance equations, *wSp*,*in* represents the mass fraction of species *S<sup>p</sup>* in the reactor feed and *MS<sup>p</sup>* is the molecular weight of *S<sup>p</sup>* (except for *M<sup>T</sup>* that is total mass in the reactor). Differential equations representing the balances are written in molar units, except where noted. They are written for the initiator *Iw*, reducing agents *Y*<sup>1</sup> and *Y*2, surfactant *S*, and chain transfer agent *T<sup>r</sup>* , and they are summarized below. More details on the assumptions and derivation of theses equations are included in the Supplementary Materials.

Initiator *Iw*:

$$\frac{dI\_w}{dt} = -k\_{d1}Y\_1^r \ I\_w / V\_w \ + F\_{\text{in}} \frac{w\_{I,\text{in}}}{M\_I} - \frac{I\_w \ F\_{\text{out}}}{M\_T} \tag{12}$$

Reducing agents *Y*<sup>1</sup> and *Y*2:

$$\frac{dY\_1}{dt} = \frac{F\_{\text{in}} w\_{y1,\text{in}}}{M\_{y1}} - \frac{Y\_1 \ F\_{\text{out}}}{M\_T} \tag{13}$$

*Processes* **2020**, *8*, 1508

$$\frac{dY\_2}{dt} = \frac{F\_{\text{in}}W\_{y2,\text{in}}}{M\_{y2}} - \frac{Y\_2 F\_{\text{out}}}{M\_T} - k\_{d2} Y\_1^o Y\_2 / V\_w \tag{14}$$

For further details on the use of these equations, see Supplementary Material Equations (S1a–S1r) in the Supplementary Materials.

Surfactant S:

$$\frac{d\mathbf{S}}{dt} = \frac{F\_{\text{in}} w\_{\text{s,in}}}{M\_{\text{s}}} - F\_{\text{out}} \frac{\mathbf{S}}{M\_{\text{T}}} \tag{15}$$

The surfactant balance is connected and solved simultaneously with the corresponding adsorption equilibrium (Langmuir isotherm) described by the Supplementary Material Equations (S1u)–(S1x) in the Supplementary Materials.

Given the total amount of surfactants by the mass balance, the partitioning of the surfactant adsorbed in particles (*Sa*), and free in water (*S<sup>f</sup>* ) can be defined; *S<sup>f</sup>* is given by a quadratic equation and includes the surfactant in solution and the surfactant in micelles. The micelle concentration is obtained by Equation (16) or (17):

$$M\_{\rm ic} = \frac{\left(\mathcal{S}\_f / V\_w - \, [\mathcal{S}]^{cm}\right) \mathcal{N}\_A a\_{\rm cm}}{4\pi \, r\_m^2} \text{ if } \mathcal{S}\_f / V\_w \ge \left[\mathcal{S}\right]^{cmc} \tag{16}$$

$$M\_{\rm ic} = 0 \text{ if } \mathbb{S}\_f / V\_w < \ [\mathbb{S}]^{\rm cunc} \tag{17}$$

where [*S*] *cmc* is the critical micelle concentration of the surfactant (mol L−<sup>1</sup> ), *aem* is the surface area per surfactant molecule, and *r<sup>m</sup>* is the radius of a micelle.

Chain transfer agent (Tr):

$$\frac{dT\_r}{dt} = \frac{F\_{\text{in}} w\_{Tr, \text{in}}}{M\_{Tr}} - \frac{F\_{\text{out}} T\_r}{M\_T} - k\_{\text{tr}T} [T\_r]\_p \; \frac{\sum\_{n=1}^2 nF\_n}{N\_A} \; V\_w + \frac{F\_{\text{ri}, Tr}}{M\_{Tr}} \tag{18}$$

where *Fr*,*Tr* is a side feed stream (mass flow) of CTAs to reactor *ri*. The partitioning equilibrium of the CTAs defines their concentration in the particles. This is based on the assumption that the component partitions between the particles and the aqueous phase according to the partition coefficient are defined as:

$$\mathcal{K}\_T = \frac{[Tr]\_p}{[Tr]\_w} \tag{19}$$

where [*Tr*] *p* and [*Tr*]*<sup>w</sup>* are the concentrations of CTAs in the particle and aqueous phase, respectively. The simultaneous solution of the mass balance and the CTA partitioning is described in detail by

Supplementary Material Equations (S1y)–(S1aa) in the Supplementary Materials.

Water (in mass units):

$$\frac{d\mathcal{W}}{dt} = F\_{\text{in}} w\_{\text{w1 in}} - \frac{F\_{\text{out}} \,\mathcal{W}}{M\_T} \tag{20}$$

Differential mass balance equations were also derived for primary radicals *R<sup>w</sup>* and polymeric radicals *P<sup>w</sup>* (aqueous phase), but the quasi-steady state approximation (QSSA) was assumed for these species, resulting in a simple quadratic equation for *Pw*:

$$-k\_{\rm t} \, P\_w^2 / V\_w - k\_{mp} 4\pi \, r^2 \, \text{P}\_w \sum\_{n=0}^2 \, \text{F}\_n - k\_{mm} a\_{\rm m} M\_{\rm i\bar{x}} \, \text{P}\_w + \frac{K\_{\rm des}}{N\_{\rm A}} \sum\_{n=0}^2 \, \text{F}\_n + k\_{\rm d1} f Y\_1^{\tau} \, I\_w / V\_w = 0 \tag{21}$$

where *kmp* is an entry coefficient for radicals in particles (m/s −1 ); *r* is the radius of an average particle; *<sup>k</sup>mm* is the entry coefficient for radicals in micelles; *<sup>a</sup><sup>m</sup>* is the total micellar surface area; and *<sup>M</sup>ic* (L−<sup>1</sup> ) is the micelle concentration.

#### *2.6. Monomer Partitioning*

The modeling approach used here was relatively simple but still effective. Madhuranthakam and Penlidis [30] concluded that the use of different models for monomer partitioning (e.g., partition coefficients and models based on equations of state) does not make a significant difference in the case of the more complex NBR system in which the water solubility of AN is important. For the simpler SBR case, an equipartition approach (similar proportions of the monomers in the droplets and particle phases) equivalent to that proposed by Dougherty [38], is used. Additionally, the approach used here allows for the a priori estimation of key parameters of the partitioning model based on the known physicochemical properties of the system components. The used modeling approach was based on the simple limit conversion for saturation, *xsat*, used by Gardon [39] as the conversion at which the monomer droplets disappear (end of interval 2). Up to this point, in the presence of excess monomers, the polymer particles are saturated with monomers and maintain a nearly constant concentration ([*M*] *p* ).

[*M*] *p* is calculated differently depending on the stage (interval) of the polymerization. The decision is made based on the equivalent conversion, *x*, defined as:

$$\chi = \frac{P\_1^m + P\_2^m}{M\_1^m + M\_2^m - M\_{2w} + P\_1^m + P\_2^m} \tag{22}$$

where *M*2*<sup>w</sup>* stands for the mass of monomer 2 in water. This latter quantity is nearly zero for SBR, but it is significant for NBR due to the partial solubility of AN in water.

If *x* ≤ *xsat* (intervals 1 and 2), [*M*] *p* is calculated by the Equation (23). The details of how to arrive to this equation and how to calculate each of the involved quantities can be found in the Supplementary Materials Section 2, Supplementary Material Equations (S2a)–(S2c) of the Supplementary Materials.

$$[M]\_p = \ [M]\_{psat} = \frac{(M\_1/M\_{w1}) + (M\_2/M\_{w2})}{V\_p'} \tag{23}$$

where *V* 0 *p* is the volume of this hypothetical particle, which can be calculated assuming volume additivity:

$$V\_p' = \frac{M\_1}{\rho\_1} + \frac{M\_2}{\rho\_2} + \frac{P\_1}{\rho\_{p1}} + \frac{P\_2}{\rho\_{p2}}\tag{24}$$

where *Mwi* = molecular weight of monomer *i* and ρ*<sup>i</sup>* and ρ*pi* (*i* = 1,2) are the densities of monomer *i* and homopolymer *i*, respectively.

If *x* > *xsat* (*mass conversion x*, interval 3), [*M*] *p* is calculated assuming that all the remaining monomers, except for the possible presence of AN in the water phase (NBR case), are in the particles:

$$[M]\_p = \frac{\left(M\_1^m / M\_{w1}\right) + \left(\left(M\_2^m - M\_{2w}^m\right) / M\_{w2}\right)}{V\_p} \tag{25}$$

where:

$$V\_p = \frac{M\_1^m}{\rho\_1} + \frac{M\_2^m - M\_{2w}^m}{\rho\_2} + \frac{H\_1^m}{\rho\_{p1}} + \frac{H\_2^m}{\rho\_{p2}} \tag{26}$$

and the superindex *m* indicates mass units.

An approximate way of estimating the amount of AN in the water phase (*M*2*w*) in the NBR case is based on a partition coefficient, defined as the ratio of mass concentrations of monomer 2 in the particles and the water phase:

$$\chi\_{AN} = \frac{\left[\left(M\_2^m\right]\_p\right]}{\left[\left(M\_2^m\right]\_w\right]}\tag{27}$$

As shown in Supplementary Materials Section 2 (see Supplementary Material Equation (S2d)), *M*2*<sup>w</sup>* can be obtained by solving the following quadratic equation:

$$\frac{\chi\_{AN}}{\rho\_2} M\_{2w}^2 - \left[ \chi\_{AN} \left( \frac{M\_1^m}{\rho\_1} + \frac{M\_2^m}{\rho\_2} + \frac{P\_1^m}{\rho\_{p1}} + \frac{P\_2^m}{\rho\_{p2}} \right) + V\_w \right] M\_{2w} + M\_2^m V\_w = 0 \tag{28}$$

#### *2.7. Total Mass Balance*

Since there might be side-feed streams of monomers and/or CTAs along the train, a total mass balance per reactor (in mass units) is necessary:

$$\frac{dM\_{T,ri}}{dt} = F\_{in,ri} - F\_{out,ri} + F\_{ri,tr} + F\_{ri,M1} + F\_{ri,M2} \tag{29}$$

where *MT*,*ri* is the total mass present in reactor *ri* (reactors of different volumes can be considered along the train). Notice that in previous equations, the second subindex *ri* was omitted when it was clear that the balances were written for any reactor in the train. The same convention is used for the total mass flows entering and leaving the reactor (*Fin*,*ri* and *Fout*,*ri*, respectively). As an approximation, it is assumed that the time derivative is zero; that is, a constant mass is maintained in each reactor. In reality, a constant volume operation is used instead (level control); nonetheless, since the density variations are relatively small, this approximation introduces little error during the transient calculations and no error when the steady-state is achieved. In this way, the total output mass flow can be explicitly calculated from Equation (29) because the other terms are known (naturally, *Fin*,*ri* = *Fout*,*ri*−<sup>1</sup> for *ri* > 1).

#### *2.8. Molecular Weight Distribution*

The quantities of two different distributions, one for the live polymer and the other for the dead polymer, are defined as follows:

*Nr <sup>n</sup>* = number of length-*r* radicals in particles having *n* radicals, per L of water.

*Dr <sup>n</sup>* = number of length-*r* inactive polymer chains radicals in particles having *n* radicals, per L of water.

For the assumptions and details of the derivation of the corresponding PBEs, the reader is referred to [9]. Though the derivations used here follow the general assumptions in [9], the model was adapted to comply with the 0-1-2 scheme. Some of the terms in the balances represent a shift in the classification of radicals in the presence of a given kinetic or physicochemical event.

A general balance for the live polymer *r* = 1, 2, 3 . . . and *n* = 1, 2, 3 . . . , in any reactor, where the identification of the reactor number (*ri*) is kept only for the inlet and outlet mass flows, is:

*d V<sup>w</sup> N<sup>r</sup> n dt* <sup>=</sup> *<sup>k</sup>pNr*−<sup>1</sup> *n* [*M*] *<sup>p</sup>V<sup>w</sup>* <sup>−</sup> *<sup>k</sup>pN<sup>r</sup> n* [*M*] *<sup>p</sup>V<sup>w</sup>* <sup>−</sup> *<sup>k</sup>trMN<sup>r</sup> n* [*M*] *<sup>p</sup>Vw*+ *kdes n N<sup>r</sup> n*+1 −*kdes N<sup>r</sup> n* (*n* − 1) − *KtrTN<sup>r</sup> n* [*Tr*] *<sup>p</sup>Vw*+ρ*N<sup>o</sup>* δ (*n* = 1) δ (*r* = 1) + ρ*mic Mic* δ(*n* = 1)δ(*r* = 1) − *kdes Nn*δ (*r* = 1)δ(*n* = 1)+ρ *N<sup>r</sup> n*−1 − ρ*N<sup>r</sup> <sup>n</sup>*/*n* + *ktrTNn*[*Tr*] *<sup>p</sup> Vw*δ(*r* = 1)+*ktrMNn*[*M*] *<sup>p</sup>Vw*δ(*r* = 1) + *ktVw* 2*NAv* h *n*(*n* + 1) *N<sup>r</sup> n*+2 − *n*(*n* − 1)*N<sup>r</sup> n* i + ρ*N<sup>r</sup>* 2 2 δ (*n* = 1) +*kdbF*<sup>1</sup> h − P<sup>∞</sup> *j*=1 *jDj n Nr <sup>n</sup>* + P*r j*=1 *jDj <sup>n</sup>N r*−*j n* i *nVw <sup>N</sup>nNA<sup>v</sup>* + h <sup>−</sup>*ktrp*P<sup>∞</sup> *j*=1 *jDj n Nr n* +*ktrprD<sup>r</sup> n* P<sup>∞</sup> *<sup>j</sup>*=<sup>1</sup> *N j n* i *nVw <sup>N</sup>nNA<sup>v</sup>* + *<sup>F</sup>in*,*ri<sup>w</sup> f <sup>w</sup>N r*, *f <sup>n</sup>* /ρ*<sup>w</sup>* <sup>−</sup> *Nr <sup>n</sup>VwFout*,*ri M<sup>T</sup>* (30)

where δ (*j* = *j*0) represents the Kronecker delta function, which is 1 when a given integer variable *j* = *j*<sup>0</sup> and is 0 elsewhere. For a dead polymer, the balance is, in principle, valid for *r* = 1, 2, 3 . . . and = 0, 1, 2, 3 . . .; the resulting PBE is:

*dVwD<sup>r</sup> n dt* = ρ *Dr n*−1 − *D<sup>r</sup> n* + *kdes*(*n* + 1)*D<sup>r</sup> n*+1 − *kdesnD<sup>r</sup> <sup>n</sup>* + *VwktrM*[*M*] *pN<sup>r</sup> n* +*VwktrT*[*Tr*] *pN<sup>r</sup> <sup>n</sup>* + *Vw* 2*NAv kt*(*n* + 2)(*n* + 1)*D<sup>r</sup> n*+2 − *Vw* 2*NAv k<sup>t</sup> n*(*n* − 1)*D<sup>r</sup> <sup>n</sup>*+ *ktdV<sup>w</sup>* 2*NAv* (*n* + 1) *Nr <sup>n</sup>*+<sup>2</sup> + *ktcVw* 2*NAv* (*n*+1) *Nn*+2 *r*P−1 *j*=1 *N j n*+2 *N r*−*j n*+2 + <sup>−</sup>*ktrprD<sup>r</sup> n* P∞ *j*=1 *N j <sup>n</sup>* + *ktrpN<sup>r</sup> n* P∞ *j*=1 *jDj n nVw NnNAv* −*kdbF*1*rD<sup>r</sup> n* P∞ *j*=1 *N j n nVw <sup>N</sup>nNA<sup>v</sup>* + *<sup>F</sup>in*,*ri<sup>w</sup> f <sup>w</sup>D r*, *f <sup>n</sup>* /ρ*<sup>w</sup>* <sup>−</sup> *Dr <sup>n</sup>VwFout*,*ri M<sup>T</sup>* (31)

where *w f <sup>w</sup>* is the mass fraction of water at the inlet and *P r*, *f <sup>n</sup>* and *D r*, *f <sup>n</sup>* are the values of the distribution in the inlet stream for the live and dead polymers, respectively.

Since a 0-1-2 model was to be used, Equation (30) was applied only to *n* = 1,2, while Equation (31) was applied to *n* = 0,1,2. Then, the method of moments was applied to derive equations for the following moments:

$$\mu\_{K,n} = \sum\_{r=1}^{\infty} r^K n = 1, \text{ 2, } K = 0, \text{ 1, 2} \tag{32}$$

$$\lambda\_{K,n} = \sum\_{r=1}^{\infty} r^K D\_n^r \text{ } n = 0\text{, 1, 2, } K = 0\text{, 1, 2} \tag{33}$$

where µ*K*,*<sup>n</sup>* and λ*K*,*<sup>n</sup>* are the *K*-th moments of the live and dead polymers, respectively, in particles having *n* radicals. Only the first three moments (0, 1, and 2) of both populations are needed to calculate the number- and weight-average molecular weights. The final equations for these moments, six for the live polymer (Supplementary Material Equations (S3d)–(S3i)) and nine for the dead polymer (Supplementary Material Equations (S3j)–(S3r)), are provided in the Supplementary Materials, both in their original and dimensionless versions. The modeling approach followed here involved much more detail than previous modeling efforts since a full set of moments (0–2) had to be calculated for particles with 0, 1, and 2 radicals and then added together, taking compartmentalization effects into account more faithfully.

#### *2.9. Number of Branches*

The number of branches, *B*, per L of water in any reactor, can be estimated simply by the following ODE:

$$\frac{dBV\_w}{dt} = \frac{F\_{\text{in}}w\_w^f B^f}{\rho\_w} - \frac{F\_{\text{out}}V\_w B}{M\_T} + \frac{k\_{lrp} + k\_{db}}{N\_A V\_p} \sum\_{n=1}^2 n\lambda\_{1,n} V\_w \tag{34}$$

In steady-state, this can be written as:

$$0 = \frac{F\_{in}B^f}{M\_T} - \frac{F\_{out}B}{M\_T} + \frac{k\_{trp} + k\_{db}}{N\_A V\_p} \sum\_{n=1}^2 n\lambda\_{1,n} \tag{35}$$

or, equivalently:

$$B = \frac{F\_{in}B^f}{F\_{out}} + \frac{M\_T}{F\_{out}} \frac{k\_{trp} + k\_{db}}{N\_A V\_p} \sum\_{n=1}^{2} n \lambda\_{1,n} \tag{36}$$

#### *2.10. Strategy of Solution*

The model described in the previous section resulted in a first set of 13 ODEs per reactor (*ri*) associated to the states *M*1, *M*2, *H*1, *H*2, *S*, *W*, *Tr*, *Iw*, *Y*1, *Y*2, *F*0,*ri*, *F*1,*ri*, and *F*2,*ri*, as well as 15 additional ODEs per reactor for the moments of live (6) and dead polymers (9) (Supplementary Material Equations (S2d)–(S2r)). Note, however, that the equations for the live polymer moments µ0,1 and µ0,2 did not need to be solved, since equivalent information could be obtained from the solution of Equations (8) and (9) by realizing that µ0, 1 ≡ *F*1,*ri* and µ0, 2 ≡ 2*F*2,*ri*. Though the number of branches *B* could be dynamically calculated using Equation (34), it was decided to calculate its value only at the final steady-state by using Equation (36). In summary, the model resulted in a set of 26 ODEs per reactor. A summary of the ODEs and the most important auxiliary equations, which are the set of working equations implemented in the solution code, is included in Table 2.

To reduce the stiffness of the ODE system, the QSSA was applied to the remaining 4 live polymer moments. Simple explicit algebraic expressions were obtained for these variables, reducing the effective number of ODEs per reactor to 22.

Given the relative complexity of the equations (especially those related to the MWD) and to facilitate the programming and debugging steps of the implementation of the numerical solution of the ODEs, this task was implemented in two stages. In the first stage (computer program), only the first 13 ODEs (for each reactor) were integrated by using a given set of initial conditions in all the reactors in the chain until the steady-state was reached. The set of initial conditions was designed to reach the steady-state fast; linear profiles of the species along the reactor train, according to the approximate profiles expected at the steady-state, were used as the initial conditions. The evolution of the states was recorded at regular time intervals in a file, as illustrated in Figure 1. The variables in the first set of ODEs are independent of the MWD, and, therefore, the corresponding equations could be solved first. Once the solution of this set was obtained, the set of ODEs corresponding to the MWD was numerically solved in a dynamic way in a second computer program by using the dynamic information of the states obtained from the previous integration stage. To obtain more accuracy in the second stage, linear interpolation was used as needed to obtain values of the variables in the first set at integration times between the recorded time intervals (every ~20 s). The numerical routine DDASSL (Differential-Algebraic System Solver) [40] was used in the two stages (programs) for the integration of the equations, which were coded in Fortran. DDASSL is capable of solving either a system of differential-algebraic equations or a system of pure ordinary differential equations, and, since the present system is numerically a system of ODEs (all the coupled algebraic equations are explicit), any ODE integrator for stiff equations could have been used. In this case, all the variables were declared as double-precision in Fortran, and the used numerical relative tolerance was 10−7–10−<sup>6</sup> . These numerical parameters allowed for safe and robust convergence in all cases. The execution of each stage took ~1–2 s from time zero to the achievement of the steady-state in a standard laptop, even in the case of simulations for the maximum number of implemented reactors (16 in a train). The simulated time needed to achieve the steady-state was around 2–3 times the overall residence time of the reaction in the reactor train, which is usually close to 5 h. The implemented solution strategy turned out to be efficient, robust, and relatively easy to debug.


**Table 2.** Working equations for implementation in the solution code.

**Figure 1.** Data flow and computer programs for the implementation of the numerical solution of the mathematical model. Each stage (darker blocks) represents a computer program that is run separately and sequentially. See the text for more details.

#### **3. Results**

#### *3.1. Parameter Values*

The first problem to solve for the practical use of the developed computer programs was the selection of the physicochemical parameters of the simulated systems (SBR and NBR). Many of them were found in literature sources, but in several cases, the information was scarce or nonexistent. In the latter cases, most of the parameters were estimated a priori based on reasonable assumptions or indirect information; for example, the kinetic constants of the redox initiation system were estimated a priori based on data on the concentration of initiator at the last reactor exit. Notice though that these should not be considered adjustable parameters, since they were not fitted after comparison with the experimental data corresponding to the model outputs but estimated before the simulations were performed. A possible concern with this approach is how significant is the effect of these parameters on the final calculations and, therefore, if the error introduced with these estimations is significant. To answer this, the discussion can be better organized by classifying the prior-estimated parameters into three categories that are discussed next:


Essentially, only one parameter was left as an adjustable value (entry rate coefficients to micelles and particles which were assumed equal, *kmm* = *kmp*) to fit the kinetic data (mainly conversion), and a second parameter (transfer to polymer rate coefficient *ktrp*) was used to fit the MWD dispersity. This was justified since the value of the last parameter is usually model-dependent, as it is quite difficult to estimate it from independent experiments; see, for example, [28]. Since both the *ktrp* and the *kdb* affect the MWD dispersity and the number of branches *B* but it is practically impossible to discern their individual values based solely on the experimental values of the dispersity and *B*, it was decided to arbitrarily set *kdb* = 0 and to fit only the value of *ktrp*, which should then be considered as an effective value. In principle, it should be possible to measure by <sup>13</sup>C NMR techniques the number of trifunctional and tetrafunctional branches, mainly associated with transfer to polymer (one tertiary carbon) and internal double bond (two tertiary carbons) reactions, respectively, and from these values to independently obtain a rough estimate of both kinetic coefficients; however, this fell out of the scope of the present work. Regarding propagation and termination coefficients, it is well-known that the most reliable values for these parameters are those obtained by pulse laser polymerization (PLP) (when available), and there have been some reported for these coefficients for styrene and AN (see for example [41,42] for propagation coefficients); however, it was also interesting to compare the model results with kinetic parameters that had been used before by other groups working with emulsion systems, in particular for styrene polymerization, [43] or specifically for NBR polymerization, [24,28], even though they were not determined by PLP. Both types of parameters were tested in the present model; however, for comparison with previous work, the values shown in the simulations were not necessarily obtained with PLP-determined parameters (see details in Table 3). Nonetheless, the sensitivity of the model outputs to the used parameter type (determination method) was quite low. For example, for the SBR system, by changing the value of the propagation coefficient *kp*,22 for styrene from the non-PLP parameter value (151 L mol−<sup>1</sup> s <sup>−</sup><sup>1</sup> at 10 ◦C) [43] to the PLP parameter value (43 L mol−<sup>1</sup> s <sup>−</sup><sup>1</sup> at 10 ◦C) [41] for simulated conditions very similar to those in Table 4, there were only minor changes in the final conversion in a 10-reactor train (from 0.740 to 0.732, ~1%) and in the number of particles (from 2.74 × 10<sup>18</sup> to 2.76 × 10<sup>18</sup> L −1 , less than 0.5%), as well as even smaller changes in all the other outputs. This low sensitivity can be explained in part due to the relatively low amount of styrene used in the SBR formulation (~17% on a molar basis), but other mechanistic features of this kind of system could also play a role. For the NBR system, a similar situation occurs. A change of the propagation coefficient value *<sup>k</sup>p*,22 for AN from the non-PLP parameter (6630 L mol−<sup>1</sup> s <sup>−</sup><sup>1</sup> at 10 ◦C) [24] to the PLP parameter value (2604 L mol−<sup>1</sup> s <sup>−</sup><sup>1</sup> at 10 ◦C) [42] for an initial monomer feed of *<sup>f</sup>* <sup>1</sup> <sup>=</sup> 0.68 (wt. basis (without additional side stream of AN)), results in changes smaller than 0.2% in all the model outputs. Unfortunately, for butadiene, the available PLP-determined expression for *kp*,11 [44] is not applicable in the temperature range of interest for these processes; the same is true for the termination coefficient for butadiene, *kt*,11, but in this case, it was decided to extrapolate the Arrhenius expression in [45] out of its valid temperature range because the activation energy and, therefore, the temperature effect is much smaller for the termination than for the propagation step. Table 3 contains the values or range of values used in the simulations. For proprietary reasons, some of the details of the formulation and specific parameters used are not disclosed.


**Table 3.** Parameter values tested or used in the simulations for the styrene–butadiene rubber (SBR) and acrylonitrile–butadiene rubber (NBR) systems.


**Table 3.** *Cont.*

1 In all kinetic constants *k*<sup>x</sup> where an Arrhenius expression was available, the used notation is: *k<sup>x</sup>* = *A* exp − *E RT* with R = 1.987 cal mol−<sup>1</sup> K−<sup>1</sup> and the temperature T in K. <sup>2</sup> Values used in all the simulations in figures. <sup>3</sup> Only the used values that are different from those of the SBR system are listed.


**Table 4.** Typical recipe in parts per hundred monomer (pphm). Approximately the medium value of the range for each component was used in the simulations.

<sup>1</sup> Based on ferrous sulfate. <sup>2</sup> Based on sodium formaldehyde sulfoxylate.

#### *3.2. Model Validation*

Figure 2 shows a comparison of the model predictions with plant data for conversion (18 experimental points) in the SBR case for a 10-reactor train. The steady-state operation of each reactor (except for reactor 2 (R2) in the train) was sampled at two different times (data labeled as experiments 1 and 2); slight temperature variations (± 1 *C*) were recorded for the temperatures of the reactors at these different times, except for reactor 10, in which the temperature difference was +3 C for experiment 2. Details of the used temperature profile, as well as the formulation data used in the plant, are provided in Table 4. A similar formulation has been published before for the NBR case [29,32]. The agreement between model and experiment was quite reasonable using a single value of *kmm* = *kmp*.

**Figure 2.** Comparison of plant conversion data with model predictions for two operating conditions. The two simulated operating conditions are indicated as Model 1 and Model 2, See the conditions in Table 4.

Figure 3 shows the model predictions with a few experimental data from the plant available for the number average molecular weight *Mn*, the MWD dispersity, and the accumulated copolymer composition *F*1. The copolymer composition was measured by refractive index (ASTM D5775-95(2019)), while *M<sup>n</sup>* and the MWD dispersity were measured by size exclusion chromatography (estimating absolute values from those relative to a polystyrene standard). The reactivity ratios were slightly adjusted (~10%) with respect to those reported in the open literature to get a better agreement with the copolymer composition data (Table 3), but this was justified since the literature data were obtained at a different temperature. The agreement of the model with the experimental data was reasonable but unfortunately limited to the few available experimental data points.

Figures 4–6 show the effect of variation of surfactant on the most important outputs and process parameters with respect to the base case, which corresponded to the conditions of Experiment 1 and included a side stream of CTAs in R5. The predicted effects were all expected and easy to explain. As the surfactant increased with respect to the base value, more particles were produced (Figure 5A) and the reaction rate was higher (Figure 5C), as was the conversion (Figure 4A). At a higher conversion, the composition drift increased, as shown in Figure 4B–D. A higher number of particles resulted in smaller particle diameters (Figure 5D), although the effect is barely perceptible in the plot. Finally, higher reaction rates resulted in higher molecular weights and dispersities. The jump observed in *R*<sup>p</sup> in R6 was due to the abrupt increase in temperature from R5 to R6 (+10 C). Figure 5B shows the value of the average number of radicals per particle (e*n*) according to the model. It must be highlighted that some of the previous models ([3] and partially in [4,5]) assumed a constant value of <sup>e</sup>*<sup>n</sup>* <sup>=</sup> 0.5 (Smith–Ewart case II kinetics), but, as seen in the figure, this may have deviated from the presumably true value by 20–30%, thus having a direct impact on the estimation of the polymerization rate (see Equation (5)).

**Figure 3.** Comparison of model prediction with plant experimental data for conversion (**A**), number average molecular weight instantaneous (*M*n) (**B**), MWD dispersity (**C**), and accumulated copolymer composition *F*<sup>1</sup> (**D**). The parameter values were those in Table 3, and the reaction conditions were similar to those in Table 4.

**Figure 4.** Effect of variation of surfactant (from −45% to +45%) on conversion (**A**), remaining monomer feed composition *f* <sup>1</sup> (**B**), instantaneous copolymer composition (*F*1,inst) (**C**), and accumulated copolymer composition *F*<sup>1</sup> (**D**) along the reactor train with respect to the base case in SBR production. All compositions are weight fractions.

**Figure 5.** Effect of variation of surfactant (from −45% to +45%) on number of particles *N*<sup>p</sup> (**A**), average number of radicals <sup>e</sup>*<sup>n</sup>* (**B**), reaction rate *<sup>R</sup>*<sup>p</sup> (**C**), and particle diameter *<sup>D</sup>*<sup>p</sup> (**D**) along the reactor train with respect to the base case in SBR production.

Figure 6 shows how an increase in surfactant concentration also caused an increase in the average molecular weight (Figure 6A,B). This was due to an increased number of particles that delayed the entry of a second radical to the particle (more competition for radicals), extending its life and, therefore, its molecular weight. The increased competition for radicals also resulted in an increased dispersion in the probability of capturing a second radical, leading to increased MWD dispersity (Figure 6C). Notice that this explanation has also been put forward by other researchers, albeit in the frame of miniemulsion polymerization. [52,53]

**Figure 6.** Effect of variation of surfactant (from −45% to +45%) on the number-average molecular weight *M*n (**A**), weight-average molecular weight *M*w (**B**), and MWD dispersity (**C**) along the reactor train with respect to the base case in SBR production.

Figure 7 shows the evolution of the number of particles with zero, one, and two radicals for the simulated SBR base case. At the first reactor, due to the relative abundance of initiator radicals, the number of particles with one radical (*F*1) was higher than that with zero radicals (*F*0); however, down the train, this situation was reversed as the initiator was depleted. From the plots, it is also clear that the number of particles with two radicals (*F*2) was negligible in all the reactors, pointing to the conclusion that a 0-1 model should suffice to represent these systems. Notice that this behavior does not mean that an average value of <sup>e</sup>*<sup>n</sup>* <sup>=</sup> 0.5 would provide a good representation of the system as assumed in early models; it rather means that <sup>e</sup>*<sup>n</sup>* <sup>&</sup>lt; 0.5. However, the whole situation could change if a gel effect was included in the model.

**Figure 7.** Number of particles with zero (F\_0), one (F\_1), and two (F\_2) radicals per L of water at each reactor for the SBR base case. (**A**) Linear scale; (**B**) log scale.

Figure 8 shows the effect of the variation of the amount of CTAs on conversion, average molecular weights, and MWD dispersity. The effect was almost negligible on the conversion since the desorption of CTAs was not significant; however, the effect on the molecular weight averages and dispersity was as expected, with large increases in these quantities as the amount of CTAs decreased. The plot for a 45% decrease of CTAs is not shown since the calculation diverged for the last reactor as the system reached the gelation point (the second moment of the MWD tended to infinity).

**Figure 8.** Effect of variation of chain transfer agents (CTAs) (from −30% to +45%) on the conversion (**A**), number-average molecular weight *M*n (**B**), weight-average molecular weight *M*w (**C**), and MWD dispersity (**D**) along the reactor train with respect to the base case in SBR production.

Finally, Figure 9 shows the effect of the side stream of CTAs on the MWD. As mentioned above, in the base case, a (correction) side stream of CTAs was fed to R5. Figure 9 compares the behavior of the

train with and without the intermediate CTA feed. As previously discussed, the effect was negligible for the progress of the conversion along the train (Figure 9A), but the final *M*<sup>n</sup> and the dispersity increased from ~70,000 to ~90,000 and from ~7 to ~9, respectively, when the CTA side stream was suppressed. Similarly, the final number of branches increased from ~1 to ~1.2 in the absence of the side stream. Clearly, the intermediate CTA stream avoided an excessive increase in the molecular weight and the eventual formation of gel in the last reactors of the train that operate at high conversions.

The effects of the changes in the redox initiator system were milder than those found for the surfactant and CTAs and are shown in the Supplementary Materials (Figure S1).

Similar behavior was found for the NBR system upon changes in the surfactant, CTA, and initiator. The effects of changes in the surfactant for an NBR system are shown in Figures 10–12 for the same important outputs and process parameters shown for the SBR system in Figures 4–6. The variations of the surfactant were made with respect to a base case, which corresponded to a composition feed of *f* <sup>1</sup> = 0.68 (wt. basis) and weight feed ratios of 1.25/0.082/0.35/100 for S/Iw/CTAs/monomers. The base case also included a side stream of an additional 20% flow rate of CTAs at R6. The used temperature was 10 ◦C in all reactors except in R1 (12 ◦C) and R10 (15 ◦C). The effects were very similar to those found for the SBR system, with a few exceptions. In Figure 10, the butadiene composition of the remaining monomer *f* <sup>1</sup>, as well as those of the copolymer (accumulated and instantaneous (*F*<sup>1</sup> and

*F*1,inst, respectively)), increased instead of decreasing like in the SBR case. This was expected given the differences in the reactivity ratios between the two systems. In this region of the monomer feed composition (above the azeotropic point for NBR) in the NBR case, acrylonitrile was consumed faster than butadiene (see the Mayo–Lewis plot for this system in Figure S2), and the composition of the remaining monomer therefore increased in butadiene with conversion, also leading to a richer butadiene composition in the copolymer.

**Figure 10.** Effect of variation of surfactant (from −45% to +45%) on conversion (**A**), remaining monomer feed composition *f* <sup>1</sup> (**B**), instantaneous copolymer composition (*F*1,inst) (**C**), and accumulated copolymer composition *F*<sup>1</sup> (**D**) along the reactor train with respect to the base case in NBR production. All compositions are weight fractions.

**Figure 11.** Effect of variation of surfactant (from −45% to +45%) on number of particles *N*<sup>p</sup> (**A**), average number of radicals <sup>e</sup>*<sup>n</sup>* (**B**), reaction rate *<sup>R</sup>*<sup>p</sup> (**C**), and particle diameter *<sup>D</sup>*<sup>p</sup> (**D**) along the reactor train with respect to the base case in NBR production.

**Figure 12.** Effect of variation of surfactant (from −45% to +45%) on the number-average molecular weight *M*n (**A**), weight-average molecular weight *M*w (**B**), and MWD dispersity (**C**) along the reactor train with respect to the base case in NBR production.

The effects of the surfactant changes on *<sup>N</sup>p*, *<sup>R</sup>p*, <sup>e</sup>*n*, and *<sup>D</sup><sup>p</sup>* shown in Figure <sup>11</sup> for the NBR system were qualitatively very similar to those observed and already discussed for the SBR case (Figure 5). The only minor difference was that the effects were smaller for <sup>e</sup>*<sup>n</sup>* in the NBR than in the SBR case. Another difference is that in the NBR case, there was not a jump in the reaction rate as in the SBR case; however, this was simply because in the former case, there was not a jump in an intermediate reactor temperature as in the latter. Figure 12 shows similar effects of the surfactant concentration for the NBR system compared to those in the SBR system (Figure 6). A difference was that in the NBR case, the molecular weights reached a peak in R5 to decrease in R6 due to the side stream of CTAs. This effect was less important in the SBR base case, but in that system, the higher temperatures in the last reactors also significantly contributed to higher molecular weights. Additionally, the molecular weights and dispersities were lower for NBR compared to those exhibited by the SBR system.

The effects of changes in the CTA feed with respect to the NBR base case are shown in Figure 13, which is comparable to Figure 8 for the SBR system. Again, the behavior of both systems was qualitatively similar with the only exception of the peak exhibited by the molecular weights in R5 for the NBR case due to the side stream of CTAs in R6, as in Figure 12, and the different temperature profiles used for the SBR and the NBR cases.

**Figure 13.** Effect of variation of CTAs (from −45% to +45%) on the conversion (**A**), number-average molecular weight *M*n (**B**), weight-average molecular weight *M*w (**C**), and MWD dispersity (**D**) along the reactor train with respect to the base case in NBR production.

More interesting behavior for the NBR system was related to the drift in copolymer composition along the train. Given the reactivity ratios of this system (0.30 and 0.04), the system exhibited an azeotropic composition at *f* 1,az = *F*1,az = 0.578 (molar basis, 0.582 weight basis) (see the *F*<sup>1</sup> vs. *f* <sup>1</sup> Mayo–Lewis curve in the Supplementary Materials, Figure S2) Therefore, composition drift depends on what side of the azeotrope the composition of the feed is located. If the feed composition is above *f* 1,az, as in most commercial grades of NBR, without correction, the instantaneous product composition will initially deviate towards copolymers with a lower Bd composition than the feed (R1); however, as the reaction mixture proceeds along the reactor train, the Bd-enriched environment will tend to increase the Bd content in the copolymer. This can be seen in Figure 14 for a product grade in which *f* <sup>1</sup> in the feed was 0.68 (wt. basis). Since composition drift may lead to a product with heterogeneous composition or a product deviated from the target average composition, it is common to introduce a side feed stream of AN in an intermediate reactor of the train to "correct" or attenuate the composition drift. Figure 10 also shows the composition profile when an additional AN side stream (20% of the initial feed) was fed to reactor R7. Without the correction, the final *F*<sup>1</sup> was above 0.66; on the other hand, with the side feed stream at R7, the final *F*<sup>1</sup> was slightly below 0.65. A second example is shown

in Figure 15; in this case, the composition of the feed was *f* <sup>1</sup> = 0.5 on a wt. basis, a composition which was located on the other side of the azeotrope (see the *F*<sup>1</sup> vs. *f* <sup>1</sup> curve in the Supplementary Materials). The behavior of this system was the opposite of that in the previous example, exhibiting a composition drift that initially produced copolymers with a higher Bd composition than the feed (R1, *F*<sup>1</sup> ~0.58) and a decreasing Bd content in the copolymer along the reactor train to end up with an *F*<sup>1</sup> slightly above 0.56, although the instantaneous *F*<sup>1</sup> fell below 0.54 at R10. Figure 15 also shows how this drift could be attenuated by feeding a side stream of Bd (20% of the initial feed). This correction worked very well, resulting in a final average composition near the target (*F*<sup>1</sup> ~0.58) while avoiding large deviations in the instantaneous composition produced in each reactor (better composition homogeneity).

**Figure 14.** Effect of the operation of the reactor train for NBR production with (corrected) and without a side stream of additional AN in R7 (20% with respect to the AN feed in R1) on conversion (**A**) and instantaneous butadiene copolymer composition (*F*1,inst) and accumulated butadiene copolymer composition (*F*<sup>1</sup> ) (**B**). The initial monomer feed was *f* <sup>1</sup> = 0.68 (wt. basis).

**Figure 15.** Effect of the operation of the reactor train for NBR production with (corrected) and without a side stream of additional Bd in R7 (20% with respect to the Bd feed in R1) on conversion (**A**) and instantaneous butadiene copolymer composition (*F*1,inst) and accumulated butadiene copolymer composition (*F*<sup>1</sup> ) (**B**). The initial monomer feed is *f* <sup>1</sup> = 0.5 (wt. basis).

#### **4. Conclusions**

A single mechanistic model was built and used to describe the production of both SBR and NBR elastomers in trains of emulsion polymerization CSTRs. The aim was to balance the level of complexity of the model with its capability of describing the effects observed at the plant on the process parameters in a quantitative way and the product quality upon variable changes. Some novel aspects of this model compared to previous models are summarized below, along with additional comments:


The model is presently used in the plant to design changes of configurations of the reactor train, such as the number of reactors in the train (up to 16 reactors in the existing facility), the reactor volume, the location and flows of side streams of CTAs and monomers, and grade recipes, among other variables, to maximize the plant operation versatility maintaining high productivity and product quality. The model is quite flexible, and it has been implemented with the help of a friendly graphical interphase. Current work involves the generation of product quality correlations that link model outputs, e.g., *M*<sup>n</sup> and *M*w, and number of branches, with product performance parameters, such as Mooney viscosity.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2227-9717/8/11/1508/s1, File with figures and additional equations.

**Author Contributions:** Conceptualization, E.S.-G. and J.M.I.-M.; methodology, E.S.-G., R.I.-M., and J.M.I.-M; software, E.S.-G. and R.I.-M.; validation, E.S.-G., R.I.-M., and J.M.I.-M.; writing—original draft preparation, E.S.-G. and J.M.I.-M; writing—review and editing, E.S.-G. and J.M.I.-M.; project administration, E.S.-G. and J.M.I.-M; funding acquisition, E.S.-G. and J.M.I.-M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Dynasol Elastomers and the publication fee by CIQA.

**Acknowledgments:** The authors thank Dynasol Elastomers for industrial data and CIQA for administrative support.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
